text:
export DAG, NodeType, Leaf, Inner,
       children, has_children, num_children, isleaf, isinner,
       foreach, foreach_down, filter, foldup, foldup_aggregate,
       num_nodes, num_edges, tree_num_nodes, tree_num_edges, in,
       innernodes, leafnodes, num_innernodes, num_leafnodes, linearize,
       left_most_descendent, right_most_descendent,
       num_parents, label_nodes, feedforward_layers,
       node_stats, innernode_stats, leaf_stats,
       parent_stats

#####################
# types and traits
#####################

"A node in a directed acyclic graph (of which it is the root)"
abstract type DAG end

"""
A trait hierarchy denoting types of `DAG` nodes
`NodeType` defines an orthogonal type hierarchy of node types, so we can dispatch on node type regardless of the graph type.
See @ref{https://docs.julialang.org/en/v1/manual/methods/#Trait-based-dispatch-1}
"""
abstract type NodeType end

"The trait of leaf nodes (nodes without children)"
struct Leaf <: NodeType end

"The trait of inner nodes (nodes that have children)"
struct Inner <: NodeType end

#####################
# basic fields and methods
#####################

# Each `DAG` is required to provide a specialized method for the following functions.

"Get the node type trait of the given `Node`"
@inline NodeType(node::DAG) = NodeType(typeof(node))

"Get the children of a given inner node"
function children end

#####################
# derived node functions
#####################

"Does the DAG node have children?"
@inline has_children(n::DAG)::Bool = has_children(NodeType(n), n)
@inline has_children(::Inner, n::DAG)::Bool = true
@inline has_children(::Leaf, n::DAG)::Bool = false

"Get the number of children of a given inner DAG node"
@inline num_children(n::DAG)::Int = num_children(NodeType(n), n)
@inline num_children(::Inner, n::DAG)::Int = length(children(n))
@inline num_children(::Leaf, n::DAG)::Int = 0

"Is the DAG node a leaf node?"
@inline isleaf(n::DAG) = NodeType(n) isa Leaf

"Is the DAG node an inner node?"
@inline isinner(n::DAG) = NodeType(n) isa Inner

#####################
# traversal
#####################

import Base.foreach #extend

"Apply a function to each node in a graph, bottom up"
foreach(f::Function, node::DAG, ::Nothing=nothing) =
    foreach(f, node, Dict{DAG,Nothing}())

function foreach(f::Function, node::DAG, seen)
    get!(seen, node) do
        if isinner(node)
            for c in children(node)
                foreach(f, c, seen)
            end
        end
        f(node)
        nothing
    end
    nothing
end

function foreach(node::DAG, f_leaf::Function, f_inner::Function, seen=nothing)
    foreach(node, seen) do n
        isinner(n) ? f_inner(n) : f_leaf(n)
    end
    nothing
end

"Apply a function to each node in a graph, top down"
function foreach_down(f::Function, node::DAG)
    # naive implementation
    lin = linearize(node)
    foreach(f, Iterators.reverse(lin))
end

import Base.filter #extend

"""Retrieve list of nodes in DAG matching predicate `p`"""
function filter(p::Function, root::DAG, seen=nothing, ::Type{T} = Union{})::Vector where T
    results = Vector{T}()
    foreach(root, seen) do n
        if p(n)
            if !(n isa eltype(results))
                results = collect(typejoin(eltype(results), typeof(n)), results)
            end
            push!(results, n)
        end
    end
    results
end

"""
    foldup(node::DAG,
        f_leaf::Function,
        f_inner::Function,
        ::Type{T})::T where {T}

Compute a function bottom-up on the graph.
`f_leaf` is called on leaf nodes, and `f_inner` is called on inner nodes.
Values of type `T` are passed up the circuit and given to `f_inner` as a function on the children.
"""
function foldup(node::DAG, f_leaf::Function, f_inner::Function, ::Type{T}, ::Nothing=nothing) where {T}
    foldup(node, f_leaf, f_inner, T, Dict{DAG,T}())
end

function foldup(node::DAG, f_leaf::Function, f_inner::Function, ::Type{T}, cache) where {T}
    get!(cache, node) do
        if isinner(node)
            callback(c) = foldup(c, f_leaf, f_inner, T, cache)::T
            f_inner(node, callback)::T
        else
            f_leaf(node)::T
        end
    end
end

"""
Compute a function bottom-up on the circuit.
`f_leaf` is called on leaf nodes, and `f_inner` is called on inner nodes.
Values of type `T` are passed up the circuit and given to `f_inner` in aggregate
as a vector from the children.
"""
function foldup_aggregate(node::DAG, f_leaf::Function, f_inner::Function, ::Type{T}, ::Nothing=nothing) where {T}
    foldup_aggregate(node, f_leaf, f_inner, T, Dict{DAG,T}())
end

function foldup_aggregate(node::DAG, f_leaf::Function, f_inner::Function, ::Type{T}, cache) where {T}
    get!(cache, node) do
        if isinner(node)
            child_values = Vector{T}(undef, num_children(node))
            map!(c -> foldup_aggregate(c, f_leaf, f_inner, T, cache)::T,
                                        child_values, children(node))
            f_inner(node, child_values)::T
        else
            f_leaf(node)::T
        end
    end
end

#####################
# methods using circuit traversal
#####################

import Base.iterate #extend

function iterate(node::DAG, state=nothing)
    if isnothing(state)
        lin = linearize(node)
        liniter = iterate(lin)
    else
        lin = first(state)
        liniter = iterate(lin, last(state))
    end
    if isnothing(liniter)
        nothing
    else
        el, outstate = liniter
        el, (lin, outstate)
    end
end

"""
    num_nodes(node::DAG)

Count the number of nodes in the `DAG`
"""
function num_nodes(node::DAG, seen=nothing)
    count::Int = 0
    foreach(node, seen) do n
        count += 1
    end
    count
end

"Number of edges in the `DAG`"
function num_edges(node::DAG, seen=nothing)
    count::Int = 0
    foreach(node, seen) do n
        count += num_children(n)
    end
    count
end

"""
    tree_num_nodes(node::DAG)::BigInt

Compute the number of nodes in of a tree-unfolding of the `DAG`.
"""
function tree_num_nodes(node::DAG, cache=nothing)
    @inline f_leaf(n) = one(BigInt)
    @inline f_inner(n, call) = (1 + mapreduce(call, +, children(n)))
    foldup(node, f_leaf, f_inner, BigInt, cache)
end

"""
    tree_num_edges(node::DAG)::BigInt

Compute the number of edges in the tree-unfolding of the `DAG`.
"""
function tree_num_edges(node::DAG, cache=nothing)::BigInt
    @inline f_leaf(n) = zero(BigInt)
    @inline f_inner(n, call) = (num_children(n) + mapreduce(c -> call(c), +, children(n)))
    foldup(node, f_leaf, f_inner, BigInt, cache)
end

"Is the node contained in the `DAG`?"
function Base.in(needle::DAG, root::DAG, seen=nothing)
    contained::Bool = false
    foreach(root, seen) do n
        contained |= (n == needle)
    end
    contained
end

"Get the list of inner nodes in a given graph"
innernodes(c::DAG) =
    filter(isinner, c)

"Get the list of leaf nodes in a given graph"
leafnodes(c::DAG) =
    filter(isleaf, c)

"Count the number of leaf nodes in a given graph"
num_leafnodes(c::DAG) =
    length(leafnodes(c))

"Count the number of inner nodes in a given graph"
num_innernodes(c::DAG) =
    length(innernodes(c))

"Order the `DAG`'s nodes bottom-up in a list (with optional element type)"
@inline linearize(r::DAG, ::Type{T} = Union{}) where T =
    filter(x -> true, r, nothing, typejoin(T, typeof(r)))

"""
    left_most_descendent(root::DAG)::DAG

Return the left-most descendent.
"""
function left_most_descendent(root::DAG)::DAG
    while isinner(root)
        root = children(root)[1]
    end
    root
end

"""
    right_most_descendent(root::DAG)::DAG

Return the right-most descendent.
"""
function right_most_descendent(root::DAG)::DAG
    while isinner(root)
        root = children(root)[end]
    end
    root
end

"""
Assign an integer label to each circuit node, bottom up, starting at `1`
"""
function label_nodes(root::DAG)
    labeling = Dict{DAG,Int}()
    i = 0
    f_inner(n, call) = begin
        foreach(call, children(n))
        (i += 1)
    end
    f_leaf(n) = (i += 1)
    foldup(root, f_leaf, f_inner, Int, labeling)
    labeling
end

"""
    num_parents(root::DAG)

Count the number of parents for each node in `DAG` under `root`
"""
function num_parents(root::DAG)
    count = Dict{DAG,Int}()
    count[root] = 0
    foreach(root) do n
        if isinner(n)
            for c in children(n)
                count[c] = get(count, c, 0) + 1
            end
        end
    end
    count
end

"""
    feedforward_layers(root::DAG)

Assign a layer id with each node, starting from leafs to root
"""
function feedforward_layers(root::DAG)
    node2layer = Dict{DAG, Int}()
    f_inner(n, call) =
        1 + mapreduce(call, max, children(n))
    f_leaf(n) = 1
    num_layers = foldup(root, f_leaf, f_inner, Int, node2layer)
    node2layer, num_layers
end

#####################
# debugging methods (not performance critical)
#####################

# When you suspect there is a bug but execution halts, it may be because of
# pretty printing a huge recursive graph structure.
# To safeguard against that case, we set a default show:
Base.show(io::IO, c::DAG) = print(io, "$(typeof(c))($(hash(c)))")

"""
    node_stats(c::DAG)

Give count of types and fan-ins of all nodes in the graph
"""
node_stats(c::DAG) = merge(leaf_stats(c), innernode_stats(c))

"""
    innernode_stats(c::DAG)

Give count of types and fan-ins of inner nodes in the graph
"""
function innernode_stats(c::DAG)
    groups = groupby(e -> (typeof(e), num_children(e)), innernodes(c))
    map_values(v -> length(v), groups, Int)
end

"""
    leaf_stats(c::DAG)

Give count of types of leaf nodes in the graph
"""
function leaf_stats(c::DAG)
    groups = groupby(e -> typeof(e), leafnodes(c))
    map_values(v -> length(v), groups, Int)
end

"Group the elements of `list` by their values according to function `f`"
function groupby(f::Function, list::Union{Vector{E},Set{E}})::Dict{Any,Vector{E}} where E
    groups = Dict{Any,Vector{E}}()
    for v in list
        push!(get!(groups, f(v), []), v)
    end
    groups
end

"Map the values in the dictionary, retaining the same keys"
function map_values(f::Function, dict::AbstractDict{K}, vtype::Type)::AbstractDict{K,vtype} where K
    mapped_dict = Dict{K,vtype}()
    for key in keys(dict)
        mapped_dict[key] = f(dict[key])
    end
    mapped_dict
end

"""
    parent_stats(c::DAG)

Give number of nodes grouped by (type, parent_count)
"""
function parent_stats(c::DAG)
    par_count = num_parents(c)
    groups = groupby(e -> (typeof(e[1]), e[2]), collect(par_count));
    map_values(v->length(v), groups, Int)
end


text:
# This file transformed from Lazy Foo' Productions "Playing Sounds" tutorial:
# http://lazyfoo.net/SDL_tutorials/lesson11/index.php

using SimpleDirectMediaLayer.LibSDL2

SDL_Init(SDL_INIT_AUDIO)

#Load the music
aud_files = dirname(@__FILE__)
music = Mix_LoadMUS("$aud_files/beat.wav")

if (music == C_NULL)
    error("$aud_files/beat.wav not found.")
end

scratch = Mix_LoadWAV_RW(SDL_RWFromFile("$aud_files/scratch.wav", "rb"), 1)
high = Mix_LoadWAV_RW(SDL_RWFromFile("$aud_files/high.wav", "rb"), 1)
med = Mix_LoadWAV_RW(SDL_RWFromFile("$aud_files/medium.wav", "rb"), 1)
low = Mix_LoadWAV_RW(SDL_RWFromFile("$aud_files/low.wav", "rb"), 1)
Mix_PlayChannelTimed(-1, med, 0, -1)

Mix_PlayMusic(music, -1)
sleep(1)
Mix_PauseMusic()
sleep(1)
Mix_ResumeMusic()
sleep(1)
Mix_HaltMusic()

ext: jl


text:
function tree(tg, av, T)
    if av["txt"] == "Atgal | Повернутися"
        if av["path"][end] == "0"
            tbegin(tg, av)
            return
        else
            av["path"] = av["path"][1:end-1]
        end
    end
    pathkey = join(av["path"])
    ch = T[pathkey]["children"]
    children = []
    ktree = [[]]
    for i in 1:length(ch)
        if any(T[string(pathkey,ch[i])]["dav_id"] .== av["id"])
            push!(children,string(ch[i],"✓"))
        else
            push!(children,ch[i])
        end
    end
    for i in 1:2:length(children)
        if (length(children)-i) > 0
            push!(ktree,children[i:i+1])
        else
            push!(ktree,[children[i]])
        end
    end
    if any(T[pathkey]["dav_id"] .== av["id"]) && av["step"] == "enter"
        push!(ktree,["Atgal | Повернутися","Nuimti žymę","Namo | Додому"])
    else
        push!(ktree,["Atgal | Повернутися","Pasirinkti","Namo | Додому"])
    end
    popfirst!(ktree)
    if length(ktree) > 1
        msg = string(T[pathkey]["descript"]," Išsirinkite sritį.")
    else
        field = av["path"][end]
        msg = "Eikite atgal arba pasirinkite šią $field sritį."
    end
    kb = Dict(:keyboard => ktree, :one_time_keyboard => true, :resize_keyboard=>true)
    sendMessage(tg, chat_id = av["id"], text = msg, reply_markup = kb)

    return
end

function maketree(filename = "tree.txt")
    mtxt = readlines(filename, keep = false)
    steps = ["0"]
    E = Dict("field"=>"0","children"=>[],"steps"=>deepcopy(steps),"dav_id"=>[])
    T = Dict("0"=>Dict("field"=>"0","children"=>[],"steps"=>deepcopy(steps),"dav_id"=>[],"descript"=>""))
    lv = 0
    for i in 1:length(mtxt)
        it = findall("\t", mtxt[i])
        eile = mtxt[i][length(it)+1:end]
        if '/' in eile
            (eil,descript) = split(eile,'/')
        else
            eil = eile
            descript = ""
        end
        y = length(it)+1
        if y>lv
            step = join(steps)
            T[step]["children"] = eil
            append!(steps,[eil])
        elseif y<= lv
            steps = steps[1:end+y-lv]
            steps[end] = eil
            T[join(steps[1:end-1])]["children"] = vcat(T[join(steps[1:end-1])]["children"],eil)
        end
        step = join(steps)
        T[step] = deepcopy(E)
        T[step]["field"] = eil
        T[step]["steps"] = deepcopy(steps)
        T[step]["descript"] = descript
        lv = y
    end
    return T
end

function kas2nkas(name)
    d = Dict("as"=>"ai",'ė'=>"e","us"=>"au","ys"=>"y","is"=>"i","inkas"=>"inke")
    if length(name)>5 && haskey(d,name[end-4:end])
        nname=name[1:end-5]*d[name[end-4:end]]
    elseif haskey(d,name[end-1:end])
        nname=name[1:end-2]*d[name[end-1:end]]
    elseif haskey(d,name[end])
        nname=name[1:end-1]*d[name[end]]
    else
        nname=name
    end
end

function tik(k) #talka inline keyboard
    kb=[]
    for ki in k
        kl=[]
        for i in 1:2:length(ki)
            push!(kl,Dict(:text=>ki[i],:callback_data =>ki[i+1]))
        end
        push!(kb,kl)
    end
    d=Dict(:inline_keyboard => kb)
    return d
end

function getgroup(groups)
    for group in groups
        if group.second["state"]=="free"
            return group
        end
    end

    print("not enought free rooms")
end

function taskai(tg)
    sendMessage(tg, chat_id = av["id"], text = "Jūsų taškų skaičius $(round(10*av["token"])/10)")
end

function cleanK(d::DataBase)
    (T, Av, updateId, K, groups) = deserialize(DATAFILE[])
    d.K = Dict(now() => Dict("getter" => Av[5090964479], "txt"=>"refresh",
                             "giver"=>[1], "state"=>"nothing"))
    serialize(DATAFILE[], [T, Av, updateId, K, groups])
end

function cleanT()
    (T,Av,updateId,K,groups) = deserialize(DATAFILE[])
    ak=keys(Av)
    for t in T
        v=[]
        for i in 1:length(t.second["dav_id"])
            if !any(ak.==t.second["dav_id"][i])
                push!(v,i)
            end
        end
        deleteat!(t.second["dav_id"],v)
        t.second["dav_id"]=union(t.second["dav_id"])
    end
    serialize(DATAFILE[], [T,Av,updateId,K,groups])
    print("T išvalytas2")
end


text:
open("core/test1.txt", "r") do io
    for c in TrySource(io)
        print(c)
    end
end

open("core/test1.txt", "r") do io
    s = TrySource(io)
    i = start(s)
    (c, i) = next(s, i)
    @test c == 'a'
    @test forwards(s, i) == "bcdefghijklmnopqrstuvwxyz\n"
end

#open("test1.txt", "r") do io
#    parse_one_dbg(TrySource(io), Trace(p"[a-z]"[0:end] + e"m" > string); debug=true)
#end

for parse in (parse_try, parse_try_cache, parse_try_dbg, parse_try_cache_dbg)

    open("core/test1.txt", "r") do io
        @test_throws ParserException parse(io, Trace(p"[a-z]"[0:end] + e"m" > string))
    end

    open("core/test1.txt", "r") do io
        result = parse(io, Try(p"[a-z]"[0:end] + e"m" > string))
        println(result)
        @test result == Any["abcdefghijklm"]
    end

    open("core/test1.txt", "r") do io
        # multiple lines
        result = parse(io, Try(p"(.|\n)"[0:end] + e"5" > string))
        println(result)
        @test result == Any["abcdefghijklmnopqrstuvwxyz\n012345"]
    end

    @test_throws ParserError{LineIter} parse("?", Alt!(p"[a-z]", p"\d", Error("not letter or number")))

end

println("try ok")


text:
using EngEconomics, Roots, Plots

# Given
init = -610000
annual = 200000
N = 10
salvage = -1500000

MARR = 0.1
# Determine the better alternative
NPW(x) = init + annual * seriesPresentAmountFactor(x, N) + salvage * presentWorthFactor(x, N)
xVec = collect(0:0.01:1.0)
zVec = zeros(size(xVec))

# Two irrs that are equal to zero
plot(xVec, NPW.(xVec))
plot!(xVec, zVec)


text:
# rational number reconstruction implementation borrowed from CLUE
# and modified a bit to suit the 'Modern Computer Algebra' definitions
# returns a rational r // h of QQ field in a canonical form such that
#   r // h ≡ a (mod m)
#
# let n = max( λ(a), λ(m) ) , where λ(x) is a number of bits for x
# O(n^2)
"""
    rational_reconstruction(a, m)
Rational number reconstruction implementation borrowed from CLUE
and modified a bit to suit the 'Modern Computer Algebra' definitions.
Returns a rational of QQ field in a canonical form that
is congruent a modulo m
a, m are integers
"""
function rational_reconstruction(a::I, m::I) where {I<:Union{Int, BigInt}}
    # restricting to I<:Union{Int, BigInt} for compatibility reasons

    a = mod(a, m)
    if a == 0 || m == 0
        return Nemo.QQ(0, 1)
    end
    if m < 0
        m = -m
    end
    if a < 0
        a = m - a
    end
    if a == 1
        return Nemo.QQ(1, 1)
    end
    bnd = sqrt(float(m) / 2)

    U = I[1, 0, m]
    V = I[0, 1, a]
    while abs(V[3]) >= bnd
        q = div(U[3], V[3])
        T = U .- q .* V
        U = V
        V = T
    end

    t = abs(V[2])
    r = V[3] * sign(V[2])
    # changed from `<= bnd` to `<= m / bnd`
    if t <= m / bnd && gcd(r, t) == 1
        return Nemo.QQ(r, t)
    end

    throw(DomainError(
        :($a//$m), "rational reconstruction of $a (mod $m) does not exist"
    ))
end

### Several ad hoc functions for reconstruction

function rational_reconstruction(f::Nemo.gfp_elem, m)
    rational_reconstruction(BigInt(f.data), m)
end

function rational_reconstruction(f::Nemo.fmpz, m)
    rational_reconstruction(BigInt(f), m)
end

function rational_reconstruction(f::AbstractAlgebra.GFElem{T}, m) where {T}
    rational_reconstruction(AbstractAlgebra.lift(f), m)
end

function rational_reconstruction(f::Union{MPoly, gfp_mpoly, fmpz_mpoly}, m)
    map_coefficients(c -> rational_reconstruction(c, m), f)
end

function rational_reconstruction(f::Frac, m)
    q1 = rational_reconstruction(numerator(f), m)
    q2 = rational_reconstruction(denominator(f), m)
    q1 // q2
end

function rational_reconstruction(A::Array, m)
    map(c -> rational_reconstruction(c, m), A)
end

function rational_reconstruction(f::Singular.n_Zp, m)
    rational_reconstruction(BigInt(Int(f)), m)
end

#------------------------------------------------------------------------------

"""
    Coerces the given rational field element into the given finite field
"""
function modular_reduction(x, field)
    n, d = field(Int(numerator(x))), field(Int(denominator(x)))
    if iszero(d)
        throw(DomainError(
            :($x), "modular reduction of $x (to $field) does not exist"
        ))
    end
    n // d
end

### Reduction ad hoc functions

function modular_reduction(f::Union{MPoly, fmpq_mpoly}, field)
    map_coefficients(c -> modular_reduction(c, field), f)
end

function modular_reduction(f::Frac, field)
    modular_reduction(numerator(f), field) //
        modular_reduction(denominator(f), field)
end

function modular_reduction(x::gfp_fmpz_elem, field)
    return field(convert(BigInt, x))
end

function modular_reduction(x::gfp_elem, field)
    return field(x)
end

function modular_reduction(arr::Array, field)
    map(f -> modular_reduction(f, field), arr)
end

#------------------------------------------------------------------------------

"""
    Chinese remainder algorithm
    (Eucleadean one. Ii there a Gauss version of the algorithm btw?..)
"""
function CRT(moduli, remainders)
    """
    The function implements the CR Algorithm mentioned as Algorithm 5.4
    in  'Modern Computer Algebra', 2nd edition,
        Joachim von zur Gathen and Jurgen Gerhard
    """

    m = prod(moduli)

    ans = BigInt(0)
    for (i, (mi, vi)) in enumerate(zip(moduli, remainders))
        mi_hat = divexact(m, mi)
        _, s, t = Nemo.gcdx(mi_hat, mi)
        ans = BigInt(mod(ans + vi * s * mi_hat, m))
    end

    ans
end

"""
    Reconstructs the ultimate Groebner basis with coefficients
    obtained from applying CRT to the given bases coeffs and moduli
"""
function CRT_reconstruction(groebner_bases, moduli)
    onegb = first(groebner_bases)
    onepoly = first(onegb)

    basepolyring = base_ring(base_ring(onepoly))
    yoverx = parent(onepoly)

    ystring = string.(symbols(yoverx))
    xstring = string.(symbols(basepolyring))

    ground = Nemo.ZZ
    inner_ring, = Nemo.PolynomialRing(ground, xstring)
    outer_ring, = Nemo.PolynomialRing(FractionField(inner_ring), ystring)

    npolys = length(onegb)
    ncoeffs = map(length, onegb)

    new_basis = []

    for (polyidx, poly) in zip(1:npolys, onegb)
        outer_f = deepcopy(poly)
        outer_builder = AbstractAlgebra.MPolyBuildCtx(outer_ring)

        for (outer_i, outer_e) in zip(1:ncoeffs[polyidx], exponent_vectors(poly))
            inner_frac = coeff(poly, outer_e)
            num_and_denom = [zero(inner_ring), zero(inner_ring)]

            for (jk, locator) in enumerate([numerator, denominator])
                inner_f = locator(inner_frac)
                inner_builder = AbstractAlgebra.MPolyBuildCtx(inner_ring)

                for (inner_i, inner_e) in enumerate(exponent_vectors(inner_f))
                    oldcoeffs = [
                        coeff(locator(coeff(gb[polyidx], outer_e)), inner_e)
                        for gb in groebner_bases
                    ]

                    newcoeff = CRT(moduli, map(x -> BigInt(x.data), oldcoeffs))
                    push_term!(inner_builder, newcoeff, inner_e)
                end

                num_and_denom[jk] = finish(inner_builder)

            end

            num, denom = num_and_denom
            push_term!(outer_builder, num // denom, outer_e)
        end

        push!(new_basis, finish(outer_builder))
    end

    new_basis
end


text:
module PiecewiseInterpolation

using IntervalSets
using Dierckx
import Base: append!

export PiecewiseSpline1D, append!

immutable PiecewiseSpline1D{T <: Real}
    pieces::Vector{Spline1D}
    intervals::Vector{ClosedInterval{T}}
end

function (p::PiecewiseSpline1D)(t)
    i = findlast(interval -> t ∈ interval, p.intervals)
    return p.pieces[i](t)
end

function split_interval{T <: Real}(jumps::Vector{T},
                                   interval::ClosedInterval{T})
    isempty(jumps) && return [interval]
    intervals = Vector{typeof(interval)}(length(jumps) + 1)
    intervals[1] = interval.left..jumps[1]
    for i in 2:length(jumps)
        intervals[i] = jumps[i - 1]..jumps[i]
    end
    intervals[end] = jumps[end]..interval.right
    return intervals
end

function get_index_range{T}(times::Vector{T}, interval::ClosedInterval{T})
    f = t -> t ∈ interval
    return findfirst(f, times):findlast(f, times)
end

get_interpolation(times, values) =
    Spline1D(times, values, bc="extrapolate", k = min(5, length(times)))

function append!{T}(p::PiecewiseSpline1D{T}, new_times::Vector{T},
                    new_values::Vector, new_jumps::Vector{T})
    old_t_end = p.intervals[end].right
    if any(new_times .< old_t_end) || any(new_jumps .< old_t_end)
        error("new data overlaps with existing interpolations")
    end
    intervals = split_interval(new_jumps, old_t_end..new_times[end])
    for interval in intervals
        indices = get_index_range(new_times, interval)
        interpolation = get_interpolation(view(new_times, indices),
                                          view(new_values, indices))
        push!(p.pieces, interpolation)
        push!(p.intervals, interval)
    end
    return nothing
end

function PiecewiseSpline1D{T <: Real}(times::Vector{T}, values::Vector,
                                      jumps::Vector{T})
    assert(length(times) == length(values))
    # create object with first interval, then append remaining
    time_interval = times[1]..jumps[1]
    indices = 1:searchsortedlast(times, jumps[1])
    interpolation =
        get_interpolation(view(times, indices), view(values, indices))
    p = PiecewiseSpline1D([interpolation], [time_interval])
    rest_start = searchsortedfirst(times, jumps[1])
    append!(p, times[rest_start:end], values[rest_start:end], jumps[2:end])
    return p
end

end #module


text:
push!(LOAD_PATH, joinpath(homedir(),"GitHub/KenyaCoV/src"))
using Plots,Parameters,Distributions,DifferentialEquations,JLD2,DataFrames,StatsPlots,FileIO,MAT,RecursiveArrayTools
import KenyaCoV
using LinearAlgebra:eigen
using Statistics: median, quantile

"""
Load age structured data, define initial state and declare the KenyaCoV problem
"""

u0,P,P_dest = KenyaCoV.model_ingredients_from_data("data/data_for_age_structuredmodel.jld2",
                                            "data/flight_numbers.csv",
                                            "data/projected_global_prevelance.csv")
@load "data/detection_rates_for_different_taus.jld2" d_0 d_01 d_025 d_05 d_1

P.χ = ones(KenyaCoV.n_a)
P.rel_detection_rate = d_0
P.dt = 0.25;
P.ext_inf_rate = 0.;
P.ϵ = 0.
P.γ = 1/2.5
P.σ = 1/rand(KenyaCoV.d_incubation)
P.β = 2.5*P.γ
R₀_scale = KenyaCoV.calculate_R₀_scale(P)
P.χ = ones(KenyaCoV.n_a)/R₀_scale
KenyaCoV.calculate_R₀(P) #Before

R_kenya_tau_0 = [R_kenya_prediction(1,ones(16),P.χ,0.2,0.,x,0,1/2.5) for x in x_range]
R_kenya_tau_7 = [R_kenya_prediction(1,ones(16),P.χ,0.2,0.,x,1/7,1/2.5) for x in x_range]
R_kenya_tau_3_5 = [R_kenya_prediction(1,ones(16),P.χ,0.2,0.,x,1/3.5,1/2.5) for x in x_range]

plot((1 .- collect(x_range))*100,R_kenya_tau_0,lw=3,lab = "No isolation")
plot!((1 .- collect(x_range))*100,R_kenya_tau_7,lw=3,lab = "Av. 1 week to isolation")
plot!((1 .- collect(x_range))*100,R_kenya_tau_3_5,lw=3,lab = "Av. 1/2 week to isolation")
plot!([0,50],[1/2.5,1/2.5],lw=3,lab = "Threshold for Chinese reproductive ratio = 2.5",ls= :dash,color=:black)
xlabel!("Percentage reduction in transmissibility of clinical cases")
ylabel!("Relative reproductive ratio for Kenya")
title!("Effectiveness of intervention: SARS-like scenario")
savefig("plotting/eff_SARS_like.png")

e_range = 0:0.01:0.5
MERS_R_tau_0 = [R_kenya_prediction(1,reporting_rate,ones(16),0.2,e,1.,0,1/2.5) for e in e_range]
MERS_R_tau_max = [R_kenya_prediction(1,reporting_rate,ones(16),0.2,e,0.5,1/3.5,1/2.5) for e in e_range]

MERS_R_tau_0_alt = [R_kenya_prediction(1,reporting_rate,ones(16),0.8,e,1.,0,1/2.5) for e in e_range]
MERS_R_tau_max_alt = [R_kenya_prediction(1,reporting_rate,ones(16),0.8,e,0.5,1/3.5,1/2.5) for e in e_range]

plot(e_range,MERS_R_tau_0,lab = "No intervention",lw=3,color = :red,legend = :bottomright)
plot!(e_range,MERS_R_tau_0_alt,lab = "Baseline 80% case detection",lw=3,color=:blue)
plot!(e_range,MERS_R_tau_max,lab = "Max. intervention",lw=3,color = :blue)
plot!(e_range,MERS_R_tau_max_alt,lab = "Baseline 80% case detection - max. intervention",lw=3,color = :blue,ls = :dot)
plot!([0,0.5],[1/2.5,1/2.5],lw=3,lab = "Threshold for Chinese reproductive ratio = 2.5",ls= :dash,color=:black)
xlabel!("Relative infectiousness of undetected cases")
ylabel!("Relative reproductive ratio for Kenya")
title!("Effectiveness of intervention: MERS-like scenario")

savefig("plotting/eff_MERS_like.png")


text:
##### Beginning of file

testmodulea_filename = joinpath("TestModuleA", "TestModuleA.jl")
testmoduleb_filename  = joinpath(
    "TestModuleB", "directory1", "directory2", "directory3",
    "directory4", "directory5", "TestModuleB.jl",
    )
testmodulec_filename = joinpath(mktempdir(), "TestModuleC.jl")
rm(testmodulec_filename; force = true, recursive = true)
open(testmodulec_filename, "w") do io
    write(io, "module TestModuleC end")
end
include(testmodulea_filename)
include(testmoduleb_filename)
include(testmodulec_filename)

##### End of file


text:
import Adapt
import CUDA

Adapt.adapt_structure(to, data::IJKFVH{S, Nij, Nk}) where {S, Nij, Nk} =
    IJKFVH{S, Nij, Nk}(Adapt.adapt(to, parent(data)))

Adapt.adapt_structure(to, data::IJFH{S, Nij}) where {S, Nij} =
    IJFH{S, Nij}(Adapt.adapt(to, parent(data)))

Adapt.adapt_structure(to, data::IFH{S, Ni}) where {S, Ni} =
    IFH{S, Ni}(Adapt.adapt(to, parent(data)))

Adapt.adapt_structure(to, data::IJF{S, Nij}) where {S, Nij} =
    IJF{S, Nij}(Adapt.adapt(to, parent(data)))

Adapt.adapt_structure(to, data::IF{S, Ni}) where {S, Ni} =
    IF{S, Ni}(Adapt.adapt(to, parent(data)))

Adapt.adapt_structure(to, data::VF{S}) where {S} =
    VF{S}(Adapt.adapt(to, parent(data)))

parent_array_type(::Type{<:CUDA.CuArray{T, N, B} where {N}}) where {T, B} =
    CUDA.CuArray{T, N, B} where {N}

# Ensure that both parent array types have the same memory buffer type.
promote_parent_array_type(
    ::Type{CUDA.CuArray{T1, N, B} where {N}},
    ::Type{CUDA.CuArray{T2, N, B} where {N}},
) where {T1, T2, B} = CUDA.CuArray{promote_type(T1, T2), N, B} where {N}

# Make `similar` accept our special `UnionAll` parent array type for CuArray.
Base.similar(
    ::Type{CUDA.CuArray{T, N′, B} where {N′}},
    dims::Dims{N},
) where {T, N, B} = similar(CUDA.CuArray{T, N, B}, dims)

function knl_copyto!(dest, src)

    #=
    nij, nh = size(dest)

    thread_idx = CUDA.threadIdx().x
    block_idx = CUDA.blockIdx().x
    block_dim = CUDA.blockDim().x

    # mapping to global idx to make insensitive
    # to number of blocks / threads per device
    global_idx = thread_idx + (block_idx - 1) * block_dim

    nx, ny = nij, nij
    i = global_idx % nx == 0 ? nx : global_idx % nx
    j = cld(global_idx, nx)
    h = ((global_idx-1) % (nx*nx)) + 1
    =#

    i = CUDA.threadIdx().x
    j = CUDA.threadIdx().y

    h = CUDA.blockIdx().x

    p_dest = slab(dest, h)
    p_src = slab(src, h)

    @inbounds p_dest[i, j] = p_src[i, j]
    return nothing
end

function Base.copyto!(
    dest::IJFH{S, Nij},
    bc::Union{IJFH{S, Nij, A}, Base.Broadcast.Broadcasted{IJFHStyle{Nij, A}}},
) where {S, Nij, A <: CUDA.CuArray}
    _, _, _, _, Nh = size(bc)
    CUDA.@cuda threads = (Nij, Nij) blocks = (Nh,) knl_copyto!(dest, bc)
    return dest
end


text:
#=
Popis
  Funkcia zisti vyskyt sluciek v grafe.

Parametre
  g : datova struktura grafu :: GenericGraph

Navratova hodnota
  : parameter reprezentujuci najdenie slucky(True=najdena/False=nenajdena) :: Bool

Poznamky
  Pri najdeni prvej slucky v grafe je vykonavanie funkcie ukoncene.
=#

function contains_selfloops(g::GenericGraph)

  for e in edges(g)
    if e.source == e.target
      return true
    end
  end

  return false

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis


text:
#SimulationOptions Type
struct SimulationOptions
  ODE_SolverType::String
  DofSimType::String
  SimulationTime::Float64

  function SimulationOptions(;ODE_SolverIn::String = "tsit5",
                    DofSimTypeIn::String = "Trans", SimulationTime::Float64 = 50.0)
    #ODE Type Checker
    if !haskey(ODETYPES, lowercase(ODE_SolverIn))
      warn("ODE Solver Type Not Found, Setting Solver Type to Default (Tsit5)")
      ODE_SolverIn::String = "tsit5"
    end

    #DofSimType Checker
    if !haskey(DOFSIMTYPES, lowercase(DofSimTypeIn))
      warn("Degrees of Freedom Type Not Found, Setting Degrees
            of Freedom Type to Default (Trans)")
      DofSimTypeIn::String = "trans"
    end

    #SimulationTime Checker
    if SimulationTime <= 0
        warn("Simulation Time must be positive, Setting Simulation Time
            to Default (50.0 sec)")
        SimulationTime = 50.0
    end

    this = new(lowercase(ODE_SolverIn),
               lowercase(DofSimTypeIn));

  end
end


text:
module EscapeString

using HttpCommon

export escapestring

# Registry of available output formats and their handlers
formats = Dict{Symbol, Function}()
formats[:html_text] = escapeHTML

"Escapes string s to format fmt."
function escapestring(s::AbstractString, fmt::Symbol)
    formats[fmt](s)
end

end # module


text:
function commandHelp(ed::Editor, args::String)
    if args == ""
        setStatusMessage(ed, "Type help <command> for command specific help")
        return
    end

    if !Base.isidentifier(args)
        setStatusMessage(ed, "help: '$args' is not a valid command name")
        return
    end

    sym = Symbol(args)

    # Try getting help information from commands and params
    helptext = ""

    if configIsParam(sym)
        helptext = configDesc(sym)
    elseif sym in keys(COMMANDS)
        helptext = COMMANDS[sym].help
    end

    if helptext == ""
        setStatusMessage(ed, "No help documents for '$args'")
    else
        setStatusMessage(ed, helptext)
    end
end

addCommand(:help, commandHelp,
    help="Type help <command> for command specific help")


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

using JSON

push!(LOAD_PATH,joinpath(pwd(),"/../"))
using MDP
import MCTSdpw: DPW, State, Action, StateNode, StateActionNode, StateActionStateNode

function saveSimTree(dpw::DPW,s::State,outfile::AbstractString)

  function process(sn::StateNode,Tout)
    for (a,san) in sn.a
      node_ = Dict{ASCIIString, Any}()
      node_["action"] = hash(a)
      node_["N"] = san.n
      node_["Q"] = san.q
      node_["states"] = Dict{ASCIIString, Any}[]

      push!(Tout, node_)

      process(san, node_["states"])
    end
  end

  function process(san::StateActionNode, Tout)
    for (s, sasn) in san.s
      node_ = Dict{ASCIIString, Any}()
      node_["state"] = hash(s)
      node_["N"] = haskey(dpw.s,s) ? dpw.s[s].n : sasn.n
      node_["r"] = sasn.r
      node_["actions"] = Dict{ASCIIString, Any}[]

      push!(Tout, node_)

      if haskey(dpw.s,s)
        process(dpw.s[s], node_["actions"])
      end
    end
  end

  function process{S<:State}(ss::Vector{S}, Tout)
    for s in ss
      node_ = Dict{ASCIIString, Any}()
      node_["state"] = hash(s)
      node_["N"] = dpw.s[s].n
      node_["actions"] = Dict{ASCIIString, Any}[]

      push!(Tout, node_)

      if haskey(dpw.s,s)
        process(dpw.s[s], node_["actions"])
      end
    end
  end

  Tout = Dict{ASCIIString, Any}()
  Tout["name"] = "root"
  Tout["states"] = Dict{ASCIIString, Any}[]

  process([s],Tout["states"])
  f = open(outfile, "w")
  JSON.print(f, Tout, 2)
  close(f)

  return Tout
end

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
#!/usr/bin/env julia

function solve(case, N, A, B, C)
    filler = N - A - B + C
    if filler < 0
        println("Case #$case: IMPOSSIBLE")
        return
    end

    str = []
    if A == 1
        append!(str, Vector((N - A + C):N))
        append!(str, repeat([N], C - 1))
        append!(str, repeat([1], filler))
        append!(str, Vector((N - 1):-1:(N - B + C)))
    else
        append!(str, Vector((N - A + C):(N - 1)))
        append!(str, repeat([N], C - 1))
        append!(str, repeat([1], filler))
        append!(str, Vector(N:-1:(N - B + C)))
    end
    println("Case #$case: $(join(str, " "))")
end

function main()
    T = parse(Int, readline(stdin))
    for case in 1:T
        N, A, B, C = parse.(Int, stdin |> readline |> split)
        solve(case, N, A, B, C)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
include("./bbh_to_compare_with_sbh.jl")

s = [ p1 ]

for x in s
    solve_wave_equation_2D(x)
end


text:
using MyPackage
using Test

@testset "MyPackage.jl" begin

    @test f(10) == 100
    @test f2(10,20) == 30
    # Write your tests here.
end


text:
cd(@__DIR__)
using Pkg
Pkg.activate(".")

##

using Revise
using Statistics
using ProBeTune
using DiffEqFlux
using LaTeXStrings

##
include("PlotUtils.jl")
path = "../plots/kuramoto/" # for saving plots
mkpath(path)
save_plots = true # set to true
##
#=
Define global system parameters to use in the tune_run function
=#
global t_steps = 0.:0.1:10.
global tspan = (0., 10.)
global N_osc = 10
global size_p_sys = N_osc * (N_osc - 1) ÷ 2 + N_osc
global size_p_spec = 2
global N_samples = 10
global p_sys_init = 6. * rand(size_p_sys) .+ 1.
global p_spec_init = rand(size_p_spec) .+ 1.
global p_initial = vcat(p_sys_init, repeat(p_spec_init, N_samples))
global K_av = 1.
##
#=
Create a set of intrinsic frequencies with mean(omega)=0
In tune_run we use omega*scale with variating scales, to explore systems with different synchronisity.
=#
omega = randn(N_osc)
omega[1] = 0.
omega .-= mean(omega[1:end])

##
#=
Tuning pipeline for the system of kuramoto oscillators.
We want to tune a system of 10 oscillators with one grid connection to behave like a single oscillator (specification) in terms of input-output behavior.
We seek to find such parameters for the system, that it behaves close to some single oscillator under all possible inputs.
In this example we use a sample of size 10.

First, we estimate the distance and optimize only specification parameters by running behavioural_distance.
Then we tune the system to specification using different optimizers.
We verify the tuning by resampling and tuning to the new sample.
All systems have the same set of inputs, as we fix the random seed.
=#
function tune_run(omega, scale; save_plots = false)
    Random.seed!(1)
    kur = PBTExample.create_kuramoto_example(omega*scale, N_osc, size_p_spec, K_av, t_steps, N_samples)

    d, p_dist, = behavioural_distance(kur, p_initial; abstol=1e-4, reltol=1e-4,
                                        optimizer=DiffEqFlux.BFGS(),
                                        optimizer_options=(
                                            :maxiters => 10,
                                            :cb => PBTLibrary.basic_pbt_callback))
    d_initial = d
    println("Initial distance to specified behaviour is lower than $d\n")
    ##
    scen = 1:3
    ##
    plot_callback(kur, p_dist, d, scenario_nums = scen, title = "Scenarios $scen, s = $scale", offset = 0.5, xlims = (kur.t_span[2]/2, kur.t_span[2]))
    if save_plots
        savefig(path*string(round(scale, digits = 4))*"_initial_"*string(round(d, digits = 4))*".png")
    end
    ##
    res = pbt_tuning(kur, p_dist; abstol=1e-4, reltol=1e-4)
    p_tuned = res.minimizer
    ##
    d2, p_dist_2, = behavioural_distance(kur, p_tuned; abstol=1e-4, reltol=1e-4)

    println("After a first tuning using ADAM(0.01) the distance to specified behaviour is lower than $d2\n")

    ##

    res = pbt_tuning(kur, p_dist_2; abstol=1e-4, reltol=1e-4,
                    optimizer=DiffEqFlux.BFGS(),
                    optimizer_options=(
                        :maxiters => 5,
                        :cb => PBTLibrary.basic_pbt_callback))

    p_t2 = res.minimizer
    ##
    plot_callback(kur, p_t2, res.minimum, scenario_nums = scen, title = "Scenarios $scen, s = $scale", offset = 0.5, xlims = (kur.t_span[2]/2, kur.t_span[2]))

    if save_plots
        savefig(path*string(round(scale, digits = 4))*"_tuned_"*string(round(res.minimum, digits = 4))*".png")
    end
    ##
    d3, p_dist_3, = behavioural_distance(kur, p_t2; abstol=1e-4, reltol=1e-4, optimizer_options=(:maxiters => 1000,))

    println("After a second tuning using BFGS the distance to specified behaviour is lower than $d3\n")
    ##
    plot_callback(kur, p_dist_3, d3, scenario_nums = scen, title = "Scenarios $scen, s = $scale", offset = 0.5, xlims = (kur.t_span[2]/2, kur.t_span[2]))
    if save_plots
        savefig(path*string(round(scale, digits = 4))*"_tuned2_"*string(round(d3, digits = 4))*".png")
    end
    resample!(kur)
    d4, p_dist_4, = behavioural_distance(kur, p_dist_3; abstol=1e-4, reltol=1e-4,
                    optimizer = DiffEqFlux.BFGS(),
                    optimizer_options=(
                        :maxiters => 10,
                        :cb => PBTLibrary.basic_pbt_callback))

    d5, p_dist_5, = behavioural_distance(kur, p_dist_4; abstol=1e-4, reltol=1e-4,
                    optimizer = DiffEqFlux.AMSGrad(0.01),
                    optimizer_options=(
                        :maxiters => 1000,))

    plot_callback(kur, p_dist_5, d5, scenario_nums = scen, title = "Scenarios $scen, s = $scale", offset = 0.5, xlims = (kur.t_span[2]/2, kur.t_span[2]))
    if save_plots
        savefig(path*string(round(scale, digits = 4))*"_final_rs_"*string(round(d5, digits = 4))*".png")
    end
    println("After resampling the distance to specified behaviour is lower than $d4\n")
    d_final = d5
    p_final = p_dist_5
    d_initial, d_final, p_final
end
##
#=
Different scales represent degree of synchronisity, the higher the number the less synchronous the system is.
Non-synchronous systems are more difficult to tune, so we chose a non-uniform distribution of scales.
Thus we can explore the large range of possible systems, while keeping the computation time not too high.
=#
scales = vcat(1.0:0.1:1.5,2.0:0.5:5.)

##
#=
Run tune_run for the set of systems, defined by scales.
res[i] contains initial and final distances to the specification, as well as the final parameters of the tuned system.
=#
res = [tune_run(omega, s, save_plots = save_plots) for s in scales]

##
#=
We want to evaluate performance of the tuning on all systems.
q = d_initial/d_final
First we plot the improvement depending on scale.
Secondly, we look at the absolute value of final distance.
=#
q = [r[1]/r[2] for r in res]

plot(scales, q)
if save_plots
    savefig(path*"improvement(scale).png")
end

plot(scales, [r[2] for r in res], xlabel = "s", ylabel = "d", label = "Final distance to specification")
if save_plots
    savefig(path*"final_d(omega).png")
end
##
#=
By running tune_run on all systems, we kept tuning parameters constant for all of them.
However, the chosen tuning parameters were limited by performance on systems with larger frequency spread.
In other words, we stop tuning when the systems with larger s start to become unstable.
At that point we could still continue tuning for more synchronous systems to show the best tuning result.
We chose to continue optimization for s = 1.2
=#
kur_new = PBTExample.create_kuramoto_example(omega*1.2, N_osc, size_p_spec, K_av, t_steps, N_samples)
p_initial_new = res[3][3] # for scale = 1.2

d, p_dist, = behavioural_distance(kur_new, p_initial_new; abstol=1e-4, reltol=1e-4,
                    optimizer = DiffEqFlux.BFGS(),
                    optimizer_options=(
                        :maxiters => 10,
                        :cb => PBTLibrary.basic_pbt_callback))
##
scen = 1:3
plot_callback(kur_new, p_dist, d, scenario_nums = scen, title = "Scenarios $scen, s = 1.2", offset = 0.5, xlims = (kur_new.t_span[2]/2, kur_new.t_span[2]))

##
d, p_dist, = behavioural_distance(kur_new, p_initial_new; abstol=1e-4, reltol=1e-4,
                    optimizer = DiffEqFlux.BFGS(),
                    optimizer_options=(
                        :maxiters => 10,
                        :cb => PBTLibrary.basic_pbt_callback))

plot_callback(kur_new, p_dist, d, scenario_nums = scen, title = "Scenarios $scen, s = 1.2", offset = 0.5, xlims = (kur_new.t_span[2]/2, kur_new.t_span[2]))
##
res = pbt_tuning(kur_new, p_dist)

plot_callback(kur_new, res.minimizer, res.minimum, scenario_nums = scen, title = "Scenarios $scen, s = 1.2", offset = 0.5, xlims = (kur_new.t_span[2]/2, kur_new.t_span[2]))
if save_plots
    savefig(path*"1.2_tuned_3_"*string(round(res.minimum, digits = 4))*".png")
end

resample!(kur_new)

d_rs, p_dist_rs = behavioural_distance(kur_new, res.minimizer; abstol=1e-4, reltol=1e-4,
optimizer = DiffEqFlux.BFGS(),
optimizer_options=(
    :maxiters => 10,
    :cb => PBTLibrary.basic_pbt_callback))

plot_callback(kur_new, p_dist_rs, d_rs, scenario_nums = scen, title = "Scenarios $scen, s = 1.2", offset = 0.5, xlims = (kur_new.t_span[2]/2, kur_new.t_span[2]))
if save_plots
    savefig(path*"1.2_rs_further_"*string(round(d_rs, digits = 4))*".png")
end
##
res2 = pbt_tuning(kur_new, p_dist_rs)

plot_callback(kur_new, res2.minimizer, res2.minimum, scenario_nums = scen, title = "Scenarios $scen, s = 1.2", offset = 0.5, xlims = (kur_new.t_span[2]/2, kur_new.t_span[2]))
if save_plots
    savefig(path*"1.2_rtuned4_"*string(round(res2.minimum, digits = 4))*".png")
end
##
resample!(kur_new)

d_rs, p_dist_rs, losses = behavioural_distance(kur_new, res.minimizer; abstol=1e-4, reltol=1e-4,
optimizer = DiffEqFlux.BFGS(),
optimizer_options=(
    :maxiters => 10,
    :cb => PBTLibrary.basic_pbt_callback))

plot_callback(kur_new, p_dist_rs, d_rs, scenario_nums = scen, title = "Scenarios $scen, s = 1.2", offset = 0.5, xlims = (kur_new.t_span[2]/2, kur_new.t_span[2]))
if save_plots
    savefig(path*"1.2_rs_final_"*string(round(d_rs, digits = 4))*".png")
end


text:
"""
    ProductModel{T <: VecOrTupOrNTup{<:AbstractMillModel}, U <: ArrayModel} <: AbstractMillModel

A model node for processing [`ProductNode`](@ref)s. For each subtree of the data node it applies one
(sub)model from `ms` and then applies `m` on the concatenation of results.

# Examples
```jldoctest product_model
julia> Random.seed!(0);

julia> n = ProductNode((a=ArrayNode([0 1; 2 3]), b=ArrayNode([4 5; 6 7])))
ProductNode with 2 obs
  ├── a: ArrayNode(2×2 Array with Int64 elements)
  └── b: ArrayNode(2×2 Array with Int64 elements)

julia> m1 = ProductModel((a=ArrayModel(Dense(2, 2)), b=ArrayModel(Dense(2, 2))))
ProductModel … ↦ ArrayModel(identity)
  ├── a: ArrayModel(Dense(2, 2))
  └── b: ArrayModel(Dense(2, 2))
```

```jldoctest product_model; filter=$(DOCTEST_FILTER)
julia> m1(n)
4×2 ArrayNode{Matrix{Float32}, Nothing}:
 -1.284...  -1.254...
  1.802...   3.711...
 -4.036...  -5.013...
  0.587...   0.372...
```

```jldoctest product_model
julia> m2 = ProductModel((a=identity, b=identity))
ProductModel … ↦ ArrayModel(identity)
  ├── a: ArrayModel(identity)
  └── b: ArrayModel(identity)

julia> m2(n)
4×2 ArrayNode{Matrix{Int64}, Nothing}:
 0  1
 2  3
 4  5
 6  7
```

See also: [`AbstractMillModel`](@ref), [`AbstractProductNode`](@ref), [`ProductNode`](@ref).
"""
struct ProductModel{T <: VecOrTupOrNTup{AbstractMillModel}, U <: ArrayModel} <: AbstractMillModel
    ms::T
    m::U
end

Flux.@functor ProductModel

"""
    ProductModel(ms, m=identity_model())

Construct a [`ProductModel`](@ref) from the arguments. `ms` should an iterable
(`Tuple`, `NamedTuple` or `Vector`) of one or more [`AbstractMillModel`](@ref)s, and `m` should be
an [`ArrayModel`](@ref).

It is also possible to pass any function (`Flux.Dense`, `Flux.Chain`, `identity`...) as elements of `ms`
or as `m`. In that case, they are wrapped into an [`ArrayNode`](@ref).

If `ms` is [`AbstractMillModel`](@ref), a one-element `Tuple` is constructed from it.

# Examples
```jldoctest
julia> ProductModel((a=ArrayModel(Dense(2, 2)), b=identity))
ProductModel … ↦ ArrayModel(identity)
  ├── a: ArrayModel(Dense(2, 2))
  └── b: ArrayModel(identity)

julia> ProductModel((identity_model(), BagModel(ArrayModel(Dense(2, 2)), SegmentedMean(2), identity)))
ProductModel … ↦ ArrayModel(identity)
  ├── ArrayModel(identity)
  └── BagModel … ↦ SegmentedMean(2) ↦ ArrayModel(identity)
        └── ArrayModel(Dense(2, 2))

julia> ProductModel(identity)
ProductModel … ↦ ArrayModel(identity)
  └── ArrayModel(identity)
```

See also: [`AbstractMillModel`](@ref), [`AbstractProductNode`](@ref), [`ProductNode`](@ref).
"""
function ProductModel(ms::VecOrTupOrNTup{Union{MillFunction, AbstractMillModel}},
                                m::Union{MillFunction, ArrayModel}=identity_model())
    ProductModel(map(_make_array_model, ms), _make_array_model(m))
end
ProductModel(ms::Union{MillFunction, AbstractMillModel},
             m::Union{MillFunction, ArrayModel}=identity_model()) = ProductModel((ms,), m)

Base.getindex(m::ProductModel, i::Symbol) = m.ms[i]
Base.keys(m::ProductModel) = keys(m.ms)

function (m::ProductModel{<:Tuple})(x::ProductNode{<:Tuple})
    m.m(vcat(map((sm, sx) -> sm(sx), m.ms, getfield(x, :data))...))
end
# function (m::ProductModel{<:NamedTuple})(x::ProductNode{<:NamedTuple})
#     m.m(vcat(map((sm, sx) -> sm(sx), m.ms, getfield(x, :data))...))
# end

function (m::ProductModel{<:NamedTuple})(x::ProductNode{<:NamedTuple})
    ms = getfield(m, :ms)
    cm = getfield(m, :m)
    cm(vcat(map((sm, sx) -> sm(sx), ms, getfield(x, :data))...))
end


text:
"""
Contains abstractions for describing orbital states and bodies.
Implementations are provided in TwoBody, and NBody.
"""
module AstrodynamicsCore

export  AbstractBody, AbstractOrbitalSystem, AbstractUnitfulStructure, AbstractState, AbstractSystem, AbstractOrbit, AbstractTrajectory, CartesianState,
        getindex, setindex!, lengthunit, timeunit, velocityunit, massparameterunit, coordinateframe,
        position_vector, velocity_vector, scalar_position, scalar_velocity,
        AbstractFrame, Bodycentric, Synodic, Perifocal, size, length, convert, epoch, eltype

export NormalizedLengthUnit, NormalizedTimeUnit

using Reexport
@reexport using Unitful, UnitfulAngles, UnitfulAstro
using StaticArrays: StaticVector, MVector
using LinearAlgebra: norm

include("../Misc/DocStringExtensions.jl")
include("../Misc/UnitfulAliases.jl")

"""
Abstract type for bodies in space: both `CelestialBody`s (in
`TwoBody.jl`), and `Body`s (in `NBody.jl`).
"""
abstract type AbstractBody end

abstract type AbstractOrbitalSystem end

"""
Abstract type describing a collection of states resulting from numerical integration
"""
abstract type AbstractTrajectory end

"""
Absract type describing all frames.
"""
abstract type AbstractFrame end

abstract type AbstractUnitfulStructure{F <: AbstractFloat, LU <: Unitful.LengthUnits, TU <: Unitful.TimeUnits} end

abstract type AbstractState{F, LU, TU, FR} <: AbstractUnitfulStructure{F, LU, TU} where FR <: AbstractFrame end

abstract type AbstractSystem{F, LU, TU} <: AbstractUnitfulStructure{F, LU, TU} end

abstract type AbstractOrbit{F, LU, TU} <: AbstractUnitfulStructure{F, LU, TU} end

"""
Body centered frame.
"""
struct Bodycentric <: AbstractFrame end

"""
Synodic frame.
"""
struct Synodic <: AbstractFrame end

"""
Perifocal frame.
"""
struct Perifocal <: AbstractFrame end

"""
Cartesian state which describes a spacecraft or body's position and velocity with respect to _something_.
"""
mutable struct CartesianState{F, LU, TU, FR} <: AbstractState{F, LU, TU, FR}
    t::F
    r::SubArray{F, 1, MVector{6, F}, Tuple{UnitRange{Int64}}, true}
    v::SubArray{F, 1, MVector{6, F}, Tuple{UnitRange{Int64}}, true}
    rv::MVector{6,F}

    function CartesianState(r::AbstractVector{R}, v::AbstractVector{V}, epoch::E=0, frame=Bodycentric; lengthunit=u"km", timeunit=u"s") where {R<:Real, V<:Real, E<:Real}
        F  = promote_type(R, V, E)
        if !(F <: AbstractFloat)
            @warn "Type provided ($(string(F))) is not a float: defaulting to Float64."
            F = Float64
        end
        @assert length(r) == length(v) == 3 "Both arguments must have length equal to 3!"
        rv = MVector{6, F}(r[1], r[2], r[3], v[1], v[2], v[3])
        pos = @views rv[1:3]
        vel = @views rv[4:6]
        return new{F, typeof(lengthunit), typeof(timeunit), frame}(F(epoch), pos, vel, rv)
    end

    function CartesianState(r::AbstractVector{R}, v::AbstractVector{V}, epoch::Unitful.Time = 0u"s", frame=Bodycentric) where {R<:Unitful.Length, V<:Unitful.Velocity}

        # Get ready to commit a crime... we need the Time unit provided in velocity quantity V

        # V is a quantity with units Length / Time.  Let's reverse that so we have Time instead of Time^(-1)
        TL = inv(unit(V))

        # Now we need to find which index (1 or 2) contains the Time unit
        timeaxis = findfirst(T -> T isa Unitful.Dimension{:Time}, collect(typeof(typeof(TL).parameters[2]).parameters[1]))

        # Now that we have the proper index, let's select the time unit
        timeunit = Unitful.FreeUnits{(typeof(TL).parameters[1][timeaxis],), Unitful.𝐓, nothing}()

        # This is easy...
        lengthunit = unit(R)

        # Phew!
        return CartesianState(ustrip.(lengthunit, r), ustrip.(lengthunit/timeunit, v), ustrip(timeunit, epoch), frame; lengthunit = lengthunit, timeunit = timeunit)
    end

    function CartesianState(cart::CartesianState{F,LU,TU,FR}) where {F,LU,TU,FR}
        return CartesianState(cart.r,  cart.v, cart.t, FR; lengthunit = LU(), timeunit = TU())
    end

    function CartesianState(arr::StaticVector{6})
        return CartesianState(arr[1:3], arr[4:6])
    end
end

Base.getindex(cart::CartesianState, i::Int) = cart.rv[i]
Base.setindex!(cart::CartesianState, value, i::Int) = (cart.rv[i] = value)
Base.length(::CartesianState) = 6
Base.size(::CartesianState) = (6,)

function Base.convert(::Type{CartesianState{F,LU,TU}}, cart::CartesianState) where {F,LU,TU}
    r = F.(ustrip.(LU(), position_vector(cart)))
    v = F.(ustrip.(LU()/TU(), velocity_vector(cart)))
    t = F.(ustrip.(TU(), epoch(cart)))
    return CartesianState(r, v, t, coordinateframe(cart); lengthunit = LU(), timeunit = TU())
end

epoch(cart::CartesianState) = cart.t * timeunit(cart)

const NormalizedLengthUnit   = Unitful.FreeUnits{(), Unitful.𝐋, nothing}
const NormalizedTimeUnit     = Unitful.FreeUnits{(), Unitful.𝐓, nothing}

Base.eltype(::AbstractUnitfulStructure{F}) where F = F

"""
Returns the `Unitful.Length` unit associated with the state.
"""
lengthunit(::AbstractUnitfulStructure{F, LU, TU}) where {F,LU,TU} = LU()

"""
Returns the `Unitful.Time` unit associated with the state.
"""
timeunit(::AbstractUnitfulStructure{F, LU, TU}) where {F, LU, TU} = TU()

"""
Returns the `Unitful.Velocity` unit associated with the state.
"""
velocityunit(state::AbstractUnitfulStructure) = lengthunit(state) / timeunit(state)

"""
Returns the `MassParameter` unit associated with the state.
"""
massparameterunit(state::AbstractUnitfulStructure) = lengthunit(state)^3 / timeunit(state)^2

"""
Returns the coordinate frame.
"""
coordinateframe(::AbstractState{F, LU, TU, FR}) where {F, LU, TU, FR} = FR

function Base.show(io::IO, state::CartesianState{F,LU,TU,FR}) where {F,LU,TU,FR}
    println(io, "  ", string(FR), " Cartesian State:")
    println("")
    println(io, "    t:  ", state.t, " ", string(TU()))
    println(io, "    r = ", [state.r[1] state.r[2] state.r[3]], " ", string(LU()))
    println(io, "    v = ", [state.v[1] state.v[2] state.v[3]], " ", string(LU()/TU()))
end

Base.show(io::IO, ::MIME"text/plain", state::CartesianState{F,LU,TU,FR}) where {F,LU,TU,FR} = show(io, state)

"""
Returns the `Unitful` position vector of the `Cartesianstate`.
"""
position_vector(state::CartesianState{F,LU,TU}) where {F, LU, TU} = state.r * lengthunit(state)

"""
Returns the `Unitful` velocity vector of the `CartesianState`.
"""
velocity_vector(state::CartesianState{F,LU,TU}) where {F, LU, TU} = state.v * velocityunit(state)

"""
Returns the `Unitful` scalar position of the `CartesianState`.
"""
scalar_position(state::CartesianState{F,LU,TU}) where {F, LU, TU} = norm(position_vector(state))

"""
Returns the `Unitful` scalar velocity of the `CartesianState`.
"""
scalar_velocity(state::CartesianState{F,LU,TU}) where {F, LU, TU} = norm(velocity_vector(state))

end


text:
module AccuracyAtTopKernels

# -------------------------------------------------------------------------------
# Used packages
# -------------------------------------------------------------------------------
using Statistics, LinearAlgebra, Random, Parameters
import Convex, Roots, Mmap, ProgressMeter
import Base: convert, show

import ECOS
import ECOS: Optimizer

import Flux.Optimise
import Flux.Optimise: Descent,
                      ADAM,
                      Momentum,
                      Nesterov,
                      RMSProp,
                      ADAGrad,
                      AdaMax,
                      ADADelta,
                      AMSGrad,
                      NADAM,
                      RADAM

import MLKernels
import MLKernels: Orientation,
                  Kernel,
                  ExponentialKernel,
                  LaplacianKernel,
                  SquaredExponentialKernel,
                  GaussianKernel,
                  RadialBasisKernel,
                  GammaExponentialKernel,
                  RationalQuadraticKernel,
                  MaternKernel,
                  LinearKernel,
                  PolynomialKernel,
                  ExponentiatedKernel,
                  PeriodicKernel,
                  PowerKernel,
                  LogKernel,
                  SigmoidKernel

# -------------------------------------------------------------------------------
# Export
# -------------------------------------------------------------------------------
export
    solve,
    scores,
    predict,
    exact_threshold,

    # Solvers
    AbstractSolver,
        General,
        Gradient,
        Coordinate,
    Optimizer,

    # Surrogates
    AbstractSurrogate,
        Hinge,
        Quadratic,
        Exponential,

    # Models
    AbstractModel,
        AbstractPatMat,
          PatMat,
          PatMatNP,
        AbstractTopPushK,
            TopPushK,
            TopPush,
            τFPL,

    # Data
    AbstractData,
        Primal,
        Dual,
            DTrain,
            DValidation,
            DTest,

    # Gradient descent optimizers (reexport Flux.Optimise)
    AbstractOptimizer,
        Descent,
        ADAM,
        Momentum,
        Nesterov,
        RMSProp,
        ADAGrad,
        AdaMax,
        ADADelta,
        AMSGrad,
        NADAM,
        RADAM,

    # Kernels (reexport MLKernels)
    Kernel,
    ExponentialKernel,
    LaplacianKernel,
    SquaredExponentialKernel,
    GaussianKernel,
    RadialBasisKernel,
    GammaExponentialKernel,
    RationalQuadraticKernel,
    MaternKernel,
    LinearKernel,
    PolynomialKernel,
    ExponentiatedKernel,
    PeriodicKernel,
    PowerKernel,
    LogKernel,
    SigmoidKernel

# -------------------------------------------------------------------------------
# types
# -------------------------------------------------------------------------------
abstract type AbstractSurrogate end
abstract type AbstractData end
abstract type AbstractModel end
abstract type AbstractPatMat{AbstractSurrogate}   <: AbstractModel end
abstract type AbstractTopPushK{AbstractSurrogate} <: AbstractModel end
abstract type AbstractSolver end

AbstractOptimizer = Union{Descent, ADAM, Momentum, Nesterov, RMSProp, ADAGrad,
                          AdaMax, ADADelta, AMSGrad, NADAM, RADAM}

function show(io::IO, opt::O) where {O<:AbstractOptimizer}
    args   = [getfield(opt, field) for field in fieldnames(O) if !in(field, (:velocity, :acc, :state))]
    print(io, "$(nameof(O))($(join(args, ",")))")
end

@with_kw_noshow struct General{I<:Integer,S} <: AbstractSolver
    seed::I   = rand(1:10000)
    solver::S = Optimizer(verbose = false)
end

show(io::IO, solver::General) =
    print(io, "General($(nameof(typeof(solver.solver))))")

struct Gradient{O<:AbstractOptimizer} <: AbstractSolver
    seed::Int
    maxiter::Int
    verbose::Bool
    optimizer::O
    iters::Vector{Int}

    function Gradient(;
        seed = rand(1:10000),
        maxiter = 1000,
        verbose = true,
        optimizer::O = ADAM(),
        iters = Int[],
    ) where {O<:AbstractOptimizer}

        iters = sort(unique(vcat(iters, maxiter)))
        return new{O}(seed, maxiter, verbose, optimizer, iters)
    end
end

show(io::IO, solver::Gradient) =
    print(io, "Gradient($(solver.maxiter), $(solver.optimizer))")

struct Coordinate <: AbstractSolver
    seed::Int
    maxiter::Int
    verbose::Bool
    iters::Vector{Int}

    function Coordinate(;
        seed = rand(1:10000),
        maxiter = 1000,
        verbose = true,
        iters = Int[],
    )

        iters = sort(unique(vcat(iters, maxiter)))
        return new(seed, maxiter, verbose, iters)
    end
end

show(io::IO, solver::Coordinate) =
    print(io, "Coordinate($(solver.maxiter))")

mutable struct BestUpdate{I<:Integer, T<:Real}
    k::I
    l::I
    Δ::T
    L::T
    vars::NamedTuple
end

# includes
include("surrogates.jl")
include("dataset.jl")

include("PatMat.jl")
include("TopPushK.jl")

include("utilities.jl")
include("solver.jl")
include("projections.jl")
include("kernels.jl")
include("scores_predict.jl")

end # module


text:
export tmapreduce, maptreduce, tmaptreduce;

"""
    tmapreduce(f::Function, op::Function, itr::AbstractArray; init)

Multi-threaded version of [mapreduce(f, op, itr; init)](https://docs.julialang.org/en/v1.2/base/collections/#Base.mapreduce-Tuple{Any,Any,Any}). Only the **mapping** part is parallelised.
"""
function tmapreduce(f::Function, op::Function, itr...; init)
	return reduce(op, tmap(f, itr...); init = init);
end

"""
    tmapreduce(f::Function, op::Function, itr::AbstractArray)

Multi-threaded version of [mapreduce(f, op, itr)](https://docs.julialang.org/en/v1.2/base/collections/#Base.mapreduce-Tuple{Any,Any,Any}). Only the **mapping** part is parallelised.
"""
function tmapreduce(f::Function, op::Function, itr...)
	return reduce(op, tmap(f, itr...));
end

"""
    maptreduce(f::Function, op::Function, itr::AbstractArray; init)

Multi-threaded version of [mapreduce(f, op, itr; init)](https://docs.julialang.org/en/v1.2/base/collections/#Base.mapreduce-Tuple{Any,Any,Any}). Only the **reduction** part is parallelised. Note that function `op` must **not** change type!
"""
function maptreduce(f::Function, op::Function, itr...; init)
	return treduce(op, map(f, itr...); init = init);
end

"""
    maptreduce(f::Function, op::Function, itr::AbstractArray)

Multi-threaded version of [mapreduce(f, op, itr)](https://docs.julialang.org/en/v1.2/base/collections/#Base.mapreduce-Tuple{Any,Any,Any}). Only the **reduction** part is parallelised. Note that function `op` must **not** change type!
"""
function maptreduce(f::Function, op::Function, itr...)
	return treduce(op, map(f, itr...));
end

"""
    tmaptreduce(f::Function, op::Function, itr::AbstractArray; init)

Multi-threaded version of [mapreduce(f, op, itr; init)](https://docs.julialang.org/en/v1.2/base/collections/#Base.mapreduce-Tuple{Any,Any,Any}). Both the **mapping** and the **reduction** part is parallelised. Note that function `op` must **not** change type!
"""
function tmaptreduce(f::Function, op::Function, itr...; init)
	return treduce(op, tmap(f, itr...); init = init);
end

"""
    tmaptreduce(f::Function, op::Function, itr::AbstractArray)

Multi-threaded version of [mapreduce(f, op, itr)](https://docs.julialang.org/en/v1.2/base/collections/#Base.mapreduce-Tuple{Any,Any,Any}). Both the **mapping** and the **reduction** part is parallelised. Note that function `op` must **not** change type!
"""
function tmaptreduce(f::Function, op::Function, itr...)
	return treduce(op, tmap(f, itr...));
end


text:
module TestSetIndex

include("preamble.jl")

@testset begin
    @test setindex!!((1, 2, 3), :two, 2) === (1, :two, 3)
    @test setindex!!((a=1, b=2, c=3), :two, :b) === (a=1, b=:two, c=3)
    @test setindex!!([1, 2, 3], :two, 2) == [1, :two, 3]
    @test setindex!!(Dict{Symbol,Int}(:a=>1, :b=>2), 10, :a) ==
        Dict(:a=>10, :b=>2)
    @test setindex!!(Dict{Symbol,Int}(:a=>1, :b=>2), 3, "c") ==
        Dict(:a=>1, :b=>2, "c"=>3)
end

@testset "mutation" begin
    @testset for args in [
        ([1, 2, 3], 20, 2),
        (Dict(:a=>1, :b=>2), 10, :a),
    ]
        @test setindex!!(args...) === args[1]
    end
end

end  # module


text:
"""
    φ , λ = yx2latlon(φ₀, λ₀, N, E)

Coordinate transformation between geographic and plane rectangular coordinates on the Gauss-Krüger Projection
This function is based on the following documents:
        https://www.gsi.go.jp/common/000061216.pdf
        https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/algorithm/xy2bl/xy2bl.htm
        https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/algorithm/bl2xy/bl2xy.htm

### input arguments
- φ₀, λ₀ : latitude of origin and central meridian
- N, E   : plane rectangular coordinates northing and easting
### return values
- φ , λ  : geographic latitude and longitude

### Examples
```julia-repl
julia> # ANSWER: lat = 36.10404755, lon = 140.08539843

julia> yx2latlon(36.0, 139.83333333, 11573.375, 22694.980)
(36.104047552508895, 140.08539842726532)
```

"""
function yx2latlon(φ₀, λ₀, N, E)

    S̅φ₀ = arclength(deg2rad(φ₀))

    # ξ and η
    ξ = (N + S̅φ₀)/A̅
    η = E/A̅

    # β terms
    βcoef = [1/2  -2/3   37/96       -1/360      -81/512; # β₁
             0.0  1/48    1/15    -437/1440       46/105; # β₂
             0.0   0.0  17/480      -37/840    -209/4480; # β₃
             0.0   0.0     0.0  4397/161280      -11/504; # β₄
             0.0   0.0     0.0          0.0  4583/161280; # β₅
            ]
    βj = βcoef * [n^k for k=1:5]

    ξ′ = ξ - sum([βj[j]sin(2j*ξ)cosh(2j*η) for j=1:5])
    η′ = η - sum([βj[j]cos(2j*ξ)sinh(2j*η) for j=1:5])

    # χ
    χ = asin(sin(ξ′)/cosh(η′))

    # δ terms
    δcoef = [2.0 -2/3  -2.0    116/45     26/45     -2854/675; # δ₁
             0.0  7/3  -8/5   -227/45  2704/315      2323/945; # δ₂
             0.0  0.0 56/15   -136/35 -1262/105    73814/2835; # δ₃
             0.0  0.0   0.0  4279/630   -332/35 -339572/14175; # δ₄
             0.0  0.0   0.0       0.0  4174/315  -144838/6237; # δ₅
             0.0  0.0   0.0       0.0       0.0  601676/22275; # δ₆
            ]
    δj = δcoef * [n^k for k=1:6]

    φ = rad2deg( χ + sum([δj[j]sin(2j*χ) for j=1:6]) )
    λ = λ₀ + rad2deg( atan(sinh(η′),cos(ξ′)) )

    #=
    σ′ = 1 - sum([2j*βj[j]*cos(2j*ξ)*cosh(2j*η) for j=1:5])
    τ′ =     sum([2j*βj[j]*sin(2j*ξ)*sinh(2j*η) for j=1:5])
    γ = atan(τ′+σ′*tan(ξ′)tanh(η′), σ′-τ′*tan(ξ′)tanh(η′))
    m = m₀*A₀*(1/(1+n))*sqrt( (cos(ξ′)^2+sinh(η′)^2)/(σ′^2+τ′^2) * (1+((1-n)tan(φ)/(1+n))^2) )
    =#

    return [φ λ]
end
# -------------------------------------
"""
    λ, φ = xy2lonlat(λ₀, φ₀, E, N)

See [`yx2latlon`](@ref)
"""
xy2lonlat(λ₀, φ₀, E, N) = reverse(yx2latlon(φ₀, λ₀, N, E), dims=2)
# -------------------------------------
function yx2latlon_ja(zone::Integer, N, E)
    !(0 < zone < 20) && error("zone ID must be a integer in the range of 1 to 19.")
    return yx2latlon(Origin_LatLon_Japan[zone,:]..., N, E)
end
# -------------------------------------
xy2lonlat_ja(zone::Integer, E, N) = reverse(yx2latlon_ja(zone, N, E), dims=2)
# -------------------------------------