text:
# A special kind of container promise that performs table layout optimization.

mutable struct Table <: ContainerPromise
    # Direct children must be Contexts, and not just Containers. If
    # children[i,j] has a vector with multiple children it indicates multiple
    # possible layouts for that cell in the table.
    children::Matrix{Vector{Context}}

    # In the formulation of the table layout problem used here, we are trying
    # find a feasible solution in which the width + height of a particular
    # group of cells in the table is maximized.
    x_focus::UnitRange{Int}
    y_focus::UnitRange{Int}

    # If non-nothing, constrain the focused cells to have a proportional
    # relationship.
    x_prop::Union{Vector{Float64}, Nothing}
    y_prop::Union{Vector{Float64}, Nothing}

    # If non-nothing, constrain the focused cells to have a fixed aspect ratio.
    aspect_ratio::Union{Float64, Nothing}

    # fixed configuration
    fixed_configs::Vector

    # Coordinate system used for children
    units::Union{UnitBox, Nothing}

    # Z-order of this context relative to its siblings.
    order::Int

    # Ignore this context and everything under it if we are
    # not drawing to the javascript backend.
    withjs::Bool

    # Ignore this context if we are drawing to the SVGJS backend.
    withoutjs::Bool
end

function Table(m::Integer, n::Integer, y_focus::UnitRange{Int}, x_focus::UnitRange{Int};
               y_prop=nothing, x_prop=nothing,
               aspect_ratio=nothing,
               units=nothing, order=0, withjs=false, withoutjs=false,
               fixed_configs=Any[])

    if x_prop != nothing
        @assert length(x_prop) == length(x_focus)
        x_prop ./= sum(filter(x -> !isnan(x), x_prop))
    end

    if y_prop != nothing
        @assert length(y_prop) == length(y_focus)
        y_prop ./= sum(filter(x -> !isnan(x), y_prop))
    end

    tbl = Table(Array{Vector{Context}}(undef, (m, n)),
              x_focus, y_focus,
              x_prop, y_prop,
              aspect_ratio,
              fixed_configs,
              units, order, withjs, withoutjs)
    for i in 1:m, j in 1:n
        tbl.children[i, j] = Array{Context}(undef, 0)
    end
    return tbl
end

const table = Table

getindex(t::Table, i::Integer, j::Integer) = t.children[i, j]
setindex!(t::Table, child, i::Integer, j::Integer) = t.children[i, j] = child
size(t::Table, i::Integer) = size(t.children, i)
size(t::Table) = size(t.children)

# Adjust a table solution so that the aspect ratio matches tbl.aspect_ratio
#
# Returns:
#   true if the solution is feasibly, false if not
#
# Modifies:
#   x_solution, y_solution, w_solution, h_solution
#
function force_aspect_ratio!(tbl::Table,
                             x_solution::Vector, y_solution::Vector,
                             w_solution::Vector, h_solution::Vector)

    w0 = sum(w_solution[tbl.x_focus])
    h0 = sum(h_solution[tbl.y_focus])
    adj = (w0 / h0) / tbl.aspect_ratio

    # we can't expand either dimension (since it's presumably of maximum size)
    # so shrink one or the other.
    if adj > 1.0
        w = w0 / adj
        delta = w0 - w
        x_solution[1:tbl.x_focus.stop] .+= delta/2
        x_solution[tbl.x_focus.stop+1:end] .-= delta/2
        w_solution[tbl.x_focus] ./= adj
    else
        w = w0
        h = h0 * adj
        delta = h0 - h
        y_solution[1:tbl.y_focus.stop] .+= delta/2
        y_solution[tbl.y_focus.stop+1:end] .-= delta/2
        h_solution[tbl.y_focus] *= adj
    end
end

# Return true if the minwidth and minheight constraints on ctx are satisfied
# by the given width/height
function issatisfied(ctx::Context, width, height)
    eps = 1e-4
    return (minwidth(ctx) == nothing || width + eps >= minwidth(ctx)) &&
           (minheight(ctx) == nothing || height + eps >= minheight(ctx))
end

# Solve the table layout using a brute force approach, when a MILP isn't
# available.
function realize_brute_force(tbl::Table, drawctx::ParentDrawContext)
    m, n = size(tbl.children)

    maxobjective = 0.0
    optimal_choice = nothing
    feasible = false # is the current optimal_choice feasible

    # if the current solution is infeasible, we try to minimize badness,
    # which is basically "size needed" - "size available".
    minbadness = Inf

    focused_col_widths = Array{Float64}(undef, length(tbl.x_focus))
    focused_row_heights = Array{Float64}(undef, length(tbl.y_focus))

    # minimum sizes for each column and row
    minrowheights = Array{Float64}(undef, m)
    mincolwidths = Array{Float64}(undef, n)

    # convert tbl.fixed_configs to linear indexing
    fixed_configs = Any[
        Set([(j-1)*m + i for (i, j) in fixed_config])
        for fixed_config in tbl.fixed_configs
    ]

    # build equilavence classes of configurations basen on fixed_configs
    constrained_cells = Set()
    num_choices = [length(child) for child in tbl.children]
    num_group_choices = Any[]
    for fixed_config in fixed_configs
        push!(num_group_choices, num_choices[first(fixed_config)])
        for idx in fixed_config
            push!(constrained_cells, idx)
        end
    end

    for (idx, child) in enumerate(tbl.children)
        if length(child) > 1
            if !in(idx, constrained_cells)
                push!(num_group_choices, length(child))
                push!(fixed_configs, Set([idx]))
            end
        end
    end

    # compute the optimal column widths/row heights for fixed choice and
    # pre-computed minrowheight/mincolwidths.
    function update_focused_col_widths!(focused_col_widths)
        minwidth = sum(mincolwidths)
        total_focus_width =
            drawctx.box.a[1].value - minwidth + sum(mincolwidths[tbl.x_focus])

        if tbl.x_prop != nothing
            for k in 1:length(tbl.x_focus)
                if isnan(tbl.x_prop[k])
                    total_focus_width -= mincolwidths[tbl.x_focus[k]]
                end
            end

            for k in 1:length(tbl.x_focus)
                if !isnan(tbl.x_prop[k])
                    focused_col_widths[k] = tbl.x_prop[k] * total_focus_width
                else
                    focused_col_widths[k] = mincolwidths[tbl.x_focus[k]]
                end
            end
        else
            extra_width = total_focus_width - sum(mincolwidths[tbl.x_focus])
            for k in 1:length(tbl.x_focus)
                focused_col_widths[k] =
                        mincolwidths[tbl.x_focus[k]] + extra_width / length(tbl.x_focus)
            end
        end
    end

    function update_focused_row_heights!(focused_row_heights)
        total_focus_height =
                drawctx.box.a[2].value - sum(minrowheights) + sum(minrowheights[tbl.y_focus])

        if tbl.y_prop != nothing
            for k in 1:length(tbl.y_focus)
                if isnan(tbl.y_prop[k])
                    total_focus_height -= minrowheights(tbl.y_focus[k])
                end
            end

            for k in 1:length(tbl.y_focus)
                if !isnan(tbl.y_prop[k])
                    focused_row_heights[k] = tbl.y_prop[k] * total_focus_height
                else
                    focused_row_heights[k] = minrowheights[tbl.y_focus[k]]
                end
            end
        else
            extra_height = total_focus_height - sum(minrowheights[tbl.y_focus])
            for k in 1:length(tbl.y_focus)
                focused_row_heights[k] =
                        minrowheights[tbl.y_focus[k]] + extra_height / length(tbl.y_focus)
            end
        end
    end

    # for a given configuration, compute the minimum width for every column and
    # minimum height for every row. Return the penalty for choice.
    function update_mincolrow_sizes!(choice, minrowheights, mincolwidths)
        fill!(minrowheights, -Inf)
        fill!(mincolwidths, -Inf)
        penalty = 0.0
        for i in 1:m, j in 1:n
            isempty(tbl.children[i, j]) && continue

            choice_ij = choice[(j-1)*m + i]
            child = tbl.children[i, j][(choice_ij == 0 ? 1 : choice_ij)]
            penalty += child.penalty
            mw, mh = minwidth(child), minheight(child)
            if mw != nothing && mw > mincolwidths[j]
                mincolwidths[j] = mw
            end
            if mh != nothing && mh > minrowheights[i]
                minrowheights[i] = mh
            end
        end

        minrowheights[isfinite.(minrowheights) .== false] .= 0.0
        mincolwidths[isfinite.(mincolwidths) .== false] .= 0.0

        return penalty
    end

    it_count = 0
    group_choices = [l == 0 ? (0:0) : (1:l) for l in num_group_choices]
    choice = zeros(Int, m * n)
    optimal_choice = nothing
    if !isempty(group_choices)
        for group_choice in Iterators.product(group_choices...)
            it_count += 1
            for (l, k) in enumerate(group_choice)
                for p in fixed_configs[l]
                    choice[p] = k
                end
            end

            penalty = update_mincolrow_sizes!(choice, minrowheights, mincolwidths)

            minheightval = sum(minrowheights)
            minwidthval = sum(mincolwidths)

            update_focused_col_widths!(focused_col_widths)
            update_focused_row_heights!(focused_row_heights)

            objective = sum(focused_col_widths) + sum(focused_row_heights) - penalty

            # feasible?
            if minwidthval < drawctx.box.a[1].value && minheightval < drawctx.box.a[2].value &&
               all(focused_col_widths .>= mincolwidths[tbl.x_focus]) &&
               all(focused_row_heights .>= minrowheights[tbl.y_focus])
                if objective > maxobjective || !feasible
                    maxobjective = objective
                    minbadness = 0.0
                    optimal_choice = copy(choice)
                end
                feasible = true
            else
                badness = max(minwidthval - drawctx.box.a[1].value, 0.0) +
                          max(minheightval - drawctx.box.a[2].value, 0.0)
                if badness < minbadness && !feasible
                    minbadness = badness
                    optimal_choice = copy(choice)
                end
            end

            if optimal_choice === nothing
                optimal_choice = copy(choice)
            end
        end
    end

    if optimal_choice === nothing
        optimal_choice = zeros(Int, m * n)
    end

    update_mincolrow_sizes!(optimal_choice, minrowheights, mincolwidths)
    update_focused_col_widths!(focused_col_widths)
    update_focused_row_heights!(focused_row_heights)

    w_solution = mincolwidths
    h_solution = minrowheights

    for k in 1:length(tbl.x_focus)
        w_solution[tbl.x_focus[k]] = focused_col_widths[k]
    end

    for k in 1:length(tbl.y_focus)
        h_solution[tbl.y_focus[k]] = focused_row_heights[k]
    end

    x_solution = cumsum(mincolwidths) .- w_solution
    y_solution = cumsum(minrowheights) .- h_solution

    if tbl.aspect_ratio != nothing
        force_aspect_ratio!(tbl, x_solution, y_solution, w_solution, h_solution)
    end

    root = context(units=tbl.units, order=tbl.order)

    feasible = true
    for i in 1:m, j in 1:n
        if isempty(tbl.children[i, j])
            continue
        elseif length(tbl.children[i, j]) == 1
            ctx = copy(tbl.children[i, j][1])
        elseif length(tbl.children[i, j]) > 1
            idx = optimal_choice[(j-1)*m + i]
            ctx = copy(tbl.children[i, j][idx])
        end
        feasible == feasible && issatisfied(ctx, w_solution[j], h_solution[i])
        ctx.box = BoundingBox(
            x_solution[j]*mm, y_solution[i]*mm,
            w_solution[j]*mm, h_solution[i]*mm)
        compose!(root, ctx)
    end

    feasible || warn("Graphic may not be drawn correctly at the given size.")

    return root
end

# TODO: Enable this when we have a mechanism for making it optional
#if isinstalled("JuMP") &&
    #(isinstalled("GLPKMathProgInterface") ||
     #isinstalled("Cbc"))
    #include("table-jump.jl")
#else
realize(tbl::Table, drawctx::ParentDrawContext) = realize_brute_force(tbl, drawctx)
#end

function show(io::IO, t::Table)
    if get(io, :compact, false)
        println(io,"$(size(t.children,1))x$(size(t.children,2)) Table:")
        for i = 1:size(t.children,1)
            print(io, "  ")
            first = true
            for j = 1:size(t.children,2)
                first || print(io, ",")
                first = false
                show(io, t.children[i,j])
            end
            println(io)
        end
    else
        invoke(show, Tuple{IO, Any}, io, t)
    end
end


text:
# Elementary row operations

export row_swap!, row_scale!, row_add_mult!

"""
`row_swap!(A,i,j)` swaps rows `i` and `j` in the matrix `A`.
"""
function row_swap!(A::AbstractMatrix{T}, i::Int, j::Int) where {T}
    r, c = size(A)
    @assert (1 <= i <= r) && (1 <= j <= r) "Row index out of bounds"
    if i == j
        return
    end

    save_row = [A[i, k] for k = 1:c]

    for k = 1:c
        @inbounds A[i, k] = A[j, k]
    end

    for k = 1:c
        @inbounds A[j, k] = save_row[k]
    end
end

"""
`row_scale!(A,i,s)` multiplies all entries in row `i`
of `A` by `s`.
"""
function row_scale!(A::AbstractMatrix{T}, i::Int, s) where {T}
    r, c = size(A)
    @assert 1 <= i <= r "Row index out of bounds"

    for k = 1:c
        @inbounds A[i, k] = s * A[i, k]
    end
end

"""
`row_add_mult!(A,i,s,j)` adds `s` times row `i` to row `j`
in the matrix `A`.
"""
function row_add_mult!(A::AbstractMatrix{T}, i::Int, s, j::Int) where {T}
    r, c = size(A)
    @assert (1 <= i <= r) && (1 <= j <= r) "Row index out of bounds"

    save_row = [s * A[i, k] for k = 1:c]
    for k = 1:c
        @inbounds A[j, k] = A[j, k] + save_row[k]
    end
end

ext: jl


text:
#
# # DATE: 2011-09-15 CONTRIBUTOR: Unknown CITATION: Justo, Bazant, Kaxiras, Bulatov and Yip, Phys Rev B, 58, 2539 (1998)
#
# # EDIP parameters for various elements and mixtures
# # multiple entries can be added to this file, LAMMPS reads the ones it needs
# # these entries are in LAMMPS "metal" units
#
# # format of a single entry (one or more lines)
# #
# #   element 1, element 2, element 3,
# #     A   B   cutoffA   cutoffC   alpha   beta   eta
# #     gamma   lambda    mu   rho   sigma   Q0
# #     u1   u2   u3   u4
# #
# # units for each parameters:
# #     A , lambda are in eV
# #     B, cutoffA, cutoffC, gamma, sigma are in Angstrom
# #     alpha, beta, eta, mu, rho, Q0, u1-u4 are pure numbers
#
# # Here are the original parameters in metal units, for Silicon from:
# # J. F. Justo, M. Z. Bazant, E. Kaxiras, V. V. Bulatov, S. Yip
# #       Phys. Rev. B 58, 2539 (1998)
# #
#
# Si Si Si 7.9821730 1.5075463 3.1213820 2.5609104 3.1083847 0.0070975 0.2523244
#          1.1247945 1.4533108 0.6966326 1.2085196 0.5774108 312.1341346
#          -0.165799 32.557 0.286198 0.66

module EDIP

using JuLIP
JuLIP.Potentials.WrappedAnalyticFunction
import JuLIP
import JuLIP.Potentials: SitePotential, WrappedAnalyticFunction

struct EdIP <: SitePotential
   neigf::WrappedAnalyticFunction
end

function EdIP(;
                  A        = 7.9821730
                  B        = 1.5075463
                  cutoffA  = 3.1213820   # a = cutoffA
                  cutoffC  = 2.5609104   # c = cutoffC
                  alpha    = 3.1083847
                  beta     = 0.0070975
                  eta      = 0.2523244
                  gamma    = 1.1247945
                  lambda   = 1.4533108
                  mu       = 0.6966326
                  rho      = 1.2085196
                  sigma    = 0.5774108
                  Q0       = 312.1341346
                  u1       = -0.165799
                  u2       = 32.557
                  u3       = 0.286198
                  u4       = 0.66 )
   neigf = let a = cutoffA, c = cutoffC, α = alpha
      @analytic x -> exp(α / (1 - 1 / x^3))
   end
   moll2 = let σ = sigma, a = cutoffA, A = A
      @analytic r -> A *  exp( σ / (r - a) )
   end
   moll3 = let γ = gamma, a = cutoffA
      @analytic r -> exp( γ / (r - a) )
   end
   pf = let β = beta
      @analytic Z -> exp( - β * Z^2 )
   end
   v2in = let B = B, ρ = rho
      @analytic r -> (B / r)^ρ
   end
end

function neigf(V::EdIP, r)
   x = (r - V.cutoffC) / (V.cutoffA - V.cutoffC)
   if x <= 1e-7
      return 1.0
   elseif x >= 1 - 2.5e-2
      return 0.0
   else
      return V.neigf(x)
   end
end

function neigf_d(V::EdIP, r)
   x = (r - V.cutoffC) / (V.cutoffA - V.cutoffC)
   dxdr = 1 / (V.cutoffA - V.cutoffC)
   if (x <= 1e-7) || (x >= 1 - 2.5e-2)
      return 0.0
   else
      return (@D V.neigf(x)) * dxdr
   end
end

V2(V::EdIP, R, Z) = (V.v2in(norm(R)) - V.pf(Z)) * V.moll2(norm(R))

function  V2_d(V::EdIP, R, Z)   # -> ∂V2/∂R, ∂V2/∂Z
   r = norm(R)
   v2in = V.v2in(r)
   p = V.pf(Z)
   moll2 = V.moll2(r)
   DV2Dr = (@D V.v2in(r)) * moll2 + (v2in - p) * (@D V.moll2(r))
   return (DV2Dr / r) *  R, (- moll2) * (@D V.pf(Z))
end

function V3(V::EdIP, R1, R2, Z)
   r1 = norm(R1)
   r2 = norm(R2)
   g1 = V.moll3(r1)
   g2 = V.moll3(r2)
   l12 = dot(R1, R2) / (r1*r2)
   h = hf(V, l12, Z)
   return g1 * g2 * h
end

function evaluate!(tmp, V::EdIP, R::AbstractVector{<: JVec})

   # compute the neighbour counter + the r's along the way
   Zi = 0.0
   for j = 1:length(R)
      Zi += neigf(V, norm(R[j]))
   end

   # compute the V2's
   Es = 0.0
   for j = 1:length(R)
      Es += V2(V, R[j], Zi)
   end

   # compute the V3's
   for j1 = 2:length(R), j2 = 1:j1-1
      l12 = dot(R[j1], R[j2]) / (tmp.r[j1]*tmp.r[j2])
      Es += V3(V, R[j1], R[j2], Zi)
   end
end


text:
module CURAND

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, libraryPropertyType, DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK
using ..CUDA: libcurand, @retry_reclaim

using CEnum

# core library
include("libcurand_common.jl")
include("error.jl")
include("libcurand.jl")

# low-level wrappers
include("wrappers.jl")

# high-level integrations
include("random.jl")

# cache for created, but unused handles
const idle_curand_rngs = HandleCache{CuContext,RNG}()
const idle_gpuarray_rngs = HandleCache{CuContext,GPUArrays.RNG}()

function default_rng()
    cuda = CUDA.active_state()

    # every task maintains library state per device
    LibraryState = @NamedTuple{rng::RNG}
    states = get!(task_local_storage(), :CURAND) do
        Dict{CuContext,LibraryState}()
    end::Dict{CuContext,LibraryState}

    # get library state
    @noinline function new_state(cuda)
        new_rng = pop!(idle_curand_rngs, cuda.context) do
            RNG()
        end

        finalizer(current_task()) do task
            push!(idle_curand_rngs, cuda.context, new_rng) do
                # no need to do anything, as the RNG is collected by its finalizer
            end
        end

        Random.seed!(new_rng)
        (; rng=new_rng)
    end
    state = get!(states, cuda.context) do
        new_state(cuda)
    end

    return state.rng
end

function GPUArrays.default_rng(::Type{<:CuArray})
    cuda = CUDA.active_state()

    # every task maintains library state per device
    LibraryState = @NamedTuple{rng::GPUArrays.RNG}
    states = get!(task_local_storage(), :GPUArraysRNG) do
        Dict{CuContext,LibraryState}()
    end::Dict{CuContext,LibraryState}

    # get library state
    @noinline function new_state(cuda)
        new_rng = pop!(idle_gpuarray_rngs, cuda.context) do
            N = attribute(cuda.device, DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK)
            buf = CuArray{NTuple{4, UInt32}}(undef, N)
            GPUArrays.RNG(buf)
        end

        finalizer(current_task()) do task
            push!(idle_gpuarray_rngs, cuda.context, new_rng) do
                # no need to do anything, as the RNG is collected by its finalizer
            end
        end

        Random.seed!(new_rng)
        (; rng=new_rng)
    end
    state = get!(states, cuda.context) do
        new_state(cuda)
    end

    return state.rng
end

@deprecate seed!() CUDA.seed!()
@deprecate seed!(seed) CUDA.seed!(seed)

end


text:
using ArbitrageFreeSABR, Test

@testset "HaganSet" begin
    forward = 1.0; expiry = 1.0;
    alpha = 0.35; beta = 0.25; nu = 1.0; rho = -0.1;
    N = 500; timesteps = 5; nd= 4;
    maturity = SABRMaturity(alpha,beta,rho,nu,forward,expiry,ArbitrageFreeSABRModel())
    density = makeTransformedDensityLawsonSwayne(maturity, N, timesteps, nd)
    @test isapprox(priceTransformedDensity(density, true, forward, ArbitrageFreeSABR.midpoint), 0.149701955629, atol=1e-10)
    @test isapprox(priceTransformedDensity(density, true, forward, ArbitrageFreeSABR.linear), 0.14970195574320477, atol=1e-10)
end

@testset "ImpliedVolatilitySet" begin
    price = 271.43234885190117
    strike = 275.0
    f = 356.73063159822254
    tte = 1.5917808219178082
    isCall = false
    tolerance = 1e-8
    vol = impliedVolatilityLiSORTS(price, isCall, strike, f, tte, 1.0, 0.0, tolerance,64)
    @test isapprox(blackScholesFormula(isCall, strike, f, vol*vol*tte, 1.0, 1.0), price, atol=tolerance)
end

@testset "BachelierVolatilitySet" begin
    price = 271.43234885190117
    strike = 275.0
    f = 356.73063159822254
    tte = 1.5917808219178082
    isCall = false
    tolerance = 1e-13
    vol = bachelierVolatility(price, isCall, strike, f, tte,1.0)
    call = f - strike + price
    tv = call - max(f-strike,0.0)
    y = tv/abs(f-strike)*4*sqrt(pi)
    @test isapprox(bachelierFormula(isCall, strike, f, vol*vol*tte, 1.0), price, atol=tolerance)
end


text:
using .StaticArrays

# import TensorCast: static_slice, static_glue # because Revise doesn't track these

"""
    static_slice(A, sizes)

Slice array `A` into an `AbstractArray{StaticArray}`,
with the slices always belonging to the first few indices, i.e. `code = (:,:,...,*,*)`.
Expects `sizes = Size(size(A)[1:ncolons])` for stability.

Typically produces `reshape(reinterpret(SArray...`; copy this to get `Array{SArray}`,
or call `static_slice(A, sizes, false)` to omit the reshape.
"""
@inline function static_slice(A::AbstractArray{T,N}, sizes::Size{tup}, finalshape::Bool=true) where {T,N,tup}
    IN = length(tup)
    IT = SArray{Tuple{tup...}, T, IN, prod(tup)}
    if N-IN>1 && finalshape
        finalsize = size(A)[1+IN:end]
        reshape(reinterpret(IT, vec(A)), finalsize)
    else
        reinterpret(IT, vec(A)) # always a vector
    end
end

"""
    static_glue(A)

Glues the output of `static_slice` back into one array, again with `code = (:,:,...,*,*)`.
"""
@inline function static_glue(A::AbstractArray{IT,ON}, finalshape::Bool=true) where {IT,ON}
    IT <: StaticArray || error("static_glue needs an array of StaticArrays")
    if finalshape
        finalsize = (size(IT)..., size(A)...)
        reshape(reinterpret(eltype(IT), A), finalsize)
    else
        reinterpret(eltype(eltype(A)), A)
    end
end

function auto_glue(A::AbstractArray{IT,ON}, code::Tuple) where {IT<:StaticArray,ON}
    if iscodesorted(code)
        static_glue(A)
    else
        copy_glue(A, code)
    end
end


text:
using StatsBase

input = joinpath(@__DIR__, "input")
lines = readlines(input)

template = lines[1]
rules = lines[3:end]

rules = map(line -> split(line, " -> "), rules)

function replacements(rules)
    rs = Dict{String, Tuple{String, String}}()

    for (pair, insert) in rules
        l = first(pair) * insert
        r = insert * last(pair)
        rs[pair] = (l, r)
    end

    return rs
end

function expand(template, rules, n)
    counts = Dict{String, Int}()
    for i in 1:(length(template) - 1)
        pair = template[i:i+1]
        haskey(counts, pair) ? counts[pair] += 1 : counts[pair] = 1
    end

    for _ in 1:n
        ncounts = Dict{String, Int}()
        for (pair, ns) in collect(counts)
            l, r = rules[pair]
            haskey(ncounts, l) ? ncounts[l] += ns : ncounts[l] = ns
            haskey(ncounts, r) ? ncounts[r] += ns : ncounts[r] = ns
        end
        counts = ncounts
    end

    return counts
end

function solve(counts)
    cm = Dict{Char, Int}()
    for pair in keys(counts)
        ch = first(pair)
        haskey(cm, ch) ? cm[ch] += counts[pair] : cm[ch] = counts[pair]
    end

    cm[last(template)] += 1

    mini, maxi = extrema(values(cm))
    return maxi - mini
end

rules = replacements(rules)
counts = expand(template, rules, 10)
p1 = solve(counts)
counts = expand(template, rules, 40)
p2 = solve(counts)

println("-----------------------------------------------------------------------")
println("extended polymerization -- part one :: $p1")
println("extended polymerization -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 2027)
@assert(p2 == 2265039461737)


text:
#=
BOCF solver.

Notice that the unicode symbol "⁻" (a minus exponent, e.g. τ⁻), represents
"minus", e.g. `tau_minus` is in fact `τ⁻`.
=#

include("fast_convolutions.jl")

function makesim(sts::Union{STS{:BuenoOrovioCherryFenton}, STS{:bocf}})
    @unpack N, T, Δt, S, ic = sts
    Ttr = S
    Nx, Ny = typeof(N) <: Tuple ? (N[1], N[2]) : (N, N)

    # TODO: Change the model to use Float32 numbers instead

    u, v, w, s = [zeros(Nx, Ny) for i in 1:4]

    pset = get(sts.p, :pset, "tnpp")

    # here we get the initial condition
    if ic isa Union{Integer, Nothing} # random seed
        initialize_random!(u, ic isa Integer ? ic : 42)
    elseif ic == "spiral_virtheart"
        initialize_spiral_virtheart!(u, v)
    elseif ic == "spiral"
        initialize_spiral!(u, v, w)
    else error("unknown i.c.")
    end

    t, U, V, W, S = bocf(u, v, s, w, Ttr, T, Δt; sts.p...)
    @strdict t sts U V W S
end

################################################
# PDE Stepping
################################################
"""
    bocf(u, v, s, w, Ttr, T, Δt; kwargs...) → t, U, V, W, S
Simulate the Bueno-Orovio, Cherry & Fenton model for cardiac dynamics
with initial conditions `u, v, s, w`. Simulate for total time `T`
with sampling time `Δt`, while transiently integrating for `Ttr`.

A note on the initial condition, `sts.ic`.
The initial condition is by default `ic = 42`.
If number, random initial condition with `ic` seed is generated.
Otherwise `ic` can be `"spiral_virtheart"` or `"spiral"`, which initializes
a specific pattern.

The keywords of the BOCF model are (given to `sts.p` as a `NamedTuple`)
* `D = 0.2` the diffusion constant
* `pset = "tnpp"` the parameter set. Can be `"tnpp", "thomas", "epi", "virtheart", "pb"`.
  See the source code for the actual parameters (too many to list).
* `dt = Δt` : integration timestep (NOT the sampling time `Δt`).
"""
function bocf(u, v, s, w, Ttr, T, Δt;
        pset = "tnpp", D = 0.2, dt = Δt)
        # TODO: what is a good default dt?

    @assert Δt ≥ dt
    simpars = bocf_parameters(pset)
    pushfirst!(simpars, D) # add diffusion constant

    Nx, Ny = size(u)
    olds = [copy(u), copy(v), copy(w), copy(s)]
    derivs = [copy(u), copy(v), copy(w), copy(s)]
    aux = [zeros(Nx, Ny) for i in 1:7]
    currents = [zeros(Nx, Ny) for i in 1:3]

    TU = typeof(u)
    U, V, W, S = TU[], TU[], TU[], TU[]

    for j in 1:Ttr÷dt
        bocf_step!(u, v, w, s, olds, aux, currents, derivs, simpars, dt)
    end
    push!(U, copy(w)); push!(V, copy(v)); push!(W, copy(w)); push!(S, copy(s))
    for j in 1:Int(T÷Δt) # how many data to save
        for i in 1:Δt÷dt # how often to sample
            bocf_step!(u, v, w, s, olds, aux, currents, derivs, simpars, dt)
        end
        push!(U, copy(w)); push!(V, copy(v)); push!(W, copy(w)); push!(S, copy(s))
    end
    return 0:Δt:T, U, V, W, S

end

function bocf_step!(u, v, w, s, olds, aux, currents, derivs, simpars, dt)

    D,u_o,u_u,θ_v,θ_w,θ_v⁻,θ_o,τ_v1⁻,τ_v2⁻,
    τ_v_plus,τ_w1⁻,τ_w2⁻,k_w⁻,u_w⁻,τ_w_plus,
    τ_fi,τ_o1,τ_o2,τ_so1,τ_so2,k_so,u_so,τ_s1,τ_s2,k_s,
    u_s,τ_si,τ_winfinity,w_inf_star = simpars

    # update old fields (only necessary for boundary conditions)
    # TODO: this can be optimized to only update the values
    # used in the _set_boundaries! function
    for (new, old) in zip((u, v, w, s), olds)
        old .= new
    end

    _update_constants!(aux, u, v, w, s, simpars)
    τ_v⁻, τ_w⁻, τ_so, τ_s, τ_o, v_infinity, w_infinity = aux
    J_fi, J_so, J_si = currents

    # Compute currents
    _current_fi!(J_fi, u, v, θ_v, u_u, τ_fi)
    _current_so!(J_so, u, u_o, θ_w, τ_o, τ_so)
    _current_si!(J_si, u, θ_w, w, s, τ_si)

    dudt, dvdt, dwdt, dsdt = derivs
    lap = laplace(u)
    @. dudt = D * lap - (J_fi + J_so + J_si)
    @. dvdt = (1.0 - H(u - θ_v)) * (v_infinity - v) / τ_v⁻ - H(u - θ_v) * v / τ_v_plus
    @. dwdt = (1.0 - H(u - θ_w)) * (w_infinity - w) / τ_w⁻ - H(u - θ_w) * w / τ_w_plus
    @. dsdt = ((1.0 + tanh(k_s * (u - u_s))) / 2.0 - s) / τ_s

    @. u += dt * dudt
    @. v += dt * dvdt
    @. w += dt * dwdt
    @. s += dt * dsdt

    _set_boundaries!((u, v, w, s), olds)
    return u, v, w, s
end

function _set_boundaries!(fields, old_fields)
    # TODO: We have to check if this is correct, because I took this from
    # Python code, which interchanges columns with rows
    # TODO: What is happening in this function seems really weird for me.
    # Is it really correct? The end values before the step go the the near
    # end values after the step? That's so weird.
    for (field, old_field) in zip(fields, old_fields)
        field[:, 1] .= @view old_field[:, 2]
        field[:, end] .= @view old_field[:, end-1]
        field[1, :] .= @view old_field[2, :]
        field[end, :] .= @view old_field[end-1, :]
    end
end

# TODO: `laplace` creates arrays. It could be optmized by replacing the
# convolution with an explicit loop, because
# laplace = -4*copy(field)
# laplace += roll(field, +1, axis=0)
# laplace += roll(field, -1, axis=0)
# laplace += roll(field, +1, axis=1)
# laplace += roll(field, -1, axis=1)
# notice however that our fastconv is using an optimized FFT algorithm.
# I doubt that normal loops will be faster.
# Laplacian matrix
const Λ = [1  4  1;
           4 -20 4;
           1  4  1]./6
laplace(u) = fastconv(u, Λ)[2:end-1, 2:end-1]

H(y) = Float64(y > 0)

# TODO: All currents can be hugely optimized by normal loop over
# the operations, and if statements on the result of H(), so that
# all other operations are not done

function _current_fi!(cfi, u, v, θ_v, u_u, τ_fi)
    @. cfi = -v * H(u - θ_v) * (u - θ_v) * (u_u - u) / τ_fi
end

function _current_so!(cso, u, u_o, θ_w, τ_o, τ_so)
    @. cso = (u - u_o) * (1 - H(u - θ_w)) / τ_o + H(u - θ_w) / τ_so
end

function _current_si!(csi, u, θ_w, w, s, τ_si)
    @. csi = -H(u - θ_w) * w * s / τ_si
end

function _update_constants!(aux, u, v, w, s, simpars)
    τ_v⁻, τ_w⁻, τ_so, τ_s, τ_o, v_infinity, w_infinity = aux

    D,u_o,u_u,θ_v,θ_w,θ_v⁻,θ_o,τ_v1⁻,τ_v2⁻,
    τ_v_plus,τ_w1⁻,τ_w2⁻,k_w⁻,u_w⁻,τ_w_plus,
    τ_fi,τ_o1,τ_o2,τ_so1,τ_so2,k_so,u_so,τ_s1,τ_s2,k_s,
    u_s,τ_si,τ_winfinity,w_inf_star = simpars

    @. τ_v⁻ = (1.0 - H(u - θ_v⁻)) * τ_v1⁻ + H(u - θ_v⁻) * τ_v2⁻
    @. τ_w⁻ = τ_w1⁻ + (τ_w2⁻ - τ_w1⁻) * (1 + tanh(k_w⁻ * (u - u_w⁻))) / 2.0
    @. τ_so      = τ_so1      + (τ_so2      - τ_so1)      * (1 + tanh(k_so      * (u - u_so)))      / 2.0
    @. τ_s = (1.0 - H(u - θ_w)) * τ_s1 + H(u - θ_w) * τ_s2
    @. τ_o = (1.0 - H(u - θ_o)) * τ_o1 + H(u - θ_o) * τ_o2

    @. v_infinity = Float64(u .< θ_v⁻)
    @. w_infinity = (1.0 - H(u - θ_o)) * (1.0 - u / τ_winfinity) + H(u - θ_o) * w_inf_star
    return
end

################################################
# Initial conditions
################################################
function initialize_random!(u, ic)
    Nx, Ny = size(u)
    rng = MersenneTwister(ic)
    n = 1 # ceil(Int, 1/deltaX), what is deltaX ?
    u .= rand(rng, Nx, Ny)
end

function initialize_spiral_virtheart!(u, v)
    Nx, Ny = size(u)
    for i in 1:Ny
        for j in (Nx÷2):Nx
            t  = (i - Ny÷2)*0.1
            t2 = (i - Ny÷2 + 20)*0.05
            u[j, i] = 1.5*exp(-t*t)
            v[j, i] = 1.0 - 0.9*exp(-t2*t2)
        end
    end
end

function initialize_spiral!(u, v, w)
    Nx, Ny = size(u)
    for x in 1:(Nx÷4)
        for y in 1:Ny
            u[x, y] = 0.5
        end
    end
    for x in (3Nx÷4):Nx
        for y in 1:Ny
            u[x, y] = 0.5
        end
    end
    for x in 1:(Nx÷2)
        for y in 1:(2Ny÷5)
            v[x, y] = 0.5
        end
    end
    for x in (Nx÷2):Nx
        for y in (3Ny÷5):Ny
            v[x, y] = 0.5
        end
    end
    fill!(w, 1.0)
end

################################################
# Parameter sets
################################################
function bocf_parameters(ic)
    if ic == "tnpp"
        u_o = 0.0
        u_u = 1.58
        θ_v = 0.3
        θ_w = 0.015
        θ_v⁻ = 0.015
        θ_o = 0.006
        τ_v1⁻ = 60
        τ_v2⁻ = 1150
        τ_v_plus = 1.4506
        τ_w1⁻ = 70
        τ_w2⁻ = 20
        k_w⁻ = 65
        u_w⁻ = 0.03
        τ_w_plus = 280
        τ_fi = 0.11
        τ_o1 = 6
        τ_o2 = 6
        τ_so1 = 43
        τ_so2 = 0.2
        k_so = 2
        u_so = 0.65
        τ_s1 = 2.7342
        τ_s2 = 3
        k_s = 2.0994
        u_s = 0.9087
        τ_si = 2.8723
        τ_winfinity = 0.07
        w_inf_star = 0.94
    elseif ic == "virtheart"
        u_o = 0
        u_u = 1.58
        θ_v = 0.3
        θ_w = 0.015
        θ_v⁻ = 0.015
        θ_o = 0.006
        τ_v1⁻ = 60
        τ_v2⁻ = 60
        τ_v_plus = 1.4506
        τ_w1⁻ = 170
        τ_w2⁻ = 120
        k_w⁻ = 65
        u_w⁻ = 0.03
        τ_w_plus = 280
        τ_fi = 0.2
        τ_o1 = 6
        τ_o2 = 6
        τ_so1 = 43
        τ_so2 = 0.2
        k_so = 2
        u_so = 0.65
        τ_s1 = 2.7342
        τ_s2 = 3
        k_s = 2.0994
        u_s = 0.9087
        τ_si = 3.8723
        τ_winfinity = 0.07
        w_inf_star = 0.94
    elseif ic == "thomas"
        u_o = 0#13.03
        τ_o1 = 33.25
        θ_w = 800
        τ_so2 = 0.85
        τ_v1⁻ = 0.45
        τ_s1 = 0.04
        τ_w1⁻ = 0.45
        u_s = 0.04
        u_w⁻ = 0.45
        w_inf_star = 0.04
        τ_fi = 12.5
        θ_v = 1250
        τ_so1 = 0.13
        θ_o = 0.45
        u_so = 0.04
        τ_v_plus = 0.45
        k_s = 0.04
        k_w⁻ = 0.45
        τ_winfinity = 0.04
        u_u = 19.6
        τ_o2 = 29
        θ_v⁻ = 40
        k_so = 0.04
        τ_v2⁻ = 0.45
        τ_s2 = 0.04
        τ_w2⁻ = 0.45
        τ_si = 0.04
        τ_w_plus = 0.45
    elseif ic == "epi"
        u_o = 0
        u_u = 1.55
        θ_v = 0.3
        θ_w = 0.13
        θ_v⁻ = 0.006
        θ_o = 0.006
        τ_v1⁻ = 60
        τ_v2⁻ = 1150
        τ_v_plus = 1.4506
        τ_w1⁻ = 60
        τ_w2⁻ = 15
        k_w⁻ = 65
        u_w⁻ = 0.03
        τ_w_plus = 200
        τ_fi = 0.11
        τ_o1 = 400
        τ_o2 = 6
        τ_so1 = 30.0181
        τ_so2 = 0.9957
        k_so = 2.0458
        u_so = 0.65
        τ_s1 = 2.7342
        τ_s2 = 16
        k_s = 2.0994
        u_s = 0.9087
        τ_si = 1.8875
        τ_winfinity = 0.07
        w_inf_star = 0.94
    elseif ic == "pb"
        u_o = 0
        u_u = 1.45
        θ_v = 0.35
        θ_w = 0.13
        θ_v⁻ = 0.175
        θ_o = 0.006
        τ_v1⁻ = 10
        τ_v2⁻ = 1150
        τ_v_plus = 1.4506
        τ_w1⁻ = 140
        τ_w2⁻ = 6.25
        k_w⁻ = 65
        u_w⁻ = 0.015
        τ_w_plus = 326
        τ_fi = 0.105
        τ_o1 = 400
        τ_o2 = 6
        τ_so1 = 30.0181
        τ_so2 = 0.9957
        k_so = 2.0458
        u_so = 0.65
        τ_s1 = 2.7342
        τ_s2 = 16
        k_s = 2.0994
        u_s = 0.9087
        τ_si = 1.8875
        τ_winfinity = 0.175
        w_inf_star = 0.9
    else
        error("Unknown parameter set")
    end

    simpars =
    [u_o,u_u,θ_v,θ_w,θ_v⁻,θ_o,τ_v1⁻,τ_v2⁻,
    τ_v_plus,τ_w1⁻,τ_w2⁻,k_w⁻,u_w⁻,τ_w_plus,
    τ_fi,τ_o1,τ_o2,τ_so1,τ_so2,k_so,u_so,τ_s1,τ_s2,k_s,
    u_s,τ_si,τ_winfinity,w_inf_star]
 end


text:
using TcpInstruments
using Plots
scope = initialize(AgilentDSOX4034A,  "10.1.30.32")
data = get_data(scope, 1)
plot(data)


text:
module RandomFourierFeatures

using KernelFunctions
using BayesianLinearRegressors
using AbstractGPs
using Distributions
using LinearAlgebra
using FillArrays
using Random

export sample_rff_basis, build_rff_weight_space_approx

include("random_fourier_features.jl")

end


text:
export Logit

"""
Description:

    Logit mapping between a real variable in1 ∈ R and a binary variable out ∈ {0, 1}.

    f(out, in1, xi) =  Ber(out | σ(in1))
                    >= exp(in1*out) σ(xi) exp[-(in1 + xi)/2 - λ(xi)(in1^2 - xi^2)], where
               σ(x) =  1/(1 + exp(-x))
               λ(x) =  (σ(x) - 1/2)/(2*x)

Interfaces:

    1. out (binary)
    2. in1 (real)
    3. xi (auxiliary variable)

Construction:

    Logit(out, in1, xi)
"""
mutable struct Logit <: SoftFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol, Interface}

    function Logit(out, in1, xi; id=generateId(Logit))
        @ensureVariables(out, in1, xi)
        self = new(id, Vector{Interface}(undef, 3), Dict{Symbol,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:in1] = self.interfaces[2] = associate!(Interface(self), in1)
        self.i[:xi] = self.interfaces[3] = associate!(Interface(self), xi)

        return self
    end
end

slug(::Type{Logit}) = "σ"

logisticSigmoid(x::Float64) = 1/(1 + exp(-x))
logLogisticSigmoid(x::Float64) = -log(1 + exp(-x))
logisticLambda(x::Float64) = (logisticSigmoid(x) - 0.5)/(2*x)

# Average energy functionals
function averageEnergy(::Type{Logit}, marg_out::ProbabilityDistribution{Univariate}, marg_in1::ProbabilityDistribution{Univariate}, marg_xi::ProbabilityDistribution{Univariate})
    xi_hat = unsafeMode(marg_xi)

    logisticLambda(xi_hat)*(unsafeMean(marg_in1)^2 + unsafeCov(marg_in1) - xi_hat^2) +
    0.5*(unsafeMean(marg_in1) + xi_hat) -
    logLogisticSigmoid(xi_hat) -
    unsafeMean(marg_in1)*unsafeMean(marg_out)
end


text:
module tracingtriangle # --------------------------------------------------------------

export longest_path

function longest_path(tree, row=1, col=1)
    row == length(tree) && return tree[row][col]

    left = longest_path(tree, row + 1, col)
    right = longest_path(tree, row + 1, col + 1)
    return max(left, right) + tree[row][col]
end

end # module tracingtriangle

module test # -------------------------------------------------------------------------

using Test
using ..tracingtriangle

@testset "Sample Test" begin
    input_one = [[1], [2, 3], [1, 5, 1]]
    @test longest_path(input_one) == 9

    input_two = [[6], [4, 4], [1, 2, 1], [5, 4, 3, 2]]
    @test longest_path(input_two) == 16

    input_three = [[5]]
    @test longest_path(input_three) == 5

    input_four = [[1], [1, 1], [1, 1, 1], [2, 1, 1, 1]]
    @test longest_path(input_four) == 5

    input_five = [[0], [0, 0], [0, 0, 0]]
    @test longest_path(input_five) == 0
end

end # module test

max_stars_repo_name: JenniferToops/CodingDojo

max_issues_repo_name: JenniferToops/CodingDojo

max_forks_repo_name: JenniferToops/CodingDojo


text:
using Base.Test
using MPI

@test !MPI.Initialized()
MPI.Init()
@test MPI.Initialized()

@test MPI.Comm(MPI.CComm(MPI.COMM_WORLD)).val == MPI.COMM_WORLD.val

@test !MPI.Finalized()
MPI.Finalize()
@test MPI.Finalized()


text:
using RBNF
using MLStyle
using PrettyPrint
using DataStructures

struct ReMLLang
end

second((a, b)) = b
second(vec::V) where V <: AbstractArray = vec[2]

join_token_as_str = xs -> join(x.str for x in xs)

maybe_to_bool(::Some{Nothing}) = false
maybe_to_bool(::Some{T}) where T = true

"""
nonempty sequence
"""
struct GoodSeq{T}
    head :: T
    tail :: Vector{T}
end

get_list(::Some{Nothing}) = []
get_list(a::Some) = let (a, b) = a.value; [a, b...] end
join_rule(sep, x) = begin
    :([$x, ([$sep, $x] % second){*}].? % get_list)
end

Token = RBNF.Token

get_loc(token::Token{T}) where T = (t=T, lineno=token.lineno, colno=token.colno)
get_str(token::Token{T}) where T = token.str
RBNF.typename(name:: Symbol, ::Type{ReMLLang}) = Symbol(:R, name)

RBNF.@parser ReMLLang begin
    # define the ignores
    ignore{space, comment}

    # define keywords
    reserved = [true, false]

    @grammar
    # necessary

    Bind      := [name=id %get_str, '=', value=Exp]
    Let       := [loc=:let %get_loc, rec=:rec.? % maybe_to_bool,
                        binds=join_rule(:and, Bind),
                     :in, body=Exp]

    Fun       := [loc=:fn %get_loc, "(", args=join_rule(",", (id % get_str)), ")", "->", body=Exp]
    Match     := [loc=:match %get_loc, sc=Exp, :with,
                        cases=(['|', Comp, "->", Exp] % ((_, case, _, body), ) -> (case, body)){*},
                    :end.?] # end for nested match
    If        := [loc=:if %get_loc, cond=Exp, :then,
                        br1=Exp,
                     :else,
                        br2=Exp]
    NestedExpr = ['(', value=Exp, ')'] => _.value
    Var       := value=id %get_str
    Block     := [loc='{' %get_loc, stmts=Stmt{*}, '}']
    Atom      =  Nil | NestedExpr | Num | Str | Boolean | Var | List
    Attr      := [value=Atom, attrs=(['.', id % get_str] % second){*}]
    Call      := [fn=Attr, ['(', args = join_rule(",", Exp), ')'].?]
    List      := [loc='[', elts=join_rule(',', Exp), ']']
    Comp      = Call | Let | Fun | Match | If | Block
    Op        := ['`', name=_, '`']
    Top       := [hd=Comp, tl=[Op, Comp]{*}]
    Custom    := [value=Exp, [',', next=Custom].?] | [kw=id, value=Exp, [',', next=Custom].?]
    Exp       := [top=Top, [do_custom::Bool='{' => true, [custom=Custom].?, '}'].?]

    id_str    = id%get_str
    Define    := [loc=:def %get_loc, name=id_str, '=', value=Exp]
    Infix     := [loc=:infix %get_loc, name=id_str, prec=integer %get_str, is_right="right".? % maybe_to_bool]
    Foreign   := [loc=:foreign %get_loc, paths=join_rule('.', id_str)]
    Import    := [loc=:import %get_loc, is_qual=:qualified.? %maybe_to_bool, paths=join_rule('.', id_str)]
    Ops       := [loc=:export %get_loc, names=id_str{*}]

    Stmt      =  [Exp | Import | Module, ';'.?] % first
    TopStmt   = Define | Import | Infix | Stmt | Foreign | Ops
    Module    := [loc=:module %get_loc, name=id_str, params=id_str{*}, :where, stmts=TopStmt{*}, :end.?]

    Str       :=  value=str
    Nil       := ["()"]
    Num       := [neg="-".? % maybe_to_bool, (int=integer) | (float=float)]
    Boolean   := value=("true" | "false")

    @token
    comment   := @quote ("(*", "\\*)", "*)")
    str       := @quote ("\"", "\\\"", "\"")
    id        := r"\G[A-Za-z_]{1}[A-Za-z0-9_\!]*"
    float     := r"\G([0-9]+\.[0-9]*|[0-9]*\.[0.9]+)([eE][-+]?[0-9]+)?"
    integer   := r"\G([1-9]+[0-9]*|0)"
    space     := r"\G\s+"
end

function runparser(a, v :: Symbol)
    runparser(a, Val(v))
end

function runparser(a, v :: String)
    runparser(a, Symbol(v))
end

struct ParserFailed <: Exception
    msg :: String
end

function runparser(parser :: F, source_code :: String, ::Val{:rexp}) where F <: Function
    tokens = RBNF.runlexer(ReMLLang, source_code)
    ast, ctx = RBNF.runparser(parser, tokens)
    if ctx.maxfetched >= ctx.tokens.length
        return ast
    end
    token = tokens[ctx.maxfetched]
    lineno = token.lineno
    colno = token.colno
    str   = token.str
    throw(ParserFailed("parsing error at $(repr(str)) at lineno $lineno, colno $colno"))
end

function runparser(source_code :: String, repr_form::Val, mode::Val{:module})
    runparser(Module, source_code, repr_form)
end

function runparser(source_code :: String, repr_form::Val, mode::Val{:stmt})
    runparser(TopStmt, source_code, repr_form)
end

function runparser(source_code :: String, repr_form::Symbol, mode::Symbol)
    runparser(Module, source_code, Val(repr_form), Val(mode))
end

max_forks_repo_name: thautwarm/Rupy


text:
module ExpressionUtils

using Compat

export is_funcdef_expr, get_funcdef_expr, funcdef_longform
export funcdef_name, funcdef_params, funcdef_args, funcdef_argnames, funcdef_argtypeexprs
export walk, expr_replace, expr_bind, expr_bindings

include("blocks.jl")
include("functions.jl")
include("utils.jl")

end


text:
using Logging
using LarSurf
using SparseArrays
using Pio3d
using BenchmarkTools

larmodel = LarSurf.Lar.cuboidGrid([1,1,2], true)
V, topology = larmodel
VV, EVorig, FV, CV = topology
deleteat!(FV,10)
# FVsurf = cat([FV[1:9], FV[11:end]])
EV = LarSurf.Smoothing.get_EV_quads(FV;return_unfiltered=true)
function ev_half(EV)
    sortedEV = sort(EV)
    EV = LarSurf.Smoothing.get_EV_quads(FV;return_unfiltered=false)
    return EV
end

function ev_set(EV)
    EV = collect(Set(EV))
    return EV
end
@benchmark ev_set(EV)
@benchmark ev_half(EV)


text:
# PDB Types
# =========

@auto_hash_equals immutable PDBResidueIdentifier
    PDBe_number::ASCIIString # PDBe
    number::ASCIIString # PDB
    name::ASCIIString
    group::ASCIIString
    model::ASCIIString
    chain::ASCIIString
end

@auto_hash_equals immutable Coordinates <: FixedVectorNoTuple{3, Float64}
    x::Float64
    y::Float64
    z::Float64
    function Coordinates(a::NTuple{3, Float64})
        new(a[1], a[2], a[3])
    end
end

@auto_hash_equals immutable PDBAtom
    coordinates::Coordinates
    atom::ASCIIString
    element::ASCIIString
    occupancy::Float64
    B::ASCIIString
end

@auto_hash_equals type PDBResidue
    id::PDBResidueIdentifier
    atoms::Vector{PDBAtom}
end

length(res::PDBResidue) = length(res.atoms)

# Coordinates
# ===========

vec(a::Coordinates) = Float64[a.x, a.y, a.z]

# Distances and geometry
# ----------------------

distance(a::Coordinates, b::Coordinates) = sqrt((a.x - b.x)^2 + (a.y - b.y)^2 + (a.z - b.z)^2)

"""
`contact(a::Coordinates, b::Coordinates, limit::AbstractFloat)`

This simply returns `distance(a,b) <= limit`.
"""
contact(a::Coordinates, b::Coordinates, limit::AbstractFloat) = distance(a,b) <= limit

"""
`angle(a::Coordinates, b::Coordinates, c::Coordinates)`

Angle (in degrees) at `b` between `a-b` and `b-c`
"""
function angle(a::Coordinates, b::Coordinates, c::Coordinates)
    A = b - a
    B = b - c
    norms = (norm(A)*norm(B))
    if norms != 0
        return( acosd(dot(A,B) / norms) )
    else
        return(0.0)
    end
end

distance(a::PDBAtom, b::PDBAtom) = distance(a.coordinates, b.coordinates)

contact(a::PDBAtom, b::PDBAtom, limit::Float64) = contact(a.coordinates, b.coordinates, limit)

angle(a::PDBAtom, b::PDBAtom, c::PDBAtom) = angle(a.coordinates, b.coordinates, c.coordinates)

cross(a::PDBAtom, b::PDBAtom) = cross(a.coordinates, b.coordinates)

# Find Residues/Atoms
# ===================

isobject(res::PDBResidue, tests::AbstractTest...) = isobject(res.id, tests...)

findobjects(res::PDBResidue, tests::AbstractTest...) = findobjects(res.atoms, tests...)

# @residues
# =========

_test_stringfield(field::Symbol, test::Real) = Is(field, string(test))
_test_stringfield(field::Symbol, test::Union{Char, Symbol}) = Is(field, string(test))
_test_stringfield(field::Symbol, test::Union{ASCIIString, Regex, Function}) = Is(field, test)
_test_stringfield(field::Symbol, test::Union{UnitRange, IntSet, Set, Array, Base.KeyIterator}) = In(field, test)

_is_wildcard(test::ASCIIString) = test == "*"
_is_wildcard(test::Char) = test == '*'
_is_wildcard(test) = false

function _add_test_stringfield!(tests::Vector{TestType}, field::Symbol, test)
    if !_is_wildcard(test)
        push!(tests, _test_stringfield(field, test))
    end
    tests
end

function _residues_tests(model, chain, group, residue)
    args = TestType[]
    _add_test_stringfield!(args, :model, model)
    _add_test_stringfield!(args, :chain, chain)
    _add_test_stringfield!(args, :group, group)
    _add_test_stringfield!(args, :number, residue)
    args
end

"""
`residues(residue_list, model, chain, group, residue)`

These return a new vector with the selected subset of residues from a list of residues.
"""
residues(residue_list, model, chain, group, residue) = collectobjects(residue_list, _residues_tests(model, chain, group, residue)...)

"""
`@residues ... model ... chain ... group ... residue ...`

These return a new vector with the selected subset of residues from a list of residues.
"""
macro residues(residue_list,
               model::Symbol, m, #::Union(Int, Char, ASCIIString, Symbol),
               chain::Symbol, c, #::Union(Char, ASCIIString, Symbol),
               group::Symbol, g,
               residue::Symbol, r)
    if model == :model && chain == :chain && group == :group && residue == :residue
        return :(residues($(esc(residue_list)), $(esc(m)), $(esc(c)), $(esc(g)), $(esc(r))))
    else
        throw(ArgumentError("The signature is @residues ___ model ___ chain ___ group ___ residue ___"))
    end
end

"""
`residuesdict(residue_list, model, chain, group, residue)`

These return a dictionary (using PDB residue numbers as keys) with the selected subset of residues.
"""
function residuesdict(residue_list, model, chain, group, residue)
    res_index = findobjects(residue_list, _residues_tests(model, chain, group, residue)...)
    dict = sizehint!( OrderedDict{ASCIIString, PDBResidue}() , length(res_index) )
    for i in res_index
        dict[ residue_list[i].id.number ] = residue_list[i]
    end
    dict
end

"""
`@residuesdict ... model ... chain ... group ... residue ...`

These return a dictionary (using PDB residue numbers as keys) with the selected subset of residues.
"""
macro residuesdict(residue_list,
                   model::Symbol, m, #::Union(Int, Char, ASCIIString, Symbol),
                   chain::Symbol, c, #::Union(Char, ASCIIString, Symbol),
                   group::Symbol, g,
                   residue::Symbol, r)
    if model == :model && chain == :chain && group == :group && residue == :residue
        return :(residuesdict($(esc(residue_list)), $(esc(m)), $(esc(c)), $(esc(g)), $(esc(r))))
    else
        throw(ArgumentError("The signature is @residues ___ model ___ chain ___ group ___ residue ___"))
    end
end

# @atoms
# ======

"""
`atoms(residue_list, model, chain, group, residue, atom)`

These return a vector of `PDBAtom`s with the selected subset of atoms.
"""
function atoms(residue_list, model, chain, group, residue, atom)
    _is_wildcard(atom) ?
        vcat(Vector{PDBAtom}[ res.atoms for res in residues(residue_list, model, chain, group, residue) ]...) :
        vcat(Vector{PDBAtom}[ collectobjects(res.atoms, _test_stringfield(:atom, atom)) for res in residues(residue_list, model, chain, group, residue) ]...)
end

"""
`@atoms ... model ... chain ... group ... residue ... atom ...`

These return a vector of `PDBAtom`s with the selected subset of atoms.
"""
macro atoms(residue_list,
            model::Symbol, m, #::Union(Int, Char, ASCIIString, Symbol),
            chain::Symbol, c, #::Union(Char, ASCIIString, Symbol),
            group::Symbol, g,
            residue::Symbol, r,
            atom::Symbol, a)
    if model == :model && chain == :chain && group == :group && residue == :residue && atom == :atom
        return :(atoms($(esc(residue_list)), $(esc(m)), $(esc(c)), $(esc(g)), $(esc(r)), $(esc(a))))
    else
        throw(ArgumentError("The signature is @atoms ___ model ___ chain ___ group ___ residue ___ atom ___"))
    end
end

# Special find...
# ===============

"Returns a list with the index of the heavy atoms (all atoms except hydrogen) in the `PDBResidue`"
function findheavy(res::PDBResidue)
    N = length(res)
    indices = Array(Int,N)
    j = 0
    @inbounds for i in 1:N
        if res.atoms[i].element != "H"
            j += 1
            indices[j] = i
        end
    end
    resize!(indices, j)
end

"""
`findatoms(res::PDBResidue, atom::ASCIIString)`

Returns a index vector of the atoms with the given `atom` name.
"""
function findatoms(res::PDBResidue, atom::ASCIIString)
    N = length(res)
    indices = Array(Int,N)
    j = 0
    @inbounds for i in 1:N
        if res.atoms[i].atom == atom
            j += 1
            indices[j] = i
        end
    end
    resize!(indices, j)
end

"Returns a vector of indices for `CB` (`CA` for `GLY`)"
function findCB(res::PDBResidue)
    N = length(res)
    indices = Array(Int,N)
    atom = res.id.name == "GLY" ? "CA" : "CB"
    j = 0
    @inbounds for i in 1:N
        if res.atoms[i].atom == atom
            j += 1
            indices[j] = i
        end
    end
    resize!(indices, j)
end

# occupancy
# =========

"""
Takes a `PDBResidue` and a `Vector` of atom indices.
Returns the index value of the `Vector` with maximum occupancy.
"""
function selectbestoccupancy(res::PDBResidue, indices::Vector{Int})
    Ni = length(indices)
    if Ni == 1
        return(indices[1])
    end
    Na = length(res)
    if Ni == 0 || Ni > Na
        throw(ErrorException("There are no atom indices or they are more atom indices ($Ni) than atoms in the Residue ($Na)"))
    end
    indice = indices[1]
    occupancy = 0.0
    for i in indices
        actual_occupancy = res.atoms[i].occupancy
        if actual_occupancy > occupancy
            occupancy = actual_occupancy
            indice = i
        end
    end
    return(indice)
end

"""
Takes a `Vector` of `PDBAtom`s and returns a `Vector` of the `PDBAtom`s with best occupancy.
"""
function bestoccupancy(atoms::Vector{PDBAtom})
    N = length(atoms)
    if N == 0
        warn("There are no atoms.")
        return(atoms)
    elseif N == 1
        return(atoms)
    else
        atomdict = sizehint!(OrderedDict{ASCIIString,PDBAtom}(), N)
        for atom in atoms
            name = atom.atom
            if haskey(atomdict, name)
                if atom.occupancy > atomdict[name].occupancy
                    atomdict[name] = atom
                end
            else
                atomdict[name] = atom
            end
        end
        return(collect(values(atomdict)))
    end
end

# More Distances
# ==============

function _update_distance(a, b, i, j, dist)
    actual_dist = distance(a.atoms[i], b.atoms[j])
    if actual_dist < dist
        return(actual_dist)
    else
        return(dist)
    end
end

"""
`distance(a::PDBResidue, b::PDBResidue; criteria::ASCIIString="All")`

Returns the distance between the residues `a` and `b`.
The available `criteria` are: `Heavy`, `All`, `CA`, `CB` (`CA` for `GLY`)
"""
function distance(a::PDBResidue, b::PDBResidue; criteria::ASCIIString="All")
    dist = Inf
    if criteria == "All"
        Na = length(a)
        Nb = length(b)
        @inbounds for i in 1:Na
            for j in 1:Nb
                dist = _update_distance(a, b, i, j, dist)
            end
        end
    elseif criteria == "Heavy"
        indices_a = findheavy(a)
        indices_b = findheavy(b)
        if length(indices_a) != 0 && length(indices_b) != 0
            for i in indices_a
                for j in indices_b
                    dist = _update_distance(a, b, i, j, dist)
                end
            end
        end
    elseif criteria == "CA"
        indices_a = findatoms(a, "CA")
        indices_b = findatoms(b, "CA")
        if length(indices_a) != 0 && length(indices_b) != 0
            for i in indices_a
                for j in indices_b
                    dist = _update_distance(a, b, i, j, dist)
                end
            end
        end
    elseif criteria == "CB"
        indices_a = findCB(a)
        indices_b = findCB(b)
        if length(indices_a) != 0 && length(indices_b) != 0
            for i in indices_a
                for j in indices_b
                    dist = _update_distance(a, b, i, j, dist)
                end
            end
        end
    end
    dist
end

"""
`any(f::Function, a::PDBResidue, b::PDBResidue)`

Test if the function `f` is true for any pair of atoms between the residues `a` and `b`
"""
function any(f::Function, a::PDBResidue, b::PDBResidue)
    Na = length(a)
    Nb = length(b)
    @inbounds for i in 1:Na
        for j in 1:Nb
            if f(a.atoms[i], b.atoms[j])
                return(true)
            end
        end
    end
    return(false)
end

"""
`contact(a::PDBResidue, b::PDBResidue, limit::AbstractFloat; criteria::ASCIIString="All")`

Returns `true` if the residues `a` and `b` are at contact distance (`limit`).
The available distance `criteria` are: `Heavy`, `All`, `CA`, `CB` (`CA` for `GLY`)
"""
function contact(a::PDBResidue, b::PDBResidue, limit::AbstractFloat; criteria::ASCIIString="All")
    if criteria == "All"
        Na = length(a)
        Nb = length(b)
        @inbounds for i in 1:Na
            for j in 1:Nb
                if contact(a.atoms[i], b.atoms[j], limit)
                    return(true)
                end
            end
        end
    elseif criteria == "Heavy"
        indices_a = findheavy(a)
        indices_b = findheavy(b)
        if length(indices_a) != 0 && length(indices_b) != 0
            for i in indices_a
                for j in indices_b
                    if contact(a.atoms[i], b.atoms[j], limit)
                        return(true)
                    end
                end
            end
        end
    elseif criteria == "CA"
        indices_a = findatoms(a, "CA")
        indices_b = findatoms(b, "CA")
        if length(indices_a) != 0 && length(indices_b) != 0
            for i in indices_a
                for j in indices_b
                    if contact(a.atoms[i], b.atoms[j], limit)
                        return(true)
                    end
                end
            end
        end
    elseif criteria == "CB"
        indices_a = findCB(a)
        indices_b = findCB(b)
        if length(indices_a) != 0 && length(indices_b) != 0
            for i in indices_a
                for j in indices_b
                    if contact(a.atoms[i], b.atoms[j], limit)
                        return(true)
                    end
                end
            end
        end
    end
    false
end

# Vectorize
# ---------

"""
`contact(vec::Vector{PDBResidue}, limit::AbstractFloat; criteria::ASCIIString="All")`

If `contact` takes a `Vector{PDBResidue}`, It returns a `BitMatrix` with all the pairwise comparisons (contact map).
"""
function contact(vec::Vector{PDBResidue}, limit::AbstractFloat; criteria::ASCIIString="All")
    N = length(vec)
    cmap = trues(N,N)
    for i in 1:(N-1)
        for j in (i+1):N
            cmap[i, j] = cmap[j, i] = contact(vec[i], vec[j], limit, criteria=criteria)
        end
    end
    cmap
end

"""
`distance(vec::Vector{PDBResidue}; criteria::ASCIIString="All")`

If `distance` takes a `Vector{PDBResidue}` returns a `PairwiseListMatrix{Float64, false}` with all the pairwise comparisons (distance matrix).
"""
function distance(vec::Vector{PDBResidue}; criteria::ASCIIString="All")
    PLM = PairwiseListMatrix(Float64, length(vec), false, 0.0)
    @iterateupper PLM false list[k] = :($distance)( :($vec)[i], :($vec)[j], criteria=:($criteria))
    PLM
end

# Proximity average
# -----------------

"""
`proximitymean` calculates the proximity mean/average for each residue as the average score (from a `scores` list)
of all the residues within a certain physical distance to a given amino acid.
The score of that residue is not included in the mean unless you set `include` to `true`.
The default values are 6.05 for the distance threshold/`limit` and "Heavy" for the `criteria` keyword argument.
This function allows to calculate pMI (proximity mutual information) and pC (proximity conservation) as in Buslje et. al. 2010.

Buslje, Cristina Marino, Elin Teppa, Tomas Di Doménico, José María Delfino, and Morten Nielsen.
*Networks of high mutual information define the structural proximity of catalytic sites: implications for catalytic residue identification.*
PLoS Comput Biol 6, no. 11 (2010): e1000978.
"""
function proximitymean{T}(residues::Vector{PDBResidue}, scores::AbstractVector{T}, limit::AbstractFloat=6.05;
                          criteria::ASCIIString="Heavy", include::Bool=false)
    N = length(residues)
    if N != length(scores)
        throw(ErrorException("Vectors must have the same length."))
    end
    count = zeros(Int, N)
    sum   = zeros(T, N)
    offset = include ? 0 : 1
    @inbounds for i in 1:(N-offset)
        res_i = residues[i]
        for j in (i+offset):N
            if include && (i == j)
                count[i] += 1
                sum[i] += scores[i]
            elseif contact(res_i, residues[j], limit, criteria=criteria)
                count[i] += 1
                count[j] += 1
                sum[i] += scores[j]
                sum[j] += scores[i]
            end
        end
    end
    sum ./ count
end

# For Aromatic
# ============

function _get_plane(residue::PDBResidue)
    name = residue.id.name
    planes = Vector{PDBAtom}[]
    if name != "TRP"
        plane = PDBAtom[]
        for atom in residue.atoms
            if (name, atom.atom) in PDB._aromatic
                push!(plane, atom)
            end
        end
        push!(planes, plane)
    else
        plane1 = PDBAtom[]
        plane2 = PDBAtom[]
        for atom in residue.atoms
            if atom.atom in Set{ASCIIString}(["CE2", "CD2", "CZ2", "CZ3", "CH2", "CE3"])
                push!(plane1, atom)
            end
            if atom.atom in Set{ASCIIString}(["CG", "CD1", "NE1", "CE2", "CD2"])
                push!(plane2, atom)
            end
        end
        push!(planes, plane1)
        push!(planes, plane2)
    end
    planes
end

function _centre(planes::Vector{Vector{PDBAtom}})
    subset = Vector{PDBAtom}[ bestoccupancy(atoms) for atoms in planes ]
    polyg = Vector{Coordinates}[ Coordinates[ a.coordinates for a in atoms ] for atoms in subset ]
    Coordinates[ sum(points)./Float64(length(points)) for points in polyg ]
end

# function _simple_normal_and_centre(atoms::Vector{PDBAtom})
#   atoms = bestoccupancy(atoms)
#   points = Coordinates[ a.coordinates for a in atoms ]
#   (cross(points[2] - points[1], points[3] - points[1]), sum(points)./length(points))
# end

# For Parsing...
# ==============

"Used for parsing a PDB file into `Vector{PDBResidue}`"
function _generate_residues(residue_dict::OrderedDict{PDBResidueIdentifier, Vector{PDBAtom}}, occupancyfilter::Bool=false)
    occupancyfilter ? PDBResidue[ PDBResidue(k, bestoccupancy(v)) for (k,v) in residue_dict ] : PDBResidue[ PDBResidue(k, v) for (k,v) in residue_dict ]
end

# Show PDB* objects (using Formatting)
# ====================================

const _Format_ResidueID = FormatExpr("{:>15} {:>15} {:>15} {:>15} {:>15} {:>15}\n")
const _Format_ATOM = FormatExpr("{:>50} {:>15} {:>15} {:>15} {:>15}\n")

const _Format_ResidueID_Res = FormatExpr("\t\t{:>15} {:>15} {:>15} {:>15} {:>15} {:>15}\n")
const _Format_ATOM_Res = FormatExpr("\t\t{:<10} {:>50} {:>15} {:>15} {:>15} {:>15}\n")

function show(io::IO, id::PDBResidueIdentifier)
    printfmt(io, _Format_ResidueID, "PDBe_number", "number", "name", "group", "model", "chain")
    printfmt(io, _Format_ResidueID, string('"',id.PDBe_number,'"'), string('"',id.number,'"'), string('"',id.name,'"'),
             string('"',id.group,'"'), string('"',id.model,'"'), string('"',id.chain,'"'))
end

function show(io::IO, atom::PDBAtom)
    printfmt(io, _Format_ATOM, "coordinates", "atom", "element", "occupancy", "B")
    printfmt(io, _Format_ATOM, atom.coordinates, string('"',atom.atom,'"'), string('"',atom.element,'"'),
             atom.occupancy, string('"',atom.B,'"'))
end

function show(io::IO, res::PDBResidue)
    println(io, "PDBResidue:\n\tid::PDBResidueIdentifier")
    printfmt(io, _Format_ResidueID_Res, "PDBe_number", "number", "name", "group", "model", "chain")
    printfmt(io, _Format_ResidueID_Res, string('"',res.id.PDBe_number,'"'), string('"',res.id.number,'"'),
             string('"',res.id.name,'"'), string('"',res.id.group,'"'), string('"',res.id.model,'"'), string('"',res.id.chain,'"'))
    len = length(res)
    println(io, "\tatoms::Vector{PDBAtom}\tlength: ", len)
    for i in 1:len
        printfmt(io, _Format_ATOM_Res, "", "coordinates", "atom", "element", "occupancy", "B")
        printfmt(io, _Format_ATOM_Res, string(i,":"), res.atoms[i].coordinates, string('"',res.atoms[i].atom,'"'),
                 string('"',res.atoms[i].element,'"'), res.atoms[i].occupancy, string('"',res.atoms[i].B,'"'))
    end
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# Creates cache for time series callback
function create_cache_time_series(point_coordinates, mesh::TreeMesh{2}, dg, cache)
  # Determine element ids for point coordinates
  element_ids = get_elements_by_coordinates(point_coordinates, mesh, dg, cache)

  # Calculate & store Lagrange interpolation polynomials
  interpolating_polynomials = calc_interpolating_polynomials(point_coordinates, element_ids, mesh,
                                                             dg, cache)

  time_series_cache = (; element_ids, interpolating_polynomials)

  return time_series_cache
end

# Find element ids containing coordinates given as a matrix [ndims, npoints]
function get_elements_by_coordinates!(element_ids, coordinates, mesh::TreeMesh, dg, cache)
  if length(element_ids) != size(coordinates, 2)
    throw(DimensionMismatch("storage length for element ids does not match the number of coordinates"))
  end

  @unpack cell_ids = cache.elements
  @unpack tree = mesh

  # Reset element ids - 0 indicates "not (yet) found"
  element_ids .= 0
  found_elements = 0

  # Iterate over all elements
  for element in eachelement(dg, cache)
    # Get cell id
    cell_id = cell_ids[element]

    # Iterate over coordinates
    for index in 1:length(element_ids)
      # Skip coordinates for which an element has already been found
      if element_ids[index] > 0
        continue
      end

      # Construct point
      x = SVector(ntuple(i -> coordinates[i, index], ndims(mesh)))

      # Skip if point is not in cell
      if !is_point_in_cell(tree, x, cell_id)
        continue
      end

      # Otherwise point is in cell and thus in element
      element_ids[index] = element
      found_elements += 1
    end

    # Exit loop if all elements have already been found
    if found_elements == length(element_ids)
      break
    end
  end

  return element_ids
end

function get_elements_by_coordinates(coordinates, mesh, dg, cache)
  element_ids = Vector{Int}(undef, size(coordinates, 2))
  get_elements_by_coordinates!(element_ids, coordinates, mesh, dg, cache)

  return element_ids
end

# Calculate the interpolating polynomials to extract data at the given coordinates
# The coordinates are known to be located in the respective element in `element_ids`
function calc_interpolating_polynomials!(interpolating_polynomials, coordinates, element_ids,
                                         mesh::TreeMesh, dg::DGSEM, cache)
  @unpack tree = mesh
  @unpack nodes = dg.basis

  wbary = barycentric_weights(nodes)

  for index in 1:length(element_ids)
    # Construct point
    x = SVector(ntuple(i -> coordinates[i, index], ndims(mesh)))

    # Convert to unit coordinates
    cell_id = cache.elements.cell_ids[element_ids[index]]
    cell_coordinates_ = cell_coordinates(tree, cell_id)
    cell_length = length_at_cell(tree, cell_id)
    unit_coordinates = (x .- cell_coordinates_) * 2 / cell_length

    # Calculate interpolating polynomial for each dimension, making use of tensor product structure
    for d in 1:ndims(mesh)
      interpolating_polynomials[:, d, index] .= lagrange_interpolating_polynomials(
          unit_coordinates[d], nodes, wbary)
    end
  end

  return interpolating_polynomials
end

function calc_interpolating_polynomials(coordinates, element_ids, mesh::TreeMesh, dg, cache)
  interpolating_polynomials = Array{real(dg), 3}(undef,
                                                 nnodes(dg), ndims(mesh), length(element_ids))
  calc_interpolating_polynomials!(interpolating_polynomials, coordinates, element_ids, mesh, dg,
                                  cache)

  return interpolating_polynomials
end

# Record the solution variables at each given point
function record_state_at_points!(point_data, u, solution_variables, n_solution_variables,
                                 mesh::TreeMesh{2}, equations, dg::DG, time_series_cache)
  @unpack element_ids, interpolating_polynomials = time_series_cache
  old_length = length(first(point_data))
  new_length = old_length + n_solution_variables

  # Loop over all points/elements that should be recorded
  for index in 1:length(element_ids)
    # Extract data array and element id
    data = point_data[index]
    element_id = element_ids[index]

    # Make room for new data to be recorded
    resize!(data, new_length)
    data[(old_length+1):new_length] .= zero(eltype(data))

    # Loop over all nodes to compute their contribution to the interpolated values
    for j in eachnode(dg), i in eachnode(dg)
      u_node = solution_variables(get_node_vars(u, equations, dg, i, j, element_id), equations)

      for v in 1:length(u_node)
        data[old_length + v] += (u_node[v]
                                * interpolating_polynomials[i, 1, index]
                                * interpolating_polynomials[j, 2, index])
      end
    end
  end
end

end # @muladd


text:
"""

"""
function corey_gas_oil()

end # corey_gas_oil

"""
brooks_corey_gas_oil(labda=4.17, sor=0.2, pce_go=5000)

Original Brooks and Corey (1964)

Input parameters:

    + labda: sorting factor (pore-size distribution)
    + sor: residual oil saturation
    + pce_go: gas-oil pore entry capillary pressure
"""
function brooks_corey_gas_oil(labda=4.17, sor=0.2, pce_go=5000)
  kro(so)=((so-sor)/(1-sor)).^((2+3*labda)/labda)
  krg(so)=((1-so)/(1-sor)).^2.*(1-((so-sor)/(1-sor)).^((2+labda)/labda))
  pc_go(so)=pce_go*((so-sor)/(1-sor)).^(-1/labda)
  return kro, krg, pc_go
end # brooks_corey

"""

"""
function brooks_corey_general(kro_max=1.0, krg_max=1.0, ng=2.0, no=2.0)

end # brooks_corey

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
function gen_pawn_one_step_valid(source_square::UInt64,
    color::String="white")

    if color == "white"
        return source_square << 8
    else
        return source_square >> 8
    end
end

function gen_all_pawns_one_step_valid_moves(color::String="white")
    pawn_moves = Dict{UInt64, UInt64}()

    for i in 1:8
        pawn_moves[INT2UINT[i]] = EMPTY
    end

    for i in 9:56
        pawn_moves[INT2UINT[i]] = gen_pawn_one_step_valid(INT2UINT[i],
            color)
    end

    for i in 57:64
        pawn_moves[INT2UINT[i]] = EMPTY
    end

    return pawn_moves
end
const WHITE_PAWN_ONESTEP_MOVES = gen_all_pawns_one_step_valid_moves()
const BLACK_PAWN_ONESTEP_MOVES = gen_all_pawns_one_step_valid_moves("black")

function gen_pawn_two_steps_valid(source_square::UInt64,
    color::String="white")

    if color == "white"
        return source_square << 16
    else
        return source_square >> 16
    end
end

function gen_all_pawns_two_steps_valid_moves(color::String="white")
    pawn_moves = Dict{UInt64, UInt64}()

    if color == "white"
        for i in 1:48
            pawn_moves[INT2UINT[i]] = EMPTY
        end

        for i in 49:56
            pawn_moves[INT2UINT[i]] = gen_pawn_two_steps_valid(
                INT2UINT[i])
        end

        for i in 57:64
            pawn_moves[INT2UINT[i]] = EMPTY
        end
    else
        for i in 1:8
            pawn_moves[INT2UINT[i]] = EMPTY
        end

        for i in 9:16
            pawn_moves[INT2UINT[i]] = gen_pawn_two_steps_valid(
                INT2UINT[i], color)
        end

        for i in 17:64
            pawn_moves[INT2UINT[i]] = EMPTY
        end
    end

    return pawn_moves
end
const WHITE_PAWN_TWOSTEPS_MOVES = gen_all_pawns_two_steps_valid_moves()
const BLACK_PAWN_TWOSTEPS_MOVES = gen_all_pawns_two_steps_valid_moves("black")

function gen_pawn_attacked_valid(source_square::UInt64,
    color::String="white")

    target_squares = zeros(UInt64, 0)

    if color == "white"
        if source_square & CLEAR_FILE_A != EMPTY
            push!(target_squares, source_square << 9)
        end
        if source_square & CLEAR_FILE_H != EMPTY
            push!(target_squares, source_square << 7)
        end
    else
        if source_square & CLEAR_FILE_H != EMPTY
            push!(target_squares, source_square >> 9)
        end
        if source_square & CLEAR_FILE_A != EMPTY
            push!(target_squares, source_square >> 7)
        end
    end

    return target_squares
end

function gen_all_pawns_valid_attack(color::String="white")
    pawn_moves = Dict{UInt64, Array{UInt64,1}}()
    for i in 1:64
        pawn_moves[INT2UINT[i]] = gen_pawn_attacked_valid(INT2UINT[i], color)
    end
    return pawn_moves
end
const WHITE_PAWN_ATTACK = gen_all_pawns_valid_attack()
const BLACK_PAWN_ATTACK = gen_all_pawns_valid_attack("black")

function add_promotions(pawns_moves::Array{Move,1}, source::UInt64,
    target::UInt64, taken_piece::String="none")

    for promotion_type in ["queen", "rook", "night", "bishop"]
        push!(pawns_moves, Move(source, target, "pawn", taken_piece,
            promotion_type, EMPTY, "-"))
    end

    return pawns_moves
end

function find_pawn_pseudo!(pawns_moves::Array{Move,1}, board::Bitboard,
    pawn::UInt64, color::String="white")

    if color == "white"
        home_rank = MASK_RANK_2
        promotion_rank = MASK_RANK_7
        one_step = WHITE_PAWN_ONESTEP_MOVES
        two_steps = WHITE_PAWN_TWOSTEPS_MOVES
        attacks = WHITE_PAWN_ATTACK
        others = board.black
        enpassant_rank = MASK_RANK_6
        opponent_color = "black"
    else
        home_rank = MASK_RANK_7
        promotion_rank = MASK_RANK_2
        one_step = BLACK_PAWN_ONESTEP_MOVES
        two_steps = BLACK_PAWN_TWOSTEPS_MOVES
        attacks = BLACK_PAWN_ATTACK
        others = board.white
        enpassant_rank = MASK_RANK_3
        opponent_color = "white"
    end

    # move once
    if one_step[pawn] & board.taken == EMPTY # check front square
        if pawn & promotion_rank != EMPTY
            pawns_moves = add_promotions(pawns_moves, pawn,
                one_step[pawn], "none")
        else
            # move once
            push!(pawns_moves, Move(pawn, one_step[pawn], "pawn", "none",
                "none", EMPTY, "-"))

            # move twice and take note of the en-passant square
            if pawn & home_rank != EMPTY
                if two_steps[pawn] & board.taken == EMPTY
                    push!(pawns_moves, Move(pawn, two_steps[pawn],
                        "pawn", "none", "none", one_step[pawn], "-"))
                end
            end
        end
    end

    # attack squares
    pawn_attacks = attacks[pawn]
    for attack in pawn_attacks
        if attack & others != EMPTY
            taken_piece = find_piece_type(board, attack, opponent_color)
            if pawn & promotion_rank != EMPTY
                pawns_moves = add_promotions(pawns_moves, pawn,
                    attack, taken_piece)
            else
                push!(pawns_moves, Move(pawn, attack, "pawn", taken_piece,
                    "none", EMPTY, "-"))
            end
        elseif attack == board.enpassant_square &&
            attack & enpassant_rank != EMPTY # use enpassant square
            push!(pawns_moves, Move(pawn, attack, "pawn", "none", "none",
                EMPTY, "-"))
        end
    end
    return
end

function get_pawns_list(pawns_moves::Array{Move,1}, board::Bitboard,
    color::String="white")

    if color == "white"
        pawns = board.P
    else
        pawns = board.p
    end

    # pawns_moves = Array{Move,1}()
    for pawn in pawns
        find_pawn_pseudo!(pawns_moves, board, pawn, color)
    end

    return pawns_moves
end

function get_pawns_list(board::Bitboard, color::String="white")

    if color == "white"
        pawns = board.P
    else
        pawns = board.p
    end

    pawns_moves = Array{Move,1}()
    for pawn in pawns
        find_pawn_pseudo!(pawns_moves, board, pawn, color)
    end

    return pawns_moves
end


text:
using DelimitedFiles
using Downloads
using Tar
using CodecZlib

# ---------------------------------------------------------------------------------------- #
# MISC CONVERSION UTILITIES

if !(@isdefined GREEKS)
    const GREEKS = Dict("GAMMA" => 'Γ', "LAMBDA" => 'Λ', "SIGMA" => 'Σ', "DELTA"=>'Δ')
end
if !(@isdefined SUBSCRIPTS)
    const SUBSCRIPTS = Dict('0'=>'₀', '1'=>'₁', '2'=>'₂', '3'=>'₃', '4'=>'₄',
                            '5'=>'₅', '6'=>'₆', '7'=>'₇', '8'=>'₈', '9'=>'₉')
end

function greekify(s::String)
    for (key, substitute) in GREEKS
        if occursin(key, s)
            return replace(s, key=>substitute)
        end
    end
    return s
end

function subscriptify(s::String)
    m = match(r"_[0-9]", s)
    m === nothing && return s
    sub    = m.match
    newsub = SUBSCRIPTS[sub[end]]
    return replace(s, sub => newsub)
end

# ---------------------------------------------------------------------------------------- #
# DOWNLOAD THE SEEK-PATH REPO; FIND DIRECTORY WITH BAND_PATH_DATA
user = "giovannipizzi"
repo = "seekpath"
url  = "https://github.com/$(user)/$(repo)/tarball/master" # tarball url

# create a temporary file path to download in and then manipulate the data
mktemp(mktempdir()) do path, io
    # download a (shallow) tarball version of SeeK's source code
    Downloads.download(url, io)
    # extract compressed tarball to to tmpdir/ using Tar.jl and CodecZlib.jl
    tmpdir = open(path) do targz_io
        tar = GzipDecompressorStream(targz_io)
        Tar.extract(tar)
    end

    # find the extracted content, which is always of the form $user-$repo-$hash
    contents = readdir(tmpdir)
    i = findfirst(str -> contains(str, user), contents)::Int
    localdir = contents[i]
    basedir  = joinpath(tmpdir, localdir) # path to "root" of repo directory

    # directory containing all the extended Bravais type data we want
    datadir  = joinpath(basedir, "seekpath", "hpkot", "band_path_data")

    # license file
    license_file = joinpath(basedir, "LICENSE.txt")

    # ------------------------------------------------------------------------------------ #
    # EXTRACT ALL EXTENDED BRAVAIS TYPES
    bravaistypes = Symbol.(filter!(s->endswith(s, r"[1-3]"), readdir(datadir)))

    pathsd  = Dict{Symbol, Vector{Vector{Symbol}}}()
    pointsd = Dict{Symbol, Vector{Pair{Symbol, Expr}}}()
    paramsd = Dict{Symbol, Vector{Pair{Symbol, Expr}}}()

    for bt in bravaistypes
        btstr = string(bt)
        # -------------------------------------------------------------------------------- #
        # READ PATHS/SEGMENTS
        raw_paths = readdlm(joinpath(datadir, btstr, "path.txt"), String)
        # adjust formatting
        raw_paths .= greekify.(raw_paths)
        raw_paths .= subscriptify.(raw_paths)
        raw_paths = Symbol.(raw_paths)

        # split into connected paths
        pathsd[bt] = [[(first(raw_paths))]]
        last_stop = first(raw_paths)
        i = 1
        for (start, stop) in eachrow(raw_paths)
            if start == last_stop
                push!(pathsd[bt][i], stop) # extend path
            else
                push!(pathsd[bt], [start, stop]) # start new path
                i += 1 # continue on new path
            end
            last_stop = stop
        end

        # -------------------------------------------------------------------------------- #
        # READ POINTS
        raw_points = readdlm(joinpath(datadir, btstr, "points.txt"), String)
        labels = Symbol.(subscriptify.(greekify.(raw_points[:,1])))
        values_str  = "[".*join.(eachrow(raw_points[:,2:4]), ',').*']' # as `String`s
        values_repr = Meta.parse.(values_str) # as `Expr`s
        pointsd[bt] = [label => value for (label, value) in zip(labels, values_repr)]

        # -------------------------------------------------------------------------------- #
        # READ VARIABLES/PARAMETERS, IF ANY
        param_filepath = joinpath(datadir, btstr, "k_vector_parameters.txt")
        (read(param_filepath, String) == "") && continue # for empty files, nothing to add
        raw_params = readdlm(param_filepath, String)
        paramsd[bt] = [Symbol(label) => Meta.parse(replace(value, "beta"=>'β'))
                                                    for (label, value) in eachrow(raw_params)]
    end

    # ------------------------------------------------------------------------------------ #
    # FILTER OUT POINTS FROM `pointsd` THAT ARE NOT REFERENCED IN `pathsd`
    # NB: For a few extended Bravais types, SeeK includes points that are not referenced
    #     in the paths; I haven't found a reason/use-case for these extra points, so it
    #     seems easier just to get rid of them entirely (then we can e.g. guarantee that
    #     every point point will be referenced)
    for bt in keys(pointsd)
        points = pointsd[bt]
        paths  = pathsd[bt]
        remove_idxs = Int[]
        for (idx, (klab, kvexpr)) in enumerate(points)
            if klab ∉ Iterators.flatten(paths)
                push!(remove_idxs, idx)
            end
        end
        deleteat!(points, remove_idxs)
    end

    # ------------------------------------------------------------------------------------ #
    # WRITE DICTIONARIES TO A .JL FILE; USE JULIA'S PARSER AS FILE FORMAT
    data_dir = joinpath(@__DIR__, "..", "assets")
    isdir(data_dir) || mkpath(data_dir)
    open(joinpath(data_dir, "data-SeeK.jl"), "w") do io
        println(io,
        "# ", "-"^76, "\n",
        """
        # This file was autogenerated by /$(last(splitpath(@__DIR__)))/$(basename(@__FILE__)).
        # It contains data included in the SeeK-path package (github.com/$user/$repo).
        # The SeeK-path package is licensed under the MIT license, a copy of which is
        # included below.
        """,
        "# ", "-"^76, "\n"
        )

        for name in (:pathsd, :pointsd, :paramsd)
            # obtain the variable via its name (can't use `getfield(Main, name)` since this
            # is a local scope)
            d = Base.@locals()[name]

            print(io, "const ", name, "_3d = Dict(")
            for (key, val) in d
                println(io)
                print(io, "   :", key, " => ")
                if val isa Vector{Pair{Symbol, Expr}}
                    print(io, '[')
                    for (i, (sym, ex)) in enumerate(val)
                        print(io, ':', sym, " => :(", ex, ")")
                        i != length(val) && print(io, ", ")
                    end
                    print(io, "],")
                else
                    print(io, val, ',')
                end
            end
            println(io, "\n)\n")
        end

        # copy over licence information from SeeK-path
        license_lines = readlines(license_file)
        println(io, "# ------------------ COPY OF ", uppercase(user*"/"*repo),
                    " LICENSE ------------------")
        foreach(line -> println(io, "# ", line), license_lines)
        print(io, "# ", "-"^76)
    end
end
# ---------------------------------------------------------------------------------------- #


text:
using RMInsurance

using Distributions
using DataFrames

include("S2NonLife_Input.jl")

#################################################################
println("Start S2NonLife.jl")

lobs = Array{NLLob}(undef, 0)

for lob in select_lob_names
  push!(lobs, NLLob(df_lobs[df_lobs[!,:name] .== lob, :],
                    convert(Array, β_raw[!,lob]),
                    rf,
                    corr_prem_res, df_lobs_prem_risk))
end

## premium reserve risk #########################################
scr_prem_res = 3 * premrestotalsd(lobs, corr_lob)[1]
scr_lapse = 0.0

## catastrophe risk #############################################
## proportional reinsurance seems to be ignored for cat scr

## natural catastrophe risk -------------------------------------
scr_cat_nat = 0.0   ## no natural catastrophe risk in portfolio

## man-made catastrophe risk ------------------------------------
## only fire risk and liability risk in "man-made" portfolio

## fire
scr_cat_fire = sum(cat_fire_is) ## no risk management mechanisms

## liability
#prem_liab_vec = Array{Float64}(undef, nrow(df_cat_liability))
scr_cat_liab_vec = Array{Float64}(undef, nrow(df_cat_liability))
#cat_liab_grp = Array{Int}(undef, 0)
for 𝑙𝑜𝑏 ∈ lobs
  if 𝑙𝑜𝑏.name == :liability
    prem_liab_vec = 𝑙𝑜𝑏.prem_gross_cy * cat_liability_mix
    global scr_cat_liab_vec =
      convert(Array,
              prem_liab_vec .* df_cat_liability[:,:f_liab])
#    for 𝑖 ∈ 1:length(prem_liab_vec)
#      if prem_liab_vec[𝑖] > 0
#        push!(cat_liab_grp, 𝑖)
#      end
#    end
  end
end
scr_cat_liab =
  sqrt(scr_cat_liab_vec' * corr_liab * scr_cat_liab_vec)[1,1]

## total man-made catastrophe risk ------------------------------
scr_cat_man_made_vec = [scr_cat_fire, scr_cat_liab]
scr_cat_man_made =
  sqrt(scr_cat_man_made_vec ⋅ scr_cat_man_made_vec)

## other catastrophe risk ---------------------------------------
scr_cat_other = 0.0

## catastrophe risk for non-prop reinsurance --------------------
scr_cat_nprop = 0.0  ## XYZ does not act as reinsurer

## total catastrophe risk ---------------------------------------
scr_cat = sqrt((scr_cat_nat + scr_cat_nprop)^2 +
                 scr_cat_man_made^2 +
                 scr_cat_other^2)

## Total SCR ####################################################
scr_nl_vec = [scr_prem_res, scr_lapse, scr_cat]
scr_nl = sqrt(scr_nl_vec ⋅ (corr_scr * scr_nl_vec))

## main results #################################################
println("scr_prem_risk        :  $(round(scr_prem_res,
                                         digits = 2))")
println("scr_lapse            :  $(round(scr_lapse,
                                         digits = 2))")
println("scr_cat_man_made_fire:  $(round(scr_cat_fire,
                                         digits = 2))")
println("scr_cat_man_made_liab:  $(round(scr_cat_liab,
                                         digits = 2))")
println("scr_cat_man_made     :  $(round(scr_cat_man_made,
                                         digits = 2))")
println("scr_cat_other        :  $(round(scr_cat_other,
                                         digits = 2))")
println("scr_cat_nprob        :  $(round(scr_cat_nprop,
                                         digits = 2))")
println("scr_cat              :  $(round(scr_cat, digits = 2))")
println("scr_total            :  $(round(scr_nl, digits = 2))")

println("End S2NonLife.jl")


text:
# circles and ellipses

#Luxor
using Luxor

tiles = Tiler(500, 300, 5, 5)
width = 20
height = 25
for (pos, n) in tiles
    global width, height
    randomhue()
    ellipse(pos, width, height, :fill)
    sethue("black")
    label = string(round(width/height, digits=2))
    textcentered(label, pos.x, pos.y + 25)
    width += 2
end

#overlapping

f1 = Point(-100, 0)
f2 = Point(100, 0)
ellipsepoly = ellipse(f1, f2, 170, :none, vertices=true)
[ begin
    setgray(rescale(c, 150, 1, 0, 1))
    poly(offsetpoly(ellipsepoly, c), close=true, :fill);
    rotate(π/20)
  end
     for c in 150:-10:1 ]

# circlepath

setline(4)
tiles = Tiler(600, 250, 1, 5)
for (pos, n) in tiles
    randomhue()
    circlepath(pos, tiles.tilewidth/2, :path)
    newsubpath()
    circlepath(pos, rand(5:tiles.tilewidth/2 - 1), :fill, reversepath=true)
end

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
"""
    KLocalOperatorSum

A KLocalOperator representing the sum of several KLocalOperator-s. Internally,
the sum is stored as a vector of local operators acting on some sites.
"""
struct KLocalOperatorSum{VS<:AbstractVector,VOp} <: AbsLinearOperator
    # list of sites in this sum
    sites::VS

    # List of operators
    operators::VOp
end

KLocalOperatorSum(op::KLocalOperator) = KLocalOperatorSum([sites(op)], [op])

## Accessors
operators(op::KLocalOperatorSum) = op.operators
conn_type(op::KLocalOperatorSum) = conn_type(eltype(operators(op)))

# Copy
function duplicate(op::KLocalOperatorSum)
    new_sites = deepcopy(op.sites)
    ops = eltype(op.operators)[];
    for _op=operators(op)
        push!(ops, duplicate(_op))
    end

    KLocalOperatorSum(new_sites, ops)
end

#
function row_valdiff!(conn::OpConnection, op::KLocalOperatorSum, v::State)
    for _op=operators(op)
        row_valdiff!(conn, _op, v)
    end
    conn
end

function row_valdiff_index!(conn::OpConnectionIndex, op::KLocalOperatorSum, v::State)
    for _op=operators(op)
        row_valdiff_index!(conn, _op, v)
    end
    conn
end

function Base.sum!(op_sum::KLocalOperatorSum, op::KLocalOperator)
    id = findfirst(isequal(sites(op)), op_sum.sites)

    if isnothing(id)
        push!(op_sum.sites, sites(op))
        push!(op_sum.operators, op)
    else
        sum_samesite!(op_sum.operators[id], op)
    end

    return op_sum
end

function Base.sum!(op_l::KLocalOperatorSum, op_r::KLocalOperatorSum)
    for op=operators(op_r)
        sum!(op_l, op)
    end
    op_l
end

+(op_l::KLocalOperatorSum, op::KLocalOperator) = sum!(duplicate(op_l), op)
+(op::KLocalOperator, ops::KLocalOperatorSum) = ops + op

function Base.transpose(ops::KLocalOperatorSum)
    new_sites  = similar(ops.sites)
    new_sites .= ops.sites

    new_ops    = eltype(operators(ops))[]
    for op = operators(ops)
        push!(new_ops, transpose(op))
    end
    return KLocalOperatorSum(new_sites, new_ops)
end

function Base.conj(ops::KLocalOperatorSum)
    new_sites  = similar(ops.sites)
    new_sites .= ops.sites

    new_ops    = eltype(operators(ops))[]
    for op = operators(ops)
        push!(new_ops, conj(op))
    end
    return KLocalOperatorSum(new_sites, new_ops)
end

function Base.conj!(ops::KLocalOperatorSum)
    for op = operators(ops)
        conj!(op)
    end
    return ops
end

Base.adjoint(ops::KLocalOperatorSum) = conj!(transpose(ops))

Base.show(io::IO, ::MIME"text/plain", op::KLocalOperatorSum) = print(io,
    "KLocalOperatorSum: \n\t -sites: $(op.sites)")

Base.eltype(::T) where {T<:KLocalOperatorSum} = eltype(T)
Base.eltype(T::Type{KLocalOperatorSum{Vec,VOp}}) where {Vec,VOp} =
    eltype(eltype(VOp))


text:
srand(1)

function f(x::Vector)
    (x[1] - 5.0)^4
end

results = Optim.simulated_annealing(f, [0.0])
@assert norm(results.minimum - [5.0]) < 0.1

function rosenbrock(x::Vector)
    (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

results = Optim.simulated_annealing(rosenbrock, [0.0, 0.0])
@assert norm(results.minimum - [1.0, 1.0]) < 0.1