text:
"""
    struct GradedSpace{I<:Sector, D} <: EuclideanSpace{ℂ}
        dims::D
        dual::Bool
    end

A complex Euclidean space with a direct sum structure corresponding to labels in a set `I`,
the objects of which have the structure of a monoid with respect to a monoidal product `⊗`.
In practice, we restrict the label set to be a set of superselection sectors of type
`I<:Sector`, e.g. the set of distinct irreps of a finite or compact group, or the
isomorphism classes of simple objects of a unitary and pivotal (pre-)fusion category.

Here `dims` represents the degeneracy or multiplicity of every sector.

The data structure `D` of `dims` will depend on the result `Base.IteratorElsize(values(I))`;
if the result is of type `HasLength` or `HasShape`, `dims` will be stored in a
`NTuple{N,Int}` with `N = length(values(I))`. This requires that a sector `s::I` can be
transformed into an index via `s == getindex(values(I), i)` and
`i == findindex(values(I), s)`. If `Base.IteratorElsize(values(I))` results `IsInfinite()`
or `SizeUnknown()`, a `SectorDict{I,Int}` is used to store the non-zero degeneracy
dimensions with the corresponding sector as key. The parameter `D` is hidden from the user
and should typically be of no concern.

The concrete type `GradedSpace{I,D}` with correct `D` can be obtained as `Vect[I]`, or if
`I == Irrep[G]` for some `G<:Group`, as `Rep[G]`.
"""
struct GradedSpace{I<:Sector, D} <: EuclideanSpace{ℂ}
    dims::D
    dual::Bool
end
sectortype(::Type{<:GradedSpace{I}}) where {I<:Sector} = I

function GradedSpace{I, NTuple{N, Int}}(dims; dual::Bool = false) where {I, N}
    d = ntuple(n->0, N)
    isset = ntuple(n->false, N)
    for (c, dc) in dims
        i = findindex(values(I), convert(I, c))
        isset[i] && throw(ArgumentError("Sector $c appears multiple times"))
        isset = TupleTools.setindex(isset, true, i)
        d = TupleTools.setindex(d, dc, i)
    end
    return GradedSpace{I, NTuple{N, Int}}(d, dual)
end
GradedSpace{I, NTuple{N, Int}}(dims::Pair; dual::Bool = false) where {I, N} =
    GradedSpace{I, NTuple{N, Int}}((dims,); dual = dual)

function GradedSpace{I, SectorDict{I, Int}}(dims; dual::Bool = false) where {I<:Sector}
    d = SectorDict{I, Int}()
    for (c, dc) in dims
        k = convert(I, c)
        haskey(d, k) && throw(ArgumentError("Sector $k appears multiple times"))
        !iszero(dc) && push!(d, k=>dc)
    end
    return GradedSpace{I, SectorDict{I, Int}}(d, dual)
end
GradedSpace{I, SectorDict{I, Int}}(dims::Pair; dual::Bool = false) where {I<:Sector} =
    GradedSpace{I, SectorDict{I, Int}}((dims,); dual = dual)

GradedSpace{I,D}(; kwargs...) where {I<:Sector,D} = GradedSpace{I,D}((); kwargs...)
GradedSpace{I,D}(d1::Pair, d2::Pair, dims::Vararg{Pair}; kwargs...) where {I<:Sector,D} =
    GradedSpace{I,D}((d1, d2, dims...); kwargs...)

GradedSpace{I}(args...; kwargs...) where {I<:Sector} = Vect[I](args..., kwargs...)

GradedSpace(dims::Tuple{Vararg{Pair{I, <:Integer}}}; dual::Bool = false) where {I<:Sector} =
    Vect[I](dims; dual = dual)
GradedSpace(dims::Vararg{Pair{I, <:Integer}}; dual::Bool = false) where {I<:Sector} =
    Vect[I](dims; dual = dual)
GradedSpace(dims::AbstractDict{I, <:Integer}; dual::Bool = false) where {I<:Sector} =
    Vect[I](dims; dual = dual)
# not inferrable
GradedSpace(g::Base.Generator; dual::Bool = false) =
    GradedSpace(g...; dual = dual)
GradedSpace(g::AbstractDict; dual::Bool = false) =
    GradedSpace(g...; dual = dual)

Base.hash(V::GradedSpace, h::UInt) = hash(V.dual, hash(V.dims, h))

# Corresponding methods:
# properties
dim(V::GradedSpace) =
    reduce(+, dim(V, c) * dim(c) for c in sectors(V); init = zero(dim(one(sectortype(V)))))

dim(V::GradedSpace{I,<:AbstractDict}, c::I) where {I<:Sector} =
    get(V.dims, isdual(V) ? dual(c) : c, 0)
dim(V::GradedSpace{I,<:Tuple}, c::I) where {I<:Sector} =
    V.dims[findindex(values(I), isdual(V) ? dual(c) : c)]

sectors(V::GradedSpace{I,<:AbstractDict}) where {I<:Sector} =
    SectorSet{I}(s->isdual(V) ? dual(s) : s, keys(V.dims))
sectors(V::GradedSpace{I,NTuple{N,Int}}) where {I<:Sector, N} =
    SectorSet{I}(Iterators.filter(n->V.dims[n]!=0, 1:N)) do n
        isdual(V) ? dual(values(I)[n]) : values(I)[n]
    end

hassector(V::GradedSpace{I}, s::I) where {I<:Sector} = dim(V, s) != 0

Base.conj(V::GradedSpace) = typeof(V)(V.dims, !V.dual)
isdual(V::GradedSpace) = V.dual

# equality / comparison
Base.:(==)(V1::GradedSpace, V2::GradedSpace) =
    sectortype(V1) == sectortype(V2) && (V1.dims == V2.dims) && V1.dual == V2.dual

# axes
Base.axes(V::GradedSpace) = Base.OneTo(dim(V))
function Base.axes(V::GradedSpace{I}, c::I) where {I<:Sector}
    offset = 0
    for c′ in sectors(V)
        c′ == c && break
        offset += dim(c′)*dim(V, c′)
    end
    return (offset+1):(offset+dim(c)*dim(V, c))
end

Base.oneunit(S::Type{<:GradedSpace{I}}) where {I<:Sector} = S(one(I)=>1)

# TODO: the following methods can probably be implemented more efficiently for
# `FiniteGradedSpace`, but we don't expect them to be used often in hot loops, so
# these generic definitions (which are still quite efficient) are good for now.
function ⊕(V1::GradedSpace{I}, V2::GradedSpace{I}) where {I<:Sector}
    dual1 = isdual(V1)
    dual1 == isdual(V2) ||
        throw(SpaceMismatch("Direct sum of a vector space and a dual space does not exist"))
    dims = SectorDict{I, Int}()
    for c in union(sectors(V1), sectors(V2))
        cout = ifelse(dual1, dual(c), c)
        dims[cout] = dim(V1, c) + dim(V2, c)
    end
    return typeof(V1)(dims; dual = dual1)
end

function flip(V::GradedSpace{I}) where {I<:Sector}
    if isdual(V)
        typeof(V)(c=>dim(V, c) for c in sectors(V))
    else
        typeof(V)(dual(c)=>dim(V, c) for c in sectors(V))'
    end
end

function fuse(V1::GradedSpace{I}, V2::GradedSpace{I}) where {I<:Sector}
    dims = SectorDict{I, Int}()
    for a in sectors(V1), b in sectors(V2)
        for c in a ⊗ b
            dims[c] = get(dims, c, 0) + Nsymbol(a, b, c)*dim(V1, a)*dim(V2, b)
        end
    end
    return typeof(V1)(dims)
end

function infimum(V1::GradedSpace{I}, V2::GradedSpace{I}) where {I<:Sector}
    if V1.dual == V2.dual
        typeof(V1)(c=>min(dim(V1, c), dim(V2, c)) for c in
            union(sectors(V1), sectors(V2)), dual = V1.dual)
    else
        throw(SpaceMismatch("Infimum of space and dual space does not exist"))
    end
end

function supremum(V1::GradedSpace{I}, V2::GradedSpace{I}) where {I<:Sector}
    if V1.dual == V2.dual
        typeof(V1)(c=>max(dim(V1, c), dim(V2, c)) for c in
            union(sectors(V1), sectors(V2)), dual = V1.dual)
    else
        throw(SpaceMismatch("Supremum of space and dual space does not exist"))
    end
end

function Base.show(io::IO, V::GradedSpace{I}) where {I<:Sector}
    print(io, type_repr(typeof(V)), "(")
    seperator = ""
    comma = ", "
    io2 = IOContext(io, :typeinfo => I)
    for c in sectors(V)
        if isdual(V)
            print(io2, seperator, dual(c), "=>", dim(V, c))
        else
            print(io2, seperator, c, "=>", dim(V, c))
        end
        seperator = comma
    end
    print(io, ")")
    V.dual && print(io, "'")
    return nothing
end

struct SpaceTable end
"""
    const Vect

A constant of a singleton type used as `Vect[I]` with `I<:Sector` a type of sector, to
construct or obtain the concrete type `GradedSpace{I,D}` instances without having to
specify `D`.
"""
const Vect = SpaceTable()
Base.getindex(::SpaceTable) = ComplexSpace
Base.getindex(::SpaceTable, ::Type{Trivial}) = ComplexSpace
function Base.getindex(::SpaceTable, I::Type{<:Sector})
    if Base.IteratorSize(values(I)) isa Union{HasLength, HasShape}
        N = length(values(I))
        return GradedSpace{I, NTuple{N, Int}}
    else
        return GradedSpace{I, SectorDict{I, Int}}
    end
end

Base.getindex(::ComplexNumbers, I::Type{<:Sector}) = Vect[I]
struct RepTable end
"""
    const Rep

A constant of a singleton type used as `Rep[G]` with `G<:Group` a type of group, to
construct or obtain the concrete type `GradedSpace{Irrep[G],D}` instances without having to
specify `D`. Note that `Rep[G] == Vect[Irrep[G]]`.

See also [`Irrep`](@ref) and [`Vect`](@ref).
"""
const Rep = RepTable()
Base.getindex(::RepTable, G::Type{<:Group}) = Vect[Irrep[G]]

type_repr(::Type{<:GradedSpace{I}}) where {I<:Sector} =
    "Vect[" * type_repr(I) * "]"
type_repr(::Type{<:GradedSpace{<:AbstractIrrep{G}}}) where {G<:Group} =
    "Rep[" * type_repr(G) * "]"
function type_repr(::Type{<:GradedSpace{ProductSector{T}}}) where
                                                        {T<:Tuple{Vararg{AbstractIrrep}}}
    sectors = T.parameters
    s = "Rep["
    for i in 1:length(sectors)
        if i != 1
            s *= " × "
        end
        s *= type_repr(supertype(sectors[i]).parameters[1])
    end
    s *= "]"
    return s
end

# Specific constructors for Z_N
const ZNSpace{N} = GradedSpace{ZNIrrep{N}, NTuple{N,Int}}
ZNSpace{N}(dims::NTuple{N, Int}; dual::Bool = false) where {N} = ZNSpace{N}(dims, dual)
ZNSpace{N}(dims::Vararg{Int, N}; dual::Bool = false) where {N} = ZNSpace{N}(dims, dual)
ZNSpace(dims::NTuple{N, Int}; dual::Bool = false) where {N} = ZNSpace{N}(dims, dual)
ZNSpace(dims::Vararg{Int, N}; dual::Bool = false) where {N} = ZNSpace{N}(dims, dual)

# TODO: Do we still need all of those
# ASCII type aliases
const ZNSpace{N} = GradedSpace{ZNIrrep{N}, NTuple{N,Int}}
const Z2Space = ZNSpace{2}
const Z3Space = ZNSpace{3}
const Z4Space = ZNSpace{4}
const U1Space = Rep[U₁]
const CU1Space = Rep[CU₁]
const SU2Space = Rep[SU₂]

# Unicode alternatives
const ℤ₂Space = Z2Space
const ℤ₃Space = Z3Space
const ℤ₄Space = Z4Space
const U₁Space = U1Space
const CU₁Space = CU1Space
const SU₂Space = SU2Space


text:
#=
 Explicit MDP model of an urban environment with one car
(extension of the single intersection model)
=#

# State type
struct CarMDPState
    crash::Bool
    ego::VehicleState
    car::VehicleState
    route::SVector{2, LaneTag}
end

# copy b to a
function Base.copyto!(a::CarMDPState, b::CarMDPState)
    a.crash = b.crash
    a.ego = b.ego
    a.car = b.car
    a.route = b.route
end

function Base.hash(s::CarMDPState, h::UInt64 = zero(UInt64))
    return hash(s.crash, hash(s.ego, hash(s.car, hash(s.route, h))))
end

function Base.:(==)(a::CarMDPState, b::CarMDPState)
    return a.crash == b.crash && a.ego == b.ego && a.car == b.car && a.route == b.route
end

# Action type
const CarMDPAction = UrbanAction

@with_kw mutable struct CarMDP <: MDP{CarMDPState, CarMDPAction}
    env::UrbanEnv = UrbanEnv(params=UrbanParams(nlanes_main=1,
                     crosswalk_pos =  [VecSE2(6, 0., pi/2), VecSE2(-6, 0., pi/2), VecSE2(0., -5., 0.)],
                     crosswalk_length =  [10.0, 10., 10.0],
                     crosswalk_width = [4.0, 4.0, 3.1],
                     stop_line = 22.0))
    ego_type::VehicleDef = VehicleDef()
    car_type::VehicleDef = VehicleDef()
    car_model::DriverModel = RouteFollowingIDM(route=random_route(env, Base.GLOBAL_RNG), σ=1.0)
    max_acc::Float64 = 2.0
    pos_res::Float64 = 3.0
    vel_res::Float64 = 2.0
    ego_start::Float64 = env.params.stop_line - ego_type.length/2
    ego_goal::LaneTag = LaneTag(2,1)
    off_grid::VecSE2 = VecSE2(UrbanEnv().params.x_min+VehicleDef().length/2, env.params.y_min+VehicleDef().width/2, 0)
    ΔT::Float64 = 0.5
    car_birth::Float64 = 0.3
    collision_cost::Float64 = -1.0
    action_cost::Float64 = 0.
    goal_reward::Float64 = 1.
    γ::Float64 = 0.95
    _ped_grid::Dict{LaneTag, RectangleGrid{3}} = init_ped_grid(env, pos_res, vel_ped_res)
    _car_grid::Dict{LaneTag, RectangleGrid{2}} = init_car_grid(env, pos_res, vel_res)
    _l_grid::Dict{LaneTag, RectangleGrid{1}} = init_l_grid(env, pos_res)
    _v_grid::RectangleGrid{1} = init_v_grid(env, vel_res)
    _collision_checker::Dict{Tuple{Vararg{VehicleState, 3}}, Bool} = Dict{Tuple{Vararg{VehicleState, 3}}, Bool}()
end

### REWARD MODEL ##################################################################################

function POMDPs.reward(mdp::CarMDP, s::CarMDPState, a::CarMDPAction, sp::CarMDPState)
    r = mdp.action_cost
    ego = sp.ego
    if sp.crash
        r += mdp.collision_cost
    end
    if ego.posF.s >= get_end(mdp.env.roadway[mdp.ego_goal]) &&
       get_lane(mdp.env.roadway, ego).tag == mdp.ego_goal
        r += mdp.goal_reward
    end
    return r
end

function POMDPs.isterminal(mdp::CarMDP, s::CarMDPState)
    if s.crash
        return true
    elseif s.ego.posF.s >= get_end(mdp.env.roadway[mdp.ego_goal]) &&
       get_lane(mdp.env.roadway, s.ego).tag == mdp.ego_goal
       return true
   end
   return false
end

## Helpers

function POMDPs.discount(mdp::CarMDP)
    return mdp.γ
end

POMDPs.actions(mdp::CarMDP) = [CarMDPAction(-4.0), CarMDPAction(-2.0), CarMDPAction(0.0), CarMDPAction(2.0)]

function POMDPs.actionindex(mdp::CarMDP, action::CarMDPAction)
    if action.acc == -4.0
        return 1
    elseif action.acc == -2.0
        return 2
    elseif action.acc == 0.
        return 3
    else
        return 4
    end
end

function get_off_the_grid(mdp::CarMDP)
    return VehicleState(mdp.off_grid, mdp.env.roadway, 0.)
end

function AutomotivePOMDPs.random_route(env::UrbanEnv, rng::AbstractRNG)
    start_lanes = get_start_lanes(env.roadway)
    ego_lanes = Set([env.roadway[tag] for tag in get_ego_route(env)])
    start_lanes = setdiff(start_lanes, ego_lanes)
    start_lane = rand(rng, start_lanes)
    return random_route(rng, env.roadway, start_lane)
end

# create labels for model checking
function labeling(mdp::CarMDP)
    labels = Dict{CarMDPState, Vector{String}}()
    for s in ordered_states(mdp)
        if s.crash
            labels[s] = ["crash"]
        elseif s.ego.posF.s >= get_end(mdp.env.roadway[mdp.ego_goal]) &&
            get_lane(mdp.env.roadway, s.ego).tag == mdp.ego_goal
            labels[s] = ["goal"]
        end
    end
    return labels
end

ext: jl

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
using RecipesBase, Unitful, Measurements

@recipe function f(spec::Spectrum)
    seriestype --> :step
    xlabel --> "wave"
    ylabel --> "flux density"
    label --> ""
    ustrip.(spec.wave), Measurements.value.(ustrip.(spec.flux))
end

@recipe function f(spec::EchelleSpectrum)
    seriestype --> :step
    xlabel --> "wave"
    ylabel --> "flux density"
    label --> ["Order $i" for i in 1:size(spec, 1)]
    spec.wave', spec.flux'
end


text:
using KeyTuples
using Base.Test
import DataFrames

macro capture_io(io, body)
    quote
        $io = IOBuffer()
        $(body)
        takebuf_string($io)
    end
end

k = KeyTuple(:a => 1, :b => "string")

@testset "constructors" begin
    @test k == KeyTuple((:a,:b),(1,"string"))
end

@testset "accessors" begin
    @test keys(k) ≡ (:a,:b)
    @test values(k) == (1,"string")
    @test k[:a] == 1
    @test k[:b] == "string"
    @test k[end] == "string"
    @test_throws BoundsError k[:c]
end

@testset "comparison" begin
    @test k == KeyTuple(:a => 1, :b => "string")
    @test k != KeyTuple(:a => 1, :c => "string")
    @test k != KeyTuple(:a => 1, :b => 2)
end

@testset "show" begin
    string1 = @capture_io io show(io, k)
    string2 = @capture_io io begin
        print(io, "KeyTuple")
        show(io, (:a => 1, :b => "string"))
    end
    @test string1 == string2
end

@testset "iteration" begin
    @test first(k) == (:a => 1)
    @test last(k) == (:b => "string")
    @test collect(k) == [:a => 1, :b => "string"]
end

@testset "dataframe conversion" begin
    kts = [KeyTuple(:a => i, :b => string(i)) for i in 1:5]
    df = keytuples_to_df(kts)
    df2 = DataFrames.DataFrame(a=collect(1:5), b=string.(1:5))
    @test isequal(df, df2)
end


text:
using NCDatasets

# issue #96

filename = tempname()

ds = Dataset(filename, "c")

data = defGroup(ds, "data")
st = defGroup(data, "stats")

defDim(data,"n_levels",10)
defDim(data,"n_fov",5)
defDim(data,"n_for",4)
defDim(data,"n_lines",2)

defVar(st, "air_temperature", Float32, ("n_levels","n_fov","n_for","n_lines"))

close(ds)

ds = Dataset(filename, "a")
data = ds.group["data"]
oe = defGroup(data, "optimal_estimation")

defVar(oe, "other_variable", Float32, ("n_levels","n_fov","n_for","n_lines"))

@test "other_variable" in keys(oe)

close(ds)

rm(filename)


text:
module DataAPI

"""
    defaultarray(T, N)

For a given element type `T` and number of dimensions `N`, return the appropriate array
type.

The default definition returns `Array{T, N}`. This function is useful for custom types that
have a more efficient vectorized representation (usually using SOA optimizations).

This generic function is owned by DataAPI.jl itself, which is the sole provider of the
default definition.
"""
function defaultarray end
defaultarray(::Type{T}, N) where {T} = Array{T, N}

"""
    refarray(A::AbstractArray)

For a given array `A`, potentially return an optimized "ref array" representation of the
original array, which can allow for faster comparison and sorting.

The default definition just returns the input array. This function is useful for custom
array types which already store a "hashed"-like representation of elements where testing
equality or permuting elements in place can be much faster than the original scalar value,
like pooled arrays.

This generic function is owned by DataAPI.jl itself, which is the sole provider of the
default definition.
"""
function refarray end
refarray(A::AbstractArray) = A

"""
    refvalue(A, x)

For the *original* array `A`, and a "ref value" `x` taken from `refarray(A)`, return the
appropriate *original* value. `refvalue(A, refarray(A)[I...])` must be equal to `A[I...]`.

By default, `refvalue(A, x)` returns `x` (since `refarray(A)` returns `A` by default).
This allows recovering an original array element after operating on the "ref values".

This generic function is owned by DataAPI.jl itself, which is the sole provider of the
default definition.
"""
function refvalue end
refvalue(A::AbstractArray, x) = x

"""
    refpool(A)

Whenever available, return an indexable object `pool` such that, given the *original* array `A` and
a "ref value" `x` taken from `refarray(A)`, `pool[x]` is the appropriate *original* value. Return
`nothing` if such object is not available. If `refpool(A)` is not `nothing`, then
`refpool(A)[refarray(A)[I...]]` must be equal to `A[I...]`.

By default, `refpool(A)` returns `nothing`.

This generic function is owned by DataAPI.jl itself, which is the sole provider of the
default definition.
"""
function refpool end
refpool(A::AbstractArray) = nothing

"""
    describe(io::IO, x)

For an object `x`, print descriptive statistics to `io`.

This generic function is owned by StatsBase.jl, which is the sole provider of the default
definition.
"""
function describe end

"""
    levels(x)

Return a vector of unique values which occur or could occur in collection `x`,
omitting `missing` even if present. Values are returned in the preferred order
for the collection, with the result of [`sort`](@ref) as a default.

Contrary to [`unique`](@ref), this function may return values which do not
actually occur in the data, and does not preserve their order of appearance in `x`.
"""
function levels(x)
    T = Base.nonmissingtype(eltype(x))
    levs = convert(AbstractArray{T}, filter!(!ismissing, unique(x)))
    try
        sort!(levs)
    catch
    end
    levs
end

"""
    Between(first, last)

Select the columns between `first` and `last` from a table.
"""
struct Between{T1 <: Union{Int, Symbol}, T2 <: Union{Int, Symbol}}
    first::T1
    last::T2
end

"""
    All(cols...)

Select the union of the selections in `cols`. If `cols == ()`, select all columns.
"""
struct All{T<:Tuple}
    cols::T
    All(args...) = new{typeof(args)}(args)
end

end # module


text:
@testset "DeviceModel Tests" begin
    @test_throws ArgumentError DeviceModel(PSY.ThermalGen, PSI.ThermalStandardUnitCommitment)
    @test_throws ArgumentError DeviceModel(PSY.ThermalStandard, PSI.AbstractDeviceFormulation)
    @test_throws ArgumentError PSI.CanonicalModel(JuMP.Model(),
                                                nothing,
                                                false,
                                                true,
                                                1:24,
                                                Dates.Hour(1),
                                                Dates.now(),
                                                Dict{Symbol, JuMP.Containers.DenseAxisArray}(),
                                                Dict{Symbol, JuMP.Containers.DenseAxisArray}(),
                                                JuMP.AffExpr(0),
                                                Dict{Symbol, JuMP.Containers.DenseAxisArray}(),
                                                nothing,
                                                Dict{PSI.ICKey, Array{InitialCondition}}(),
                                                nothing)
end

@testset "OperationModel Tests" begin
    for p in [true, false]
        t = OperationModel(TestOptModel, CopperPlatePowerModel, c_sys5; parameters = p)
        @test PSI.model_has_parameters(t.canonical) == p
    end
end


text:
using Distributions
using JSON

mu = -2
s = 3
d = Logistic( mu,s )

x = linspace( -0.3, 0.3, 400 )

dmgf(t) = mgf(d, t )
y = map( dmgf, x )
println( y )

data = Dict([
	("mu", mu),
	("s", s),
	("data", x),
	("expected", y)
])

outfile = open("./test/fixtures/typedarray.json", "w")
JSON.json(data)

write( outfile, JSON.json(data) )

max_stars_repo_name: distributions-io/logistic-mgf

max_issues_repo_name: distributions-io/logistic-mgf

max_forks_repo_name: distributions-io/logistic-mgf


text:
using Plots
plotlyjs()

# Rossler space equation
function rossler(t, state, a, b, c)
    x, y, z = state
    dx = -y - z
    dy = x + a*y
    dz = b + z*(x-c)
    return [dx, dy, dz]
end

"""
    rk4(f, x0, time_span[, h, args])

4-th order Runge-Kutta. ``f`` is a handler to the
space state funtion, ``x0`` is a vector of initial
conditions, ``time_span`` is a tuple containing
the start and end simulation times and ``h`` is the
integration interval.

The function ``f`` is called as:

    dx = f(t, x, args...)
"""
function rk4_old(f, x0, time_span; h=(time_span[2]-time_span[1])/1e4, args=())
    t = time_span[1]:h:time_span[2]
    N = length(t)
    n = length(x0)
    # Allocate buffers
    x = zeros(N, n)
    xd = zeros(4, n)
    # Simulate
    x[1, :] = x0
    for i = 1:N-1
        # 1-st evaluation
        xd[1, :] = f(t[i], x[i, :], args...)
        # 2-nd evaluation
        xd[2, :] = f(t[i]+h/2, x[i, :]+h/2*xd[1, :], args...)
        # 3-rd evaluation
        xd[3, :] = f(t[i]+h/2, x[i, :]+h/2*xd[2, :], args...)
        # 4-th evaluation
        xd[4, :] = f(t[i]+h, x[i, :]+h*xd[3, :], args...)
        # Next point
        x[i+1, :] = x[i, :]
        x[i+1, :] += (xd[1, :] + 2*xd[2, :] + 2*xd[3, :] + xd[4, :]) * h/6
    end
    return x
end

# Simulate Rossler
time_span = (0, 1000)
x0 = [1.0, 1.0, 2.0]
args = (0.2, 0.2, 5.7)

@time x = rk4_old(rossler, x0, time_span, h=0.001, args=args)

plot(x)

plot(x[:, 1], x[:,2], x[:, 3])

max_stars_repo_name: MACSIN-UFMG/julia-rossler-example

max_issues_repo_name: MACSIN-UFMG/julia-rossler-example

max_forks_repo_name: MACSIN-UFMG/julia-rossler-example


text:
# abstract definition
abstract type AbstractMarkovDecisionProcess end

## some constants
# transition probs
const αts = 1/2
const αtw = 4/7
const αew = 4/5
const αewo = 1/2
# transition rewards
const rtse = -0.1
const rtst = -0.1
const rtwt = 1.0
const rtwoe = -0.5
const rtwoet = -0.3
const rewt = 1.4
const rewe = 1.6
const rewoh = -0.6
const rewoe = -0.2
const rhwt = 10

# our problem
struct HumanMDP <: AbstractMarkovDecisionProcess
    initial_state::String
    states::Set{String}
    rewards::Dict
    transitions::Dict
    γ::Float64
    function HumanMDP(initial_state::String, γ::Float64)
        # states (s)
        states = Set(["tired", "energized","healthier"])
        # transition
        trans = Dict(
            "tired" => Dict("sleep"=> [(αts,"energized"),(1-αts,"tired")],
                            "work" => [(1, "tired")],
                            "workout" => [(αtw,"energized"),(1-αtw,"tired")]),
            "energized" => Dict("work"=> [(αew,"tired"),(1-αew,"energized")],
                                "workout"=> [(αewo,"healthier"),(1-αewo,"energized")]),
            "healthier" => Dict("work"=> [(1, "tired")]))
        # rewards R(s,a,s')
        rewards = Dict(
            "tired" => Dict("sleep"=> [(rtse,"energized"),(rtst,"tired")],
                            "work" => [(rtwt, "tired")],
                            "workout" => [(rtwoe,"energized"),(rtwoet,"tired")]),
            "energized" => Dict("work"=> [(rewt,"tired"),(rewe,"energized")],
                                "workout"=> [(rewoh,"healthier"),(rewoe,"energized")]),
            "healthier" => Dict("work"=> [(rhwt, "tired")]))
        return new(initial_state, states, rewards, trans, γ)
    end
end

function actions(mdp::HumanMDP, state)
    return keys(mdp.transitions[state])
end

function transition_model(mdp::HumanMDP, state::String, action::String)
    return mdp.transitions[state][action]
end

function reward(mdp::HumanMDP, state::String, action::String, ns::String)
    idx = findfirst(x-> x[2]==ns, mdp.rewards[state][action])
    return mdp.rewards[state][action][idx][1]
end

"""
    value_iteration(mdp::HumanMDP; epsilon::Float64=0.001)

Return the utilities of the MDP's states as a Dict by applying the value iteration algorithm
on the given Markov decision process 'mdp' and a arbitarily small positive number 'epsilon'.

Modified version that works with R(s,a,s')
"""
function value_iteration(mdp::HumanMDP; ϵ::Float64=0.001)
    # initialise utility to zero for all states
    U_prime = Dict(Pair(state, 0.0) for state in mdp.states)
    while true
        U = copy(U_prime)
        δ = 0.0
        # go over all possible state
        for state in mdp.states
            # possible outcomes: Array[ Pair(a => [(T(s,a,s', s'), (T(s,a,s', s')]) ]
            outcomes = [Pair(action,transition_model(mdp, state, action)) for action in actions(mdp, state)]
            expectations = map(x-> sum( T*(reward(mdp, state, x.first, ns) + mdp.γ * U[ns]) for (T,ns) in x.second), outcomes)
            U_prime[state] = maximum(expectations)
            δ = max(δ, abs(U_prime[state] - U[state]))
        end
        if (δ < ((ϵ * (1 - mdp.γ))/mdp.γ))
            return U_prime
        end
    end
end

"""
    policy_extraction(mdp::HumanMDP; epsilon::Float64=0.001)

Return the optimal policy given the computed utilities

Modified version that works with R(s,a,s')
"""
function policy_extraction(mdp::HumanMDP, U::Dict)
	pi = Dict()
	for state in mdp.states
        # possible outcomes: Array[ Pair(a => [(T(s,a,s', s'), (T(s,a,s', s')]) ]
        outcomes = [Pair(action,transition_model(mdp, state, action)) for action in actions(mdp, state)]
        expectations = map(x-> sum( T*(reward(mdp, state, x.first, ns) + mdp.γ * U[ns]) for (T,ns) in x.second), outcomes)
		pi[state] = outcomes[argmax(expectations)].first
	end
	return pi
end

h = HumanMDP("tired", 0.)
s = "tired"
a = "workout"
ns = "tired"
actions(h,s)
transition_model(h, s, a)
reward(h, s, a, ns)

# γ = 0 => only look at immediate (local) gains
@info "running the model using γ = 0"
h = HumanMDP("tired", 0.)
V_star = value_iteration(h)
π_star = policy_extraction(h,V_star)
display(π_star)
# γ > 0 => also consider neighbors
γ = 0.85
@info "running the model using γ = $(γ)"
H = HumanMDP("tired", γ)
VV_star = value_iteration(H)
ππ_star = policy_extraction(H,VV_star)
display(ππ_star)


text:
mutable struct LL_approx
    models::Vector{AbstractModel}
    λ::Vector{Float64}
    indexes::Matrix{Bool}

    X_train::Matrix{Float64}
    Y_train::Matrix{Float64}

    X_test::Matrix{Float64}
    Y_test::Matrix{Float64}

    ys::Vector{Function}

    approx_errors::Vector{Float64}
    scale_values::Matrix{Float64}

end

function LL_approx( X_train,
                    Y_train;
                    models = AbstractModel[],
                    λ = [1e-5 for i=1:size(Y_train, 2)],
                    indexes = zeros(Bool, 0, 0),
                    X_test  = zeros(Float64, 0, 0),
                    Y_test  = zeros(Float64, 0, 0),
                    ys = Function[],
                    p_train = 0.8,
                    approx_errors = [Inf],
                    scale_values = ones(1,size(X_train, 2)))

    if size(X_train, 1) != size(Y_train, 1)
        @error("Check size of the train set (X_train, Y_train)")
    end
    if length(size(scale_values)) == 1
        scale_values = Array(scale_values')
    end

    X_train = X_train .* scale_values

    if (isempty(X_test) || isempty(Y_test))
        n = round(Int, p_train*size(X_train, 1))
        idx = randperm(size(X_train, 1))

        X_test = X_train[idx[n+1:end], :]
        Y_test = Y_train[idx[n+1:end], :]

        X_train = X_train[idx[1:n], :]
        Y_train = Y_train[idx[1:n], :]
    else
        X_test = X_test .* scale_values
    end

    LL_approx(models, λ, indexes, X_train, Y_train, X_test, Y_test, ys,
                approx_errors, scale_values)
end

function fit!(ll_approx::LL_approx)
    empty!(ll_approx.models)
    empty!(ll_approx.ys)

    for i = 1:size(ll_approx.Y_train, 2)
        Yx = view(ll_approx.Y_train, :, i)
        if sum(ll_approx.indexes[i,:]) == 0 || var(Yx) < 1e-8
            y =  x -> mean(Yx)
            model = KernelInterpolation(zeros(1), zeros(1,1))
        else
            Xtr = view(ll_approx.X_train, :, ll_approx.indexes[i,:])
            model = KernelInterpolation(Yx, Xtr, λ = ll_approx.λ[i], kernel=PeriodicKernel())
            train!(model)
            y = approximate(model)

        end
        # approximate(m) returns a function ŷ_i
        push!(ll_approx.models, model)
        push!(ll_approx.ys, y)

    end
end

function update_approximation_errors!(ll_approx)
    empty!(ll_approx.approx_errors)

    s = zeros(size(ll_approx.Y_train, 2))
    for i = 1:size(ll_approx.X_test, 1)
        y_approx = Ψ_approx(view(ll_approx.X_test, i,:), ll_approx, false)
        s += (y_approx - view(ll_approx.Y_test, i, :)).^2
    end

    ll_approx.approx_errors = sqrt.(s)

end

function optimize_indexes!(model::LL_approx, debug=false)
    # D_ul = size(model.X_test, 1)
    # if D_ul > 10
        # indexes_initial = rand(Bool[0, 1], size(X_train, 2))
        # indexes_new, res = hillClimbing(ff, indexes_initial, iters=100, debug=true, r = 2)
    # else
    # end
    bruteforce!(model, debug=debug)

end

function Ψ_approx(x, ll_approx, scale_data=true)
    y = zeros(length(ll_approx.ys))

    if scale_data
        x = x .* ll_approx.scale_values
    end
    # x = xx[1,:]

    for i = 1:length(ll_approx.ys)
        y[i] = ll_approx.ys[i](x[view(ll_approx.indexes, i,:)])
    end
    y
end

function bruteforce!(ll_approx; debug=false)
    D = size(ll_approx.X_train, 2)
    D_ll = size(ll_approx.Y_train, 2)

    approx_errors_old = Inf*ones(D_ll)
    indexes_best = zeros(D_ll, D)

    for t = 1:2^D
        x = parse.(Int, collect(bitstring(t)), base=10)
        reverse!(x)

        x_new = Bool.(x[1:D])
        indexes = repeat(x_new', D_ll)

        ll_approx.indexes = indexes

        fit!(ll_approx)
        update_approximation_errors!(ll_approx)

        for i = 1:D_ll
            if ll_approx.approx_errors[i] < approx_errors_old[i]
                approx_errors_old[i] = ll_approx.approx_errors[i]
                indexes_best[i,:] = indexes[i,:]
            end
        end

        if debug
            println(">>>>>>>>>>>> iter: ", t)
            display(indexes_best)
            println("err: ", approx_errors_old)
        end

    end

    ll_approx.indexes = indexes_best

    fit!(ll_approx)
    update_approximation_errors!(ll_approx)

    indexes_best
end


text:
#=--------------------------------------------------------------------

cornerfeatures - Functions for the detection of 'corner' features

Copyright (c) 2015-2016 Peter Kovesi
peterkovesi.com

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

The Software is provided "as is", without warranty of any kind.

PK August 2015
   September 2018 updated for v0.7/v1.0

---------------------------------------------------------------------=#

export structuretensor, harris, noble, shi_tomasi, coherence
export hessianfeatures
export fastradial

import LinearAlgebra, Images

#----------------------------------------------------------------------
"""
harris - Harris corner detector
```
Usage:                cimg = harris(img, sigma=1; k=0.04)
              (cimg, r, c) = harris(img, sigma=1; k=0.04, args...)

Arguments:
           img    - Image to be processed.
                    ::Array{T,2} or ::ImageMeta{T,2}
           sigma  - Standard deviation of Gaussian summation window.
                    Typical values to use might be 1-3. Default is 1.

Keyword arguments:
                k - The scaling factor to apply to the trace of the
                    structure tensor. Defaults to the traditional value
                    of 0.04

          args... - A sequence of keyword arguments for performing non-maximal suppression.
                    These are: (see nonmaxsuppts() for more details)
           radius - Radius of region considered in non-maximal
                    suppression. Typical values to use might
                    be 1-3 pixels. Default (and minimum value) is 1.
           thresh - Threshold, only features with value greater than
                    threshold are returned. Default is 0.
                N - Maximum number of features to return.  In this case the
                    N strongest features with value above 'thresh' are
                    returned. Default is Inf.
         subpixel - If set to true features are localised to subpixel
                    precision. Default is false.
             img  - Optional image data.  If an image is supplied the
                    detected corners are overlayed on this image. This can
                    be useful for parameter tuning. Default is [].

Returns:
           cimg   - Corner strength image.
           r      - Row coordinates of corner points.
           c      - Column coordinates of corner points.
```
If the input is an ImageMeta cimg will be an ImageMeta with
property spatialorder set to ["y","x"] so that it is consistent with
the coordinates returned in r, c

With none of the optional nonmaximal keyword arguments specified only
`cimg` is returned as a raw corner strength image.  You may then want
to look at the values within 'cimg' to determine the appropriate
threshold value to use. Note that the Harris corner strength varies
with the intensity gradient raised to the 4th power.  Small changes in
input image contrast result in huge changes in the appropriate
threshold.

The Harris measure is det(M) - k*trace^2(M), where k is a parameter
you have to set (traditionally k = 0.04) and M is the structure
tensor.  Use Noble's measure if you wish to avoid the need to set a
parameter k.  However the Shi-Tomasi measure is probably what you
really want.

See also: noble(), shi_tomasi(), hessianfeatures(), nonmaxsuppts(), derivative5()
"""
function  harris(img::Array{T,2}, sigma::Real=1; k::Real=0.04, args...) where T <: Real
    #=
    References:
    C.G. Harris and M.J. Stephens. "A combined corner and edge detector",
    Proceedings Fourth Alvey Vision Conference, Manchester.
    pp 147-151, 1988.
    =#

    (Ix2, Iy2, Ixy) = structuretensor(img, sigma)
    cimg = (Ix2.*Iy2 .- Ixy.^2) .- k*(Ix2 + Iy2).^2

    if isempty(args)
        return cimg
    else
        (r,c) = nonmaxsuppts(cimg; args...)
        return cimg, r, c
    end

end

# Case for img::ImageMeta
function harris(img::Images.ImageMeta{T,2}, sigma::Real=1; k::Real=0.04, args...) where T <: Real

    # Extract 2D float array in y x spatial order
    (dimg, prop) = floatyx(img)

    if isempty(args)
        cimg = harris(dimg, sigma, k=k)
        return Images.Image(cimg, prop)
    else
        (cimg, r, c) = harris(dimg, sigma, k=k; args...)
        return Images.Image(cimg, prop), r, c
    end
end

#----------------------------------------------------------------------
"""
noble - Noble's variant of the Harris corner detector
```

Usage:                cimg = noble(img, sigma=1)
              (cimg, r, c) = noble(img, sigma=1; args...)

Arguments:
           img    - Image to be processed.
                    ::Array{T,2} or ::ImageMeta{T,2}
           sigma  - Standard deviation of Gaussian summation window.
                    Typical values to use might be 1-3. Default is 1.
Keyword arguments:
          args... - A sequence of keyword arguments for performing non-maximal suppression.
                    These are: (see nonmaxsuppts() for more details)
           radius - Radius of region considered in non-maximal
                    suppression. Typical values to use might
                    be 1-3 pixels. Default (and minimum value) is 1.
           thresh - Threshold, only features with value greater than
                    threshold are returned. Default is 0.
                N - Maximum number of features to return.  In this case the
                    N strongest features with value above 'thresh' are
                    returned. Default is Inf.
         subpixel - If set to true features are localised to subpixel
                    precision. Default is false.
             img  - Optional image data.  If an image is supplied the
                    detected corners are overlayed on this image. This can
                    be useful for parameter tuning. Default is [].

Returns:
           cimg   - Corner strength image.
           r      - Row coordinates of corner points.
           c      - Column coordinates of corner points.
```
If the input is an ImageMeta cimg will be an ImageMeta with
property spatialorder set to ["y","x"] so that it is consistent with
the coordinates returned in r, c

With none of the optional keyword arguments specified only `cimg` is
returned as a raw corner strength image.  You may then want to look at
the values within 'cimg' to determine the appropriate threshold value
to use.

Noble's variation of the Harris operator avoids the need to set the
parameter `k`.  However the Shi-Tomasi measure is probably what you
really want.

See also: harris(), shi_tomasi(), hessianfeatures(), nonmaxsuppts(),
structuretensor(), derivative5()
"""
function  noble(img::Array{T,2}, sigma::Real=1; args...) where T <: Real
    #=
    References:
    C.G. Harris and M.J. Stephens. "A combined corner and edge detector",
    Proceedings Fourth Alvey Vision Conference, Manchester.
    pp 147-151, 1988.

    Alison Noble, "Descriptions of Image Surfaces", PhD thesis, Department
    of Engineering Science, Oxford University 1989, p45.
    =#

    (Ix2, Iy2, Ixy) = structuretensor(img, sigma)
    cimg = (Ix2.*Iy2 .- Ixy.^2)./(Ix2 .+ Iy2 .+ eps())

    if !isempty(args)
        (r,c) = nonmaxsuppts(cimg; args...)
        return cimg, r, c
    else
        return cimg
    end

end

# Case for img::ImageMeta
function noble(img::Images.ImageMeta{T,2}, sigma::Real=1; args...) where T <: Real

    # Extract 2D float array in y x spatial order
    (dimg, prop) = floatyx(img)

    if isempty(args)
        cimg = noble(dimg, sigma)
        return Images.Image(cimg, prop)
    else
        (cimg, r, c) = noble(dimg, sigma; args...)
        return Images.Image(cimg, prop), r, c
    end
end

#----------------------------------------------------------------------
"""
shi_tomasi - Shi - Tomasi corner detector
```
Usage:               cimg = shi_tomasi(img, sigma=1)
             (cimg, r, c) = shi_tomasi(img, sigma=1; args...)

Arguments:
           img    - Image to be processed.
                    ::Array{T,2} or ::ImageMeta{T,2}
           sigma  - Standard deviation of Gaussian summation window.
                    Typical values to use might be 1-3. Default is 1.
Keyword arguments:
          args... - A sequence of keyword arguments for performing non-maximal suppression.
                    These are: (see nonmaxsuppts() for more details)
           radius - Radius of region considered in non-maximal
                    suppression. Typical values to use might
                    be 1-3 pixels. Default (and minimum value) is 1.
           thresh - Threshold, only features with value greater than
                    threshold are returned. Default is 0.
                N - Maximum number of features to return.  In this case the
                    N strongest features with value above 'thresh' are
                    returned. Default is Inf.
         subpixel - If set to true features are localised to subpixel
                    precision. Default is false.
             img  - Optional image data.  If an image is supplied the
                    detected corners are overlayed on this image. This can
                    be useful for parameter tuning. Default is [].

Returns:
           cimg   - Corner strength image.
           r      - Row coordinates of corner points.
           c      - Column coordinates of corner points.
```
If the input is an ImageMeta cimg will be an ImageMeta with
property spatialorder set to ["y","x"] so that it is consistent with
the coordinates returned in r, c

With no keyword arguments supplied only 'cimg' is returned as a raw
corner strength image.  You may then want to look at the values within
'cimg' to determine the appropriate threshold value to use.

The Shi - Tomasi measure returns the minimum eigenvalue of the
structure tensor.  This represents the ideal that the Harris and Noble
detectors attempt to approximate.  Back in 1988 Harris wanted to avoid
the computational cost of taking a square root when computing
features.  This is no longer relevant today!

See also: harris(), noble(), hessianfeatures(), structuretensor(),
nonmaxsuppts(), derivative5()
"""
function shi_tomasi(img::Array{T1,2}, sigma::Real=1; args...) where T1 <: Real
    #=
    Reference:
    J. Shi and C. Tomasi. "Good Features to Track,". 9th IEEE Conference on
    Computer Vision and Pattern Recognition. 1994.
    =#

    (Ix2, Iy2, Ixy) = structuretensor(img, sigma)

    T = Ix2 .+ Iy2                 # trace
    D = Ix2.*Iy2 .- Ixy.^2         # determinant

    # The two eigenvalues of the 2x2 structure tensor are:
    # L1 = T/2 + sqrt(T.^2/4 - D)
    # L2 = T/2 - sqrt(T.^2/4 - D)

    # We just want the minimum eigenvalue
    cimg = T/2 .- sqrt.(T.^2/4 .- D .+ eps())

    if !isempty(args)
        (r,c) = nonmaxsuppts(cimg; args...)
        return cimg, r, c
    else
        return cimg
    end

end

# Case for img::ImageMeta
function shi_tomasi(img::Images.ImageMeta{T,2}, sigma::Real=1; args...) where T <: Real

    # Extract 2D float array in y x spatial order
    (dimg, prop) = floatyx(img)

    if isempty(args)
        cimg = shi_tomasi(dimg, sigma)
        return Images.Image(cimg, prop)
    else
        (cimg, r, c) = shi_tomasi(dimg, sigma; args...)
        return Images.Image(cimg, prop), r, c
    end
end

#----------------------------------------------------------------------
"""
coherence - Compute image coherence from structure tensor
```
Usage:      cimg = coherence(img, sigma=1)

Arguments:
           img    - Image to be processed.
                    ::Array{T,2} or ::ImageMeta{T,2}
           sigma  - Standard deviation of Gaussian summation window.
                    Typical values to use might be 1-3. Default is 1.
Returns:
           cimg   - Coherence image.
```

Coherence is defined as the difference between the eigenvalues of the
structure tensor divided by the sum of the eigenvalues, all squared.
```
     ((L1-L2)./(L1+L2)).^2
```
A small value indicates that the eigenvalues are similar in magnitude
and hence the local structure has a strong 2D component.

If the input is an ImageMeta cimg will be an ImageMeta with
property spatialorder set to ["y","x"].

See also: structuretensor()
"""
function coherence(img::Array{T1,2}, sigma::Real=1) where T1 <: Real

    (Ix2, Iy2, Ixy) = structuretensor(img, sigma)

    T = Ix2 .+ Iy2                 # trace
    D = Ix2.*Iy2 .- Ixy.^2         # determinant

    # The two eigenvalues of the 2x2 structure tensor are:
    L1 = T/2 .+ sqrt(T.^2/4 .- D)
    L2 = T/2 .- sqrt(T.^2/4 .- D)

    # Coherence is defined as ((L1-L2)./(L1+L2)).^2
    # this can be reduced to the following
    return 4*(T.^2/4 .- D) ./ T.^2

end

# Case for img::ImageMeta
function coherence(img::Images.ImageMeta{T,2}, sigma::Real=1) where T <: Real

    # Extract 2D float array in y x spatial order
    (dimg, prop) = floatyx(img)

    cimg = coherence(dimg, sigma)
    return Images.Image(cimg, prop)
end

#----------------------------------------------------------------------
"""
structuretensor - Compute structure tensor values over an image
```
Usage:  (Ix2, Iy2, Ixy) = structuretensor(img, sigma=1)

Arguments:
           img   - Image to be processed  ::Array{T,2}
           sigma - Standard deviation of Gaussian summation window. Typical
                   values to use might be 1-3. Default is 1.
Returns:
           Ix2   - 2nd moment of gradients in x
           Iy2   - 2nd moment of gradients in y
           Ixy   - 2nd moment of gradients in xy
```
See also: shi_tomasi(), harris(), noble(), coherence(), derivative5()
"""
function structuretensor(img::Array{T,2}, sigma::Real=1) where T <: Real

    # Convert to float if needed
    if isa(img[1,1], Integer)
        fimg = float(img)
    else
        fimg = img
    end

    # Compute derivatives and elements of the structure tensor.
    (Ix, Iy) = derivative5(fimg, ("x", "y"))
    Ix2 = gaussfilt(Ix.^2,  sigma)
    Iy2 = gaussfilt(Iy.^2,  sigma)
    Ixy = gaussfilt(Ix.*Iy, sigma)

    return Ix2, Iy2, Ixy
end

#----------------------------------------------------------------------
"""
hessianfeatures  - Computes determiant of hessian features in an image.
```
Usage: hdet = hessianfeatures(img, sigma)

Arguments:
            img      - Greyscale image to be processed.
                       ::Array{T,2} or ::ImageMeta{T,2}
            sigma    - Defines smoothing scale.

Returns:    hdet     - Matrix of determinants of Hessian
```
The local maxima of hdet tend to mark the centres of dark or light blobs.
However, the point that gets localised can be dependent on scale.  If the
blobs you wish to detect are large you will need to use a value of sigma that
is comparable in magnitude.

The local minima of hdet is useful for marking the intersection points of a
camera calibration checkerboard pattern.  These saddle features seem to be
more stable under scale variations.

For example to get the 100 strongest saddle features in image `img` use:
```
 > hdet = hessianfeatures(img, 1)      # sigma = 1
 > (r, c) = nonmaxsuppts(-hdet, N=100)
```
If the input is an ImageMeta cimg will be an ImageMeta with
property spatialorder set to ["y","x"].

See also: harris(), noble(), shi_tomasi(), derivative5(), nonmaxsuppts()
"""
function  hessianfeatures(img::Array{T,2}, sigma::Real=1) where T <: Real

    if sigma > 0    # Convolve with Gaussian at desired sigma
        Gimg = gaussfilt(img,  sigma)
    else           # No smoothing
        Gimg = img
        sigma = 1  # Needed for normalisation later
    end

    # Take 2nd derivatives in x and y
    (Lxx, Lxy, Lyy) = derivative5(Gimg, ("xx", "xy", "yy"))

    # Apply normalizing scaling factor of sigma^2 to 2nd derivatives
    Lxx = Lxx*sigma^2
    Lyy = Lyy*sigma^2
    Lxy = Lxy*sigma^2

    # Determinant
    return hdet = Lxx.*Lyy .- Lxy.^2
end

# Case for img::ImageMeta
function hessianfeatures(img::Images.ImageMeta{T,2}, sigma::Real=1) where T <: Real

    # Extract 2D float array in y x spatial order
    (dimg, prop) = floatyx(img)

    hdet = hessianfeatures(dimg, sigma)
    return Images.Image(cimg, prop)
end

#----------------------------------------------------------------------
"""
fastradial - Loy and Zelinski's fast radial feature detector
```
Usage: (S, So) = fastradial(img, radii, alpha=2, beta=0)

Arguments:
           img   - Image to be analysed
                   ::Array{T,2} or ::ImageMeta{T,2}
           radii - Vector of integer radius values to be processed
                   suggested radii might be [1, 3, 5]
           alpha - Radial strictness parameter.
                   1 - slack, accepts features with bilateral symmetry.
                   2 - a reasonable compromise.
                   3 - strict, only accepts radial symmetry.
                       ... and you can go higher
           beta  - Gradient threshold.  Gradients below this threshold do
                   not contribute to symmetry measure, defaults to 0.

Returns    S     - Symmetry map.  Bright points with high symmetry are
                   marked with large positive values. Dark points of
                   high symmetry marked with large -ve values.
           So    - Symmetry map based on orientation only.
```
To localize points use nonmaxsuppts() on S, -S or abs(S) depending on
what you are seeking to find.

A good algorithm for detecting eyes in face images.

If the input is an ImageMeta cimg will be an ImageMeta with
property spatialorder set to ["y","x"].

Reference:
Loy, G.  Zelinsky, A.  Fast radial symmetry for detecting points of
interest.  IEEE PAMI, Vol. 25, No. 8, August 2003. pp 959-973.
"""
function fastradial(img::Array{T,2}, radii::Vector, alpha::Real=2, beta::Real=0) where T <: Real

    if !isa(radii, Vector{Int}) && !isa(radii, Int) || minimum(radii) < 1
        error("radii must be a vector of integers and > 1")
    end

    (rows,cols) = size(img)

    # Compute derivatives in x and y via Farid and Simoncelli's 5 tap
    # derivative filters
    (imgx, imgy) = derivative5(img, ("x", "y"))
    mag = sqrt.(imgx.^2 .+ imgy.^2) .+ eps() # (+eps to avoid division by 0)

    # Normalise gradient values so that [imgx imgy] form unit
    # direction vectors.
    imgx = imgx./mag
    imgy = imgy./mag

    S = zeros(rows,cols)  # Symmetry matrix
    So = zeros(rows,cols) # Orientation only symmetry matrix

    for n in radii
        M = zeros(rows,cols)  # Magnitude projection image
        O = zeros(rows,cols)  # Orientation projection image

        # Form the orientation and magnitude projection matrices
        for r = 1:rows
            for c = 1:cols
                if mag[r,c] > beta

                    # Coordinates of 'positively' and 'negatively' affected pixels
                    posx = c + round(Int, n*imgx[r,c])
                    posy = r + round(Int, n*imgy[r,c])

                    negx = c - round(Int, n*imgx[r,c])
                    negy = r - round(Int, n*imgy[r,c])

                    # clamp to image limits
                    posx = max(min(posx,cols),1)
                    posy = max(min(posy,rows),1)
                    negx = max(min(negx,cols),1)
                    negy = max(min(negy,rows),1)

                    O[posy, posx] += 1
                    O[negy, negx] -= 1

                    M[posy, posx] += mag[r,c]
                    M[negy, negx] -= mag[r,c]
                end
            end
        end

        # Clamp Orientation projection matrix values to a maximum of
        # +/-kappa,  but first set the normalization parameter kappa to the
        # values suggested by Loy and Zelinski
        if n == 1
            kappa = 8
        else
            kappa = 9.9
        end

        O[O .>  kappa] .=  kappa
        O[O .< -kappa] .= -kappa

        # Unsmoothed symmetry measure at this radius value
        F = M./kappa .* (abs.(O)/kappa).^alpha
        Fo = sign.(O) .* (abs.(O)/kappa).^alpha   # Orientation only based measure

        # Smooth and spread the symmetry measure with a Gaussian
        # proportional to n.  Also scale the smoothed result by n so
        # that large scales do not lose their relative weighting.
        # Note we use the function gaussfilt because it seems that
        # small sigma values (< 1) cause problems for
        # Images.imfilt_gaussian()
        # ** But if sigma < 1 do we need to smooth anyway? **
        S  +=  gaussfilt(F,  0.25*n) * n
        So +=  gaussfilt(Fo, 0.25*n) * n

    end  # for each radius

    S  = S /length(radii)  # Average
    So = So/length(radii)

    return S, So
end

# Case for img::ImageMeta
function fastradial(img::Images.ImageMeta{T,2}, radii::Vector, alpha::Real=2, beta::Real=0) where T <: Real

    # Extract 2D float array in y x spatial order
    (dimg, prop) = floatyx(img)

    (S, So) = fastradial(dimg, radii, alpha, beta)
    return Images.Image(S, prop), Images.Image(So, prop)
end


text:
## Initializes GUI
# Toolbar
win1 = Window("GUI", 450, 200)
t = Grid()
bz = Button("Zoom Mode")
t[1,1] = bz
br = Button("Rectangle Mode")
t[1,2] = br
bf = Button("Freehand Mode")
t[1,3] = bf
bp = Button("Polygon Mode")
t[1,4] = bp
g = Grid()
g[1,1] = t

# Button actions for mode switches
zoomid = signal_connect(bz, "clicked") do z
  set_mode_all(zoom_mode)
end

rectid = signal_connect(br, "clicked") do r
  set_mode_all(rectangle_mode)
end

freeid = signal_connect(bf, "clicked") do f
  set_mode_all(freehand_mode)
end

polyid = signal_connect(bp, "clicked") do p
  set_mode_all(polygon_mode)
end

nb = Notebook()
g[2,1] = nb
# Measurements show for every open window
function display_measure(ctx::ImageContext)
  txt = TextView() # turn off ability for user to enter text
  setproperty!(txt,:expand,true)
  tblsig = map(ctx.canvas.mouse.buttonrelease) do btn
    df = init_measure_tbl()
    push!(df, measure(ctx.image))
    push!(df, measure(value(ctx.rectview)))
    #@show df
    df
  end
  push!(txt,string(value(tblsig))) # make this update w/ every buttonrelease
  push!(nb, txt, "Context $(ctx.id)")
  showall(win1)
  append!(ctx.canvas.preserved, [tblsig])
end

push!(win1,g)
showall(win1)


text:
using Test

# Customize the testset slightly to line up the output.
struct MBTestSet <: Test.AbstractTestSet end
MBTestSet(description) = Test.DefaultTestSet(rpad(description, 50))

include("tutorial.jl")


text:
# Sieve of Atkin for generating primes:
#     http://en.wikipedia.org/wiki/Sieve_of_Atkin
# Code very loosely based on this:
#     http://thomasinterestingblog.wordpress.com/2011/11/30/generating-primes-with-the-sieve-of-atkin-in-c/
#     http://dl.dropboxusercontent.com/u/29023244/atkin.cpp
#
function primesmask(s::AbstractVector{Bool})
    n = length(s)
    n < 2 && return s; s[2] = true
    n < 3 && return s; s[3] = true
    r = ifloor(sqrt(n))
    for x = 1:r
        xx = x*x
        for y = 1:r
            yy = y*y
            i, j, k = 4xx+yy, 3xx+yy, 3xx-yy
            i <= n && (s[i] $= (i%12==1)|(i%12==5))
            j <= n && (s[j] $= (j%12==7))
            1 <= k <= n && (s[k] $= (x>y)&(k%12==11))
        end
    end
    for i = 5:r
        s[i] && (s[i*i:i*i:n] = false)
    end
    return s
end
primesmask(n::Int) = primesmask(falses(n))
primesmask(n::Integer) = n <= typemax(Int) ? primesmask(int(n)) :
    error("you want WAY too many primes ($n)")

primes(n::Union(Integer,AbstractVector{Bool})) = find(primesmask(n))

# Miller-Rabin for primality testing:
#     http://en.wikipedia.org/wiki/Miller–Rabin_primality_test
#
function isprime(n::Integer)
    n == 2 && return true
    (n < 2) | iseven(n) && return false
    s = trailing_zeros(n-1)
    d = (n-1) >>> s
    for a in witnesses(n)
        a < n || break
        x = powermod(a,d,n)
        x == 1 && continue
        t = s
        while x != n-1
            (t-=1) <= 0 && return false
            x = oftype(n, widemul(x,x) % n)
            x == 1 && return false
        end
    end
    return true
end

# Miller-Rabin witness choices based on:
#     http://mathoverflow.net/questions/101922/smallest-collection-of-bases-for-prime-testing-of-64-bit-numbers
#     http://primes.utm.edu/prove/merged.html
#     http://miller-rabin.appspot.com
#
witnesses(n::Union(Uint8,Int8,Uint16,Int16)) = (2,3)
witnesses(n::Union(Uint32,Int32)) = n < 1373653 ? (2,3) : (2,7,61)
witnesses(n::Union(Uint64,Int64)) =
        n < 1373653         ? (2,3) :
        n < 4759123141      ? (2,7,61) :
        n < 2152302898747   ? (2,3,5,7,11) :
        n < 3474749660383   ? (2,3,5,7,11,13) :
                              (2,325,9375,28178,450775,9780504,1795265022)

isprime(n::Uint128) =
    n <= typemax(Uint64) ? isprime(uint64(n)) : isprime(BigInt(n))
isprime(n::Int128) = n < 2 ? false :
    n <= typemax(Int64)  ? isprime(int64(n))  : isprime(BigInt(n))

# TODO: faster factorization algorithms?

const PRIMES = primes(10000)

function factor{T<:Integer}(n::T)
    0 < n || error("number to be factored must be positive")
    h = (T=>Int)[]
    n == 1 && return h
    n <= 3 && (h[n] = 1; return h)
    local s::T, p::T
    s = isqrt(n)
    for p in PRIMES
        if p > s
            h[n] = 1
            return h
        end
        if n % p == 0
            while n % p == 0
                h[p] = get(h,p,0)+1
                n = oftype(n, div(n,p))
            end
            n == 1 && return h
            s = isqrt(n)
        end
    end
    p = PRIMES[end]+2
    while p <= s
        if n % p == 0
            while n % p == 0
                h[p] = get(h,p,0)+1
                n = oftype(n, div(n,p))
            end
            if n == 1
                return h
            end
            s = isqrt(n)
        end
        p += 2
    end
    h[n] = 1
    return h
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
struct BackTracking{T1, T2} <: LineSearch
    ratio::T1
	c::T1
	maxiter::T2
	verbose::Bool
end
BackTracking(; ratio=0.5, c=1e-4, maxiter=50, verbose=false) = BackTracking(ratio, c, maxiter, verbose)
function (ls::BackTracking)(f∇f::T, d, x, f_0, ∇f_0, α_0) where T
	ratio, c, maxiter, verbose = ls.ratio, ls.c, ls.maxiter, ls.verbose

    if verbose
        println("Entering line search with step size: ", α_0)
        println("Initial value: ", f_0)
        println("Value at first step: ", f∇f(nothing, x+α_0*d))
    end

	m = dot(d, ∇f_0)
    t = -c*m

    α, β = α_0, α_0

    iter = 0

    f_α = f∇f(nothing, x + α*d) # initial function value

	is_solved = isfinite(f_α) && f_α <= f_0 + c*α*t
    while !is_solved && iter <= maxiter
        iter += 1
        β, α = α, α*ratio # backtrack according to specified ratio
        f_α = f∇f(nothing, x + α*d) # update function value
		is_solved = isfinite(f_α) && f_α <= f_0 + c*α*t
    end

	ls_success = iter >= maxiter ? false : true
	#
    # if verbose
	# 	!ls_success && println("maxiter exceeded in backtracking")
    #     println("Exiting line search with step size: ", α)
    #     println("Exiting line search with value: ", f_α)
    # end
    return α, f_α, ls_success
end
function backtracking!(objective, x, d, f_0, ∇f_0, α_0, opt)
	ratio, c, maxiter, verbose = opt.ratio, opt.c, opt.maxiter, opt.verbose

    # if verbose
    #     println("Entering line search with step size: ", α_0)
    #     println("Initial value: ", f(x))
    #     println("Value at first step: ", f(x+α_0*d))
    # end

	m = dot(d, ∇f_0)
    t = -c*m

    α, β = α_0, α_0

    iter = 0

    f_α = objective(nothing, nothing, x + α*d) # initial function value
	is_solved = false
    while !is_solved && iter <= maxiter
        iter += 1
        β, α = α, α*ratio # backtrack according to specified ratio
        f_α = objective(nothing, nothing, x + α*d) # update function value
		is_solved = isfinite(f_α) && f_0 - f_α >= α*t
    end

	ls_success = iter >= maxiter ? false : true
	#
    # if verbose
	# 	!ls_success && println("maxiter exceeded in backtracking")
    #     println("Exiting line search with step size: ", α)
    #     println("Exiting line search with value: ", f_α)
    # end
    return α, f_α, ls_success
end


text:
"""
Measures Q
"""
function Q(M::Modular)
   m = sum(M.A)
   B = M.A - kron(sum(M.A, 1), sum(M.A, 2)) ./ m
   inner_sum = 0.0
   nr, nc = size(M.A)
   for col = 1:nc
      for row = 1:nr
         if M.S[col,:] == M.S[nc+row,:]
            inner_sum += B[row,col]
         end
      end
   end
   return inner_sum / m
end

"""
Measures Qr'

Qr' is the proportion of arvs established between nodes belonging to the same
modules. Values of Qr' below 0 indicate that there are *more* arcs between than
within modules.

Poisot T. An a posteriori measure of network modularity. *F1000Research* 2013,
2:130, `10.12688/f1000research.2-130.v3`.
"""
function Qr(M::Modular)
   m = sum(M.A)
   within = 0
   nr, nc = size(M.A)
   for col = 1:nc
      for row = 1:nr
         if M.A[row,col] > 0
            if M.S[col,:] == M.S[nc+row,:]
               within += 1
            end
         end
      end
   end
   return 2*(within/m)-1
end


text:
using TMI, Test

@testset "TMI.jl" begin

    TMIversion = "modern_90x45x33_GH10_GH12"
    #TMIversion = "modern_180x90x33_GH10_GH12"
    #TMIversion = "modern_90x45x33_unpub12"

    A, Alu, γ, TMIfile, L, B = config_from_nc(TMIversion)

    @testset "steadyinversion" begin

        yPO₄ = readfield(TMIfile,"PO₄",γ)
        bPO₄ = getsurfaceboundary(yPO₄)
        PO₄pre = steadyinversion(Alu,bPO₄,γ)
        qPO₄ = readfield(TMIfile,"qPO₄",γ)
        b₀ = zerosurfaceboundary(γ)
        PO₄ᴿ = steadyinversion(Alu,b₀,γ,q=qPO₄)
        PO₄total = PO₄ᴿ + PO₄pre
        PO₄direct = steadyinversion(Alu,bPO₄,γ,q=qPO₄)

        ## how big is the maximum difference?
        # could replace with abs function
        @test maximum(PO₄direct - PO₄total) < 0.1
        @test minimum(PO₄direct - PO₄total) > -0.1

        ## oxygen distribution, just be sure it runs
        yO₂ = readfield(TMIfile,"O₂",γ)
        bO₂ = getsurfaceboundary(yO₂)
        O₂ = steadyinversion(Alu,bO₂,γ,q=qPO₄,r=-170.0)

    end

    ############################
    ## trackpathways
    @testset "trackpathways" begin

        @test isapprox(maximum(sum(A,dims=2)),1.0)
        @test minimum(sum(A,dims=2))> -1e-14

        latbox = [50,60]; # 50 N -> 60 N, for example.
        lonbox = [-50, 0]; # 50 W -> prime meridian

        c = trackpathways(Alu,latbox,lonbox,γ)

        @test maximum(c) ≤ 1.0
        @test minimum(c) ≥ 0.0
    end

    @testset "watermassdistribution" begin
        list = ("GLOBAL","ANT","SUBANT",
                "NATL","NPAC","TROP","ARC",
                "MED","ROSS","WED","LAB","GIN",
                "ADEL","SUBANTATL","SUBANTPAC","SUBANTIND",
                "TROPATL","TROPPAC","TROPIND")
        region = list[2]
        g = watermassdistribution(TMIversion,Alu,region,γ)
        @test maximum(g) ≤ 1.0
        @test minimum(g) ≥ 0.0

    end

    #######################################
    ## example: regeneration
    @testset "regeneration" begin

        PO₄ᴿ = regeneratedphosphate(TMIversion,Alu,γ)
        @test maximum(PO₄ᴿ) < 5
        @test minimum(PO₄ᴿ) ≥ 0
    end

    #######################################
    ## example: mean age
    @testset "meanage" begin

        a = meanage(TMIversion,Alu,γ)
        @test maximum(a) < 3000
        @test minimum(a) ≥ 0
    end

    ########################################
    ## example: howoceansfilled
    @testset "howoceansfilled" begin

        # are the areas and volumes consistent?
        v = cellvolume(γ)
        area = cellarea(γ)
        @test sum(0. .< v./area .< 1000.)/length(γ.I) == 1
        volume = volumefilled(TMIversion,Alu,γ)
        # volumefill no smaller than smallest box?
        @test exp10(minimum(volume)) ≥ 5.0
    end

    ####################################
    ## example: surfaceorigin
    @testset "surfaceorigin" begin

        # choose linear or nearest neighbor interpolation
        linearinterp = true #true
        # randomized location.
        loc = wetlocation(γ)

        # choose linear or nearest neighbor interpolation

        if linearinterp
            # temporary tracer to initialize linear interp
            ctmp = zeros(γ.wet)
            δ = interpweights(loc,γ)
        else
            # alternate choice
            δ = nearestneighbormask(loc,γ)
        end
        @test isapprox(sum(filter(!isnan,δ)),1.0)

        origin = surfaceorigin(loc, Alu, γ)

        #@test isapprox(sum(filter(!isnan,origin)),1)
        #@test isapprox(sum(filter(!isnan,origin)),1)
        #@test minimum(origin) ≥ -20
        @test maximum(origin) ≤ 0 # log10(1) = 0
    end

    #########################################
    ## formerly filterdata
    @testset "steadyclimatology" begin

        for ii = 1:2
            y, W⁻, ctrue = synthetic_observations(TMIversion,"θ",γ)
            if ii == 1
                println("NamedTuple type")
                u = (;surface = zerosurfaceboundary(γ))
                b = (;surface = getsurfaceboundary(y))
            else
                println("BoundaryCondition type")
                u = zerosurfaceboundary(γ)
                b = getsurfaceboundary(y)
            end

            uvec = vec(u)
            F,G = costfunction_gridded_obs(uvec,Alu,b,u,y,W⁻,γ)
            fg!(F,G,x) = costfunction_gridded_obs!(F,G,x,Alu,b,u,y,W⁻,γ)
            fg(x) = costfunction_gridded_obs(x,Alu,b,u,y,W⁻,γ)
            f(x) = fg(x)[1]
            J₀,gJ₀ = fg(uvec)

            iterations = 10
            out = steadyclimatology(uvec,fg!,iterations)

            # check with forward differences
            ϵ = 1e-3
            ii = rand(1:sum(γ.wet[:,:,1]))
            println("Location for test =",ii)
            δu = copy(uvec); δu[ii] += ϵ
            ∇f_finite = (f(δu) - f(uvec))/ϵ
            println(∇f_finite)

            fg!(J₀,gJ₀,(uvec+δu)./2) # J̃₀ is not overwritten
            ∇f = gJ₀[ii]
            println(∇f)

            # error less than 10 percent?
            println("Percent error ",100*abs(∇f - ∇f_finite)/abs(∇f + ∇f_finite))
            @test abs(∇f - ∇f_finite)/abs(∇f + ∇f_finite) < 0.1

            # was cost function decreased?
            @test out.minimum < J₀

            # reconstruct by hand to double-check.
            ũ = out.minimizer
            J,gJ = fg(ũ)
            @test J < J₀
        end
    end

    #########################################
    ## globalmap
    @testset "sparsedatamap" begin

        using Statistics, Interpolations

        N = 20
        y, W⁻, ctrue, ytrue, locs, wis = synthetic_observations(TMIversion,"θ",γ,N)

        for ii = 1:2
            if ii == 1
                println("NamedTuple type")
                u = (;surface = zerosurfaceboundary(γ))
                b = (;surface = mean(y) * onesurfaceboundary(γ))
            else
                println("BoundaryCondition type")
                u = zerosurfaceboundary(γ)
                b = mean(y) * onesurfaceboundary(γ)
            end
            uvec = vec(u)
            σb = 5.0
            Q⁻ = 1.0/(σb^2)
            fg(x) = costfunction_point_obs(x,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ)
            f(x) = fg(x)[1]
            J0 = f(uvec)
            J̃₀,gJ₀ = fg(uvec)
            fg!(F,G,x) = costfunction_point_obs!(F,G,x,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ)

            ϵ = 1e-3 # size of finite perturbation
            ii = rand(1:sum(γ.wet[:,:,1]))
            println("gradient check location=",ii)
            δu = copy(uvec); δu[ii] += ϵ
            ∇f_finite = (f(δu) - f(uvec))/ϵ
            println("∇f_finite=",∇f_finite)

            fg!(J̃₀,gJ₀,(uvec+δu)./2) # J̃₀ is not overwritten
            ∇f = gJ₀[ii]
            println("∇f=",∇f)

            # error less than 10 percent?
            println("Percent error=",100*abs(∇f - ∇f_finite)/abs(∇f + ∇f_finite))
            @test abs(∇f - ∇f_finite)/abs(∇f + ∇f_finite) < 0.1
            iterations = 5
            out = sparsedatamap(uvec,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ,iterations)
            # was cost function decreased?
            @test out.minimum < J̃₀
            # reconstruct by hand to double-check.
            ũ = out.minimizer
            J̃,gJ̃ = fg(ũ)
            @test J̃ < J̃₀
        end
    end
end


text:
using DataStructures # needed for PriorityQueue

s = readlines("../data/input15.txt")
N =  length(s)
A = zeros(Int, N, N)
@inbounds for i in 1:N
    for j in 1:N
        A[j, i] = parse(Int, s[j][i])
    end
end

function part1(grid)
    row, col = size(grid)
    dis = fill(typemax(Int), row, col)
    dx = (-1, 0, 1, 0)
    dy = (0, 1, 0, -1)

    st = PriorityQueue{Tuple{Int, Int}, Int}()
    st[(1, 1)] = 0
    dis[1, 1] = grid[1, 1]
    while !isempty(st)
        (kx, ky), _ = peek(st)
        dequeue!(st)

        @inbounds for i in 1:4
            x = kx + dx[i]
            y = ky + dy[i]
            (1 <= x <= col && 1 <= y <= row) || continue

            dcand = dis[ky, kx] + grid[y, x]
            if dis[y, x] > dcand
                dis[y, x] = dcand
                st[(x, y)] = dis[y, x]
            end
        end
    end
    return dis[end, end] - grid[1, 1]
end

function part2(grid)
    row, col = size(grid)
    row2 = 5*row
    col2 = 5*col
    dis = fill(typemax(Int), row2, col2)
    dx = (-1, 0, 1, 0)
    dy = (0, 1, 0, -1)

    st = PriorityQueue{Tuple{Int, Int}, Int}()
    st[(1, 1)] = 0
    dis[1, 1] = grid[1, 1]
    while !isempty(st)
        (kx, ky), _ = peek(st)
        dequeue!(st)

        @inbounds for i in 1:4
            x = kx + dx[i]
            y = ky + dy[i]
            (1 <= x <= col2 && 1 <= y <= row2) || continue

            dcand = mod1(grid[mod1(y, row), mod1(x, col)] + (x-1)÷col + (y-1)÷row, 9)
            dcand += dis[ky, kx]
            if dis[y, x] > dcand
                dis[y, x] = dcand
                st[(x, y)] = dis[y, x]
            end
        end
    end
    return dis[end, end] - grid[1, 1]
end

max_stars_repo_name: lucaferranti/AoC

max_issues_repo_name: lucaferranti/AoC

max_forks_repo_name: lucaferranti/AoC


text:
using Plots

export DefaultVals,Bags,Pens,Turtles

mutable struct DefaultVals
    turn_rad::Real
    step::Real
    do_till_escaped::Bool
    N::Int
    increase_step::Real
    ang::Real
    lsys_data::Tuple
    iterations::Int

    function DefaultVals()
        turn_rad = 0
        step = 100
        do_till_escaped = false
        N = 10
        increase_step = 50
        ang = deg2rad(120)
        lsys_data = (Dict("L" => "L+RF+",
                    "R" => "-FL-R"),
                    "L")
        iterations = 3
        new(turn_rad,step,do_till_escaped,N,increase_step,ang,lsys_data,iterations)
    end
end

mutable struct Bags
    vertexes::Union{Matrix{T} where T<:Real,Tuple}

    function Bags(x::Union{Array{T} where T<:Real,Tuple} = vec([-350. -350. 350. 350.]), y::Union{Array{T} where T<:Real,Tuple} = vec([350. -350. -350. 350.]))
        if x isa Tuple
            x = collect(x)
        end

        if !(x isa Vector)
            # @warn "x is converted from $(typeof(x)) to Array{$(eltype(x)),1} by columns"
            x = vec(x)
        end

        if y isa Tuple
            y = collect(y)
        end

        if !(y isa Vector)
            # @warn "y is converted from $(typeof(y)) to Array{$(eltype(y)),1} by columns"
            y = vec(y)
        end

        if size(x) != size(y)
            DimensionMismatch("vector x,y dimensions mismatch")
        end

        new([x y])
    end
end

mutable struct Pens
    color::Symbol
    size::Int
    up::Bool

    function Pens(;color::Union{String,Symbol,AbstractFloat} = :blue, size = 5, up::Bool = false)
        new(color,size,up)
    end
end

mutable struct Turtles
    bag::Bags
    pen::Pens
    icon_shape::Symbol
    icon_size::Int
    pos::Vector
    heading::AbstractFloat
    plt::Plots.Plot
    history::Array{Any,2}
    junctions::Vector

    function Turtles(;  bag::Bags = Bags(),
                        pen::Pens = Pens(),
                        icon_shape::Symbol = :rect,
                        icon_size::Int = 5,
                        pos::Union{Tuple,Array{T} where T<:Real} = vec([0.0,0.0]),
                        heading::Real = 0.0,
                        world_size::Union{Tuple,Vector{T} where T<:Real} = vec([600., 600.]))

        if pos isa Tuple
            pos = collect(pos)
        end

        if !(pos isa Vector)
            pos = vec(pos)
        end

        if length(pos) != 2
            DimensionMismatch("The new point to goto must have 2 elements [x,y]")
        end

        if world_size isa Tuple
            world_size = vec(world_size)
        end

        if length(world_size) != 2
            DimensionMismatch("The worldSize must have 2 elements [width,hight]")
        end
        plt = plot();
        plot!(plt,  bag.vertexes[:,1],bag.vertexes[:,2],
                    xlims = (-world_size[1],world_size[1]),
                    ylims = (-world_size[2],world_size[2]),
                    leg = false, xgrid = false, ygrid = false,
                    linewidth = pen.size, linecolor = pen.color);
        scatter!(plt,[pos[1]],[pos[2]],marker = (icon_shape,icon_size));
        plot!(plt,[pos[1]],[pos[2]]);
        history = [pos;false]'
        junctions = []
        new(bag,pen,icon_shape,icon_size,pos,heading,plt,history,junctions);
    end
end


text:
#
#      SO(n) - The manifold of the special orthogonal group
#  A point on the manifold is a orthogonal matrix with determinant +1.
#  A point in an tangential space is a skew-symmetric matrix.
#
import LinearAlgebra: Diagonal, norm, dot, nullspace, det, tr, qr, triu, eigvals, diag, svd, sylvester
import Base: exp, log, show, rand, Matrix
export Rotations, SOPoint, SOTVector, getValue
export addNoise, distance, dot, exp, log, manifoldDimension, norm, parallelTransport, randomTVector, randomMPoint, retractionQR, retractionPolar, inverseRetractionPolar, inverseRetractionQR, retraction, inverseRetraction
export zeroTVector
#
# Type definitions
#

@doc doc"""
    Rotations <: Manifold

The manifold $\mathcal M = \mathrm{SO}(n)$ represented by $n\times n$
real-valued orthogonal matrices with determinant $+1$.

# Abbreviation
`SO`

# Constructor
    Rotations(n)

generates the manifold $\mathrm{SO}(n)$ where the integer `n` is the
number of rows or columns of the matrices.
"""
struct Rotations <: Manifold
  name::String
  dimension::Int
  abbreviation::String
  Rotations(n::Int) = new("of the set of rotations in R$n",n, "SO($n)")
end

@doc doc"""
    SOPoint <: MPoint

A point $x$ on the manifold $\mathcal M = \mathrm{SO}(n)$ is represented by an
orthogonal matrix with determinant $+1$ from $\mathbb R^{n\times n}.$

# Constructor
    SOPoint(x)

where `x` is an orthogonal matrix with determinant $+1$ of dimension $n×n$.
"""
struct SOPoint <: MPoint
  value::Matrix
  SOPoint(value::Matrix) = new(value)
end
getValue(x::SOPoint) = x.value;

@doc doc"""
    SOTVector <: TVector

A tangent vector $\xi \in T_x\mathcal M$ on the manifold
$\mathcal M = \mathrm{SO}(n)$. The tangent space is given by as

$T_x\mathrm{SO}(n) = \bigl\{x\xi \in \mathbb R^{n\times n}
\big| \xi + \xi^T = 0 \bigr\}.$

Since the manifold of rotations is a Lie group, it suffices to store just the
skew-symmetric matrix $\xi$. This has to be taken into account in all formulae.

# Constructor
    SOTVector(ξ)

where `ξ` is an $n\times n$ `Matrix` that is skew-symmetric.
"""
struct SOTVector <: TVector
  value::Matrix
  SOTVector(value::Matrix) = new(value)
end
getValue(ξ::SOTVector) = ξ.value;

# Traits
# ---
# (a) SO(n) is a MatrixManifold
@traitimpl IsMatrixM{Rotations}
@traitimpl IsMatrixP{SOPoint}
@traitimpl IsMatrixTV{SOTVector}
# (b) SO(n) is Embedded
@traitimpl IsEmbeddedM{Rotations}
@traitimpl IsEmbeddedP{SOPoint}
@traitimpl IsEmbeddedV{SOTVector}
# (c) SO(n) is a LieGroup
@traitimpl IsLieGroupM{Rotations}
@traitimpl IsLieGroupP{SOPoint}
@traitimpl IsLieGroupV{SOTVector}
#
LieGroupOp(x::SOPoint, y::SOPoint) = SOPoint( transpose(getValue(x))*getValue(y) )

@doc doc"""
    dot(M,x,ξ,ν)

compute the Riemannian inner product for two [`SOTVector`](@ref)s `ξ` and `ν`
from $T_x\mathcal M$ of the [`Rotations`](@ref) manifold `M` given by

$\langle \xi, \nu \rangle_x = \operatorname{tr}(\xi^T\nu)$

i.e. the inner product in the embedded space $\mathbb R^{n\times n}$.
"""
dot(M::Rotations, x::SOPoint, ξ::SOTVector, ν::SOTVector)::Float64 = tr( transpose(getValue(ξ))*getValue(ν) )

@doc doc"""
    distance(M,x,y)

compute the Riemannian distance on [`Rotations`](@ref) manifold `M`
$= \mathrm{SO}(n)$ embedded in $\mathbb R^{n\times n}$, which is given by

$d(x,y) = \lVert \operatorname{log}_{x}y \rVert_x$

where $\operatorname{log}_{\cdot}\cdot$ denotes the logarithmic map on
the [`Rotations`](@ref) $\mathcal M=\mathrm{SO}(n)$.
"""
function distance(M::Rotations,x::SOPoint,y::SOPoint)
  U = transpose(getValue(x)) * getValue(y)
  if (getValue(x) != getValue(y)) && (abs.(U) == one(U))
    return π
  end
  return norm(M, x, log(M,x,y))
end

@doc doc"""
    exp(M,x,ξ,[t=1.0])

compute the exponential map on the [`Rotations`](@ref) manifold `M`$=\mathrm{SO}(n)$ with
respect to the [`SOPoint`](@ref) `x` and the [`SOTVector`](@ref) `ξ`, which can
be shortened with `t` to `tξ`. The formula reads

```math
\operatorname{exp}_{x}(tξ) = x \cdot \operatorname{Exp}(tξ)
```

where $\operatorname{Exp}$ denotes matrix exponential.
"""
exp(M::Rotations,x::SOPoint,ξ::SOTVector,t::Float64=1.0) = SOPoint(getValue(x) * exp( t * getValue(ξ)))

@doc doc"""
    log(M,x,y)

compute the logarithmic map on the [`Rotations`](@ref) manifold
`M`$=\mathrm{SO}(n)$, which is given by

```math
\operatorname{log}_{x} y =
  \frac{1}{2} \bigl(\operatorname{Log}(x^{\mathrm{T}}y)
  - (\operatorname{Log} x^{\mathrm{T}}y)^{\mathrm{T}}),
```

where $\operatorname{Log}$ denotes the matrix logarithm.
"""
function log(M::Rotations,x::SOPoint,y::SOPoint)
  U = transpose(getValue(x)) * getValue(y)
  if (getValue(x) != getValue(y)) && (abs.(U) == one(U))
    throw( ErrorException("The points $x and $y are antipodal, thus these input parameters are invalid.") )
  end
  U1 = real(log(U))
  U2 = 0.5 * (U1 - transpose(U1))
  return SOTVector(U2)
end

@doc doc"""
    manifoldDimension(x)

return the dimension of the [`Rotations`](@ref) manifold `M`$= \mathrm{SO}(n)$, the
[`SOPoint`](@ref) `x`, itself embedded in $\mathbb R^{n\times n}$, belongs to.
The dimension is defined by

$\frac{n(n-1)}{2}.$
"""
manifoldDimension(x::SOPoint)::Integer = size(getValue(x), 1)*(size(getValue(x), 2)-1)/2

@doc doc"""
    manifoldDimension(M)

return the dimension of the [`Rotations`](@ref) manifold `M`$= \mathrm{SO}(n)$.
The dimension is defined by

$\frac{n(n-1)}{2}.$
"""
manifoldDimension(M::Rotations)::Integer = (M.dimension*(M.dimension-1))/2

@doc doc"""
    norm(M,x,ξ)

compute the norm of the [`SOTVector`](@ref) `ξ` in the tangent space
$T_x\mathcal M$ at [`SOPoint`](@ref) `x` of the [`Rotations`](@ref) manifold `M`.

$\lVert \xi \rVert_x = \sqrt{\sum_{i,j=0}^n \xi_{ij}^2}$

where $\xi_{ij}$ are the entries of the skew-symmetric matrix `ξ`, i.e. the norm
is given by the Frobenius norm of `ξ`.
"""
norm(M::Rotations, x::SOPoint, ξ::SOTVector) = norm(getValue(ξ))

@doc doc"""
    parallelTransport(M,x,y,ξ)

compute the parallel transport of the [`SOTVector`](@ref) `ξ` from
the tangent space $T_x\mathcal M$ at [`SOPoint`](@ref) `x` to
$T_y\mathcal M$ at [`SOPoint`](@ref) `y` on the [`Rotations`](@ref) `M` along
$g$ provided that the corresponding [`geodesic`](@ref) $g(\cdot;x,y)$ is unique.
Since we have only stored the skew-symmetric matrix as a
[`SOTVector`](@ref) `ξ`, the function returns the the [`SOTVector`](@ref) `ξ`.
"""
parallelTransport(M::Rotations,x::SOPoint,y::SOPoint,ξ::SOTVector) = ξ

@doc doc"""
    randomTVector(M,x[, type=:Gaussian, σ=1.0])

return a random [`SOTVector`](@ref) in the tangent space
$T_x\mathrm{SO}(n)$ of the [`SOPoint`](@ref) `x` on the [`Rotations`](@ref)
manifold `M` by generating a random skew-symmetric matrix. The function
takes the real upper triangular matrix of a (Gaussian) random matrix $A$ with
dimension $n\times n$ and subtracts its transposed matrix.
Finally, the matrix is normalized.
"""
function randomTVector(M::Rotations,x::SOPoint, ::Val{:Gaussian}, σ::Real=1.0)
  if M.dimension==1
    SOTVector(zeros(1,1))
  else
    A = randn(Float64, M.dimension, M.dimension)
    A = triu(A,1) - transpose(triu(A,1))
    A = (1/norm(A))*A
    SOTVector(A)
  end
end

@doc doc"""
    randomMPoint(M[, type=:Gaussian, σ=1.0])

return a random [`SOPoint`](@ref) `x` on the manifold [`Rotations`](@ref) `M`
by generating a (Gaussian) random orthogonal matrix with determinant $+1$. Let

$QR = A$

be the QR decomposition of a random matrix $A$, then the formula reads

$x = QD$

where $D$ is a diagonal matrix with the signs of the diagonal entries of $R$,
i.e.

$D_{ij}=\begin{cases} \operatorname{sgn}(R_{ij}) & \text{if} \; i=j \\ 0 & \, \text{otherwise} \end{cases}.$

It can happen that the matrix gets -1 as a determinant. In this case, the first
and second columns are swapped.
"""
function randomMPoint(M::Rotations, ::Val{:Gaussian}, σ::Real=1.0)
  if M.dimension==1
    SOPoint(ones(1,1))
  else
    A=randn(Float64, M.dimension, M.dimension)
    s=diag(sign.(qr(A).R))
    D=Diagonal(s)
    C = qr(A).Q*D
    if det(C)<0
      C[:,[1,2]] = C[:,[2,1]]
    end
    SOPoint(C)
  end
end

@doc doc"""
    retractionPolar(M,x,ξ [,t=1.0])

move the [`SOPoint`](@ref) `x` in the direction of the [`SOTVector`](@ref) `ξ`
on the  [`Rotations`](@ref) manifold `M`. This SVD-based retraction is a second-order
approximation of the exponential map. Let

$USV = x + txξ$

be the singular value decomposition, then the formula reads

$\operatorname{retr}_x\xi = UV^\mathrm{T}$
"""
function retractionPolar(M::Rotations, x::SOPoint, ξ::SOTVector, t::Float64=1.0)
  y = getValue(x) + t * getValue(x) * getValue(ξ)
  S = svd(y)
  A = S.U * transpose(S.V)
  SOPoint(A)
end

@doc doc"""
    retractionQR(M,x,ξ [,t=1.0])

move the [`SOPoint`](@ref) `x` in the direction of the [`SOTVector`](@ref) `ξ`
on the [`Rotations`](@ref) manifold `M`. This QR-based retraction is a
first-order approximation of the exponential map. Let

$QR = x + txξ$

be the QR decomposition, then the formula reads

$\operatorname{retr}_x\xi = QD$

where the matrix $D$ is given by

```math
D_{ij}=\begin{cases}
\operatorname{sgn}(R_{ij}+0,5) & \text{if} \; i=j \\
0 & \, \text{otherwise.}
\end{cases}
```
"""
function retractionQR(M::Rotations, x::SOPoint, ξ::SOTVector, t::Float64=1.0)
  y = getValue(x) + t * getValue(x)*getValue(ξ)
  QRdecomp = qr(y)
  d = diag(QRdecomp.R)
  D = Diagonal( sign.( sign.(d .+ 0.5) ) )
  A = QRdecomp.Q * D
  SOPoint(A)
end
retraction(M::Rotations, x::SOPoint, ξ::SOTVector, t::Float64=1.0) = retractionQR(M,x,ξ,t)

@doc doc"""
    inverseRetractionPolar(M,x,y)

return a [`SOTVector`](@ref) `ξ` of the tagent space $T_x\mathrm{SO}(n)$
of the [`SOPoint`](@ref) `x` on the [`Rotations`](@ref) manifold `M`
with which the [`SOPoint`](@ref) `y` can be reached by the
[`retractionPolar`](@ref) after time 1. The formula reads

$ξ = -\frac{1}{2}(x^{\mathrm{T}}ys - (x^{\mathrm{T}}ys)^{\mathrm{T}})$

where $s$ is the solution to the Sylvester equation

$x^{\mathrm{T}}ys + s(x^{\mathrm{T}}y)^{\mathrm{T}} + 2\mathrm{I}_n = 0.$
"""
function inverseRetractionPolar(M::Rotations, x::SOPoint, y::SOPoint)
  A = transpose(getValue(x)) * getValue(y)
  H = 2 * one(getValue(x))
  K = convert(Array{Float64,2}, transpose(A))
  B = sylvester(A, K, H)
  C =  A * B
  SOTVector( -0.5 * ( C - transpose(C) ) )
end

@doc doc"""
    inverseRetractionQR(M,x,y)

return a [`SOTVector`](@ref) `ξ` of the tagent space $T_x\mathrm{SO}(n)$
of the [`SOPoint`](@ref) `x` on the [`Rotations`](@ref) manifold `M`
with which the [`SOPoint`](@ref) `y` can be reached by the
[`retractionQR`](@ref) from the [`SOPoint`](@ref) `x` after time 1.
"""
function inverseRetractionQR(M::Rotations, x::SOPoint, y::SOPoint)
  A = transpose(getValue(x)) * getValue(y)
  R = zeros(M.dimension, M.dimension)
  for i = 1:M.dimension
    b = zeros(i)
    b[end] = 1
    b[1:(end-1)] = - transpose(R[1:(i-1), 1:(i-1)]) * A[i, 1:(i-1)]
    R[1:i, i] = A[1:i, 1:i] \ b
  end
  C =  A * R
  SOTVector( 0.5 * ( C - transpose(C) ) )
end
inverseRetraction(M::Rotations, x::SOPoint, y::SOPoint) = inverseRetractionQR(M,x,y)

@doc doc"""
    zeroTVector(M,x)

return a zero [`SOTVector`](@ref) $\xi$ from the tagent space $T_x\mathrm{SO}(n)$
of [`SOPoint`](@ref) `x` on the [`Rotations`](@ref) manifold `M`, i.e. a zero
matrix.
"""
zeroTVector(M::Rotations, x::SOPoint) = SOTVector( zero(getValue(x)) )

function validateMPoint(M::Rotations, x::SOPoint)
  if det(getValue(x)) ≉ 1
    throw( ErrorException("The determinant of $x has to be +1 but it is $(det(getValue(x))).") )
  end
  if transpose(getValue(x))*getValue(x) ≉ one(getValue(x))
    throw( ErrorException("$x has to be orthogonal but it's not.") )
  end
  return true
end

function validateTVector(M::Rotations, x::SOPoint, ξ::SOTVector)
  if transpose(getValue(ξ))+getValue(ξ) ≉ zero(getValue(ξ))
    throw( ErrorException("$ξ has to be skew-symmetric but it's not.") )
  end
  return true
end


text:
module FinancialMath

export fv, monthly_payment, monthly_rate

include("core.jl")

end # module


text:
"""
    generate(filename; ofolder="", failure=false)

Generates the julia tests from the file filename. The tests in filename must be written
using the ITL domain specific language. The tests are written in a .jl file with the same
name of filename. The folder where to save the output file is specifified

If failure=true, than each test is also executed before printing to the target file. If the
test fails, then the test is generated as ´@test_broken`.
"""
function generate(filename; failure=true, test_warn=true, output="test_ITF1788")

    # read file
    src = joinpath(@__DIR__, "itl", filename)
    f = open(src)
    lines = readlines(f)
    close(f)

    # file where to Write
    dest = joinpath(output, filename[1:end-4]*".jl")
    mkpath(output)
    f = open(dest; write=true)

    # where testcase blocks start
    rx_start = r"^\s*testcase\s+\S+\s+\{\s*$"
    rx_end = r"^\s*\}\s*$"
    block_start = findall(x -> occursin(rx_start, x), lines)
    block_end = findall(x -> occursin(rx_end, x), lines)

    # check opening and closing blocks match
    length(block_start) == length(block_end) || throw(ArgumentError("opening and closing braces not not much in $filename"))

    for (bstart, bend) in zip(block_start, block_end)
        testset = parse_block(lines[bstart:bend]; failure=failure, test_warn=test_warn)
        write(f, testset)
    end

    close(f)
    nothing
end

function generate(; failure=true, test_warn=true, output="test_ITF1788")

    files = ("atan2.itl",
            "c-xsc.itl",
            "fi_lib.itl",
            "ieee1788-constructors.itl",
            "ieee1788-exceptions.itl",
            "libieeep1788_bool.itl",
            "libieeep1788_cancel.itl",
            "libieeep1788_class.itl",
            "libieeep1788_elem.itl",
            "libieeep1788_num.itl",
            "libieeep1788_overlap.itl",
            "libieeep1788_rec_bool.itl",
            "libieeep1788_reduction.itl",
            "libieeep1788_set.itl",
            "mpfi.itl",
            "abs_rev.itl",
            "libieeep1788_rev.itl",
            "libieeep1788_mul_rev.itl",
            "pow_rev.itl"
    )
    for file in files
        generate(file; failure=failure, test_warn=test_warn, output=output)
    end

    f = open(joinpath(output, "run_ITF1788.jl"); write=true)
    for file in files
        write(f, """include("$(file[1:end-4]).jl")\n""")
    end
    close(f)
end


text:
using StackOverflow

"""
    searchtag(intitle::String)
Adding the ability to search for an tag through StackOverflow.jl while in the Julia Terminal.
"""
function searchtag(intitle::String)
    data = replace(intitle, " " => "%20")
    r = HTTP.request("GET", "https://api.stackexchange.com/2.2/tags?order=desc&sort=popular&inname=$(data)&site=stackoverflow")
    json_obj = StackOverflow.convert_HTTP_Response_To_JSON(r)
    for (k, v) in json_obj
        if occursin("items", k)
            for item in v
                for (key, value) in item
                    if key == "name" || key == "count"
                        println("$(key) $(value)")
                    end
                end
                println("\n")
            end
        end
    end
end


text:
module TestPyCallBuild

include(joinpath(dirname(@__FILE__), "..", "deps", "depsutils.jl"))
include(joinpath(dirname(@__FILE__), "..", "deps", "buildutils.jl"))

using Test

@testset "find_libpython" begin
    for python in ["python", "python2", "python3"]
        # TODO: In Windows, word size should also be checked.
        Sys.iswindows() && break
        if Sys.which(python) === nothing
            @info "$python not available; skipping test"
        else
            @test isfile(find_libpython(python)[2])
        end
    end

    # Test the case `find_libpython.py` does not print anything.  We
    # use the command `true` to mimic this case.
    if Sys.which("true") === nothing
        @info "no `true` command; skipping test"
    else
        let err, msg
            @test try
                find_libpython("true")
                false
            catch err
                err isa ErrorException
            end
            msg = sprint(showerror, err)
            @test occursin("Couldn't find libpython", msg)
            @test occursin("ENV[\"PYCALL_DEBUG_BUILD\"] = \"yes\"", msg)
        end
    end

    # Test the case `dlopen` failed to open the library.
    let err, msg
        @test try
            find_libpython("python"; _dlopen = (_...) -> error("dummy"))
            false
        catch err
            err isa ErrorException
        end
        msg = sprint(showerror, err)
        @test occursin("Couldn't find libpython", msg)
        @test occursin("ENV[\"PYCALL_DEBUG_BUILD\"] = \"yes\"", msg)
    end
end

end  # module


text:
include("utils.jl")

function Ctx(env::AlgorithmicEnv, mode::Symbol=:human)
    if mode == :human
        HumanCtx()
    elseif mode == :ansi
        ANSICtx()
    elseif mode == :no_render
        NoCtx()
    else
        error("Unrecognized mode in Ctx(): $(mode)")
    end
end

function colorize_bold(string::Union{AbstractString, Char}, color::Symbol)
    colormap = Dict(
    :gray => 30, # highlights to black
    :red => 31,
    :green => 32,
    :yellow => 34,
    :magenta => 35,
    :cyan => 36,
    :white => 37)

    num = colormap[color]
    num += 10 # For highlighting
    return "\x1b[$num;1m" * string * "\x1b[0m"
end

_get_str_target(env::AlgorithmicEnv, pos) =
    pos < 1 || length(env.target) < pos ? " " : env.charmap[env.target[pos]]

function render_observation(env::TapeAlgorithmicEnv)
    x = env.read_head_position
    x_str = "Observation Tape    : "
    for i=-1:length(env.input_data)+2
        x_str *= i == x ? colorize_bold(env.charmap[_get_obs(env, i)[1]], :green) :
                          env.charmap[_get_obs(env, i)[1]]
    end
    x_str *= "\n"
    return x_str
end

function render_observation(env::GridAlgorithmicEnv)
    x = env.read_head_position
    label = "Observation Grid    : "
    x_str = ""
    for j=0:env.rows
        if j != 0
            for _=1:length(label)
                x_str *= " "
            end
        end
        for i=0:length(env.input_data)
            println("Okay for i=$i and j=$j, _get_obs=$(_get_obs(env, [i, j]))")
            x_str *= i == x[1] && j == x[2] ? colorize_bold(env.charmap[_get_obs(env, [i, j])[1]], :green) :
                                              env.charmap[_get_obs(env, [i, j])[1]]
        end
        x_str *= "\n"
    end
    x_str = label * x_str
    return x_str
end

function drawcanvas!(env::AlgorithmicEnv)
    inp = "Total length of input instance: $(length(env.input_data)), step: $(env.time)\n"
    output = ""
    output *= inp
    x = env.read_head_position
    y = env.write_head_position
    action = env.last_action
    !isnothing(action) &&
        ((inp_act, out_act, pred) = action)
    equal_sign = ""
    for i=1:length(inp)+1
        equal_sign *= '='
    end
    output *= equal_sign * '\n'
    y_str =      "Output Tape         : "
    target_str = "Targets             : "
    !isnothing(action) &&
        (pred_str = env.charmap[pred])
    x_str = render_observation(env)
    for i=-1:length(env.target)+2
        target_str *= _get_str_target(env, i)
        if i < y - 1
            y_str *= _get_str_target(env, i)
        elseif i == (y - 1)
            if !isnothing(action) && out_act == 2
                color = pred == env.target[i] ? :green : :red
                y_str *= colorize_bold(pred_str, color)
            else
                y_str *= _get_str_target(env, i)
            end
        end
    end
    output *= x_str * y_str * '\n' * target_str * "\n\n"

    if !isnothing(action)
        output *= "Current reward      :   $(env.last_reward)\n"
        output *= "Cumulative reward   :   $(env.episode_total_reward)\n"
        output *= "Action              :   Tuple(move over input: $(String(env.MOVEMENTS[inp_act])),\n"
        output *= "                              write to the output tape: $(out_act==2)\n"
        output *= "                              prediction: $pred_str)\n"
    else
        output *= "\n\n\n\n\n"
    end
    return output
end


text:
using NumericalDifferentiation, PyPlot

x = range(-1, stop=1, length=101)
y = abs.(x) + (rand(length(x)) .- .5)/5

u_TV = differentiate(x, y, TotalVariation(), 0.2, 1e-3, maxit=200)
u_TK = differentiate(x, y, Tikhonov(), 0.0002)

figure()
plot(x, y .- y[1],
     x, integrationoperator(x,u_TV),
     x, integrationoperator(x,u_TK))
plot(x, u_TV,
     x, u_TK)

z = sin.(π*x) + (rand(length(x)) .- .5)/5

u_TV = differentiate(x, z, TotalVariation(), 0.2, 1e-3, maxit=200)
u_TK = differentiate(x, z, Tikhonov(), 0.0002)

figure()
plot(x, z .- z[1],
     x, integrationoperator(x,u_TV),
     x, integrationoperator(x,u_TK))

plot(x, u_TV,
     x, u_TK)


text:
# Test of accumulators

ct = counter(Compat.ASCIIString)
@assert isa(ct, Accumulator{Compat.ASCIIString,Int})

@test ct["abc"] == 0
@test !haskey(ct, "abc")
@test isempty(collect(keys(ct)))

push!(ct, "a")
@test haskey(ct, "a")
@test ct["a"] == 1

push!(ct, "b", 2)
@test haskey(ct, "b")
@test ct["b"] == 2

# Test convert
push!(ct, "b", 0x3)
@test ct["b"] == 5

@test !haskey(ct, "abc")
@test ct["abc"] == 0

@test length(ct) == 2
@test length(collect(ct)) == 2
@test length(collect(keys(ct))) == 2
@test length(collect(values(ct))) == 2
@test sum(ct) == 6

ct2 = counter(["a", "a", "b", "b", "a", "c", "c"])
@test isa(ct2, Accumulator{Compat.ASCIIString,Int})
@test haskey(ct2, "a")
@test haskey(ct2, "b")
@test haskey(ct2, "c")
@test ct2["a"] == 3
@test ct2["b"] == 2
@test ct2["c"] == 2

push!(ct, ct2)
@test ct["a"] == 4
@test ct["b"] == 7
@test ct["c"] == 2

ct3 = counter(Dict([("a",10), ("b",20)]))
@test isa(ct3, Accumulator{Compat.ASCIIString,Int})
@test haskey(ct3, "a")
@test haskey(ct3, "b")
@test ct3["a"] == 10
@test ct3["b"] == 20

ctm = merge(ct2, ct3)
@test isa(ctm, Accumulator{Compat.ASCIIString,Int})
@test haskey(ctm, "a")
@test haskey(ctm, "b")
@test haskey(ctm, "c")
@test ctm["a"] == 13
@test ctm["b"] == 22
@test ctm["c"] == 2

@test pop!(ctm, "b") == 22
@test !haskey(ctm, "b")
@test ctm["b"] == 0


text:
#
# Date created: 2022-03-14
# Author: aradclif
#
#
############################################################################################
function sizeblock(N::Int)
    block = Expr(:block)
    for d = 1:N
        ex = Expr(:(=), Symbol(:D_, d), Expr(:call, :size, :A, d))
        push!(block.args, ex)
    end
    block
end

function sizeproductsblock(N::Int)
    block = Expr(:block)
    for k = 3:N
        ex = Expr(:(=), Symbol(:D_, ntuple(identity, k - 1)...),
                  Expr(:call, :*, ntuple(d -> Symbol(:D_, d), k - 1)...))
        push!(block.args, ex)
    end
    block
end

function sumprodprecomputed2(N::Int)
    e = Expr(:call, :+, Symbol(:i_, 1))
    for k = 2:N
        if k == 2
            ex = Expr(:call, :*, Symbol(:D_, 1), Symbol(:i_, 2))
            push!(e.args, ex)
        else
            ex = Expr(:call, :*, Symbol(:D_, ntuple(identity, k - 1)...), Symbol(:i_, k))
            push!(e.args, ex)
        end
    end
    e
end

function sumprodconstant(N::Int)
    Expr(:(=), :D_sp,
         Expr(:call, :+, ntuple(d -> Expr(:call, :*, ntuple(i -> Symbol(:D_, i), d)...), N - 1)...))
end

function outerloopgen(N::Int, D)
    block = Expr(:block)
    params = D.parameters
    for d = N:-1:1
        if params[d] !== Static.One
            ex = Expr(:(=), Symbol(:i_, d), Expr(:call, :axes, :A, d))
            push!(block.args, ex)
        end
    end
    Expr(:for, block)
end
function innerloopgen(N::Int, D)
    block = Expr(:block)
    params = D.parameters
    for d = N:-1:1
        if params[d] === Static.One
            ex = Expr(:(=), Symbol(:i_, d), Expr(:call, :axes, :A, d))
            push!(block.args, ex)
        end
    end
    Expr(:for, block)
end

function innerpost(N::Int, D)
    params = D.parameters
    b = Expr(:ref, :B, ntuple(d -> params[d] === Static.One ? 1 : Symbol(:i_, d), N)...)
    c = Expr(:ref, :C, ntuple(d -> params[d] === Static.One ? 1 : Symbol(:i_, d), N)...)
    Expr(:block, Expr(:(=), b, :m), Expr(:(=), c, :j))
end

function maxcompareblock(N::Int)
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    d = sumprodprecomputed2(N)
    push!(d.args, :D_sp)
    yₑ = Expr(:(=), :y, Expr(:call, :(>), a, :m)) # f should only be > or <
    # mₑ = Expr(:(=), :m, Expr(:if, :y, a, :m))
    # jₑ = Expr(:(=), :j, Expr(:if, :y, d, :j))
    mₑ = Expr(:(=), :m, Expr(:call, :ifelse, :y, a, :m))
    jₑ = Expr(:(=), :j, Expr(:call, :ifelse, :y, d, :j))
    Expr(:block, yₑ, mₑ, jₑ)
end

function lvfindmax_quote(N::Int, D)
    b1 = sizeblock(N)
    b2 = sizeproductsblock(N)
    b3 = Expr(:block, sumprodconstant(N), Expr(:(=), :D_sp, Expr(:call, :-, :D_sp)))
    outer = outerloopgen(N, D)
    b4 = Expr(:block, Expr(:(=), :j, 1), Expr(:(=), :m, Expr(:call, :typemin, :T)))
    inner = innerloopgen(N, D)
    b5 = maxcompareblock(N)
    push!(inner.args, b5)
    push!(b4.args, inner)
    b6 = innerpost(N, D)
    push!(b4.args, b6.args...)
    push!(outer.args, b4)
    return quote
        $b1
        $b2
        $b3
        @turbo $outer
    end
end
@generated function _lvfindmax!(C::AbstractArray{Tₒ, N}, A::AbstractArray{T, N},
                             B::AbstractArray{T, N}, dims::D) where {Tₒ, T, N, D}
    lvfindmax_quote(N, D)
end
function lvfindmax(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    if 1 ∈ dims
        return bfindmax(A, dims)
    else
        Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : size(A, d), N)
        B = similar(A, Dᴮ′)
        C = similar(B, Int)
        _lvfindmax!(C, A, B, Dᴮ′)
        return B, CartesianIndices(A)[C]
    end
end

A = rand(4, 3, 5);
N = ndims(A)
dims = (2,)
DD = typeof(ntuple(d -> d ∈ dims ? StaticInt(1) : size(A, d), N))
lvfindmax_quote(N, DD)
for d₂ = 1:ndims(A), d₁ = 1:ndims(A)
    dims = (d₁, d₂)
    @assert lvfindmax(A, dims) == findmax(A, dims=dims)
end
@benchmark findmax(A, dims=dims)
@benchmark lvfindmax(A, dims)