text:
## CPItree.jl - Tipos y métodos para operar la estructura jerárquica del IPC

##  ----------------------------------------------------------------------------
#   Main type definitions
#   ----------------------------------------------------------------------------

"""
    Item{T}(code::String, name::String, weight::T<:AbstractFloat)

Representa un gasto básico en la estructura de nodos del IPC. Es el nivel más
bajo de la estructura. Almacena el código del gasto básico, su nombre o
descripción y su ponderación en el IPC. Los datos de este nodo deben estar
disponibles en algún [`FullCPIBase`](@ref).

Posee los campos:
- `code`: que almacena el código del gasto básico como un `String`.
- `name`: que almacena el nombre del gasto básico como un `String`.
- `weight::T`: que almacena la ponderación del gasto básico como un valor
  flotante de tipo `T`.

Aunque este nodo será usualmente creado de manera automática por métodos como
[`CPITree`](@ref), se pueden crear estructuras jerárquicas manualmente. Por
ejemplo, para crear un nodo del nivel inferior:
```julia-repl
julia> Item("_011101", "Item A", 7.352945f0)
Item{Float32}("_011101", "Item A", 7.352945f0)
```

Ver también: [`Group`](@ref), [`CPITree`](@ref).
"""
struct Item{T}
    code::String
    name::String
    weight::T
end

"""
    Group{S,T}

    Group(code, name, children::Vector{S}) where S
    Group(code, name, children...)

Representa un nodo de agrupación de cualquier nivel en la estructura de nodos
del IPC. Puede almacenar gastos básicos u otros grupos de mayor jerarquía.
Almacena el código del grupo, su nombre o descripción y su ponderación en el
IPC.

Posee los campos:
- `code`: almacena el código del grupo como un `String`.
- `name`: almacena el nombre del grupo como un `String`.
- `weight::T`: almacena la ponderación del grupo como un valor flotante de tipo
  `T`.
- `children::Vector{S}`: almacena el vector de nodos "hijos" de la estructura.
  Por ejemplo, este vector podría ser un vector de elementos `Item` para agrupar
  un conjunto de gastos básicos.

Aunque este nodo será usualmente creado de manera automática por métodos como
[`CPITree`](@ref), se pueden crear estructuras jerárquicas manualmente. Por
ejemplo, para crear un grupo:
```julia-repl
julia> a = Item("_011201", "Item B", 6.7442417f0)
Item{Float32}("_011201", "Item B", 6.7442417f0)

julia> b = Item("_011202", "Item C", 7.394718f0)
Item{Float32}("_011202", "Item C", 7.394718f0)

julia> g = Group("_0112", "Subgr._0112", a, b)
Group{Item{Float32}, Float32}
_0112: Subgr._0112 [14.13896]
├─ _011201: Item B [6.7442417]
└─ _011202: Item C [7.394718]
```

Ver también: [`Item`](@ref), [`CPITree`](@ref).
"""
struct Group{S,T}
    code::String
    name::String
    weight::T
    children::Vector{S}

    function Group(code, name, children...)
        sum_weights = sum(child.weight for child in children)
        T = eltype(sum_weights)
        S = eltype(children)
        new{S,T}(code, name, sum_weights, S[children...])
    end
    Group(code, name, children::Vector{S}) where {S} = Group(code, name, children...)
end

# Redefine methods for getting children
children(::Item) = ()
children(g::Group) = g.children

# Redefine how to print a node in the print_tree function
printnode(io::IO, g::Item) = print(io, g.code * ": " * g.name * " [" * string(g.weight) * "] ")
printnode(io::IO, g::Group) = print(io, g.code * ": " * g.name * " [" * string(g.weight) * "] ")

# How to show a Group
function Base.show(io::IO, g::Group)
    println(io, typeof(g))
    print_tree(io, g)
end

##  ----------------------------------------------------------------------------
#   CPI tree functions
#   ----------------------------------------------------------------------------

"""
    cpi_tree_nodes(codes::Vector{<:AbstractString};
        characters::(NTuple{N, Int} where N), depth::Int=1, chars::Int=characters[depth], prefix::AbstractString="",
        base::FullCPIBase,
        group_names::Vector{<:AbstractString},
        group_codes::Vector{<:AbstractString})

Construye y devuelve una lista de nodos a partir de la lista de códigos `codes`
o de la especificación jerárquica de caracteres en `characters`. Los nombres y
las ponderaciones del nivel inferior (nivel de gasto básico) son
obtenidas de la estructura `base`. Se debe proveer el vector de códigos y
nombres de todas las jerarquías superiores en la estructura de códigos en los
vectores `group_names` y `group_codes`.

Esta función permite crear únicamente la estructura jerárquica de nodos. Para construir de manera automática una estructura del IPC, se recomienda utilizar preferentemente [`CPITree`](@ref).

Vea también: [`CPITree`](@ref).
"""
function cpi_tree_nodes(codes::Vector{<:AbstractString};
    characters::(NTuple{N, Int} where N), depth::Int=1, chars::Int=characters[depth], prefix::AbstractString="",
    base::FullCPIBase,
    group_names::Vector{<:AbstractString},
    group_codes::Vector{<:AbstractString})

    # Get available starting codes
    available = filter(code -> startswith(code, prefix), codes)

    # Base case:
    # If code length is the last available then we reached a leaf node
    if chars == last(characters)
        # With available codes construct leaf CPI nodes
        children = map(available) do code
            # Find the code index
            icode = findfirst(==(code), codes)
            Item(code, base.names[icode], base.w[icode])
        end
        return children
    end

    # Get possible prefixes values from the available list. Possible prefixes
    # are the ones from the beginning of the string to the number of chars,
    # which depends on the depth we are on.
    possibles = unique(getindex.(available, Ref(1:chars)))

    # For each available code, call the function itself with the next downward hierarchy
    groups = map(possibles) do prefixcode
        # Go get the children in the next downward level
        children = cpi_tree_nodes(codes;
            characters, depth=depth+1, prefix=prefixcode,
            base, group_names, group_codes
        )
        # Get the group name
        @debug "Prefix code:" prefixcode
        i = findfirst(==(prefixcode), group_codes)
        if i === nothing
            @warn "Código de grupo para $(prefixcode) no encontrado en `group_codes`. Utilizando nombre genérico."
            gname = "Group: $prefixcode"
        else
            gname = group_names[i]
        end
        # With the children create a group
        group = Group(prefixcode, gname, children)
        group
    end

    # Return the group for the upward parents
    return groups
end

"""
    get_cpi_tree(;
        base::FullCPIBase,
        group_names::Vector{<:AbstractString},
        group_codes::Vector{<:AbstractString},
        characters::(NTuple{N, Int} where N),
        upperlevel_code = "_0",
        upperlevel_name = "IPC")

Función superior para obtener estructura jerárquica del IPC. Devuelve el nodo
superior del árbol jerárquico. Utiliza la función de más bajo nivel
`cpi_tree_nodes` para construir los nodos del nivel más alto y hacia abajo en la
estructura jerárquica. Se debe proveer el vector de códigos y nombres de todas
las jerarquías superiores en la estructura de códigos en los vectores
`group_names` y `group_codes`.

Esta función permite crear únicamente la estructura jerárquica de nodos. Para construir de manera automática una estructura del IPC, se recomienda utilizar preferentemente [`CPITree`](@ref).

Vea también: [`CPITree`](@ref).
"""
function get_cpi_tree(;
    base::FullCPIBase,
    group_names::Vector{<:AbstractString},
    group_codes::Vector{<:AbstractString},
    characters::(NTuple{N, Int} where N),
    upperlevel_code = "_0",
    upperlevel_name = "IPC")

    length(characters) >= 2 || throw(ArgumentError("`characters` debe ser una tupla con al menos dos valores"))
    for i in 2:length(characters)
        characters[i] > characters[i-1] || throw(ArgumentError("Valores en `characters` deben ser ascendentes hasta el largo de los códigos en base"))
    end
    # Get the codes list from the FullCPIBase object
    codes = base.codes
    @debug "Codes:" codes

    # Call lower level tree building function
    upper_nodes = cpi_tree_nodes(codes; characters, base, group_names, group_codes)

    # Build upper level tree node
    tree = Group(upperlevel_code, upperlevel_name, upper_nodes)
    tree
end

## Build functions to find the inner tree of a given code
function find_tree(code, tree::Item)
    code == tree.code && return tree
    nothing
end

function find_tree(code, tree::Group)
    # Most basic case: the code is the same as the tree in which to search
    code == tree.code && return tree

    # Search in the tree's nodes
    for child in tree.children
        # If code searched is one of the children, return the child
        code == child.code && return child

        # Look if the code starts the same as the child's code, i.e the code
        # contains the child's code. If so, find in the inner tree and break out
        # of this level's search
        contains(code, child.code) && return find_tree(code, child)
    end

    # Code not found at any level
    nothing
end

# Redefine getindex to search for specific nodes within the tree
Base.getindex(tree::Union{Item, Group}, code::AbstractString) = find_tree(code, tree)

##  ----------------------------------------------------------------------------
#   Functions to compute any code's price index from the lower level data
#   ----------------------------------------------------------------------------

# Basic case: compute index of an Item, which is stored in the `base` structure
function compute_index(good::Item, base::FullCPIBase)
    i = findfirst(==(good.code), base.codes)
    # Maybe add a warning here if code not found
    base.ipc[:, i]
end

# Recursive function to compute inner price indexes for groups
function compute_index(group, base::FullCPIBase)
    # Get the indexes of the children. At the lowest level the dispatch will
    # select compute_index(::Item, ::FullCPIBase) to return the Goods indices
    ipcs = mapreduce(c -> compute_index(c, base), hcat, group.children)

    # If there exists only one good in the group, that is the group's index
    size(ipcs, 2) == 1 && return ipcs

    # Get the weights
    weights = map(c -> c.weight, group.children)

    # Return normalized sum product
    ipcs * weights / sum(weights)
end

# Edge case called when searching for a node that doest not exist. For example,
# if called compute_index(tree["_0101101"], base) and node with code "_0101101"
# does not exist, returns nothing and raises a warning
function compute_index(::Nothing, ::FullCPIBase)
    @warn "Nodo no disponible en la estructura"
    nothing
end

##  ----------------------------------------------------------------------------
#   In-place functions to compute any code's price index from the lower level data
#   ----------------------------------------------------------------------------

function compute_index!(cache::Dict, good::Item, base::FullCPIBase)
    i = findfirst(==(good.code), base.codes)
    cache[good.code] = base.ipc[:, i] # save a copy in the cache
    cache[good.code]
end

function compute_index!(cache::Dict, group::Group, base::FullCPIBase)
    # If code is available in cache, just return it
    group.code in keys(cache) && return cache[group.code]

    # Else, compute the index and store it
    # Get the indexes of the children. At the lowest level the dispatch will
    # select compute_index(::Item, ::FullCPIBase) to return the Goods indices
    ipcs = mapreduce(c -> compute_index!(cache, c, base), hcat, group.children)

    # If there exists only one good in the group, that is the group's index
    if size(ipcs, 2) == 1
        cache[group.code] = ipcs
        return ipcs
    end

    # Get the weights
    weights = map(c -> c.weight, group.children)

    # Store normalized sum product
    cache[group.code] = ipcs * weights / sum(weights)
    cache[group.code]
end

function compute_index!(::Dict, ::Nothing, ::FullCPIBase)
    @warn "Nodo no disponible en la estructura"
    nothing
end

##  ----------------------------------------------------------------------------
#   CPITree: estructura contenedora del árbol jerárquico y los datos necesarios
#   para computar cualquier nodo en la estructura del IPC.
#   ----------------------------------------------------------------------------

"""
    CPITree

    CPITree(base::FullCPIBase, tree::Union{Group, Item}, group_names::Vector{String}, group_codes::Vector{String})
    CPITree(; base::FullCPIBase, groupsdf::DataFrame, characters::(NTuple{N, Int} where N), upperlevel_code = "_0", upperlevel_name = "IPC")

Contenedor envolvente de un árbol jerárquico del IPC y los datos necesarios de
los gastos básicos para computar cualquier jerarquía dentro del árbol. Permite
visualizar y explorar la composición del IPC de un país, así como computar los
índices de precios de las diferentes jerarquías de la estructura del IPC. Está
compuesto por:
- Un objeto `base`, de tipo [`FullCPIBase`](@ref), el cual almacena las series
  de tiempo de los índices de los gastos básicos.
- Un objeto `tree` que contiene la estructura de nodos del IPC.
- El vector de nombres `group_names` de los grupos del árbol `tree`.
- El vector de códigos `group_codes` de los grupos del árbol `tree`.

El constructor simple requiere una estructura jerárquica de nodos como la
devuelta por [`get_cpi_tree`](@ref). Al utilizar el constructor con `groupsdf` y
`characters`, se construye automáticamente un `CPITree` utilizando los códigos
como indicadores de la estructura jerárquica. Los códigos de los gastos básicos
contenidos en `base` describen cómo se agrupan las jerarquías. Por ejemplo, el
siguiente `FullCPIBase` contiene 10 gastos básicos:
```julia-repl
julia> base
FullCPIBase{Float32, Float32}: 36 períodos × 10 gastos básicos Jan-01-Dec-03
┌─────┬────────────┬─────────┬─────────┬─────────┬─────────┬─────────┬─────────┬─────────┬─────────┬─────────┬─────────┐
│ Row │      Dates │ _011101 │ _011201 │ _011202 │ _021101 │ _022101 │ _031101 │ _041101 │ _041201 │ _041202 │ _041301 │
│     │            │  Item A │  Item B │  Item C │  Item D │  Item E │  Item F │  Item G │  Item H │  Item I │  Item J │
│     │            │ 11.8947 │ 2.39931 │ 7.80836 │ 1.58585 │  20.141 │ 12.0526 │ 9.56901 │ 14.6256 │ 15.0276 │ 4.89602 │
├─────┼────────────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┤
│   1 │ 2001-01-01 │  100.00 │  100.00 │  100.00 │  100.00 │  100.00 │  101.13 │  100.00 │  100.00 │  100.51 │  100.00 │
│   2 │ 2001-02-01 │  100.00 │  100.00 │  100.00 │  100.00 │  100.00 │  102.28 │  100.00 │  100.00 │  101.02 │  100.00 │
│   3 │ 2001-03-01 │  100.00 │  100.00 │  100.00 │  100.00 │  100.00 │  103.44 │  100.00 │  100.00 │  101.53 │  100.00 │
│   4 │ 2001-04-01 │  100.00 │  100.00 │  100.00 │  100.00 │  100.00 │  104.61 │  100.00 │  100.00 │  102.05 │  100.00 │
│   5 │ 2001-05-01 │  100.00 │  100.00 │  100.00 │  100.00 │  100.00 │  105.79 │  100.00 │  100.00 │  102.56 │  100.00 │
│  ⋮  │     ⋮      │    ⋮    │    ⋮    │    ⋮    │    ⋮    │    ⋮    │    ⋮    │    ⋮    │    ⋮    │    ⋮    │    ⋮    │
│  32 │ 2003-08-01 │  100.00 │  100.00 │  100.00 │  100.00 │  100.00 │  143.39 │  100.00 │  100.00 │  117.59 │  100.00 │
│  33 │ 2003-09-01 │  100.00 │  100.00 │  100.00 │  100.00 │  100.00 │  145.02 │  100.00 │  100.00 │  118.19 │  100.00 │
│  34 │ 2003-10-01 │  100.00 │  100.00 │  100.00 │  100.00 │  100.00 │  146.66 │  100.00 │  100.00 │  118.79 │  100.00 │
│  35 │ 2003-11-01 │  100.00 │  100.00 │  100.00 │  100.00 │  100.00 │  148.32 │  100.00 │  100.00 │  119.39 │  100.00 │
│  36 │ 2003-12-01 │  100.00 │  100.00 │  100.00 │  100.00 │  100.00 │  150.00 │  100.00 │  100.00 │  120.00 │  100.00 │
└─────┴────────────┴─────────┴─────────┴─────────┴─────────┴─────────┴─────────┴─────────┴─────────┴─────────┴─────────┘
                                                                                26 rows omitted
```

En este ejemplo, el argumento debe especificarse como `characters = (3, 4, 5,
7)`, pues los códigos indican las jerarquías en el IPC de manera siguiente:
- Los primeros 3 caracteres indican la jerarquía de *división de gasto*.
- El siguiente caracter indica la jerarquía de *agrupación de gasto*.
- El siguiente caracter indica la jerarquía de *subgrupo de gasto*.
- Los siguientes 2 caracteres indican el *número de gasto básico dentro de su
  grupo*.

Por su parte, el DataFrame `groupsdf` debe tener la estructura mínima siguiente:
- La primera columna debe ser de tipo `String` y contiene los códigos de los
  grupos disponibles en la estructura del IPC.
- La segunda columna debe ser de tipo `String` y contiene las descripciones o
nombres de los grupos disponibles en la estructura del IPC. Por ejemplo, el
DataFrame `groupsdf` puede verse de esta forma:
```
17×2 DataFrame
 Row │ code    description
     │ String  String
─────┼─────────────────────
   1 │ _01     Div._01
   2 │ _02     Div._02
   3 │ _03     Div._03
   4 │ _04     Div._04
   5 │ _011    Agr._011
   6 │ _021    Agr._021
   7 │ _022    Agr._022
   8 │ _031    Agr._031
   9 │ _041    Agr._041
  10 │ _0111   Subgr._0111
  11 │ _0112   Subgr._0112
  12 │ _0211   Subgr._0211
  13 │ _0221   Subgr._0221
  14 │ _0311   Subgr._0311
  15 │ _0411   Subgr._0411
  16 │ _0412   Subgr._0412
  17 │ _0413   Subgr._0413
```

Por ejemplo, al construir un `CPITree` con la estructura de códigos indicada anteriormente y el DataFrame de ejemplo, el árbol del IPC se puede ver de esta forma:
```julia-repl
julia> tree = CPITree(; base, groupsdf, characters=(3,4,5,7))
CPITree{Group{Group{Group{Group{Item{Float32}, Float32}, Float32}, Float32}, Float32}} con datos
└─→ FullCPIBase{Float32, Float32}: 36 períodos × 10 gastos básicos Jan-01-Dec-03
_0: IPC [100.0]
├─ _01: Div._01 [21.491905]
│  └─ _011: Agr._011 [21.491905]
│     ├─ _0111: Subgr._0111 [7.352945]
│     │  └─ _011101: Item A [7.352945]
│     └─ _0112: Subgr._0112 [14.13896]
│        ├─ _011201: Item B [6.7442417]
│        └─ _011202: Item C [7.394718]
├─ _02: Div._02 [3.0530455]
│  ├─ _021: Agr._021 [1.1036392]
│  │  └─ _0211: Subgr._0211 [1.1036392]
│  │     └─ _021101: Item D [1.1036392]
│  └─ _022: Agr._022 [1.9494063]
│     └─ _0221: Subgr._0221 [1.9494063]
│        └─ _022101: Item E [1.9494063]
├─ _03: Div._03 [11.68543]
│  └─ _031: Agr._031 [11.68543]
│     └─ _0311: Subgr._0311 [11.68543]
│        └─ _031101: Item F [11.68543]
└─ _04: Div._04 [63.769615]
   └─ _041: Agr._041 [63.769615]
      ├─ _0411: Subgr._0411 [16.103952]
      │  └─ _041101: Item G [16.103952]
      ├─ _0412: Subgr._0412 [28.824577]
      │  ├─ _041201: Item H [11.367162]
      │  └─ _041202: Item I [17.457417]
      └─ _0413: Subgr._0413 [18.841085]
         └─ _041301: Item J [18.841085]
```
"""
struct CPITree{G}
    base::FullCPIBase
    tree::G
    group_names::Vector{String}
    group_codes::Vector{String}
    function CPITree(base::FullCPIBase, tree::Union{Group, Item}, group_names::Vector{String}, group_codes::Vector{String})
        G = typeof(tree)
        new{G}(base, tree, group_names, group_codes)
    end
end

function CPITree(;
    base::FullCPIBase,
    groupsdf::DataFrame,
    characters::(NTuple{N, Int} where N),
    upperlevel_code = "_0",
    upperlevel_name = "IPC")

    # Obtener códigos y nombres de los grupos en las jerarquías
    group_codes = convert.(String, groupsdf[!, 1])
    group_names = groupsdf[!, 2]

    tree = get_cpi_tree(;
        base,
        characters,
        group_names,
        group_codes,
        upperlevel_code,
        upperlevel_name
    )

    CPITree(base, tree, group_names, group_codes)
end

function Base.show(io::IO, cpitree::CPITree)
    println(io, typeof(cpitree), " con datos")
    println(io, "└─→ ", sprint(summary, cpitree.base))
    print_tree(io, cpitree.tree)
end

"""
    Base.getindex(cpitree::CPITree, code::AbstractString)

Este método se utiliza para indexar el árbol jerárquico `cpitree` y obtener una
estructura similar cuyo nodo superior sea el nodo con el código provisto `code`. Por ejemplo, si tenemos el siguiente árbol:
```julia-repl
julia> tree
CPITree{Group{Group{Group{Group{Item{Float32}, Float32}, Float32}, Float32}, Float32}} con datos
└─→ FullCPIBase{Float32, Float32}: 36 períodos × 10 gastos básicos Jan-01-Dec-03
_0: IPC [100.0]
├─ _01: Div._01 [21.491905]
│  └─ _011: Agr._011 [21.491905]
│     ├─ _0111: Subgr._0111 [7.352945]
│     │  └─ _011101: Item A [7.352945]
│     └─ _0112: Subgr._0112 [14.13896]
│        ├─ _011201: Item B [6.7442417]
│        └─ _011202: Item C [7.394718]
├─ _02: Div._02 [3.0530455]
│  ├─ _021: Agr._021 [1.1036392]
│  │  └─ _0211: Subgr._0211 [1.1036392]
│  │     └─ _021101: Item D [1.1036392]
│  └─ _022: Agr._022 [1.9494063]
│     └─ _0221: Subgr._0221 [1.9494063]
│        └─ _022101: Item E [1.9494063]
├─ _03: Div._03 [11.68543]
│  └─ _031: Agr._031 [11.68543]
│     └─ _0311: Subgr._0311 [11.68543]
│        └─ _031101: Item F [11.68543]
└─ _04: Div._04 [63.769615]
   └─ _041: Agr._041 [63.769615]
      ├─ _0411: Subgr._0411 [16.103952]
      │  └─ _041101: Item G [16.103952]
      ├─ _0412: Subgr._0412 [28.824577]
      │  ├─ _041201: Item H [11.367162]
      │  └─ _041202: Item I [17.457417]
      └─ _0413: Subgr._0413 [18.841085]
         └─ _041301: Item J [18.841085]
```

Al indexar por un código, como `_041`, obtenemos una estructura similar a partir de ese nodo:
```julia-repl
julia> tree["_041"]
CPITree{Group{Group{Item{Float32}, Float32}, Float32}} con datos
└─→ FullCPIBase{Float32, Float32}: 36 períodos × 10 gastos básicos Jan-01-Dec-03
_041: Agr._041 [63.769615]
├─ _0411: Subgr._0411 [16.103952]
│  └─ _041101: Item G [16.103952]
├─ _0412: Subgr._0412 [28.824577]
│  ├─ _041201: Item H [11.367162]
│  └─ _041202: Item I [17.457417]
└─ _0413: Subgr._0413 [18.841085]
   └─ _041301: Item J [18.841085]
```
"""
function Base.getindex(cpitree::CPITree, code::AbstractString)
    node = cpitree.tree[code]
    node === nothing && return nothing
    CPITree(cpitree.base, node, cpitree.group_names, cpitree.group_codes)
end

"""
    compute_index(cpitree::CPITree [, code::AbstractString])

Permite computar el índice de precios de la jerarquía provista en `code`. Si se
omite `code`, se computa la jerarquía padre de la estructura `cpitree`. Si el
nodo no se encuentra en el árbol, devuelve `nothing`.

```julia-repl
julia> tree
CPITree{Group{Group{Group{Group{Item{Float32}, Float32}, Float32}, Float32}, Float32}} con datos
└─→ FullCPIBase{Float32, Float32}: 36 períodos × 10 gastos básicos Jan-01-Dec-03
_0: IPC [100.0]
├─ _01: Div._01 [21.491905]
│  └─ _011: Agr._011 [21.491905]
│     ├─ _0111: Subgr._0111 [7.352945]
│     │  └─ _011101: Item A [7.352945]
│     └─ _0112: Subgr._0112 [14.13896]
│        ├─ _011201: Item B [6.7442417]
│        └─ _011202: Item C [7.394718]
├─ _02: Div._02 [3.0530455]
│  ├─ _021: Agr._021 [1.1036392]
│  │  └─ _0211: Subgr._0211 [1.1036392]
│  │     └─ _021101: Item D [1.1036392]
│  └─ _022: Agr._022 [1.9494063]
│     └─ _0221: Subgr._0221 [1.9494063]
│        └─ _022101: Item E [1.9494063]
├─ _03: Div._03 [11.68543]
│  └─ _031: Agr._031 [11.68543]
│     └─ _0311: Subgr._0311 [11.68543]
│        └─ _031101: Item F [11.68543]
└─ _04: Div._04 [63.769615]
   └─ _041: Agr._041 [63.769615]
      ├─ _0411: Subgr._0411 [16.103952]
      │  └─ _041101: Item G [16.103952]
      ├─ _0412: Subgr._0412 [28.824577]
      │  ├─ _041201: Item H [11.367162]
      │  └─ _041202: Item I [17.457417]
      └─ _0413: Subgr._0413 [18.841085]
         └─ _041301: Item J [18.841085]

julia> compute_index(tree, "_041")
36-element Vector{Float32}:
 100.13899
 100.2787
 100.41912
 100.560234
 100.70207
 100.84464
 100.98792
   ⋮
 104.65377
 104.81638
 104.97984
 105.14412
 105.309235
 105.47519

julia> compute_index(tree["_041"])
36-element Vector{Float32}:
 100.13899
 100.2787
 100.41912
 100.560234
 100.70207
 100.84464
 100.98792
   ⋮
 104.65377
 104.81638
 104.97984
 105.14412
 105.309235
 105.47519

julia> a = tree["_041302"]

julia> a === nothing
true
```
"""
function compute_index(cpitree::CPITree, code::AbstractString)
    node = cpitree.tree[code]
    node === nothing && return nothing
    compute_index(node, cpitree.base)
end

# When single argument, compute the top level node
function compute_index(cpitree::CPITree)
    node = cpitree.tree
    compute_index(node, cpitree.base)
end

children(cpitree::CPITree) = children(cpitree.tree)


text:
import LinearAlgebra: Schur, GeneralizedSchur

export Near, SmallestMagnitude, LargestMagnitude,
       LargestRealPart, SmallestRealPart,
       LargestImaginaryPart, SmallestImaginaryPart,
       schur_permutation, schur_sort!

abstract type Target end

# For finding eigenvalues near a target
struct Near{T} <: Target
  τ::Complex{T}
end

# For finding eigenvalues with the smallest magnitude
struct SmallestMagnitude{T} <: Target
  τ::Complex{T}
end

# For finding eigenvalues with the largest magnitude
struct LargestMagnitude{T} <: Target
  τ::Complex{T}
end

# For finding eigenvalues with the largest real part
struct LargestRealPart{T} <: Target
  τ::Complex{T}
end

# For finding eigenvalues with the smallest real part
struct SmallestRealPart{T} <: Target
  τ::Complex{T}
end

# For finding eigenvalues with the largest imaginary part
struct LargestImaginaryPart{T} <: Target
  τ::Complex{T}
end

# For finding eigenvalues with the smallest imaginary part
struct SmallestImaginaryPart{T} <: Target
  τ::Complex{T}
end

eigval(F::GeneralizedSchur, idx::Int) = F.alpha[idx] / F.beta[idx]

schur_permutation(target::Near, F::GeneralizedSchur, k::OrdinalRange) =
  partialsortperm(1 : length(F.alpha), k, by = idx -> abs(eigval(F, idx) - target.τ))

schur_permutation(target::LargestMagnitude, F::GeneralizedSchur, k::OrdinalRange) =
  partialsortperm(1 : length(F.alpha), k, by = idx -> abs(eigval(F, idx)), rev = true)

schur_permutation(target::SmallestMagnitude, F::GeneralizedSchur, k::OrdinalRange) =
  partialsortperm(1 : length(F.alpha), k, by = idx -> abs(eigval(F, idx)))

schur_permutation(target::LargestRealPart, F::GeneralizedSchur, k::OrdinalRange) =
  partialsortperm(1 : length(F.alpha), k, by = idx -> real(eigval(F, idx)), rev = true)

schur_permutation(target::SmallestRealPart, F::GeneralizedSchur, k::OrdinalRange) =
  partialsortperm(1 : length(F.alpha), k, by = idx -> real(eigval(F, idx)))

schur_permutation(target::LargestImaginaryPart, F::GeneralizedSchur, k::OrdinalRange) =
  partialsortperm(1 : length(F.alpha), k, by = idx -> imag(eigval(F, idx)), rev = true)

schur_permutation(target::SmallestImaginaryPart, F::GeneralizedSchur, k::OrdinalRange) =
  partialsortperm(1 : length(F.alpha), k, by = idx -> imag(eigval(F, idx)))

function schur_sort!(target::Target, F::GeneralizedSchur, k::OrdinalRange)
  best = vec(schur_permutation(target, F, k))
  perm = falses(length(F.alpha))
  perm[best] .= true
  ordschur!(F, perm)
end

schur_sort!(target::Target, F::GeneralizedSchur, k::Int) =
    schur_sort!(target, F, k:k)

ext: jl

max_forks_repo_name: haampie/jacobi-davidson-julia


text:
using ReproductionNumbers

# setup
fig_name    =   "reported"
file_name   =   "../data/raw/reported-cases-data.csv" # this file is huge and is hence added to .gitignore
days_col    =   "Refdatum"
data_col    =   "AnzahlFall_sum"
ylabel_R    =   "Reproduction number R"
ylabel_N    =   "Reported cases"
# The reported cases by RKI are changed and updated daily
# That means, cases that were not submitted on time will be added the day after etc.
# This leads to today's and yesterday's values of R and N being rather useless
# Hence, the variable `pop_days` allows to go back as many days from today as we believe are credible
# Looking at the RKI dashboard we find that today - 2 days is fine.
pop_days    =   2

N, R, df_cases = main_reported(fig_name, file_name, days_col, data_col, ylabel_R, ylabel_N, "country", pop_days)


text:
using Intensity, Random
using Test

@testset "Intensity.jl" begin
    @testset "TermStructure" begin
        terms = [0.5, 1.0, 2.0]
        values = [1.0, 2.0, 3.0]
        test_terms = [0.0, 1.5, 2.0, 3.0]
        test_pwl_values = [0.0, 2.5
        , 3.0, 3.0]
        test_const_values = [1.0, 3.0, 3.0, 3.0]
        ts_pwl = Intensity.TermStructure(terms, values, Intensity.PWL)
        ts_const = Intensity.TermStructure(terms, values, Intensity.Constant)
        combs = zip(
            [Intensity.PWL, Intensity.Constant],
            [test_pwl_values, test_const_values],
            [ts_pwl, ts_const],
        )
        for (interptype, test_vals, ts) in combs
            @testset "$(ts.interpolation)" begin
                for (i, term) in enumerate(test_terms)
                    pred = Intensity.interextrapolate(ts, term)
                    @test pred ≈ test_vals[i] atol=1e-10
                end
            end
        end

        @testset "SurvivalProbability" begin
            terms = [1.0, 0.0, 2.0]
            probs = [1.0, 0.5, 0.4]
            @test_throws DomainError Intensity.SurvivalProbability(terms, probs)
            terms = [0.0, 1.0, 2.0]
            @test isa(Intensity.SurvivalProbability(terms, probs), Intensity.AbstractTermStructure)
            for probs in [[1.0, 0.5], [1.0, 0.4, 0.5], [1.0, 0.5, -0.5], [1.0, 1.1, 0.4]]
                @test_throws DomainError Intensity.SurvivalProbability(terms, probs)
            end
        end

        @testset "DeterministicIntensity" begin
            terms = [0.0, 1.0, 2.0]
            intensities = [0.5, 0.4, 0.3]
            valid = Intensity.DeterministicIntensity(terms, intensities)
            @test isa(valid, Intensity.DeterministicIntensity)
            @test valid.interpolation === Intensity.PWL
            for intensities in [[0.5, 0.4], [0.5, 0.4, -0.3]]
                @test_throws DomainError Intensity.DeterministicIntensity(terms, intensities)
            end
            pred = Intensity.interextrapolate(valid, 0.0)
            @test pred == 0.5
            pred = Intensity.interextrapolate(valid, terms)
            @test pred == intensities
            @test Intensity.interextrapolate(valid, 3.0) == 0.3

            @testset "probabilities" begin
                terms = [0.0, 1.0, 2.0]
                intensities = [0.5, 0.4, 0.3]
                const_probs = exp.(-cumsum([0.0, 0.4, 0.3]))
                pwl_probs = exp.(-cumsum([0.0, 0.45, 0.35]))
                test_terms = [0.0, 0.5, 1.0, 1.5, 2.0]
                for (probs, itype) in zip([const_probs, pwl_probs], [Intensity.Constant, Intensity.PWL])
                    @testset "$itype" begin
                        di = Intensity.DeterministicIntensity(terms, intensities, itype)
                        pred = Intensity.get_rolling_probabilities(di)
                        @test probs == pred
                        t_probs = Intensity.get_survival_probability(di, test_terms)
                        @test probs == t_probs[[1, 3, 5]]
                        sp = Intensity.SurvivalProbability(test_terms, t_probs)
                        @test isa(sp, Intensity.AbstractTermStructure)
                        di_2 = Intensity.calibrate_intensity(sp, itype)
                        @test isa(di_2, Intensity.DeterministicIntensity)
                        t_probs2 = Intensity.get_survival_probability(di, test_terms)
                        @test t_probs == t_probs2
                    end
                end
            end

            @testset "InhomogeneousExponential" begin
                IE = Intensity.InhomogeneousExponential
                ts = [0.0, 0.50, 1.00, 1.50, 2.00, 2.50, 3.00, 4.00, 5.00]
                obs_surv_prob = [1.0, 0.99, 0.95, 0.88, 0.80, 0.75, 0.70, 0.60, 0.45]
                sp = Intensity.SurvivalProbability(ts, obs_surv_prob)
                dc = Intensity.calibrate_intensity(sp, Intensity.PWL)
                ie = IE(dc, 0:0.01:5)
                rng = MersenneTwister(31)
                shape = (4, 3)
                dts = rand(rng, ie, shape)
                rng = MersenneTwister(31)
                dts2 = rand(rng, ie, shape)
                @test size(dts) == shape
                @test dts == dts2
                #surv_prob = Intensity.get_survival_probability(ie.intensity, ie.grid)
                #emp_surv_prob = [1 - mean(dts .< t) for t in ie.grid]
            end
        end
    end
end

max_stars_repo_name: lsindoni/Intensity

max_issues_repo_name: lsindoni/Intensity

max_forks_repo_name: lsindoni/Intensity


text:
bottles(n) = n==0 ? "No more bottles" :
             n==1 ? "1 bottle" :
             "$n bottles"

for n = 99:-1:1
    println("""
        $(bottles(n)) of beer on the wall
        $(bottles(n)) of beer
        Take one down, pass it around
        $(bottles(n-1)) of beer on the wall
    """)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module NonPromoting

export NP
struct NP{T} <: AbstractFloat where {T <: AbstractFloat}
    val::T
    NP{T}(x::T) where {T <: AbstractFloat} = new{T}(x)
end

NP{T}(x::I) where {T <: AbstractFloat, I <: Integer} = NP{T}(T(x))
NP{T}(x::Rational{I}) where {T <: AbstractFloat, I <: Integer} = NP{T}(T(x))
NP{T}(x::Irrational{I}) where {T <: AbstractFloat, I} = NP{T}(T(x))
NP(x::T) where {T <: AbstractFloat} = NP{T}(x)

function Base.convert(::Type{T}, x::NP{T})::T where {T <: AbstractFloat}
    x.val
end

function Base.show(io::IO, x::NP{T}) where {T}
    show(io, x.val)
end

# Prevent promotion
Base.promote_rule(::Type{NP{T}}, ::Type{NP{T}}) where {T} = NP{T}
Base.promote_rule(::Type{NP{T}}, ::Type{NP{U}}) where {T, U} = error
Base.promote_rule(::Type{NP{T}}, ::Type) where {T} = error
Base.promote_rule(::Type, ::Type{NP{T}}) where {T} = error

# Nullary functions (which take a type as argument)
for fun in [:one, :zero]
    @eval begin
        function Base.$fun(::Type{NP{T}})::NP{T} where {T <: AbstractFloat}
            NP{T}($fun(T))
        end
    end
end

# Unary functions
for fun in [:(-),
            :abs, :acos, :acosh, :asin, :asinh, :atan, :atanh, :cbrt,
            :cos, :cosh, :exp, :exp10, :exp2,
            :inv, :log, :log10, :log2,
            :sign,:sin, :sinh, :sqrt, :tan, :tanh]
    @eval begin
        function Base.$fun(x::NP{T})::NP{T} where {T <: AbstractFloat}
            NP{T}($fun(x.val))
        end
    end
end

# Unary functions that return another type
for fun in [:isfinite, :isinf, :isnan, :issubnormal, :signbit]
    @eval begin
        function Base.$fun(x::NP{T}) where {T <: AbstractFloat}
            $fun(x.val)
        end
    end
end

# Binary functions
for fun in [:(-), :(/), :(\), :(^),
            :atan, :copysign, :hypot, :modf, :rem]
    @eval begin
        function Base.$fun(x::NP{T}, y::NP{T})::NP{T} where {T <: AbstractFloat}
            NP{T}($fun(x.val, y.val))
        end
    end
end

# Binary functions that return another type
for fun in [:(==), :(<)]
    @eval begin
        function Base.$fun(x::NP{T}, y::NP{T}) where {T <: AbstractFloat}
            $fun(x.val, y.val)
        end
    end
end

# n-ary functions
for fun in [:(+), :(*), :max, :min]
    @eval begin
        function Base.$fun(x::NP{T}, ys::NP{T}...)::NP{T} where
                {T <: AbstractFloat}
            NP{T}($fun(x.val, map(y -> y.val, ys)...))
        end
    end
end

# Functions that take additional arguments
# function Base.norm(x::NP{T}, p::Real=2)::NP{T} where {T <: AbstractFloat}
#     NP{T}(norm(x.val, p))
# end

end

max_stars_repo_name: eschnett/NonPromoting

max_issues_repo_name: eschnett/NonPromoting

max_forks_repo_name: eschnett/NonPromoting


text:
type NfToNfMor <: Map{AnticNumberField, AnticNumberField}
  header::MapHeader{AnticNumberField, AnticNumberField}
  prim_img::nf_elem

  function NfToNfMor()
    z = new()
    z.header = MapHeader()
    return r
  end

  function NfToNfMor(K::AnticNumberField, L::AnticNumberField, y::nf_elem)
    z = new()
    z.prim_img = y

    function image(x::nf_elem)
      g = parent(K.pol)(x)
      return evaluate(g, y)
    end

    z.header = MapHeader(K, L, image)
    return z
  end
end

type NfMaxOrdToFqNmodMor <: Map{NfMaxOrd, FqNmodFiniteField}
  header::MapHeader{NfMaxOrd, FqNmodFiniteField}
  poly_of_the_field::nmod_poly

  function NfMaxOrdToFqNmodMor()
    r = new()
    r.header = MapHeader{NfMaxOrd, FqNmodFiniteField}()
    return r
  end

  function NfMaxOrdToFqNmodMor(O::NfMaxOrd, F::FqNmodFiniteField, g::nmod_poly)
    # assume that F = F_p[X]/(g) and g is a factor of f mod p

    z = new()
    p = characteristic(F)
    a = gen(nf(O))
    tmp_nmod_poly = parent(g)()

    z.poly_of_the_field = g

    function _image(x::NfOrdElem)
      u = F()
      gg = parent(nf(O).pol)(elem_in_nf(x))::fmpq_poly
      fmpq_poly_to_nmod_poly_raw!(tmp_nmod_poly, gg)
      ccall((:fq_nmod_set, :libflint), Void, (Ptr{fq_nmod}, Ptr{nmod_poly}, Ptr{FqNmodFiniteField}), &u, &tmp_nmod_poly, &F)
      ccall((:fq_nmod_reduce, :libflint), Void, (Ptr{fq_nmod}, Ptr{FqNmodFiniteField}), &u, &F)
      return u
    end

    function _preimage(x::fq_nmod)
      zz = nf(O)()

      for i in 0:degree(F)-1
        zz = zz + _get_coeff_raw(x, i)*a^i
      end

      return O(zz, false)
    end

    z.header = MapHeader{NfMaxOrd, FqNmodFiniteField}(O, F, _image, _preimage)

    return z
  end

  function NfMaxOrdToFqNmodMor(O::NfMaxOrd, F::FqNmodFiniteField, y::fq_nmod)
    z = new()

    p = characteristic(F)
    Zx = PolynomialRing(FlintIntegerRing(), "x")[1]
    a = gen(nf(O))
    h = Zx()
    t_fq_nmod = F()
    tt_fq_nmod = F()
    t_fmpz = fmpz()

    function _image(x::NfOrdElem)
      g = parent(nf(O).pol)(elem_in_nf(x))

      #pseudocode:
      #u = inv(F(den(g)))
      #return u*evaluate(num(g), y)

      ccall((:fmpq_poly_get_denominator, :libflint), Void,
            (Ptr{fmpz}, Ptr{fmpq_poly}), &t_fmpz, &g)

      ccall((:fq_nmod_set_fmpz, :libflint), Void,
            (Ptr{fq_nmod}, Ptr{fmpz}, Ptr{FqNmodFiniteField}),
            &tt_fq_nmod, &t_fmpz, &F)

      ccall((:fq_nmod_inv, :libflint), Void,
            (Ptr{fq_nmod}, Ptr{fq_nmod}, Ptr{FqNmodFiniteField}),
            &tt_fq_nmod, &tt_fq_nmod, &F)

      ccall((:fmpq_poly_get_numerator, :libflint), Void,
                  (Ptr{fmpz_poly}, Ptr{fmpq_poly}), &h, &g)

      evaluate!(t_fq_nmod, h, y)
      #@assert t_fq_nmod == evaluate(h, y)
      return tt_fq_nmod * t_fq_nmod
    end

    function _preimage(x::fq_nmod)
      z = nf(O)()

      for i in 0:degree(F)-1
        z = z + _get_coeff_raw(x, i)*a^i
      end

      return O(z, false)
    end

    z.header = MapHeader{NfMaxOrd, FqNmodFiniteField}(O, F, _image, _preimage)

    return z
  end
end

type NfMaxOrdQuoMap <: Map{NfMaxOrd, NfMaxOrdQuoRing}
  header::MapHeader{NfMaxOrd, NfMaxOrdQuoRing}

  function NfMaxOrdQuoMap(O::NfMaxOrd, Q::NfMaxOrdQuoRing)
    z = new()

    _image = function (x::NfOrdElem)
      return Q(x)
    end

    _preimage = function (x::NfMaxOrdQuoRingElem)
      return x.elem
    end

    z.header = MapHeader(O, Q, _image, _preimage)
    return z
  end
end

function Mor(O::NfMaxOrd, F::FqNmodFiniteField, y::fq_nmod)
  return NfMaxOrdToFqNmodMor(O, F, y)
end

function Mor(O::NfMaxOrd, F::FqNmodFiniteField, h::nmod_poly)
  return NfMaxOrdToFqNmodMor(O, F, h)
end

type NfToFqNmodMor <: Map{AnticNumberField, FqNmodFiniteField}
  header::MapHeader

  function NfToFqNmodMor()
    r = new()
    r.header = MapHeader{AnticNumberField, FqNmodFiniteField}()
    return r
  end
end

function extend(f::NfMaxOrdToFqNmodMor, K::AnticNumberField)
  nf(domain(f)) != K && error("Number field is not the number field of the order")

  z = NfToFqNmodMor()
  z.header.domain = K
  z.header.codomain = f.header.codomain

  p = characteristic(z.header.codomain)
  Zx = PolynomialRing(FlintIntegerRing(), "x")[1]
  y = f(NfOrdElem{NfMaxOrd}(domain(f), gen(K)))

  function _image(x::nf_elem)
    g = parent(K.pol)(x)
    u = inv(z.header.codomain(den(g)))

    g = Zx(den(g)*g)
    return u*evaluate(g, y)
  end

  z.header.image = _image
  return z
end

function evaluate(f::fmpz_poly, r::fq_nmod)
  #Horner - stolen from Claus

  if length(f) == 0
    return parent(r)()
  end

  l = f.length-1
  s = parent(r)(coeff(f, l))
  for i =l-1:-1:0
    s = s*r + parent(r)(coeff(f, i))
  end
  return s
end

function evaluate!(z::fq_nmod, f::fmpz_poly, r::fq_nmod)
  #Horner - stolen from Claus

  zero!(z)

  if length(f) == 0
    return z
  end

  l = f.length-1
  set!(z, parent(r)(coeff(f, l)))
  #s = parent(r)(coeff(f, l))
  for i =l-1:-1:0
    mul!(z, z, r)
    add!(z, z, parent(r)(coeff(f, i)))
    #s = s*r + parent(r)(coeff(f, i))
  end
  return z
end

function Mor(K::AnticNumberField, L::AnticNumberField, y::nf_elem)
  z = NfToNfMor(K, L, y)
  return z
end

function _get_coeff_raw(x::fq_nmod, i::Int)
  u = ccall((:nmod_poly_get_coeff_ui, :libflint), UInt, (Ptr{fq_nmod}, Int), &x, i)
  return u
end

function _get_coeff_raw(x::fq, i::Int)
  t = FlintIntegerRing()
  ccall((:fmpz_poly_get_coeff_fmpz, :libflint), Void, (Ptr{fmpz}, Ptr{fq}, Int), &t, &x, i)
  return t
end

function call(f::NfMaxOrdToFqNmodMor, p::PolyElem{NfOrdElem{NfMaxOrd}})
  F = codomain(f)
  Fx,_ = PolynomialRing(F, "_\$")

  ar = NfOrdElem[ coeff(p, i) for i in 0:degree(p) ]

  z = Fx(map(f, ar))

  return z
end

# this is super slow
# improve
function call(f::NfMaxOrdQuoMap, I::NfMaxOrdIdl)
  O = domain(f)
  Q = codomain(f)
  B = Q.ideal + I
  b = basis(B)

  z = O()

  while true
    z = rand(fmpz(1):norm(Q.ideal)^2) * b[1]

    for i in 2:degree(O)
      z = z + rand(fmpz(1):norm(Q.ideal)^2) * b[i]
    end

    if norm(ideal(O, z) + ideal(O, O(norm(Q.ideal)))) == norm(B)
      break
    end
  end

  return Q(z)
end

function call{T}(f::NfMaxOrdQuoMap, p::PolyElem{NfOrdElem{T}})
  F = codomain(f)
  Fx,_ = PolynomialRing(F, "_\$")

  ar = NfOrdElem[ coeff(p, i) for i in 0:degree(p) ]

  z = Fx(map(f, ar))

  return z
end

base_ring(::NfMaxOrd) = Union{}

Nemo.needs_parentheses(::NfOrdElem) = true

Nemo.is_negative(::NfOrdElem) = false


text:
using GymSpaces: Box, Discrete
using ArcadeLearningEnvironment

setLoggerMode!(:error)  # Only log errors

struct AtariEnv <: AbstractEnv
    game_path::String
    ale::Ptr{Nothing}
    action_set::Array{UInt, 1}
    action_space::Discrete
    observation_space::Box
    frameskip::Union{Int, UnitRange{Int}}
    obs_type::Symbol
end

include("vis/utils.jl")

function AtariEnv( ; game_path::String, obs_type::Symbol=:ram, frameskip::Union{Int, UnitRange{Int}}=2:5,
    repeat_action_probability::AbstractFloat=0f0, full_action_space::Bool=false)

    @assert obs_type ∈ (:ram, :color, :grey) "Invalid observation type. Choose either :ram, :image or :grey"
    @assert isa(repeat_action_probability, AbstractFloat) || isa(repeat_action_probability, Integer) "Invalid repeat_action_probability"

    ale = ALE_new()
    loadROM(ale, game_path)  # This is will also check if the loaded file is valid or not.
    action_set = full_action_space ? getLegalActionSet(ale) : getMinimalActionSet(ale)
    #action_set .+= 1  # Converting to one indexed array
    action_space = full_action_space ? Discrete(getLegalActionSize(ale)) : Discrete(getMinimalActionSize(ale))

    setFloat(ale, "repeat_action_probability", repeat_action_probability)

    screen_width, screen_height = getScreenWidth(ale), getScreenHeight(ale)

    observation_space = obs_type==:ram ? Box(0, 255, (128,), UInt8) : Box(0, 255, (screen_height, screen_width, 3), UInt8)

    AtariEnv(game_path, ale, action_set, action_space, observation_space, frameskip, obs_type)
end

function step!(env::AtariEnv, action)
    @assert action ∈ env.action_space "Action $action is invalid."
    reward = 0f0

    num_steps = isa(env.frameskip, UnitRange) ? rand(env.frameskip) : env.frameskip

    for _ in 1:num_steps
        reward += act(env.ale, env.action_set[action])
    end
    ob = _get_obs(env)
    return ob, reward, ArcadeLearningEnvironment.game_over(env.ale), Dict(:ale_lives => lives(env.ale))
end

function get_preprocessed_RGB(env::AtariEnv)
    # getScreenRGB returns an `Array{UInt8, 1}`, which needs to be converted to
    # an `Array{UInt8, 3}`.
    screen_grab = getScreenRGB(env.ale)
    w, h = getScreenWidth(env.ale), getScreenHeight(env.ale)

    screen_grab = reshape(screen_grab, 3, w, h)

    # This step is neccessay because the image we recieve is row major
    screen_grab = permutedims(screen_grab, [3, 2, 1])
    return screen_grab
end

function get_preprocessed_greyscale(env::AtariEnv)
    screen_grab = getScreenGrayscale(env.ale) ./ 255f0
    w, h = getScreenWidth(env.ale), getScreenHeight(env.ale)

    screen_grab = reshape(screen_grab, Int(w), Int(h))

    return permutedims(screen_grab, [2, 1])
end

_get_obs(env::AtariEnv) = env.obs_type == :ram ? getRAM(env.ale) :
                            (env.obs_type == :color ? get_preprocessed_RGB(env) :
                                get_preprocessed_greyscale(env))

function reset!(env::AtariEnv)
    reset_game(env.ale)
    _get_obs(env)
end

function clone_state(env::AtariEnv)
    state_ref = cloneState(env.ale)
    state = encodeState(state_ref)
    deleteState(state_ref)
    return state
end

function restore_state(env::AtariEnv, state::Array{Int8, 1})
    state_ref = decodeState(state)
    restoreState(env.ale, state_ref)
    deleteState(state_ref)
end

ACTION_MEANING = [
    (:NOOP,),
    (:FIRE,),
    (:UP,),
    (:RIGHT,),
    (:LEFT,),
    (:DOWN,),
    (:UP, :RIGHT),
    (:UP, :LEFT),
    (:DOWN, :RIGHT),
    (:DOWN, :LEFT),
    (:UP, :FIRE),
    (:RIGHT, :FIRE),
    (:LEFT, :FIRE),
    (:DOWN, :FIRE),
    (:UP, :RIGHT, :FIRE),
    (:UP, :LEFT, :FIRE),
    (:DOWN, :RIGHT, :FIRE),
    (:DOWN, :LEFT, :FIRE)
]

get_action_meanings(env::AtariEnv) = [ACTION_MEANING[i] for i in env.action_set]

function get_keys_to_action(env::AtariEnv)
    KEYWORD_TO_KEY = Dict(
        :UP    => Int('w'),
        :DOWN  => Int('s'),
        :LEFT  => Int('a'),
        :RIGHT => Int('d'),
        :FIRE  => Int(' ')
    )

    keystoaction = Dict()

    for (action_id, action_meaning) ∈ enumerate(get_action_meanings(env))
        keys_ = [key for (keyword, key) ∈ collect(zip(keys(KEYWORD_TO_KEY), values(KEYWORD_TO_KEY)))
                    if keyword ∈ action_meaning] |> sort |> Tuple
        keystoaction[keys_] = action_id
    end
    keystoaction
end

Base.show(io::IO, env::AtariEnv) = print(io, "AtariEnv($(env.game_path), $(env.obs_type))")
export get_keys_to_action, restore_state, clone_state


text:
#
# Protein - TMAO (compare naive and linkedcells)
#

using ComplexMixtures ; const CM = ComplexMixtures
using PDBTools
using DelimitedFiles

using Plots
ENV["GKSwstype"] = "nul"

# Here we use the PDBTools package to read the pdb file (from http://github.com/m3g/PDBTools)
atoms = readPDB("../../test/data/NAMD/structure.pdb")

# The solute is a single protein molecule (infinte dilution case). In this case,
# use the option nmols=1
solute_indexes = select(atoms,"protein")
solute = CM.Selection( solute_indexes, nmols=1 )

# The solvent is TMAO, which has 14 atoms. Use the natomspermol to indicate how many
# atoms each molecule has, such that there is no ambiguity on how to split the coordinates
# of the selection into individual molecules.
solvent_indexes = select(atoms,"resname TMAO")
solvent = CM.Selection( solvent_indexes, natomspermol=14 )

# Input options for the calcualtion
options = CM.Options(binstep=0.2)

# Run MDDF calculation, and get the resutls in the R structure
trajectory = CM.Trajectory("../../test/data/NAMD/trajectory.dcd",solute,solvent)
@time R = CM.mddf(trajectory,options)

plot(layout=(6,1))

sp=1
plot!(ylabel="MDDF",subplot=sp)
plot!(R.d,R.mddf,subplot=sp,label="new - mddf")
plot!(R.d,R.rdf,subplot=sp,label="new - rdf")
plot!(legend=:topright,subplot=sp)

sp=2
plot!(ylabel="KB",subplot=sp)
plot!(R.d,R.kb,subplot=sp,label="new - mddf")
plot!(R.d,R.kb_rdf,subplot=sp,label="new - rdf")
plot!(legend=:topright,subplot=sp)

sp=3
plot!(ylabel="Count",subplot=sp)
plot!(R.d,R.md_count,subplot=sp,label="new")
plot!(R.d,R.md_count_random,subplot=sp,label="new -rand")

sp=4
plot!(ylabel="Shell vol", subplot=sp)
plot!(R.d,R.volume.shell,subplot=sp,label="new")

sp=5
plot!(ylabel="Sum MD", subplot=sp)
plot!(R.d,R.sum_md_count,subplot=sp,label="new - md")

sp=6
plot!(ylabel="Sum RAND", subplot=sp)
plot!(R.d,R.sum_md_count_random,subplot=sp,label="new - rand")
plot!(R.d,R.sum_rdf_count,subplot=sp,label="new - rdf")
plot!(legend=:topleft,subplot=sp)

plot!(size=(800,1300))
savefig("./example4.pdf")


text:
using Revise

using BNPRegressionGGA2021
using CSV
using DataStructures
using DataFrames
using Distributions
using Gadfly
using Random
using Statistics

const BNP = BNPRegressionGGA2021

#=================================#
# Example 1 - 1 trivial predictor #
#=================================#

function simulate_sample(N0, N1)
    dy = MixtureModel(LogNormal, [(0.0, 1.0), (1.0, 0.25)], [0.6, 0.4]);
    X0 = ones(N0, 1)
    X1 = ones(N1, 1)
    y0 = rand(dy, N0)
    c0 = rand(Exponential(20), N0)
    y1 = LinRange(0, 7, N1) |> collect
    z0 = deepcopy(y0)
    return dy, z0, c0, y0, X0, y1, X1
end;

Random.seed!(1);
N0, N1 = 500, 50;
dy, z0, c0, y0, X0, y1, X1 = simulate_sample(N0, N1);

m = BNP.DGPMErlang(; c0, y0, X0, y1, X1);
chainf, chainS, chainβ, chaing = BNP.sample!(m; mcmcsize = 2000, burnin = 1000);
plot(
    layer(x = y1, y = mean(chainf), Geom.line, color=["bnp"]),
    layer(x = y1, y = pdf.(dy, y1), Geom.line, color=["true"]),
)
mean(m.event)

#=====================================#
# Example 2 - 1 non-trivial predictor #
#=====================================#

function simulate_sample(N0, N1)
    dy1 = LogNormal(0, 0.5)
    dy2 = MixtureModel(LogNormal, [(0.0, 1.0), (1.0, 0.25)], [0.6, 0.4]);
    X0 = [ones(N0) rand([1, 2], N0)]
    X1 = kron(ones(N1), [1 1; 1 2])
    y0 = [X0[i, 2] == 2 ? rand(dy2) : rand(dy1) for i in 1:N0]
    y1 = LinRange(0, 6, N1) |> collect  |> x -> repeat(x, 2)
    c0 = rand(Exponential(10), N0)
    z0 = deepcopy(y0)
    return dy1, dy2, z0, c0, y0, X0, y1, X1
end

Random.seed!(1);
N0, N1 = 500, 50;
dy1, dy2, z0, c0, y0, X0, y1, X1 = simulate_sample(N0, N1);
grid = LinRange(0, 6, N1) |> collect;
plot(
    layer(x = grid, y = pdf.(dy2, grid), Geom.line, color=["kdensity2"]),
    layer(x = grid, y = pdf.(dy1, grid), Geom.line, color=["kdensity1"]),
)
plot(
    layer(x = grid, y = pdf.(dy2, grid) ./ (1.0 .- cdf.(dy2, grid)), Geom.line, color=["pdf2"]),
    layer(x = grid, y = pdf.(dy1, grid) ./ (1.0 .- cdf.(dy1, grid)), Geom.line, color=["pdf1"]),
)
m = BNP.DGPMErlang(; c0, y0, X0, y1, X1);
chainf, chainβ = BNP.sample!(m; mcmcsize = 20000, burnin = 10000);
plot(x = y1, y = mean(chainf), color = string.(X1[:, 2]), Geom.line)
mean(m.event)

#===================================================#
# Example 3 - 1 trivial and 1 non-trivial predictor #
#===================================================#

function simulate_sample(N0, N1)
    dy1 = LogNormal(0, 0.5)
    dy2 = MixtureModel(LogNormal, [(0.0, 1.0), (1.0, 0.25)], [0.6, 0.4]);
    X0 = [ones(N0) rand([1, 2], N0, 2)]
    X1 = kron(ones(N1), [1 1 1; 1 2 1; 1 1 2; 1 2 2])
    y0 = [X0[i, 2] == 2 ? rand(dy2) : rand(dy1) for i in 1:N0]
    c0 = rand(Exponential(10), N0)
    z0 = deepcopy(y0)
    y1 = LinRange(0, 6, N1) |> collect  |> x -> repeat(x, 4)
    return dy1, dy2, z0, c0, y0, X0, y1, X1
end

Random.seed!(1);
N0, N1 = 500, 50;
dy1, dy2, z0, c0, y0, X0, y1, X1 = simulate_sample(N0, N1);
ygrid = LinRange(0, 6, N1) |> collect;
plot(
    layer(x = ygrid, y = pdf.(dy2, ygrid), Geom.line, color=["pdf2"]),
    layer(x = ygrid, y = pdf.(dy1, ygrid), Geom.line, color=["pdf1"]),
)
m = BNP.DGPMErlang(; c0, y0, X0, y1, X1);
chainf, chainβ = BNP.sample!(m; mcmcsize = 20000, burnin = 10000);
plot(x = y1, y = mean(chainf), color = string.(X1[:, 2]), Geom.line)
mean(m.event)

#===================================================#
# Example 4 - 1 trivial and 4 non-trivial predictor #
#===================================================#

function simulate_sample(N0, N1)
    dy1 = LogNormal(0, 0.5)
    dy2 = MixtureModel(LogNormal, [(0.0, 1.0), (1.0, 0.25)], [0.6, 0.4]);
    X0 = [ones(N0) rand([0, 1], N0, 5)]
    X1 = kron(ones(N1), [ones(2, 5) [1, 0]])
    y0 = [X0[i, end] == 1 ? rand(dy2) : rand(dy1) for i in 1:N0]
    c0 = rand(Exponential(10), N0)
    y1 = LinRange(0, 6, N1) |> collect  |> x -> repeat(x, size(X1, 1) ÷ N1)
    z0 = deepcopy(y0)
    return dy1, dy2, z0, c0, y0, X0, y1, X1
end

Random.seed!(1);
N0, N1 = 500, 50;
dy1, dy2, z0, c0, y0, X0, y1, X1 = simulate_sample(N0, N1);

# Save the data for replication purposes
data = DataFrame(X0, :auto)
data[!, :c0] = c0
data[!, :y0] = y0
CSV.write("data/example_survival.csv", data)

# Plot the true density / survival / hazard functions
ygrid = LinRange(0, 6, N1) |> collect;
plot(
    layer(x = ygrid, y = pdf.(dy2, ygrid), Geom.line, color=["pdf2"]),
    layer(x = ygrid, y = pdf.(dy1, ygrid), Geom.line, color=["pdf1"]),
)
m = BNP.DGPMErlang(; c0, y0, X0, y1, X1);
chainf, chainβ = BNP.sample!(m; mcmcsize = 40000, burnin = 20000);
plot(x = y1, y = mean(chainf), color = string.(X1[:, end]), Geom.line)
mean(m.event)

#====================================================#
# Example 5 - 1 non-trivial and 4 trivial predictors #
#====================================================#

function simulate_sample(N0, N1)
    # dy1 = LogNormal(0, 0.3)
    # dy2 = MixtureModel(LogNormal, [(-1.0, 1.0), (0.3, 1.0)], [0.6, 0.4]);
    dy1 = LogNormal(0.3, 0.7);
    dy2 = MixtureModel(LogNormal, [(0.0, 1.0), (0.3, 1.2)], [0.6, 0.4]);
    X0 = [ones(N0) rand([0, 1], N0, 5)]
    X1 = kron([ones(2, 5) [1, 0]], ones(N1))
    y0 = [X0[i, end] == 1 ? rand(dy2) : rand(dy1) for i in 1:N0]
    c0 = rand(Exponential(10), N0)
    y1 = LinRange(0, 6, N1) |> collect |> x -> repeat(x, size(X1, 1) ÷ N1)
    z0 = deepcopy(y0)
    return dy1, dy2, z0, c0, y0, X0, y1, X1
end

Random.seed!(1);
N0, N1 = 500, 50;
dy1, dy2, z0, c0, y0, X0, y1, X1 = simulate_sample(N0, N1);
ygrid = LinRange(0, 6, N1) |> collect;

# Densities
plot(
    layer(x = ygrid, y = pdf.(dy2, ygrid), Geom.line, color=["x = 1"]),
    layer(x = ygrid, y = pdf.(dy1, ygrid), Geom.line, color=["x = 0"]),
)

# Hazard curves
plot(
    layer(x = ygrid, y = pdf.(dy2, ygrid) ./ (1.0 .- cdf.(dy2, ygrid)), Geom.line, color=["x = 1"]),
    layer(x = ygrid, y = pdf.(dy1, ygrid) ./ (1.0 .- cdf.(dy1, ygrid)), Geom.line, color=["x = 0"]),
)

m = BNP.DGPMErlang(; c0, y0, X0, y1, X1);
chainf, chainS, chainβ, chaing = BNP.sample!(m; mcmcsize = 20000, burnin = 10000);

# Show the marginal "significance" of each variable
mean(chaing)

# Generate a dataframe for plotting the posterior survival curves
chainS_mat = hcat(chainS...)
chainS_v = [chainS_mat[i, :] for i in 1:size(chainS_mat, 1)]
plot_data = DataFrame(
    z = y1,
    x5 = X1[:, end],
    St = [1 - cdf(X1[i, end] == 1 ? dy2 : dy1, y1[i]) for i in 1:size(X1, 1)],
    Sh = mean.(chainS_v),
    lb = quantile.(chainS_v, 0.025),
    ub = quantile.(chainS_v, 0.975),

)
CSV.write("data/crossing-survival-curves-fitted.csv", plot_data)

# plot(
#     plot(x = y1, y = mean.(chainS_v), color = Symbol.(X1[:, end]), Geom.line)
#     plot(x = y1, y = quantile.(chainS_v, 0.05), color = Symbol.(X1[:, end]), Geom.line)
#     plot(x = y1, y = quantile.(chainS_v, 0.95), color = Symbol.(X1[:, end]), Geom.line)
# )
# plot(x = y1, y = mean.(chainS_v), color = Symbol.(X1[:, end]), Geom.line)
# plot(x = y1, y = quantile.(chainS_v, 0.05), color = Symbol.(X1[:, end]), Geom.line)
# plot(x = y1, y = quantile.(chainS_v, 0.95), color = Symbol.(X1[:, end]), Geom.line)
# plot(x = z0, color = string.(X0[:, end]), Geom.density)
# a = mean(chaing)
# mean(m.event)
most_common_γ =
    chaing |>
    DataStructures.counter |>
    collect |>
    x -> sort(x, by = x -> x[2], rev = true)


text:
@testset "fsm_stream_reassembler_many.jl" begin
    @testset begin
        rs = StreamReassembler(1000)
        push_substring!(rs, "a", 0)
        push_substring!(rs, "ab", 0)

        @test rs.cur_idx == 2
        @test read!(rs.output) == "ab"
        @test !eof(rs.output)
    end

    @testset begin
        rs = StreamReassembler(1000)
        push_substring!(rs, "a", 0)
        @test read!(rs.output) == "a"

        push_substring!(rs, "ab", 0)

        @test rs.cur_idx == 2
        @test read!(rs.output) == "b"
    end
    @testset begin
        rs = StreamReassembler(1000)
        push_substring!(rs, "b", 1)
        @test read!(rs.output) == ""

        push_substring!(rs, "ab", 0)

        @test rs.cur_idx == 2
        @test read!(rs.output) == "ab"

        @test unassembled_bytes(rs) == 0
    end

    @testset begin
        rs = StreamReassembler(1000)
        push_substring!(rs, "b", 1)
        @test read!(rs.output) == ""

        push_substring!(rs, "bc", 1)
        @test read!(rs.output) == ""

        @test rs.cur_idx == 0
        @test rs.unassembled_bytes == 2
    end

    @testset begin
        rs = StreamReassembler(1000)
        push_substring!(rs, "c", 2)
        @test read!(rs.output) == ""

        push_substring!(rs, "bcd", 1)
        @test read!(rs.output) == ""

        @test rs.cur_idx == 0
        @test rs.unassembled_bytes == 3
    end

    @testset begin
        rs = StreamReassembler(1000)
        push_substring!(rs, "b", 1)
        push_substring!(rs, "d", 3)
        @test read!(rs.output) == ""

        push_substring!(rs, "bcde", 1)
        @test read!(rs.output) == ""

        @test rs.cur_idx == 0
        @test rs.unassembled_bytes == 4
    end

    @testset begin
        rs = StreamReassembler(1000)
        push_substring!(rs, "c", 2)
        push_substring!(rs, "bcd", 1)
        @test read!(rs.output) == ""
        @test rs.cur_idx == 0
        @test rs.unassembled_bytes == 3

        push_substring!(rs, "a", 0)
        @test read!(rs.output) == "abcd"

        @test rs.cur_idx == 4
        @test rs.unassembled_bytes == 0
    end

    @testset begin
        rs = StreamReassembler(1000)
        push_substring!(rs, "bcd", 1)
        push_substring!(rs, "c", 2)
        @test read!(rs.output) == ""
        @test rs.cur_idx == 0
        @test rs.unassembled_bytes == 3

        push_substring!(rs, "a", 0)
        @test read!(rs.output) == "abcd"

        @test rs.cur_idx == 4
        @test rs.unassembled_bytes == 0
    end

end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
# ------------------- GLOBAL ----------------------------
rdfFmt=[:key,:model,
        :unadj_avg_expsd_hh_pre,:unadj_avg_cntrl_hh_pre,:unadj_avg_expsd_hh_pst,:unadj_avg_cntrl_hh_pst,:adj_mean_expsd_grp,:adj_mean_cntrl_grp,
        :adj_dod_effct,:twotail_pval,:onetail_pval,
        :onetail_80_pct_intrvl_lb,:onetail_80_pct_intrvl_ub,:twotail_80_pct_intrvl_lb,:twotail_80_pct_intrvl_ub,
        :onetail_90_pct_intrvl_lb,:onetail_90_pct_intrvl_ub,:twotail_90_pct_intrvl_lb,:twotail_90_pct_intrvl_ub
       ]

SDF = DataFrame(key=String[],
                model=String[],
                unadj_avg_expsd_hh_pre=Float64[],
                unadj_avg_cntrl_hh_pre=Float64[],
                unadj_avg_expsd_hh_pst=Float64[],
                unadj_avg_cntrl_hh_pst=Float64[],
                adj_mean_expsd_grp=Float64[],
                adj_mean_cntrl_grp=Float64[],
                adj_dod_effct=Float64[],
                twotail_pval=Float64[],
                onetail_pval=Float64[],
                onetail_80_pct_intrvl_lb=Float64[],
                onetail_80_pct_intrvl_ub=Float64[],
                twotail_80_pct_intrvl_lb=Float64[],
                twotail_80_pct_intrvl_ub=Float64[],
                onetail_90_pct_intrvl_lb=Float64[],
                onetail_90_pct_intrvl_ub=Float64[],
                twotail_90_pct_intrvl_lb=Float64[],
                twotail_90_pct_intrvl_ub=Float64[],
                mean_score0=Float64[],
                mean_score1=Float64[],
                B=Float64[],
                SE=Float64[],
                P=Float64[],
                M=Float64[],
                Mt=Float64[],
                Mc=Float64[]
              )

function pushSDFrow!(sdf::DataFrame, key::AbstractString, v_model::AbstractString)
    push!(sdf, [key; v_model; 0.0; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA])
    idx=length(sdf[1])
    sdf[idx,:unadj_avg_expsd_hh_pre] = NA
    return sdf
end

function pushSDFrow!(sdf::DataFrame, iDict::OrderedDict)
    key=get(iDict, :key, NA)
    model=get(iDict, :model, NA)
    if isna(key) | isna(model)
        error("Error: pushSDFrow! - missing key or model values")
    else
        B=get(iDict, :B, NA)
        SE=get(iDict, :SE, NA)
        P=get(iDict, :P, NA)
        M=get(iDict, :M, NA)
        Mt=get(iDict, :Mt, NA)
        Mc=get(iDict, :Mc, NA)
        push!(sdf, [key; model; 0.0; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; B; SE; P; M; Mt; Mc])

        #push!(sdf, [key; v_model; 0.0; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA])
        idx=length(sdf[1])
        sdf[idx,:unadj_avg_expsd_hh_pre] = NA
    end
    return sdf
end

#this.sdf=pushSDFrow!(deepcopy(SDF),"Total Campaign","occ")

#push!(SDF, ["Total Campaign"; ""; 0.0; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA; NA])
#SDF[:unadj_avg_expsd_hh_pre]=NA
#SDF[:model]=NA

function df2dict(idf::DataFrame)
    dout=OrderedDict()
    for n in names(idf)
       dout[n] = idf[1,Symbol(n)]
    end
    return dout
end

function df2dict(row::DataFrameRow)
    dout=OrderedDict()
    for n in names(row)
       dout[n] = row[Symbol(n)]
    end
    return dout
end

function write2disk(idf::DataFrame,fname::AbstractString, repvalues::Bool=true)
    ofile=fname
    if repvalues
        writetable(ofile*"tmp", idf, separator = ',', header = true)
        outfile = open(ofile, "w")
        open(ofile*"tmp") do filehandle
            for line in eachline(filehandle)
                write(outfile, replace(line,"NA",""))
            end
        end
        close(outfile)
        rm(ofile*"tmp")
    else
        writetable(ofile, idf, separator = ',', header = true)
    end
    println("Complete !!!! - output written to ",ofile)
end

function Base.lowercase(df::DataFrame)
    for (i,v) in enumerate(names(df))   if v != Symbol(lowercase(string(v)))  rename!(df,v,Symbol(lowercase(string(v)))) end end
end

"""
app.cfg
P2_Competitor = true
pvalue_lvl = 0.20
excludedKeys = key1,key2,key3
exposed_flag_var = new_exposed_flag
sigLevel = 0.2
random_demos =
random_campaigns =
dropvars = goop,snop
scoring_vars =
TotalModelsOnly=false
"""

function loadCFG(d::OrderedDict,fn::AbstractString)
    function cvt(v::AbstractString, t::DataType)
        if t in [Bool]
            return lowercase(v) =="true" ? true : lowercase(v) =="false" ? false : NA
        elseif t in [nothing]
            return NA
        elseif t in [Array{AbstractString,1}]
            #return AbstractString[string(x) split(v,",")]
            return AbstractString[string(x) for x in split(replace(v,"\"",""),",")]
        elseif t in [Array{Symbol,1}]
            return Symbol[Symbol(x) for x in split(replace(v,"\"",""),",")]
        elseif t in [Symbol]
            return Symbol(v)
        elseif t in [String]
            return string(v)
        elseif t in [Float64]
            return parse(Float64,v)
        else
            return "__ooops__"
        end
    end
    if isfile(fn)
        f = open(fn)
        for ln in eachline(f)
            ln=replace(ln,"\n","")
            sidx=search(ln,"=").start-1
            eidx=search(ln,"=").stop+1
            if sidx>0
                k=strip(ln[1:sidx])
                v=strip(ln[eidx:end])
                if (length(k)>0) & (length(v)>0)
                    #println("LINE:: ~"*k*"~"*v*"~")
                    #println("DEF::       ______ ",d[Symbol(k)],"  ----> ",typeof(d[Symbol(k)]) )
                    newk=Symbol(k)
                    newv = cvt(v,typeof(d[Symbol(k)]))
                    #println("IN::        ______ ",Symbol(k),"  ----> ",ct, " ....... ", typeof(ct))
                    #println( d[Symbol(k)],"|",v," {",typeof(d[Symbol(k)]),"}  ----->  ",newk,"|",newv,"  {",typeof(newv),"}")
                    d[Symbol(k)] = newv
                else
                    #println("CFG ERROR :: invalid length : "*"~"*k*"~"*v*"~",typeof(d[Symbol(k)]))
                end
            else
                if length(ln)>0   println("CFG ERROR :: invalid cfg parameter : ~",ln,"~") end
            end
        end
        close(f)
    else
        println("Using Default configuration")
    end
    return d
end

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
using Base.Test
using DataFrames
import WCS

const rcf = RunCamcolField(3900, 6, 269)

"""
Get the PSF located at a particular world location in an image.

Args:
 - world_loc: A location in world coordinates.
 - img: An TiledImage

Returns:
 - An array of PsfComponent objects that represents the PSF as a mixture
     of Gaussians.
"""
function get_source_psf(world_loc::Vector{Float64}, img::TiledImage)
    # Some stamps or simulated data have no raw psf information.    In that case,
    # just use the psf from the image.
    if size(img.raw_psf_comp.rrows) == (0, 0)
        # Also return a vector of empty psf params
        return img.psf, fill(fill(NaN, length(PsfParams)), psf_K)
    else
        pixel_loc = WCS.world_to_pix(img.wcs, world_loc)
        psfstamp = img.raw_psf_comp(pixel_loc[1], pixel_loc[2])
        return PSF.fit_raw_psf_for_celeste(psfstamp)
    end
end

function test_blob()
    # A lot of tests are in a single function to avoid having to reload
    # the full image multiple times.
    blob = SDSSIO.load_field_images(rcf, datadir)

    dir = "$datadir/$(rcf.run)/$(rcf.camcol)/$(rcf.field)"
    fname = @sprintf "%s/photoObj-%06d-%d-%04d.fits" dir rcf.run rcf.camcol rcf.field
    cat_entries = SDSSIO.read_photoobj_celeste(fname)

    ea = make_elbo_args(blob, cat_entries,
                       patch_radius_pix=1e-6, fit_psf=false, tile_width=20)

    # Just check some basic facts about the catalog.
    @test length(cat_entries) == 805
    @test 0 < sum([ce.is_star for ce in cat_entries]) < 805

    # Find an object near the middle of the image.
    ctr = WCS.pix_to_world(blob[3].wcs, [blob[3].H / 2, blob[3].W / 2])
    dist = [sqrt((ce.pos[1] - ctr[1])^2 + (ce.pos[2] - ctr[2])^2)
                    for ce in cat_entries]
    obj_index = findmin(dist)[2]    # index of closest object
    obj_loc = cat_entries[obj_index].pos    # location of closest object

    # Test get_source_psf at point while we have the blob loaded.
    test_b = 3
    img = ea.images[test_b]
    ea_obj = make_elbo_args(blob, cat_entries[obj_index:obj_index])
    pixel_loc = WCS.world_to_pix(img.wcs, obj_loc);
    original_psf_val = img.raw_psf_comp(pixel_loc[1], pixel_loc[2])

    original_psf_celeste = PSF.fit_raw_psf_for_celeste(original_psf_val)[1];
    fit_original_psf_val = PSF.get_psf_at_point(original_psf_celeste);

    obj_psf = get_source_psf(ea_obj.vp[1][ids.u], img)[1];
    obj_psf_val = PSF.get_psf_at_point(obj_psf);

    # The fits should match exactly.
    @test_approx_eq_eps(obj_psf_val, fit_original_psf_val, 1e-6)

    # The raw psf will not be as good.
    @test_approx_eq_eps(obj_psf_val, original_psf_val, 1e-2)

    cat_several = [cat_entries[1]; cat_entries[obj_index]]
    ea_several = make_elbo_args(ea_obj.images, cat_several)

    # The second set of vp is the object of interest
    point_patch_psf = PSF.get_psf_at_point(ea_several.patches[2, test_b].psf);
    # The threshold for the test below was formerly 1e-6.
    # Is it a problem I needed to increase it?
    @test_approx_eq_eps(obj_psf_val, point_patch_psf, 1e-4)
end

function test_stamp_get_object_psf()
    stamp_blob, stamp_mp, body = gen_sample_star_dataset();
    img = TiledImage(stamp_blob[3]);
    obj_index =    stamp_mp.vp[1][ids.u]
    pixel_loc = WCS.world_to_pix(img.wcs, obj_index)
    original_psf_val = PSF.get_psf_at_point(img.psf);

    obj_psf_val =
        PSF.get_psf_at_point(get_source_psf(stamp_mp.vp[1][ids.u], img)[1])
    @test_approx_eq_eps(obj_psf_val, original_psf_val, 1e-6)
end

function test_get_tiled_image_source()
    # Test that an object only occurs the appropriate tile's local sources.
    blob, ea, body = gen_sample_star_dataset();

    ea = make_elbo_args(blob, body; patch_radius_pix=1e-6)

    tiled_img = TiledImage(blob[3], tile_width=10);
    for hh in 1:size(tiled_img.tiles, 1), ww in 1:size(tiled_img.tiles, 2)
        tile = tiled_img.tiles[hh, ww]
        loc = Float64[mean(tile.h_range), mean(tile.w_range)]
        for b = 1:5
            ea.vp[1][ids.u] = loc
            pixel_center = WCS.world_to_pix(blob[b].wcs, loc)
            wcs_jacobian = Model.pixel_world_jacobian(blob[b].wcs, pixel_center)
            radius_pix = maxabs(eigvals(wcs_jacobian)) * 1e-6
            ea.patches[1, b] = SkyPatch(loc,
                                        radius_pix,
                                        blob[b].psf,
                                        wcs_jacobian,
                                        pixel_center)
        end
        patches = vec(ea.patches[:, 3])
        local_sources = Model.get_sources_per_tile(tiled_img.tiles,
                                                   Model.patch_ctrs_pix(patches),
                                                   Model.patch_radii_pix(patches))
        @test local_sources[hh, ww] == Int[1]
        for hh2 in 1:size(tiled_img.tiles, 1), ww2 in 1:size(tiled_img.tiles, 2)
            if (hh2 != hh) || (ww2 != ww)
                @test local_sources[hh2, ww2] == Int[]
            end
        end
    end
end

function test_local_source_candidate()
    blob, ea, body = gen_n_body_dataset(100);

    for b=1:ea.N
        # Get the sources by iterating over everything.
        patches = vec(ea.patches[:,b])

        tile_source_map = Model.get_sources_per_tile(ea.images[b].tiles,
                        Model.patch_ctrs_pix(patches),
                        Model.patch_radii_pix(patches))

        # Check that all the actual sources are candidates and that this is the
        # same as what is returned by make_elbo_args
        HH, WW = size(tile_source_map)
        for h=1:HH, w=1:WW
            # Get a set of candidates.
            candidates = Model.local_source_candidates(
                                                ea.images[b].tiles[h, w],
                                                Model.patch_ctrs_pix(patches),
                                                Model.patch_radii_pix(patches))
            @test setdiff(tile_source_map[h, w], candidates) == []
            @test tile_source_map[h, w] == ea.tile_source_map[b][h, w]
        end
    end
end

function test_set_patch_size()
    # Test that the patch size gets most of the light from a variety of
    # galaxy shapes.
    # This shows that the current patch size is actually far too conservative.

    function gal_catalog_from_scale(gal_scale::Float64, flux_scale::Float64)
        CatalogEntry[CatalogEntry(world_location, false,
                                   flux_scale * fluxes, flux_scale * fluxes,
                                   0.1, .01, pi/4, gal_scale, "sample", 0) ]
    end

    srand(1)
    blob0 = SampleData.load_stamp_blob(datadir, "164.4311-39.0359_2kpsf");
    img_size = 150
    for b in 1:5
            blob0[b].H, blob0[b].W = img_size, img_size
    end
    fluxes = [4.451805E+03,1.491065E+03,2.264545E+03,2.027004E+03,1.846822E+04]

    world_location = WCS.pix_to_world(blob0[3].wcs,
                                      Float64[img_size / 2, img_size / 2])

    for gal_scale in [1.0, 10.0], flux_scale in [0.1, 10.0]
        cat = gal_catalog_from_scale(gal_scale, flux_scale);
        blob = Synthetic.gen_blob(blob0, cat);
        ea = make_elbo_args(blob, cat, tile_width=typemax(Int));

        for b=1:length(blob)
            @assert size(ea.images[b].tiles) == (1, 1)
            tile_image = DeterministicVI.tile_predicted_image(
                ea.images[b].tiles[1,1], ea, ea.tile_source_map[b][1,1]);

            pixel_center = WCS.world_to_pix(blob[b].wcs, cat[1].pos)
            radius = Model.choose_patch_radius(
                pixel_center, cat[1], blob[b].psf, ea.images[b])

            circle_pts = fill(false, blob[b].H, blob[b].W);
            in_circle = 0.0
            for x=1:size(tile_image)[1], y=1:size(tile_image)[2]
                if ((x - pixel_center[1]) ^ 2 + (y - pixel_center[2]) ^ 2) < radius ^ 2
                    in_circle += tile_image[x, y]
                    circle_pts[x, y] = true
                end
            end
            @test in_circle / sum(tile_image) > 0.95

            # Convenient for visualizing:
            # using PyPlot
            # in_circle / sum(tile_image)
            # imshow(tile_image); colorbar()
            # imshow(circle_pts, alpha=0.4)
        end
    end
end

test_blob()
test_stamp_get_object_psf()
test_get_tiled_image_source()
test_local_source_candidate()
test_set_patch_size()


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    PointProcess

A spatial point process.
"""
abstract type PointProcess end

"""
    ishomogeneous(p)

Tells whether or not the spatial point process `p` is homogeneous.
"""
ishomogeneous(p::PointProcess) = false

"""
    rand(p, g, n=1; [algo])

Generate `n` realizations of spatial point process `p`
inside geometry `g`. Optionally specify sampling
algorithm `algo`.
"""
Base.rand(p::PointProcess, g::Geometry, n::Int;
          algo=default_sampling_algorithm(p)) =
  [rand_single(p, g, algo) for i in 1:n]

Base.rand(p::PointProcess, g::Geometry;
          algo=default_sampling_algorithm(p)) =
  rand_single(p, g, algo)

"""
    rand_single(p, g, algo)

Generate a single realization of spatial point process
`p` inside geometry `g` with sampling `algo`.
"""
rand_single(::PointProcess, ::Geometry, algo) =
  @error "not implemented"

"""
    default_sampling_algorithm(p)

Default sampling algorithm for spatial point process `p`.
"""
default_sampling_algorithm(::PointProcess) = @error "not implemented"

"""
    p₁ ∪ p₂

Union (or superposition) of spatial point processes `p₁` and `p₂`.
"""
Base.union(p₁::PointProcess, p₂::PointProcess) = UnionProcess(p₁, p₂)

#-----------------
# IMPLEMENTATIONS
#-----------------
include("processes/binomial.jl")
include("processes/poisson.jl")
include("processes/union.jl")


text:
# ucnv.jl - Wrapper for ICU (International Components for Unicode) library

# Some content of the documentation strings was derived from the ICU header files ucnv.h
# (Those portions copyright (C) 1996-2015, International Business Machines Corporation and others)

"""
"""
module ucnv end

macro libcnv(s)     ; _libicu(s, iculibi18n, "ucnv_")     ; end

export UConverter

mutable struct UConverter
    p::Ptr{Cvoid}

    function UConverter(name::ASCIIStr)
        err = Ref{UErrorCode}(0)
        p = ccall(@libcnv(open), Ptr{Cvoid}, (Cstring, Ptr{UErrorCode}), name, err)
        SUCCESS(err[]) || error("ICU: $(err[]), could not open converter ", name)
        self = new(p)
        finalizer(self, close)
        self
    end
end

close(c::UConverter) =
    c.p == C_NULL || (ccall(@libcnv(close), Cvoid, (Ptr{Cvoid},), c.p); c.p = C_NULL)

struct UConverterPivot
    buf::Vector{UCS2Chr}
    pos::Vector{Ptr{UCS2Chr}}

    function UConverterPivot(n::Int)
        buf = create_vector(UCS2Chr, n)
        p = pointer(buf)
        new(buf, [p,p])
    end
end

function convert!(dstcnv::UConverter, srccnv::UConverter,
                  dst::IOBuffer, src::IOBuffer, pivot::UConverterPivot,
                  reset::Bool=false, flush::Bool=true)
    @preserve dst src pivot begin
        p = Ptr{UInt8}[pointer(dst.data, position(dst)+1),
                       pointer(src.data, position(src)+1)]
        p0 = copy(p)
        err = Ref{UErrorCode}(0)
        ccall(@libcnv(convertEx), Cvoid,
              (Ptr{Cvoid}, Ptr{Cvoid},
               Ptr{Ptr{UInt8}}, Ptr{UInt8}, Ptr{Ptr{UInt8}}, Ptr{UInt8},
               Ptr{UCS2Chr}, Ptr{Ptr{UCS2Chr}}, Ptr{Ptr{UCS2Chr}}, Ptr{UCS2Chr},
               UBool, UBool, Ptr{UErrorCode}),
              dstcnv.p, srccnv.p,
              pointer(p, 1), pointer(dst.data, length(dst.data)+1),
              pointer(p, 2), pointer(src.data, src.size+1),
              pointer(pivot.buf, 1),
              pointer(pivot.pos, 1),
              pointer(pivot.pos, 2),
              pointer(pivot.buf, length(pivot.buf)+1),
              reset, flush, err)
        dst.size += p[1] - p0[1]
        dst.ptr += p[1] - p0[1]
        src.ptr += p[2] - p0[2]
    end
    err[] == U_BUFFER_OVERFLOW_ERROR && return true
    @assert SUCCESS(err[])
    false
end

function to_uchars(cnv::UConverter, src::Vector{UInt8})
    @preserve src begin
        srclen = length(src)
        dstlen = 2 * srclen + 1
        buf, pnt = _allocate(UInt16, dstlen)
        err = Ref{UErrorCode}(0)
        n = ccall(@libcnv(toUChars), Int32,
                  (Ptr{Cvoid}, Ptr{UCS2Chr}, Int32, Ptr{UInt8}, Int32, Ptr{UErrorCode}),
                  cnv.p, pnt, dstlen, src, srclen, err)
        SUCCESS(err[]) || error("ICU: $(err[]), could not convert string")
        Str(UTF16CSE, u[1:n])
    end
end

"""
    Determines if the converter contains ambiguous mappings of the same character or not.

    Arguments:
    cnv - the converter to be tested

    Return `true` if the converter contains ambiguous mapping of the same character, `false` otherwise.
"""
function isambiguous(cnv::UConverter)
    err = Ref{UErrorCode}(0)
    v = ccall(@libcnv(isAmbiguous), Bool, (Ptr{Cvoid}, Ptr{UErrorCode}), cnv.p, err)
    SUCCESS(err[]) || error("ICU: internal error $(err[]) in ucnv_isFixedWidth")
    v
end

"""
    Detects Unicode signature byte sequences at the start of the byte stream
    and returns the charset nameof the indicated Unicode charset.
    NULL is returned when no Unicode signature is recognized.
    The number of bytes in the signature is output as well.

    The caller can ucnv_open() a converter using the charset name.
    The first code unit (UCS2Chr) from the start of the stream will be U+FEFF
    (the Unicode BOM/signature character) and can usually be ignored.

    For most Unicode charsets it is also possible to ignore the indicated
    number of initial stream bytes and start converting after them.
    However, there are stateful Unicode charsets (UTF-7 and BOCU-1) for which
    this will not work. Therefore, it is best to ignore the first output UCS2Chr
    instead of the input signature bytes.

    Arguments:
    source - The source string in which the signature should be detected.

    Returns: the name of the encoding detected, and the length of the signature
"""
function detect_unicode_signature(src::Vector{UInt8})
    err = Ref{UErrorCode}(0)
    sig = Ref{Int32}(0)
    p = ccall(@libcnv(detectUnicodeSignature), Ptr{UInt8},
              (Ptr{UInt8}, Int32, Ptr{Int32}, Ptr{UErrorCode}),
              src, sizeof(src), sig, err)
    SUCCESS(err[]) || error("ICU: internal error $(err[]) in ucnv_detectUnicodeSignature")
    return (p == C_NULL ? UTF8Str() : UTF8Str(p), sig[])
end

"""
    Returns the number of chars held in the converter's internal state
    because more input is needed for completing the conversion. This function is
    useful for mapping semantics of ICU's converter interface to those of iconv,
    and this information is not needed for normal conversion.

    Arguments:
    cnv - The converter in which the input is held as internal state

    Returns: the number of chars in the state or -1 if an error is encountered.
"""
function from_ucount_pending(cnv::UConverter)
    err = Ref{UErrorCode}(0)
    v = ccall(@libcnv(toUCountPending), Int32, (Ptr{Cvoid}, Ptr{UErrorCode}), cnv.p, err)
    SUCCESS(err[]) || error("ICU: internal error $(err[]) in ucnv_fromUCountPending")
    v
end

"""
    Returns the number of chars held in the converter's internal state
    because more input is needed for completing the conversion. This function is
    useful for mapping semantics of ICU's converter interface to those of iconv,
    and this information is not needed for normal conversion.

    Arguments:
    cnv - The converter in which the input is held as internal state

    Returns: the number of chars in the state or -1 if an error is encountered.
"""
function to_ucount_pending(cnv::UConverter)
    err = Ref{UErrorCode}(0)
    v = ccall(@libcnv(toUCountPending), Int32, (Ptr{Cvoid}, Ptr{UErrorCode}), cnv.p, err)
    SUCCESS(err[]) || error("ICU: internal error $(err[]) in ucnv_toUCountPending")
    v
end

"""
    Returns whether or not the charset of the converter has a fixed number of bytes
    per charset character.
    An example of this are converters that are of the type UCNV_SBCS or UCNV_DBCS.
    Another example is UTF-32 which is always 4 bytes per character.
    A Unicode code point may be represented by more than one UTF-8 or UTF-16 code unit
    but a UTF-32 converter encodes each code point with 4 bytes.
    Note: This method is not intended to be used to determine whether the charset has a
    fixed ratio of bytes to Unicode codes <i>units</i> for any particular Unicode encoding form.
    `false` is returned with the UErrorCode if error occurs or cnv is NULL.

    Arguments:
    cnv       The converter to be tested
    status    ICU error code in/out paramter

    Returns: `true` if the converter is fixed-width
"""
function is_fixed_width(cnv::UConverter)
    err = Ref{UErrorCode}(0)
    v = ccall(@libcnv(isFixedWidth), Bool, (Ptr{Cvoid}, Ptr{UErrorCode}), cnv.p, err)
    SUCCESS(err[]) || error("ICU: internal error $(err[]) in ucnv_isFixedWidth")
    v
end


text:
abstract type AbstractServiceFormulation end

abstract type AbstractReservesForm <: AbstractServiceFormulation end

abstract type AbstractRegulationReserveForm <: AbstractReservesForm end

struct RampLimitedReserve <: AbstractReservesForm end

struct LoadProportionalReserve <: AbstractReservesForm end


text:
[
     "T4" => "T4 - en_Latn",
     "T5" => "T5 - en_Latn",
     "T6" => "T6 - en_Latn",
]


text:
using Test, NES

@testset "Build Test" begin
    env = SMBEnv("smb", :RIGHT_ONLY)
    @test true
end


text:
"""
```
initialize!(m::Li2020)
```

sets up all initial conditions for solving Li2020, such as the grid and boundary conditions.
"""
function initialize!(m::Li2020)

    # Use no jump solution for p as initial guess
    stategrid, funcvar, derivs, endo = solve(m; nojump = true, nojump_method = :ode) # This calls initialize_nojump! already

    # Interpolate w/SLM, first over the "ODE" part, then interpolate once more
    p₀, p₁ = get_setting(m, :boundary_conditions)[:p]
    ψ_is_1 = findfirst(funcvar[:p] .>= p₁)
    p_SLM = SLM(stategrid[:w][1:ψ_is_1],  funcvar[:p][1:ψ_is_1];  concave_down = true, left_value = p₀,
                right_value = p₁, increasing = true, knots = floor(Int, ψ_is_1 / 4))
    funcvar[:p][1:ψ_is_1] .= eval(p_SLM, stategrid[:w][1:ψ_is_1])
    p_SLM = SLM(stategrid[:w], funcvar[:p], concave_down = true, increasing = true, left_value = p₀,
                right_value = p₁, knots = floor(Int, ψ_is_1 / 8))
    funcvar[:p] = eval(p_SLM, stategrid[:w])

    # Guess for xg and Q̂
    Q = get_setting(m, :gov_bond_gdp_level) * get_setting(m, :avg_gdp)
    funcvar[:xg]     .= Q ./ (2. .* stategrid[:w] .* funcvar[:p])
    funcvar[:xg][1]   = 2. * funcvar[:xg][2]
    funcvar[:xg][end] = m[:ϵ]
    funcvar[:Q̂]      .= fill(.05 * Q, length(stategrid))

    # Some settings for functional iteration
    m <= Setting(:κp_grid, vcat(0., exp.(range(log(1e-3), stop = log((p₁ - p₀) / p₁), length = 19))))

    return stategrid, funcvar, derivs, endo
end

"""
```
initialize_nojump!(m::Li2020)
```

initializes the no-jump equilibrium for the Li (2020) model.
"""
function initialize_nojump!(m::Li2020)
    model_type = eltype(m)
    N          = get_setting(m, :N)

    # Create StateGrid object
    stategrid_init = OrderedDict{Symbol, Vector{model_type}}()
    gen_grid(l, u, n) = exp.(range(log(l), stop = log(u), length = n))
    stategrid_init[:w] = vcat(0., gen_grid(get_setting(m, :stategrid_dimensions)[:w][1],
                                           get_setting(m, :stategrid_splice), Int(round(N / 2))),
                              gen_grid(get_setting(m, :stategrid_splice) + .01, get_setting(m, :stategrid_dimensions)[:w][2],
                              Int(N - round(N / 2) - 1)))
    stategrid = StateGrid(stategrid_init)

    # Construct dictionary of functional variables
    funcvar = OrderedDict{Symbol, Vector{model_type}}()
    for k in keys(get_functional_variables(m))
        funcvar[k] = Vector{model_type}(undef, N)
    end

    # Construct dictionary of derivatives
    derivs = OrderedDict{Symbol, Vector{model_type}}()
    for k in keys(get_derivatives(m))
        derivs[k] = Vector{model_type}(undef, N)
    end

    # Construct dictionary of endogenous variables
    endo = OrderedDict{Symbol, Vector{model_type}}()
    for k in keys(get_endogenous_variables(m))
        endo[k] = Vector{model_type}(undef, N)
    end

    # Establish boundary conditions
    p₀ = find_zero(p -> get_setting(m, :Φ)(p, m[:χ].value, m[:δ].value) + m[:ρ] * p - m[:AL], 1.) # bisection search for p(0) and p(1)
    p₁ = find_zero(p -> get_setting(m, :Φ)(p, m[:χ].value, m[:δ].value) + m[:ρ] * p - m[:AH], 1.)
    p₀ *= 1 + get_setting(m, :p₀_perturb)
    ∂p∂w0 = (m[:AH] - m[:AL]) / (get_setting(m, :∂Φ)(p₀, m[:χ].value) + m[:ρ]) *
                              ((m[:AH] - m[:AL]) / (p₀ * m[:σK])^2 * p₀ + 1)
    ∂p∂wN = 0.
    set_boundary_conditions!(m, :p, [p₀, p₁])
    set_boundary_conditions!(m, :∂p_∂w, [∂p∂w0, ∂p∂wN])

    return stategrid, funcvar, derivs, endo
end


text:
@testset "textprocessing.jl" begin
    text = "So dim, so dark, So dense, so dull, So damp, so dank, so dead! The weather, now warm, now cold, Makes it harder Than ever to forget!"
    c = countwords(text)
    @test c["so"] == 4
    @test c["Makes"] == 1
    WordCloud.TextProcessing.casemerge!(c)
    @test "so" in keys(c)
    @test !("So" in keys(c)) # casemerge!
    words, weights = WordCloud.TextProcessing.processtext(c)
    @test !("to" in words) # stopwords
    @test lemmatize("Cars") == "Car"
    @test lemmatize("monkeys") == "monkey"
    @test lemmatize("politics") == "politics"
    @test lemmatize("less") == "less"
    @test lemmatize("makes") == "make"
    @test lemmatize("does") == "do"
    @test lemmatize("shoes") == "shoe"
    @test lemmatize("buses") == "bus"
    @test lemmatize("classes") == "class"
    @test lemmatize("closes") == "close"
    @test lemmatize("licenses") == "license"
    @test lemmatize("babies") == "baby"
    @test lemmatize("studies") == "study"
    @test lemmatize("dies") == "die"
    @test lemmatize("series") == "series"
    @test lemmatize("wolves") == "wolf"
    @test lemmatize("knives") == "knife"
    @test lemmatize("ourselves") == "ourselves"
    @test lemmatize("loves") == "love"
    @test lemmatize("lives") in ("life", "live")
    @test lemmatize("cos") == "cos"
    @test lemmatize("中文") == "中文"

    @test lemmatize!(Dict("dog" => 1, "dogs" => 2)) == Dict("dog" => 3)
    @test lemmatize!(Dict("cat" => 1, "dogs" => 2)) == Dict("dog" => 2, "cat" => 1)
    @test length(processtext(["cat" => 1, "dog" => 1, "dogs" => 3, "Dogs" => 2, "Dog" => 1])[1]) == 2
    @test processtext(["cat" => 3, "Dog" => 1, "dogs" => 2])[2] |> diff |> only |> iszero
    @test processtext("word cloud") == processtext(["word","cloud"], [12,12]) == processtext([("word", 3), ("cloud", 3)])
end


text:
function getTableauHPARK(name, q::Tableau{T}, p::Tableau{T}, d=[]) where {T}

    @assert q.s == p.s

    o = min(q.o, p.o)

    a_q = q.a
    a_p = p.a
    b_q = q.b
    b_p = p.b
    c_q = q.c
    c_p = p.c

    if length(d) == 0
        return TableauHPARK(name, o,
                            a_q, a_p, a_q, a_p,
                            a_q, a_p, a_q, a_p,
                            b_q, b_p, b_q, b_p,
                            c_q, c_p, c_q, c_p)
    else
        @assert length(d) == q.s == p.s

        return TableauHPARK(name, o,
                            a_q, a_p, a_q, a_p,
                            a_q, a_p, a_q, a_p,
                            b_q, b_p, b_q, b_p,
                            c_q, c_p, c_q, c_p,
                            d)
    end
end

"SPARK tableau for Gauss-Lobatto IIIA-IIIB HPARK method with s stages."
function TableauHPARKLobattoIIIAIIIB(s)
    getTableauHPARK(Symbol("HPARKLobattoIIIAIIIB($s)"), TableauLobattoIIIA(s), TableauLobattoIIIB(s))
end

"SPARK tableau for Gauss-Lobatto IIIB-IIIA  method with s stages."
function TableauHPARKLobattoIIIBIIIA(s)
    getTableauHPARK(Symbol("HPARKLobattoIIIBIIIA($s)"), TableauLobattoIIIB(s), TableauLobattoIIIA(s))
end

"Tableau for Gauss-Legendre HPARK method with s stages."
function TableauHPARKGLRK(s)
    glrk = TableauGauss(s)
    getTableauHPARK(Symbol("HPARKGLRK", s), glrk, glrk)
end


text:
using ValueInferenceTools
using ValueInferenceTools: to_sign
using Test
using Random

@testset "inferenceaccumulator" begin
    prob = TaskStats(0.9, 0.3)
    rews = [true, true, false, false, true]
    side = [true, true ,true, false, false]

    ratios = zeros(length(rews))
    for i in 2:length(rews)
        term = (ratios[i-1]+prob.dpl)/(1-prob.dpl)
        side[i] == side[i-1] || (term = 1/term)
        ratios[i] = rews[i] ? 0.0 : term/(1-prob.rwd)
    end
    vals = @. $(prob.rwd) * (1 - 2 / (ratios ^ (-to_sign(side)) + 1))

    @test vals ≈ accumulatefromdata(InferenceAccumulator(prob), rews, side)

    Random.seed!(4)
    rews = rand(Bool, 100)
    side = rand(Bool, 100)

    ratios = zeros(length(rews))
    for i in 2:length(rews)
        term = (ratios[i-1]+prob.dpl)/(1-prob.dpl)
        side[i] == side[i-1] || (term = 1/term)
        ratios[i] = rews[i] ? 0.0 : term/(1-prob.rwd)
    end
    vals = @. $(prob.rwd) * (1 - 2 / (ratios ^ (-to_sign(side)) + 1))
    @test vals ≈ accumulatefromdata(InferenceAccumulator(prob), rews, side)
end


text:
@testset "convert" begin
    @test convert(NanoDate, daytime) == nd0
    @test convert(NanoDate, dayt) == nd00
end

@testset "promote" begin
    @test promote(nd, daytime) == (nd, NanoDate(daytime))
    @test promote(nd, dayt) == (nd, NanoDate(dayt))
end

@testset "splice" begin
    @test NanoDate(nd, latertime) == NanoDate(Date(nd), latertime)
    @test NanoDate(nd, laterdate) == NanoDate(laterdate, Time(nd))
    @test NanoDate(nd, Nanosecond(444)) ==
        NanoDate(nd.datetime, nanosecs(Microsecond(nd), Nanosecond(444)))
    @test NanoDate(nd, Microsecond(888)) ==
        NanoDate(nd.datetime, nanosecs(Microsecond(888),Nanosecond(nd)))
    @test NanoDate(nd, Microsecond(888), Nanosecond(444)) ==
        NanoDate(nd.datetime, nanosecs(Microsecond(888),Nanosecond(444)))
end

@testset "substitute" begin
    @test NanoDate(nd, Year(nd) - Year(2)) + Year(2) == nd
    @test NanoDate(nd, Quarter(nd) - Quarter(2)) + Quarter(2) == nd
    @test NanoDate(nd, Month(nd) - Month(2)) + Month(2) == nd
    @test NanoDate(nd, Week(nd) - Week(2)) + Week(2) == nd
    @test nd - NanoDate(nd, Dy+Day(2)) == -Day(2)

    @test nd - NanoDate(nd, Hr-Hour(2)) == Hour(2)
    @test nd - NanoDate(nd, Mi+Minute(122)) == -Minute(122)
    @test nd - NanoDate(nd, Sc-Second(64)) == Second(64)
    @test nd - NanoDate(nd, Ms-Millisecond(324)) == Millisecond(324)

    @test nd - NanoDate(nd, Cs-Microsecond(800)) == Microsecond(800)
    @test nd - NanoDate(nd, Ns-Nanosecond(999)) == Nanosecond(999)
end


text:
using Dates
using Printf

struct Clock

    hr::Integer
    min::Integer

    function Clock(hr::Integer, min::Integer)

        # Add/subtract extra minutes
        if min % 60 == 0
            hr += div(min, 60)
            min = 0
        elseif min < 0
            hr += div(min, 60) - 1
            min = min % 60 + 60
        else
            hr += div(min, 60)
            min = min % 60
        end

        if hr % 24 == 0
            hr = 0
        elseif hr < 0
            hr = hr % 24 + 24
        else
            hr = hr % 24
        end

        new(hr, min)
    end
end

function Base.:+(clock::Clock, time::Dates.Minute)
    delta_mins = time.value
    return Clock(clock.hr, clock.min + delta_mins)
end

function Base.:-(clock::Clock, time::Dates.Minute)
    delta_mins = time.value
    return Clock(clock.hr, clock.min - delta_mins)
end

function Base.show(io::IO, clock::Clock)
    @printf(io, "%2d:%2d", clock.hr, clock.min)
end

max_stars_repo_name: jfcliang/jl-exercism-solutions

max_issues_repo_name: jfcliang/jl-exercism-solutions

max_forks_repo_name: jfcliang/jl-exercism-solutions


text:
# ToDo: Do crazy stuff for build here
# Use e.g. Binary Dependency Builder for Julia

max_stars_repo_name: JKRT/JuliaFMI

max_issues_repo_name: JKRT/JuliaFMI

max_forks_repo_name: JKRT/JuliaFMI


text:
# adding Inverse to have efficient ternary dot product and multiply
# allows for especially efficient Ma
import LinearAlgebra: dot, *, diag
function dot(x::AbstractVector, A::Inverse{<:Any, <:UpdatableCholesky}, y::AbstractVector)
	xp = copy(x)
	yp = x ≡ y ? xp : copy(y)
	dot!!(xp, A, yp)
end
# IDEA: could have parallel_threshold argument here
import LazyInverses: dot!!, inverse_cholesky_dot!!
function dot!!(x::AbstractVector, A::Inverse{<:Any, <:UpdatableCholesky}, y::AbstractVector)
	C = A.parent
	inverse_cholesky_dot!!(x, C, y)
end

using LazyInverses: inverse_cholesky_diag, inverse_cholesky_mul, inverse_cholesky_diag_mul
function diag(A::Inverse{<:Any, <:UpdatableCholesky})
	inverse_cholesky_diag(A.parent)
end

import LazyInverses: diag_mul
function diag_mul(X, A::Inverse{<:Any, <:UpdatableCholesky}, Y)
	C = A.parent
	inverse_cholesky_diag_mul(X, C, Y)
end

function *(X, A::Inverse{<:Any, <:UpdatableCholesky}, Y)
	inverse_cholesky_mul(X, A.parent, Y)
end


text:
"""
celsius_to_fahrenheit(celsius, ndigits::Int = 2)

Converts celsius to fahrenheit and round to 2 decimal places

# Example

```julia
celsius_to_fahrenheit(273.354, 3) == 524.037  # returns true
celsius_to_fahrenheit(273.354, 0) == 524.0    # returns true
celsius_to_fahrenheit(-40.0) == -40.0         # returns true
celsius_to_fahrenheit(-20.0) == -4.0          # returns true
celsius_to_fahrenheit(0) == 32.0              # returns true
celsius_to_fahrenheit(20) == 68.0             # returns true
```
"""
function celsius_to_fahrenheit(celsius, ndigits::Int = 2)
    return round((float(celsius) * 1.8) + 32.0; digits = ndigits)
end

"""
function celsius_to_kelvin(celsius, ndigits::Int = 2)

Converts celsius to kelvin and round to 2 decimal places

# Example

```julia
celsius_to_kelvin(273.354, 3) == 546.504    # returns true
celsius_to_kelvin(273.354, 0) == 547.0      # returns true
celsius_to_kelvin(0.0) == 273.15            # returns true
celsius_to_kelvin(20.0) == 293.15           # returns true
```

"""
function celsius_to_kelvin(celsius, ndigits::Int = 2)
    return round(float(celsius) + 273.15; digits = ndigits)
end

"""
fahrenheit_to_celsius(fahrenheit, ndigits::Int = 2)

Converts fahrenheit to celsius and round to 2 decimal places

# Example

```julia
fahrenheit_to_celsius(273.354, 3) == 134.086 # returns true
fahrenheit_to_celsius(273.354, 0) == 134.0   # returns true
fahrenheit_to_celsius(0.0) == -17.78         # returns true
fahrenheit_to_celsius(20.0) == -6.67         # returns true
fahrenheit_to_celsius(40.0) == 4.44          # returns true
fahrenheit_to_celsius(60.0) == 15.56         # returns true
fahrenheit_to_celsius(80.0) == 26.67         # returns true
```
"""
function fahrenheit_to_celsius(fahrenheit, ndigits::Int = 2)
    return round((float(fahrenheit) - 32.0) * 5 / 9; digits = ndigits)
end

"""
fahrenheit_to_kelvin(fahrenheit, ndigits::Int = 2)

Converts fahrenheit to kelvin and round to 2 decimal places

# Example

```julia
fahrenheit_to_kelvin(273.354, 3) == 407.236 # returns true
fahrenheit_to_kelvin(273.354, 0) == 407.0   # returns true
fahrenheit_to_kelvin(0) == 255.37           # returns true
fahrenheit_to_kelvin(20.0) == 266.48        # returns true
fahrenheit_to_kelvin(40.0) == 277.59        # returns true
fahrenheit_to_kelvin(60.0) == 288.71        # returns true
fahrenheit_to_kelvin(80.0) == 299.82        # returns true
```
"""
function fahrenheit_to_kelvin(fahrenheit, ndigits::Int = 2)
    return round(
        ((float(fahrenheit) - 32.0) * 5 / 9) + 273.15;
        digits = ndigits,
    )
end

"""
function kelvin_to_celsius(kelvin, ndigits::Int = 2)

Converts kelvin to celsius and round to 2 decimal places

# Example

```julia
kelvin_to_celsius(273.354, 3) == 0.204   # returns true
kelvin_to_celsius(273.354, 0) == 0.0     # returns true
kelvin_to_celsius(273.15) == 0.0         # returns true
kelvin_to_celsius(300) == 26.85          # returns true
```
"""
function kelvin_to_celsius(kelvin, ndigits::Int = 2)
    return round(float(kelvin) - 273.15; digits = ndigits)
end

"""
function kelvin_to_fahrenheit(kelvin, ndigits::Int = 2)

Converts kelvin to fahrenheit and round to 2 decimal places

# Example

```julia
kelvin_to_fahrenheit(273.354, 3) == 32.367  # returns true
kelvin_to_fahrenheit(273.354, 0) == 32.0    # returns true
kelvin_to_fahrenheit(273.15) == 32.0        # returns true
kelvin_to_fahrenheit(300) == 80.33          # returns true
```
"""
function kelvin_to_fahrenheit(kelvin, ndigits::Int = 2)
    return round(((float(kelvin) - 273.15) * 1.8) + 32.0; digits = ndigits)
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia