text:
module Neighborhoods

using ConstructionBase, StaticArrays

export Neighborhood, Window, AbstractKernelNeighborhood, Kernel,
       Moore, VonNeumann, AbstractPositionalNeighborhood, Positional, LayeredPositional

export neighbors, neighborhood, kernel, kernelproduct, offsets, positions, radius, distances

export setwindow, updatewindow, unsafe_updatewindow

export pad_axes, unpad_axes

export broadcast_neighborhood, broadcast_neighborhood!

"""
    Neighborhood

Neighborhoods define the pattern of surrounding cells in the "neighborhood"
of the current cell. The `neighbors` function returns the surrounding
cells as an iterable.

The main kinds of neighborhood are demonstrated below:

![Neighborhoods](https://raw.githubusercontent.com/cesaraustralia/DynamicGrids.jl/media/Neighborhoods.png)

Neighborhoods can be used in `NeighborhoodRule` and `SetNeighborhoodRule` -
the same shapes with different purposes. In a `NeighborhoodRule` the neighborhood specifies
which cells around the current cell are returned as an iterable from the `neighbors` function.
These can be counted, summed, compared, or multiplied with a kernel in an
`AbstractKernelNeighborhood`, using [`kernelproduct`](@ref).

In `SetNeighborhoodRule` neighborhoods give the locations of cells around the central cell,
as [`offsets`] and absolute [`positions`](@ref) around the index of each neighbor. These
can then be written to manually.
"""
abstract type Neighborhood{R,N,L} end

ConstructionBase.constructorof(::Type{<:T}) where T <: Neighborhood{R,N,L} where {R,N,L} =
    T.name.wrapper{R,N,L}

"""
    kernelproduct(rule::NeighborhoodRule})
    kernelproduct(hood::AbstractKernelNeighborhood)
    kernelproduct(hood::Neighborhood, kernel)

Returns the vector dot product of the neighborhood and the kernel,
although differing from `dot` in that the dot product is not take for
vector members of the neighborhood - they are treated as scalars.
"""
function kernelproduct end

"""
    radius(rule, [key]) -> Int

Return the radius of a rule or ruleset if it has one, otherwise zero.
"""
function radius end
radius(hood::Neighborhood{R}) where R = R

"""
    neighbors(x::Union{Neighborhood,NeighborhoodRule}}) -> iterable

Returns an indexable iterator for all cells in the neighborhood,
either a `Tuple` of values or a range.

Custom `Neighborhood`s must define this method.
"""
function neighbors end
neighbors(hood::Neighborhood) = begin
    map(i -> _window(hood)[i], window_indices(hood))
end

"""
    offsets(x) -> iterable

Returns an indexable iterable over all cells, containing `Tuple`s of
the index offset from the central cell.

Custom `Neighborhood`s must define this method.
"""
function offsets end
offsets(hood::Neighborhood) = offsets(typeof(hood))
_window(hood::Neighborhood) = hood._window

"""
    positions(x::Union{Neighborhood,NeighborhoodRule}}, cellindex::Tuple) -> iterable

Returns an indexable iterable, over all cells as `Tuple`s of each
index in the main array. Useful in `SetNeighborhoodRule` for
setting neighborhood values, or for getting values in an Aux array.
"""
function positions end
@inline positions(hood::Neighborhood, I::CartesianIndex) = positions(hood, Tuple(I))
@inline positions(hood::Neighborhood, I::Int...) = positions(hood, I)
@inline positions(hood::Neighborhood, I::Tuple) = map(o -> o .+ I, offsets(hood))

"""
    distances(hood::Neighborhood)

Get the center-to-center distance of each neighborhood position from the central cell,
so that horizontally or vertically adjacent cells have a distance of `1.0`, and a
diagonally adjacent cell has a distance of `sqrt(2.0)`.

Vales are calculated at compile time, so `distances` can be used inside rules with little
overhead.
"""
@generated function distances(hood::Neighborhood{N,R,L}) where {N,R,L}
    expr = Expr(:tuple, ntuple(i -> :(bd[bi[$i]]), L)...)
    return quote
        bd = window_distances(hood)
        bi = window_indices(hood)
        $expr
    end
end

@generated function window_indices(H::Neighborhood{R,N}) where {R,N}
    # Offsets can be anywhere in the window, specified with
    # all dimensions. Here we transform them to a linear index.
    bi = map(offsets(H)) do o
        # Calculate strides
        S = 2R + 1
        strides = ntuple(i -> S^(i-1), N)
        # Offset indices are centered, we want them as regular array indices
        # Return the linear index in the square window
        sum(map((i, s) -> (i + R) * s, o, strides)) + 1
    end
    return Expr(:tuple, bi...)
end

@generated function window_distances(hood::Neighborhood{R,N}) where {R,N}
    values = map(CartesianIndices(ntuple(_ -> SOneTo{2R+1}(), N))) do I
        sqrt(sum((Tuple(I) .- (R + 1)) .^ 2))
    end
    x = SArray(values)
    quote
        return $x
    end
end

Base.eltype(hood::Neighborhood) = eltype(_window(hood))
Base.length(hood::Neighborhood{<:Any,<:Any,L}) where L = L
Base.ndims(hood::Neighborhood{<:Any,N}) where N = N
# Note: size is radial, and may not relate to `length` in the same way
# as in an array. A neighborhood does not have to include all cells includeding
# in the area covered by `size` and `axes`.
Base.size(hood::Neighborhood{R,N}) where {R,N} = ntuple(_ -> 2R+1, N)
Base.axes(hood::Neighborhood{R,N}) where {R,N} = ntuple(_ -> SOneTo{2R+1}(), N)
Base.iterate(hood::Neighborhood, args...) = iterate(neighbors(hood), args...)
Base.getindex(hood::Neighborhood, i) = begin
    # @show _window(hood) window_indices(hood)
    getindex(_window(hood), window_indices(hood)[i])
end

"""
    Moore <: Neighborhood

    Moore(radius::Int=1; ndims=2)
    Moore(; radius=1, ndims=2)
    Moore{R}(; ndims=2)
    Moore{R,N}()

Moore neighborhoods define the neighborhood as all cells within a horizontal or
vertical distance of the central cell. The central cell is omitted.

Radius `R = 1`:

```
N = 1   N = 2

 ▄ ▄     █▀█
         ▀▀▀
```

Radius `R = 2`:

```
N = 1   N = 2

        █████
▀▀ ▀▀   ██▄██
        ▀▀▀▀▀
```

Using `R` and `N` type parameters removes runtime cost of generating the neighborhood,
compated to passing arguments/keywords.
"""
struct Moore{R,N,L,W} <: Neighborhood{R,N,L}
    _window::W
end
Moore(radius::Int=1; ndims=2) = Moore{radius,ndims}()
Moore(args...; radius=1, ndims=2) = Moore{radius,ndims}(args...)
Moore{R}(_window=nothing; ndims=2) where R = Moore{R,ndims,}(_window)
Moore{R,N}(_window=nothing) where {R,N} = Moore{R,N,(2R+1)^N-1}(_window)
Moore{R,N,L}(_window::W=nothing) where {R,N,L,W} = Moore{R,N,L,W}(_window)

@generated function offsets(::Type{<:Moore{R,N}}) where {R,N}
    exp = Expr(:tuple)
    for I in CartesianIndices(ntuple(_-> -R:R, N))
        if !all(map(iszero, Tuple(I)))
            push!(exp.args, :($(Tuple(I))))
        end
    end
    return exp
end
@inline setwindow(n::Moore{R,N,L}, win::W2) where {R,N,L,W2} = Moore{R,N,L,W2}(win)

"""
    VonNeumann(radius=1; ndims=2) -> Positional
    VonNeumann(; radius=1, ndims=2) -> Positional
    VonNeumann{R,N}() -> Positional

A Von Neuman neighborhood is a damond-shaped, omitting the central cell:

Radius `R = 1`:

```
N = 1   N = 2

 ▄ ▄     ▄▀▄
          ▀
```

Radius `R = 2`:

```
N = 1   N = 2

         ▄█▄
▀▀ ▀▀   ▀█▄█▀
          ▀
```

In 1 dimension it is identical to [`Moore`](@ref).

Using `R` and `N` type parameters removes runtime cost of generating the neighborhood,
compated to passing arguments/keywords.
"""
struct VonNeumann{R,N,L,W} <: Neighborhood{R,N,L}
    _window::W
end
VonNeumann(; radius=1, ndims=2) = VonNeumann(radius; ndims)
VonNeumann(radius, _window=nothing; ndims=2) = VonNeumann{radius,ndims}(_window)
VonNeumann{R}(_window=nothing; ndims=2) where R = VonNeumann{R,ndims}(_window)
function VonNeumann{R,N}(_window=nothing) where {R,N}
    L = 2sum(1:R) + 2R
    VonNeumann{R,N,L}(_window)
end
VonNeumann{R,N,L}(_window::W=nothing) where {R,N,L,W} = VonNeumann{R,N,L,W}(_window)

@inline setwindow(n::VonNeumann{R,N,L}, win::W2) where {R,N,L,W2} = VonNeumann{R,N,L,W2}(win)

@generated function offsets(::Type{T}) where {T<:VonNeumann{R,N}} where {R,N}
    offsets_expr = Expr(:tuple)
    rngs = ntuple(_ -> -R:R, N)
    for I in CartesianIndices(rngs)
        manhatten_distance = sum(map(abs, Tuple(I)))
        if manhatten_distance in 1:R
            push!(offsets_expr.args, Tuple(I))
        end
    end
    return offsets_expr
end

"""
    Window <: Neighborhood

    Window(; radius=1, ndims=2)
    Window{R}(; ndims=2)
    Window{R,N}()

A neighboorhood of radius R that includes the central cell.

Radius `R = 1`:

```
N = 1   N = 2

 ▄▄▄     ███
         ▀▀▀
```

Radius `R = 2`:

```
N = 1   N = 2

        █████
▀▀▀▀▀   █████
        ▀▀▀▀▀
```
"""
struct Window{R,N,L,W} <: Neighborhood{R,N,L}
    _window::W
end
Window(; radius=1, ndims=2) = Window{radius,ndims}(args...)
Window(R::Int, args...; ndims=2) = Window{R,ndims}(args...)
Window{R}(_window=nothing; ndims=2) where {R} = Window{R,ndims}(_window)
Window{R,N}(_window=nothing) where {R,N} = Window{R,N,(2R+1)^N}(_window)
Window{R,N,L}(_window::W=nothing) where {R,N,L,W} = Window{R,N,L,W}(_window)
Window(A::AbstractArray) = Window{(size(A, 1) - 1) ÷ 2,ndims(A)}()

# The central cell is included
@inline function offsets(::Type{<:Window{R,N}}) where {R,N}
    D = 2R + 1
    ntuple(i -> (rem(i - 1, D) - R, (i - 1) ÷ D - R), D^N)
end

distances(hood::Window) = Tuple(window_distances(hood))

@inline setwindow(::Window{R,N,L}, win::W2) where {R,N,L,W2} = Window{R,N,L,W2}(win)

window_indices(hood::Window{R,N}) where {R,N} = SOneTo{(2R + 1)^N}()

neighbors(hood::Window) = _window(hood)

"""
    AbstractKernelNeighborhood <: Neighborhood

Abstract supertype for kernel neighborhoods.

These can wrap any other neighborhood object, and include a kernel of
the same length and positions as the neighborhood.
"""
abstract type AbstractKernelNeighborhood{R,N,L,H} <: Neighborhood{R,N,L} end

neighbors(hood::AbstractKernelNeighborhood) = neighbors(neighborhood(hood))
offsets(::Type{<:AbstractKernelNeighborhood{<:Any,<:Any,<:Any,H}}) where H = offsets(H)
positions(hood::AbstractKernelNeighborhood, I::Tuple) = positions(neighborhood(hood), I)

"""
    kernel(hood::AbstractKernelNeighborhood) => iterable

Returns the kernel object, an array or iterable matching the length
of the neighborhood.
"""
function kernel end
kernel(hood::AbstractKernelNeighborhood) = hood.kernel

"""
    neighborhood(x) -> Neighborhood

Returns a neighborhood object.
"""
function neighborhood end
neighborhood(hood::AbstractKernelNeighborhood) = hood.neighborhood

"""
    kernelproduct(hood::AbstractKernelNeighborhood)
    kernelproduct(hood::Neighborhood, kernel)

Take the vector dot produce of the neighborhood and the kernel,
without recursion into the values of either. Essentially `Base.dot`
without recursive calls on the contents, as these are rarely what is
intended.
"""
function kernelproduct(hood::AbstractKernelNeighborhood)
    kernelproduct(neighborhood(hood), kernel(hood))
end
function kernelproduct(hood::Neighborhood{<:Any,<:Any,L}, kernel) where L
    sum = zero(first(hood))
    @simd for i in 1:L
        @inbounds sum += hood[i] * kernel[i]
    end
    return sum
end
function kernelproduct(hood::Window{<:Any,<:Any,L}, kernel) where L
    sum = zero(first(hood))
    @simd for i in 1:L
        @inbounds sum += _window(hood)[i] * kernel[i]
    end
    return sum
end

"""
    Kernel <: AbstractKernelNeighborhood

    Kernel(neighborhood, kernel)

Wrap any other neighborhood object, and includes a kernel of
the same length and positions as the neighborhood.
"""
struct Kernel{R,N,L,H,K} <: AbstractKernelNeighborhood{R,N,L,H}
    neighborhood::H
    kernel::K
end
Kernel(A::AbstractMatrix) = Kernel(Window(A), A)
function Kernel(hood::H, kernel::K) where {H<:Neighborhood{R,N,L},K} where {R,N,L}
    length(hood) == length(kernel) || _kernel_length_error(hood, kernel)
    Kernel{R,N,L,H,K}(hood, kernel)
end
function Kernel{R,N,L}(hood::H, kernel::K) where {R,N,L,H<:Neighborhood{R,N,L},K}
    Kernel{R,N,L,H,K}(hood, kernel)
end

function _kernel_length_error(hood, kernel)
    throw(ArgumentError("Neighborhood length $(length(hood)) does not match kernel length $(length(kernel))"))
end

function setwindow(n::Kernel{R,N,L,<:Any,K}, win) where {R,N,L,K}
    hood = setwindow(neighborhood(n), win)
    return Kernel{R,N,L,typeof(hood),K}(hood, kernel(n))
end

"""
    AbstractPositionalNeighborhood <: Neighborhood

Positional neighborhoods are tuples of coordinates that are specified in relation
to the central point of the current cell. They can be any arbitrary shape or size,
but should be listed in column-major order for performance.
"""
abstract type AbstractPositionalNeighborhood{R,N,L} <: Neighborhood{R,N,L} end

const CustomOffset = Tuple{Vararg{Int}}
const CustomOffsets = Union{AbstractArray{<:CustomOffset},Tuple{Vararg{<:CustomOffset}}}

"""
    Positional <: AbstractPositionalNeighborhood

    Positional(coord::Tuple{Vararg{Int}}...)
    Positional(offsets::Tuple{Tuple{Vararg{Int}}})
    Positional{O}()

Neighborhoods that can take arbitrary shapes by specifying each coordinate,
as `Tuple{Int,Int}` of the row/column distance (positive and negative)
from the central point.

The neighborhood radius is calculated from the most distant coordinate.
For simplicity the window read from the main grid is a square with sides
`2r + 1` around the central point.

The dimensionality `N` of the neighborhood is taken from the length of
the first coordinate, e.g. `1`, `2` or `3`.

Example radius `R = 1`:

```
N = 1   N = 2

 ▄▄      ▀▄
          ▀
```

Example radius `R = 2`:

```
N = 1   N = 2

         ▄▄
 ▀ ▀▀   ▀███
           ▀
```

Using the `O` parameter e.g. `Positional{((1, 2), (1, 1))}()` removes any
runtime cost of generating the neighborhood.
"""
struct Positional{O,R,N,L,W} <: AbstractPositionalNeighborhood{R,N,L}
    "A tuple of tuples of Int, containing 2-D coordinates relative to the central point"
    _window::W
end
Positional(args::CustomOffset...) = Positional(args)
function Positional(offsets::CustomOffsets, _window=nothing)
    Positional{offsets}(_window)
end
function Positional(offsets::O, _window=nothing) where O
    Positional{offsets}(_window)
end
function Positional{O}(_window=nothing) where O
    R = _absmaxcoord(O)
    N = length(first(O))
    L = length(O)
    Positional{O,R,N,L}(_window)
end
function Positional{O,R,N,L}(_window::W=nothing) where {O,R,N,L,W}
    Positional{O,R,N,L,W}(_window)
end

# Calculate the maximum absolute value in the offsets to use as the radius
function _absmaxcoord(offsets::Union{AbstractArray,Tuple})
    maximum(map(x -> maximum(map(abs, x)), offsets))
end

function ConstructionBase.constructorof(::Type{Positional{O,R,N,L,W}}) where {O,R,N,L,W}
    Positional{O,R,N,L}
end

offsets(::Type{<:Positional{O}}) where O = O

@inline function setwindow(n::Positional{O,R,N,L}, win::W2) where {O,R,N,L,W2}
    Positional{O,R,N,L,W2}(win)
end

"""
    LayeredPositional <: AbstractPositional

    LayeredPositional(layers::Positional...)

Sets of [`Positional`](@ref) neighborhoods that can have separate rules for each set.

`neighbors` for `LayeredPositional` returns a tuple of iterators
for each neighborhood layer.
"""
struct LayeredPositional{R,N,L,La,W} <: AbstractPositionalNeighborhood{R,N,L}
    "A tuple of custom neighborhoods"
    layers::La
    _window::W
end
LayeredPositional(layers::Positional...) = LayeredPositional(layers)
function LayeredPositional(layers::Tuple{Vararg{<:Positional}}, _window=nothing)
    R = maximum(map(radius, layers))
    N = ndims(first(layers))
    L = map(length, layers)
    LayeredPositional{R,N,L}(layers, _window)
end
function LayeredPositional{R,N,L}(layers, _window::W) where {R,N,L,W}
    # Child layers must have the same _window, and the
    # same R and L parameters. So we rebuild them all here.
    layers = map(l -> Positional{offsets(l),R,N,L}(_window), layers)
    LayeredPositional{R,N,L,typeof(layers),W}(layers, _window)
end

@inline neighbors(hood::LayeredPositional) = map(l -> neighbors(l), hood.layers)
@inline offsets(::Type{<:LayeredPositional{R,N,L,La}}) where {R,N,L,La} =
    map(p -> offsets(p), tuple_contents(La))
@inline positions(hood::LayeredPositional, args::Tuple) = map(l -> positions(l, args...), hood.layers)
@inline function setwindow(n::LayeredPositional{R,N,L}, win) where {R,N,L}
    LayeredPositional{R,N,L}(n.layers, win)
end

@inline Base.sum(hood::LayeredPositional) = map(sum, neighbors(hood))

function _subwindow(l::Neighborhood{R,N}, window) where {R,N}

    window_inds = ntuple(_-> R + 1, R + 1, N)

    vals = map(ps) do p
        window[p...]
    end
    L = (2R + 1) ^ N
    S = Tuple{ntuple(_ -> 2R + 1, N)...}
    return SArray{S,eltype(window),N,L}(vals)
end

# Utils

# Get the size of a neighborhood dimension from its radius,
# which is always 2r + 1.
@inline hoodsize(hood::Neighborhood{R}) where R = hoodsize(R)
@inline hoodsize(radius::Integer) = 2radius + 1

# Copied from StaticArrays. If they can do it...
Base.@pure function tuple_contents(::Type{X}) where {X<:Tuple}
    return tuple(X.parameters...)
end
tuple_contents(xs::Tuple) = xs

"""
    unsafe_readwindow(hood::Neighborhood, A::AbstractArray, I) => SArray

Get a single window square from an array, as an `SArray`, checking bounds.
"""
readwindow(hood::Neighborhood, A::AbstractArray, I::Int...) = readwindow(hood, A, I)
@inline function readwindow(hood::Neighborhood{R,N}, A::AbstractArray, I) where {R,N}
    for O in ntuple(_ -> (-R, R), N)
        edges = Tuple(I) .+ O
        map(I -> checkbounds(A, I...), edges)
    end
    return unsafe_readwindow(hood, A, I...)
end

"""
    unsafe_readwindow(hood::Neighborhood, A::AbstractArray, I) => SArray

Get a single window square from an array, as an `SArray`, without checking bounds.
"""
@inline unsafe_readwindow(hood::Neighborhood, A::AbstractArray, I::CartesianIndex) =
    unsafe_readwindow(hood, A, Tuple(I))
@inline unsafe_readwindow(hood::Neighborhood, A::AbstractArray, I::Int...) =
    unsafe_readwindow(hood, A, I)
@generated function unsafe_readwindow(
    ::Neighborhood{R,N}, A::AbstractArray{T,N}, I::NTuple{N,Int}
) where {T,R,N}
    R = 1
    S = 2R+1
    L = S^N
    sze = ntuple(_ -> S, N)
    vals = Expr(:tuple)
    nh = CartesianIndices(ntuple(_ -> -R:R, N))
    for i in 1:L
        Iargs = map(Tuple(nh[i]), 1:N) do nhi, n
            :(I[$n] + $nhi)
        end
        Iexp = Expr(:tuple, Iargs...)
        exp = :(@inbounds A[$Iexp...])
        push!(vals.args, exp)
    end

    sze_exp = Expr(:curly, :Tuple, sze...)
    return :(SArray{$sze_exp,$T,$N,$L}($vals))
end
@generated function unsafe_readwindow(
    ::Neighborhood{R,N1}, A::AbstractArray{T,N2}, I::NTuple{N3,Int}
) where {T,R,N1,N2,N3}
    throw(DimensionMismatch("neighborhood has $N1 dimensions while array has $N2 and index has $N3"))
end

# Reading windows without padding
# struct Padded{S,V}
#     padval::V
# end
# Padded{S}(padval::V) where {S,V} = Padded{S,V}(padval)

# @generated function unsafe_readwindow(
#     ::Neighborhood{R,N}, bounds::Padded{V}, ::AbstractArray{T,N}, I::NTuple{N,Int}
# ) where {T,R,N,V}
#     R = 1
#     S = 2R+1
#     L = S^N
#     sze = ntuple(_ -> S, N)
#     vals = Expr(:tuple)
#     for X in CartesianIndices(ntuple(_ -> -R:R, N))
#         if X in CartesianIndices(V)
#             # Generate indices for this position
#             Iargs = map(Tuple(X), 1:N) do x, n
#                 :(I[$n] + $x)
#             end
#             Iexp = Expr(:tuple, Iargs...)
#             push!(vals.args, :(@inbounds A[$Iexp...]))
#         else
#             push!(vals.args, :(bounds.padval))
#         end
#     end

#     sze_exp = Expr(:curly, :Tuple, sze...)
#     return :(SArray{$sze_exp,$T,$N,$L}($vals))
# end

"""
    updatewindow(x, A::AbstractArray, I...) => Neighborhood

Set the window of a neighborhood to values from the array A around index `I`.

Bounds checks will reduce performance, aim to use `unsafe_setwindow` directly.
"""
@inline function updatewindow(x, A::AbstractArray, i, I...)
    setwindow(x, readwindow(x, A, i, I...))
end

"""
    unsafe_setwindow(x, A::AbstractArray, I...) => Neighborhood

Set the window of a neighborhood to values from the array A around index `I`.

No bounds checks occur, ensure that A has padding of at least the neighborhood radius.
"""
@inline function unsafe_updatewindow(h::Neighborhood, A::AbstractArray, i, I...)
    setwindow(h, unsafe_readwindow(h, A, i, I...))
end

"""
    broadcast_neighborhood(f, hood::Neighborhood, As...)

Simple neighborhood application, where `f` is passed
each neighborhood in `A`, returning a new array.

The result is smaller than `A` on all sides, by the neighborhood radius.
"""
function broadcast_neighborhood(f, hood::Neighborhood, sources...)
    checksizes(sources...)
    ax = unpad_axes(first(sources), hood)
    sourceview = view(first(sources), ax...)
    broadcast(sourceview, CartesianIndices(ax)) do _, I
        applyneighborhood(f, sources, I)
    end
end

"""
    broadcast_neighborhood!(f, hood::Neighborhood{R}, dest, sources...)

Simple neighborhood broadcast where `f` is passed each neighborhood
of `src` (except padding), writing the result of `f` to `dest`.

`dest` must either be smaller than `src` by the neighborhood radius on all
sides, or be the same size, in which case it is assumed to also be padded.
"""
function broadcast_neighborhood!(f, hood::Neighborhood, dest, sources)
    checksizes(sources...)
    if axes(dest) === axes(src)
        ax = unpad_axes(src, hood)
        destview = view(dest, ax...)
        broadcast!(destview, CartesianIndices(ax)) do I
            applyneighborhood(f, sources, I)
        end
    else
        broadcast!(dest, CartesianIndices(unpad_axes(src, hood))) do I
            applyneighborhood(f, sources, I)
        end
    end
end

function checksizes(sources...)
    map(sources) do s
        size(s) === size(first(sources)) || throw(ArgumentError("Source array sizes must match"))
    end
end

function applyneighborhood(f, sources, I)
    hoods = map(sources) do s
        unsafe_updatewindow(hood, s, I)
    end
    f(hoods...)
end

"""
    pad_axes(A, hood::Neighborhood{R})
    pad_axes(A, radius::Int)

Add padding to axes.
"""
pad_axes(A, hood::Neighborhood{R}) where R = pad_axes(A, R)
function pad_axes(A, radius::Int)
    map(axes(A)) do axis
        firstindex(axis) - radius:lastindex(axis) + radius
    end
end
"""
    unpad_axes(A, hood::Neighborhood{R})
    unpad_axes(A, radius::Int)

Remove padding from axes.
"""
unpad_axes(A, hood::Neighborhood{R}) where R = unpad_axes(A, R)
function unpad_axes(A, radius::Int)
    map(axes(A)) do axis
        firstindex(axis) + radius:lastindex(axis) - radius
    end
end

end # Module Neighborhoods


text:
module MyModule
    export f

    g() = "Internal function"
    f() = println(g())
end

ext: jl

max_stars_repo_name: oist/skillpill-julia

max_issues_repo_name: oist/skillpill-julia

max_forks_repo_name: oist/skillpill-julia


text:
using InterPol
using Base.Test

@test interpolateAt([0], [0], 1) == 0

@test interpolateAt([0], [1], 1) == 1

@test interpolateAt([0,1,2,3], [0,1,2,3], 0) == 0

@test interpolateAt([0,1,2,3], [0,1,2,3], 4) == 4

@test interpolateAt([0,1,2], [0,1,4], 0) == 0

@test interpolateAt([0,1,2], [0,1,4], 1) == 1

@test interpolateAt([0,1,2], [0,1,4], 3) == 9

@test interpolateDerivativeAt([0], [0], 1) == 0

@test interpolateDerivativeAt([0], [1], 1) == 0

@test interpolateDerivativeAt([0,1,2,3], [0,1,2,3], 0) == 1

@test interpolateDerivativeAt([0,1,2,3], [0,1,2,3], 4) == 1

@test interpolateDerivativeAt([0,1,2], [0,1,4], 0) == 0

@test interpolateDerivativeAt([0,1,2], [0,1,4], 1) == 2

@test interpolateDerivativeAt([0,1,2], [0,1,4], 3) == 6


text:
function all_devices(sys, filter::Array)
    dev = Array{PowerSystems.PowerSystemDevice}([])

    for source in sys.generators
        if typeof(source) <: Array{<:PowerSystems.Generator}
            for d in source
                d.name in filter ? push!(dev,d) : continue
            end
        end
    end

    for d in sys.loads
        d.name in filter ? push!(dev,d) : continue
    end

    return dev
end

function all_devices(sys)
    dev = Array{PowerSystems.PowerSystemDevice}([])

    for source in sys.generators
        if typeof(source) <: Array{<:PowerSystems.Generator}
            for d in source
                push!(dev,d)
            end
        end
    end

    for d in sys.loads
        push!(dev,d)
    end

    return dev
end

#TODO: Make additional methods to handle other device types
function get_pg(m::JuMP.AbstractModel, gen::G, t::Int64) where G <: PowerSystems.ThermalGen
    return m.obj_dict[:p_th][gen.name,t]
end

function get_pg(m::JuMP.AbstractModel, gen::G, t::Int64) where G <: PowerSystems.RenewableCurtailment
    return m.obj_dict[:p_re][gen.name,t]
end

# Methods for accessing jump, moi, and optimizer variables
function get_all_vars(obj_dict)
    # get all variables in a jump model
    var_arays = [v.data for (k,v) in obj_dict if isa(v,JuMP.JuMPArray{JuMP.VariableRef}) ];
    vars = [i for arr in var_arays for i in arr]
end

function map_moi_opt_variables(model::JuMP.AbstractModel)
    # maps moi variablerefs (keys) to optimizer variable refs (values)
    vmap = model.moi_backend.model.optimizer.variable_mapping
    moi_optimizer_vmap = Dict()
    for v in vmap
        moi_optimizer_vmap[model.moi_backend.model.optimizer_to_model_map.varmap[v[1]]]=v[2]
    end
    return moi_optimizer_vmap
end

function map_jump_vars(model::JuMP.AbstractModel)
    # maps jump variablerefs (keys) to optimizer variable refs (values)
    vars = get_all_vars(model.obj_dict)
    moi_optimizer_vmap = map_moi_opt_variables(model)
    moivariables = [moi_optimizer_vmap[v.index] for v in vars]
    moivars = [moi_optimizer_vmap[v.index] for v in vars]
    return Dict(zip(vars,moivars))
end

function map_optimizer_vars(model::JuMP.AbstractModel)
    # maps optimizer variablerefs (keys) to jump variable refs (values)
    vars = get_all_vars(model.obj_dict)
    moi_optimizer_vmap = map_moi_opt_variables(model)
    moivariables = [moi_optimizer_vmap[v.index] for v in vars]
    moivars = [moi_optimizer_vmap[v.index] for v in vars]
    return Dict(zip(moivars,vars))
end

# Methods for accessing jump, moi, and optimizer constraintrefs
function get_all_constraints(obj_dict)
    # get all constraints in a jump model
    constraint_arrays = [v.data for (k,v) in obj_dict if isa(v,JuMP.JuMPArray{JuMP.ConstraintRef}) ];
    constraints = [i for arr in constraint_arrays for i in arr]
end

function map_moi_opt_constraints(model::JuMP.AbstractModel)
    # maps moi constraintrefs (keys) to optimizer constraint refs (values)
    cmap = model.moi_backend.model.optimizer.constraint_mapping
    moi_optimizer_map = Dict()
    for f in fieldnames(typeof(cmap))
        ctype = getfield(cmap,f)
        for c in ctype
            moi_optimizer_map[model.moi_backend.model.optimizer_to_model_map.conmap[c[1]]]=c[2]
        end
    end
    return moi_optimizer_map
end

function map_jump_constraints(model::JuMP.AbstractModel)
    # maps jump constraints (keys) to optimizer constraint refs (values)
    constraints = get_all_constraints(model.obj_dict)
    moi_optimizer_map = map_moi_opt_constraints(model)
    moiconstraints = [moi_optimizer_map[c.index] for c in constraints]
    return Dict(zip(constraints,moiconstraints))
end

function map_optimizer_constraints(model::JuMP.AbstractModel)
    # maps optimizer constraint refs (keys) to jump constraints (values)
    constraints = get_all_constraints(model.obj_dict)
    moi_optimizer_map = map_moi_opt_constraints(model)
    moiconstraints = [moi_optimizer_map[c.index] for c in constraints]
    return Dict(zip(moiconstraints,constraints))
end

function create_result_dict(jump_array, ::Type{Int64})

    d = Dict{String, Array{Int64}}()

    for var in jump_array.axes[1]
        arr = Int64[]

        for t in jump_array.axes[2]
            push!(arr, Int64(round(JuMP.result_value(jump_array[var, t]))))
        end

        d[var] = arr
    end

    return d

end

function create_result_dict(jump_array, ::Type{Float64})

    d = Dict{String, Array{Float64}}()

    for var in jump_array.axes[1]
        arr = Float64[]

        for t in jump_array.axes[2]
            push!(arr, Float64(JuMP.result_value(jump_array[var, t])))
        end

        d[var] = arr
    end

    return d

end

function create_result_dict(jump_array, k)

    if k in [:on_th, :start_th, :stop_th]
        return create_result_dict(jump_array, Int64)
    elseif k in [:p_th, :fbr, :p_re, :p_cl, :p_rsv, :q_th, :q_cl, :q_re]
        return create_result_dict(jump_array, Float64)
    end

end

function get_model_result(pspom::PS.PowerOperationModel)

    d = Dict{Symbol, DataFrame}()
    for (k, v) in pspom.model.obj_dict
        if typeof(v) <: JuMPArray{VariableRef}
            d[k] = create_result_dict(v, k)
        end
    end

    return d

end

function get_previous_value_df(res::DataFrame)
    res[:period] = 1:size(res, 1)
    var_res = melt(res[end,:], :period, variable_name = :Device)
    return var_res
end

function get_previous_value(res::DataFrame)
    var_res = get_previous_value_df(res)
    prev_val = Dict(zip(map(String,var_res[:Device]),var_res[:value]))
    return prev_val
end

function commitment_duration(res::Dict, initial,  transition::Symbol, minutes_per_step = 60)
    # res is the results df from the previous solution
    # initial is the dict that defines the initial status (e.g. initialonduration)
    # transition is the variable of commitment transition {:start_th or :stop_th}

    last_period = size(res[:on_th],1)

    on_devices = get_previous_value_df(res[:on_th])
    if transition == :start_th
        status = 1
    elseif transition == :stop_th
        status = 0
    end
    off_devices = copy(on_devices)
    off_devices = off_devices[off_devices.value.!=status,[:Device]]
    off_devices.value = 0.0
    on_devices = on_devices[on_devices.value.==status,[:Device]]

    initial = melt(DataFrame(initial), variable_name = :Device)
    initial.value = initial.value .+ (last_period * minutes_per_step/60)

    # for devices that have changed status in the last step, calculate how long they have been at their current status
    res_df = copy(res[transition])
    res_df[:period] = 1:size(res_df, 1)

    res_df = melt(res_df, :period, variable_name = :Device)

    res_df = join(res_df,on_devices, on = :Device)
    res_df = by(res_df[res_df[:value] .== 1 ,[:Device,:period]], :Device, df -> DataFrames.tail(df[[:period]],1))

    if size(res_df,1) > 0
        res_df.value  = ((last_period + 1) .- res_df.period) .* minutes_per_step/60
        res_df = join(on_devices,res_df[[:Device,:value]], on = :Device, kind = :outer)
        res_df[findall(ismissing,res_df.value),:] = join(initial, res_df[findall(ismissing,res_df.value),[:Device]], on=:Device)
    else
        # for everything else, add the current step periods to initial status
        res_df = copy(on_devices)
        res_df = join(initial,res_df,on=:Device)
    end
    res_df = vcat(res_df,off_devices)

    return Dict(zip(map(String,res_df.Device),res_df.value))
end

function collect_results(simulation_results)
    dfs = []
    vars = []

    for (d,step) in simulation_results
        for (v,df) in step
            df.Date = d
            df.period = 1:size(df,1)
            df.DateTime = [(DateTime(df.Date[ix])+Hour(df.period[ix]-1)) for ix in df.period]
            push!(dfs,df)
            push!(vars,v)
        end
    end

    sim_res = Dict()
    for var in unique(vars)
        ids = findall(vars -> vars == var,vars)
        sim_res[var] = copy(dfs[ids[1]])
        for id in ids[2:end]
            sim_res[var] = vcat(sim_res[var],dfs[id])
        end
    end
    return(sim_res)
end


text:
"""
`pandoc2html(formatted::AbstractString, doc::WeaveDoc)`

Convert output from pandoc markdown to html using Weave.jl template
"""
function pandoc2html(formatted::AbstractString, doc::WeaveDoc, outname::AbstractString)
  html_template = joinpath(Pkg.dir("Weave"), "templates/pandoc_skeleton.html")
  css_template = joinpath(Pkg.dir("Weave"), "templates/pandoc_skeleton.css")

  path, wsource = splitdir(abspath(doc.source))
  wversion = string(Pkg.installed("Weave"))
  wtime =  string(Date(now()))

  #Change path for pandoc
  old_wd = pwd()
  cd(doc.cwd)
  html =""
  outname = basename(outname)

  try
    pandoc_out, pandoc_in, proc = readandwrite(`pandoc -R -s --mathjax="" --self-contained --highlight-style=tango
    --template $html_template -c $css_template
     -V wversion=$wversion -V wtime=$wtime -V wsource=$wsource
     -o $outname`)
    println(pandoc_in, formatted)

    close(pandoc_in)
    proc_output = readall(pandoc_out)
    cd(old_wd)
  catch e
    cd(old_wd)
    warn("Error converting document to HTML")
    throw(e)
  end
end

"""
`pandoc2pdf(formatted::AbstractString, doc::WeaveDoc)`

Convert output from pandoc markdown to pdf using Weave.jl template
"""
function pandoc2pdf(formatted::AbstractString, doc::WeaveDoc, outname::AbstractString)

  header_template = joinpath(Pkg.dir("Weave"), "templates/pandoc_header.txt")

  path, wsource = splitdir(abspath(doc.source))
  wversion = string(Pkg.installed("Weave"))
  wtime =  Date(now())
  outname = basename(outname)

  #Change path for pandoc
  old_wd = pwd()
  cd(doc.cwd)
  html =""

  info("Done executing code. Running xelatex")
  try
    pandoc_out, pandoc_in, proc = readandwrite(`pandoc -R -s  --latex-engine=xelatex --highlight-style=tango
     --include-in-header=$header_template
     -V fontsize=12pt -o $outname`)
    println(pandoc_in, formatted)

    close(pandoc_in)
    proc_output = readall(pandoc_out)
    cd(old_wd)
  catch e
    cd(old_wd)
    warn("Error converting document to pdf")
    throw(e)
  end
end


text:
type CircleShape <: Drawable
    ptr::Ptr{Void}
    _texture::Texture

    function CircleShape(ptr::Ptr{Void})
        c = new(ptr)
        finalizer(c, destroy)
        c
    end
end

function CircleShape()
    CircleShape(ccall((:sfCircleShape_create, libcsfml_graphics), Ptr{Void}, ()))
end

function copy(shape::CircleShape)
    return CircleShape(ccall((:sfCircleShape_copy, libcsfml_graphics), Ptr{Void}, (Ptr{Void},), shape.ptr))
end

function destroy(shape::CircleShape)
    ccall((:sfCircleShape_destroy, libcsfml_graphics), Void, (Ptr{Void},), shape.ptr)
end

function set_position(shape::CircleShape, position::Vector2f)
    ccall((:sfCircleShape_setPosition, libcsfml_graphics), Void, (Ptr{Void}, Vector2f,), shape.ptr, position)
end

function set_rotation(shape::CircleShape, angle::Real)
    ccall((:sfCircleShape_setRotation, libcsfml_graphics), Void, (Ptr{Void}, Cfloat,), shape.ptr, angle)
end

function set_scale(shape::CircleShape, scale::Vector2f)
    ccall((:sfCircleShape_setScale, libcsfml_graphics), Void, (Ptr{Void}, Vector2f,), shape.ptr, scale)
end

function set_origin(shape::CircleShape, origin::Vector2f)
    ccall((:sfCircleShape_setOrigin, libcsfml_graphics), Void, (Ptr{Void}, Vector2f,), shape.ptr, origin)
end

function set_radius(shape::CircleShape, radius::Real)
    ccall((:sfCircleShape_setRadius, libcsfml_graphics), Void, (Ptr{Void}, Cfloat,), shape.ptr, radius)
end

function set_fillcolor(shape::CircleShape, color::Color)
    ccall((:sfCircleShape_setFillColor, libcsfml_graphics), Void, (Ptr{Void}, Color,), shape.ptr, color)
end

function set_outlinecolor(shape::CircleShape, color::Color)
    ccall((:sfCircleShape_setOutlineColor, libcsfml_graphics), Void, (Ptr{Void}, Color,), shape.ptr, color)
end

function set_outline_thickness(shape::CircleShape, thickness::Real)
    ccall((:sfCircleShape_setOutlineThickness, libcsfml_graphics), Void, (Ptr{Void}, Cfloat,), shape.ptr, thickness)
end

function get_origin(shape::CircleShape)
    ccall((:sfCircleShape_getOrigin, libcsfml_graphics), Vector2f, (Ptr{Void},), shape.ptr)
end

function get_scale(shape::CircleShape)
    ccall((:sfCircleShape_getScale, libcsfml_graphics), Vector2f, (Ptr{Void},), shape.ptr)
end

function get_rotation(shape::CircleShape)
    ccall((:sfCircleShape_getRotation, libcsfml_graphics), Cfloat, (Ptr{Void},), shape.ptr)
end

function get_position(shape::CircleShape)
    ccall((:sfCircleShape_getPosition, libcsfml_graphics), Vector2f, (Ptr{Void},), shape.ptr)
end

function get_radius(shape::CircleShape)
    ccall((:sfCircleShape_getRadius, libcsfml_graphics), Cfloat, (Ptr{Void},), shape.ptr)
end

function set_pointcount(shape::CircleShape, count::Int)
    ccall((:sfCircleShape_setPointCount, libcsfml_graphics), Void, (Ptr{Void}, UInt32,), shape.ptr, count)
end

function get_fillcolor(shape::CircleShape)
    ccall((:sfCircleShape_getFillColor, libcsfml_graphics), Color, (Ptr{Void},), shape.ptr)
end

function get_outlinecolor(shape::CircleShape)
    ccall((:sfCircleShape_getOutlineColor, libcsfml_graphics), Color, (Ptr{Void},), shape.ptr)
end

function get_outline_thickness(shape::CircleShape)
    ccall((:sfCircleShape_getOutlineThickness, libcsfml_graphics), Cfloat, (Ptr{Void},), shape.ptr)
end

function get_pointcount(shape::CircleShape)
    Int(ccall((:sfCircleShape_getPointCount, libcsfml_graphics), Csize_t, (Ptr{Void},), shape.ptr))
end

function get_point(shape::CircleShape, index::Integer)
    ccall((:sfCircleShape_getPoint, libcsfml_graphics), Vector2f, (Ptr{Void}, Csize_t), shape, index)
end

function move(shape::CircleShape, offset::Vector2f)
    ccall((:sfCircleShape_move, libcsfml_graphics), Void, (Ptr{Void}, Vector2f,), shape.ptr, offset)
end

function scale(shape::CircleShape, factors::Vector2f)
    ccall((:sfCircleShape_scale, libcsfml_graphics), Void, (Ptr{Void}, Vector2f,), shape.ptr, factors)
end

function rotate(shape::CircleShape, angle::Real)
    ccall((:sfCircleShape_rotate, libcsfml_graphics), Void, (Ptr{Void}, Cfloat,), shape.ptr, angle)
end

function set_texture(shape::CircleShape, texture::Texture)
    ccall((:sfCircleShape_setTexture, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Void}, Int32,), shape.ptr, texture.ptr, false)
    shape._texture = texture
end

function set_texture_rect(shape::CircleShape, rect::IntRect)
    ccall((:sfCircleShape_setTextureRect, libcsfml_graphics), Void, (Ptr{Void}, IntRect,), shape.ptr, rect)
end

function get_texture(shape::CircleShape)
    Texture(call((:sfCircleShape_getTexture, libcsfml_graphics), Ptr{Void}, (Ptr{Void},), shape.ptr))
end

function get_texture_rect(shape::CircleShape)
    ccall((:sfCircleShape_getTextureRect, libcsfml_graphics), IntRect, (Ptr{Void},), shape.ptr)
end

function get_localbounds(shape::CircleShape)
    ccall((:sfCircleShape_getLocalBounds, libcsfml_graphics), FloatRect, (Ptr{Void},), shape.ptr)
end

function get_globalbounds(shape::CircleShape)
    ccall((:sfCircleShape_getGlobalBounds, libcsfml_graphics), FloatRect, (Ptr{Void},), shape.ptr)
end


text:
# Normal distribution

import StatsFuns
export Normal

import Base: eltype

@measure Normal(μ,σ) ≃ (1/sqrt2π) * Lebesgue(X)

function logdensity(d::Normal{P} , x::X) where {P <: NamedTuple{(:μ, :σ)}, X}
    return - log(d.par.σ)  - (x - d.par.μ)^2 / (2 * d.par.σ^2)
end

# Standard normal

function logdensity(d::Normal{EmptyNamedTuple} , x::X) where {X}
    return - x^2 / 2
end

# Normal() = Normal{EmptyNamedTuple,Real}(NamedTuple())

sampletype(::Normal{NamedTuple{(),Tuple{}}}) = Real

Base.rand(μ::Normal{EmptyNamedTuple}) = randn()

≪(::Normal, ::Lebesgue{X}) where X <: Real = true
representative(::Normal) = Lebesgue(Real)


text:
# 1-dim version of the code

using Winston

#basic parameters
gamma = 1.4
cfl = 0.8
dt = 1.0e-5
dtp = dt

nzones = 100
tend = 0.2

type data1d
    x::Array{Float64,1} #cell centers
    xi::Array{Float64,1} #cell LEFT interfaces

    rho::Array{Float64,1}
    rhop::Array{Float64,1}
    rhom::Array{Float64,1}

    vel::Array{Float64,1}
    velp::Array{Float64,1}
    velm::Array{Float64,1}

    eps::Array{Float64,1}
    epsp::Array{Float64,1}
    epsm::Array{Float64,1}

    press::Array{Float64,1}
    pressp::Array{Float64,1}
    pressm::Array{Float64,1}

    q::Array{Float64,2} #conserved quantities
    qp::Array{Float64,2}
    qm::Array{Float64,2}

    n::Int64
    g::Int64 #ghost cells

    function data1d(nzones::Int64)

        x = zeros(nzones)
        xi = zeros(nzones)

        rho = zeros(nzones)
        rhop = zeros(nzones)
        rhom = zeros(nzones)

        vel = zeros(nzones)
        velp = zeros(nzones)
        velm = zeros(nzones)

        eps = zeros(nzones)
        epsp = zeros(nzones)
        epsm = zeros(nzones)

        press = zeros(nzones)
        pressp = zeros(nzones)
        pressm = zeros(nzones)

        q = zeros(nzones, 3)
        qp = zeros(nzones, 3)
        qm = zeros(nzones, 3)
        n = nzones
        g = 3

        new(x,
            xi,
            rho,
            rhop,
            rhom,
            vel,
            velp,
            velm,
            eps,
            epsp,
            epsm,
            press,
            pressp,
            pressm,
            q,
            qp,
            qm,
            n,
            g)
    end
end

#1-dim grid
function grid_setup(self::data1d, xmin, xmax)
    dx = (xmax - xmin) / (self.n -self.g*2 - 1)
    xmin = xmin - self.g*dx
    xmax = xmax + self.g*dx

    for i = 1:self.n
        self.x[i] = xmin + i*dx
    end

    for i = 1:self.n
        self.xi[i] = self.x[i] - 0.5dx
    end

    return self
end

#Shoctube initial data
function setup_tube(self::data1d)
    rchange = 0.5(self.x[self.n - self.g] - self.x[self.g + 1])

    rho1 = 1.0
    rho2 = 0.125
    press1 = 1.0
    press2 = 0.1

    for i = 1:self.n
        if self.x[i] < rchange
            self.rho[i] = rho1
            self.press[i] = press1
            self.eps[i] = press1/rho1/(gamma - 1.0)
            self.vel[i] = 0.0
        else
            self.rho[i] = rho2
            self.press[i] = press2
            self.eps[i] = press2/rho2/(gamma - 1.0)
            self.vel[i] = 0.0
        end
    end

    return self
end

#Shoctube initial data
function setup_blast(self::data1d)

    rho1 = 0.4
    press1 = 0.2
    press2 = 0.99

    self.rho[:] = rho1*ones(self.n)
    self.press[:] = press1*ones(self.n)
    self.eps[:] = press1./rho1./(gamma - 1.0)
    self.vel[:] = zeros(self.n)

    mid = int(self.n/2)
    self.press[mid-2:mid+2] = press2
    self.eps[mid-2:mid+2] = press2./rho1./(gamma - 1.0)

    return self
end

function apply_bcs(hyd::data1d)

    #arrays starting from zero
    #       |g                  |n-g #
    #[0 1 2 x x x  .....  x x x 7 8 9]

    #arrays starting from 1
    #     |g                  |n-g    #
    #[1 2 3 x x x  .....  x x x 8 9 10]
    hyd.rho[1:hyd.g] = hyd.rho[hyd.g+1]
    hyd.vel[1:hyd.g] = hyd.vel[hyd.g+1]
    hyd.eps[1:hyd.g] = hyd.eps[hyd.g+1]
    hyd.press[1:hyd.g] = hyd.press[hyd.g+1]

    hyd.rho[(hyd.n-hyd.g+1) : hyd.n] = hyd.rho[hyd.n-hyd.g]
    hyd.vel[(hyd.n-hyd.g+1) : hyd.n] = hyd.vel[hyd.n-hyd.g]
    hyd.eps[(hyd.n-hyd.g+1) : hyd.n] = hyd.eps[hyd.n-hyd.g]
    hyd.press[(hyd.n-hyd.g+1) : hyd.n] = hyd.press[hyd.n-hyd.g]

    return hyd
end

#equation of state (ideal gas)
function eos_press(rho, eps, gamma)
    press = (gamma - 1.0) .* rho .* eps
    return press
end

#soundspeed
function eos_cs2(rho, eps, gamma)
    prs = (gamma - 1.0) .* rho .* eps
    dpde = (gamma - 1.0) .* rho
    dpdrho = (gamma - 1.0) .* eps
    cs2 = dpdrho .+ dpde .* prs ./ (rho + 1.0e-30).^2.0
    return cs2
end

#1-dim
function prim2con(rho::AbstractVector,
                  vel::AbstractVector,
                  eps::AbstractVector)

    q = zeros(size(rho, 1), 3)
    q[:,1] = rho
    q[:,2] = rho .* vel
    q[:,3] = rho .* eps .+ 0.5rho .* vel.^2.0

    return q
end

#1-dim
function con2prim(q)
    rho = vec(q[:,1])
    vel = vec(q[:,2] ./ rho)
    eps = vec(q[:,3] ./ rho - 0.5vel.^2.0)
    press = eos_press(rho, eps, gamma)

    return rho, eps, press, vel
end

#time step calculation
function calc_dt(hyd::data1d, dtp)
    cs = sqrt(eos_cs2(hyd.rho, hyd.eps, gamma))
    dtnew = 1.0
    for i = (hyd.g+1):(hyd.n-hyd.g+1)
        dtnew = min(dtnew, (hyd.x[i+1] - hyd.x[i]) / max(abs(hyd.vel[i]+cs[i]), abs(hyd.vel[i]-cs[i])))
    end

    dtnew = min(cfl*dtnew, 1.05*dtp)

    return dtnew
end

function minmod(a,b)
    if a*b < 0.0
        return 0.0
    elseif abs(a) < abs(b)
        return a
    else
        return b
    end
end

signum(x,y) = y >= 0.0 ? abs(x) : -abs(x)

function tvd_mc_reconstruction(n, g, f, x, xi)
    fp = zeros(n)
    fm = zeros(n)

    for i = g:(n-g+2)
        dx_up = x[i] - x[i-1]
        dx_down = x[i+1] - x[i]
        dx_m = x[i] -xi[i]
        dx_p = xi[i+1] - x[i]
        df_up = (f[i]-f[i-1]) / dx_up
        df_down = (f[i+1]-f[i]) / dx_down

        if df_up*df_down < 0.0
            delta = 0.0
        else
            delta = signum(min(2.0abs(df_up), 2.0abs(df_down), 0.5(abs(df_up)+abs(df_down))), df_up + df_down)
        end

        fp[i] = f[i] + delta*dx_p
        fm[i] = f[i] - delta*dx_m
    end

    return fp, fm
end

function reconstruct(hyd::data1d)

    hyd.rhop, hyd.rhom = tvd_mc_reconstruction(hyd.n,
                                               hyd.g,
                                               hyd.rho,
                                               hyd.x,
                                               hyd.xi)
    hyd.epsp, hyd.epsm = tvd_mc_reconstruction(hyd.n,
                                               hyd.g,
                                               hyd.eps,
                                               hyd.x,
                                               hyd.xi)
    hyd.velp, hyd.velm = tvd_mc_reconstruction(hyd.n,
                                               hyd.g,
                                               hyd.vel,
                                               hyd.x,
                                               hyd.xi)

    hyd.pressp = eos_press(hyd.rhop, hyd.epsp, gamma)
    hyd.pressm = eos_press(hyd.rhom, hyd.epsm, gamma)

    hyd.qp = prim2con(hyd.rhop, hyd.velp, hyd.epsp)
    hyd.qm = prim2con(hyd.rhom, hyd.velm, hyd.epsm)

    return hyd
end

#Load solvers
include("rsolvers1d.jl")

function calc_rhs(hyd::data1d)
    #reconstruction and prim2con
    hyd = reconstruct(hyd)

    #compute flux difference
    fluxdiff = hlle(hyd)

    #return RHS = -fluxdiff
    return -fluxdiff
end

###############
# main program
###############

#initialize
hyd = data1d(nzones)

#set up grid
hyd = grid_setup(hyd, 0.0, 1.0)

#set up initial data
hyd = setup_tube(hyd)

#get initial timestep
dt = calc_dt(hyd, dt)

#initial prim2con
hyd.q = prim2con(hyd.rho, hyd.vel, hyd.eps)

t = 0.0
i = 1

tsave = collect(0.0:0.01:0.2)

k = 1
while t < tend

    if i % 10 == 0
        #sleep(1.0)
        #output
        println("$i $t $dt")
        p=plot(hyd.x, hyd.rho, "r-")
        p=oplot(hyd.x, hyd.vel, "b-")
        p=oplot(hyd.x, hyd.press, "g-")
        display(p)
    end

    if (t >= tsave[k])
        println("t = $t")
        wmatr = zeros(nzones, 5)

        wmatr[:,1] = hyd.x
        wmatr[:,2] = hyd.rho
        wmatr[:,3] = hyd.vel
        wmatr[:,4] = hyd.eps
        wmatr[:,5] = hyd.press

        writecsv("step_$k.csv", wmatr)

        k += 1
    end

    #calculate new timestep
    dt = calc_dt(hyd, dt)

    #save old state
    hydold = hyd
    qold = hyd.q

    #calc rhs
    k1 = calc_rhs(hyd)
    #calculate intermediate step
    hyd.q = qold + 0.5dt*k1
    #con2prim
    hyd.rho, hyd.eps, hyd.press, hyd.vel = con2prim(hyd.q)
    #boundaries
    hyd = apply_bcs(hyd)

    #calc rhs
    k2 = calc_rhs(hyd)
    #apply update
    hyd.q = qold + dt*(0.5k1 + 0.5k2)
    #con2prim
    hyd.rho, hyd.eps, hyd.press, hyd.vel = con2prim(hyd.q)
    #apply bcs
    hyd = apply_bcs(hyd)

    #update time
    t += dt
    i += 1

end

max_stars_repo_name: gzahl/hydro

max_issues_repo_name: gzahl/hydro

max_forks_repo_name: gzahl/hydro


text:
using DSP
using Unitful

export samplerate, set_default_samplerate!, mix, envelope, silence, noise,
  highpass, lowpass, bandpass, bandstop, tone, ramp, harmonic_complex, amplify,
  rampon, rampoff, fadeto, irn, normpower, dc_offset

"""
    mix(x,...)

Add several sounds together so that they play at the same time.

Unlike normal addition, this acts as if each sound is padded with
zeros at the end so that the lengths of all sounds match.

With one argument `x`, this returns a function `f(y)` that calls mix(x,y)
"""
mix(xs...) = soundop(+,xs...)
mix(x) = y -> mix(x,y)

"""
    envelope([x],y)

Mutliply several sounds together. Typically used to apply an
amplitude envelope x to sound y.

Unlike normal multiplication, this acts as if each sound is padded with
ones at the end so that the lengths of all sounds match.

With one argument `x`, this returns a function `f(y)` that calls envelope(x,y)
"""
envelope(x,y) = soundop(*,x,y)
envelope(x) = y -> envelope(x,y)

function soundop(op,x_in...)
  xs = promote(x_in...)
  len = maximum(nframes.(xs))

  sorted = sortperm(collect(nframes.(xs)),rev=true)
  y = similar(xs[1],(len,nchannels(xs[1])))
  y .= xs[sorted[1]]

  for i in sorted[2:end]
    y[1:nframes(xs[i]),:] .= op.(y[1:nframes(xs[i]),:],xs[i][:,:])
  end

  y
end

"""
    silence(length;[rate=samplerate()])

Creates period of silence of the given length (in seconds).
"""
function silence(length;rate=samplerate())
  Sound(t -> fill(0.0,size(t)),length,asseconds=false,rate=rate)
end

"""
    dc_offset(length;[rate=44100Hz])

Creates a DC offset of unit value.

In other words, this just returns samples all with the value 1: `silence` is
to `zeros` just as `dc_offset` is to `ones`.

Though this technically constructs a sound it is normally only used
in combination with `envelope`, because it produces no audible sound.
For example, it could be used to transition from a unmodulated to an amplitude
modulated noise:

    env = dc_offset(2s) |> fadeto(tone(8Hz,2s)) |> ramp(250ms)
    sound = noise(duration(env)) |> envelope(env) |> normpower
"""
function dc_offset(length;rate=samplerate())
  Sound(t -> fill(1.0,size(t)),length,asseconds=false,rate=rate)
end

"""
    noise(length;[rate=samplerate()],[rng=RandomDevice()])

Creates a period of white noise of the given length (in seconds).

"""
function noise(len;rate=samplerate(),rng=RandomDevice())
  Sound(i -> 1.0-2.0rand(rng,length(i)),len,asseconds=false,rate=rate)
end

"""
    tone(freq,length;[rate=samplerate()],[phase=0])

Creates a pure tone of the given frequency and length (in seconds).

"""
function tone(freq,len=Inf;rate=samplerate(),phase=0.0)
  Sound(t -> sin.(2π.*t .* ustrip(inHz(freq)) .+ phase),len,rate=rate)
end

function complex_cycle(f0,harmonics,amps,rate,phases)
  @assert all(0 .<= phases) && all(phases .< 2π)
	n = maximum(harmonics)+1

  # generate single cycle of complex
  cycle = silence((1/f0),rate=rate)

	highest_freq = tone(f0,2n*length(cycle)*frames;rate=rate)

	for (amp,harm,phase) in zip(amps,harmonics,phases)
		phase_offset = round(Int,n*phase/2π*rate/f0)
    wave = highest_freq[(1:length(cycle)) * (n-harm) + phase_offset]
		cycle += amp*wave[1:length(cycle)]
	end

  cycle
end

"""
    harmonic_complex(f0,harmonics,amps,length,
                     [rate=samplerate()],[phases=zeros(length(harmonics))])

Creates a harmonic complex of the given length, with the specified harmonics
at the given amplitudes. This implementation is somewhat superior
to simply summing a number of pure tones generated using `tone`, because
it avoids beating in the sound that may occur due floating point errors.

"""
function harmonic_complex(f0,harmonics,amps,len=Inf;
						              rate=samplerate(),
                          phases=zeros(length(harmonics)))
  cycle = complex_cycle(inHz(f0),harmonics,amps,
                        inHz(Int,rate),phases)

  N = size(cycle,1)
  Sound(i -> cycle[(i.-1) .% N + 1],len,asseconds=false,rate=rate)
end

"""
    irn(n,freq,length;[g=1],[rate=samplerate()],[rng=Base.GLOBAL_RNG])

Creates an iterated ripple ``y_n(t)`` for a noise ``y_0(t)`` according to
the following formula.

``
y_n(t) = y_{n-1}(t) + g⋅y_{n-1}(t-1/freq)
``
"""
function irn(n,λ,length;g=1,rate=samplerate(),rng=Base.GLOBAL_RNG)
  irn_helper(noise(length,rate=rate,rng=rng),n,λ,g,rng)
end

function irn_helper(source,n,λ,g,rng)
  if n == 0
    source
  else
    irn_helper(mix(source,[silence(λ); g * source]),n-1,λ,g,rng)
  end
end

"""
    bandpass([x],low,high;[order=5])

Band-pass filter the sound at the specified frequencies.

Filtering uses a butterworth filter of the given order.

With two positional arguments this returns a function f(x)
that calls bandpass(x,low,hight;order=order)
"""
bandpass(x,low,high;order=5) = filter_helper(x,low,high,Bandpass,order)
bandpass(low,high;order=5) = x -> bandpass(x,low,high,order=order)

"""
    bandstop([x],low,high,[order=5])

Band-stop filter of the sound at the specified frequencies.

Filtering uses a butterworth filter of the given order.

With two positional arguments this returns a function f(x)
that calls bandstop(x,low,hight;order=order)
"""
bandstop(x,low,high;order=5) = filter_helper(x,low,high,Bandstop,order)
bandstop(low,high;order=5) = x -> bandstop(x,low,high,order=order)

"""
    lowpass([x],low,[order=5])

Low-pass filter the sound at the specified frequency.

Filtering uses a butterworth filter of the given order.

With one positional argument this returns a function f(x)
that calls lowpass(x,low;order=order)
"""
lowpass(x,low;order=5) = filter_helper(x,low,0,Lowpass,order)
lowpass(low;order=5) = x -> lowpass(x,low,order=order)

"""
    highpass([x],high,[order=5])

High-pass filter the sound at the specified frequency.

Filtering uses a butterworth filter of the given order.

With one positional argument this returns a function f(x)
that calls highpass(x,high;order=order)
"""
highpass(x,high;order=5) = filter_helper(x,0,high,Highpass,order)
highpass(high;order=5) = x -> highpass(x,high,order=order)

function buildfilt(samplerate,low,high,kind)
  if kind == Bandpass
	  Bandpass(float(ustrip(inHz(low))),float(ustrip(inHz(high))),fs=samplerate)
  elseif kind == Lowpass
    Lowpass(float(ustrip(inHz(low))),fs=samplerate)
  elseif kind == Highpass
    Highpass(float(ustrip(inHz(high))),fs=samplerate)
  elseif kind == Bandstop
    Bandstop(float(ustrip(inHz(low))),float(ustrip(inHz(high))),fs=samplerate)
  end
end

function filter_helper(sound,low,high,kind,order)
  ftype = buildfilt(ustrip(samplerate(sound)),low,high,kind)
  f = digitalfilter(ftype,Butterworth(order))
  Sound(DSP.filt(f,sound),rate=samplerate(sound))
end

"""
    ramp([sound],[length=5ms])

Applies a half cosine ramp to start and end of the sound.

Ramps prevent clicks at the start and end of sounds.

When there is no sound argument, this returns a function f(x) that
ramp(x,length)
"""
function ramp(x::AbstractArray,length=5ms)
  ramp_n = inframes(length,samplerate(x))
	if nframes(x) < 2ramp_n
    error("Cannot apply two $(rounded_time(length,samplerate(x))) ramps to ",
          "$(rounded_time(duration(x),samplerate(x))) sound.")
  end

  n = nframes(x)
	r = Sound(n*frames,asseconds=false,rate=samplerate(x)) do t
    ifelse.(t .< ramp_n,
      -0.5.*cos.(π.*t./ramp_n).+0.5,
    ifelse.(t .< n .- ramp_n,
      1.0,
      -0.5.*cos.(π.*(t .- n .+ ramp_n)./ramp_n.+π).+0.5))
	end
	envelope(x,r)
end
ramp(length=5ms) = x -> ramp(x,length)

"""
    rampon([sound],[length=5ms])

Applies a half consine ramp to start of the sound.

When passed no sound argument, this returns a function f(x)
which calls rampon(x,length)
"""
function rampon(x::AbstractArray,length=5ms)
  ramp_n = inframes(length,samplerate(x))
	if nframes(x) < ramp_n
    error("Cannot apply a $(rounded_time(length,samplerate(x))) ramp to ",
          "$(rounded_time(duration(x),samplerate(x))) sound.")
  end

	r = Sound(ramp_n*frames,asseconds=false,rate=samplerate(x)) do t
    -0.5.*cos.(π.*t./ramp_n).+0.5
	end
	envelope(x,r)
end
rampon(length=5ms) = x -> rampon(x,length)

"""
    rampoff([sound],[length=5ms])

Applies a half consine ramp to the end of the sound.

When passed no sound argument, this returns a function f(x)
which calls rampoff(x,length)
"""
function rampoff(x::AbstractArray,length=5ms)
  len_s = inframes(length,samplerate(x))
  after = nframes(x)

  R = ustrip(samplerate(x))
  if len_s > nframes(x)
    error("Cannot apply $(rounded_time(len_s/R,R)) ramp to",
          " $(rounded_time(duration(x),R)) of audio.")
  end

  rampstart = (after - len_s)
	r = Sound(after*frames,asseconds=false,rate=samplerate(x)) do t
    ifelse.(t .< rampstart,1.0,-0.5.*cos.(π.*(t.-rampstart)./len_s.+π).+0.5)
	end
	envelope(x,r)
end
rampoff(length=5ms) = x -> rampon(x,length)

"""
    fadeto([a],[b],[transition=50ms])

A smooth transition from a to b, overlapping the end of one
with the start of the other by `transition`.

When passed a single sound, `a`, this returns a function `f(b)` that
calls fadeto(a,b,transition).
"""
function fadeto(a::AbstractArray,b::AbstractArray,transition=50ms)
  @assert(samplerate(a) == samplerate(b),
          "Sounds must have the same sample rate.")
  mix(rampoff(a,transition),
      [silence(duration(a) - transition,rate=samplerate(a));
       rampon(b,transition)])
end
fadeto(b::AbstractArray,transition=50ms) = a -> fadeto(a,b,transition)

"""
    amplify([x],ratio)

Amplify (positive) or attenuate (negative) the sound by a given ratio, typically
specified in decibels (e.g. amplify(x,10dB)).

*Note*: you can also directly multiply by a factor, e.g. x * 10dB,
which has the same effect as this function.

With one positional argument, `ratio`, this returns a function f(x) which calls
amplify(x,ratio).
"""
amplify(x,ratio) = x*ratio
amplify(ratio) = x -> amplify(x,ratio)

"""
    normpower(x)

Normalize the sound so it has a power of 1.
"""
normpower(x) = x ./ sqrt.(mean(x.^2,1))
normpower(x::Sound{R,T,C}) where {R,T,C} = Sound(R,C,x ./ sqrt.(mean(x.^2,1)))


text:
struct GenericVariableBridge{T, S} <: AbstractVariableBridge
    gram_matrix::Vector{MOI.VariableIndex}
    gram_constraint::MOI.ConstraintIndex{MOI.VectorOfVariables, S}
end

function add_variable_bridge(::Type{GenericVariableBridge{T, S}},
                             model::MOI.ModelLike, set::S) where {T, S}
    gram_matrix = MOI.add_variables(model, MOI.dimension(set))
    ci = MOI.add_constraint(model, MOI.VectorOfVariables(gram_matrix), set)
    func = MOI.SingleVariable[MOI.SingleVariable(v) for v in gram_matrix]
    # Need to specify `T` in the constructor as it cannot be deduced
    return func, GenericVariableBridge{T, S}(gram_matrix, ci)
end

function MOIB.added_constraint_types(::Type{GenericVariableBridge{T, S}}) where {T, S}
    return [(MOI.VectorOfVariables, S)]
end

function variable_bridge_type(S::Type{<:MOI.AbstractVectorSet}, T::Type)
    return GenericVariableBridge{T, S}
end

# Attributes, VariableBridge acting as an model
function MOI.get(bridge::GenericVariableBridge, ::MOI.NumberOfVariables)
    return length(bridge.gram_matrix)
end
function MOI.get(bridge::GenericVariableBridge{T, S},
                 ::MOI.NumberOfConstraints{MOI.VectorOfVariables, S}) where {T, S}
    return bridge.gram_constraint isa MOI.ConstraintIndex{MOI.VectorOfVariables, S} ? 1 : 0
end
function MOI.get(bridge::GenericVariableBridge{T, S},
                 ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables, S}) where {T, S}
    return [bridge.gram_constraint]
end

function MOI.delete(model::MOI.ModelLike, bridge::GenericVariableBridge)
    # First delete the constraint in which the Gram matrix appears
    MOI.delete(model, bridge.gram_constraint)
    # Now we delete the Gram matrix
    for variable in bridge.gram_matrix
        MOI.delete(model, variable)
    end
end

function MOI.get(model::MOI.ModelLike, ::MomentMatrixAttribute,
                 bridge::GenericVariableBridge)
    return MOI.get(model, MOI.ConstraintDual(), bridge.gram_constraint)
end

function MOI.get(model::MOI.ModelLike, ::GramMatrixAttribute,
                 bridge::GenericVariableBridge)
    return MOI.get(model, MOI.VariablePrimal(), bridge.gram_matrix)
end


text:
using DelimitedFiles

function ReadData(fname::AbstractString)
    data = readdlm(fname, '|', String, '\n')

    entries = []
    for i = 1:size(data)[1]
        push!(entries, (split(data[i,1], " ", keepempty=false), split(data[i,2], " ", keepempty=false)))
    end

    return entries
end

# 1 :   c  f   -> 2
# 7 : a c  f   -> 3
# 4 :  bcd f   -> 4
# 2 : a cde g  -> 5
# 3 : a cd fg  -> 5
# 5 : ab d fg  -> 5
# 0 : abc efg  -> 6
# 6 : ab defg  -> 6
# 9 : abcd fg  -> 6
# 8 : abcdefg  -> 7

# intersection of 0, 6, 9, 1 = "f"
# 0 : abc efg  -> 6
# 1 :   c  f   -> 2
# 6 : ab defg  -> 6
# 9 : abcd fg  -> 6

# difference between 7 and 1 = "a"
# 7 : a c  f   -> 3
# 1 :   c  f   -> 2

# difference between 1 and "f" = "c"

# intersection of all length >= 5 codes  = "ag"
# 2 : a cde g  -> 5
# 3 : a cd fg  -> 5
# 5 : ab d fg  -> 5
# 0 : abc efg  -> 6
# 6 : ab defg  -> 6
# 9 : abcd fg  -> 6
# 8 : abcdefg  -> 7

# difference between "ag" and "a" = g

# intersection of all length 5 codes with 4 = "d"
# 4 :  bcd f   -> 4
# 2 : a cde g  -> 5
# 3 : a cd fg  -> 5
# 5 : ab d fg  -> 5

# intersection of all length 6 codes with 4 = "bf"
# 4 :  bcd f   -> 4
# 0 : abc efg  -> 6
# 6 : ab defg  -> 6
# 9 : abcd fg  -> 6

# difference between "bf" and "f" = b

# e is just the other character

digitLengths = Dict([(0,6), (1,2), (2,5), (3,5), (4,4), (5,5), (6,6), (7,3), (8,7), (9,6)])

function SortString(s::AbstractString)
    return join(sort(split(s, "")), "")
end

function DecodeSignal(signal::Array)
    # dictionary of the coded digits
    code = Dict([(i,"") for i in [1,4,7,8]])

    for i in [1,4,7,8]
        code[i] = signal[findfirst(x -> length(x)==digitLengths[i], signal)]
    end

    charCode = Dict()

    charCode["a"] = StringSetDiff([code[7], code[1]])[1]
    charCode["f"] = StringIntersection(push!(filter(x -> length(x) == 6, signal), code[1]))[1]
    charCode["c"] = StringSetDiff([code[1], charCode["f"]])[1]

    # temp is the code for "ag"
    temp = join(StringIntersection(filter(x -> length(x) >= 5, signal)), "")
    charCode["g"] = StringSetDiff([temp, charCode["a"]])[1]

    charCode["d"] = StringIntersection(push!(filter(x -> length(x) == 5, signal), code[4]))[1]

    # temp is the code for "bf"
    temp = join(StringIntersection(push!(filter(x -> length(x) == 6, signal), code[4])), "")
    charCode["b"] = StringSetDiff([temp, charCode["f"]])[1]

    charCode["e"] = StringSetDiff(["abcdefg", join(values(charCode), "")])[1]

    code[0] = join([charCode[i] for i in ["a","b","c","e","f","g"]], "")
    code[2] = join([charCode[i] for i in ["a","c","d","e","g"]], "")
    code[3] = join([charCode[i] for i in ["a","c","d","f","g"]], "")
    code[5] = join([charCode[i] for i in ["a","b","d","f","g"]], "")
    code[6] = join([charCode[i] for i in ["a","b","d","e","f","g"]], "")
    code[9] = join([charCode[i] for i in ["a","b","c","d","f","g"]], "")

    # decoded digits
    decoded = Dict([(SortString(code[i]),i) for i=0:9])

    return decoded
end

# get the union of an array of strings
function StringUnion(strings::Array)
    return map(String, union([split(s, "") for s in strings]...))
end

# get the intersection of an array of strings
function StringIntersection(strings::Array)
    return map(String, intersect([split(s, "") for s in strings]...))
end

# get the set difference of an array of strings
function StringSetDiff(strings::Array)
    return map(String, setdiff([split(s, "") for s in strings]...))
end

function Part1(entries::AbstractArray)
    counter = 0
    for e in entries
        counter += count(x -> length(x) in [2,3,4,7], e[2])
    end

    return counter
end

function Part2(entries::AbstractArray)
    total = 0

    for e in entries
        decoded = DecodeSignal(e[1])

        for i=1:4
            total += decoded[SortString(String(e[2][i]))]*10^(4-i)
        end
    end

    return total
end

#data = ReadData("../data/test8.txt")
data = ReadData("../data/day8.txt")

println(Part1(data))
println(Part2(data))


text:
"""
    analytical_coefficients(setup, coeffs)

Get coefficients for the analytical module.
"""
function analytical_coefficients(setup, coeffs)
    rmean = (setup.rmin + setup.rmax) / 2
    include_in = setup.include_in
    include_ecs = setup.ecs_shape != :no_ecs
    dim = radial_dimension(setup)
    n = unique([1:include_in+include_ecs; length(coeffs.κ)])

    # Get OUT parameters
    r_out = (setup.rmin + setup.rmax) / 2

    # Get IN parameters
    if setup.include_in
        r_in = setup.in_ratio * r_out
    else
        r_in = []
    end

    # Get ECS parameters
    if include_ecs
        # Include ECS
        r_ecs = r_out + setup.ecs_ratio * rmean
    else
        # Empty ECS
        r_ecs = []
    end

    ρ = coeffs.ρ
    r = vcat(r_in, r_out, r_ecs) * 1e-6
    D = tr.(coeffs.D) ./ 3 .* 1e-6
    T₂ = coeffs.T₂ * 1e-6
    W = coeffs.κ[n]
    γ = coeffs.γ

    (; ρ, r, D, W, T₂, γ, dim)
end


text:
using Test
using Spandex

@testset "range 1" begin
    local u = SegmentTree{Int64}([1], max, typemin(Int64))

    @test 1 == top(u)
    @test 1 == range(u, 1, 1)
    @test 1 == u[1]
end

@testset "range 2" begin
    local u =
        SegmentTree{Int64}([3, 8, 6, 4, 2, 5, 9, 0, 7, 1], max, typemin(Int64))

    @test 9 == top(u)
    @test 9 == range(u, 6, 10)
    @test 9 == range(u, 1, 10)
    @test 8 == range(u, 1, 5)
    @test 8 == range(u, 2, 2)
    @test 8 == range(u, 2, 5)
    @test 7 == range(u, 8, 9)
    @test 5 == range(u, 5, 6)
    @test 6 == range(u, 3, 6)
end

@testset "range 3" begin
    local u = SegmentTree{Int64}([5, 4, 3, 2, 1], max, typemin(Int64))

    @test 5 == top(u)
    @test 3 == range(u, 3, 5)
    @test 5 == u[1]
    @test 4 == u[2]
    @test 3 == u[3]
end

@testset "set! 1" begin
    local values = [3, 8, 6, 4, 2, 5, 9, 0, 7, 1]
    local u = SegmentTree{Int64}(values, max, typemin(Int64))
    local v = SegmentTree{Int64}(10, max, typemin(Int64))

    @test !equals(u, v)

    for i = 1:10
        v[i] = values[i]
    end

    @test equals(u, v)
end

@testset "set! 2" begin
    local u =
        SegmentTree{Int64}([3, 8, 6, 4, 2, 5, 9, 0, 7, 1], max, typemin(Int64))

    u[7] = 0

    @test 8 == top(u)
    @test 8 == range(u, 1, 10)
    @test 7 == range(u, 6, 10)
    @test 8 == range(u, 1, 5)
    @test 8 == u[2]
    @test 8 == range(u, 2, 5)
    @test 7 == range(u, 8, 9)
    @test 5 == range(u, 5, 6)
    @test 6 == range(u, 3, 6)
end


text:
# Seed random generator to avoid random effects in testing
srand(1)

X = [1 1; 2 2;  1 -1]
y = [-1; -1; 1]

@test_approx_eq_eps optimize(RidgeReg(X, y), optimizer=:closed_form).theta [-0.0666666666666667,-0.23333333333333334] 1e-5
@test_approx_eq_eps optimize(RidgeReg(X, y, λ=1.0), optimizer=:closed_form).theta [0.06060606060606055,-0.606060606060606] 1e-5
@test_approx_eq_eps optimize(RidgeReg(X, y, λ=10.0), optimizer=:closed_form).theta [0.1791607732201792,-0.7732201791607733] 1e-5
@test_approx_eq_eps optimize(RidgeReg(X, y), optimizer=:sgd).theta [-0.0666666666666667,-0.23333333333333334] 5e-2
@test_approx_eq_eps optimize(RidgeReg(X, y, λ=1.0), optimizer=:sgd).theta [0.06060606060606055,-0.606060606060606] 5e-2
@test_approx_eq_eps optimize(RidgeReg(X, y, λ=10.0), optimizer=:sgd).theta [0.1791607732201792,-0.7732201791607733] 5e-2
@test_approx_eq_eps optimize(RidgeReg(X, y), optimizer=:l_bfgs).theta [-0.0666666666666667,-0.23333333333333334] 1e-5
@test_approx_eq_eps optimize(RidgeReg(X, y, λ=1.0), optimizer=:l_bfgs).theta [0.06060606060606055,-0.606060606060606] 1e-5
@test_approx_eq_eps optimize(RidgeReg(X, y, λ=10.0), optimizer=:l_bfgs).theta [0.1791607732201792,-0.7732201791607733] 1e-5
show(IOBuffer(), RidgeReg(X, y))

@test_throws DimensionMismatch RidgeReg(X', y)
@test_throws ArgumentError optimize(RidgeReg(X, y), optimizer=:blubs)
@test_throws ArgumentError optimize(RidgeReg(X, y, λ=0.0))


text:
# Translating
# https://proplot.readthedocs.io/en/latest/subplots.html#Automatic-subplot-spacing
# to Julia

# Import proplot using PyCall
using PyCall
proplot = pyimport("proplot")

f, axs = proplot.subplots(
    ncols=4, nrows=3, wspace=(0, 0, nothing), hspace=(0, nothing), # auto spacing not working AFAIK
    bottom="5em", right="5em", span=false,
    axwidth=1.1,
)

# Must reshape axs (and transpose it...)
axs = permutedims(reshape(axs, 4, 3), [2, 1])

[ax.format(
    suptitle="Tight layout with user overrides",
    rowlabels=["Row 1", "Row 2", "Row 3"],
    collabels=["Column 1", "Column 2", "Column 3", "Column 4"]
) for ax in axs]
[ax.format(xtickloc="top") for ax in axs[1, :]]
[ax.format(xtickloc="both") for ax in axs[3, :]]
[ax.format(ytickloc="neither") for ax in axs[:, 2]]
[ax.format(ytickloc="right") for ax in axs[:, 3]]
[ax.format(ytickloc="both") for ax in axs[:, 4]]
[ax.format(title="Title\nTitle\nTitle", xlabel="xlabel") for ax in axs[end, :]]
[ax.format(ylabel="ylabel\nylabel") for ax in axs[:, 1]]
f.savefig("figs/Automatic-subplot-spacing-2.svg", transparent=false)

max_stars_repo_name: natgeo-wong/JuliaProPlotExamples

max_issues_repo_name: natgeo-wong/JuliaProPlotExamples

max_forks_repo_name: natgeo-wong/JuliaProPlotExamples


text:
@testset "Violations" begin

    # Test dynamics violation
    T = Float64
    N = 10
    dt = 0.1
    i = 1
    model = UnicycleGame(p=3)
    probsize = ProblemSize(N, model)
    pi = probsize.pz[i]
    pdtraj = PrimalDualTraj(probsize, dt)

    init_traj!(pdtraj, x0=zeros(SVector{model.n,T}), f=zeros)
    dyn_vio = dynamics_violation(model, pdtraj)
    @test dyn_vio.N == pdtraj.probsize.N
    @test dyn_vio.vio == zeros(N-1)
    @test dyn_vio.max == 0.0
    init_traj!(pdtraj, x0=ones(SVector{model.n,T}), f=ones, amplitude=1.0)
    @test dynamics_violation(model, pdtraj).max - maximum(abs.(dynamics_residual(model, pdtraj, 1))) < 1e-10
    @test dynamics_violation(model, pdtraj, i).max - maximum(abs.(dynamics_residual(model, pdtraj, 1)[pi])) < 1e-10

    # Test control violation
    game_con = GameConstraintValues(probsize)
    u_max =  0.1*ones(model.m)
    u_min = -0.1*ones(model.m)
    add_control_bound!(game_con, u_max, u_min)
    init_traj!(pdtraj, x0=zeros(SVector{model.n,T}), f=ones, amplitude=1.0)
    con_vio = control_violation(game_con, pdtraj)
    @test con_vio.N == pdtraj.probsize.N
    @test con_vio.vio == 0.9*ones(N-1)
    @test con_vio.max == 0.9
    con_vio = control_violation(game_con, pdtraj, i)
    @test con_vio.N == pdtraj.probsize.N
    @test con_vio.vio == 0.9*ones(N-1)
    @test con_vio.max == 0.9

    # Test state violation
    game_con = GameConstraintValues(probsize)
    walls = [Wall([0.,1], [1,0], [1,1]/sqrt(2))]
    add_wall_constraint!(game_con, walls)
    init_traj!(pdtraj, x0=zeros(SVector{model.n,T}), f=ones, amplitude=1.0)
    sta_vio = state_violation(game_con, pdtraj)
    @test sta_vio.N == pdtraj.probsize.N
    @test norm(sta_vio.vio - [0; sqrt(2)/2*ones(N-1)...], 1) <= 1e-10
    @test norm(sta_vio.max - sqrt(2)/2,1) < 1e-10
    sta_vio = state_violation(game_con, pdtraj, i)
    @test sta_vio.N == pdtraj.probsize.N
    @test norm(sta_vio.vio - [0; sqrt(2)/2*ones(N-1)...], 1) <= 1e-10
    @test norm(sta_vio.max - sqrt(2)/2,1) < 1e-10

    # Test optimality violation
    core = NewtonCore(probsize)
    opt_vio = optimality_violation(core)
    @test opt_vio.N == core.probsize.N
    @test norm(opt_vio.vio - zeros(N), 1) <= 1e-10
    @test norm(opt_vio.max - 0.0,1) < 1e-10
    opt_vio = optimality_violation(core, i)
    @test opt_vio.N == core.probsize.N
    @test norm(opt_vio.vio - zeros(N), 1) <= 1e-10
    @test norm(opt_vio.max - 0.0,1) < 1e-10

    stamp = stampify(:opt, 2, :u, 2, 5)
    add2sub(core.res_sub[stamp], 1e2*ones(model.mi[2]))
    opt_vio = optimality_violation(core)
    vio = zeros(N)
    vio[5] = 1e2
    @test norm(opt_vio.vio - vio, 1) <= 1e-10
    @test norm(opt_vio.max - 1e2,1) < 1e-10

end


text:
"""
    Design{T1<:Integer, T2<:Real, PType<:Parameters}

Binary two-stage design object defining a discrete sampe size function
``n(x_1)`` and ``c(x_1)``, where ``x_1`` is the observed number of stage-one events
and the test rejects the superiority null hypothesis whenever ``x_1 + x_2 > c(x_1)``.
"""
mutable struct Design{T1<:Integer, T2<:Real, PType<:Parameters}

  n::Vector{T1}
  c::Vector{T2} # rejected if x1 + x2 > c(x1), must allow real to hold Infinity
  params::PType

  function Design{T1,T2,PType}(
    n::Vector{T1}, c::Vector{T2}, params::PType
  ) where {T1<:Integer, T2<:Real, PType<:Parameters}

    any(n .< (length(n) - 1)) ?
      error("n must always be larger than n1 (= length(n) - 1") : nothing
    length(n) != length(c) ?
      error("n and c must be of equal length") : nothing
    new(n, c, params)

  end

end # Design

function Design(
  n::Vector{T1}, c::Vector{T2}, params::PType
) where {T1<:Integer,T2<:Real,PType<:Parameters}

  return Design{T1,T2,PType}(n, c, params)

end

"""
    Design(n, c, params::PType) where {PType<:Parameters}

Construct a [`Design`](@ref) object from `n` and `c` directly.
Here, `n` must be convertable to an integer vector and `c` to a real vector
(both of same length ``n_1``).
Then, `n[x1 + 1]` is the final sample size after observing `x1` responses in
stage one and `c[x1 + 1]` the corresponding critical value.
A constructor without explicitly specifying a set of parameters with respect to
which thee design is optimal also exists.

# Parameters

| Parameter  | Description |
| ---------: | :---------- |
| n          | sample size function (must be convertable to integer vector) |
| c          | critical value function (must be convertable to real vector) |
| params     | parameter object |

"""
function Design(n, c, params::PType) where {PType<:Parameters}

  n = convert(Vector{Int}, n)
  c = convert(Vector{Float64}, c)
  return Design{Int,Float64,Ptype}(n, c, params)

end

"""
    Design(n, c)

Construct a [`Design`](@ref) object from `n` and `c` directly.
Here, `n` must be convertable to an integer vector and `c` to a real vector
(both of same length ``n_1``).

# Parameters

| Parameter  | Description |
| ---------: | :---------- |
| n          | sample size function (must be convertable to integer vector) |
| c          | critical value function (must be convertable to real vector) |

"""
function Design(n, c)

  n = convert(Vector{Int}, n)
  c = convert(Vector{Float64}, c)
  return Design{Int,Float64,NoParameters}(n, c, NoParameters())

end

"""
    convert(::Type{DataFrames.DataFrame}, design::Design)

Convert a design to a `DataFrame` with column names `x1`, `n`, and `c`.
"""
convert(::Type{DataFrames.DataFrame}, design::Design) =
    DataFrames.DataFrame(
        x1 = 0:interimsamplesize(design),
        n  = samplesize(design),
        c  = criticalvalue(design)
    )

# these enable array broadcasting of all methods!
Base.size(::Design) = ()

Base.length(design::Design) = 1

Base.start(design::Design) = true

Base.done(design::Design, state) = true

Base.getindex(design::Design, i) = design

Base.show(io::IO, design::Design) = print("Design")

"""
    interimsamplesize(design::Design)

Return interim sample size ``n_1`` of a design.
"""
interimsamplesize(design::Design) = length(design.n) - 1

"""
    parameters(design::Design)

Return stored parameter object of a design.
"""
parameters(design::Design) = design.params

"""
    samplesize(design::Design)

Return vector of final sample sizes.
"""
samplesize(design::Design) = design.n

"""
    samplesize(design::Design, x1::T) where {T<:Integer}

Return final sample size for `x1`.
"""
samplesize(design::Design, x1::T) where {T<:Integer} =
    (checkx1(x1, design); design.n[x1 + 1])

"""
    criticalvalue(design::Design)

Return vector of critical values.
"""
criticalvalue(design::Design) = design.c

"""
    criticalvalue(design::Design, x1::T) where {T<:Integer}

Return critical value for `x1`.
"""
criticalvalue(design::Design, x1::T) where {T<:Integer} =
    (checkx1(x1, design); design.c[x1 + 1])

"""
    pdf{T1<:Integer, T2<:Real}(design::Design, x1::T1, x2::T1, p::T2)

Probability density function of ``(x1, x2)`` responses in stage one and two,
respectively, under `design` given response rate ``p``.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| design       | a Design |
| x1           | number of stage-one responses |
| x2           | number of stage-two responses |
| p            | response probability |

"""
function pdf(
  design::Design,
  x1::T1, x2::T1, p::T2
) where {T1<:Integer, T2<:Real}

  @checkprob p
  try
    checkx1x2(x1, x2, design)
  catch
    return 0.0
  end
  n1 = interimsamplesize(design)
  n  = samplesize(design, x1)
  try
    return p^(x1 + x2)*(1 - p)^(n - x1 - x2)*binomial(n1, x1)*binomial(n - n1, x2) # is 0 if impossible
  catch
    return p^(x1 + x2)*(1 - p)^(n - x1 - x2)*binomial(BigInt(n1), BigInt(x1))*binomial(BigInt(n - n1), BigInt(x2)) # is 0 if impossible
  end

end

"""
    power(design::Design, x1::T1, p::T2) where {T1<:Integer, T2<:Real}

conditional power of a design for a given response rate ``p`` and
the observed number of stage-one responses `x1`.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| design       | a [`Design`](@ref) |
| p            | response probability |
| x1           | number of stage-one responses |

"""
function power(
  design::Design, x1::T1, p::T2
) where {T1<:Integer, T2<:Real}

  @checkprob p
  checkx1(x1, design)
  res = _cpr(
    x1,
    interimsamplesize(design),
    samplesize(design, x1),
    criticalvalue(design, x1),
    p
  )
  return min(1, max(0, res)) # guarantee bounds!

end

"""
    power(design::Design, p::T) where {T<:Real}

Power of a design for a given response rate ``p`` .

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| design       | a [`Design`](@ref) |
| p            | response probability |

"""
function power(design::Design, p::T) where {T<:Real}

  @checkprob p
  n1      = interimsamplesize(design)
  X1      = Distributions.Binomial(n1, p) # stage one responses
  x1range = collect(0:n1)
  return min(1, max(0, vecdot(Distributions.pdf.(X1, x1range), power.(design, x1range, p))))

end

"""
    expectedpower(
      design::Design, x1::T, prior::Function; mcrv::Real = mcrv(parameters(design))
    ) where {T<:Integer}

Compute the conditional expected power of a design given `x1`.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| design       | a Design |
| x1           | stage one number of responses |
| prior        | prior function prior(p) for response probability p |
| mcrv         | minimal clinically relevant value for expected power calculation |

"""
function expectedpower(
    design::Design, x1::T, prior::Function; mcrv::Real = mcrv(parameters(design))
) where {T<:Integer}

  checkx1(x1, design)
  z   = QuadGK.quadgk(
      p -> prior(p),  # f(p)
      mcrv,           # p_min
      1,              # p_max
      atol    = 0.001  # tolerance
  )[1]
  res = QuadGK.quadgk(
      p -> prior(p)*power(design, x1, p)/z, # f(p)
      mcrv,         # p_min
      1,             # p_max
      atol    = 0.001 # tolerance
  )[1]
  return min(1, max(0, res)) # guarantee bounds!

end

"""
    expectedpower(design::Design, prior::Function; mcrv::Real = mcrv(parameters(design)))

Compute the expected power of a given design.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| design       | a Design |
| prior        | prior function prior(p) for response probability p |

"""
function expectedpower(design::Design, prior::Function; mcrv::Real = BinaryTwoStageDesigns.mcrv(parameters(design)) )

  z   = QuadGK.quadgk(
      p -> prior(p),             # f(p)
      mcrv,  # p_min
      1,                         # p_max
      atol    = 0.001             # tolerance
  )[1]
  res = QuadGK.quadgk(
      p -> prior(p)*power(design, p)/z, # f(p)
      mcrv,    # p_min
      1,                           # p_max
      atol    = 0.001               # tolerance
  )[1]
  return min(1, max(0, res)) # guarantee bounds!

end

"""
    stoppingforfutility{T<:Real}(design::Design, p::T) where {T<:Real}

Compute probability of stopping early for futility of a given design.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| design       | a [`Design`](@ref) |
| p            | response probability |

"""
function stoppingforfutility(design::Design, p::T) where {T<:Real}

  @checkprob p
  n1  = interimsamplesize(design)
  X1  = Distributions.Binomial(n1, p) # stage one responses
  res = 0.0
  c   = criticalvalue(design)
  for x1 in 0:n1
    if c[x1 + 1] == Inf
      res += Distributions.pdf(X1, x1)
    end
  end
  return res

end

"""
    stoppingforefficacy{T<:Real}(design::Design, p::T) where {T<:Real}

Compute probability of stopping early for futility of a given design.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| design       | a [`Design`](@ref) |
| p            | response probability |

"""
function stoppingforefficacy(design::Design, p::T) where {T<:Real}

  @checkprob p
  n1  = interimsamplesize(design)
  X1  = Distributions.Binomial(n1, p) # stage one responses
  res = 0.0
  c   = criticalvalue(design)
  for x1 in 0:n1
    if c[x1 + 1] == -Inf
      res += Distributions.pdf(X1, x1)
    end
  end
  return res

end

"""
    test(design::Design, x1::T, x2::T)::Bool where {T<:Integer}

Binary test decision of `design` when observing `x1` responses in stage one and
`x2` responses in stage two.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| design       | a [`Design`](@ref) |
| x1           | number of stage-one responses |
| x2           | number of stage-two responses |

"""
function test(design::Design, x1::T, x2::T)::Bool where {T<:Integer}

  checkx1x2(x1, x2, design)
  return x1 + x2 > criticalvalue(design, x1)

end

"""
    simulate(design::Design, p::T2, nsim::T1) where {T1<:Integer, T2<:Real}

Simulate `nsim` trials of design `design` with true response probability `p`.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| design       | a [`Design`](@ref) |
| p            | true response probability |
| nsim         | number of simulated trials |

# Return Value

A `DataFrame` with columns `x1` (stage one responses) `n` (overall sample size)
`c` (critical value), `x2` (stage-two responses), and `rejectedH0` (test decision).

"""
function simulate(design::Design, p::T2, nsim::T1) where {T1<:Integer, T2<:Real}

  x2    = SharedArray{Int}(nsim)
  n     = SharedArray{Int}(nsim)
  c     = SharedArray{Float64}(nsim)
  rej   = SharedArray{Bool}(nsim)
  n1    = interimsamplesize(design)
  rv_x1 = Distributions.Binomial(n1, p)
  x1    = SharedArray{Int}(nsim)
  @sync @parallel for i in 1:nsim
      x1[i]  = rand(rv_x1)
      n2     = samplesize(design, x1[i]) - n1
      n[i]   = n2 + n1
      c[i]   = criticalvalue(design, x1[i])
      x2[i]  = rand(Distributions.Binomial(n2, p))
      rej[i] = test(design, x1[i], x2[i])
  end
  return DataFrames.DataFrame(
      x1 = convert(Vector{Int}, x1),
      n  = convert(Vector{Int}, n),
      c  = convert(Vector{Float64}, c),
      x2 = convert(Vector{Int}, x2),
      rejectedH0 = convert(Vector{Bool}, rej)
  )

end

"""
    score(design::Design, params::Parameters)::Real

Evaluates the score of a design for a given set of parameters.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| design       | a [`Design`](@ref) |
| params       | a parameters object holding the required information about the score |

# Details

A method `score(design::Design)::Real` is also implemented which
uses the parameter object stored within the design object after optimization.
Note that this is only possible if the design was created as result of calling
`getoptimaldesign()`.

"""
score(design::Design, params::Parameters)::Real = error("not implemented")

score(design::Design)::Real = score(design, parameters(design))

"""
    jeffreysprior(design::Design)

Computes the Jeffreys prior of any given design.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| design       | a Design |

# Return Value

A function with signature `prior{T<:Real}(p::T)::Real` where `p` is the response
probability and `prior(p)` the PDF of the Jeffres prior at `p`.

# Examples
```julia-repl
julia> ss = SimpleSampleSpace(10:25, 100, n2min = 5)
julia> params = SimpleMinimalExpectedSampleSize(ss, .2, .4, .05, .2, .4)
julia> design, res = getoptimaldesign(params, solver = Gurobi.GurobiSolver())
julia> f = jeffreysprior(design)
julia> f(.5)
```
"""
function jeffreysprior(design::Design)

  function sqrtfi(p::Float64)::Float64

    res  = 0.0
    n1   = interimsamplesize(design)
    nmax = maximum(samplesize(design))
    supp = support(design)
    for i in 1:size(supp, 1)
        x1, x2 = supp[i, 1], supp[i, 2]
        x   = x1 + x2
        n   = samplesize(design, x1)
        res += binomial(BigInt(n1), BigInt(x1))*binomial(BigInt(n - n1), BigInt(x2))*BigFloat(p^x*(1 - p)^(n - x)*(x/p - (n - x)/(1 - p))^2)
    end
    return sqrt(res)

  end

  z = QuadGK.quadgk(sqrtfi, 0, 1, atol    = 0.001)[1] # exact integration from 0 to 1 is expensive!

  function prior{T<:Real}(p::T)::Real

    @checkprob p
    return sqrtfi(p)/z

  end

  return prior

end

"""
    save(filename::String, design::Design)

Save design as .jls file.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| filename     | filename for sving design |
| design       | a [`Design`](@ref) |

"""
function save(filename::String, design::Design)

  file = open(filename, "w")
  serialize(file, design)
  close(file)

end

"""
    writecsv(filename::String, design::Design; label::String = "")

Save design as .csv file.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| filename     | filename for sving design |
| design       | a [`Design`](@ref) |
| label        | string name of design |
"""
function writecsv(filename::String, design::Design; label::String = "")

  df = DataFrames.DataFrame(design)
  label != "" ? df[:design] = label : nothing
  CSV.write(filename, df)
  return df

end

function writepropertiescsv(
  filename::String, design::Design;
  label::String = "", pvec = linspace(0, 1, 101), nmax::Int = Int(round(1.1*maximum(samplesize(design))))
)

  df = DataFrame(
    p      = Real[],
    n      = Integer[],
    PDF    = Real[],
    ESS    = Real[],
    Q1SS   = Real[],
    Q3SS   = Real[],
    Power  = Real[]
  )
  for p in pvec
    ss = SampleSize(design, p)
    for n in 0:nmax
      push!(df, [p n pdf(ss, n ) mean(ss) quantile(ss, .25) quantile(ss, .75) power(design, p)])
    end
  end
  label != "" ? df[:design] = label : nothing
  CSV.write(filename, df)
  return df

end

# utility functions
function checkx1(x1::T, design::Design) where {T<:Integer}

  x1 < 0 ? throw(InexactError("x1 must be non-negative")) : nothing
  x1 > interimsamplesize(design) ? throw(InexactError("x1 smaller or equal to n1")) : nothing

end

function checkx1x2(x1::T, x2::T, design::Design) where {T<:Integer}

  checkx1(x1, design)
  x2 < 0 ? throw(InexactError("x2 must be non-negative")) : nothing
  n1 = interimsamplesize(design)
  n2 = samplesize(design, x1) - n1
  x2 > n2 ? throw(InexactError("x2 must be smaller or equal to n2")) : nothing

end

function ispossible(design::Design, x1::T, x2::T) where {T<:Integer}

  res = x1 < 0 ? false : true
  res = x1 > interimsamplesize(design) ? false : true
  res = x2 < 0 ? false : true
  res = x2 > samplesize(design, x1) - interimsamplesize(design) ? false : true

end

function support(design::Design)

    n1     = interimsamplesize(design)
    nmax   = maximum(samplesize(design))
    return [[x1, x2] for x1 in 0:n1, x2 in 0:(nmax - n1) if
        (x2 <= samplesize(design, x1) - n1) & ((samplesize(design, x1) > n1) | (x2 == 0))
    ] |> x-> hcat(x...)'

end

function _cpr(x1, n1, n, c, p)
    # conditional probability to reject
    if x1 > c
        return(1.0)
    end
    if n - n1 <= c - x1
        return(0.0)
    end # TODO: simplify
    return 1 - Distributions.cdf(Distributions.Binomial(n - n1, p), convert(Int64, c - x1))
end

function _cprprior(x1, n1, n, c, p0, prior)
    if x1 > c
        return 1.0
    end
    if n - n1 + x1 <= c
        return 0.0
    end
    z = QuadGK.quadgk(p -> prior(p)*dbinom(c - x1, n - n1, p), p0, 1.0, atol = 1e-4)[1]
    cposterior(p) = prior(p)*dbinom(c - x1, n - n1, p)/z
    return QuadGK.quadgk(p -> cposterior(p)*_cpr(x1, n1, n, c, p), p0, 1.0, atol = 1e-4)[1]
end

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Frames (Modia3D/Frames/_module.jl)
#

"""
    const Modia3D.Vector3D = SVector{3,Float64}

Type of a vector in 3D space (e.g. position vector of the origin of a frame)
"""
const Vector3D = SVector{3,Float64}

"""
    const Modia3D.ZeroVector3D = Vector3D(0.0, 0.0, 0.0)

Constant of a Vector3D where all elements are zero
"""
const ZeroVector3D = Vector3D(0.0, 0.0, 0.0)

"""
    vec = Modia3D.axisValue(axis, value)
    vec = Modia3D.axisValue(axis, positive, value)

Return `vec::Vector3D` where all elements are zero with exception of `vec[axis] = value` or
`vec[axis] = positive ? value : -value`.
"""
@inline axisValue(axis::Int, value::Number)::Vector3D = axis==1 ? @SVector([value, 0.0, 0.0]) : (axis==2 ? @SVector([0.0, value, 0.0]) :
                                                       (axis==3 ? @SVector([0.0, 0.0, value]) : error("Bug in Modia3D: axisValue($axis, ...) - argument needs to be 1,2 or 3.")))
@inline axisValue(axis::Int, positive::Bool, value::Number)::Vector3D = positive ? axisValue(axis,value) : axisValue(axis,-value)

"""
    M = Modia3D.skew(e::AbstractVector)

Return the skew symmetric matrix `M::SMatrix{3,3,Float64,9}` of vector `e` (`length(e) = 3`)
"""
skew(e::AbstractVector) = @SMatrix([  0.0   -e[3]    e[2];
                                      e[3]   0.0    -e[1];
                                     -e[2]   e[1]    0.0 ])


text:
using EconPDEs, Distributions, LookupTables

mutable struct BansalYaronModel
    # consumption process parameters
    μbar::Float64
    vbar::Float64
    κμ::Float64
    νμ::Float64
    κv::Float64
    νv::Float64

    # utility parameters
    ρ::Float64
    γ::Float64
    ψ::Float64

    # lookup
    μn::Int64
    vn::Int64
    τ::LookupTable2D
end

function BansalYaronModel(;μbar = 0.018, vbar = 0.00073, κμ = 0.252, νμ = 0.528, κv = 0.156, νv = 0.00354, ρ = 0.024, γ = 7.5, ψ = 1.5, μn=30, vn=30, τ=LookupTable2D(30, 30))
    BansalYaronModel(μbar, vbar, κμ, νμ, κv, νv, ρ, γ, ψ, μn, vn, τ)
end

function initialize_stategrid(m::BansalYaronModel)
    μbar = m.μbar ; vbar = m.vbar ; κμ = m.κμ ; νμ = m.νμ ; κv = m.κv ; νv = m.νv ; ρ = m.ρ ; γ = m.γ ; ψ = m.ψ
    μn = m.μn ; vn = m.vn

    σ = sqrt(νμ^2 * vbar / (2 * κμ))
    μmin = quantile(Normal(μbar, σ), 0.025)
    μmax = quantile(Normal(μbar, σ), 0.975)
    μs = range(μmin, stop = μmax, length = μn)

    α = 2 * κv * vbar / νv^2
    β = νv^2 / (2 * κv)
    vmin = quantile(Gamma(α, β), 0.025)
    vmax = quantile(Gamma(α, β), 0.975)
    vs = range(vmin, stop = vmax, length = vn)
    OrderedDict(:μ => μs, :v => vs)
end

function initialize_y(m::BansalYaronModel, stategrid::OrderedDict)
    OrderedDict(:p => ones(length(stategrid[:μ]), length(stategrid[:v])))
end

function update_τ(m::BansalYaronModel, stategrid::OrderedDict, τmat::Array{Float64, 2})
    m.τ = LookupTable2D(stategrid[:μ], stategrid[:v], τmat)
    m.τ.xvar = stategrid[:μ]
    m.τ.yvar = stategrid[:v]
end

function (m::BansalYaronModel)(state::NamedTuple, y::NamedTuple)
    μbar = m.μbar ; vbar = m.vbar ; κμ = m.κμ ; νμ = m.νμ ; κv = m.κv ; νv = m.νv ; ρ = m.ρ ; γ = m.γ ; ψ = m.ψ ; τ = m.τ
    μ, v = state.μ, state.v
    p, pμ, pv, pμμ, pμv, pvv = y.p, y.pμ, y.pv, y.pμμ, y.pμv, y.pvv

    τvar = lookup(τ, μ, v)
    τnew = μ^2 - v * μ

    # drift and volatility of c, μ, σ, p
    μc = μ
    σc = sqrt(v)
    μμ = κμ * (μbar - μ)
    σμ = νμ * sqrt(v)
    μv = κv * (vbar - v)
    σv = νv * sqrt(v)
    σp_Zμ = pμ / p * σμ
    σp_Zv = pv / p * σv
    σp2 = σp_Zμ^2 + σp_Zv^2
    μp = pμ / p * μμ + pv / p * μv + 0.5 * pμμ / p * σμ^2 + 0.5 * pvv / p * σv^2

    # Market price of risk κ
    κ_Zc = γ * σc
    κ_Zμ = - (1 - γ * ψ) / (ψ - 1) * σp_Zμ
    κ_Zv = - (1 - γ * ψ) / (ψ - 1) * σp_Zv
    κ2 = κ_Zc^2 + κ_Zμ^2 + κ_Zv^2

    # Risk free rate r
    r = ρ + μc / ψ - (1 + 1 / ψ) / 2 * γ * σc^2 - (γ * ψ - 1) / (2 * (ψ - 1)) * σp2

    # Market Pricing
    # pt = p * (1 / p + μc + μp - r - κ_Zc * σc - κ_Zμ * σp_Zμ - κ_Zv * σp_Zv)
    pt = p * (1 / p - ρ + (1 - 1 / ψ) * (μc - 0.5 * γ * σc^2) + μp + 0.5 * (1 / ψ - γ) / (1 - 1 / ψ) * σp2)

    # return (pt,), (μμ, μv), (p = p, r = r, κ_Zc = κ_Zc, κ_Zμ = κ_Zμ, κ_Zv = κ_Zv, σμ = σμ, σμ_Zv = 0.0, σv_Zμ = 0.0, σv = σv, μ = μ, v = v, σμ2 = σμ^2, σv2 = σv^2, σμv = 0.0, μμ = μμ, μv = μv, σp2 = σp2, σp_Zμ = σp_Zμ, σp_Zv = σp_Zv)
    return (pt,), (μμ, μv), (p = p, μμ = μμ, σμ = σμ, μv = μv, σv = σv, τnew = τnew, μ = μ, v = v)
end

# Bansal Yaron (2004)
m = BansalYaronModel()
stategrid = initialize_stategrid(m)
y0 = initialize_y(m, stategrid)
update_τ(m, stategrid, zeros(m.μn, m.vn))
y, result, distance = pdesolve(m, stategrid, y0)

update_τ(m, stategrid, result[:τnew])


text:
"""
draw(ctx::CairoContext, board, w, h)

draw `board` on `ctx` with width `w` and height `h`
"""
function draw end

"""
defaultsize(::Type{Game})

returns the default size of the drawing board in pixels for `Game`
"""
function defaultsize end

defaultsize(::Type{<:Game}) = (500,500)

"""
correspondingmove(board, x, y)

take normalized coordinates (between 0 and 1) `x` and `y` and return the corresponding move

TODO: can return missing for no move?
"""
function correspondingmove end


text:
module FeaturesTest

    using FactCheck
    using DataFrames
    include(Pkg.dir("MethylUtils", "src","Features.jl"))

   # GenomePosition

    g1 = GenomePos("chr1","3","15")
    g2 = GenomePos("chr1",3,15)
    g3 = GenomePos(1,3,15)

    facts("GenomicPosition") do
        @fact g1.chr   => int16(1)
        @fact g1.start => int16(3)
        @fact g1.stop  => int16(15)
        @fact g2.chr   => int16(1)
        @fact g2.start => int16(3)
        @fact g2.stop  => int16(15)
        @fact g3.chr   => int16(1)
    end

   # CytosineCount ?

   # DO SOME TESTS

   # FeatureTable

   df = DataFrame(name = ["A","B","C","D"], chr = ["chr1","chr2","chr1", "chr3"], start = 0:3, stop = 0:3 )
   ft = FeatureTable( ["organism"=>"mm10"], df )

   # FeatureRegions
   fr = FeatureRegions(ft)
   println(fr)
   regions = getoverlaps(fr, "chr2", 2, 2)
   println(regions)

   # load ensGene rows into DF
   # names_str = [ "bin", "name", "chrom", "strand", "txStart", "txEnd", "cdsStart",
   #            "cdsEnd", "exonCount", "exonStarts", "exonEnds", "score", "name2",
   ##           "cdsStartStat", "cdsEndStat", "exonFrames" ]
   # names = map(symbol, names_str)
   #DF = readtable( Pkg.dir("MethylUtils","testdata","ensGene1000lines.txt.gz"),
   #                separator='\t',
   #                names=names
   #                )
   # make a feature table
   #FT = FeatureTable({"organism" => "superman" }, DF)
   #FR = FeatureRegions(FT)

   # SequenceFeatures
   SF = SequenceFeatures( Dict(), Dict())
   feature = CytosineCount( SequencePos(100,200),true,15,90)
   add_sequence_feature( SF, "chr1", feature  )
   feature = CytosineCount( SequencePos(105,205),true,45,90)
   add_sequence_feature( SF, "chr1", feature  )
   println(SF)
end


text:
using DelimitedFiles

function solve(input)
    claims = readdlm(input,'\n')
    fabric = fill([],(1000,1000))
    mapreduce(parseClaim, fillGrid, claims, init=fabric)
    # part 1
    surface = count(x->length(x)>1,fabric)
    println("square inches of fabric within 2 or more claims: $surface")
    # part 2
    mapreduce(parseClaim, checkClaim, claims, init=fabric)

end

function parseClaim(s)
    m = match(r"#(?<id>\d+) @ (?<x>\d+),(?<y>\d+): (?<w>\d+)x(?<h>\d+)", s)
    return tuple(map(parseInt,m.captures)...)
end

function parseInt(s)
    return tryparse(Int16,s)
end

function fillGrid(fabric,(id,x,y,w,h))
    for i in x+1:(x+w)
        for j in y+1:(y+h)
            if length(fabric[i,j])==0
                fabric[i,j]=[id]
            else
                push!(fabric[i,j],id)
            end
        end
    end
    return fabric
end

function checkClaim(fabric,(id,x,y,w,h))
    nonoverlapping=true
    for i in x+1:(x+w)
        for j in y+1:(y+h)
            nonoverlapping=nonoverlapping&&(length(fabric[i,j])==1)
        end
    end
    if nonoverlapping
        println("Found claim $id")
    end
    return fabric
end

solve("input.txt")


text:
using Random, Plots
include("plotting.jl")

rng = MersenneTwister(1)
x, f = random_probability_distribution(rng, 3, 5, -2.0, 3.0)
α = 0.25

@assert issorted(x)
@assert all(f .≥ 0)
@assert sum(f) ≈ 1

dir = "images"

@info "Distributions"
plt1 = plot_distributions(x, f, α)
savefig(plt1, joinpath(dir, "distributions.svg"))

@info "VaR"
plt2 = plot_VaR(x, f, α)

@info "CVaR"
plt3 = plot_CVaR(x, f)

plt4 = plot(plt2, plt3, layout=(1, 2), legend=false, size=(720, 480))
savefig(plt4, joinpath(dir, "conditional-value-at-risk.svg"))

max_stars_repo_name: jaantollander/ConditionalValueAtRisk

max_issues_repo_name: jaantollander/ConditionalValueAtRisk

max_forks_repo_name: jaantollander/ConditionalValueAtRisk