text:
"""
Thermometer state machine:

@startuml
state ThermometerStateMachine {
    state OffState
    state OnState {
        state CelsiusState
        state FahrenheitState
        state KelvinState

        [*] --> CelsiusState
        CelsiusState : entry / set converter to celsius.
        CelsiusState -> FahrenheitState : UnitsEvent
        FahrenheitState: entry / set converter to fahrenheit.
        FahrenheitState --> KelvinState : UnitsEvent
        KelvinState : entry / set converter to kelvin.
        KelvinState --> CelsiusState : UnitsEvent
    }

    [*] --> OffState

    OffState : entry / set display to "off".
    OffState --> OnState[H*] : PowerEvent
    OnState : entry / set display to "--".
    OnState : TemperatureEvent / set display to converted temperature value.
    OnState --> OffState : PowerEvent
}
@enduml
"""

using Logging, Test, HSM

#####
# Temperature converter functions.
#####

celsius(t::Float64) = string(round(t, digits=1)) * "°C"
fahrenheit(t::Float64) = string(round((t * 1.8) + 32.0, digits=1)) * "°F"
kelvin(t::Float64) = string(round(t + 273.15, digits=1)) * "K"

#####
# Thermometer model.
#####

mutable struct Thermometer
    display::String
    temperature::Union{Float64, Nothing}
    convert::Function

    function Thermometer(display="", temperature=nothing, convert=celsius)
        new(display, temperature, convert)
    end
end

thermometer = Thermometer()

#####
# Events.
#####

struct PowerEvent <: HSM.AbstractHsmEvent
end

struct UnitsEvent <: HSM.AbstractHsmEvent
end

struct TemperatureEvent <:HSM.AbstractHsmEvent
    temperature::Float64
end

#####
# State machine.
#####

macro thermometer_state(name)
    return :(
        struct $name <: HSM.AbstractHsmState
            state_info::HSM.HsmStateInfo
            thermometer::Thermometer

            $name(parent, thermometer) = new(HSM.HsmStateInfo(parent), thermometer)
        end
    )
end

@thermometer_state(ThermometerStateMachine)
@thermometer_state(OffState)
@thermometer_state(OnState)
@thermometer_state(CelsiusState)
@thermometer_state(FahrenheitState)
@thermometer_state(KelvinState)

thermometer_state_machine = ThermometerStateMachine(nothing, thermometer)
off_state = OffState(thermometer_state_machine, thermometer)
on_state = OnState(thermometer_state_machine, thermometer)
celsius_state = CelsiusState(on_state, thermometer)
fahrenheit_state = FahrenheitState(on_state, thermometer)
kelvin_state = KelvinState(on_state, thermometer)

# Root state machine state event handlers.
# The root state machine state must handle all possible user-defined events --
# defining `on_event!()` functions for each -- which return `true` to avoid
# UnhandledExceptionErrors being thrown.

function HSM.on_initialize!(state::ThermometerStateMachine)
    @debug "on_initialize!(ThermometerStateMachine)"
    HSM.transition_to_state!(thermometer_state_machine, off_state)
end

function HSM.on_event!(state::ThermometerStateMachine, event::PowerEvent)
    @debug "on_event!(ThermometerStateMachine, PowerEvent)"
    return true
end

function HSM.on_event!(state::ThermometerStateMachine, event::UnitsEvent)
    @debug "on_event!(ThermometerStateMachine, UnitsEvent)"
    return true
end

function HSM.on_event!(state::ThermometerStateMachine, event::TemperatureEvent)
    @debug "on_event!(ThermometerStateMachine, TemperatureEvent)"
    return true
end

# Off state event handlers.

function HSM.on_entry!(state::OffState)
    @debug "on_entry!(OffState)"
    state.thermometer.display = "off"
    state.thermometer.temperature = nothing
end

function HSM.on_event!(state::OffState, event::PowerEvent)
    @debug "on_event!(OffState, PowerEvent)"
    HSM.transition_to_deep_history!(thermometer_state_machine, on_state)
    return true
end

# On state event handlers.

function HSM.on_initialize!(state::OnState)
    @debug "on_initialize!(OnState)"
    HSM.transition_to_state!(thermometer_state_machine, celsius_state)
end

function HSM.on_entry!(state::OnState)
    @debug "on_entry!(OnState)"
    state.thermometer.display = "--"
end

function HSM.on_event!(state::OnState, event::PowerEvent)
    @debug "on_event!(OnState, PowerEvent)"
    HSM.transition_to_state!(thermometer_state_machine, off_state)
    return true
end

function HSM.on_event!(state::OnState, event::TemperatureEvent)
    @debug "on_event!(OnState, TemperatureEvent)"
    state.thermometer.display = state.thermometer.convert(event.temperature)
    state.thermometer.temperature = event.temperature
    return true
end

# Celsius state event handlers.

function HSM.on_entry!(state::CelsiusState)
    @debug "on_entry!(CelsiusState)"
    state.thermometer.convert = celsius
    if !isnothing(state.thermometer.temperature)
        state.thermometer.display = state.thermometer.convert(state.thermometer.temperature)
    end
end

function HSM.on_event!(state::CelsiusState, event::UnitsEvent)
    @debug "on_event!(CelsiusState, UnitsEvent)"
    HSM.transition_to_state!(thermometer_state_machine, fahrenheit_state)
    return true
end

# Fahrenheit state event handlers.

function HSM.on_entry!(state::FahrenheitState)
    @debug "on_entry!(FahrenheitState)"
    state.thermometer.convert = fahrenheit
    if !isnothing(state.thermometer.temperature)
        state.thermometer.display = state.thermometer.convert(state.thermometer.temperature)
    end
end

function HSM.on_event!(state::FahrenheitState, event::UnitsEvent)
    @debug "on_event!(FahrenheitState, UnitsEvent)"
    HSM.transition_to_state!(thermometer_state_machine, kelvin_state)
    return true
end

# Kelvin state event handlers.

function HSM.on_entry!(state::KelvinState)
    @debug "on_entry!(KelvinState)"
    state.thermometer.convert = kelvin
    if !isnothing(state.thermometer.temperature)
        state.thermometer.display = state.thermometer.convert(state.thermometer.temperature)
    end
end

function HSM.on_event!(state::KelvinState, event::UnitsEvent)
    @debug "on_event!(KelvinState, UnitsEvent)"
    HSM.transition_to_state!(thermometer_state_machine, celsius_state)
    return true
end

#####
# Exercise / Test state machine.
#####

@testset "State machine initialization" begin
    # Initialize the state machine.
    HSM.transition_to_state!(thermometer_state_machine, thermometer_state_machine)

    @test HSM.active_state(thermometer_state_machine) == off_state
    @test thermometer.display == "off"
end

@testset "Off State tests" begin
    # Test ignored events.
    HSM.handle_event!(thermometer_state_machine, UnitsEvent())

    @test HSM.active_state(thermometer_state_machine) == off_state
    @test thermometer.display == "off"

    HSM.handle_event!(thermometer_state_machine, TemperatureEvent(10.0))

    @test HSM.active_state(thermometer_state_machine) == off_state
    @test thermometer.display == "off"

    # Test power event.
    HSM.handle_event!(thermometer_state_machine, PowerEvent())

    @test HSM.active_state(thermometer_state_machine) == celsius_state
    @test thermometer.display == "--"
end

@testset "On State tests" begin
    # Test units event.
    HSM.handle_event!(thermometer_state_machine, UnitsEvent())

    @test HSM.active_state(thermometer_state_machine) == fahrenheit_state
    @test thermometer.display == "--"

    HSM.handle_event!(thermometer_state_machine, UnitsEvent())

    @test HSM.active_state(thermometer_state_machine) == kelvin_state
    @test thermometer.display == "--"

    HSM.handle_event!(thermometer_state_machine, TemperatureEvent(0.0))

    @test thermometer.display == "273.2K"

    HSM.handle_event!(thermometer_state_machine, UnitsEvent())

    @test HSM.active_state(thermometer_state_machine) == celsius_state
    @test thermometer.display == "0.0°C"

    HSM.handle_event!(thermometer_state_machine, UnitsEvent())

    @test HSM.active_state(thermometer_state_machine) == fahrenheit_state
    @test thermometer.display == "32.0°F"

    HSM.handle_event!(thermometer_state_machine, PowerEvent())

    @test HSM.active_state(thermometer_state_machine) == off_state
    @test thermometer.display == "off"
end

@testset "On - Off - On retains memory" begin
    HSM.handle_event!(thermometer_state_machine, PowerEvent())

    @test HSM.active_state(thermometer_state_machine) == fahrenheit_state
    @test thermometer.display == "--"
end


text:
function Base.show(io::IO, A::AbstractGeoArray)
    l = nameof(typeof(A))
    printstyled(io, nameof(typeof(A)); color=:blue)
    if label(A) != ""
        print(io, " (named ")
        printstyled(io, label(A); color=:blue)
        print(io, ")")
    end
    print(io, " with dimensions:\n")
    for d in dims(A)
        print(io, " ", d, " (length ", length(d), ") \n")
    end
    if !isempty(refdims(A))
        print(io, "and referenced dimensions:\n")
        for d in refdims(A)
            print(io, " ", d, "\n")
        end
    end
    A isa DiskGeoArray && print(io, "\n  From file: $(filename(A))")
end

function Base.show(io::IO, stack::AbstractGeoStack)
    n_fields = length(keys(stack))
    fields_str = n_fields == 1 ? "field" : "fields"
    printstyled(io, "$(Base.typename(typeof(stack)))", color=:blue)
    print(io, " with $n_fields $fields_str")
    stack isa DiskGeoStack && print(io, ": $(filename(stack))")
    print(io, '\n')

    for var in keys(stack)
        printstyled(io, " $var", color=:green)

        field_dims = dims(stack, var)
        n_dims = length(field_dims)
        dims_str = n_dims == 1 ? "dimension" : "dimensions"
        print(io, " with $n_dims $dims_str: ")
        if n_dims > 0
            for (d, dim) in enumerate(field_dims)
                printstyled(io, "$(name(dim))", color=:red)
                d != length(field_dims) && print(io, ", ")
            end
            print(io, " (size ")
            for (d, dim) in enumerate(field_dims)
                print(io, "$(length(dim))")
                d != length(field_dims) && print(io, '×')
            end
            print(io, ')')
        end
        print(io, '\n')
    end

    n_windows = length(window(stack))
    if n_windows > 0
        print(io, "and with window:\n")
        for window in window(stack)
            print(io, ' ')
            show(window)
            print(io, '\n')
        end
    end

    if !isnothing(metadata(stack))
        n_metadata = length(metadata(stack))
        entries_str = n_metadata == 1 ? "entry" : "entries"
        if n_metadata > 0
            print(io, "and $n_metadata metadata $entries_str:\n")
            display(stack.metadata)
        end
    end
end

ext: jl


text:
# Contains constants and constructors for SAC file precision, and defaults
# for routines and endianness

# SAC types
const SACFloat = Float32
const SACChar = String
const SACInt = Int32
const SACBool = Bool
# Constructors
sacfloat(x) = map(Float32, x)
sacint(x) = map(Int32, x)
const sacchar = ascii
sacbool(x) = x != 0
# Length of SAC floats and ints
const sac_byte_len = 4
# Length of SAC character headers (except kevnm, which is twice the length)
const saccharlen = 8
# SAC file version number
const sac_ver_num = SACInt(6)
# Whether this machine is big- or little endian.  SAC files are meant to be big-endian,
# so this determines whether a file is 'native-endian' or not.
const machine_is_little_endian = Base.ENDIAN_BOM == 0x04030201

# Convert a number into a SACChar
sacstring(x, maxlen=saccharlen) = sacchar(string(x)[1:minimum((length(string(x)),maxlen))]*" "^(maximum((0,maxlen-length(string(x))))))

# SAC unset values
const sac_rnull = -12345.
const sac_inull = -12345
const sac_cnull = "-12345"

# Default values for filtering
const sac_npoles = 2
const sac_passes = 1

# For SAC/BRIS (MacSAC), files are always big-endian, so set this appropriately
const sac_force_swap = machine_is_little_endian

# Flag for verbosity
const sac_verbose = Ref(true)

# Lists of SAC headers as symbols
const sac_float_hdr = [:delta, :depmin, :depmax, :scale, :odelta, :b,
    :e, :o, :a, :internal0, :t0, :t1,
    :t2, :t3, :t4, :t5, :t6, :t7,
    :t8, :t9, :f, :resp0, :resp1, :resp2,
    :resp3, :resp4, :resp5, :resp6, :resp7, :resp8,
    :resp9, :stla, :stlo, :stel, :stdp, :evla,
    :evlo, :evel, :evdp, :mag, :user0, :user1,
    :user2, :user3, :user4, :user5, :user6, :user7,
    :user8, :user9, :dist, :az, :baz, :gcarc,
    :internal1, :internal2, :depmen, :cmpaz, :cmpinc, :xminimum,
    :xmaximum, :yminimum, :ymaximum, :unused1, :unused2, :unused3,
    :unused4, :unused5, :unused6, :unused7]
const sac_int_hdr = [:nzyear, :nzjday, :nzhour, :nzmin, :nzsec, :nzmsec,
    :nvhdr, :norid, :nevid, :npts, :internal3, :nwfid,
    :nxsize, :nysize, :unused8, :iftype, :idep, :iztype,
    :unused9, :iinst, :istreg, :ievreg, :ievtyp, :iqual,
    :isynth, :imagtyp, :imagsrc, :unused10, :unused11, :unused12,
    :unused13, :unused14, :unused15, :unused16, :unused17]
const sac_bool_hdr = [:leven, :lpspol, :lovrok, :lcalda, :unused18]
const sac_char_hdr = [:kstnm, :kevnm, :khole, :ko, :ka, :kt0,
    :kt1, :kt2, :kt3, :kt4, :kt5, :kt6,
    :kt7, :kt8, :kt9, :kf, :kuser0, :kuser1,
    :kuser2, :kcmpnm, :knetwk, :kdatrd, :kinst]
const sac_all_hdr = [sac_float_hdr; sac_int_hdr; sac_bool_hdr; sac_char_hdr]

# Where in the file the NVHDR value is
const sac_nvhdr_pos = length(sac_float_hdr) + findfirst(sac_int_hdr .== :nvhdr) - 1

# Length in bytes of total SAC header, accounting for double-length kevnm
const sac_header_len = sac_byte_len*(length(sac_float_hdr) + length(sac_int_hdr) +
    length(sac_bool_hdr)) + saccharlen*(length(sac_char_hdr) + 1)


text:
module Display

using UUIDs
import LibGit2

using Pkg3.Types
import Pkg3: @info, Nothing

const colors = Dict(
    ' ' => :white,
    '+' => :light_green,
    '-' => :light_red,
    '↑' => :light_yellow,
    '~' => :light_yellow,
    '↓' => :light_magenta,
    '?' => :red,
)
const color_dark = :light_black

function git_file_stream(repo::LibGit2.GitRepo, spec::String; fakeit::Bool=false)::IO
    blob = try LibGit2.GitBlob(repo, spec)
    catch err
        err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND || rethrow(err)
        fakeit && return devnull
    end
    return IOBuffer(LibGit2.rawcontent(blob))
end

function status(ctx::Context, mode::PackageMode, use_as_api=false)
    env = ctx.env
    project₀ = project₁ = env.project
    manifest₀ = manifest₁ = env.manifest
    diff = nothing

    if env.git != nothing
        git_path = LibGit2.path(env.git)
        project_path = relpath(env.project_file, git_path)
        manifest_path = relpath(env.manifest_file, git_path)
        project₀ = read_project(git_file_stream(env.git, "HEAD:$project_path", fakeit=true))
        manifest₀ = read_manifest(git_file_stream(env.git, "HEAD:$manifest_path", fakeit=true))
    end
    if mode == PKGMODE_PROJECT || mode == PKGMODE_COMBINED
        # TODO: handle project deps missing from manifest
        m₀ = filter_manifest(in_project(project₀["deps"]), manifest₀)
        m₁ = filter_manifest(in_project(project₁["deps"]), manifest₁)
        use_as_api || @info("Status $(pathrepr(env, env.project_file))")
        diff = manifest_diff(m₀, m₁)
        use_as_api || print_diff(diff)
    end
    if mode == PKGMODE_MANIFEST
        use_as_api || @info("Status $(pathrepr(env, env.manifest_file))")
        diff = manifest_diff(manifest₀, manifest₁)
        use_as_api || print_diff(diff)
    elseif mode == PKGMODE_COMBINED
        p = not_in_project(merge(project₀["deps"], project₁["deps"]))
        m₀ = filter_manifest(p, manifest₀)
        m₁ = filter_manifest(p, manifest₁)
        c_diff = filter!(x->x.old != x.new, manifest_diff(m₀, m₁))
        if !isempty(c_diff)
            use_as_api || @info("Status $(pathrepr(env, env.manifest_file))")
            use_as_api || print_diff(c_diff)
            diff = Base.vcat(c_diff, diff)
        end
    end
    return diff
end

function print_project_diff(env₀::EnvCache, env₁::EnvCache)
    pm₀ = filter_manifest(in_project(env₀.project["deps"]), env₀.manifest)
    pm₁ = filter_manifest(in_project(env₁.project["deps"]), env₁.manifest)
    diff = filter!(x->x.old != x.new, manifest_diff(pm₀, pm₁))
    if isempty(diff)
        printstyled(color = color_dark, " [no changes]\n")
    else
        print_diff(diff)
    end
end

function print_manifest_diff(env₀::EnvCache, env₁::EnvCache)
    diff = manifest_diff(env₀.manifest, env₁.manifest)
    diff = filter!(x->x.old != x.new, diff)
    if isempty(diff)
        printstyled(color = color_dark, " [no changes]\n")
    else
        print_diff(diff)
    end
end

struct VerInfo
    hash::Union{SHA1,Nothing}
    path::Union{String,Nothing}
    ver::Union{VersionNumber,Nothing}
    pinned::Bool
end

vstring(a::VerInfo) =
    string(a.ver == nothing ? "[$(string(a.hash)[1:16])]" : "v$(a.ver)",
           a.pinned == true ? "⚲" : "",
           a.path != nothing ? " [$(a.path)]" : "")

Base.:(==)(a::VerInfo, b::VerInfo) =
    a.hash == b.hash && a.ver == b.ver && a.pinned == b.pinned

≈(a::VerInfo, b::VerInfo) = a.hash == b.hash &&
    (a.ver == nothing || b.ver == nothing || a.ver == b.ver) &&
    (a.pinned == b.pinned)

struct DiffEntry
    uuid::UUID
    name::String
    old::Union{VerInfo,Nothing}
    new::Union{VerInfo,Nothing}
end

function print_diff(io::IO, diff::Vector{DiffEntry})
    same = all(x.old == x.new for x in diff)
    for x in diff
        warnings = String[]
        if x.old != nothing && x.new != nothing
            if x.old ≈ x.new
                verb = ' '
                vstr = vstring(x.new)
            else
                if x.old.hash != x.new.hash && x.old.ver != x.new.ver
                    verb = x.old.ver == nothing || x.new.ver == nothing ||
                           x.old.ver == x.new.ver ? '~' :
                           x.old.ver < x.new.ver  ? '↑' : '↓'
                elseif x.old.ver == x.new.ver && x.old.pinned != x.new.pinned
                    verb = '~'
                else
                    verb = '?'
                    msg = x.old.hash == x.new.hash ?
                        "hashes match but versions don't: $(x.old.ver) ≠ $(x.new.ver)" :
                        "versions match but hashes don't: $(x.old.hash) ≠ $(x.new.hash)"
                    push!(warnings, msg)
                end
                vstr = (x.old.ver == x.new.ver && x.old.pinned == x.new.pinned) ?
                       vstring(x.new) :
                       vstring(x.old) * " ⇒ " * vstring(x.new)
            end
        elseif x.new != nothing
            verb = '+'
            vstr = vstring(x.new)
        elseif x.old != nothing
            verb = '-'
            vstr = vstring(x.old)
        else
            verb = '?'
            vstr = "[unknown]"
        end
        v = same ? "" : " $verb"
        printstyled(color = color_dark, " [$(string(x.uuid)[1:8])]")
        printstyled(color = colors[verb], "$v $(x.name) $vstr\n")
    end
end
print_diff(diff::Vector{DiffEntry}) = print_diff(stdout, diff)

function manifest_by_uuid(manifest::Dict)
    entries = Dict{UUID,Dict}()
    for (name, infos) in manifest, info in infos
        uuid = UUID(info["uuid"])
        haskey(entries, uuid) && @warn("Duplicate UUID in manifest: $uuid")
        entries[uuid] = merge(info, Dict("name" => name))
    end
    return entries
end

function name_ver_info(info::Dict)
    name = info["name"]
    hash = haskey(info, "git-tree-sha1") ? SHA1(info["git-tree-sha1"])    : nothing
    ver  = haskey(info, "version")       ? VersionNumber(info["version"]) : nothing
    path =  get(info, "path", nothing)
    pin  =  get(info, "pinned", false)
    name, VerInfo(hash, path, ver, pin)
end

function manifest_diff(manifest₀::Dict, manifest₁::Dict)
    diff = DiffEntry[]
    entries₀ = manifest_by_uuid(manifest₀)
    entries₁ = manifest_by_uuid(manifest₁)
    for uuid in union(keys(entries₀), keys(entries₁))
        name₀ = name₁ = v₀ = v₁ = nothing
        haskey(entries₀, uuid) && ((name₀, v₀) = name_ver_info(entries₀[uuid]))
        haskey(entries₁, uuid) && ((name₁, v₁) = name_ver_info(entries₁[uuid]))
        name₀ == nothing && (name₀ = name₁)
        name₁ == nothing && (name₁ = name₀)
        if name₀ == name₁
            push!(diff, DiffEntry(uuid, name₀, v₀, v₁))
        else
            push!(diff, DiffEntry(uuid, name₀, v₀, nothing))
            push!(diff, DiffEntry(uuid, name₁, nothing, v₁))
        end
    end
    sort!(diff, by=x->(x.name, x.uuid))
end

function filter_manifest!(predicate, manifest::Dict)
    empty = String[]
    for (name, infos) in manifest
        filter!(infos) do info
            predicate(name, info)
        end
        isempty(infos) && push!(empty, name)
    end
    for name in empty
        pop!(manifest, name)
    end
    return manifest
end
filter_manifest(predicate, manifest::Dict) =
    filter_manifest!(predicate, deepcopy(manifest))

# This is precompilable, an anonymous function is not.
struct InProject{D <: Dict}
    deps::D
    neg::Bool
end
function (ip::InProject)(name::String, info::Dict)
    v = haskey(ip.deps, name) && haskey(info, "uuid") && ip.deps[name] == info["uuid"]
    return ip.neg ? !v : v
end
in_project(deps::Dict) = InProject(deps, false)
not_in_project(deps::Dict) = InProject(deps, true)

end # module

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text: Dict{String,Any}("multiline" => Dict{String,Any}("value" => "This string\r\nhas ' a quote character\r\nand more than\r\none newline\r\nin it.","type" => "string"),"firstnl" => Dict{String,Any}("value" => "This string has a ' quote character.","type" => "string"),"oneline" => Dict{String,Any}("value" => "This string has a ' quote character.","type" => "string"))

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using FactCheck
using ReadWriteLocks

facts("Single-threaded tests") do
    context("Initialization") do
        rwlock = ReadWriteLock()

        @fact rwlock.read_lock.rwlock --> rwlock
        @fact rwlock.write_lock.rwlock --> rwlock
        @fact rwlock.readers --> 0
        @fact rwlock.writer --> false
        @fact read_lock(rwlock) --> rwlock.read_lock
        @fact write_lock(rwlock) --> rwlock.write_lock
    end
end

facts("Two-threaded tests") do
    context("Read locks") do
        NUM_LOCKS = 10

        context("$NUM_LOCKS locks") do
            rwlock = ReadWriteLock()
            rlock = read_lock(rwlock)

            @async begin
                @fact rwlock.writer --> false
                @fact rwlock.readers --> 0
                for i = 1:NUM_LOCKS
                    @fact lock!(rlock) --> nothing
                    @fact rwlock.readers --> i
                end
            end

            sleep(1)

            @fact rwlock.readers --> NUM_LOCKS

            @async begin
                @fact rwlock.writer --> false
                @fact rwlock.readers --> NUM_LOCKS
                for i = NUM_LOCKS:-1:1
                    @fact unlock!(rlock) --> nothing
                    @fact rwlock.readers --> i - 1
                end
                @fact rwlock.readers --> 0
            end

            sleep(1)

            @fact rwlock.readers --> 0
        end
    end

    context("Write locks") do
        context("two locks") do
            rwlock = ReadWriteLock()
            wlock = write_lock(rwlock)

            c = Channel{Symbol}(1)
            put!(c, :pretest)

            @async begin
                @fact rwlock.writer --> false
                @fact rwlock.readers --> 0
                @fact lock!(wlock) --> nothing
                @fact rwlock.writer --> true
                @fact rwlock.readers --> 0
                @fact take!(c) --> :pretest
                put!(c, :prelock)
                @fact lock!(wlock) --> nothing

                # this code should never be reached
                @fact take!(c) --> :prelock
                put!(c, :postlock)
            end

            sleep(1)

            @fact take!(c) --> :prelock
            put!(c, :posttest)

            @fact rwlock.writer --> true
            @fact rwlock.readers --> 0
        end

        context("unlock") do
            rwlock = ReadWriteLock()
            wlock = write_lock(rwlock)

            c = Channel{Symbol}(1)
            put!(c, :pretest)

            @async begin
                @fact rwlock.writer --> false
                @fact rwlock.readers --> 0
                @fact lock!(wlock) --> nothing
                @fact rwlock.writer --> true
                @fact rwlock.readers --> 0
                @fact take!(c) --> :pretest
                put!(c, :preunlock)
                @fact unlock!(wlock) --> nothing
                @fact take!(c) --> :preunlock
                put!(c, :postunlock)
            end

            sleep(1)

            @fact take!(c) --> :postunlock
            put!(c, :posttest)

            @fact rwlock.writer --> false
            @fact rwlock.readers --> 0
        end
    end

    context("read and write locks") do
        context("write then read") do
            rwlock = ReadWriteLock()
            wlock = write_lock(rwlock)
            rlock = read_lock(rwlock)

            c = Channel{Symbol}(1)
            put!(c, :pretest)

            @async begin
                @fact rwlock.writer --> false
                @fact rwlock.readers --> 0
                @fact lock!(wlock) --> nothing
                @fact rwlock.writer --> true
                @fact rwlock.readers --> 0
                @fact take!(c) --> :pretest
                put!(c, :prelock)
                @fact lock!(rlock) --> nothing

                # this code should never be reached
                @fact take!(c) --> :prelock
                put!(c, :postlock)
            end

            sleep(1)

            @fact take!(c) --> :prelock
            put!(c, :posttest)

            @fact rwlock.writer --> true
            @fact rwlock.readers --> 0
        end

        context("read then write") do
            rwlock = ReadWriteLock()
            wlock = write_lock(rwlock)
            rlock = read_lock(rwlock)

            c = Channel{Symbol}(1)
            put!(c, :pretest)

            @async begin
                @fact rwlock.writer --> false
                @fact rwlock.readers --> 0
                @fact lock!(rlock) --> nothing
                @fact rwlock.writer --> false
                @fact rwlock.readers --> 1
                @fact take!(c) --> :pretest
                put!(c, :prelock)
                @fact lock!(wlock) --> nothing

                # this code should never be reached
                @fact take!(c) --> :prelock
                put!(c, :postlock)
            end

            sleep(1)

            @fact take!(c) --> :prelock
            put!(c, :posttest)

            @fact rwlock.writer --> false
            @fact rwlock.readers --> 1
        end
    end

end


text:
using Equate
using Statistics, StatsBase

# methods = [
#     "Tucker"
#     "LevineCongeneric"
#     "ChainedLinear"
#     "FrequencyEstimation"
#     "BraunHolland"
#     "ChainedEquipercentile"
# ]

@testset "Tucker" begin
    res = Tucker(neattabX, neattabY)
    @test round(coef(res)[1]; digits = 4) == 1.0292
    @test round(coef(res)[2]; digits = 4) == .5378
end

@testset "LevineObservedScore (external, not show in K&B) " begin
    res = LevineObservedScore(neattabX, neattabY; common = :external)
    @test round(coef(res)[1]; digits = 4) == 1.0168
    @test round(coef(res)[2]; digits = 4) == .3257
end

@testset "LevineObservedScore (internal) " begin
    res = LevineObservedScore(neattabX, neattabY; common = :internal)
    @test round(coef(res)[1]; digits = 4) == 1.0110
    @test round(coef(res)[2]; digits = 4) == .2514
end

@testset "LevineTrueScore (external, not show in K&B) " begin
    res = LevineTrueScore(neattabX, neattabY; common = :external)
    @test round(coef(res)[1]; digits = 4) == 1.0165
    @test round(coef(res)[2]; digits = 4) == .331
end

@testset "LevineTrueScore (internal) " begin
    res = LevineTrueScore(neattabX, neattabY; common = :internal)
    @test round(coef(res)[1]; digits = 4) == 1.0086
    @test round(coef(res)[2]; digits = 4) == .2912
end

@testset "ChainedLinear" begin
    res = ChainedLinear(neattabX, neattabY)
    @test round(coef(res)[1]; digits = 4) == 1.0213
    @test round(coef(res)[2]; digits = 4) == .3937
end

# Invoke the percentile rank methods
# The results in Kolen & Brennan seems to be wrong.
@testset "BraunHolland (R-equate)" begin
    res = BraunHolland(neattabX, neattabY)
    @test round(coef(res)[1]; digits = 4) == 1.0067
    @test round(coef(res)[2]; digits = 4) == 0.8976
    # popXonY = ExpandTable(res.table.lYx, neattabX.tabX.freq)
    # @test round(mean(popXonY); digits = 4) == 16.8329
    # @test round(std(popXonY; corrected = false); digits = 4) == 6.6017
    # @test round(skewness(popXonY); digits = 4) == .5799
end

# @testset "BraunHolland (Kolen & Brennan)" begin
#     res = BraunHolland(dummytabX, dummytabY)
#     @test round(coef(res)[1]; digits = 4) == 1.0331
#     @test round(coef(res)[2]; digits = 4) == -.1927
#     popXonY = ExpandTable(res.table.lYx, Int64.(dummytabX.tabX.freq))
#     @test round(mean(popXonY); digits = 4) == 16.8329
#     @test round(std(popXonY; corrected = false); digits = 4) == 6.6017
#     @test round(skewness(popXonY); digits = 4) == .5799
# end

@testset "FrequencyEstimation" begin
    res = FrequencyEstimation(dummytabX, dummytabY; case = :upper, w₁ = 1.0)
    @test round.(res.table.eYx; digits = 2) == [-0.02, 0.83, 1.76, 2.92, 3.98, 5.00]
end

@testset "ChainedEquipercentile" begin
    res = ChainedEquipercentile(neattabX, neattabY; case = :upper)
    # Largest 3 values in the concordance table are matched.
    @test round.(res.table.eYx[end-2:end]; digits = 4) == [34.3125, 35.4125, 36.0938]
end


text:
# Frames ================================================================================= #
"""
    offsetframe(man::AbstractManifold{D, C}, args::NTuple{D, T})

Return a matrix that makes vectors perpendicular to `man(args...)` when multiplied with
`*ₚ` and the value of `man(args...)`.

Should return a `SVector{D+C}` and a matrix.
"""
offsetframe(man::AbstractManifold{D, C}, t::Vararg{T, D}) where {D, C, T} =
    man(t...), vcat(zeros(T, D+C, 1), one(T))
offsetframe(man::AbstractManifold{0}) =
    man(), zeros(0, 0)

function offsetframe(c::AbstractManifold{1, C}, t::T) where {C, T}
    # This only works for codimension up to 2.
    if C > 2
        return invoke(offsetframe, Tuple{AbstractManifold{1, C}, T}, c, t)
    end

    val = c(Dual{:d1}(Dual{:d2}(t, one(T)), one(T)))
    if length(val) < 3
        val = append!([val...], fill(zero(eltype(val)), 3 - length(val)))
    elseif length(val) > 3
        # This should never happen, except if C wrong.
        throw(ArgumentError("Curve has codimension $(length(val) - 1)"))
    end

    offset  = value.(value.(val))
    tangent = value.(partials.(val, 1))
    normal  = partials.(partials.(val, 1), 1)
    if iszero(normal × tangent)
        normal = !iszero(tangent × @SVector[0, 0, 1]) ?
            tangent × @SVector[0, 0, -1] : tangent × @SVector[0, -1, 0]
    end
    binormal = normalize(tangent × normal)
    normal   = normalize(binormal × tangent)

    offset[1:1+C], [normal binormal]
end

dimincrease(::AbstractManifold{0}) = 0
dimincrease(::AbstractManifold{D, C}) where {D, C} = D+C
dimincrease(::AbstractManifold{1, C}) where C = C ≤ 2 ? 1 : 1+C

"""
    u +ₚ v

Promote vectors `u` and `v` to common size and add them.
"""
+ₚ(u::SVector{N, T}, v::SVector{M, U}) where {N, T, M, U} =
    vcat(u, @SVector zeros(T, max(N, M) - N)) +
    vcat(v, @SVector zeros(U, max(N, M) - M))

function +ₚ(u::AbstractVector{T}, v::AbstractVector{U}) where {T, U}
    n = length(u)
    m = length(v)
    vcat(u, zeros(T, max(n, m) - n)) + vcat(v, zeros(U, max(n, m) - m))
end

"""
    M *ₚ v

Promote matrix `M` and vector/matrix `v` to compatible sizes and multiply them.
TODO: better explanation.

# Examples

    julia> [0 0; 1 0; 0 1] *ₚ [1,1,1]
    4-element Array{Int64, 1}:
     0
     1
     1
     1
"""
function *ₚ(M::AbstractMatrix{T}, v::AbstractVector{U}) where {T, U}
    n = length(v)
    m1, m2 = size(M)
    if m2 ≥ n
        M * vcat(v, zeros(U, m2 - n))
    else
        d = n - m2
        [M zeros(T, m1, d);
         zeros(T, d, m2) I] * v
    end
end

function *ₚ(A::AbstractMatrix{T}, B::AbstractMatrix{U}) where {T, U}
    n1, n2 = size(A)
    m1, m2 = size(B)
    if m1 > n2
        d = m1 - n2
        [A zeros(T, n1, d);
         zeros(T, d, n2) I] * B
    elseif n2 > m1
        d = n2 - m1
        A * [B zeros(U, m1, d);
             zeros(U, d, m2) I]
    else
        A * B
    end
end

# Product ================================================================================ #
struct ProductSpace{D, C, T<:NTuple{2, AbstractManifold}} <: AbstractManifold{D, C}
    spaces::T
end
Base.show(io::IO, ps::ProductSpace) =
    print(io, "(", ps.spaces[1], " × ", ps.spaces[2], ")")

# D, C calculation is wrong.
# TODO: add a trait dimensionincrease??
function LinearAlgebra.cross(m1::AbstractManifold{D1, C1},
                             m2::AbstractManifold{D2, C2}) where {D1,D2,C1,C2}
    D = D1 + D2
    #C = max(D1 + C1, D1 + D2 + C2) - D
    C = dimincrease(m1) - D1 + C2

    spaces = (m1, m2)
    ProductSpace{D, C, typeof(spaces)}((m1, m2))
end
LinearAlgebra.cross(::PointSpace, ::PointSpace) = PointSpace()
LinearAlgebra.cross(m::AbstractManifold, ::PointSpace) = m
LinearAlgebra.cross(::PointSpace, m::AbstractManifold) = m

function offsetframe(ps::ProductSpace{D, C}, args::Vararg{T, D}) where {D, C, T}
    m1, m2 = ps.spaces

    val1, mat1 = offsetframe(m1, args[1:dim(m1)]...)
    val2, mat2 = offsetframe(m2, args[dim(m1) .+ (1:dim(m2))]...)
    val2 *= scaling(m2)(args[dim(m1)]...)

    SVector{D+C, T}((val1 +ₚ mat1 *ₚ val2)[1:D+C]), mat1 *ₚ mat2
end

(ps::ProductSpace{D})(args::Vararg{T, D}) where {T, D} = offsetframe(ps, args...)[1]


text:
function Log_LikeliHood(m::AbstractModel, om, leaves)
    if leaves
        #= if the animal has left we want the probability of this to happen
        at that omission, given that it has not left before.
        Since we assume that the events are independent the product of each likelihood
        should give me the total probability. Using the log we can simply sum them =#
        return log(cdf(m,om) - cdf(m,om-1))
    else
        #= if the animal hasn't left we want the probability of not leaving until om =#
        # this is equivalent log(1-cdf(dist,om)) but it's safe when cdf is near 1
        return log(1-cdf(m, om))
    end
end

function  nll_data(m::AbstractModel,cm::Dict)
    -sum(v*Log_LikeliHood(m, k...) for (k, v) in cm)
end

function Distributions.fit(::Type{T}, b::DataFrames.AbstractDataFrame) where T<:AbstractModel
    cm = countmap(StructArray((Omissions_plus_one = b.Omissions_plus_one, Leave = b.Leave)))
    p = init(T, cm)
    res = optimize(params(p)) do param
        all(param .> 1e-10) || return 1e10
        params(p) .= param
        return nll_data(p, cm)
    end
    params(p) .= Optim.minimizer(res)
    return p
end

###
function simulate(m::AbstractModel,n_samples=1)
    w = [cdf(m,x)-cdf(m,x-1) for x in 1:maximum(m)]
    weights = StatsBase.weights(w)
    [sample(1:maximum(m),weights) for _ in 1:n_samples]
end


text:
using Bio.Seq

reader1 = open(FASTA.Reader, "hg19.fasta")

for i in reader1
    println(FASTA.identifier(i))
    println(FASTA.sequence(i)[1:10])
end

close(reader1)

## ref https://biojulia.net/Bio.jl/latest/man/reading/


text:
const term_c = (Tokens.RPAREN, Tokens.RSQUARE, Tokens.RBRACE, Tokens.END, Tokens.ELSE, Tokens.ELSEIF, Tokens.CATCH, Tokens.FINALLY, Tokens.ENDMARKER)

"""
Continue parsing statements until an element of `closers` is hit (usually
`end`). Statements are grouped in a `Block` EXPR.
"""
function parse_block(ps::ParseState, ret::Vector{EXPR}=EXPR[], closers=(Tokens.END,), docable=false)
    prevpos = position(ps)
    while kindof(ps.nt) ∉ closers # loop until an expected closer is hit
        if kindof(ps.nt) ∈ term_c # error handling if an unexpected closer is hit
            if kindof(ps.nt) === Tokens.ENDMARKER
                break
            elseif kindof(ps.nt) === Tokens.RPAREN
                push!(ret, mErrorToken(ps, INSTANCE(next(ps)), UnexpectedToken))
            elseif kindof(ps.nt) === Tokens.RBRACE
                push!(ret, mErrorToken(ps, INSTANCE(next(ps)), UnexpectedToken))
            elseif kindof(ps.nt) === Tokens.RSQUARE
                push!(ret, mErrorToken(ps, INSTANCE(next(ps)), UnexpectedToken))
            else
                push!(ret, mErrorToken(ps, INSTANCE(next(ps)), UnexpectedToken))
            end
        else
            if docable
                a = parse_doc(ps)
            else
                a = parse_expression(ps)
            end
            push!(ret, a)
        end
        prevpos = loop_check(ps, prevpos)
    end
    return ret
end

"""
Parses an iterator, allowing for the preceding keyword `outer`. Returns an
error expression if an invalid expression is parsed (anything other than
`=`, `in`, `∈`).
"""
function parse_iterator(ps::ParseState, outer=parse_outer(ps))
    arg = @closer ps :range @closer ps :ws parse_expression(ps)
    if !is_range(arg)
        arg = mErrorToken(ps, arg, InvalidIterator)
    end
    if outer !== nothing
        arg.args[1] = setparent!(EXPR(Outer, EXPR[outer, arg.args[1]]), arg)
        arg.fullspan += outer.fullspan
        arg.span = outer.fullspan + arg.span
    end
    return arg
end

function parse_outer(ps)
    if kindof(ps.nt) === Tokens.OUTER && kindof(ps.nws) !== EmptyWS && !Tokens.isoperator(kindof(ps.nnt))
        outer = INSTANCE(next(ps))
    end
end

"""
    parse_iterators(ps::ParseState, allowfilter = false)

Parses a group of iterators e.g. used in a `for` loop or generator. Can allow
for a succeeding `Filter` expression.
"""
function parse_iterators(ps::ParseState, allowfilter=false)
    arg = parse_iterator(ps)
    if iscomma(ps.nt) # we've hit a comma separated list of iterators.
        arg = EXPR(Block, EXPR[arg])
        prevpos = position(ps)
        while iscomma(ps.nt)
            accept_comma(ps, arg)
            nextarg = parse_iterator(ps)
            push!(arg, nextarg)
            prevpos = loop_check(ps, prevpos)
        end
    end

    if allowfilter
        arg = parse_filter(ps, arg)
    end
    return arg
end

"""
parse_filter(ps::ParseState, arg)

Parse a conditional filter following a generator.
"""
function parse_filter(ps::ParseState, arg)
    if kindof(ps.nt) === Tokens.IF # assumes we're inside a generator
        if typof(arg) === Block
            arg = EXPR(Filter, arg.args)
        else
            arg = EXPR(Filter, EXPR[arg])
        end
        push!(arg, mKEYWORD(next(ps)))
        cond = @closer ps :range parse_expression(ps)
        push!(arg, cond)
    end
    return arg
end

"""
    parse_call(ps, ret)

Parses a function call. Expects to start before the opening parentheses and is passed the expression declaring the function name, `ret`.
"""
function parse_call(ps::ParseState, ret::EXPR, ismacro=false)
    if is_minus(ret) || is_not(ret)
        arg = @closer ps :unary @closer ps :inwhere @precedence ps PowerOp parse_expression(ps)
        if istuple(arg)
            pushfirst!(arg.args, ret)
            ret = EXPR(Call, arg.args)
        elseif iswherecall(arg) && istuple(arg.args[1])
            ret = mWhereOpCall(EXPR(Call, EXPR[ret; arg.args[1].args]), arg.args[2], arg.args[3:end])
        else
            ret = mUnaryOpCall(ret, arg)
        end
    elseif is_and(ret) || is_decl(ret) || is_exor(ret)
        arg = @precedence ps 20 parse_expression(ps)
        if is_exor(ret) && istuple(arg) && length(arg) == 3 && is_splat(arg.args[2])
            arg = EXPR(InvisBrackets, arg.args)
        end
        ret = mUnaryOpCall(ret, arg)
    elseif is_issubt(ret) || is_issupt(ret)
        arg = @precedence ps PowerOp parse_expression(ps)
        ret = EXPR(Call, EXPR[ret; arg.args])
    else
        !ismacro && typof(ret) === MacroName && (ismacro = true)
        args = EXPR[ret, mPUNCTUATION(next(ps))]
        @closeparen ps @default ps parse_comma_sep(ps, args, !ismacro)
        accept_rparen(ps, args)
        ret = EXPR(ismacro ? MacroCall : Call, args)
    end
    return ret
end

"""
Parses a comma separated list, optionally allowing for conversion of
assignment (`=`) expressions to `Kw`.
"""
function parse_comma_sep(ps::ParseState, args::Vector{EXPR}, kw=true, block=false, istuple=false)
    prevpos = position(ps)
    @nocloser ps :inwhere @nocloser ps :newline @closer ps :comma while !closer(ps)
        a = parse_expression(ps)
        if kw && _do_kw_convert(ps, a)
            a = _kw_convert(a)
        end
        push!(args, a)
        if iscomma(ps.nt)
            accept_comma(ps, args)
        else# if kindof(ps.ws) == SemiColonWS
            break
        end
        prevpos = loop_check(ps, prevpos)
    end
    if istuple && length(args) > 2
        block = false
    end

    if kindof(ps.ws) == SemiColonWS
        if @nocloser ps :newline @closer ps :comma @nocloser ps :semicolon closer(ps)
            if block && !(length(args) == 1 && ispunctuation(args[1])) && !(typof(last(args)) === UnaryOpCall && is_dddot(last(args).args[2]))
                push!(args, EXPR(Block, EXPR[pop!(args)]))
            elseif kw && kindof(ps.nt) === Tokens.RPAREN
                push!(args, EXPR(Parameters, EXPR[], 0, 0))
            end
        else
            a = @nocloser ps :newline @closer ps :comma @nocloser ps :inwhere parse_expression(ps)
            if block && !(length(args) == 1 && ispunctuation(args[1])) && !is_splat(last(args)) && !(istuple && iscomma(ps.nt))
                args1 = EXPR[pop!(args), a]
                prevpos = position(ps)
                @nocloser ps :newline @closer ps :comma while @nocloser ps :semicolon !closer(ps)
                    a = parse_expression(ps)
                    push!(args1, a)
                    prevpos = loop_check(ps, prevpos)
                end
                body = EXPR(Block, args1)
                push!(args, body)
                args = body
            else
                parse_parameters(ps, args, EXPR[a])
            end
        end
    end
    return
end

"""
    parse_parameters(ps::ParseState, args::Vector{EXPR}, args1::Vector{EXPR} = EXPR[]; usekw = true)

Parses parameter arguments for a function call (e.g. following a semicolon).
"""
function parse_parameters(ps::ParseState, args::Vector{EXPR}, args1::Vector{EXPR}=EXPR[]; usekw=true)
    isfirst = isempty(args1)
    prevpos = position(ps)
    @nocloser ps :inwhere @nocloser ps :newline  @closer ps :comma while !isfirst || (@nocloser ps :semicolon !closer(ps))
        if isfirst
            a = parse_expression(ps)
        else
            a = first(args1)
        end
        if usekw && _do_kw_convert(ps, a)
            a = _kw_convert(a)
        end
        if isfirst
            push!(args1, a)
        else
            pop!(args1)
            push!(args1, a)
        end
        if iscomma(ps.nt)
            accept_comma(ps, args1)
        end
        if kindof(ps.ws) == SemiColonWS
            parse_parameters(ps, args1; usekw=usekw)
        end
        if isfirst
            prevpos = loop_check(ps, prevpos)
        else
            prevpos = position(ps)
        end
        isfirst = true
    end
    if !isempty(args1)
        paras = EXPR(Parameters, args1)
        push!(args, paras)
    end
    return
end

"""
    parse_macrocall(ps)

Parses a macro call. Expects to start on the `@`.
"""
function parse_macrocall(ps::ParseState)
    at = mPUNCTUATION(ps)
    if !isemptyws(ps.ws)
        mname = mErrorToken(ps, INSTANCE(next(ps)), UnexpectedWhiteSpace)
    else
        mname = EXPR(MacroName, EXPR[at, mIDENTIFIER(next(ps))])
    end

    # Handle cases with @ at start of dotted expressions
    if kindof(ps.nt) === Tokens.DOT && isemptyws(ps.ws)
        prevpos = position(ps)
        while kindof(ps.nt) === Tokens.DOT
            op = mOPERATOR(next(ps))
            nextarg = mIDENTIFIER(next(ps))
            mname = mBinaryOpCall(mname, op, EXPR(Quotenode, EXPR[nextarg]))
            prevpos = loop_check(ps, prevpos)
        end
    end

    if iscomma(ps.nt)
        return EXPR(MacroCall, EXPR[mname], mname.fullspan, mname.span)
    elseif isemptyws(ps.ws) && kindof(ps.nt) === Tokens.LPAREN
        return parse_call(ps, mname, true)
    else
        args = EXPR[mname]
        insquare = ps.closer.insquare
        prevpos = position(ps)
        @default ps while !closer(ps)
            if insquare
                a = @closer ps :insquare @closer ps :inmacro @closer ps :ws @closer ps :wsop parse_expression(ps)
            else
                a = @closer ps :inmacro @closer ps :ws @closer ps :wsop parse_expression(ps)
            end
            push!(args, a)
            if insquare && kindof(ps.nt) === Tokens.FOR
                break
            end
            prevpos = loop_check(ps, prevpos)
        end
        return EXPR(MacroCall, args)
    end
end

"""
parse_generator(ps)

Having hit `for` not at the beginning of an expression return a generator.
Comprehensions are parsed as SQUAREs containing a generator.
"""
function parse_generator(ps::ParseState, ret::EXPR)
    kw = mKEYWORD(next(ps))
    ret = EXPR(Generator, EXPR[ret, kw])
    ranges = @closesquare ps parse_iterators(ps, true)

    if typof(ranges) === Block
        append!(ret, ranges)
    else
        push!(ret, ranges)
    end

    if typof(ret.args[1]) === Generator || typof(ret.args[1]) === Flatten
        ret = EXPR(Flatten, EXPR[ret])
    end

    return ret
end

"""
Helper function for parsing import/using statements.
"""
function parse_dot_mod(ps::ParseState, is_colon=false)
    args = EXPR[]

    prevpos = position(ps)
    while kindof(ps.nt) === Tokens.DOT || kindof(ps.nt) === Tokens.DDOT || kindof(ps.nt) === Tokens.DDDOT
        d = mOPERATOR(next(ps))
        trailing_ws = d.fullspan - d.span
        if is_dot(d)
            push!(args, mOPERATOR(1 + trailing_ws, 1, Tokens.DOT, false))
        elseif is_ddot(d)
            push!(args, mOPERATOR(1, 1, Tokens.DOT, false))
            push!(args, mOPERATOR(1 + trailing_ws, 1, Tokens.DOT, false))
        elseif is_dddot(d)
            push!(args, mOPERATOR(1, 1, Tokens.DOT, false))
            push!(args, mOPERATOR(1, 1, Tokens.DOT, false))
            push!(args, mOPERATOR(1 + trailing_ws, 1, Tokens.DOT, false))
        end
        prevpos = loop_check(ps, prevpos)
    end

    # import/export ..
    # TODO: Not clear what this is for?
    # if iscomma(ps.nt) || kindof(ps.ws) == NewLineWS || kindof(ps.nt) === Tokens.ENDMARKER
    #     if length(args) == 2
    #         return EXPR[INSTANCE(ps)]
    #     end
    # end

    prevpos = position(ps)
    while true
        if kindof(ps.nt) === Tokens.AT_SIGN
            at = mPUNCTUATION(next(ps))
            a = INSTANCE(next(ps))
            push!(args, EXPR(MacroName, EXPR[at, a]))
        elseif kindof(ps.nt) === Tokens.LPAREN
            a = EXPR(InvisBrackets, EXPR[mPUNCTUATION(next(ps))])
            push!(a, @closeparen ps parse_expression(ps))
            accept_rparen(ps, a)
            push!(args, a)
        elseif kindof(ps.nt) === Tokens.EX_OR
            a = @closer ps :comma parse_expression(ps)
            push!(args, a)
        elseif !is_colon && isoperator(ps.nt)
            next(ps)
            push!(args, mOPERATOR(ps.nt.startbyte - ps.t.startbyte,  1 + ps.t.endbyte - ps.t.startbyte, kindof(ps.t), false))
        elseif VERSION > v"1.3.0-" && isidentifier(ps.nt) && isemptyws(ps.nws) && (kindof(ps.nnt) === Tokens.STRING || kindof(ps.nnt) === Tokens.TRIPLE_STRING)
            push!(args, EXPR(NONSTDIDENTIFIER, EXPR[INSTANCE(next(ps)), INSTANCE(next(ps))]))
        else
            push!(args, INSTANCE(next(ps)))
        end

        if kindof(ps.nt) === Tokens.DOT
            push!(args, mPUNCTUATION(next(ps)))
        elseif isoperator(ps.nt) && (ps.nt.dotop || kindof(ps.nt) === Tokens.DOT)
            push!(args, mPUNCTUATION(Tokens.DOT, 1, 1))
            ps.nt = RawToken(kindof(ps.nt), ps.nt.startpos, ps.nt.endpos, ps.nt.startbyte + 1, ps.nt.endbyte, ps.nt.token_error, false, ps.nt.suffix)
        else
            break
        end
        prevpos = loop_check(ps, prevpos)
    end
    args
end


text:
#=--------------------------------------
This is the CPU Simululator....RUN THIS FILE
IF YOU NEED TO EDIT ALU FUNCTIONS ONLY USE: "functions.jl".
Functions which are not operations on registers or constants (like jump, return)
are located here.
----------------------------------------=#

module CPU

include("AssemblyParser.jl")
include("ALU_functions.jl")
include("Stack.jl")

import .AssemblyParser
import .functions
import .Stack

#=-----------------------------------------
CPU ON - Innitialize power state
------------------------------------------=#
#println("on") #Test/Debug print
called = false
PC = 1 #"zero" the program counter (counting starts at 1 in Julia)
 #The zero flag is initialized to false by Flags.jl
 #The carry flag is initialized to false by Flags.jl
functions.regbank("A") #Register bank ‘A’ is selected and therefore ‘A’ is the default bank of registers.

#= ------------------------
RESET FUNCTION DECLARATION
----------------------------=#
function reset()
  global PC
  global numJumps
  global jumped
  PC = 1
  numJumps = 0
  jumped = false

  functions.Flags.set("C", false)
  functions.Flags.set("Z", false)
  functions.regbank("A")

  Stack.reset()
end

reset() #= pointer in the program counter stack is reset to ensure that the program is able to execute programs
in which up to 30 nested subroutine calls can be made.=#

#= ------------------------
GLOBAL VARIABLE DECLARATION
----------------------------=#
global numJumps = 0 #The number of recursive jump calls to the same address
global jumped = false #Is true if "jump" is the last instruction executed
labelDict = Dict() #Dictionary which, maps labels to their address in the program memory

#=-----------------------------------------------
FUNCTION DECLARATION:
(ALU functions are located in ALU_functions.jl)
------------------------------------------------=#

function jump(label)
  global jumped
  jumped = true #The last instruction was a jump instruction
  label = chomp(label) #Cleans up garbage
  PC_current = PC
  PC_new = labelDict["$(label)"] #Change the program counter to point to label
  if PC_new == PC_current
    global numJumps+=1 #The last instruction jumped to itself
  end
  global PC = PC_new #Change the program coutner to jump to label
end

function jump(inCondition, label)
  global jumped
  conditions = Dict("Z" => functions.getFlag("Z"), "NZ" => !functions.getFlag("Z")
  , "C" => functions.getFlag("C"), "NC" => !functions.getFlag("C"))

  if conditions[inCondition]
    jumped = true #The last instruction was a jump instruction
    label = chomp(label) #Cleans up garbage
    PC_current = PC
    PC_new = labelDict["$(label)"] #Change the program counter to point to label
    global PC = PC_new #Change the program coutner to jump to label
  else
    #Do nothing, the CPU will increment the PC
  end
end

function jumpAt(sX, sY)
  global jumped
  global PC
  newAddr = ((UInt8(functions.get(sX)) & (0x0F)) * (0x100)) + UInt8(functions.get(sY)) #New address is lower 4 bits of sX:sY
  jumped = true #The last instruction was a jump instruction
  PC_current = PC
#  println("The address is:$(newAddr)")
  PC_new = newAddr #Change the program counter to point to new address
  if PC_new == PC_current
    global numJumps+=1 #The last instruction jumped to itself
  end
  global PC = PC_new #Change the program coutner to jump to label
end

function thisCall(k)
  try
    global PC
    global jumped
    jumped = true
    Stack.push(PC)
    PC = k
  catch # If the Stack overflows then reset
    reset()
  end
end

function thisCall(flag, k)
  addr = labelDict["$k"]
  if flag == "Z" && functions.getFlag("Z")
    thisCall(addr)
  elseif flag == "NZ" && ! functions.getFlag("Z")
    thisCall(addr)
  elseif flag == "C" && functions.getFlag("C")
    thisCall(addr)
  elseif flag == "NC" && ! functions.getFlag("C")
    thisCall(addr)
  end
end

function thisCallAt(sX, sY)
  a = convert(UInt16, sX)
  b = convert(UInt16, sY)

  # take the lower 4 bits of a
  a = a & 0xF
  # take the lower 8 bits of b
  b = b & 0xFF
  # shift a left by 8 bits
  a = a << 8

  addr = a | b

  thisCall(Int64(addr))
end

function thisReturn()
  try
    global jumped = true
    global PC = Stack.pop()
  catch
    reset()
  end
end

function thisReturn(flag)
  if flag == "Z" && functions.getFlag("Z")
    thisReturn()
  elseif flag == "NZ" && ! functions.getFlag("Z")
    thisReturn()
  elseif flag == "C" && functions.getFlag("C")
    thisReturn()
  elseif flag == "NC" && ! functions.getFlag("C")
    thisReturn()
  end
end

function input(sX)
  nIn = parse(UInt8,readline(STDIN)) #Reads in a number from 0 to 255 to write to the target register
  functions.setReg(sX,nIn)
end

function output(sK,pp)
  secondOpRegister = (pp[1] == '(')

  if secondOpRegister
    portNumber = functions.get(pp[2:length(pp)-1])
  else
    portNumber = parse(Int,pp)
  end

  value_of_sK = functions.get(sK)
 println("$(hex(portNumber,2)) $(hex(value_of_sK,2))") #Prints output of the format "portNumber registerValue"
end

function outputk(kk,p)
  portNumber = parse(Int,p)
  value = parse(Int,kk)
  println("$(hex(portNumber,2)) $(hex(value,2))") #Prints output of the format "portNumber Value"
end

function loadAndReturn(sX, kk)
  functions.load(sX, kk)
  thisReturn()
end

#=-------------------------------------------------------------------
These dictionaries map instruction codeword syntax to Julia functions
of ALU_functions.jl and CPU.jl
(This is a co-req of parsing the assembly code)
--------------------------------------------------------------------=#
instructionsOneArgDict = Dict("SL0" => functions.sl0, "SL1" => functions.sl1, "SLX" => functions.slx, "SLA" => functions.sla
, "RL" => functions.rl, "SR0" => functions.sr0, "SR1" => functions.sr1, "SRX" => functions.srx, "SRA" => functions.sra, "RR" => functions.rr, "JUMP" => jump
, "CALL" => thisCall, "STORE" => functions.store, "FETCH" => functions.fetch, "INPUT" => input
, "RETURN" => thisReturn, "REGBANK" => functions.regbank, "HWBUILD" => functions.hwbuild) #Creates dictionary of functions w/ one arguments

instructionsTwoArgDict = Dict("LOAD" => functions.load,"STAR" => functions.star,"AND" => functions.and,"OR" => functions.or,"XOR" => functions.xor,
"ADD" => functions.add,"ADDCY" => functions.addcy,"SUB" => functions.sub,"SUBCY" => functions.subcy,"TEST" => functions.test,"TESTCY" => functions.testcy,
"COMPARE" => functions.compare,"INPUT" => input,"OUTPUT" => output,"OUTPUTK" => outputk,"STORE" => functions.store,
"FETCH" => functions.fetch,"JUMP" => jump, "JUMP@" => jumpAt, "CALL" => thisCall, "CALL@" => thisCallAt, "LOAD&RETURN" => loadAndReturn) #Creates dictionary of functions w/ two arguments

#=------------------------------------------
BEGIN READING THE INPUT FILE OF INSTRUCTIONS
--------------------------------------------=#
#Runs the parser and returns a list of instructions
#for the program memory and a dictionary containing
#all labels and their locations in the program memory
(instructions,labelDict) = AssemblyParser.Parse(ARGS) #Note that this is the function Parse from "AssemblyParser.jl", NOT the built in Julia function "parse()"

#=----------------------------------------
BEGIN EXECUTION OF PARSED ADDEMBLY CODE
-----------------------------------------=#

  while PC <= size(instructions)[1] && numJumps<10 #Checks whether the Program Counter has reached the end of the program memory or the program has jumped to the same label 10 times
    currentInst = instructions[PC,2] #Gets the next instruction from the program memory
  # println("$PC")
    # Is this redundant?
    if isempty(currentInst) #Checks if there are no more instructions to execute
      break #End execution immediately
    end

    firstArg = instructions[PC,3] #Gets the first argument of the current instruction(might be empty string)
    secondArg = instructions[PC,4] #Gets the second argument of the current instruction(might be empty string)

    hasArg1 = !isempty(firstArg) #Fails on first argument being empty
    hasArg2 = !isempty(secondArg) #Fails on the second argument being empty

    if hasArg1 == false && hasArg2 == false #1st and 2nd arguments are empty
      thisReturn() #Return is the only function which can be called with no arguments
    elseif hasArg1 == true && hasArg2 == false #1st argument is non empty, 2nd is empty
      instructionsOneArgDict["$(currentInst)"](firstArg) #Call the function with one argument
    elseif hasArg1 == true && hasArg2 == true #1st argument is non empty, 2nd arument is non empty
      instructionsTwoArgDict["$(currentInst)"](firstArg,secondArg) #Call the function with two arguments
    end

    if !jumped
      PC+=1
    end #Increment the program counter if the last instruction executed is not a jump

    jumped = false #"Jump" was not the last instruction
  end
  #=------------------------------
  END EXECUTION OF ASSEMBLY CODE
  ------------------------------=#
  #=------------
  CPU OFF
  ------------=#

  #println("off") #Test/Debug print

end

max_stars_repo_name: kinoroy/KCPSM-Simulator

max_issues_repo_name: kinoroy/KCPSM-Simulator

max_forks_repo_name: kinoroy/KCPSM-Simulator


text:
module Rewrite

using DataStructures

export xsub, xmatch, XMatch, rewrite, flatten, distribute_over, distribute_over2

import Base: (==)

# package code goes here
function xsub(e::Expr, subs::Dict)
  es = get(subs, e, nothing)
  if es != nothing
    return es
  end
  nh =  xsub(e.head, subs)
  nargs = [xsub(arg, subs) for arg=e.args]
  re = Expr(nh)
  re.args = nargs
  re
end

function xsub(e, subs::Dict)
  get(subs, e, e)
end

type Wildcard
  symbol::Symbol
#   minlength::Int64
#   maxlength::Int64
#   Wildcard(s) = new(symbol(s), 1, 1)
#   Wildcard(s,j,k) = new(symbol(s), int(j), int(k))
end

function matchrange(w::Wildcard)
  Int64[w.minlength, w.maxlength]
end

function matchrange(v::Vector)
  rng = zeros(Int, 2)
  for vj=v
    rng += matchrange(vj)
  end
  rng
end

function matchrange(x)
  Int64[1, 1]
end

type XMatch
  success::Bool
  subs::Dict{Symbol,Any}
  path::Vector{Int64}
  XMatch(success, subs, path) = new(success, subs, path)
  XMatch() = new(false, Dict{Symbol,Any}(),Int64[])
  XMatch(success::Bool) = new(success, Dict{Symbol,Any}(),Int64[])
  XMatch(success::Bool, path::Vector{Int64}) = new(success, Dict{Symbol,Any}(), path)
  XMatch(subs::Dict) = new(true, convert(Dict{Symbol,Any}, subs),Int64[])
  XMatch(subs::Dict, path::Vector{Int64}) = new(true, convert(Dict{Symbol,Any}, subs),path)
end

(==)(x1::XMatch, x2::XMatch) = (x1.success==x2.success) && (x1.subs == x2.subs)

function xmatch(p, e)
  XMatch(p==e)
end

function xmatch(p::Expr, e::Expr)
  mh = xmatch(p.head, e.head)
  if !mh.success
    return mh
  end
  p = xsub(p, mh.subs)
  ma = xmatch_args(p.args, e.args)
  if ma.success
    merge!(ma.subs, mh.subs)
    return ma
  else
    return XMatch()
  end
end

function xmatch(p::Symbol, e)
  ps = string(p)
  if endswith(ps, "_")
    if length(ps) == 1
      return XMatch(true)
    end
    return XMatch([p=>e])
  end
  return XMatch(p==e)
end

# AnyDict = Union(Dict{Any, Any}, OrderedDict{Any, Any})
# SymbolDict = Union(Dict{Symbol, AnyDict}, OrderedDict{Symbol, AnyDict})
# type Rules
#   atomic::AnyDict
#   ByHead::SymbolDict
#   CallByHead::SymbolDict
#   PartialCallByHead::SymbolDict
#   Rules() = new(OrderedDict{Any,AnyDict}(),OrderedDict{Symbol,AnyDict}(),OrderedDict{Symbol,AnyDict}(),OrderedDict{Symbol,AnyDict}())
# end

function xmatch_args(pargs::AbstractVector, eargs::AbstractVector)
  lp = length(pargs)
  le = length(eargs)
  if lp != le
    return XMatch(false)
  end
  m = XMatch(true)
  for k=1:lp
    mk = xmatch(xsub(pargs[k], m.subs), eargs[k])
    if !mk.success
      return mk
    end
    merge!(m.subs, mk.subs)
  end
  return m
end

# function xmatch_partial(callhead::Symbol, rules::Rules, eargs::AbstractVector)
#   prules = rules.PartialCallByHead[callhead]
#   lhs = keys(prules)
#   rhs = values(prules)

# end

# function rewrite(rules::Union(Dict, OrderedDict), e, recursive=true)
#   from = keys(rules)
#   to = values(rules)
#   success = true
#   while success
#     success = false
#     for (fk,tk) = zip(from, to)
#       xm = xmatch(fk, e)
#       if xm.success
#         if isa(tk, Function)
#           e = tk(;xm.subs...)
#         else
#           e = xsub(tk, xm.subs)
#         end
#         success = true
# #         println(e)
#       end
#     end
#   end
#   return e
# end

function rewrite_top(rules, e)
  from = keys(rules)
  to = values(rules)
  for (fk,tk) = zip(from, to)
    xm = xmatch(fk, e)
    if xm.success
      if isa(tk, Function)
        return XRewrite(true, tk(;xm.subs...))
      else
        return XRewrite(true, xsub(tk, xm.subs))
      end
    end
  end
  XRewrite(false, e)
end

type XRewrite
  success::Bool
  result
end

function rewrite_step(rules, e, recursive=true)
  rewrite_top(rules, e)
end

function rewrite_step(rules, e::Expr, recursive=true)
  success = false
  if recursive
    rh = rewrite_step(rules, e.head)
    success = rh.success
    nargs = {}
    for arg in e.args
      rarg = rewrite_step(rules, arg, recursive)
      success = success || rarg.success
      push!(nargs, rarg.result)
    end
    e = Expr(rh.result)
    e.args = nargs
  end
  rt = rewrite_top(rules, e)
  return XRewrite(rt.success || success, rt.result)
end

function xrewrite(rules, e, recursive=true)
#   print(e, " [ ")
  success = true
  anysuccess = false
  while success
    er = rewrite_step(rules, e, recursive)
    success = er.success
    anysuccess = anysuccess || success
    e = er.result
  end
#   println(" ] -> ", e)
  XRewrite(anysuccess, e)
end

function rewrite(rules, e, recursive=true)
  xres = xrewrite(rules, e, recursive)
  xres.result
end

function rewrite_partial(rules, e::Expr, recursive=true)

end

function flatten(e, args...)
  e
end

function operands(op, e::Expr)
  if e.head === :call && first(e.args) === op
    e.args[2:end]
  else
    Any[e]
  end
end

function operands(op, e)
  Any[e]
end

function flatten(e::Expr, flat_ops=[:+,:*], recursive::Bool=true)
  if isa(flat_ops, Symbol)
    flat_ops = [flat_ops]
  end
  fargs = e.args
  if recursive
    fargs = [flatten(arg, flat_ops, recursive) for arg in e.args]
  end
  nargs = {}
  if e.head === :call
    op = first(fargs)
    push!(nargs, op)
    for a in fargs[2:end]
      append!(nargs, operands(op, a))
    end
  else
    nargs = fargs
  end
  re = Expr(e.head)
  re.args = nargs
  re
end

function distribute_over(ex::Expr, operations, over::Symbol=:(+), recursive=true, flatten_expr=true)
  if recursive
      args = [distribute_over(ea, operations, over, recursive, flatten_expr) for ea=ex.args]
  else
      args = copy(ex.args)
  end

  if flatten_expr
    ex = flatten(Expr(ex.head, args...), [operations, over], recursive)
    args = ex.args
  end

  if ex.head === :call && (args[1] in operations)
    factor_args = {}
    dims = Int64[]
    for kk=2:length(args)
      arg = args[kk]
      if (isa(arg, Expr) && arg.head === :call && arg.args[1] === over)
        push!(factor_args, arg.args[2:end])
        push!(dims, length(arg.args)-1)
      else
        push!(factor_args, [arg])
        push!(dims, 1)
      end
    end
    ret = Expr(:call, over)
    for idxs in cartesianiter(dims...)
      part = Expr(:call, args[1])
      for (k, idx) = enumerate(idxs)
        current = factor_args[k][idx]

        if flatten_expr && isa(current, Expr) && current.head === :call && current.args[1] === args[1]

          append!(part.args, current.args[2:end])
        else
          push!(part.args, current)
        end
      end

      push!(ret.args, part)
    end
    ret
  else
    re = Expr(ex.head)
    re.args = args
    re
  end
end

function distribute_over(ex, operations, over::Symbol=:(+), recursive=true, flatten=true)
    ex
end

function cartesianiter(dims)
    @task cartesianmap(((x...)->produce(x)), dims)
end

function cartesianiter(dims::Int...)
    cartesianiter(dims)
end

end # module


text:
"""
    symband(x, y, Δy; bandscale=1.0, kwargs...)
    symband(xy, Δy; bandscale=1.0, kwargs...)

Plot a symmetric band of radius `bandscale * Δy` around `(x, y)`.

## Attributes

$(Makie.ATTRIBUTES)
"""
@recipe(SymBand, xy, Δy) do scene
    Attributes(;
        bandscale=Node(1.0),
        color=theme(scene, :patchcolor),
        colormap=theme(scene, :colormap),
        colorrange=Makie.automatic,
        linecolor=Makie.automatic,
        linestyle=theme(scene, :linestyle),
        linewidth=theme(scene, :linewidth),
        cycle=[:color => :patchcolor],
        inspectable=theme(scene, :inspectable),
    )
end

function Makie.plot!(plot::SymBand)
    @extract plot (xy, Δy)
    lower_upper = lift(xy, Δy, plot.bandscale) do xy, Δy, scale
        return (
            Point2f.(first.(xy), last.(xy) .- abs.(scale .* Δy)),
            Point2f.(first.(xy), last.(xy) .+ abs.(scale .* Δy)),
        )
    end
    lower = lift(first, lower_upper)
    upper = lift(last, lower_upper)

    linecolor = lift(plot.color, plot.linecolor) do color, linecolor
        linecolor === Makie.automatic || return linecolor
        c = to_color(color)
        return Makie.RGB(Makie.red(c), Makie.green(c), Makie.blue(c))
    end

    band!(
        plot,
        lower,
        upper;
        color=plot.color,
        colormap=plot.colormap,
        colorrange=plot.colorrange,
    )
    lines!(
        plot,
        xy;
        color=linecolor,
        linestyle=plot.linestyle,
        linewidth=plot.linewidth,
        colormap=plot.colormap,
        colorrange=plot.colorrange,
    )

    return plot
end

Makie.convert_arguments(::Type{<:SymBand}, x, y, Δy) = (map(Point2f, x, y), Δy)


text:
module DownloadArgo

using Printf, Dates, YAML, NetCDF, NCDatasets, CSV, DataFrames, Dierckx

"""
    GDAC_FTP(b::String,y::Int,m::Int)

Download Argo data files for one regional domain (b), year (y), and
month (m) from the `GDAC` FTP server (`ftp://ftp.ifremer.fr/ifremer/argo`
or, equivalently, `ftp://usgodae.org/pub/outgoing/argo`).

```
b="atlantic"; yy=2009:2009; mm=8:12;
for y=yy, m=mm;
    println("\$b/\$y/\$m"); DownloadArgo.GDAC_FTP(b,y,m)
end
```
"""
function GDAC_FTP(b::String,y::Int,m::Int)
    yy = @sprintf "%04d" y
    mm = @sprintf "%02d" m
    c=`wget --quiet -r ftp://ftp.ifremer.fr/ifremer/argo/geo/"$b"_ocean/$yy/$mm`
    run(c)
end

"""
    mitprof_interp_setup(fil::String)

Get parameters etc from yaml file (`fil`).
"""
function mitprof_interp_setup(fil::String)

    meta=YAML.load(open(fil))

    #1. file list

    d=meta["dirIn"]
    b=meta["subset"]["basin"]
    y=meta["subset"]["year"]

    list0=Array{Array,1}(undef,12)
    for m=1:12
        sd="$b"*"_ocean/$y/"*Printf.@sprintf("%02d/",m)
        tmp=readdir(d*sd)
        list0[m]=[sd*tmp[i] for i=1:length(tmp)]
    end

    nf=sum(length.(list0))
    list1=Array{String,1}(undef,nf)
    f=0
    for m=1:12
        for ff=1:length(list0[m])
            f+=1
            list1[f]=list0[m][ff]
        end
    end

    meta["fileInList"]=list1;

    #2. coordinate

    z_std=meta["depthLevels"]
    if length(z_std)>1
        tmp1=(z_std[2:end]+z_std[1:end-1])/2
        z_top=[z_std[1]-(z_std[2]-z_std[1])/2;tmp1]
        z_bot=[tmp1;z_std[end]+(z_std[end]-z_std[end-1])/2]
    else
        z_top=0.9*z_std
        z_bot=1.1*z_std
    end

    meta["z_std"]=z_std
    meta["z_top"]=z_top
    meta["z_bot"]=z_bot

    #3. various other specs

    meta["inclZ"] = false
    meta["inclT"] = true
    meta["inclS"] = true
    meta["inclU"] = false
    meta["inclV"] = false
    meta["inclPTR"] = false
    meta["inclSSH"] = false
    meta["TPOTfromTINSITU"] = true

    meta["doInterp"] = true
    meta["addGrid"] = true
    meta["outputMore"] = false
    meta["method"] = "interp"
    meta["fillval"] = -9999.0
    meta["buffer_size"] = 10000

    meta["fileOut"]=meta["name"]*"_"*meta["subset"]["basin"]*"_"*string(meta["subset"]["year"])*".nc"
    meta["var_out"]=meta["variables"]

    return meta
end

"""
    GetOneProfile(m)

Get one profile from netcdf file.
"""
function GetOneProfile(ds,m)

    #
    t=ds["JULD"][m]
    ymd=Dates.year(t)*1e4+Dates.month(t)*1e2+Dates.day(t)
    hms=Dates.hour(t)*1e4+Dates.minute(t)*1e2+Dates.second(t)

    lat=ds["LATITUDE"][m]
    lon=ds["LONGITUDE"][m]
    lon < 0.0 ? lon=lon+360.0 : nothing

    direction=ds["DIRECTION"][m]
    direc=0
    direction=='A' ? direc=1 : nothing
    direction=='D' ? direc=2 : nothing

    #
    pnum_txt=ds["PLATFORM_NUMBER"][:,m]
    #pnum_txt=pnum_txt[findall((!ismissing).(pnum_txt))]
    ii=findall(skipmissing(in.(pnum_txt,"0123456789")))
    ~isempty(ii) ? pnum_txt=String(vec(Char.(pnum_txt[ii]))) : pnum_txt="9999"
    pnum=parse(Int,pnum_txt)

    #
    p=ds["PRES_ADJUSTED"][:,m]
    p_QC=ds["PRES_ADJUSTED_QC"][:,m]
    if sum((!ismissing).(p))==0
        p=ds["PRES"][:,m]
        p_QC=ds["PRES_QC"][:,m]
    end

    #set qc to 5 if missing
    p_QC[findall(ismissing.(p))].='5'
    #avoid potential duplicates
    for n=1:length(p)-1
        if ~ismissing(p[n])
            tmp1=( (!ismissing).(p[n+1:end]) ).&( p[n+1:end].==p[n] )
            tmp1=findall(tmp1)
            p[n.+tmp1].=missing
            p_QC[n.+tmp1].='5'
        end
    end

    #position and date
    isBAD=0
    ~in(ds["POSITION_QC"][m],"1258") ? isBAD=1 : nothing
    ~in(ds["JULD_QC"][m],"1258") ? isBAD=1 : nothing

    #pressure
    tmp1=findall( (!in).(p_QC,"1258") )
    if (length(tmp1)<=5)&&(length(tmp1)>0)
        #omit these few bad points but keep the profile
        p[tmp1].=missing
    elseif length(tmp1)>5
        #flag the profile (will be omitted later)
        #but keep the bad points (for potential inspection)
        isBAD=1
    end

    #temperature
    t=ds["TEMP_ADJUSTED"][:,m]
    t_QC=ds["TEMP_ADJUSTED_QC"][:,m]
    t_ERR=ds["TEMP_ADJUSTED_ERROR"][:,m]
    t_ERR[findall( (ismissing).(t_ERR) )].=0.0

    if sum((!ismissing).(t))==0
        t=ds["TEMP"][:,m]
        t_QC=ds["TEMP_QC"][:,m]
    end

    #salinity
    if haskey(ds,"PSAL")
        s=ds["PSAL_ADJUSTED"][:,m]
        s_QC=ds["PSAL_ADJUSTED_QC"][:,m]
        s_ERR=ds["PSAL_ADJUSTED_ERROR"][:,m]
        s_ERR[findall( (ismissing).(t_ERR) )].=0.0
        if sum((!ismissing).(t))==0
            s=ds["PSAL"][:,m]
            s_QC=ds["PSAL_QC"][:,m]
        end
    else
        s=fill(missing,size(t))
        s_QC=Char.(32*ones(size(t_QC)))
    end

    if ismissing(t[1]) #this file does not contain temperature data...
        t=fill(missing,size(t))
        t_ERR=fill(0.0,size(t))
    else #apply QC
        tmp1=findall(skipmissing((!in).(t_QC,"1258")))
        t[tmp1].=missing
    end

    if ismissing(s[1]) #this file does not contain salinity data...
        s=fill(missing,size(s))
        s_ERR=fill(0.0,size(s))
    else #apply QC
        tmp1=findall(skipmissing((!in).(s_QC,"1258")))
        s[tmp1].=missing
    end

    prof=Dict()
    prof["pnum_txt"]=pnum_txt
    prof["ymd"]=convert(Union{Int,Missing},ymd)
    prof["hms"]=convert(Union{Int,Missing},hms)
    prof["lat"]=convert(Union{Float64,Missing},lat)
    prof["lon"]=convert(Union{Float64,Missing},lon)
    prof["direc"]=convert(Union{Int,Missing},direc)
    prof["T"]=convert(Array{Union{Float64,Missing}},t)
    prof["S"]=convert(Array{Union{Float64,Missing}},s)
    prof["p"]=convert(Array{Union{Float64,Missing}},p)
    prof["T_ERR"]=convert(Array{Union{Float64,Missing}},t_ERR)
    prof["S_ERR"]=convert(Array{Union{Float64,Missing}},s_ERR)
    prof["isBAD"]=isBAD
    prof["DATA_MODE"]=ds["DATA_MODE"][m]

    return prof
end

"""
    prof_PtoZ!(prof,meta)

Convert prof["p"] to prof["depth"]
"""
function prof_PtoZ!(prof,meta)
    l=prof["lat"]
    v=meta["var_out"][1]
    prof[v]=similar(prof["p"],Union{Missing,Float64})
    prof[v].=missing
    k=findall((!ismissing).(prof["p"]))
    prof[v][k]=[DownloadArgo.sw_dpth(Float64(prof["p"][kk]),Float64(l)) for kk in k]
end

"""
    prof_TtoΘ!(prof,meta)

Convert prof["T"] to potential temperature
"""
function prof_TtoΘ!(prof,meta)
    T=prof[meta["var_out"][2]]
    P=0.981*1.027*prof[meta["var_out"][1]]
    S=35.0*ones(size(T))
    k=findall( (!ismissing).(T) )
    T[k]=[DownloadArgo.sw_ptmp(Float64(S[kk]),Float64(T[kk]),Float64(P[kk])) for kk in k]
end

"""
    prof_interp!(prof,meta)

Interpolate from prof["depth"] to meta["z_std"]
"""
function prof_interp!(prof,meta)
    for ii=2:length(meta["var_out"])
        v=meta["var_out"][ii]
        v_e=v*"_ERR"

        z_std=meta["z_std"]
        t_std=similar(z_std,Union{Missing,Float64})
        e_std=similar(z_std,Union{Missing,Float64})

        z=prof["depth"]
        t=prof[v]
        do_e=haskey(prof,v_e)
        do_e ? e=prof[v_e] : e=[]

        kk=findall((!ismissing).(z.*t))
        if (meta["doInterp"])&&(length(kk)>1)
            z_in=z[kk]; t_in=t[kk]
            do_e ? e_in=e[kk] : nothing
            k=sort(1:length(kk),by= i -> z_in[i])
            z_in=z_in[k]; t_in=t_in[k]
            do_e ? e_in=e_in[k] : nothing
            #omit values outside observed range:
            D0=minimum(skipmissing(z_in))
            D1=maximum(skipmissing(z_in))
            msk1=findall( (z_std.<D0).|(z_std.>D1) )
            #avoid duplicates:
            msk2=findall( ([false;(z_in[1:end-1]-z_in[2:end]).==0.0]).==true )
            if length(kk)>5
                spl = Spline1D(z_in, t_in)
                t_std[:] = spl(z_std)
                t_std[msk1].=missing
                t_std[msk2].=missing
                if do_e
                    spl = Spline1D(z_in, e_in)
                    e_std[:] = spl(z_std)
                    e_std[msk1].=missing
                    e_std[msk2].=missing
                end
            else
                t_std = []
                e_std = []
            end
            prof[v]=t_std
        end
    end
end

"""
    sw_dpth(P,LAT)

Calculate depth in meters from pressure (P; in decibars) and
latitude (LAT; in °N)

```
d = sw_dpth(100.0,20.0)
```
"""
function sw_dpth(P,LAT)
    # Original author:  Phil Morgan 92-04-06  (morgan@ml.csiro.au)
    # Reference: Unesco 1983. Algorithms for computation of fundamental properties of
    # seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp. Eqn 25, p26

    c1 = 9.72659
    c2 = -2.2512E-5
    c3 = 2.279E-10
    c4 = -1.82E-15
    gam_dash = 2.184e-6

    LAT = abs.(LAT)
    X   = sin.(deg2rad.(LAT))
    X   = X.*X
    bot_line = 9.780318*(1.0+(5.2788E-3+2.36E-5*X).*X) + gam_dash*0.5*P
    top_line = (((c4*P+c3).*P+c2).*P+c1).*P
    DEPTHM   = top_line./bot_line

    return DEPTHM
end

"""
    sw_ptmp(S,T,P,PR)

Calculate potential temperature as per UNESCO 1983 report from salinity (S;
in psu), in situ temperature (T; in °C), and pressure (P; in decibar)
relative to PR (in decibar; 0 by default).

```
ptmp = sw_ptmp(S,T,P,PR=missing)
```
"""
function sw_ptmp(S,T,P,PR=0.0)
# Original author:  Phil Morgan
# Reference: Fofonoff, P. and Millard, R.C. Jr
#    Unesco 1983. Algorithms for computation of fundamental properties of
#    seawater. _Unesco Tech. Pap. in Mar. Sci._, No. 44. Eqn.(31) p.39
#    Bryden, H. 1973.
#    "New Polynomials for thermal expansion, adiabatic temperature gradient
#    and potential temperature of sea water." DEEP-SEA RES., 1973, Vol 20

# theta1
del_P  = PR - P
del_th = del_P.*sw_adtg(S,T,P);
th     = T + 0.5*del_th;
q      = del_th;

# theta2
del_th = del_P.*sw_adtg(S,th,P+0.5*del_P);
th     = th + (1 - 1/sqrt(2))*(del_th - q);
q      = (2-sqrt(2))*del_th + (-2+3/sqrt(2))*q;

# theta3
del_th = del_P.*sw_adtg(S,th,P+0.5*del_P);
th     = th + (1 + 1/sqrt(2))*(del_th - q);
q      = (2 + sqrt(2))*del_th + (-2-3/sqrt(2))*q;

# theta4
del_th = del_P.*sw_adtg(S,th,P+del_P);
PT     = th + (del_th - 2*q)/6;

return PT
end

"""
    sw_adtg(S,T,P)

Calculate adiabatic temperature gradient as per UNESCO 1983 routines from salinity
(S; in psu), in situ temperature (T; in °C), and pressure (P; in decibar)
```
adtg = sw_adtg(S,T,P)
```
"""
function sw_adtg(S,T,P)
    # Original author:  Phil Morgan
    # Reference: Fofonoff, P. and Millard, R.C. Jr
    #    Unesco 1983. Algorithms for computation of fundamental properties of
    #    seawater. _Unesco Tech. Pap. in Mar. Sci._, No. 44. Eqn.(31) p.39
    #    Bryden, H. 1973.
    #    "New Polynomials for thermal expansion, adiabatic temperature gradient
    #    and potential temperature of sea water." DEEP-SEA RES., 1973, Vol 20

    a0 =  3.5803E-5;
    a1 =  8.5258E-6;
    a2 = -6.836E-8;
    a3 =  6.6228E-10;

    b0 =  1.8932E-6;
    b1 = -4.2393E-8;

    c0 =  1.8741E-8;
    c1 = -6.7795E-10;
    c2 =  8.733E-12;
    c3 = -5.4481E-14;

    d0 = -1.1351E-10;
    d1 =  2.7759E-12;

    e0 = -4.6206E-13;
    e1 =  1.8676E-14;
    e2 = -2.1687E-16;

    ADTG =      a0 + (a1 + (a2 + a3.*T).*T).*T + (b0 + b1.*T).*(S-35.0) +
    ( (c0 + (c1 + (c2 + c3.*T).*T).*T) + (d0 + d1.*T).*(S-35.0) ).*P +
    (  e0 + (e1 + e2.*T).*T ).*P.*P

    return ADTG
end

end


text:
module StructsOfArrays

using Compat

export StructOfArrays

immutable StructOfArrays{T,N,U<:Tuple} <: AbstractArray{T,N}
    arrays::U
end

function gather_eltypes(T, visited = Set{Type}())
    (!isleaftype(T) || T.mutable) && throw(ArgumentError("can only create an StructOfArrays of leaf type immutables"))
    isempty(T.types) && throw(ArgumentError("cannot create an StructOfArrays of an empty or bitstype"))
    types = Type[]
    push!(visited, T)
    for S in T.types
        sizeof(S) == 0 && continue
        (S in visited) && throw(ArgumentError("Recursive types are not allowed for SoA conversion"))
        if isempty(S.types)
            push!(types, S)
        else
            append!(types, gather_eltypes(S, copy(visited)))
        end
    end
    types
end

@generated function StructOfArrays{T}(::Type{T}, dims::Integer...)
    N = length(dims)
    types = gather_eltypes(T)
    arrtuple = Tuple{[Array{S,N} for S in types]...}
    :(StructOfArrays{T,$N,$arrtuple}(($([:(Array{$(S)}(dims)) for S in types]...),)))
end
StructOfArrays(T::Type, dims::Tuple{Vararg{Integer}}) = StructOfArrays(T, dims...)

@compat Base.IndexStyle{T<:StructOfArrays}(::Type{T}) = IndexLinear()

@generated function Base.similar{T}(A::StructOfArrays, ::Type{T}, dims::Dims)
    if isbits(T) && length(T.types) > 1
        :(StructOfArrays(T, dims))
    else
        :(Array{T}(dims))
    end
end

Base.convert{T,S,N}(::Type{StructOfArrays{T,N}}, A::AbstractArray{S,N}) =
    copy!(StructOfArrays(T, size(A)), A)
Base.convert{T,S,N}(::Type{StructOfArrays{T}}, A::AbstractArray{S,N}) =
    convert(StructOfArrays{T,N}, A)
Base.convert{T,N}(::Type{StructOfArrays}, A::AbstractArray{T,N}) =
    convert(StructOfArrays{T,N}, A)

Base.size(A::StructOfArrays) = size(A.arrays[1])
Base.size(A::StructOfArrays, d) = size(A.arrays[1], d)

function generate_getindex(T, arraynum)
    members = Expr[]
    for S in T.types
        sizeof(S) == 0 && push!(members, :($(S())))
        if isempty(S.types)
            push!(members, :(A.arrays[$arraynum][i...]))
            arraynum += 1
        else
            member, arraynum = generate_getindex(S, arraynum)
            push!(members, member)
        end
    end
    Expr(:new, T, members...), arraynum
end

@generated function Base.getindex{T}(A::StructOfArrays{T}, i::Integer...)
    strct, _ = generate_getindex(T, 1)
    Expr(:block, Expr(:meta, :inline), strct)
end

function generate_setindex(T, x, arraynum)
    s = gensym()
    exprs = Expr[:($s = $x)]
    for (el,S) in enumerate(T.types)
        sizeof(S) == 0 && push!(members, :($(S())))
        if isempty(S.types)
            push!(exprs, :(A.arrays[$arraynum][i...] = getfield($s, $el)))
            arraynum += 1
        else
            nexprs, arraynum = generate_setindex(S, :(getfield($s, $el)), arraynum)
            append!(exprs, nexprs)
        end
    end
    exprs, arraynum
end

@generated function Base.setindex!{T}(A::StructOfArrays{T}, x, i::Integer...)
    exprs = Expr(:block, generate_setindex(T, :x, 1)[1]...)
    quote
        $(Expr(:meta, :inline))
        v = convert(T, x)
        $exprs
        x
    end
end
end # module


text:
"""
################################################################################
#																			   #
#		Methods for conversion between field types, interpolation, etc.		   #
#																			   #
################################################################################
"""

export unflat, _d2ẽ!, _H2d!, E⃗, E⃗x, E⃗y, E⃗z, H⃗, H⃗x, H⃗y, H⃗z, S⃗, S⃗x, S⃗y, S⃗z
export normE!, Ex_norm, Ey_norm, mn
# export _cross, _cross_x, _cross_y, _cross_z, _sum_cross, _sum_cross_x, _sum_cross_y, _sum_cross_z
# export _outer
 #, slice_inv 	# stuff to get rid of soon

"""
In-place/mutating methods
"""
function _H2d!(d::AbstractArray{Complex{T},N}, Hin::AbstractArray{Complex{T},N},
	m::AbstractArray{T,N}, n::AbstractArray{T,N}, mag::AbstractArray{T},
	𝓕!::FFTW.cFFTWPlan)::AbstractArray{Complex{T},N} where {T<:Real,N}
    kx_tc!(d,Hin,m,n,mag);
    mul!(d.data,𝓕!,d.data);
	return d
end

function _d2ẽ!(e::AbstractArray{Complex{T},N}, d::AbstractArray{Complex{T},N},
	ε⁻¹,m::AbstractArray{T,N}, n::AbstractArray{T,N}, mag::AbstractArray{T},
	𝓕⁻¹!::FFTW.cFFTWPlan)::AbstractArray{Complex{T},N} where {T<:Real,N}
    eid!(e,ε⁻¹,d);
    mul!(e.data,𝓕⁻¹!,e.data);
	return e
end
##
function _H2d!(d::AbstractArray{Complex{T},4}, Hin::AbstractArray{Complex{T},4},
	M̂::HelmholtzMap{3,T})::AbstractArray{Complex{T},4} where T<:Real
    kx_tc!(d,Hin,M̂.m,M̂.n,M̂.mag);
    mul!(d.data,M̂.𝓕!,d.data);
	return d
end

function _d2ẽ!(e::AbstractArray{Complex{T},4}, d::AbstractArray{Complex{T},4},
	M̂::HelmholtzMap{3,T})::AbstractArray{Complex{T},4} where T<:Real
    eid!(e,M̂.ε⁻¹,d);
    mul!(e.data,M̂.𝓕⁻¹!,e.data);
	return e
end

function _H2d!(d::AbstractArray{Complex{T},3}, Hin::AbstractArray{Complex{T},3},
	M̂::HelmholtzMap{2,T})::AbstractArray{Complex{T},3} where T<:Real
    kx_tc!(d,Hin,M̂.m,M̂.n,M̂.mag);
    mul!(d.data,M̂.𝓕!,d.data);
	return d
end

function _d2ẽ!(e::AbstractArray{Complex{T},3}, d::AbstractArray{Complex{T},3},
	M̂::HelmholtzMap{2,T})::AbstractArray{Complex{T},3} where T<:Real
    eid!(e,M̂.ε⁻¹,d);
    mul!(e.data,M̂.𝓕⁻¹!,e.data);
	return e
end

##
function _H2d!(d::AbstractArray{Complex{T},4}, Hin::AbstractArray{Complex{T},4},
	ms::ModeSolver{3,T})::AbstractArray{Complex{T},4} where T<:Real
    kx_tc!(d,Hin,ms.M̂.m,ms.M̂.n,ms.M̂.mag);
    mul!(d.data,ms.M̂.𝓕!,d.data);
	return d
end

function _d2ẽ!(e::AbstractArray{Complex{T},4}, d::AbstractArray{Complex{T},4},
	ms::ModeSolver{3,T})::AbstractArray{Complex{T},4} where T<:Real
    eid!(e,ms.M̂.ε⁻¹,d);
    mul!(e.data,ms.M̂.𝓕⁻¹!,e.data);
	return e
end

function _H2d!(d::AbstractArray{Complex{T},3}, Hin::AbstractArray{Complex{T},3},
	ms::ModeSolver{2,T})::AbstractArray{Complex{T},3} where T<:Real
    kx_tc!(d,Hin,ms.M̂.m,ms.M̂.n,ms.M̂.mag);
    mul!(d.data,ms.M̂.𝓕!,d.data);
	return d
end

function _d2ẽ!(e::AbstractArray{Complex{T},3}, d::AbstractArray{Complex{T},3},
	ms::ModeSolver{2,T})::AbstractArray{Complex{T},3} where T<:Real
    eid!(e,ms.M̂.ε⁻¹,d);
    mul!(e.data,ms.M̂.𝓕⁻¹!,e.data);
	return e
end

@inline function flat(f::AbstractArray{SVector{3,T}}) where T
	reinterpret(reshape,T,f)
end

@inline function flat(f::AbstractArray{T}) where {T<:Number}
	return f #reinterpret(reshape,T,f)
end

@inline function unflat(f,nvec::Int,Ns)
	reshape(f,(nvec,Ns...))
end

@inline function unflat(f; ms::ModeSolver)
	Ns = size(ms.grid)
	nev = size(f,2)
	ratio = length(f) // ( nev * N(ms.grid) ) |> Int # (length of vector) / (number of grid points)
	[reshape(f[:,i],(ratio,Ns...)) for i=1:nev]
end

@inline function unflat(f,grid::Grid)
	Ns = size(grid)
	# nev = 1 # size(f,2)
	# ratio = length(f) // ( nev * N(grid) ) |> Int # (length of vector) / (number of grid points)
	# [reshape(f[:,i],(ratio,Ns...)) for i=1:nev]
	ratio = length(f) //  N(grid) |> Int # (length of vector) / (number of grid points)
	reshape(f,(ratio,Ns...))
end

mn(ms::ModeSolver) = ms.M̂.mn # vcat(reshape(ms.M̂.m,(1,3,size(ms.grid)...)),reshape(ms.M̂.n,(1,3,size(ms.grid)...)))

function H⃗(ms::ModeSolver{ND,T}; svecs=true) where {ND,T<:Real}
	Harr = [ fft( tc(unflat(ms.H⃗;ms)[eigind],mn(ms)), (2:1+ND) ) for eigind=1:size(ms.H⃗,2) ]#.* ms.M̂.Ninv
	# svecs ? [ reinterpret(reshape, SVector{3,Complex{T}},  Harr[eigind]) for eigind=1:size(ms.H⃗,2) ] : Harr
	return Harr
end

function H⃗(k,H⃗::AbstractArray{Complex{T}},ω::T,geom::Geometry,grid::Grid{ND}; svecs=true, normalized=true) where {ND,T<:Real}
	Ns = size(grid)
	mag,m⃗,n⃗ = mag_m_n(k,grid)
	Hₜ = reshape(H⃗,(2,Ns...))
	mns = vcat(reshape(flat(m⃗),1,3,Ns...),reshape(flat(n⃗),1,3,Ns...))
	# ε⁻¹ = εₛ⁻¹(ω,geom,grid)
	ε⁻¹, nng⁻¹, ngvd⁻¹ = εₛ⁻¹_nngₛ⁻¹_ngvdₛ⁻¹(ω,geom,grid)
	H = fft( tc( Hₜ,mns ), (2:1+ND) )
	# if normalized
	# 	imagmax = argmax(abs2.(E0))
	# 	E1 = E0 / E0[imagmax]
	# 	E1s = reinterpret(reshape, SVector{3,Complex{T}},  E1)
	# 	E1norm = sum(dot.(E1s, inv.(nng⁻¹) .* E1s )) * δ(grid)
	# 	E = E1 / sqrt(E1norm)
	# else
	# 	E = E0
	# end
	# svecs ?  reinterpret(reshape, SVector{3,Complex{T}},  H) : H
	return H
end

function H⃗(k,Hv::AbstractArray{Complex{T}},ω::T,ε⁻¹,nng,grid::Grid{ND}; normalized=true, nnginv=false)::Array{Complex{T},ND+1} where {ND,T<:Real}
	Ns = size(grid)
	mag,m⃗,n⃗ = mag_m_n(k,grid)
	Hₜ = reshape(Hv,(2,Ns...))
	mns = vcat(reshape(flat(m⃗),1,3,Ns...),reshape(flat(n⃗),1,3,Ns...))
	# H0 = fft( tc( Hₜ,mns ), (2:1+ND) )
	if normalized
		E0 = 1im * ε⁻¹_dot( fft( kx_tc( Hₜ,mns,mag), (2:1+ND) ), ε⁻¹)
		imagmax = argmax(abs2.(E0))
		E0magmax = copy(E0[imagmax])
		E1 = E0 / E0[imagmax]
		# E1s = reinterpret(reshape, SVector{3,Complex{T}},  E1)
		if nnginv
			E1norm = dot(E1,_dot(slice_inv(copy(nng.data)),E1)) * δ(grid)
		else
			E1norm = dot(E1,_dot(nng,E1)) * δ(grid)
		end
		E = E1 / sqrt(E1norm)
		# fieldnorm = E0magmax * sqrt(E1norm)
		H1 = fft(tc(kx_ct( ifft( E, (2:1+ND) ), mns,mag), mns),(2:1+ND) ) / ω
		# H1 = -1.0 * fft( tc( Hₜ,mns ), (2:1+ND) ) / E0magmax / sqrt(E1norm) / ω
	else
		# E = E0
		# fieldnorm = 1.0
		H1 = fft( tc( Hₜ,mns ), (2:1+ND) )
	end

	return H1
end

H⃗x(ms::ModeSolver) = H⃗(ms;svecs=false)[1,eachindex(ms.grid)...]
H⃗y(ms::ModeSolver) = H⃗(ms;svecs=false)[2,eachindex(ms.grid)...]
H⃗z(ms::ModeSolver) = H⃗(ms;svecs=false)[3,eachindex(ms.grid)...]

function E⃗(ms::ModeSolver{ND,T}; svecs=true) where {ND,T<:Real}
	Earr = [ 1im * ε⁻¹_dot( fft( kx_tc( unflat(ms.H⃗; ms)[eigind],mn(ms),ms.M̂.mag), (2:1+ND) ), copy(flat( ms.M̂.ε⁻¹ ))) for eigind=1:size(ms.H⃗,2) ]
	# svecs ? [ reinterpret(reshape, SVector{3,Complex{T}},  Earr[eigind]) for eigind=1:size(ms.H⃗,2) ] : Earr
	return Earr
end

function E⃗(ms::ModeSolver{ND,T}, eigind::Int; svecs=true) where {ND,T<:Real}
	E = 1im * ε⁻¹_dot( fft( kx_tc( unflat(ms.H⃗; ms)[eigind],mn(ms),ms.M̂.mag), (2:1+ND) ), copy(flat( ms.M̂.ε⁻¹ )))
	# svecs ?  reinterpret(reshape, SVector{3,Complex{T}},  E) : E
	return E
end

function E⃗(k,H⃗::AbstractArray{Complex{T}},ω::T,geom::Geometry,grid::Grid{ND}; normalized=true)::Array{Complex{T},ND+1} where {ND,T<:Real}
	Ns = size(grid)
	mag,m⃗,n⃗ = mag_m_n(k,grid)
	H = reshape(H⃗,(2,Ns...))
	mns = vcat(reshape(flat(m⃗),1,3,Ns...),reshape(flat(n⃗),1,3,Ns...))
	# ε⁻¹ = εₛ⁻¹(ω,geom,grid)
	ε⁻¹, nng⁻¹, ngvd⁻¹ = εₛ⁻¹_nngₛ⁻¹_ngvdₛ⁻¹(ω,geom,grid)
	E0 = 1im * ε⁻¹_dot( fft( kx_tc( H,mns,mag), (2:1+ND) ), flat(ε⁻¹))
	if normalized
		imagmax = argmax(abs2.(E0))
		E1 = E0 / E0[imagmax]
		E1s = reinterpret(reshape, SVector{3,Complex{T}},  E1)
		E1norm = sum(dot.(E1s, inv.(nng⁻¹) .* E1s )) * δ(grid)
		E = E1 / sqrt(E1norm)
	else
		E = E0
	end
	#svecs ?  reinterpret(reshape, SVector{3,Complex{T}},  E) : E
	return E
end

function E⃗(k,H⃗::AbstractArray{Complex{T}},ω::T,ε⁻¹,nng,grid::Grid{ND}; normalized=true, nnginv=false)::Array{Complex{T},ND+1} where {ND,T<:Real}
	Ns = size(grid)
	mag,m⃗,n⃗ = mag_m_n(k,grid)
	H = reshape(H⃗,(2,Ns...))
	mns = vcat(reshape(flat(m⃗),1,3,Ns...),reshape(flat(n⃗),1,3,Ns...))
	# ε⁻¹ = εₛ⁻¹(ω,geom,grid)
	# ε⁻¹, nng⁻¹, ngvd⁻¹ = εₛ⁻¹_nngₛ⁻¹_ngvdₛ⁻¹(ω,geom,grid)
	E0 = 1im * ε⁻¹_dot( fft( kx_tc( H,mns,mag), (2:1+ND) ), ε⁻¹)
	if normalized
		imagmax = argmax(abs2.(E0))
		E1 = E0 / E0[imagmax]
		# E1s = reinterpret(reshape, SVector{3,Complex{T}},  E1)
		if nnginv
			E1norm = dot(E1,_dot(slice_inv(copy(nng.data)),E1)) * δ(grid)
		else
			E1norm = dot(E1,_dot(nng,E1)) * δ(grid)
		end
		E = E1 / sqrt(E1norm)
	else
		E = E0
	end
	#svecs ?  reinterpret(reshape, SVector{3,Complex{T}},  E) : E
	return E
end

E⃗x(ms::ModeSolver) = [ E[1,eachindex(ms.grid)...] for E in E⃗(ms;svecs=false) ]
E⃗y(ms::ModeSolver) = [ E[2,eachindex(ms.grid)...] for E in E⃗(ms;svecs=false) ]
E⃗z(ms::ModeSolver) = [ E[3,eachindex(ms.grid)...] for E in E⃗(ms;svecs=false) ]

import Base: abs2
abs2(v::SVector) = real(dot(v,v))

function S⃗(ms::ModeSolver{ND,T}; svecs=true) where {ND,T<:Real}
	# Ssvs = real.( cross.( conj.(E⃗(ms)), H⃗(ms) ) )
	Ssvs = map((E,H)->real.( cross.( conj.(E), H) ), E⃗(ms), H⃗(ms) )
	svecs ? Ssvs : [ reshape( reinterpret( Complex{T},  Ssvs[i]), (3,size(ms.grid)...)) for i in length(Ssvs) ]
end

S⃗x(ms::ModeSolver) = map((E,H)->real( getindex.( cross.( conj.(E), H), 1)), E⃗(ms), H⃗(ms) ) #real.( getindex.( cross.( conj.(E⃗(ms)), H⃗(ms) ), 1) )
S⃗y(ms::ModeSolver) = map((E,H)->real( getindex.( cross.( conj.(E), H), 2)), E⃗(ms), H⃗(ms) ) #real.( getindex.( cross.( conj.(E⃗(ms)), H⃗(ms) ), 2) )
S⃗z(ms::ModeSolver) = map((E,H)->real( getindex.( cross.( conj.(E), H), 3)), E⃗(ms), H⃗(ms) ) #real.( getindex.( cross.( conj.(E⃗(ms)), H⃗(ms) ), 3) )

function normE!(ms)
	E = E⃗(ms;svecs=false)
	Eperp = view(E,1:2,eachindex(ms.grid)...)
	imagmax = argmax(abs2.(Eperp))
	# Enorm = E / Eperp(imagmax)
	ms.H⃗ *= inv(Eperp[imagmax])
end

function normE(E)
	Eperp = view(E,1:2,eachindex(ms.grid)...)
	imagmax = argmax(abs2.(Eperp))
	# Enorm = E / Eperp(imagmax)
	ms.H⃗ *= inv(Eperp[imagmax])
end

function Ex_norm(ms)
	E = E⃗(ms;svecs=false)
	Eperp = view(E,1:2,eachindex(ms.grid)...)
	imagmax = argmax(abs2.(Eperp))
	# Enorm = E / Eperp(imagmax)
	view(E,1,eachindex(ms.grid)...) * inv(Eperp[imagmax])
end

function Ey_norm(ms)
	E = E⃗(ms;svecs=false)
	Eperp = view(E,1:2,eachindex(ms.grid)...)
	imagmax = argmax(abs2.(Eperp))
	# Enorm = E / Eperp(imagmax)
	view(E,2,eachindex(ms.grid)...) * inv(Eperp[imagmax])
end

# function ngdisp1(ω,H,geom,mag,m,n)
# 	ω = sqrt(real(ms.ω²[eigind]))
# 	ω / HMₖH(ms.H⃗[:,eigind],nngₛ(ω,geom;ms),ms.M̂.mag,ms.M̂.m,ms.M̂.n)
# end
#
# function ngdisp1(ms,eigind)
# 	ω = sqrt(real(ms.ω²[eigind]))
# 	ω / HMₖH(ms.H⃗[:,eigind],nngₛ(ω,geom;ms),ms.M̂.mag,ms.M̂.m,ms.M̂.n)
# end
# ( sum( dot.( ( inv.(ms.M̂.ε⁻¹) .* E⃗(ms)[1] ), E⃗(ms)[1]))*δ(ms.grid) ) ./ ( sum.(S⃗z(ms))*δ(ms.grid) )

# sum( real.( getindex.( cross.( conj.(E⃗(ms)), H⃗(ms) ), 3) ) )

# ∫(intgnd;ms::ModeSolver) = sum(ingnd)*δ(ms.grid)

# function 𝓘(ms::ModeSolver{ND,T}) where {ND,T<:Real}
#
# end

"""
################################################################################
#																			   #
#							   Plotting methods					   			   #
#																			   #
################################################################################
"""

using Makie: heatmap
export plot_field, plot_field!

function plot_field!(pos,F,grid;cmap=:diverging_bkr_55_10_c35_n256,label_base=["x","y"],label="E",xlim=nothing,ylim=nothing)
	xs = x(grid)
	ys = y(grid)
	xlim = isnothing(xlim) ? Tuple(extrema(xs)) : xlim
	ylim = isnothing(ylim) ? Tuple(extrema(ys)) : ylim

		# ax = [Axis(pos[1,j]) for j=1:2]
	# ax = [Axis(pos[j]) for j=1:2]
	labels = label.*label_base
	Fs = [view(F,j,:,:) for j=1:2]
	magmax = maximum(abs,F)
	ax1 = pos[1]
	# hms = [heatmap!(pos[j], xs, ys, real(Fs[j]),colormap=cmap,label=labels[j],colorrange=(-magmax,magmax)) for j=1:2]
	hm1 = heatmap!(ax1, xs, ys, real(Fs[1]),colormap=cmap,label=labels[1],colorrange=(-magmax,magmax))
	ax2 = pos[2]
	hm1 = heatmap!(ax2, xs, ys, real(Fs[2]),colormap=cmap,label=labels[2],colorrange=(-magmax,magmax))
	hms = [hm1,hm2]
	# cbar = Colorbar(pos[1,3], heatmaps[2],  width=20 )
	# wfs_E = [wireframe!(ax_E[j], xs, ys, Es[j], colormap=cmap_E,linewidth=0.02,color=:white) for j=1:2]
	# map( (axx,ll)->text!(axx,ll,position=(-1.4,1.1),textsize=0.7,color=:white), ax, labels )
	# hideydecorations!.(ax[2])
	# [ax[1].ylabel= "y [μm]" for axx in ax[1:1]]
	# for axx in ax
	# 	axx.xlabel= "x [μm]"
	# 	xlims!(axx,xlim)
	# 	ylims!(axx,ylim)
	# 	axx.aspect=DataAspect()
	# end
	# linkaxes!(ax...)
	return hms
end

function plot_field(F,grid;cmap=:diverging_bkr_55_10_c35_n256,label_base=["x","y"],label="E",xlim=nothing,ylim=nothing)
	fig=Figure()
	hms = plot_field!(fig[1,1],F,grid;cmap,label_base,label,xlim,ylim)
	fig
end


text:
"""
Make prediction for the next `t` times. Assumes that `y_past` is already ordered.

"""
@traitfn function predict_ar(
    m::TGP, p::Int=3, n::Int=1; y_past=get_y(m)
) where {T,TGP<:AbstractGP{T};!IsMultiOutput{TGP}}
    @assert y_past isa AbstractVector{<:Real}
    @assert length(y_past) >= p
    Xtest = reshape(y_past[(end - p + 1):end], 1, :)
    y_new = zeros(T, n)
    for i in 1:n
        y_new[i] = first(first(first(_predict_f(m, Xtest; covf=false))))
        Xtest = hcat(Xtest[:, 2:end], y_new[i])
    end
    return y_new
end

ks = []

@traitfn function predict_ar(
    m::TGP, p::Int=3, n::Int=1; y_past=get_y(m)
) where {T,TGP<:AbstractGP{T};IsMultiOutput{TGP}}
    @assert y_past isa AbstractVector{<:AbstractVector}
    @assert length(y_past) == m.nTask
    @assert all(length.(y_past) .>= p)
    @assert all(m.nDim .== p)
    Xtest = [reshape(y[(end - p + 1):end], 1, :) for y in y_past]
    y_new = [zeros(T, n) for _ in 1:(m.nTask)]
    for i in 1:n
        setindex!.(y_new, first.(first.(first(_predict_f(m, Xtest; covf=false)))), i)
        Xtest = [hcat(Xtest[j][:, 2:end], y_new[j][i]) for j in 1:(m.nTask)]
    end
    return y_new
end

@traitfn function sample_ar(
    m::TGP, p::Int, n::Int=1; y_past=get_y(m)
) where {T,TGP<:AbstractGP{T};!IsMultiOutput{TGP}}
    @assert y_past isa AbstractVector{<:Real}
    @assert length(y_past) >= p
    Xtest = reshape(y_past[(end - p):end], 1, :)
    y_new = zeros(T, n)
    for i in 1:n
        μ, σ² = _predict_f(m, Xtest; covf=true)
        y_new[i] = rand(Normal(first(μ), sqrt(first(σ²))))
        Xtest = hcat(Xtest[:, 2:end], y_new[i])
    end
    return y_new
end

@traitfn function sample_ar(
    m::TGP, p::Int, n::Int=1; y_past=get_y(m)
) where {T,TGP<:AbstractGP{T};IsMultiOutput{TGP}}
    @assert y_past isa AbstractVector{<:AbstractVector}
    @assert length(y_past) == m.nTask
    @assert all(length.(y_past) .>= p)
    Xtest = [reshape(y[(end - p + 1):end], 1, :) for y in y_past]
    y_new = [zeros(T, n) for _ in 1:(m.nTask)]
    for i in 1:n
        μ, σ² = _predict_f(m, Xtest; covf=true)
        μ = first.(first.(μ))
        σ² = first.(first.(σ²))
        setindex!.(y_new, rand.(Normal.(μ, sqrt.(σ²))), i)
        Xtest = [hcat(Xtest[j][:, 2:end], y_new[j][i]) for j in 1:(m.nTask)]
    end
    return y_new
end


text:
# parametric(like generics) composite
struct Point{T}
   x::T
   y::T
end

Point{Float64}(0.1, 33.3)
Point{String}("0.0", "-10")
# type error
# Point{Int64}(10, 10.1)

# parametric type is not covariant, and not contravariant
# https://docs.julialang.org/en/v1/manual/types/#man-parametric-composite-types
println((Float64 <: Real, Point{Float64} <: Point{Real})) # -> true, false

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
using OVERTVerify
using LazySets
using Dates
using JLD2
#ENV["JULIA_DEBUG"] = Main

controller_name = ARGS[1]
controller = "nnet_files/jmlr/car_"*controller_name*"_controller.nnet"
println("Controller is: ", controller)

query = OvertQuery(
	SimpleCar,  # problem
	controller,    # network file
	Id(),      	# last layer activation layer Id()=linear
	"MIP",     	# query solver, "MIP" or "ReluPlex"
	10,        	# ntime
	0.2,       	# dt
	-1,        	# N_overt
	)

input_set = Hyperrectangle(low=[9.5, -4.5, 2.1, 1.5], high=[9.55, -4.45, 2.11, 1.51])
concretization_intervals = [5, 5]
t1 = Dates.time()
concrete_state_sets, symbolic_state_sets, concrete_meas_sets, symbolic_meas_sets = symbolic_reachability_with_concretization(query, input_set, concretization_intervals)
t2 = Dates.time()
dt = (t2-t1)
print("elapsed time= $(dt) seconds")

# clean up sets
init_set0, reachable_state_sets = clean_up_sets(concrete_state_sets, symbolic_state_sets, concretization_intervals, dims=[1,2])

# In this example, our property is the following:
# We want the car to reach the box [-.6, .6] [-.2,.2]
# at SOME point in the time history
# constraint 1  x1 >= -0.6  --> -x1 <= 0.6
c1 = HalfSpace([-1.0, 0.0, 0.0, 0.0], 0.6)
# constraint 2 x1 <= 0.6
c2 = HalfSpace([1.0, 0.0, 0.0, 0.0], 0.6)
# constraint 3  x2 >= -0.2   -> -x2 <= 0.2
c3 = HalfSpace([-1.0, 0.0, 0.0, 0.0], 0.2)
# constraint 4  x2 <= 0.2
c4 = HalfSpace([1.0, 0.0, 0.0, 0.0], 0.2)
goal_set = HPolyhedron([c1, c2, c3, c4])
t1 = time()
goal_reached_steps = reachable_state_sets .⊆ Ref(goal_set)
goal_reached = any(goal_reached_steps)
dt_check = time() - t1
println("Goal reached: $goal_reached")
println("Goal reached at step: $(findfirst(goal_reached_steps))")

JLD2.@save "src/examples/jmlr/data/car_reachability_"*string(controller_name)*"_controller_data.jld2" query input_set concretization_intervals goal_set concrete_state_sets symbolic_state_sets concrete_meas_sets symbolic_meas_sets reachable_state_sets dt goal_reached goal_reached_steps dt_check