text:
normal_kurtosis(x) = kurtosis(x, false)

function ks_test(a,b)
    x = ApproximateTwoSampleKSTest(a[2:end],b[2:end])
    return sqrt(x.n_x*x.n_y/(x.n_x+x.n_y))*x.δ
end

function generalized_hurts_exp(log_prices)
    q = 1
    max_τ = 19

    k = 0
    H = zeros(15)
    L = size(log_prices, 1)

    for iter_τ = 5:max_τ
        k += 1

        x = collect(1:iter_τ)
        k_q_t = zeros(iter_τ)
        for τ = 1:iter_τ
            numer = log_prices[(τ+1):τ:L] - log_prices[1:τ:(L-τ)]
            denom = log_prices[1:τ:L]

            # Determine Drift
            N = size(numer, 1) + 1
            X = collect(1:N)
            Y = denom
            mx = sum(X)/N
            SSxx = sum(X.^2) - N * mx^2
            my = sum(Y) / N
            SSxy = sum(X .* Y) - N * mx * my
            cc = [SSxy/SSxx, my - SSxy/SSxx*mx]

            # Subtract Drift
            numer = numer .- cc[1]
            denom = denom .- cc[1] * collect(1:N) .- cc[2]

            k_q_t[τ] = mean(abs.(numer).^q) / mean(abs.(denom).^q)
        end

        # Calculate Hurst Exponent for current iteration
        log_10_x = log10.(x)
        mx = mean(log_10_x)
        SSxx = sum(log_10_x.^2) - iter_τ * mx^2
        log_10_k_q_t = log10.(k_q_t)
        my = mean(log_10_k_q_t)
        SSxy = sum(log_10_x .* log_10_k_q_t) - iter_τ * mx * my
        H[k] = SSxy/SSxx
    end

    H = mean(H) / q
    return H
end

function get_summary_stats(simualted_log_prices, log_prices)

    num_replications = size(simualted_log_prices,2)
    stats_mat = zeros(num_replications,5)

    stats_mat[:,1] = mean(simualted_log_prices, dims=1)
    stats_mat[:,2] = std(simualted_log_prices, dims=1)
    stats_mat[:,3] = mapslices(normal_kurtosis, simualted_log_prices, dims=1)

    ks_test_stat(x) = ks_test(log_prices, x)
    stats_mat[:,4] = mapslices(ks_test_stat, simualted_log_prices, dims=1)
    stats_mat[:,5] = mapslices(generalized_hurts_exp, simualted_log_prices,
        dims=1)
    return stats_mat
end


text:
module JultraDark

using Base.Threads: @threads
using LinearAlgebra
using Statistics
using AbstractFFTs: fftfreq, rfftfreq
using PencilFFTs, MPI

using FFTW

export simulate
export Grids, PencilGrids
export Config, OutputConfig

include("grids.jl")
include("pencil_grids.jl")
include("output.jl")
include("config.jl")

"""
    max_time_step(grids, a)

Calculate an upper bound on the time step
"""
function max_time_step(grids, a)
    max_time_step_gravity = 2π / maximum(grids.Φx)
    max_time_step_pressure = 2π * 2 / maximum(grids.k)^2 * a^2  # TODO: cache k_max

    @assert isfinite(max_time_step_gravity)
    @assert isfinite(max_time_step_pressure)

    time_step = min(max_time_step_gravity, max_time_step_pressure)

    time_step
end

"""
    actual_time_step(max_timestep::Real, time_interval::Real, n::Integer)

Actual size and number of time steps that should be taken if the maximum
is `max_timestep`, no more than `n` steps should be taken, and they should
fit in `time_interval`.

# Examples

```jldoctest
julia> using JultraDark: actual_time_step

julia> actual_time_step(0.11, 1, 20)
(0.1, 10)
```
"""
function actual_time_step(max_timestep, time_interval, n)::Tuple{Real, Integer}
    if max_timestep * n > time_interval
        num_steps = ceil(time_interval / max_timestep)
        time_interval / num_steps, num_steps
    else
        max_timestep, n
    end
end

"""
Take `n` steps with time step `Δt`

# Examples

```jldoctest
julia> using JultraDark: take_steps!, PencilGrids, OutputConfig

julia> take_steps!(PencilGrids(1.0, 16), 0, 0.5, 10, OutputConfig(mktempdir(), []), t->1)
5.0

```
"""
function take_steps!(grids, t_start, Δt, n, output_config, a)

    t = t_start

    half_step = true

    for step in 1:n
        if half_step
            psi_half_step!(Δt, grids)
            t += Δt / 2
            half_step = false
        else
            psi_whole_step!(Δt, grids)
            t += Δt
        end

        phi_whole_step!(Δt, grids, a=a(t))

        output_summary_row(grids, output_config, t, a(t), Δt)
    end

    psi_half_step!(Δt, grids)
    t += Δt / 2

    output_summary_row(grids, output_config, t, a(t), Δt)

    t
end

"""
Evolve `grids` forward from `t_start` to `t_end`
"""
function evolve_to!(t_start, t_end, grids, output_config, config::Config.SimulationConfig)

    @assert t_start < t_end

    t = t_start

    while (t < t_end) && ~(t ≈ t_end)

        if max_phase_grad(grids.ψx) > π / 4
            output_grids(grids, output_config, -1)
            throw("Phase gradient is too large to continue")
        end

        Δt, n_steps = actual_time_step(
            max_time_step(grids, config.a(t)),
            t_end - t,
            config.time_step_update_period,
        )
        @debug "t = $t, max_time_step = $(max_time_step(grids, config.a(t))), Δt = $Δt"

        t = take_steps!(grids, t, Δt, n_steps, output_config, config.a)
    end

    t
end

function simulate(grids, options::Config.SimulationConfig, output_config::OutputConfig)

    # Setup output
    mkpath(output_config.directory)
    output_summary_header(output_config)

    t_begin = output_config.output_times[1]

    # Initialise vars other than ψx
    grids.ψk .= (grids.fft_plan * grids.ψx)
    grids.ρx .= abs2.(grids.ψx)
    grids.Φk .= -4 * π * (grids.rfft_plan * grids.ρx) ./ (options.a(t_begin) * grids.rk.^2)
    grids.Φk[1, 1, 1] = 0
    grids.Φx .= grids.rfft_plan \ grids.Φk

    output_grids(grids, output_config, 1)

    if max_phase_grad(grids.ψx) > π / 4
        throw("Phase gradient is too large to start")
    end

    for (index, t_end) in enumerate(output_config.output_times[2:end])
        t_begin = evolve_to!(
            t_begin,
            t_end,
            grids,
            output_config,
            options,
        )
        @info "Reached time $t_begin"
        output_grids(grids, output_config, index + 1)
    end

    if max_phase_grad(grids.ψx) > π / 4
        throw("Phase gradient is too large to end")
    end

end

end # module

ext: jl


text:
"""StanfordSentimentTreebank Type
having dict and senti fields for
phrases and corresponding sentiments respectively"""
struct StanfordSentimentTreebank
    filepath :: String
end

function init_datadeps(::Type{StanfordSentimentTreebank})
    register(DataDep(
        "Stanford Sentiment Treebank",
        """
        Stanford Sentiment Treebank V1.0

        It includes fine grained sentiment labels for 215,154 phrases in the parse trees of 11,855 sentences.

        title       :   Parsing With Compositional Vector Grammars
        author      :   Richard Socher and Alex Perelygin and Jean Wu and Jason Chuang and Christopher Manning and Andrew Ng and Christopher Potts
        booktitle   :   EMNLP
        year        :   2013
        Website     :   https://nlp.stanford.edu/sentiment/index.html

        Cite this paper if this dataset is used:
        https://nlp.stanford.edu/~socherr/EMNLP2013_RNTN.pdf
        """,
        "http://nlp.stanford.edu/~socherr/stanfordSentimentTreebank.zip",
        "3f5209483b46bbf129cacbbbe6ae02fe780407034f61cf6342b7833257c3f1db",
        post_fetch_method =function (fn)
            unpack(fn)
            dir = "stanfordSentimentTreebank"
            files = readdir(dir)
            mv.(joinpath.(dir, files), files)
            rm(dir), rm("__MACOSX", recursive=true)

            local dict
            open("dictionary.txt", "r") do fdict
                dict = split.(split(read(fdict, String), '\n'), '|')[1:end-1]
                dict = permutedims(reshape(hcat(dict...), (length( dict[1]), length(dict))))
            end

            local sentiments
            open("sentiment_labels.txt", "r") do fsenti
                sentiments = split.(split(read(fsenti, String), '\n'), '|')[2:end-1]
                sentiments = permutedims(reshape(hcat(sentiments...), (length(sentiments[1]), length(sentiments))))
            end

            open("sa_dataset.txt", "w") do fdataset
                for id=1:size(dict)[1]
                    line = join([dict[id, 1], sentiments[parse(Int, dict[id, 2])+1, 2]], '|')
                    write(fdataset, line, "\n")
                end
            end
        end
    ))
end

function StanfordSentimentTreebank()
    path = joinpath(datadep"Stanford Sentiment Treebank", "sa_dataset.txt")
    StanfordSentimentTreebank(path)
end

MultiResolutionIterators.levelname_map(::Type{StanfordSentimentTreebank}) = [
    :documents => 1,
    :sentences => 2,
    :words => 3, :tokens => 3,
    :characters => 4
]

function load(dataset::StanfordSentimentTreebank)
    data = map(eachline(dataset.filepath)) do line
        doc_raw, sentiment_raw = split(line, '|')
        sents = [intern.(tokenize(sent)) for sent in split_sentences(doc_raw)]
        sentiment = parse(Float64, sentiment_raw)
        [sents, sentiment]
    end
    permutedims(hcat(data...))
end


text:
################################################################################
# util.jl
#
# Utilities used in the Air library that don't depend on other components of
# Air.
#
# @author Noah C. Benson
#
# MIT License
# Copyright (c) 2020-2021 Noah C. Benson

# #Delay #######################################################################
struct _DelayPending
    _mux::ReentrantLock
    _fn::Function
end
"""
    Delay{T}

`Delay` objects can be used to lazily calculate a single value the first time
it is requested. They act like `Ref`s in that you access a delay `d` via `d[]`.
`Delay` objects are thread-safe and are functionally immutable.

See also: [`@delay`](@ref), [`LazyDict`](@ref)

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest
julia> # Create a Delay with a long run time.
       d = Delay{Int64}(() -> (println("Running."); sleep(2); 10))
Delay{Int64}(<...>)

julia> # Start a few threads, each of which attempt to read it. The function
       # will only run once.
       for th in [(Threads.@spawn d[]) for _ in 1:5]; wait(th) end
Running.

julia> # Ensure that it produced the correct value and doesn't run again.
       d[]
10
```
"""
mutable struct Delay{T} <: Base.Ref{T}
    _val::Union{_DelayPending, Some{T}}
    function Delay{T}(t) where {T}
        return new{T}(Some{T}(t))
    end
    # We want it to be possible to share a delayed value.
    function Delay{T}(d::Delay) where {T}
        dval = d._val
        if isa(dval, _DelayPending)
            return new{T}(_DelayPending(ReentrantLock(), () -> d[]))
        else
            return new{T}(Some{T}(dval.value))
        end
    end
    function Delay{T}(d::Delay{T}) where {T}
        return d
    end
    function Delay{T}(f::Function) where {T}
        return new{T}(_DelayPending(ReentrantLock(), f))
    end
end
Delay(t) = Delay{typeof(t)}(t)
Delay(d::Delay) = d
Delay(f::Function) = Delay{Any}(f)
"""
    @delay expression

Yields a `Delay` object that matches the given expression. The expression may be
one of the following:
 1. A function of no arguments, such as `() -> 10`; in this case the delay is made
    from this function (i.e., the RHS is the `->` expression that is delayed).
 2. A function with a set of symbol arguments evaluates the RHS but uses a let
    statement to wrap all the symbols in the LHS into a closure.
 3. An expression, which is treated as equivalent to `() -> expression`.
Optionally, the expression or LHS may be tagged with a type T. In this case, a
`Delay{T}` object is yielded instead of a `Delay{Any}`.

See also: [`Delay`](@ref), [`LazyDict`](@ref)

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest
julia> # Create a Delay with a long run time.
       d = (@delay (println("Running."); sleep(2); 10)::Int64)
Delay{Int64}(<...>)

julia> # Start a few threads, each of which attempt to read it. The function will
       # only run once.
       for th in [(Threads.@spawn d[]) for _ in 1:5]; wait(th) end
Running.

julia> # Ensure that it produced the correct value and doesn't run again.
       d[]
10

julia> # The display now shows the realized value also.
       d
Delay{Int64}(10)

julia> # Create a Delay with a locally-bound symbol.
       d2 = (@delay (d) -> (d[] / 20.0)::Float64)
Delay{Float64}(<...>)

julia> # We can rebind d without affecting d2.
       d = 10
10

julia> d2[]
0.5
```
"""
macro delay(e::Expr)
    # First, parse the expression-- is it a function or an expression?
    if e.head == :->
        syms = e.args[1]
        expr = e.args[2]
    else
        syms = ()
        expr = e
    end
    tmp = expr
    while tmp.head == :block; tmp = tmp.args[2] end
    if tmp.head == :(::)
        T = tmp.args[2]
    else
        T = :Any
    end
    # If there are symbols, convert them over to expressions
    isa(syms, Symbol) && (syms = :(($syms,)))
    if isa(syms, Expr)
        syms = [
            ( isa(sym, Symbol) ? :($sym = $sym)
              : sym.head == :(::) ? :($(sym.args[1]) = $sym)
              : throw(ArgumentError("closures must be symbols or tagged-symbols")))
            for sym in syms.args]
    end
    # Generate the code
    if length(syms) == 0
        return esc(:(Air.Delay{$T}(() -> $expr)))
    else
        return esc(:(let $(syms...); Air.Delay{$T}(() -> $expr) end))
    end
end
# Some base functions.
Base.getindex(d::Delay{T}) where {T} = begin
    v = d._val
    if isa(v, _DelayPending)
        lock(v._mux)
        try
            if d._val === v
                u = (v._fn())::T
                d._val = Some{T}(u)
            else
                u = ((d._val)::Some{T}).value
            end
            return u
        finally
            unlock(v._mux)
        end
    else
        return v.value
    end
end
Base.isready(d::Delay{T}) where {T} = isa(d._val, Some{T})
Base.setindex!(d::Delay{T}, x...) where {T} = throw(ArgumentError("setindex!: Delays are immutable"))
Base.isequal(a::Delay{T}, b::Delay{S}) where {T,S} = (a === b) || isequal(a[], b[])
Base.hash(d::Delay{T}) where {T} = 0x26c850a2957fa577 + hash(d[])
Base.show(io::IO, ::MIME"text/plain", d::Delay{T}) where {T} = begin
    v = d._val
    if isa(v, Some)
        print(io, "$(typeof(d))($(v.value))")
    else
        print(io, "$(typeof(d))(<...>)")
    end
end
export Delay, @delay

# #memoize #####################################################################
# First, this is a helper functionthat makes sure than a function-arg
# declaration has a name.
_memoize_fixarg(arg::Expr) = (arg.head == :(::) && length(arg.args) == 1
                              ? Expr(:(::), gensym(), arg.args[1])
                              : arg)
# Now the memoize macro itself.
"""
    @memoize name(args...) = expr
    @memoize name(args...) where {...} = expr

`@memoize` is a macro for declaring that the function declaration that follows
should be memoized in a private dictionary and any pre-calculated value should
be returned from that dictionary instead of being recalculated. All memoization
is thread-safe: `expr` is only ever evaluated by one thread at a time, and is
only ever evaluated once per unique set of arguments.

Note that arguments are memoized according to equality, so the use of mutable
arguments can result in undefined behavior of those arguments are later changed.

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest
julia> @memoize fib(n::Int) = begin
           println("Calculating fib(\$n)...")
           if n < 1
               return 0
           elseif n == 1
               return 1
           else
               return fib(n-1) + fib(n - 2)
           end
       end::Int
fib (generic function with 1 method)

julia> fib(5)
Calculating fib(5)...
Calculating fib(4)...
Calculating fib(3)...
Calculating fib(2)...
Calculating fib(1)...
Calculating fib(0)...
5

julia> fib(6)
Calculating fib(6)...
8
```
"""
macro memoize(assgn::Expr)
    (assgn.head == :(=)) || throw(
        ArgumentError("memconst must be given an assignment expression"))
    # Parse the assignment statement.
    lhs  = assgn.args[1]
    expr = assgn.args[2]
    if lhs.head == :call
        fsym = lhs.args[1]
        args = [_memoize_fixarg(a) for a in lhs.args[2:end]]
        lhs = Expr(:call, fsym, args...)
    elseif lhs.head == :where
        fsig = lhs.args[1]
        fsym = fsig.args[1]
        args = [_memoize_fixarg(a) for a in fsig.args[2:end]]
        fsig = Expr(:call, fsym, args...)
        lhs = Expr(:where, fsig, lhs.args[2:end]...)
    else
        ArgumentError(
            "memconst assignment LHS must be a call or where expression"
        ) |> throw
    end
    # See if the expr is tagged; if so, we have a particular type we can use in
    # the memoization dict.
    MT = expr.head === :(::) ? expr.head : :Any
    # Make an expression for the tuple of arguments.
    argtup = Expr(:tuple, args...)
    # Symbols we will need in the generated code.
    s_cache = gensym("cache")
    s_delay = gensym("delay")
    s_lock = gensym("lock")
    s_tmp  = gensym("tmp")
    s_val  = gensym("val")
    quote
        let $s_lock  = ReentrantLock(),
            $s_cache = Dict{Tuple, $MT}(),
            $s_delay = Dict{Tuple, Delay{$MT}}(),
            $s_tmp, $s_val
            global $lhs = begin
                lock($s_lock)
                try
                    $s_tmp = get($s_cache, $argtup, $s_cache)
                    if $s_tmp === $s_cache
                        $s_tmp = get!(() -> Delay{$MT}(() -> $expr),
                                      $s_delay, $argtup)
                    else
                        return $s_tmp
                    end
                finally
                    unlock($s_lock)
                end
                # If we get here, we've created or grabbed a delay for the
                # arguments; go ahead and wait on it (outside of the lock so
                # that we don't prevent other argument tuples from computing at
                # the same time).
                $s_val = $s_tmp[]
                # Now, re-grab the lock and update the dictionaries.
                lock($s_lock)
                try
                    # Possibly another thread updated things before we got to it.
                    $s_tmp = get($s_cache, $argtup, $s_cache)
                    if $s_tmp === $s_cache
                        # We're the first task to finish the calculation and/or
                        # the first to grab the lock. Fix the cache.
                        $s_cache[$argtup] = $s_val
                        delete!($s_delay, $argtup)
                        return $s_val
                    else
                        return $s_tmp
                    end
                finally
                    unlock($s_lock)
                end
            end
            $fsym
        end
    end |> esc
end
export @memoize

# #Promise #####################################################################
"""
    Promise{T}

Promise objects represent placeholders for values that may or may not have been
delivered yet. This is effectively a `Channel` object that can only be `put!` to
a single time and all `take` calls on the promise will return that value.
Promise values can be accessed via the `take()` function (not the `take!()`
function) or via `p[]` for promise `p`. In both cases, the running thread is
suspended until a value is delivered.

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest
julia> p = Promise{Symbol}()
Promise{Symbol}(<...>)

julia> isready(p)
false

julia> put!(p, :value)
:value

julia> isready(p)
true

julia> take(p)
:value

julia> p[]
:value
```
"""
mutable struct Promise{T}
    _val::Union{Threads.Condition, Some{T}}
end
Promise{T}() where {T} = Promise{T}(Threads.Condition())
Promise() = Promise{Any}()
# Some base functions.
"""
    take(promise)

Yields the value delivered to the given `Promise` object after suspending the
current thread to wait for the value if necessary.

See also [`Promise`](@ref), `Base.put!`.

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest
julia> p = Promise()
Promise{Any}(<...>)

julia> put!(p, :done)
:done

julia> p
Promise{Any}(:done)

julia> take(p)
:done
```
"""
take(d::Promise{T}) where {T} = begin
    v = d._val
    if isa(v, Threads.Condition)
        lock(v)
        try
            (d._val === v) && wait(v)
            return ((d._val)::Some{T}).value
        finally
            unlock(v)
        end
    else
        return v.value
    end
end
Base.getindex(d::Promise{T}) where {T} = take(d)
Base.put!(d::Promise{T}, x) where {T} = begin
    v = d._val
    if isa(v, Threads.Condition)
        lock(v)
        try
            (d._val === v) || throw(
                ArgumentError("given promise is already fulfilled"))
            d._val = Some{T}(x)
            notify(v, all=true)
            return d._val.value
        finally
            unlock(v)
        end
    else
        throw(ArgumentError("given promise is already fulfilled"))
    end
end
Base.isready(d::Promise{T}) where {T} = isa(d._val, Some{T})
Base.isequal(a::Promise{T}, b::Promise{S}) where {T,S} = (a === b) || isequal(a[], b[])
Base.hash(d::Promise{T}) where {T} = 0x767127451c1e402a + hash(d[])
Base.show(io::IO, ::MIME"text/plain", d::Promise{T}) where {T} = begin
    v = d._val
    if isa(v, Some)
        print(io, "$(typeof(d))(")
        show(io, v.value)
        print(io, ")")
    else
        print(io, "$(typeof(d))(<...>)")
    end
end
export Promise, take

# #lockall #####################################################################
_lockall(locks::Vector{T}) where {T} = begin
    locked = 0
    try
        for l in locks
            lock(l)
            locked += 1
        end
    finally
        (locked == length(locks)) || for l in locks
            unlock(l)
            locked -= 1
            (locked == 0) && break
        end
    end
    return locked
end
_lockall(f::Function, locks::Vector{T}) where {T} = begin
    # This only gets called once we have our own copy of the vector
    sort!(locks, by=objectid)
    # This will either lock them all or raise an exception.
    _lockall(locks)
    # Now we can run the function
    try
        return f()
    finally
        for l in locks
            unlock(l)
        end
    end
end
"""
    lockall(f, l1, l2, ...)

Locks all of the lockable objects `l1`, `l2`, etc. then runs `f`, unlocks the
objects, and returns the return value of `f()`.

See also: `Threads.ReentrantLock`

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest
julia> (r1, r2, r3) = [ReentrantLock() for _ in 1:3]
       lockall(r1, r2, r3) do; :success end
:success

julia> lockall([r1, r2, r3]) do; :success end
:success

julia> lockall((r1, r2, r3)) do; :success end
:success
```
"""
function lockall end
lockall(f::Function, locks::Vector{T}) where {T,N} = _lockall(f, copy(locks))
lockall(f::Function, locks::NTuple{N,T}) where {T,N} = _lockall(f, [locks...])
lockall(f::Function, locks::Vararg{T,N}) where {T,N} = _lockall(f, [locks...])
export lockall

# #_to_pairs ###################################################################
# A utility function for turning a list of pairs/tuples into a list of pairs.
_to_pairs(kvs) = begin
    if length(kvs) == 0
        K = Any
        V = Any
        if Base.IteratorEltype(kvs) isa Base.HasEltype
            ET = Base.eltype(itr)
            if isa(ET, DataType)
                if ET <: Pair
                    K = ET.parameters[1]
                    V = ET.parameters[2]
                elseif ET <: Tuple && length(ET.parameters) == 2
                    K = ET.parameters[1]
                    V = ET.parameters[2]
                end
            end
        end
        return Pair{K,V}[]
    else
        ks = []
        vs = []
        for kv in kvs
            if kv isa Pair || (kv isa Tuple && length(kv) == 2)
                push!(ks, kv[1])
                push!(vs, kv[2])
            else
                msg = "EquivDict: arg must be iterator of tuples or pairs"
                throw(ArgumentError(msg))
            end
        end
        K = typejoin(map(typeof, ks)...)
        V = typejoin(map(typeof, vs)...)
        return Pair{K,V}[Pair{K,V}(k,v) for (k,v) in zip(ks,vs)]
    end
end


text:
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.develop(path="$(homedir())/github/SimpleSDMLayers.jl")
using Revise
using SimpleSDMLayers
using Test

# Import QGIS-cropped layer as example
ref = geotiff(SimpleSDMPredictor, "../../rasters/wc1_cropped.tif")

## Step 2 : Fix getindex bounding coordinates

# Load things
layer = SimpleSDMPredictor(WorldClim, BioClim, 1)
layer_subset = layer[coords]
size(layer_subset) # should be 330, 570
left, right, bottom, top = coords

# Preliminary tests
layer_subset
ref
# Latitudes should be (20.083333333333332, 74.91666666666667)
# Longitudes should be (-144.91666666666666, -50.083333333333336)
isequal(latitudes(ref), latitudes(layer_subset)) # false
isapprox(latitudes(ref), latitudes(layer_subset)) # true

coords_subset = (layer_subset.left, layer_subset.right, layer_subset.bottom, layer_subset.top)
coords_ref = (ref.left, ref.right, ref.bottom, ref.top)
isequal.(coords_ref, coords_subset) # false
isapprox.(coords_ref, coords_subset) # true

## Investigate the bug
# 1. getindex(layer, coords)
imax = _match_longitude(layer, isnothing(right) ? layer.right : right; side=:right) # ok
imin = _match_longitude(layer, isnothing(left) ? layer.left : left; side=:left) # ok
jmax = _match_latitude(layer, isnothing(top) ? layer.top : top; side=:top) # ok
jmin = _match_latitude(layer, isnothing(bottom) ? layer.bottom : bottom; side=:bottom) # ok
layer[jmin:jmax, imin:imax] # problem is here, so in other getindex call
@which layer[jmin:jmax, imin:imax] # let's look at getindex(layer::T, i::R, j::R) where {T<:SimpleSDMLayer, R<:UnitRange}
# 2. getindex(layer, jmin:jmax, imin:imax)
T = typeof(layer)
i = jmin:jmax
j = imin:imax
i_min = isempty(i) ? max(i.start-1, 1) : i.start
i_max = isempty(i) ? max(i.stop+2, size(layer, 1)) : i.stop
j_min = isempty(j) ? max(j.start-1, 1) : j.start
j_max = isempty(j) ? max(j.stop+2, size(layer, 2)) : j.stop
i_fix = i_min:i_max
j_fix = j_min:j_max
RT = T <: SimpleSDMResponse ? SimpleSDMResponse : SimpleSDMPredictor
# RT(
    layer.grid[i_fix,j_fix] # ok
    minimum(longitudes(layer)[j_fix])-stride(layer,1) # ok
    maximum(longitudes(layer)[j_fix])+stride(layer,1) # nope
    minimum(latitudes(layer)[i_fix])-stride(layer,2) # nope
    maximum(latitudes(layer)[i_fix])+stride(layer,2) # ok
# )
# Compare allocation for range & LinRange
@time LinRange(layer.bottom+stride(layer, 2), layer.top-stride(layer, 2), size(layer,1))
@time range(layer.bottom+stride(layer, 2), layer.top-stride(layer, 2), length=size(layer,1))
# Let's just switch it to range instead of LinRange


text:
"""
    IS()

Importance sampling algorithm.

Note that this method is particle-based, and arrays of variables
must be stored in a [`TArray`](@ref) object.

Usage:

```julia
IS()
```

Example:

```julia
# Define a simple Normal model with unknown mean and variance.
@model gdemo(x) = begin
    s ~ InverseGamma(2,3)
    m ~ Normal(0,sqrt.(s))
    x[1] ~ Normal(m, sqrt.(s))
    x[2] ~ Normal(m, sqrt.(s))
    return s, m
end

sample(gdemo([1.5, 2]), IS(), 1000)
```
"""
struct IS{space} <: InferenceAlgorithm end

IS() = IS{()}()

DynamicPPL.initialsampler(sampler::Sampler{<:IS}) = sampler

function DynamicPPL.initialstep(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:IS},
    vi::AbstractVarInfo;
    kwargs...
)
    return Transition(vi), nothing
end

function AbstractMCMC.step(
    rng::Random.AbstractRNG,
    model::Model,
    spl::Sampler{<:IS},
    ::Nothing;
    kwargs...
)
    vi = VarInfo(rng, model, spl)
    return Transition(vi), nothing
end

# Calculate evidence.
function getlogevidence(samples::Vector{<:Transition}, ::Sampler{<:IS}, state)
    return logsumexp(map(x -> x.lp, samples)) - log(length(samples))
end

function DynamicPPL.assume(rng, spl::Sampler{<:IS}, dist::Distribution, vn::VarName, vi)
    if haskey(vi, vn)
        r = vi[vn]
    else
        r = rand(rng, dist)
        push!(vi, vn, r, dist, spl)
    end
    return r, 0
end

function DynamicPPL.observe(spl::Sampler{<:IS}, dist::Distribution, value, vi)
    return logpdf(dist, value)
end


text:
# Example 1
using Dates

# @time macro
#=
julia> @time sum(rand(10_000_000))
  0.043366 seconds (7 allocations: 76.294 MiB, 18.50% gc time)
5.000372959218527e6
=#

# manual version
#=
begin
    t1 = now()
    result = sum(rand(10_000_000))
    t2 = now()
    elapsed = t2 - t1
    println("It took ", elapsed)
    result
end
=#

# if implemented as a function
#=
function timeit(func)
    t1 = now()
    result = func()
    t2 = now()
    elapsed = t2 - t1
    println("It took ", elapsed)
    result
end
=#

# how to use?
#=
julia> mycode() = sum(rand(10_000_000))
mycode (generic function with 1 method)

julia> timeit(mycode)
It took 45 milliseconds
5.002266841783055e6
=#

#------------------------------------------------------------------------------
# loop unrolling

# sample loop
#=
for i in 1:3
    println("hello: ", i)
end
=#

# unrolled version
#=
println("hello: ", 1)
println("hello: ", 2)
println("hello: ", 3)
=#

# using the Unrolled.jl package
#=
using Unrolled

@unroll function hello(xs)
    @unroll for i in xs
        println("hello: ", i)
    end
end
=#

# usage
#=
julia> seq = tuple(1:3...)
(1, 2, 3)

julia> hello(seq)
hello: 1
hello: 2
hello: 3
=#

# lowered code
#=
julia> @code_lowered(hello(seq))
CodeInfo(
1 ─ i@_5 = (Base.getindex)(xs, 1)
│ (Main.println)("hello: ", i@_5)
│ i@_4 = (Base.getindex)(xs, 2)
│ (Main.println)("hello: ", i@_4)
│ i@_3 = (Base.getindex)(xs, 3)
│ (Main.println)("hello: ", i@_3)
│ %7 = Main.nothing
└── return %7
)
=#


text:
#
# Load DFTK-compatible structural information from an external file
# Relies on ASE and other external libraries to do the parsing
#
using PyCall

function load_lattice(T, file::AbstractString)
    ase = pyimport_e("ase")
    ispynull(ase) && error("Install ASE to load data from exteral files")
    load_lattice(T, pyimport("ase.io").read(file))
end

function load_atoms(T, file::AbstractString)
    ase = pyimport_e("ase")
    ispynull(ase) && error("Install ASE to load data from exteral files")
    load_atoms(T, pyimport("ase.io").read(file))
end

function load_magnetic_moments(file::AbstractString)
    ase = pyimport_e("ase")
    ispynull(ase) && error("Install ASE to load data from exteral files")
    load_magnetic_moments(pyimport("ase.io").read(file))
end

"""
Return a DFTK-compatible `lattice` object loaded from an ASE `Atoms`, a pymatgen `Structure`
or a compatible file (e.g. xyz file, cif file etc.)
"""
load_lattice(args...; kwargs...) = load_lattice(Float64, args...; kwargs...)

"""
Return a DFTK-compatible `atoms` object loaded from an ASE `Atoms`, a pymatgen `Structure`
or a compatible file (e.g. xyz file, cif file etc.)
"""
load_atoms(  args...; kwargs...) = load_atoms(  Float64, args...; kwargs...)

load_model(  args...; kwargs...) = load_model(  Float64, args...; kwargs...)
load_basis(  args...; kwargs...) = load_basis(  Float64, args...; kwargs...)
load_density(args...; kwargs...) = load_density(Float64, args...; kwargs...)


text:
function set_all_start_values(m)
    jm = m._jump_model
    for s in m._sectors
        if s isa ScalarSector
            Complementarity.set_start_value(jm[s.name], s.benchmark)
        else
            for i in Iterators.product(s.indices...)
                Complementarity.set_start_value(jm[s.name][i...], s.benchmark[i...])
            end
        end
    end

    for c in m._commodities
        if c isa ScalarCommodity
            Complementarity.set_start_value(jm[c.name], c.benchmark)
        else
            for i in Iterators.product(c.indices...)
                Complementarity.set_start_value(jm[c.name][i...], c.benchmark[i...])
            end
        end
    end

    for s in m._productions
        for i in s.inputs
            compensated_input1_demand_name = get_comp_demand_name(i)
            Complementarity.set_start_value(jm[compensated_input1_demand_name], eval(swap_our_param_with_val(i.quantity)))
        end
    end

    for c in m._consumers
        Complementarity.set_start_value(jm[c.name], c.benchmark)
    end

    # Add compensated supply variables
    for s in m._productions
        for o in s.outputs
            Complementarity.set_start_value(jm[get_comp_supply_name(o)], eval(swap_our_param_with_val(o.quantity)))
        end
    end

    # Add final demand variables
    for demand_function in m._demands
        for demand in demand_function.demands
            Complementarity.set_start_value(jm[get_final_demand_name(demand)], eval(swap_our_param_with_val(demand.quantity)))
        end
    end
end

function set_all_bounds(m)
    jm = m._jump_model

    for c in m._commodities
        if c isa ScalarCommodity
            jump_var = jm[c.name]

            if c.fixed
                JuMP.fix(jump_var, c.benchmark, force=true)
            else
                if JuMP.is_fixed(jump_var)
                    JuMP.unfix(jump_var)
                end

                JuMP.set_lower_bound(jump_var, 0.001)
            end
        else
            for i in Iterators.product(c.indices...)
                jump_var = jm[c.name][i...]

                if c.fixed[i...]
                    JuMP.fix(jump_var, c.benchmark[i...], force=true)
                else
                    if JuMP.is_fixed(jump_var)
                        JuMP.unfix(jump_var)
                    end

                    JuMP.set_lower_bound(jump_var, 0.001)
                end
            end
        end

    end
    for cs in m._consumers
        jump_var = jm[cs.name]
        if cs.fixed
            JuMP.fix(jump_var, cs.benchmark, force=true)
        else
            if JuMP.is_fixed(jump_var)
                JuMP.unfix(jump_var)
            end
        end
    end
end

function set_all_parameters(m)
    jm = m._jump_model

    for p in m._parameters
        if p isa ScalarParameter
            JuMP.set_value(jm[p.name], p.value)
        else
            for pp in p.indices
                JuMP.set_value(jm[pp.name], pp.value)
            end
        end
    end
end


text:
# Script to generate figures in Figure 4.
using Plots, SparseArrays, JLD2, LinearAlgebra, MultiscaleGraphSignalTransforms

# Set up the resolution and display size
gr(dpi=200, size=(800,600))

# Load the Toronto street network (vehicular volume counts as its graph signal)
JLD2.@load "Toronto_new.jld2"
tmp1 = toronto["G"]
G=GraphSig(tmp1["W"],xy=tmp1["xy"],f=tmp1["f"],name =tmp1["name"],
    plotspecs = tmp1["plotspecs"])

# Assign correct edge weights via 1/Euclidean distance
G = Adj2InvEuc(G)
# Perform the hierarchical graph partitioning using the Fiedler vectors
GP = partition_tree_fiedler(G) # Lrw by default
# Compute the expansion coefficients of the full GHWT c2f dictionary
dmatrix = ghwt_analysis!(G, GP=GP)
N = length(G.f)

# Fig. 4a: Sacling vector ψ^1_{0,0}, i.e., (j,k,l)=(1,0,0).
j = 1
loc = 1
BS = bs_level(GP, j)
dvec = zeros(N, 1)
dvec[loc,1] = 1.0
(f, GS) = ghwt_synthesis(dvec, GP, BS, G)
GraphSig_Plot(GS, linewidth = 1., markersize = 4., markerstrokealpha = 0.,
    markercolor = :viridis, notitle = true, nocolorbar = true, clim = (-0.01, 0.01))
plot!(showaxis = false, ticks = false)
k, l = BS.levlist[loc][1], BS.levlist[loc][2]
print(j," ",(rs_to_region(GP.rs, GP.tag))[k,l]," ",GP.tag[k,l])
savefig("toronto_psi100.pdf")
savefig("toronto_psi100.png")

# Fig. 4b: Haar vector ψ^2_{0,1}, i.e., (j,k,l)=(2,0,1).
j = 2
loc = 2
BS = bs_level(GP, j)
dvec = zeros(N, 1)
dvec[loc,1] = 1.0
(f, GS) = ghwt_synthesis(dvec, GP, BS, G)
GraphSig_Plot(GS, linewidth = 1., markersize = 4., markerstrokealpha = 0.,
    markercolor = :viridis, notitle = true, nocolorbar = true, clim = (-0.01, 0.01))
plot!(showaxis = false, ticks = false)
k, l = BS.levlist[loc][1], BS.levlist[loc][2]
print(j," ",(rs_to_region(GP.rs, GP.tag))[k,l]," ",GP.tag[k,l])
savefig("toronto_psi201.pdf")
savefig("toronto_psi201.png")

# Fig. 4c: Walsh vector ψ^3_{0,9}, i.e., (j,k,l)=(3,0,9).
#j = 5
#loc = 1000
j = 3
loc = 10
BS = bs_level(GP, j)
dvec = zeros(N, 1)
dvec[loc,1] = 1
(f, GS) = ghwt_synthesis(dvec, GP, BS, G)
GraphSig_Plot(GS, linewidth = 1., markersize = 4., markerstrokealpha = 0.,
    markercolor = :viridis, notitle = true, nocolorbar = true, clim = (-0.01, 0.01))
plot!(showaxis = false, ticks = false)
k, l = BS.levlist[loc][1], BS.levlist[loc][2]
print(j," ",(rs_to_region(GP.rs, GP.tag))[k,l]," ",GP.tag[k,l])
savefig("toronto_psi309.pdf")
savefig("toronto_psi309.png")


text:
# Autogenerated wrapper script for LibPQ_jll for armv6l-linux-musleabihf
export LIBPQ_HANDLE

using OpenSSL_jll
using Kerberos_krb5_jll
JLLWrappers.@generate_wrapper_header("LibPQ")
JLLWrappers.@declare_library_product(LIBPQ_HANDLE, "libpq.so.5")
function __init__()
    JLLWrappers.@generate_init_header(OpenSSL_jll, Kerberos_krb5_jll)
    JLLWrappers.@init_library_product(
        LIBPQ_HANDLE,
        "lib/libpq.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
create(c::Client, ::Type{DiscordChannel}, g::AbstractGuild; kwargs...) = create_guild_channel(c, g.id; kwargs...)
create(c::Client, ::Type{DiscordChannel}, u::User; kwargs...) = create_dm(c; kwargs..., recipient_id = u.id)

retrieve(c::Client, ::Type{DiscordChannel}, channel::Integer) = get_channel(c, channel)
retrieve(c::Client, ::Type{DiscordChannel}, g::AbstractGuild) = get_guild_channels(c, g.id)

update(c::Client, ch::DiscordChannel; kwargs...) = modify_channel(c, ch.id; kwargs...)

delete(c::Client, ch::DiscordChannel) = delete_channel(c, ch.id)


text:
#include("./computer.jl")
import AoC2019

mutable struct BoardState
  dims
  my_display::Array{Int8, 2}
  current_score
end

function update_board_state!(state, x, y, arg)
  if x == -1 && y == 0
      state.current_score = arg
  else
    state.my_display[y+1, x+1] = arg
  end
end

function move_controller(state)
  the_ball = findfirst(x -> (x == 4), state.my_display)
  the_pad = findfirst(x -> (x == 3), state.my_display)
  #println(the_ball)
  #println(the_pad)
  if the_ball === nothing || the_pad === nothing
    return 0
  end
  direction = AoC2019.pseudoNormalize(the_ball[2] - the_pad[2])
  #println(stderr, direction)
  return direction
end

function draw_state(state)
  color_map = Dict(0 => ' ', 1 => 'w', 2 => 'b', 3 => '-', 4 => 'o')
  println(AoC2019.visualizeArray(state.my_display, color_map))
end

function arcade_cabinet(input_data)

  game_input = Channel{Int64}(Inf)
  game_output = Channel{Int64}(Inf)
  #dims = (24, 42)
  dims = (60, 60)
  my_display = zeros(Int8, dims...)
  current_score = 0

  state = BoardState((50,50), my_display, current_score)
  input_data[1] = 2
  condition=Condition()

  @sync begin
    my_program = @async AoC2019.run_program_async(input_data, game_input, game_output)

    @async while true
      #println(stderr, "Reading X")

      x = take!(game_output)
      #println(stderr, "Reading Y")
      y = take!(game_output)
      #println(stderr, "Reading Arg")
      arg = take!(game_output)
      #println(stderr, "X: $x, Y: $y, arg: $arg")
      #println(stderr, "X == -1: $(x == -1)")
      update_board_state!(state, x, y, arg)
      #      #display(plot!(state.my_display))
    end
    @async while true
      try
        #wait(condition)
        #println(stderr, "woken up")
        draw_state(state)
        sleep(0.033)
        #println(stderr, "Broken")
        direction = move_controller(state)
        #println(stderr, "Inputting $direction")
        put!(game_input, direction)
        #println(stderr, "Input $direction")
        yield()
      catch
        break
      end
    end
  end

  println("Game finished")
  draw_state(state)

  return state.current_score
end

input_data = AoC2019.read_numbers("input13.txt")
println("Good luck")
out = arcade_cabinet(input_data)
println(out)

max_stars_repo_name: buotex/aoc

max_issues_repo_name: buotex/aoc

max_forks_repo_name: buotex/aoc


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function barfunc1(low, high, v; extremes=extrema(values), barnumber=0, bartotal=0, mvalue=0)
    @layer begin
        sethue(rescale(v, 0, extremes[2], 0, 1), rand(), rand())
        circle(low, 2, :fill)
        box(Point(low.x - 5, low.y), Point(high.x + 5, high.y), :fill)
        fontsize(6)
        sethue("black")
        textbox([string(round(v, digits=2))], Point(low.x, 100), leading=5, alignment=:center)
        sethue("gray80")
        setline(0.5)
        setdash("dot")
        line(low, Point(low.x, 100), :stroke)
    end
end

function barfunc2(low, high, value; extremes=extrema(values), barnumber=0, bartotal=0)
    @layer begin
        sethue("red")
        randomhue()
        setline(4)
        setlinecap("round")
        line(low, high, :stroke)
        circle(high, 2, :fill)
    end
end

emptylabelfunction(args...; extremes=[], barnumber=0, bartotal=0) = nothing

function mylabelfunction(bottom::Point, top::Point, value;
        extremes=[], barnumber=0, bartotal=0)
    t = string(value)
    textoffset = textextents(t)[4]
    if top.y < 0
        tp = Point(top.x, min(top.y, bottom.y) - textoffset)
    else
        tp = Point(top.x, max(top.y, bottom.y) + textoffset)
    end
    @layer begin
        translate(tp)
        rotate(-pi/2)
        text(t, O, halign=:left, valign=:middle)
    end
end

function test_bars(fname)
    pagewidth, pageheight = 1200, 1400
    Drawing(pagewidth, pageheight, fname)
    background("antiquewhite")
    fontsize(5)
    sethue("black")
    origin()
    translate(-200, 0)
    tiles = Tiler(pagewidth, pageheight, 6, 2, margin=35)
    for (pos, n) in tiles
        @layer begin
            translate(pos)
            if n % 5 == 1
                v = randn(15)
                bars(v, labelfunction = (args...; extremes=[], barnumber=0, bartotal=0) ->  setgray(rand()))
            elseif n % 5 == 2
                # no labels
                v = rand(1:10, 100)
                bars(v, xwidth=rand(2:6), yheight=rand(50:100), labelfunction=emptylabelfunction)
            elseif n % 5 == 3
                # a custom label function
                v = rand(20)
                bars(v, yheight=rand(10:80), xwidth=rand(5:20), labelfunction=mylabelfunction)
            elseif n % 5 == 4
                # a custom barfunction
                v = rand(5:0.1:10, 20)
                bars(v, yheight=30, xwidth=20, labels=false, barfunction = barfunc1)
            elseif n % 5 == 5
                # another custom barfunction
                v = rand(-20:2:20, 30)
                bars(v, xwidth=14, barfunction=barfunc2)
            else
                v = randn(15)
                bars(v)
            end
        end
    end
    @test finish() == true
    println("...finished barstest, saved in $(fname)")
end

fname = "bars-test.pdf"
test_bars(fname)


text:
using Metatheory

@testset "Reduction Basics" begin
  t = @theory begin
    ~a + ~a --> 2 * (~a)
    ~x / ~x --> 1
    ~x * 1 --> ~x
  end

  # basic theory to check that everything works
  @test rewrite(:(a + a), t) == :(2a)
  @test rewrite(:(a + (x * 1)), t) == :(a + x)
  @test rewrite(:(a + (a * 1)), t; order = :inner) == :(2a)
end

import Base.(+)
@testset "Extending Algebra Operators" begin
  t = @theory begin
    ~a + ~a --> 2(~a)
  end

  # Let's extend an operator from base, for sake of example
  function +(x::Symbol, y)
    rewrite(:($x + $y), t)
  end

  @test (:x + :x) == :(2x)
end

## Free Monoid

@testset "Free Monoid - Overriding identity" begin
  # support symbol literals
  symbol_monoid = @theory begin
    ~a ⋅ :ε --> ~a
    :ε ⋅ ~a --> ~a
    ~a::Symbol --> ~a
    ~a::Symbol ⋅ ~b::Symbol => Symbol(String(a) * String(b))
    # i |> error("unsupported ", i)
  end

  @test rewrite(:(ε ⋅ a ⋅ ε ⋅ b ⋅ c ⋅ (ε ⋅ ε ⋅ d) ⋅ e), symbol_monoid; order = :inner) == :abcde
end

## Interpolation should be possible at runtime

@testset "Calculator" begin
  calculator = @theory begin
    ~x::Number ⊕ ~y::Number => ~x + ~y
    ~x::Number ⊗ ~y::Number => ~x * ~y
    ~x::Number ⊖ ~y::Number => ~x ÷ ~y
    ~x::Symbol --> ~x
    ~x::Number --> ~x
  end
  a = 10

  @test rewrite(:(3 ⊕ 1 ⊕ $a), calculator; order = :inner) == 14
end

## Direct rules
@testset "Direct Rules" begin
  t = @theory begin
    # maps
    ~a * ~b => ((~a isa Number && ~b isa Number) ? ~a * ~b : _lhs_expr)
  end
  @test rewrite(:(3 * 1), t) == 3

  t = @theory begin
    # maps
    ~a::Number * ~b::Number => ~a * ~b
  end
  @test rewrite(:(3 * 1), t) == 3
end

## Take advantage of subtyping.
# Subtyping in Julia has been formalized in this paper
# [Julia Subtyping: A Rational Reconstruction](https://benchung.github.io/papers/jlsub.pdf)

abstract type Vehicle end
abstract type GroundVehicle <: Vehicle end
abstract type AirVehicle <: Vehicle end
struct Airplane <: AirVehicle end
struct Car <: GroundVehicle end

airpl = Airplane()
car = Car()

t = @theory begin
  ~a::AirVehicle * ~b => "flies"
  ~a::GroundVehicle * ~b => "doesnt_fly"
end

@testset "Subtyping" begin

  sf = rewrite(:($airpl * c), t)
  df = rewrite(:($car * c), t)

  @test sf == "flies"
  @test df == "doesnt_fly"
end

@testset "Interpolation" begin
  airpl = Airplane()
  car = Car()
  t = @theory begin
    airpl * ~b => "flies"
    car * ~b => "doesnt_fly"
  end

  sf = rewrite(:($airpl * c), t)
  df = rewrite(:($car * c), t)

  @test sf == "flies"
  @test df == "doesnt_fly"
end

@testset "Segment Variables" begin
  t = @theory begin
    f(~x, ~~y) => Expr(:call, :ok, (~~y)...)
  end
  sf = rewrite(:(f(1, 2, 3, 4)), t)
  @test sf == :(ok(2, 3, 4))

  t = @theory x y begin
    f(x, y...) => Expr(:call, :ok, y...)
  end
  sf = rewrite(:(f(1, 2, 3, 4)), t)
  @test sf == :(ok(2, 3, 4))

  t = @theory x y begin
    f(x, y...) --> ok(y...)
  end
  sf = rewrite(:(f(1, 2, 3, 4)), t)
  @test sf == :(ok(2, 3, 4))
end

module NonCall
using Metatheory
t = [@rule a b (a, b) --> ok(a, b)]

test() = rewrite(:(x, y), t)
end

@testset "Non-Call expressions" begin
  @test NonCall.test() == :(ok(x, y))
end

@testset "Pattern matcher can match on both function object references and name symbols" begin
  ex = :($(+)($(sin)(x)^2, $(cos)(x)^2))
  r = @rule(sin(~x)^2 + cos(~x)^2 --> 1)

  @test r(ex) == 1
end

@testset "Pattern variable as pattern term head" begin
  foo(x) = x + 2
  ex = :(($foo)(bar, 2, pazz))
  r = @rule ((~f)(~x, 2, ~y) => (~f)(2))

  @test r(ex) == 4
end

using TermInterface

using Metatheory.Syntax: @capture
@testset "Capture form" begin
  ex = :(a^a)

  #note that @test inserts a soft local scope (try-catch) that would gobble
  #the matches from assignment statements in @capture macro, so we call it
  #outside the test macro
  ret = @capture ex (~x)^(~x)
  @test ret
  @test @isdefined x
  @test x === :a

  ex = :(b^a)
  ret = @capture ex (~y)^(~y)
  @test !ret
  @test !(@isdefined y)

  ret = @capture :(a + b) (+)(~~z)
  @test ret
  @test @isdefined z
  @test all(z .=== arguments(:(a + b)))

  #a more typical way to use the @capture macro

  f(x) =
    if @capture x (~w)^(~w)
      w
    end

  @test f(:(b^b)) == :b
  @test isnothing(f(:(b + b)))

  x = 1
  r = (@capture x x)
  @test r == true
end

using TermInterface
@testset "Matchable struct" begin
  struct qux
    args
    qux(args...) = new(args)
  end
  TermInterface.operation(::qux) = qux
  TermInterface.istree(::qux) = true
  TermInterface.arguments(x::qux) = [x.args...]

  @capture qux(1, 2) qux(1, 2)

  @test (@rule qux(1, 2) => "hello")(qux(1, 2)) == "hello"
  @test (@rule qux(1, 2) => "hello")(1) === nothing
  @test (@rule 1 => "hello")(1) == "hello"
  @test (@rule 1 => "hello")(qux(1, 2)) === nothing
  @test (@capture qux(1, 2) qux(1, 2))
  @test false == (@capture qux(1, 2) qux(3, 4))
end


text:
const VTREE_FORMAT = """c ids of vtree nodes start at 0
c ids of variables start at 1
c vtree nodes appear bottom-up, children before parents
c
c file syntax:
c vtree number-of-nodes-in-vtree
c L id-of-leaf-vtree-node id-of-variable
c I id-of-internal-vtree-node id-of-left-child id-of-right-child
c
"""

"""
Used to specify file type .vtree or .dot
"""
abstract type VtreeAbstractFile end

mutable struct VtreeConfigFile{T<:AbstractString, W<:IOStream} <: VtreeAbstractFile
    name::T
    file::W
end

mutable struct VtreeDotFile{T<:AbstractString, W<:IOStream} <: VtreeAbstractFile
    name::T
    file::W
end

VtreeConfigFile(name::AbstractString) = VtreeConfigFile(name, open(name, "w"))
VtreeDotFile(name::AbstractString) = VtreeDotFile(name, open(name, "w"))

# reorder: parents before children, as graphviz draw the picture in the order it appears
function reverse_order(dotfile::VtreeDotFile)
    close(dotfile.file)
    f = open(dotfile.name, "r")
    lines = readlines(f)
    lines[2:end-1] = reverse(lines[2:end-1])
    close(f)
    dotfile.file = open(dotfile.name, "w")
    for line in lines
        write(dotfile.file, line * "\n")
    end
end

"""
Saves a vtree in the given file path.
"""
function save(vtree::PlainVtree, file::AbstractString)

    "1. decide file type and open file"
    if endswith(file,".vtree")
        f = VtreeConfigFile(file)
    elseif endswith(file, ".dot")
        f = VtreeDotFile(file)
    else
        throw("Invalid file type")
    end

    "2. map from PlainVtreeNode to index for output"
    index_cache = Dict{PlainVtreeNode, UInt32}()
    index = -1
    node2index(n::PlainVtreeNode) =
        get!(index_cache, n) do
            index += 1
        end

    "3. saving methods for header, nodes, tailer"
    function save_vtree_header(vtree::PlainVtree, f::VtreeConfigFile)
        vtree_count = length(vtree)
        write(f.file, VTREE_FORMAT)
        write(f.file, "vtree $vtree_count\n")
    end

    function save_vtree_header(vtree::PlainVtree, f::VtreeDotFile)
        write(f.file,"strict graph vtree { node [shape=point]; splines=false; \n")
    end

    function save_vtree_node(n::PlainVtreeLeafNode, f)
        node_index = node2index(n)
        node_variable = n.var

        if f isa VtreeConfigFile
            write(f.file, "L $node_index $node_variable\n")
        elseif f isa VtreeDotFile
            write(f.file, "$node_index [label=$(node_variable), shape=\"plaintext\"]\n")
        else
            @assert 0
        end
    end

    function save_vtree_node(n::PlainVtreeInnerNode, f)
        node_index = node2index(n)
        left = node2index(n.left)
        right = node2index(n.right)

        if f isa VtreeConfigFile
            write(f.file, "I $node_index $left $right\n")
        elseif f isa VtreeDotFile
            write(f.file, "$node_index -- $right\n")
            write(f.file, "$node_index -- $left\n")
        else
            @assert 0
        end
    end

    function save_vtree_tailer(f::VtreeConfigFile)
    end #do nothing

    function save_vtree_tailer(f::VtreeDotFile)
        write(f.file, "}\n")

        reverse_order(f)

    end

    " 4. saving frame"
    order = node2dag(vtree[end])

    save_vtree_header(vtree, f)

    for node in order
        save_vtree_node(node,f)
    end

    save_vtree_tailer(f)

    close(f.file)

end


text:
using Yama
using Documenter

makedocs(;
    modules=[Yama],
    authors="Anton Smirnov",
    repo="https://github.com/pxl-th/Yama.jl/blob/{commit}{path}#L{line}",
    sitename="Yama.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://pxl-th.github.io/Yama.jl",
        assets=String[],
    ),
    pages=[
        "index.md",
        "Examples" => [
            "Surface" => "examples/surface.md",
            "Trajectory" => "examples/trajectory.md",
        ],
        "User API" => "docs.md",
    ],
)
deploydocs(repo="github.com/pxl-th/Yama.jl")


text:
# Monoid instances
# ================

# this is standard Applicative combine implementation
# there is no other sensible definition for combine and hence if this does fail, it fails correctly
TypeClasses.combine(a::Task, b::Task) = mapn(combine, a, b)

"""
    orelse(task1, task2, ...)
    task1 ⊘ task2

Runs both in parallel and collects which ever result is first. Then interrupts all other tasks and returns the found result.
"""
TypeClasses.orelse(a::Task, b::Task, more::Task...) = @async begin
  n = length(more) + 2
  c = Channel(n)  # all tasks should in principle be able to send there result without waiting
  tasks = (a, b, more...)
  tasks′ = map(tasks) do task
    @async begin
      tried = @TryCatch TaskFailedException fetch(task)
      # we need to capture errors too, so that in case everything errors out, this is still ending
      put!(c, tried[])
    end
  end

  exceptions = Vector{Exception}(undef, n)
  for i in 1:n
    result = take!(c)

    if isa(result, Thrown)  # all exceptions are Thrown{TaskFailedException}
      exceptions[i] = if isa(result.exception.task.exception, MultipleExceptions)
        # extract MultipleExceptions for better readability
        result.exception.task.exception
      else
        result
      end
      # we found an error, hence need to wait for another result
      continue
    end

    # no error, i.e. we found the one result we wanted to find
    close(c)  # we only need one result
    for t in (tasks′..., tasks...)
      # all tasks can be interrupted now that we have the first result
      @async Base.throwto(t, InterruptException())
    end
    return result  # break for loop
  end
  # only if all n tasks failed
  return throw(MultipleExceptions(exceptions...))
end

# FunctorApplicativeMonad
# =======================

function TypeClasses.foreach(f, x::Task)
  f(fetch(x))
  nothing
end

TypeClasses.map(f, x::Task) = @async f(fetch(x))
TypeClasses.pure(::Type{<:Task}, x) = @async x

# we use the default implementation of ap which follows from flatten
# TypeClasses.ap
TypeClasses.ap(f::Task, x::Task) = @async fetch(f)(fetch(x))
# we don't use convert for typesafety, as fetch is more flexible and also enables typechecks
# e.g. this works seamlessly to combine a Future into a Task
TypeClasses.flatmap(f, x::Task) = @async fetch(f(fetch(x)))

# FlipTypes
# =========

# does not make much sense as this would need to execute the Task, and map over its returned value,
# creating a bunch dummy Tasks within.


text:
function runfig4()
    p = Params(Ln = 2000, W = 2700, Ls =8.649,  Ws = 8.649, scale = 40, λ= 5, α = 0,
        EZ = SA[0, 0.6, 0], μN = 1.192, Δ = 0.3, d = 0, τ = 1);

    presets_fig = Fig4_presets(0.0, 0.0, 0, 4, false)
    data = ldosonlattice_averaged_sc(p, presets_fig)
    savepsi("fig4", p, data[1])

    presets_fig = Fig4_presets(0.0, 0.0, 0, 4, true)
    data = ldosonlattice_averaged_sc(p, presets_fig)
    savepsi("fig4", p, data[1])

    sp = spectrumsweepb(p, 0:.1:2, true)
    savespectrum("fig4", p, sp[1], sp[2])

    sp = spectrumsweepb(p, 0:.1:2, false)
    savespectrum("fig4", p, sp[1], sp[2])
end

max_stars_repo_name: fernandopenaranda/MBSinBLG

max_issues_repo_name: fernandopenaranda/MBSinBLG

max_forks_repo_name: fernandopenaranda/MBSinBLG


text:
# http://web.mit.edu/puzzle/www/2012/puzzles/william_s_bergman/behave/

@testset "behave clusters" begin
    behave = shuffle([
            "hugoweaving",
            "mountaindew",
            "mozambique",
            "sequoia",
            "annotation",
            "artificial",
            "individual",
            "omnivorous",
            "onlocation",
            "almost",
            "biopsy",
            "chimp",
            "films",
            "ghost",
            "tux",
            "balked",
            "highnoon",
            "posted"])
    cluster = Collective.best_cluster(corpus, behave, 6)
    @test sort(cluster.words) == ["almost", "biopsy", "chimp", "films", "ghost", "tux"]
    @test cluster.feature.description == "has at least 1 reverse alphabetical bigrams"
    @test !any(cluster.feature.satisfied)
end


text:
@testset "316.remove-duplicate-letters.jl" begin
    @test remove_duplicate_letters("bcabc") == "abc"
    @test remove_duplicate_letters("ecbacba") == "eacb"
    @test remove_duplicate_letters("bddbccd") == "bcd"
end


text:
export RandomPolicy

using Random

"""
    RandomPolicy(action_space=nothing; rng=Random.GLOBAL_RNG)

If `action_space` is `nothing`, then it will use the `legal_action_space` at
runtime to randomly select an action. Otherwise, a random element within
`action_space` is selected.

!!! note
    You should always set `action_space=nothing` when dealing with environments
    of `FULL_ACTION_SET`.
"""
struct RandomPolicy{S,R<:AbstractRNG} <: AbstractPolicy
    action_space::S
    rng::R
end

Random.seed!(p::RandomPolicy, seed) = Random.seed!(p.rng, seed)

RandomPolicy(s = nothing; rng = Random.GLOBAL_RNG) = RandomPolicy(s, rng)

(p::RandomPolicy{Nothing})(env) = rand(p.rng, legal_action_space(env))
(p::RandomPolicy)(env) = rand(p.rng, p.action_space)

function RLBase.prob(p::RandomPolicy{<:Union{AbstractVector,Tuple}}, env::AbstractEnv)
    prob(p, state(env))
end

function RLBase.prob(p::RandomPolicy{<:Union{AbstractVector,Tuple}}, s)
    n = length(p.action_space)
    Categorical(fill(1 / n, n); check_args = false)
end

RLBase.prob(p::RandomPolicy{Nothing}, env::AbstractEnv) = prob(p, env, ChanceStyle(env))
RLBase.prob(p::RandomPolicy{Nothing}, x) =
    @error "no I really don't know how to calculate the prob from nothing"

function RLBase.prob(
    p::RandomPolicy{Nothing},
    env::AbstractEnv,
    ::RLBase.AbstractChanceStyle,
)
    mask = legal_action_space_mask(env)
    n = sum(mask)
    prob = zeros(length(mask))
    prob[mask] .= 1 / n
    prob
end

function RLBase.prob(
    p::RandomPolicy{Nothing},
    env::AbstractEnv,
    ::RLBase.ExplicitStochastic,
)
    if current_player(env) == chance_player(env)
        prob(env, chance_player(env))
    else
        prob(p, env, DETERMINISTIC)
    end
end

RLBase.update!(p::RandomPolicy, x) = nothing

RLBase.prob(p::RandomPolicy, env_or_state, a) = 1 / length(p.action_space)

function RLBase.prob(p::RandomPolicy{Nothing}, env::AbstractEnv, a)
    # we can safely assume s is discrete here.
    s = legal_action_space(env)
    if a in s
        1.0 / length(s)
    else
        0.0
    end
end


text:
# # Planar One-Dimensional Flow

# This a WAVI.jl's simplest example:
# flow down a flat plane in one horizontal dimension. This example demonstrates
#
#   * How to load `WAVI.jl`.
#   * How to instantiate an `WAVI.jl` model.
#   * How to create simple `WAVI.jl` output.
#   * How to time-step a model forward.
#   * How to look at results.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# using Pkg
# Pkg.add("https://github.com/RJArthern/WAVI.jl"), Pkg.add(Plots)
using WAVI, Plots

# ## Instantiating and configuring a model
# We first build a WAVI `model`, by passing it a grid, information about the problem we would like to solve.
#
# Below, we build a grid with 300 grid points in the `x` direction. We use 2 grid points in the `y` direction (the minimum number of grid points in any dimension). This grid has a resolution of 12km.
grid = Grid(nx = 300, ny = 2, dx = 12000.0, dy = 12000.0);

# Next, we write a function which defines the WAVI.jl accepts both functions and arrays (of the same size as the grid) as bed inputs, but here we'll use a function for simplicity
function bed_elevation(x,y)
    B = 720 - 778.5 * x ./ (750e3);
    return B;
end
# This bed drops a height of 778.5m in every 750km, which is a typical scale for an ice sheet in Antarctica.

# Next we specify physical parameters via a `Params` object. In this case, we set the accumulation rate (the net snowfall).
params = Params(accumulation_rate = 0.3);

# Let's also set the thickness of the ice to be 300m everywhere, to start with. Initial conditions are set via `InitialConditions` objects
initial_conditions = InitialConditions(initial_thickness = 300. .* ones(grid.nx, grid.ny));

# Now we are ready to build a `Model` by assembling these pieces:
model = Model(grid = grid, bed_elevation = bed_elevation, params = params);

# ## Simple model visualization
# Before we can vizualize anything, we need to update the model so that velocities are appropriate for this particular ice thickness. We do this with the `update_state!` function:
update_state!(model)

# Now let's plot the ice profile and ice velocity, starting with the bed
ice_plot = plot(model.grid.xxh[:,1]/1e3, model.fields.gh.b[:,1],
                linewidth = 2,
                linecolor = :brown,
                label = "bed",
                xlabel = "x (km)",
                ylabel = "z (m)");

# Now add the ice surface
plot!(ice_plot, model.grid.xxh[:,1]/1e3, model.fields.gh.s[:,1],
                linewidth = 2,
                linecolor = :blue,
                label = "ice surface");

# And finally the ice base
plot!(ice_plot, model.grid.xxh[:,1]/1e3, model.fields.gh.s[:,1] .- model.fields.gh.h[:,1],
                linewidth = 2,
                linecolor = :red,
                label = "ice base");
# We see that the ice shelf goes afloat when the ice base is approximately 270m below sea level: 270m is the product of the ratio of the densities of ice (about 900 km/m^3) and ocean (about 1000 kg/m^3) with the ice thickness (300m).

# Lets also look at the velocity in the ice
vel_plot = plot(model.grid.xxh[:,1]/1e3, model.fields.gh.u[:,1],
                linewidth = 2,
                label = "ice velocity",
                xlabel = "x (km)",
                ylabel = "ice velocity (m/yr)");

p = Plots.plot(ice_plot,vel_plot, layout  = (2,1))
#display(p)
# Ice velocities are very small (but non-zero) in the grounded ice, where friction between the ice and the bed restrains the flow. In the shelf, where there is no basal friction, velocities increase linearly to a maximum of 250 m/yr at the downstream end of the shelf.

# ## Running a `Simulation`
#
# Now, let's think about advancing time. To do so, we set up a simulation, which time-steps the model forward and manages output.
#
# A `TimesteppingParams` object controls parameters related to timestepping. Let's set the model to run for 1000 years with a timestep of 0.5 years:
timestepping_params = TimesteppingParams(dt = 0.5, end_time = 10000.);

# Now we can build the `Simulation` object and then run it!
simulation = Simulation(model = model, timestepping_params = timestepping_params);
run_simulation!(simulation);
# Our simulation ran: `simulation` holds all the information about the state at time 100 years.

# ## Outputting the solution
# Our simulation ran successfully, but we don't have any information about what happened. We get around this by outputting the solution regularly.
# First, make a clean folder where solution files will go:
folder = joinpath(@__DIR__, "planar_one_dimensional_flow");
isdir(folder) && rm(folder, force = true, recursive = true);
mkdir(folder) ;

# What and when to output is specified in WAVI.jl by an instance of an `OutputtingParams` objects. Let's set one up so that the ice thickness, (unchanging) bed, ice surface and ice velocity is output every 100 years:
output_params = OutputParams(outputs = (h = model.fields.gh.h,u = model.fields.gh.u, b = model.fields.gh.b,s = model.fields.gh.s),
                            output_freq = 1000.,
                            output_path = folder);
# Note that the `outputs` keyword argument takes a named tuple, which points to the locations of fields that are to be outputted.

# Let's build a new simulation, which knows about the outputting...
simulation = Simulation(model = model, timestepping_params = timestepping_params, output_params = output_params);

# ...and run it
run_simulation!(simulation);

# ## Visualizing the results

# Let's look at how the shape of the ice sheet changes during the simulation. To do, we'll loop over the output files and put the thickness and surface info a matrix
files = [joinpath(folder, file) for file in readdir(folder) if endswith( joinpath(folder, file), ".jld2") ] ;
nout = length(files);
h_out = zeros(simulation.model.grid.nx, nout);
surface_out = zeros(simulation.model.grid.nx, nout);
base_out = zeros(simulation.model.grid.nx, nout);
t_out = zeros(1,nout);

for i = 1:nout
    d = load(files[i]);
    h_out[:,i] = d["h"][:,1];
    base_out[:,i] = d["s"][:,1] .- d["h"][:,1];
    surface_out[:,i] = d["s"][:,1];
    t_out[i] = d["t"];
end

# Now lets make the plot
pl = Plots.plot(simulation.model.grid.xxh[:,1], simulation.model.fields.gh.b[:,1],
                    linecolor = :brown,
                    xlabel = "x (km)",
                    ylabel = "z (m)",
                    legend = :none)
Plots.plot!(pl,simulation.model.grid.xxh[:,1], surface_out, legend = :none, linecolor = :blue)
Plots.plot!(pl,simulation.model.grid.xxh[:,1], base_out, legend = :none, linecolor = :red)

# Finally, we clear up the files we just outputted
rm(folder, force = true, recursive = true);


text:
using KUparser, KUnet, HDF5, JLD
KUnet.atype(Array)
@date @load "zn11oparse1.jld"
for y in (:ytrn, :ydev, :ytst); @eval $y = full($y); end
@show net = Layer[KPerceptron(size(ydev,1), KUnet.kpoly, [1,2])]
for i=1:100
    @date train(net, xtrn, ytrn; iters=100)
    @date println((i, size(net[1].s,2),
                   # accuracy(ytrn, predict(net, xtrn)),
                   accuracy(ydev, predict(net, xdev)),
                   accuracy(ytst, predict(net, xtst)),
                   ))
end
:ok


text:
"""
    compute_reduced_clique_graph!(sep::Array{Set{Int64}, 1}, snd::Array{Set{Int64}, 1})

Compute the reduced clique graph (union of all clique trees) given an initial clique tree defined by its supernodes `snd` and separator `sep` sets.

We are using the algorithm described in **Michel Habib and Juraj Stacho - Polynomial-time algorithm for the leafage ofchordal graphs (2009)**, which
computes the reduced clique graph in the following way:
1. Sort all minimal separators by size
2. Initialise graph CG(R) with cliques as nodes and no edges
3. for largest unprocessed separator S and
    |  add an edge between any two cliques C1 and C2 if they both contain S and are in different connected   components   of CG(R) and store in `edges`.
    |  Compute an edge weight used for merge decision and store in `val`.
    |  Store the index of the separator which is the intersection C1 ∩ C2 in `iter`
   end
"""
function compute_reduced_clique_graph!(sep::Array{Set{Int64}, 1}, snd::Array{Set{Int64}, 1})
    # loop over separators by decreasing cardinality
    sort!(sep, by = x -> length(x), rev = true)

    # edges = Array{Tuple{Int64, Int64}, 1}()    # a list of edges in the reduced clique graph, higher clique index first
    rows = Int64[]
    cols = Int64[]
    # inter = Array{Int64, 1}()                  # the index of the separator which corresponds to the intersection of the two cliques

    for (k, separator) in enumerate(sep)
        # find cliques that contain the separator
        clique_ind = findall(x -> separator ⊆ x, snd)

        # we compute the separator graph (see Habib, Stacho - Reduced clique graphs of chordal graphs) to analyse connectivity
        # we represent the separator graph H by a hashtable
        H = separator_graph(clique_ind, separator, snd)
        # find the connected components of H
        components = find_components(H, clique_ind)
        # for each pair of cliques that contain the separator, add an edge to the reduced clique tree if they are in unconnected components
        for pair in subsets(clique_ind, Val{2}())
            if is_unconnected(pair, components)
                # push!(edges, (max(pair...), min(pair...))) #add edge
                push!(rows, max(pair...))
                push!(cols, min(pair...))
                # push!(inter, k) # store intersection
            end
        end

    end
    return rows, cols#, inter
end

"Check whether the `pair` of cliques are in different `components`."
function is_unconnected(pair::Tuple{Int64, Int64}, components::Array{Set{Int64}, 1})
    component_ind = findfirst(x -> pair[1] ∈ x, components)
    return pair[2] ∉ components[component_ind]
end

"Find the separator graph H given a separator and the relevant index-subset of cliques."
function separator_graph(clique_ind::Array{Int64,1}, separator::Set{Int64}, snd::Array{Set{Int64}, 1})

    # make the separator graph using a hash table
    # key: clique_ind --> edges to other clique indices
    H = Dict{Int64, Array{Int64, 1}}()

    for pair in subsets(clique_ind, Val{2}())
        ca = pair[1]
        cb = pair[2]
         # if intersect_dim(snd[ca], snd[cb]) > length(separator)
        if !inter_equal(snd[ca], snd[cb], separator)
            if haskey(H, ca)
                push!(H[ca], cb)
            else
                H[ca] = [cb]
            end
            if haskey(H, cb)
                push!(H[cb], ca)
            else
                H[cb] = [ca]
            end
        end
    end
    # add unconnected cliques
    for v in clique_ind
        !haskey(H, v) && (H[v] = Int64[])
    end
    return H
end

"Find connected components in undirected separator graph represented by `H`."
function find_components(H::Dict{Int64, Array{Int64, 1}}, clique_ind::Array{Int64,1})
    visited = Dict{Int64, Bool}(v => false for v in clique_ind)
    components = Array{Set{Int64}, 1}()
    for v in clique_ind
        if visited[v] == false
            component = Set{Int64}()
            push!(components, DFS_hashtable!(component, v, visited, H))
        end
    end
    return components
end

"Depth first search on a hashtable `H`."
function DFS_hashtable!(component::Set{Int64}, v::Int64, visited::Dict{Int64, Bool}, H::Dict{Int64, Array{Int64, 1}})
    visited[v] = true
    push!(component, v)
    for n in H[v]
        if visited[n] == false
            DFS_hashtable!(component, n, visited, H)
        end
    end
    return component
end

"Check if s ∩ s2 == s3."
function inter_equal(s::Set{Int64}, s2::Set{Int64}, s3::Set{Int64})
    dim = 0
    len_s3 = length(s3)
    if length(s) < length(s2)
        sa = s
        sb = s2
    else
        sa = s2
        sb = s
    end
    for e in sa
        if e in sb
            dim += 1
            dim > len_s3 && return false
            e in s3 || return false
        end
    end
    return true
end

"Given a list of edges, return an adjacency hash-table `table` with nodes from 1 to `num_vertices`."
function compute_adjacency_table(edges::SparseMatrixCSC{Float64, Int64}, num_vertices::Int64)
    table = Dict(i => Set{Int64}() for i = 1:num_vertices)
    r = edges.rowval
    c = edges.colptr
     for col = 1:num_vertices
         for i in c[col]:c[col+1]-1
             row = r[i]
             push!(table[row], col)
            push!(table[col], row)
        end
     end
     return table
end

"Check whether `edge` is permissible for a merge. An edge is permissible if for every common neighbor N, C_1 ∩ N == C_2 ∩ N or if no common neighbors exist."
function ispermissible(edge::Tuple{Int64, Int64}, adjacency_table::Dict{Int64, Set{Int64}}, snd::Array{Set{Int64}, 1})
    c_1 = edge[1]
    c_2 = edge[2]
    common_neighbors = intersect(adjacency_table[c_1], adjacency_table[c_2])
    # N.B. This can be made faster by first checking whether the sizes of the intersection are the same before allocating anything
    for neighbor in common_neighbors
        intersect(snd[c_1], snd[neighbor]) != intersect(snd[c_2], snd[neighbor]) && return false
    end
    return true
end


text:
@testset "Quick" begin
    qm = QuickMDP(statetype=Int, actiontype=Int)
    @test statetype(qm) == Int
    @test actiontype(qm) == Int
    @test discount(qm) == 1.0
    @test isterminal(qm, 1) == false
    @test_throws MissingQuickArgument transition(qm, 1, 1)
    @test_throws MissingQuickArgument initialstate_distribution(qm)
    @test_throws MissingQuickArgument reward(qm, 1, 1, 1)
    @test_throws MissingQuickArgument initialstate(qm, MersenneTwister(2))
    @test_throws MissingQuickArgument states(qm)
    @test_throws MissingQuickArgument actions(qm)
    @test_throws MissingQuickArgument stateindex(qm, 1)
    @test_throws MissingQuickArgument actionindex(qm, 1)

    qp = QuickPOMDP(obstype=Int)
    @test obstype(qp) == Int
    @test statetype(qp) == Any
    @test actiontype(qp) == Any
    @test discount(qp) == 1.0
    @test isterminal(qp, 1) == false
    @test_throws MissingQuickArgument transition(qp, 1, 1)
    @test_throws MissingQuickArgument observation(qp, 1, 1)
    @test_throws MissingQuickArgument initialstate_distribution(qp)
    @test_throws MissingQuickArgument reward(qp, 1, 1, 1)
    @test_throws MissingQuickArgument initialstate(qp, MersenneTwister(20))
    @test_throws MissingQuickArgument initialobs(qp, 1, MersenneTwister(1))
    @test_throws MissingQuickArgument states(qp)
    @test_throws MissingQuickArgument actions(qp)
    @test_throws MissingQuickArgument observations(qp)
    @test_throws MissingQuickArgument stateindex(qp, 1)
    @test_throws MissingQuickArgument actionindex(qp, 1)
    @test_throws MissingQuickArgument obsindex(qp, 1)

    @test_throws MissingQuickArgument gen(DDNOut(:sp,:o,:r), qp, 1, 1, MersenneTwister(2))
end

@testset "Mountan Car" begin
    mountaincar = begin
        include("../examples/mountaincar.jl")
    end
    energize = FunctionPolicy(
        function (s)
            if s[2] < 0.0
                return minimum(actions(mountaincar))
            else
                return maximum(actions(mountaincar))
            end
        end
    )
    sim = HistoryRecorder(max_steps=1000)
    hist = simulate(sim, mountaincar, energize)
    @test maximum(hist[:r]) == 100.0
    @test all(isa.(hist[:s], statetype(mountaincar)))
    @test last(collect(hist[:sp]))[1] > 0.5
end

@testset "Mountan Car Visualization" begin
    m = begin
        include("../examples/mountaincar_with_visualization.jl")
    end
    energize = FunctionPolicy(
        function (s)
            if s[2] < 0.0
                return minimum(actions(m))
            else
                return maximum(actions(m))
            end
        end
    )
    sim = HistoryRecorder(max_steps=1000)
    hist = simulate(sim, m, energize)
    @test maximum(hist[:r]) == 100.0
    @test all(isa.(hist[:s], statetype(mountaincar)))
    @test last(collect(hist[:sp]))[1] > 0.5
    draw(SVG("test_render.svg", 5cm, 4cm), render(m, (s=(0.2, 0.0),)))
    @test isfile("test_render.svg")
end

@testset "Simple Light-Dark" begin
    m = begin
        include("../examples/lightdark.jl")
    end
    p = RandomPolicy(m, rng = MersenneTwister(2))
    sim = HistoryRecorder(max_steps=1000)
    hist = simulate(sim, m, p)
    @test all(-100 .<= hist[:r] .<= 100)
    @test all(isa.(hist[:s], statetype(m)))
    @test last(collect(hist[:sp])) == 61
end


text:
using Pkg
Pkg.add("JSON3")

import JSON3
using Sockets

struct Coordinate
    x::Float64
    y::Float64
    z::Float64
end

function calc(text)
    jobj = JSON3.read(text)
    coordinates = jobj["coordinates"]
    len = length(coordinates)
    x = y = z = 0

    for coord in coordinates
        x += coord["x"]
        y += coord["y"]
        z += coord["z"]
    end

    Coordinate(x / len, y / len, z / len)
end

function notify(msg)
    try
        socket = connect("localhost", 9001)
        write(socket, msg)
        close(socket)
    catch
        # standalone usage
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    right = Coordinate(2.0, 0.5, 0.25)
    for v in [
        """{"coordinates":[{"x":2.0,"y":0.5,"z":0.25}]}""",
        """{"coordinates":[{"y":0.5,"x":2.0,"z":0.25}]}""",
    ]
        left = calc(v)
        if left != right
            println(stderr, "$(left) != $(right)")
            exit(1)
        end
    end

    text = open("/tmp/1.json") do file
        read(file, String)
    end

    notify("Julia (JSON3)\t$(getpid())")
    results = calc(text)
    notify("stop")

    println(results)
end

max_stars_repo_name: clemenswasser/benchmarks

max_issues_repo_name: clemenswasser/benchmarks

max_forks_repo_name: clemenswasser/benchmarks


text:
function load(filepath::String; verbose=true, mmap = false)
    open(filepath) do io
        load(io; verbose=verbose, mmap=mmap)
    end
end

load(io::IOStream; verbose=true, mmap = false) = load(read(io, TiffFile); verbose=verbose, mmap=mmap)
function load(tf::TiffFile; verbose=true, mmap = false)
    ifds = IFD{offset(tf)}[]

    nplanes = 0
    for ifd in tf
        load!(tf, ifd)
        push!(ifds, ifd)
        nplanes += 1
    end

    if mmap
        loaded = DiskTaggedImage(tf, ifds)
    else
        if nplanes == 1
            loaded = load(tf, ifds, nothing; verbose=verbose)
        else
            loaded = load(tf, ifds, nplanes; verbose=verbose)
        end
    end
    data = fixcolors(loaded, first(ifds))

    close(tf.io)
    return DenseTaggedImage(data, ifds)
end

"""
    fixcolors(loaded, ifd)

Wrap the raw eltype of an image with color if needed, e.g. for space efficient
on-disk representations like palette-colored images and bitarrays. Otherwise,
just return the passed image.
"""
function fixcolors(loaded, ifd)
    if eltype(loaded) <: Palette
        raw = rawtype(ifd)
        loadedr = reinterpret(raw, loaded)
        maxdepth = 2^(Int(ifd[BITSPERSAMPLE].data))-1
        colors = ifd[COLORMAP].data
        color_map = vec(reinterpret(RGB{N0f16}, reshape(colors, :, 3)'))
        return IndirectArray(loadedr, OffsetArray(color_map, 0:maxdepth))
    elseif eltype(loaded) <: Bool
        return Gray.(loaded)
    else
        return loaded
    end
end

function load(tf::TiffFile, ifds::AbstractVector{<:IFD}, ::Nothing; verbose = true)
    ifd = ifds[1]
    cache = getcache(ifd)
    read!(cache, tf, ifd)

    return Matrix(cache')
end

function load(tf::TiffFile, ifds::AbstractVector{<:IFD}, N; verbose = true)
    ifd = ifds[1]

    cache = getcache(ifd)

    data = similar(cache, nrows(ifd), ncols(ifd), N)

    freq = verbose ? 1 : Inf
    @showprogress freq for (idx, ifd) in enumerate(ifds)
        read!(cache, tf, ifd)
        data[:, :, idx] .= cache'
    end

    return data
end


text:
include("data.jl")
using Flux, Statistics
using Flux: onehot, onehotbatch, throttle, logitcrossentropy, reset!, onecold
using Parameters: @with_kw

@with_kw mutable struct Args
    lr::Float64 = 1e-3    # Learning rate
    epochs::Int = 20      # Number of epochs for training
    train_len::Int = 100000  # Length of training data to be generated
    val_len::Int = 1000     # Length of Validation Data
    throttle::Int = 10    # Throttle timeout
end

function getdata(args)
    # training data of bit strings from length 2 to 50
    train = gendata(args.train_len, 1:50)
    # validation data of bit strings of length 50
    val = gendata(args.val_len, 50)
    return train, val
end

function build_model()
    scanner = LSTM(length(alphabet), 20)
    encoder = Dense(20, length(alphabet))
    return scanner, encoder
end

function model(x, scanner, encoder)
    state = scanner.(x.data)[end]
    reset!(scanner)
    encoder(state)
end

function train(; kws...)
    # Initialize the parameters
    args = Args(; kws...)

    # Load Data
    train_data, val_data = getdata(args)

    @info("Constructing Model...")
    scanner,encoder = build_model()

    loss(x, y) = logitcrossentropy(model(x, scanner, encoder), y)
    batch_loss(data) = mean(loss(d...) for d in data)

    opt = ADAM(args.lr)
    ps = params(scanner, encoder)
    evalcb = () -> @show batch_loss(val_data)

    @info("Training...")
    for i=1:args.epochs
        Flux.train!(loss, ps, train_data, opt, cb=throttle(evalcb, args.throttle))
    end

    return scanner, encoder
end

function test(scanner, encoder)
    # sanity test
    tx = map(c -> onehotbatch(c, alphabet), [
        [false, true], # 01 -> 1
        [true, false], # 10 -> 1
        [false, false], # 00 -> 0
        [true, true]]) # 11 -> 0
    @info("Test...")
    out = [onecold(model(x, scanner, encoder)) - 1 for x in tx]
    input = [[0,1],[1,0],[0,0],[1,1]]
    for i in 1:length(tx)
        print(input[i]," => ",out[i],"\n")
    end
end

cd(@__DIR__)
scanner, encoder = train()
test(scanner, encoder)

max_stars_repo_name: shawn-kern/model-zoo

max_issues_repo_name: shawn-kern/model-zoo

max_forks_repo_name: shawn-kern/model-zoo


text:
using .Zygote: Zygote, @adjoint, @nograd, pullback

using Compat: eachcol

@adjoint istraining() = true, _ -> nothing
@nograd Bijectors._debug

@adjoint function mapvcat(f, args...)
    g(f, args...) = map(f, args...)
    return pullback(g, f, args...)
end
@adjoint function eachcolmaphcat(f, x1, x2)
    function g(f, x1, x2)
        init = reshape(f(view(x1, :, 1), x2[1]), :, 1)
        return reduce(hcat, [f(view(x1, :, i), x2[i]) for i in 2:size(x1, 2)]; init = init)
    end
    return pullback(g, f, x1, x2)
end
@adjoint function eachcolmaphcat(f, x)
    function g(f, x)
        init = reshape(f(view(x, :, 1)), :, 1)
        return reduce(hcat, [f(view(x, :, i)) for i in 2:size(x, 2)]; init = init)
    end
    return pullback(g, f, x)
end
@adjoint function sumeachcol(f, x1, x2)
    g(f, x1, x2) = sum([f(view(x1, :, i), x2[i]) for i in 1:size(x1, 2)])
    return pullback(g, f, x1, x2)
end

@adjoint function logabsdetjac(b::Log{1}, x::AbstractVector)
    return -sum(log, x), Δ -> (nothing, -Δ ./ x)
end
@adjoint function logabsdetjac(b::Log{1}, x::AbstractMatrix)
    return -vec(sum(log, x; dims = 1)), Δ -> (nothing, .- Δ' ./ x)
end

# AD implementations
function jacobian(
    b::Union{<:ADBijector{<:ZygoteAD}, Inverse{<:ADBijector{<:ZygoteAD}}},
    x::Real
)
    return Zygote.gradient(b, x)[1]
end
function jacobian(
    b::Union{<:ADBijector{<:ZygoteAD}, Inverse{<:ADBijector{<:ZygoteAD}}},
    x::AbstractVector{<:Real}
)
    return Zygote.jacobian(b, x)
end
@adjoint function _logabsdetjac_scale(a::Real, x::Real, ::Val{0})
    return _logabsdetjac_scale(a, x, Val(0)), Δ -> (inv(a) .* Δ, nothing, nothing)
end
@adjoint function _logabsdetjac_scale(a::Real, x::AbstractVector, ::Val{0})
    J = fill(inv.(a), length(x))
    return _logabsdetjac_scale(a, x, Val(0)), Δ -> (transpose(J) * Δ, nothing, nothing)
end
@adjoint function _logabsdetjac_scale(a::Real, x::AbstractMatrix, ::Val{0})
    J = fill(size(x, 1) / a, size(x, 2))
    return _logabsdetjac_scale(a, x, Val(0)), Δ -> (transpose(J) * Δ, nothing, nothing)
end
@adjoint function _logabsdetjac_scale(a::AbstractVector, x::AbstractVector, ::Val{1})
    # ∂ᵢ (∑ⱼ log|aⱼ|) = ∑ⱼ δᵢⱼ ∂ᵢ log|aⱼ|
    #                 = ∂ᵢ log |aᵢ|
    #                 = (1 / aᵢ) ∂ᵢ aᵢ
    #                 = (1 / aᵢ)
    J = inv.(a)
    return _logabsdetjac_scale(a, x, Val(1)), Δ -> (J .* Δ, nothing, nothing)
end
@adjoint function _logabsdetjac_scale(a::AbstractVector, x::AbstractMatrix, ::Val{1})
    Jᵀ = repeat(inv.(a), 1, size(x, 2))
    return _logabsdetjac_scale(a, x, Val(1)), Δ -> (Jᵀ * Δ, nothing, nothing)
end
## Positive definite matrices
@adjoint function replace_diag(::typeof(log), X)
    f(i, j) = i == j ? log(X[i, j]) : X[i, j]
    out = f.(1:size(X, 1), (1:size(X, 2))')
    out, ∇ -> begin
        g(i, j) = i == j ? ∇[i, j] / X[i, j] : ∇[i, j]
        (nothing, g.(1:size(X, 1), (1:size(X, 2))'))
    end
end
@adjoint function replace_diag(::typeof(exp), X)
    f(i, j) = ifelse(i == j, exp(X[i, j]), X[i, j])
    out = f.(1:size(X, 1), (1:size(X, 2))')
    out, ∇ -> begin
        g(i, j) = ifelse(i == j, ∇[i, j] * exp(X[i, j]), ∇[i, j])
        (nothing, g.(1:size(X, 1), (1:size(X, 2))'))
    end
end

@adjoint function pd_logpdf_with_trans(
    d,
    X::AbstractMatrix{<:Real},
    transform::Bool,
)
    return pullback(pd_logpdf_with_trans_zygote, d, X, transform)
end
function pd_logpdf_with_trans_zygote(
    d,
    X::AbstractMatrix{<:Real},
    transform::Bool,
)
    T = eltype(X)
    Xcf = cholesky(X, check = false)
    if !issuccess(Xcf)
        Xcf = cholesky(X + max(eps(T), eps(T) * norm(X)) * I, check = true)
    end
    lp = getlogp(d, Xcf, X)
    if transform && isfinite(lp)
        U = Xcf.U
        @inbounds for i in 1:dim(d)
            lp += (dim(d) - i + 2) * log(U[i, i])
        end
        lp += dim(d) * log(T(2))
    end
    return lp
end

# Simplex adjoints

@adjoint function _simplex_bijector(X::AbstractVector, b::SimplexBijector{1})
    return _simplex_bijector(X, b), Δ -> (simplex_link_jacobian(X)' * Δ, nothing)
end
@adjoint function _simplex_inv_bijector(Y::AbstractVector, b::SimplexBijector{1})
    return _simplex_inv_bijector(Y, b), Δ -> (simplex_invlink_jacobian(Y)' * Δ, nothing)
end

@adjoint function _simplex_bijector(X::AbstractMatrix, b::SimplexBijector{1})
    return _simplex_bijector(X, b), Δ -> begin
        maphcat(eachcol(X), eachcol(Δ)) do c1, c2
            simplex_link_jacobian(c1)' * c2
        end, nothing
    end
end
@adjoint function _simplex_inv_bijector(Y::AbstractMatrix, b::SimplexBijector{1})
    return _simplex_inv_bijector(Y, b), Δ -> begin
        maphcat(eachcol(Y), eachcol(Δ)) do c1, c2
            simplex_invlink_jacobian(c1)' * c2
        end, nothing
    end
end

@adjoint function logabsdetjac(b::SimplexBijector{1}, x::AbstractVector)
    return logabsdetjac(b, x), Δ -> begin
        (nothing, simplex_logabsdetjac_gradient(x) * Δ)
    end
end
@adjoint function logabsdetjac(b::SimplexBijector{1}, x::AbstractMatrix)
    return logabsdetjac(b, x), Δ -> begin
        (nothing, maphcat(eachcol(x), Δ) do c, g
            simplex_logabsdetjac_gradient(c) * g
        end)
    end
end

# LocationScale fix

@adjoint function minimum(d::LocationScale)
    function _minimum(d)
        m = minimum(d.ρ)
        if isfinite(m)
            return d.μ + d.σ * m
        else
            return m
        end
    end
    return pullback(_minimum, d)
end
@adjoint function maximum(d::LocationScale)
    function _maximum(d)
        m = maximum(d.ρ)
        if isfinite(m)
            return d.μ + d.σ * m
        else
            return m
        end
    end
    return pullback(_maximum, d)
end
@adjoint function lower(A::AbstractMatrix)
    return lower(A), Δ -> (lower(Δ),)
end
@adjoint function getpd(X::AbstractMatrix)
    return LowerTriangular(X) * LowerTriangular(X)', Δ -> begin
        Xl = LowerTriangular(X)
        return (LowerTriangular(Δ' * Xl + Δ * Xl),)
    end
end
@adjoint function pd_link(X::AbstractMatrix{<:Real})
    return pullback(X) do X
        Y = cholesky(X; check = true).L
        return replace_diag(log, Y)
    end
end

@adjoint function _inv_link_chol_lkj(y)
    K = LinearAlgebra.checksquare(y)

    w = similar(y)

    z_mat = similar(y) # cache for adjoint
    tmp_mat = similar(y)

    @inbounds for j in 1:K
        w[1, j] = 1
        for i in 2:j
            z = tanh(y[i-1, j])
            tmp = w[i-1, j]

            z_mat[i, j] = z
            tmp_mat[i, j] = tmp

            w[i-1, j] = z * tmp
            w[i, j] = tmp * sqrt(1 - z^2)
        end
        for i in (j+1):K
            w[i, j] = 0
        end
    end

    function pullback_inv_link_chol_lkj(Δw)
        LinearAlgebra.checksquare(Δw)

        Δy = zero(y)

        @inbounds for j in 1:K
            Δtmp = Δw[j,j]
            for i in j:-1:2
                Δz = Δw[i-1, j] * tmp_mat[i, j] - Δtmp * tmp_mat[i, j] / sqrt(1 - z_mat[i, j]^2) * z_mat[i, j]
                Δy[i-1, j] = Δz / cosh(y[i-1, j])^2
                Δtmp = Δw[i-1, j] * z_mat[i, j] + Δtmp * sqrt(1 - z_mat[i, j]^2)
            end
        end

        return (Δy,)
    end

    return w, pullback_inv_link_chol_lkj
end

@adjoint function _link_chol_lkj(w)
    K = LinearAlgebra.checksquare(w)

    z = similar(w)

    @inbounds z[1, 1] = 0

    tmp_mat = similar(w) # cache for pullback.

    @inbounds for j=2:K
        z[1, j] = atanh(w[1, j])
        tmp = sqrt(1 - w[1, j]^2)
        tmp_mat[1, j] = tmp
        for i in 2:(j - 1)
            p = w[i, j] / tmp
            tmp *= sqrt(1 - p^2)
            tmp_mat[i, j] = tmp
            z[i, j] = atanh(p)
        end
        z[j, j] = 0
    end

    function pullback_link_chol_lkj(Δz)
        LinearAlgebra.checksquare(Δz)

        Δw = similar(w)

        @inbounds Δw[1,1] = zero(eltype(Δz))

        @inbounds for j=2:K
            Δw[j, j] = 0
            Δtmp = zero(eltype(Δz)) # Δtmp_mat[j-1,j]
            for i in (j-1):-1:2
                p = w[i, j] / tmp_mat[i-1, j]
                ftmp = sqrt(1 - p^2)
                d_ftmp_p = -p / ftmp
                d_p_tmp = -w[i,j] / tmp_mat[i-1, j]^2

                Δp = Δz[i,j] / (1-p^2) + Δtmp * tmp_mat[i-1, j] * d_ftmp_p
                Δw[i, j] = Δp / tmp_mat[i-1, j]
                Δtmp = Δp * d_p_tmp + Δtmp * ftmp # update to "previous" Δtmp
            end
            Δw[1, j] = Δz[1, j] / (1-w[1,j]^2) - Δtmp / sqrt(1 - w[1,j]^2) * w[1,j]
        end

        return (Δw,)
    end

    return z, pullback_link_chol_lkj

end


text:
# Returns two booleans where the first indicates e-box dominance and the 2nd indicates
# epsilon progress
function epsilon_box_dominates_and_epsilon_progress(u, v, epsilon)
    uindex = floor(u ./ epsilon)
    vindex = floor(v ./ epsilon)
    inner_epsilon_box_dominates_and_epsilon_progress(
        u, uindex, (epsilon .* uindex),
        v, vindex, (epsilon .* vindex))
end

# FIXME
epsilon_box_dominates_and_epsilon_progress(f1::NewFitness, f2::NewFitness, epsilon) = epsilon_box_dominates_and_epsilon_progress(fitnessvalues(f1), fitnessvalues(f2), epsilon)
epsilon_box_dominates(u, v, epsilon) = epsilon_box_dominates_and_epsilon_progress(u, v, epsilon)[1]
epsilon_box_progress(u, v, epsilon) = epsilon_box_dominates_and_epsilon_progress(u, v, epsilon)[2]

# Inner boolean expression when the intermediate results have already been calculated.
# This is to speed up processing when there is a whole archive of solutions to compare
# a new solution to.
function inner_epsilon_box_dominates_and_epsilon_progress(
    u, uindex, uindextimesepsilon,
    v, vindex, vindextimesepsilon)

    if pareto_dominates_fast(uindex, vindex)
        return (true, true)
    elseif all(uindex .== vindex)
        smaller_distance = norm(u .- uindextimesepsilon) < norm(v .- vindextimesepsilon)
        return (smaller_distance, false)
    else
        return (false, false)
    end
end


text:
using Dates, Plots

function create_directory(dirpath)
    directory = joinpath(dirpath, string(now()))
    if !ispath(directory)
        mkpath(directory)
    end
    return directory
end

hair(plt, x, y) = plot!(plt, [x, x], [0, y], linewidth=3, color=:grey, alpha=0.5, label=false)

"""Plot probability mass function.

# Examples
```julia
p1 = plot_distribution(u, p)
p2 = plot_distribution(u, cumsum(p))
directory = create_directory(joinpath("results"))
savefig(p1, joinpath(directory, "pmf.svg"))
savefig(p2, joinpath(directory, "cfd.svg"))
```
"""
function plot_distribution(u, p; plt=plot(), label="")
    x = 1:length(u)
    y = p
    plot!(plt, x, y,
        linewidth=0,
        markershape=:circle,
        ylims=(0, 1.15 * maximum(p)),
        xticks = x,
        yticks = LinRange(0, maximum(p), 11),
        xformatter = tick -> "$(u[Int(tick)])",
        label=label,
        legend=:topleft)
    for (x2, y2) in zip(x, y)
        hair(plt, x2, y2)
    end
    return plt
end


text:
using Plots
using GraphRecipes
import GraphRecipes: get_source_destiny_weight

export fixed_layout

function get_source_destiny_weight(g::Graph)
    L = ne(g)
    sources = Array{Int}(undef, L)
    destiny = Array{Int}(undef, L)
    weights = Array{Float64}(undef, L)
    i = 0
    for v in g.vertices
        for e in v
            i += 1
            sources[i] = e.v1.id
            destiny[i] = e.v2.id
            weights[i] = e.weight
        end
    end
    return sources, destiny, weights
end
function fixed_layout(locations::AbstractArray{T,2}, adjmat::AbstractMatrix; kw...) where T
    return locations[1, :], locations[2, :], zeros(Int, size(locations, 2))
end


text:
using TransformUtils

function uteOdom(x::Array{Float64,1}, vc::Float64, alpha::Float64, dt::Float64)
  L = 2.83; H=076; b=0.5;a=3.78;
  d = [(vc*cos(x[3])-vc/L*tan(x[3])*(a*sin(x[3])+b*cos(x[3])));
       (vc*sin(x[3])+vc/L*tan(x[3])*(a*cos(x[3])-b*sin(x[3])));
        vc/L*tan(alpha)]
  x += dt*d
  x[3] = wrapRad(x[3])
  return x
end

function allOdo(sensors::Array{Float64,2})
  odo = zeros(size(sensors,1),3)
  T0 = 0.0
  for i in 2:size(sensors,1)
    odo[i,:]=uteOdom(vec(odo[i-1,:]),sensors[i,2],sensors[i,3],sensors[i,1]-T0)
    T0 = sensors[i,1]
  end

  return odo
end

vc(v::Float64, alpha::Float64; L=2.80381, H=0.828329 ) = v./(1.0-tan(alpha)*H/L)
ve(v::Float64, alpha::Float64; L=2.80381, H=0.828329 ) = v.*(1.0-tan(alpha).*H/L)
dPhi(v::Float64, alpha::Float64; L=2.80381) = 1.0*v*tan(alpha)/L

function compensateRawsensors(drs::Array{Float64,1};
                    L=2.80381, H=0.828329, whlsf=0.94, strsf=1.0199, strbi=0.00159)
    return whlsf*drs[2], strsf*drs[3]+strbi
end

function uteOdomEasy(x::Array{Float64,1}, whlspd::Float64, strangl::Float64, dt::Float64; L=2.80381, H=0.828329)
  wTb = SE2(x)
  v = vc(whlspd,strangl,L=L,H=H)
  dph = dPhi(v,strangl,L=L)
  #dph = v*tan(strangl)
  bTnb = SE2(dt*[v;0.0;dph])
  wTnb = wTb*bTnb
  return vec(se2vee(wTnb))
end

# perform integration over all data, sensors is rows of data, by columns of sensors
function allOdoEasy(sensors::Array{Float64,2};
                    L=2.80381, H=0.828329, whlsf=0.94, strsf=1.0199, strbi=0.00159)
  #msensors = [sensors[:,1]';x[1]*sensors[:,2]';x[5]*sensors[:,3]'+x[2]]';
  odo = zeros(size(sensors,1),3)
  for i in 2:size(sensors,1)
    whlspd, strang = compensateRawsensors(vec(sensors[i,:]), L=L,H=H, whlsf=whlsf, strbi=strbi,strsf=strsf)
    odo[i,:] = uteOdomEasy(vec(odo[i-1,:]), whlspd, strang, (sensors[i,1]-sensors[i-1,1]), L=L, H=H )
  end

  return odo
end