text:
function optimize!(model::Model)
  @assert is_valid(model.env)
  stat = (if model.has_int
    @cpx_ccall_intercept(model, mipopt, Cint, (Ptr{Void}, Ptr{Void}), model.env.ptr, model.lp)
  elseif model.has_qc
    @cpx_ccall_intercept(model, qpopt, Cint, (Ptr{Void}, Ptr{Void}), model.env.ptr, model.lp)
  else
    @cpx_ccall_intercept(model, lpopt, Cint, (Ptr{Void}, Ptr{Void}), model.env.ptr, model.lp)
    end)
  if stat != 0
    throw(CplexError(model.env, stat))
  end
end

set_branching_priority(model::Model, priority) =
    set_branching_priority(model, Cint[1:num_var(model)], priority)

set_branching_priority(model::Model, indices, priority) =
    set_branching_priority(model, indices, priority, C_NULL)

set_branching_priority(model, indices, priority, direction) =
    set_branching_priority(model, convert(Vector{Cint},indices), convert(Vector{Cint},priority), direction)

function set_branching_priority(model::Model, indices::Vector{Cint}, priority::Vector{Cint}, direction)
    @assert (cnt = length(indices)) == length(priority)
    isa(direction,Vector) && @assert cnt == length(direction)
    stat = @cpx_ccall(copyorder, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Cint,
                      Ptr{Cint},
                      Ptr{Cint},
                      Ptr{Cint}),
                      model.env.ptr, model.lp, cnt, indices.-1, priority, direction)
    stat == 0 || throw(CplexError(model.env, stat))
    return nothing
end

function get_objval(model::Model)
  objval = Array(Cdouble, 1)
  stat = @cpx_ccall(getobjval, Cint, (
                    Ptr{Void},
                    Ptr{Void},
                    Ptr{Cdouble}
                    ),
                    model.env.ptr, model.lp, objval)
  if stat != 0
    throw(CplexError(model.env, stat))
  end
  return objval[1]
end

function get_solution(model::Model)
  nvars = num_var(model)
  x = Array(Cdouble, nvars)
  stat = @cpx_ccall(getx, Cint, (
                    Ptr{Void},
                    Ptr{Void},
                    Ptr{Cdouble},
                    Cint,
                    Cint
                    ),
                    model.env.ptr, model.lp, x, 0, nvars-1)
  if stat != 0
    throw(CplexError(model.env, stat))
  end
  return x
end

function get_reduced_costs(model::Model)
    nvars = num_var(model)
    p = Array(Cdouble, nvars)
    status = Array(Cint, 1)
    stat = @cpx_ccall(getdj, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Ptr{Cdouble},
                      Cint,
                      Cint
                      ),
                      model.env.ptr, model.lp, p, 0, nvars-1)
    if stat != 0
       throw(CplexError(model.env, stat))
   end
   return p
end

function get_constr_duals(model::Model)
    ncons = num_constr(model)
    p = Array(Cdouble, ncons)
    status = Array(Cint, 1)
    stat = @cpx_ccall(getpi, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Ptr{Cdouble},
                      Cint,
                      Cint
                      ),
                      model.env.ptr, model.lp, p, 0, ncons-1)
    if stat != 0
       throw(CplexError(model.env, stat))
   end
   return p
end

function get_constr_solution(model::Model)
  ncons = num_constr(model)
  Ax = Array(Cdouble, ncons)
  stat = @cpx_ccall(getax, Cint, (
                    Ptr{Void},
                    Ptr{Void},
                    Ptr{Cdouble},
                    Cint,
                    Cint
                    ),
                    model.env.ptr, model.lp, Ax, 0, ncons-1)
  if stat != 0
    throw(CplexError(model.env, stat))
  end
  return Ax
end

function get_infeasibility_ray(model::Model)
  ncons = num_constr(model)
  y = Array(Cdouble, ncons)
  proof_p = Array(Cdouble, 1)
  stat = @cpx_ccall(dualfarkas, Cint, (
                    Ptr{Void},
                    Ptr{Void},
                    Ptr{Cdouble},
                    Ptr{Cdouble}
                    ),
                    model.env.ptr, model.lp, y, proof_p)
  if stat != 0
    error("CPLEX is unable to grab infeasible ray; consider resolving with presolve turned off")
    throw(CplexError(model.env, stat))
  end
  return y
end

function get_unbounded_ray(model::Model)
  solve_stat = get_status(model)
  if solve_stat == :CPX_STAT_UNBOUNDED
    n = num_var(model)
    z = Array(Cdouble, n)
    stat = @cpx_ccall(getray, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Ptr{Cdouble}
                      ),
                      model.env.ptr, model.lp, z)
    if stat != 0
      throw(CplexError(model.env, stat))
    end
    return z
  else
    error("CPLEX is unable to grab unbounded ray; consider resolving with presolve turned off")
  end
end

const varmap = Dict(
    0 => :NonbasicAtLower,
    1 => :Basic,
    2 => :NonbasicAtUpper,
    3 => :Free
)

const conmap = Dict(
    0 => :NonbasicAtLower,
    1 => :Basic,
    2 => :NonbasicAtUpper
)

function get_basis(model::Model)
    cval = Array(Cint, num_var(model))
    rval = Array(Cint, num_constr(model))
    stat = @cpx_ccall(getbase, Cint, (Ptr{Void},Ptr{Void},Ptr{Cint},Ptr{Cint}),
                      model.env.ptr, model.lp, cval, rval)
    stat != 0 && throw(CplexError(model.env, stat))

    csense = get_constr_senses(model)
    cbasis = Array(Symbol, num_var(model))
    rbasis = Array(Symbol, num_constr(model))
    for it in 1:num_var(model)
        cbasis[it] = varmap[cval[it]]
    end
    for it in 1:num_constr(model)
        rbasis[it] = conmap[rval[it]]
        if (rbasis[it] == :NonbasicAtLower) && (csense[it] == convert(Cchar,'L'))
            rbasis[it] = :NonbasicAtUpper
        end
    end
    return cbasis, rbasis
end

get_node_count(model::Model) = @cpx_ccall(getnodecnt, Cint, (Ptr{Void},Ptr{Void}), model.env.ptr, model.lp)

function get_rel_gap(model::Model)
  ret = Array(Cdouble,1)
  stat = @cpx_ccall(getmiprelgap, Cint, (Ptr{Void},Ptr{Void},Ptr{Cdouble}), model.env.ptr, model.lp, ret)
  if stat != 0
    throw(CplexError(model.env, stat))
  end
  ret[1]
end

function get_num_cuts(model::Model,cuttype)
    cutcount = Array(Cint,1)

    stat = @cpx_ccall(getnumcuts, Cint, (Ptr{Void},Ptr{Void},Cint,Ptr{Void}), model.env.ptr , model.lp, cuttype, cutcount)
    if stat != 0
        error(CplexError(model.inner.env, stat).msg)
    end
    return cutcount[1]
end

const status_symbols = Dict(
    1   => :CPX_STAT_OPTIMAL,
    2   => :CPX_STAT_UNBOUNDED,
    3   => :CPX_STAT_INFEASIBLE,
    4   => :CPX_STAT_INForUNBD,
    5   => :CPX_STAT_OPTIMAL_INFEAS,
    6   => :CPX_STAT_NUM_BEST,
    7   => :CPX_STAT_FEASIBLE_RELAXED,
    8   => :CPX_STAT_OPTIMAL_RELAXED,
    10  => :CPX_STAT_ABORT_IT_LIM,
    11  => :CPX_STAT_ABORT_TIME_LIM,
    12  => :CPX_STAT_ABORT_OBJ_LIM,
    13  => :CPX_STAT_ABORT_USER,
    20  => :CPX_STAT_OPTIMAL_FACE_UNBOUNDED,
    21  => :CPX_STAT_ABORT_PRIM_OBJ_LIM,
    22  => :CPX_STAT_ABORT_DUAL_OBJ_LIM,
    101 => :CPXMIP_OPTIMAL,
    102 => :CPXMIP_OPTIMAL_TOL,
    103 => :CPXMIP_INFEASIBLE,
    104 => :CPXMIP_SOL_LIM,
    105 => :CPXMIP_NODE_LIM_FEAS,
    106 => :CPXMIP_NODE_LIM_INFEAS,
    107 => :CPXMIP_TIME_LIM_FEAS,
    108 => :CPXMIP_TIME_LIM_INFEAS,
    109 => :CPXMIP_FAIL_FEAS,
    110 => :CPXMIP_FAIL_INFEAS,
    111 => :CPXMIP_MEM_LIM_FEAS,
    112 => :CPXMIP_MEM_LIM_INFEAS,
    113 => :CPXMIP_ABORT_FEAS,
    114 => :CPXMIP_ABORT_INFEAS,
    115 => :CPXMIP_OPTIMAL_INFEAS,
    116 => :CPXMIP_FAIL_FEAS_NO_TREE,
    117 => :CPXMIP_FAIL_INFEAS_NO_TREE,
    118 => :CPXMIP_UNBOUNDED,
    119 => :CPXMIP_INForUNBD,
    120 => :CPXMIP_FEASIBLE_RELAXED,
    121 => :CPXMIP_OPTIMAL_RELAXED
)

get_status(model::Model) = status_symbols[Int(get_status_code(model))]::Symbol
get_status_code(model::Model) = @cpx_ccall(getstat, Cint, (Ptr{Void}, Ptr{Void}), model.env.ptr, model.lp)


text:
## Objects that are not AbstractArrays can also be displayed with `imshow`.
# This file demonstrates and tests the interface that such objects must satisfy.

module ConeModule

using FixedPointNumbers, Colors, ImageCore, Cairo, Reactive, ImageView
using Base: tail

export Cone

# A cone in 3d
struct Cone{C<:Colorant}
    sz3::Tuple{Int,Int,Int}
    center2::Tuple{Int,Int}
    colort::Vector{C}

    function Cone{C}(sz3, colort) where C
        center2 = ((sz3[1]+1)÷2, (sz3[2]+1)÷2)
        new{C}(sz3, center2, colort)
    end
end
Cone(sz3::Tuple{Int,Int,Int}, colort::Vector{C}) where {C<:Colorant} = Cone{C}(sz3, colort)

## These define the interface that an object needs to support in order
## to be displayable with imshow

Base.eltype(::Type{Cone{C}}) where {C} = C

# What you actually need to implement is `axes(c)`, but we can do that here via `size`
Base.size(c::Cone) = (c.sz3..., length(c.colort))

# c[x, y, z, t]
function Base.getindex(c::Cone, y::Real, x::Real, z::Real, t::Real)
    ((x-c.center2[2])^2 + (y-c.center2[1])^2 <= (c.sz3[3] - z + 1)^2)*c.colort[t]
end
function Base.getindex(c::Cone, y, x, z, t)
    # indexing with vectors or colons
    inds = map(normalize_inds, axes(c), (y, x, z, t))
    sz = map(length, filtervec(inds, inds))
    yn, xn, zn, tn = inds
    ps = filtervec((1.0, 1.0, 5.0), (yn, xn, zn))  # to get the pixelspacing correct
    out = Array{eltype(c)}(undef, sz)
    k = 0
    for ti in tn, zi in zn, xi in xn, yi in yn
        out[k+=1] = c[yi,xi,zi,ti]
    end
    SpacedArray(out, ps)
end

## Helper functions
# You don't need to implement these specific functions: they're here
# only to support the specific implementations above.

# An array with custom pixelspacing. If you don't need to customize
# ps, then you could just use a plain Array in the getindex method
# above.
struct SpacedArray{T,N} <: AbstractArray{T,N}
    a::Array{T,N}
    ps::NTuple{N,Float64}
end
ImageCore.pixelspacing(A::SpacedArray) = A.ps
@inline Base.getindex(A::SpacedArray, inds...) = A.a[inds...]
@inline Base.setindex!(A::SpacedArray, val, inds...) = A.a[inds...] = val
Base.size(A::SpacedArray) = size(A.a)

# Replace colons with explicit range vectors
normalize_inds(ref, ind) = ind
normalize_inds(ref, ::Colon) = ref

# Choose elements from the first tuple based on vector entries in the second
@inline filtervec(t1, t2::Tuple{AbstractVector,Vararg{Any}}) = (t1[1], filtervec(tail(t1), tail(t2))...)
@inline filtervec(t1, t2::Tuple{Real,Vararg{Any}}) = filtervec(tail(t1), tail(t2))
filtervec(::Tuple{}, ::Tuple{}) = ()

end

using Colors, FixedPointNumbers, GtkReactive, ImageView

## Create an object and visualize it
c = ConeModule.Cone((201, 301, 31), rand(RGB{N0f8}, 60))
imshow_now(c; name="Cone 1,2")

# Slice along axes 1 and t rather than 3 and t
imshow_now(c, axes=(2,3), name="Cone 2,3")

ext: jl


text:
using LightGraphs
using Serialization
using Statistics
import Memoize
import Base.Threads.@threads
import SharedArrays.SharedArray

get_average_nn_size(graph_data::GGraphData) = length.(graph_data.neutral_networks) |> mean
function get_nn_size_percentile(graph_data::GGraphData)
    neutral_networks = graph_data.neutral_networks
    num_neutral_networks = length(neutral_networks)
    num_genotypes = sum(length.(neutral_networks))

    percentile = 0.1*num_neutral_networks

    lengths_sorted = sort(length.(neutral_networks), rev=true)
    lengths_chosen = 0
    for (i, l) in enumerate(lengths_sorted)
        lengths_chosen += l

        if i >= percentile
            return lengths_chosen / num_genotypes
        end
    end

    return 1.0
end
get_average_diameter(graph_data::GGraphData) = analyse_neutral_networks(graph_data, diameter)
get_average_radius(graph_data::GGraphData) = analyse_neutral_networks(graph_data, radius)
get_average_clustering(graph_data::GGraphData) = analyse_neutral_networks(graph_data, global_clustering_coefficient)
get_average_phenotype_robustness(graph_data::GGraphData) = mean([get_phenotype_robustness(graph_data, nn) for nn in 1:length(graph_data.neutral_networks)])
get_average_phenotype_evolvability(graph_data::GGraphData) = mean([get_phenotype_evolvability(graph_data, nn) for nn in 1:length(graph_data.neutral_networks)])
function get_average_phenotype_evolvability_robustness_cor(graph_data::GGraphData)
    cor(
        [get_phenotype_evolvability(graph_data, i) for i in 1:length(graph_data.neutral_networks)],
        [get_phenotype_robustness(graph_data, i) for i in 1:length(graph_data.neutral_networks)]
        )
end
function get_shape_space_covering(graph_data::GGraphData, n)
    unweighted_phenotype_graph = get_simple_digraph(graph_data.phenotype_graph)
    genotypes = [g for neutral_network in graph_data.neutral_networks for g in neutral_network]

    shape_space_covering = estimate(genotypes, 0.01, 100, 10_000) do genotype_index
        return get_shape_space_covering(graph_data, unweighted_phenotype_graph, genotype_index, n)
    end

    return shape_space_covering
end

get_phenotype_robustness(graph_data::GGraphData, network_index) = mean([get_genotype_robustness(graph_data, i, network_index) for i in graph_data.neutral_networks[network_index]])

function get_phenotype_evolvability(graph_data::GGraphData, network_index)
    reachable_genotypes = []

    for genotype_index in graph_data.neutral_networks[network_index]
        current_reachable_genotypes = outneighbors(graph_data.phenotype_graph, genotype_index)
        append!(reachable_genotypes, current_reachable_genotypes)
    end

    reachable_genotypes = unique(reachable_genotypes)
    reachable_networks = get_corresponding_neutral_networks(graph_data, reachable_genotypes)
    filter!(x -> x != 0, reachable_networks)                # filter out genotypes not on any neutral network
    filter!(x -> x != network_index, reachable_networks)    # filter out genotypes on the network in question
    reachable_networks = unique(reachable_networks)

    return length(reachable_networks) / (length(graph_data.neutral_networks) - 1)
end

get_genotype_robustness(graph_data, genotype_index) = get_genotype_robustness(GGraphData, genotype_index, get_corresponding_neutral_network(graph_data, genotype_index))
function get_genotype_robustness(graph_data::GGraphData, genotype_index, network_index)
    neutral_network = graph_data.neutral_networks[network_index]

    neighbors = outneighbors(graph_data.phenotype_graph, genotype_index)
    neutral_neighbors = [n for n in neighbors if n in neutral_network]

    return length(neutral_neighbors) / length(neighbors)
end

function get_shape_space_covering(graph_data::GGraphData, genotype_index, n)
    get_shape_space_covering(graph_data, get_simple_digraph(graph_data.phenotype_graph), genotype_index, n)
end
function get_shape_space_covering(graph_data::GGraphData, unweighted_phenotype_graph::SimpleDiGraph, genotype_index, n)
    reachable_genotypes = neighborhood(unweighted_phenotype_graph, genotype_index, n)

    reachable_networks = get_corresponding_neutral_networks(graph_data, reachable_genotypes) |> unique
    filter!(x -> x != 0, reachable_networks)    # filter out genotypes not on any neutral network

    return (length(reachable_networks) - 1) / (length(graph_data.neutral_networks) - 1)   # ignore neutral network of genotype_index in question
end

function analyse_neutral_networks(graph_data::GGraphData, metric)
    results = []

    for neutral_network in graph_data.neutral_networks
        graph, vertex_mapping = induced_subgraph(graph_data.genotype_graph, neutral_network)

        metric_induced_graph = build_metric_induced_graph(graph) do v_1, v_2
            genotype_index_1 = vertex_mapping[v_1]
            genotype_index_2 = vertex_mapping[v_2]

            genotype_1 = graph_data.genotype_vertex_mapping(genotype_index_1)
            genotype_2 = graph_data.genotype_vertex_mapping(genotype_index_2)

            return Levenshtein()(genotype_1, genotype_2)
        end

        current_result = metric(metric_induced_graph)

        push!(results, current_result)
    end

    return mean(results)
end

function get_corresponding_neutral_networks(graph_data::GGraphData, genotype_indices)
    [get_corresponding_neutral_network(graph_data, genotype) for genotype in genotype_indices]
end
@memoize function get_corresponding_neutral_network(graph_data::GGraphData, genotype_index)
    for (i, neutral_network) in enumerate(graph_data.neutral_networks)
        if genotype_index in neutral_network
            return i
        end
    end

    return 0
end

function get_neutral_network_distribution(graph_data::GGraphData, snapshot)
    organisms = get_organisms(snapshot)
    genotypes = [graph_data.genotype_vertex_mapping[get_genotype(snapshot, organism)] for organism in organisms]
    neutral_networks = get_corresponding_neutral_networks(graph_data, genotypes)
    nn_distribution = get_distribution(neutral_networks)

    return nn_distribution
end

function pmean(get_value_from_item, items)
    values = SharedArray{Float64}(length(items))

    for (i, item) in enumerate(items) |> unique
        value = get_value_from_item(item)
        values[i] = value
    end

    return mean(values)
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
@testset "Polar stereographic" begin

    tm = Geodesy.TransverseMercator(wgs84)

    # Data generated independely from Charles Karney's http://geographiclib.sourceforge.net/cgi-bin/GeoConvert?input=89+1&zone=0&prec=9&option=Submit

    lat = 89.0
    lon = 1.0
    x = 2001937.679954440 - 2e6
    y = 1888990.389757469 - 2e6
    northpole = true

    (x2, y2, γ, k) = Geodesy.polarst_fwd(northpole, 0.994, tm, lat, lon)
    @test (x2, y2) ≈ (x, y)

    (lat2, lon2, γ, k) = Geodesy.polarst_inv(northpole, 0.994, tm, x, y)
    @test (lat2, lon2) ≈ (lat, lon)

    lat = 88.0
    lon = 91.0
    x = 2222035.446703544 - 2e6
    y = 2003875.643138576 - 2e6
    northpole = true

    (x2, y2, γ, k) = Geodesy.polarst_fwd(northpole, 0.994, tm, lat, lon)
    @test (x2, y2) ≈ (x, y)

    (lat2, lon2, γ, k) = Geodesy.polarst_inv(northpole, 0.994, tm, x, y)
    @test (lat2, lon2) ≈ (lat, lon)

    lat = 87.0
    lon = -179.0
    x = 1994185.827037653 - 2e6
    y = 2333093.745927457 - 2e6
    northpole = true

    (x2, y2, γ, k) = Geodesy.polarst_fwd(northpole, 0.994, tm, lat, lon)
    @test (x2, y2) ≈ (x, y)

    (lat2, lon2, γ, k) = Geodesy.polarst_inv(northpole, 0.994, tm, x, y)
    @test (lat2, lon2) ≈ (lat, lon)

    lat = 86.0
    lon = -89.0
    x = 1555799.234876368 - 2e6
    y = 1992246.446803603 - 2e6
    northpole = true

    (x2, y2, γ, k) = Geodesy.polarst_fwd(northpole, 0.994, tm, lat, lon)
    @test (x2, y2) ≈ (x, y)

    (lat2, lon2, γ, k) = Geodesy.polarst_inv(northpole, 0.994, tm, x, y)
    @test (lat2, lon2) ≈ (lat, lon)

    lat = -89.0
    lon = 1.0
    x = 2001937.679954440 - 2e6
    y = 2111009.610242531 - 2e6
    northpole = false

    (x2, y2, γ, k) = Geodesy.polarst_fwd(northpole, 0.994, tm, lat, lon)
    @test (x2, y2) ≈ (x, y)

    (lat2, lon2, γ, k) = Geodesy.polarst_inv(northpole, 0.994, tm, x, y)
    @test (lat2, lon2) ≈ (lat, lon)

    lat = -88.0
    lon = 91.0
    x = 2222035.446703544 - 2e6
    y = 1996124.356861424 - 2e6
    northpole = false

    (x2, y2, γ, k) = Geodesy.polarst_fwd(northpole, 0.994, tm, lat, lon)
    @test (x2, y2) ≈ (x, y)

    (lat2, lon2, γ, k) = Geodesy.polarst_inv(northpole, 0.994, tm, x, y)
    @test (lat2, lon2) ≈ (lat, lon)

    lat = -87.0
    lon = -179.0
    x = 1994185.827037653 - 2e6
    y = 1666906.254072543 - 2e6
    northpole = false

    (x2, y2, γ, k) = Geodesy.polarst_fwd(northpole, 0.994, tm, lat, lon)
    @test (x2, y2) ≈ (x, y)

    (lat2, lon2, γ, k) = Geodesy.polarst_inv(northpole, 0.994, tm, x, y)
    @test (lat2, lon2) ≈ (lat, lon)

    lat = -86.0
    lon = -89.0
    x = 1555799.234876368 - 2e6
    y = 2007753.553196397 - 2e6
    northpole = false

    (x2, y2, γ, k) = Geodesy.polarst_fwd(northpole, 0.994, tm, lat, lon)
    @test (x2, y2) ≈ (x, y)

    (lat2, lon2, γ, k) = Geodesy.polarst_inv(northpole, 0.994, tm, x, y)
    @test (lat2, lon2) ≈ (lat, lon)
end


text:
#################### Model Expression Operators ####################

# function modelfx(literalargs::Vector{Tuple{Symbol, DataType}}, f::Function)
#   modelfxsrc(literalargs, f)[1]
# end
#
# function modelfxsrc(literalargs::Vector{Tuple{Symbol, DataType}}, f::Function)
#   args = Expr(:tuple, map(arg -> Expr(:(::), arg[1], arg[2]), literalargs)...)
#   expr, src = modelexprsrc(f, literalargs)
#   fx = eval(Main, Expr(:function, args, expr))
#   (fx, src)
# end
#
#
# function modelexprsrc(f::Function, literalargs::Vector{Tuple{Symbol, DataType}})
#   m = first(methods(f).ms)
#   argnames = Vector{Any}(m.nargs)
#   ccall(:jl_fill_argnames, Void, (Any, Any), m.source, argnames)
#
#   fkeys = Symbol[argnames[2:end]...]
#   ftypes = DataType[m.sig.parameters[2:end]...]
#   n = length(fkeys)
#
#   literalinds = Int[]
#   for (key, T) in literalargs
#     i = findfirst(fkey -> fkey == key, fkeys)
#     if i != 0 && ftypes[i] == T
#       push!(literalinds, i)
#     end
#   end
#   nodeinds = setdiff(1:n, literalinds)
#
#   all(T -> T == Any, ftypes[nodeinds]) ||
#     throw(ArgumentError("model node arguments are not all of type Any"))
#
#   modelargs = Array{Any}(n)
#   for i in nodeinds
#     modelargs[i] = Expr(:ref, :model, QuoteNode(fkeys[i]))
#   end
#   for i in literalinds
#     modelargs[i] = fkeys[i]
#   end
#   expr = Expr(:block, Expr(:(=), :f, f), Expr(:call, :f, modelargs...))
#
#   (expr, fkeys[nodeinds])
# end

#################### Mathematical Operators ####################

isprobvec(p::AbstractVector) = isprobvec(convert(Vector{Float64}, p))

cummean(x::AbstractArray) = mapslices(cummean, x, dims = [1])

function cummean(x::AbstractVector{T}) where T<:Real
  y = similar(x, Float64)
  xs = 0.0
  for i in 1:length(x)
    xs += x[i]
    y[i] = xs / i
  end
  y
end

# dot(x) = dot(x, x)

# logit(x::Real) = log(x / (1.0 - x))
# invlogit(x::Real) = 1.0 / (exp(-x) + 1.0)

## Csorgo S and Faraway JJ. The exact and asymptotic distributions of the
## Cramer-von Mises statistic. Journal of the Royal Statistical Society,
## Series B, 58: 221-234, 1996.
function pcramer(q::Real)
  p = 0.0
  for k in 0:3
    c1 = 4.0 * k + 1.0
    c2 = c1^2 / (16.0 * q)
    p += gamma(k + 0.5) / factorial(k) * sqrt(c1) * exp(-c2) * besselk(0.25, c2)
  end
  p / (pi^1.5 * sqrt(q))
end

#################### Auxiliary Functions ####################

# ## pmap2 is a partial work-around for the pmap issue in julia 0.4.0 of worker
# ## node errors being blocked.  In single-processor mode, pmap2 calls map
# ## instead to avoid the error handling issue.  In multi-processor mode, pmap is
# ## called and will apply its error processing.
#
# function pmap2(f::Function, lsts::AbstractArray)
#   if (nprocs() > 1)
#     @everywhere importall Mamba
#     pmap(f, lsts)
#   else
#     map(f, lsts)
#   end
# end


text:
using SparseArrays, NamedArrays

function _make_incidence_matrix(applications::Vector{Application}, levelfun)
    class = reduce(vcat, classification.(applications)) |> unique
    class = levelfun.(class) |> unique
    classdict = Dict(c => i for (i, c) in enumerate(class))
    incmat = falses(length(applications), length(class))
    for (i, a) in enumerate(applications)
        for c in classification(a)
            incmat[i, classdict[levelfun(c)]] = true
        end
    end
    return incmat, class
end

function coclassification(applications::Vector{Application}; level=maingroup)
    incmat, class = _make_incidence_matrix(applications, level)
    incmat = sparse(incmat)
    coclass = incmat' * incmat
    coclass = NamedArray(Matrix(coclass))
    setnames!(coclass, string.(class), 1)
    setnames!(coclass, string.(class), 2)
    return coclass
end


text:
module Variable

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

# Definition of a variable bridge
include("bridge.jl")

# Mapping between variable indices and bridges
include("map.jl")

# Bridge optimizer bridging a specific variable bridge
include("single_bridge_optimizer.jl")

# TODO(odow): the compiler in Julia <= 1.2 (and in later versions unless
# fixed) gets stuck compiling add_constrained_variable for some inputs. This
# method seemed necessary to fix it.
# See https://github.com/JuliaLang/julia/issues/32167 for more.
function MOI.Bridges.Variable.bridge_constrained_variable(BridgeType, b, s)
    return throw(
        MOI.UnsupportedConstraint{
            MOIU.variable_function_type(typeof(s)),
            typeof(s),
        }(),
    )
end

# Variable bridges
include("zeros.jl")
const Zeros{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{ZerosBridge{T},OT}
include("free.jl")
const Free{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{FreeBridge{T},OT}
include("flip_sign.jl")
const NonposToNonneg{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{NonposToNonnegBridge{T},OT}
include("vectorize.jl")
const Vectorize{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{VectorizeBridge{T},OT}
include("soc_to_rsoc.jl")
const SOCtoRSOC{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{SOCtoRSOCBridge{T},OT}
include("rsoc_to_soc.jl")
const RSOCtoSOC{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{RSOCtoSOCBridge{T},OT}
include("rsoc_to_psd.jl")
const RSOCtoPSD{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{RSOCtoPSDBridge{T},OT}

"""
    add_all_bridges(bridged_model, T::Type)

Add all bridges defined in the `Bridges.Variable` submodule to `bridged_model`.
The coefficient type used is `T`.
"""
function add_all_bridges(bridged_model, T::Type)
    MOIB.add_bridge(bridged_model, ZerosBridge{T})
    MOIB.add_bridge(bridged_model, FreeBridge{T})
    MOIB.add_bridge(bridged_model, NonposToNonnegBridge{T})
    MOIB.add_bridge(bridged_model, VectorizeBridge{T})
    MOIB.add_bridge(bridged_model, SOCtoRSOCBridge{T})
    MOIB.add_bridge(bridged_model, RSOCtoSOCBridge{T})
    MOIB.add_bridge(bridged_model, RSOCtoPSDBridge{T})
    return
end

end


text:
module SimpleANOVA

using Distributions, Requires
include("InvertedIndices.jl")
include("AnovaEffect.jl")
include("AnovaValue.jl")
include("AnovaFactor.jl")
include("AnovaResult.jl")
include("AnovaData.jl")
include("FactorType.jl")

const totalname = "Total"
const cellsname = "Cells"
const errorname = "Error"
const remaindername = "Remainder"

function __init__()
    @require DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" include("anova_dataframes.jl")
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("anova_plots.jl")
end

"""
    anova(observations::Array{Union{Number, Vector{Number}}}, factortypes = FactorType[]; factornames = String[], hasreplicates = true)
    anova(observations::Vector{Number}, factorassignments::Vector{Vector{Any}}, factortypes = FactorType[]; factornames = String[], hasreplicates = true)
    anova(df::DataFrame, observationscolumn::Symbol, factorcolumns::Vector{Symbol}, factortypes = FactorType[]; factornames = String[])

Performs an Analysis of Variance (ANOVA) calculation.

Operates on up to 3 crossed factors (fixed or random) and arbitrarily many random nested factors, with or without
replicates, on balanced data.

# Arguments
- `observations`: Array containing the values to test. For the array, each dimension is a factor level, such that observations[2,5,3] indicates the 2nd level of the first factor, the 5th level of the second factor, and the 3rd level of the third factor. May contain values or vectors of values, where the vector contains replicates. Factors should be ordered with least significant first. For the vector, must provide `factorassignments` to specify factor levels.
- `factorassignments`: Vector of vectors of integers specifying how each observation is assigned to a factor level. Provide this when `observations` is given as a vector. Factor levels do not have to be consecutive or ordered. Nested factors must reuse factor levels currently.
- `factortypes`: Vector indicating the `FactorType` for each factor. If present, `replicates` must appear first, any `nested` after, and then `random` or `fixed` in any order. Specify `replicates` if the first dimension of the `observations` matrix contains replicate values (vs. contained in vectors). If too few values are provided, remaining are assumed to be `fixed`.
- `factornames`: Vector of names for each factor, excluding the replicate factor. If empty, will be automatically populated alphabetically.

Notes: The last index will be the top factor in the table.

Output: `AnovaData` structure containing the test results for each factor.

# Examples
```julia
anova(observations)                        # N-way fixed-effects ANOVA with replicates (vectors or first dimension)
anova(observations, hasreplicates = false) # N-way fixed-effects ANOVA without replicates (first dimension)
anova(observations, [random])              # N-way ANOVA with lower random factor and 1 or 2 upper fixed factors
anova(observations, [random])              # N-way ANOVA with lower random factor and 1 or 2 upper fixed factors
anova(observations, [fixed, random])       # N-way ANOVA with 1 lower fixed factor, 1 random factor, and 0 or 1 upper fixed factor
anova(observations, [nested, random])      # N-way fixed-effects ANOVA with 1 random nested factor, 1 random factor, and 1-2 fixed factors
```

# Glossary
- observation: The dependent variable.
- factor: An independent variable.
- factor level: A value of a factor.
- balanced: All combinations of factor levels have the same number of observations.
- crossed factor: A factor with levels that combine with the levels of all other crossed factors.
- fixed factor: A factor with fixed effects (e.g. treatment, concentration, exposure time).
- random factor: A factor with random effects (e.g. location, individual).
- nested factor: A random factor where the levels are unique to a combination of crossed factor levels (e.g. replicate).
- sum of squares (SS): A measure of variance that is dependent on sample size. Also called "sum of squared deviations."
- degrees of freedom (DF, ν): The number of bins in which the values could have been moved, if random.
- mean square (MS): SS / DF. Corrects for the larger variance expected if random values can be assigned to more bins. Also called "mean squared error" or "mean squared deviation."
- F-statistic: The division of MS values produce a result belonging to the "F distribution", the shape of which depends on the DF of the numerator and denominator. The location of this value on the distribution provides the p-value.
- p-value: The probability that, if all measurements had been drawn from the same population, you would obtain data at least as extreme as contained in your observations.
"""
function anova(observations::AbstractArray{T}, factortypes::Vector{FactorType} = FactorType[]; factornames::Vector{<:AbstractString} = String[], hasreplicates = true) where {T <: Union{Number, AbstractVector{<:Number}}}
    length(observations) > 0 || return

    firstlevelreplicates = eltype(observations) <: Number ? hasreplicates : false
    nfactors = ndims(observations) - (firstlevelreplicates ? 1 : 0)

    # defaults to assuming all unspecified factors are fixed
    if length(factortypes) < nfactors
        nremaining = nfactors - length(factortypes)
        append!(factortypes, repeat([fixed], nremaining))
    end

    validate(factortypes, factornames, nfactors)

    # automatically assigns alphabetical names if not provided.
    if isempty(factornames)
        factornames = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M", "N", "O", "P", "Q", "R", "S", "T", "U", "V", "W", "X", "Y", "Z"][1:nfactors]
        reverse!(factornames)
    end

    nnestedfactors = count(f -> f == nested, factortypes)
    crossedfactortypes = filter(f -> f ∈ [fixed, random], factortypes)
    reverse!(crossedfactortypes)
    ncrossedfactors = length(crossedfactortypes)
    ncrossedfactors < 4 || error("ANOVA with 4 or more crossed factors is not supported.")
    nreplicates = firstlevelreplicates ? size(observations, 1) : length(observations[1])
    firstlevelreplicates || all(c -> length(c) == nreplicates, observations) || throw(ErrorException("All cells must have the same number of replicates."))
    ncells = Int.(length(observations) / (firstlevelreplicates ? nreplicates : 1))
    nfactorlevels = firstlevelreplicates ? [size(observations)...][Not(1)] : [size(observations)...]

    crossedfactornames = factornames[factortypes .≠ nested]
    nestedfactornames = factornames[factortypes .== nested]

    #10kb allocated before this point
    anovakernel(observations, nreplicates, ncells, nnestedfactors, ncrossedfactors, nfactorlevels, crossedfactortypes, crossedfactornames, nestedfactornames)
end

function anova(observations::AbstractVector{T}, factorassignments::AbstractVector{<:AbstractVector}, factortypes::Vector{FactorType} = FactorType[]; factornames::Vector{<:AbstractString} = String[]) where {T <: Number}
    length(observations) > 0 || return
    nfactors = length(factorassignments)
    N = length(observations)
    N % nfactors == 0 || error("Design is unbalanced.")
    all(length.(factorassignments) .== N) || error("Each observation must have an assignment for each factor.")

    factorlevels = factorassignments .|> unique .|> sort
    nfactorlevels = length.(factorlevels)
    all(N .& nfactorlevels .== 0) || error("Design is unbalanced.")
    factorlevelcounts = [[count(l -> l == factorlevels[i][j], factorassignments[i]) for j ∈ 1:nfactorlevels[i]] for i ∈ 1:nfactors]
    nperfactorlevel = factorlevelcounts .|> unique
    all(nperfactorlevel .|> length .== 1) || error("Design is unbalanced.")
    nperfactorlevel = nperfactorlevel .|> first

    if !(isa(factorassignments, Number)) || any(maximum.(factorlevels) .> nfactorlevels)
        compressedfactorlevels = [1:i for i ∈ nfactorlevels]
        factorlevelremapping = [factorlevels[i] .=> compressedfactorlevels[i] for i ∈ 1:nfactors]
        factorassignments = [replace(factorassignments[i], factorlevelremapping[i]...) for i ∈ 1:nfactors]
    end

    nreplicates = Int(N / prod(nfactorlevels))

    nlevels = [nreplicates; nfactorlevels]
    sortorder = sortperm(repeat(1:nreplicates, Int(N / nreplicates)) .+
                         sum([factorassignments[i] .* prod(nlevels[1:i]) for i ∈ 1:nfactors]))
    observationsmatrix = reshape(observations[sortorder], nlevels...)
    anova(observationsmatrix, factortypes, factornames = factornames, hasreplicates = nreplicates > 1)
end

function validate(factortypes::Vector{FactorType}, factornames::Vector{<:AbstractString}, nfactors)
    if !isempty(factortypes)
        length(factortypes) == nfactors || error("factortypes must have an entry for each factor.")
        nested ∉ factortypes || factortypes[1:count(t -> t == nested, factortypes)] |> unique |> length == 1 || error("nested entries must come before crossed factors")
    end

    if !isempty(factornames)
        nfactors ≤ 26 || error("Can only automatically name up to 26 factors. Provide names explicitly.")
        length(factornames) == nfactors || error("factornames must have an entry for each factor.")
    end
end

function anovakernel(observations, nreplicates, ncells, nnestedfactors, ncrossedfactors, nfactorlevels, crossedfactortypes, crossedfactornames, nestedfactornames)
    N = ncells * nreplicates
    nfactors = nnestedfactors + ncrossedfactors

    # collapse replicate dimension
    cellsums = eltype(observations) <: Number && nreplicates == 1 ? observations : sumfirstdim(observations)
    C = sum(cellsums) ^ 2 / N
    total = totalcalc(observations, N, C)
    amongallnested, nestedsums, ncrossedfactorlevels, nnestedfactorlevels = amongnestedfactorscalc(cellsums, nfactorlevels, nnestedfactors, nreplicates, C)
    cells = cellscalc(cellsums, nreplicates, ncells, C)

    crossedfactors = factorscalc(nestedsums, ncrossedfactors, ncrossedfactorlevels, N, C, crossedfactornames) # 3kb allocated here, possibly can't be avoided
    interactions, interactionsmap = interactionscalc(cells, nestedsums, crossedfactors, ncrossedfactors, ncrossedfactorlevels, nnestedfactorlevels, nreplicates, C, crossedfactornames) # 9 kb allocated here!
    nestedfactors = nestedfactorscalc(amongallnested, nnestedfactors, crossedfactors, interactions, nestedfactornames)
    error = errorcalc(total, amongallnested, cells, [crossedfactors; interactions[1:end-1]], nnestedfactors, nreplicates)

    reverse!(crossedfactors)

    numerators = getnumerators(crossedfactors, ncrossedfactors, nnestedfactors, nestedfactors, interactions)
    crossedbasedenominator = nnestedfactors > 0 ? nestedfactors[end] : error;
    denominators = getdenominators(nnestedfactors, nestedfactors, nreplicates, crossedbasedenominator, error, total, crossedfactors, ncrossedfactors, crossedfactortypes, interactionsmap)# 2-3kb allocated, 50 allocations

    # drop least significant test if nreplicates == 1; either the lowest interaction level, or lowest nesting level if present
    if nreplicates == 1
        droppedfactor = pop!(numerators)
        pop!(denominators)
    end

    # perform test
    results = ftest.(numerators, denominators)

    npercell = nreplicates * prod(nnestedfactorlevels)
    nestedmeans = nestedsums ./ npercell
    data = AnovaData([total; results], crossedfactors, denominators[1:ncrossedfactors], nestedmeans, npercell)
    nnestedfactors > 0 && nreplicates == 1 && push!(data.effects, droppedfactor)
    error.df > 0 && push!(data.effects, error)

    return data
end

function sumfirstdim(observations::T) where {T <: AbstractArray{<:AbstractVector}}
    map(sumfirstdim, observations)
end

function sumfirstdim(observations::T) where {T <: AbstractArray{<:Number}}
    dropdims(sum(observations, dims = 1), dims = 1)
end

function sumfirstdim(observations::T) where {T <: AbstractVector{<:Number}}
    sum(observations, dims = 1) |> first
end

function totalcalc(observations, N, C)
    ss = sum(c -> sum(c.^2), observations) - C
    df = N - 1
    AnovaValue(totalname, ss, df)
end

function factorscalc(cellsums, nfactors, nfactorlevels, N, C, factornames)
    factorindices = 1:nfactors
    ss = map(i -> sum(sum(cellsums, dims = factorindices[Not(i)]) .^ 2) / (N / nfactorlevels[i]), factorindices) .- C
    df = nfactorlevels .- 1
    AnovaFactor.(factornames, ss, df)
end

function cellscalc(cellsums, nreplicates, ncells, C)
    ss = sum(cellsums .^ 2) / nreplicates - C
    df = ncells - 1
    AnovaValue(cellsname, ss, df)
end

function errorcalc(total, amongallnested, cells, otherfactors, nnestedfactors, nreplicates)
    if nnestedfactors > 0
        otherfactor = amongallnested[1]
        name = errorname
    elseif nreplicates > 1
        otherfactor = cells
        name = errorname
    else
        otherfactor = AnovaValue("", sum(f -> f.ss, otherfactors), sum(f -> f.df, otherfactors))
        name = remaindername
    end

    ss = total.ss - otherfactor.ss
    df = total.df - otherfactor.df
    AnovaFactor(name, ss, df)
end

function amongnestedfactorscalc(cellsums, nfactorlevels, nnestedfactors, nreplicates, C)
    nestedsums = cellsums
    nlowerfactorlevels = 1
    nupperfactorlevels = nfactorlevels
    amongallnested = Vector{AnovaValue}(undef, nnestedfactors)
    for i ∈ 1:nnestedfactors
        # collapse each nested factor
        ss = sum(nestedsums .^ 2 ./ (nreplicates * prod(nlowerfactorlevels))) - C
        df = prod(nupperfactorlevels) - 1
        amongallnested[i] = AnovaValue(ss, df)
        nlowerfactorlevels = nfactorlevels[1:i]
        nupperfactorlevels = nfactorlevels[(i+1):end]
        nestedsums = sumfirstdim(nestedsums)
    end
    amongallnested, nestedsums, nupperfactorlevels, nlowerfactorlevels
end

function interactionscalc(cells, nestedsums, crossedfactors, ncrossedfactors, ncrossedfactorlevels, nnestedfactorlevels, nreplicates, C, crossedfactorlabels)
    interactionsmap = Dict{Any,AnovaFactor}()
    interactions = Vector{AnovaFactor}()
    if ncrossedfactors > 1
        pairwise = pairwisecalc(nestedsums, crossedfactors, ncrossedfactors, ncrossedfactorlevels, nnestedfactorlevels, nreplicates, C, crossedfactorlabels)
        interactionsmap[(1,2)] = pairwise[1,2]
        interactionsmap[(2,1)] = pairwise[1,2]
        push!(interactions, interactionsmap[(1,2)])
        if ncrossedfactors > 2
            threewise = threewisecalc(cells, crossedfactors, pairwise, crossedfactorlabels)
            interactionsmap[(1,3)] = pairwise[1,3]
            interactionsmap[(3,1)] = pairwise[1,3]
            interactionsmap[(2,3)] = pairwise[2,3]
            interactionsmap[(3,2)] = pairwise[2,3]
            interactionsmap[(1,2,3)] = threewise
            push!(interactions, interactionsmap[(1,3)], interactionsmap[(2,3)])
            reverse!(interactions)
            push!(interactions, interactionsmap[(1,2,3)])
        end
    end
    return interactions, interactionsmap
end

function pairwisecalc(nestedsums, crossedfactors, ncrossedfactors, ncrossedfactorlevels, nnestedfactorlevels, nreplicates, C, crossedfactorlabels)
    pairwise = Array{AnovaFactor,2}(undef, ncrossedfactors, ncrossedfactors)
    factorindexes = 1:ncrossedfactors
    for i ∈ factorindexes
        for j ∈ (i+1):ncrossedfactors
            otherfactorindexes = intersect(factorindexes[Not(i)], factorindexes[Not(j)])
            ss = sum(sum(nestedsums, dims = otherfactorindexes) .^ 2 ./ (prod(ncrossedfactorlevels[otherfactorindexes]) * prod(nnestedfactorlevels) * nreplicates)) - C - crossedfactors[i].ss - crossedfactors[j].ss
            df = crossedfactors[i].df * crossedfactors[j].df
            pairwise[i,j] = pairwise[j,i] = AnovaFactor("$(crossedfactorlabels[j]) × $(crossedfactorlabels[i])", ss, df)
        end
    end
    return pairwise
end

function threewisecalc(cells, factors, pairwise, crossedfactorlabels)
    ss = cells.ss - sum(f -> f.ss, factors) - pairwise[1,2].ss - pairwise[1,3].ss - pairwise[2,3].ss
    df = prod(f -> f.df, factors)
    AnovaFactor("$(crossedfactorlabels[3]) × $(crossedfactorlabels[2]) × $(crossedfactorlabels[1])", ss, df)
end

function nestedfactorscalc(amongallnested, nnestedfactors, crossedfactors, interactions, nestedfactorlabels)
    nestedfactors = []
    if nnestedfactors > 0
        otherfactors = [crossedfactors; interactions]
        otherss = sum(f -> f.ss, otherfactors)
        otherdf = sum(f -> f.df, otherfactors)
        for i ∈ nnestedfactors:-1:1
            ss = amongallnested[i].ss - otherss
            df = amongallnested[i].df - otherdf
            otherss += ss
            otherdf += df
            push!(nestedfactors, AnovaFactor(nestedfactorlabels[i], ss, df))
        end
    end
    return nestedfactors
end

function getnumerators(crossedfactors, ncrossedfactors, nnestedfactors, nestedfactors, interactions)
    numerators = copy(crossedfactors)

    if ncrossedfactors > 1
        push!(numerators, interactions[1])
        if ncrossedfactors > 2
            push!(numerators, interactions[2], interactions[3], interactions[4])
        end
    end

    if nnestedfactors > 0
        append!(numerators, nestedfactors)
    end

    return numerators
end

function getdenominators(nnestedfactors, nestedfactors, nreplicates, crossedbasedenominator, error, total, crossedfactors, ncrossedfactors, crossedfactortypes, interactionsmap)
    if ncrossedfactors == 1
        denominators = [crossedbasedenominator]
    elseif ncrossedfactors == 2
        if all(f -> f == fixed, crossedfactortypes)
            crosseddenominators = repeat([crossedbasedenominator], ncrossedfactors)
        elseif all(f -> f == random, crossedfactortypes)
            crosseddenominators = repeat([interactionsmap[(1,2)]], ncrossedfactors)
        else
            crosseddenominators = map(f -> f == fixed ? interactionsmap[(1,2)] : crossedbasedenominator, crossedfactortypes)
        end

        denominators = [crosseddenominators; crossedbasedenominator]
    elseif ncrossedfactors == 3
        if all(f -> f == fixed, crossedfactortypes)
            crosseddenominators = repeat([crossedbasedenominator], ncrossedfactors)
            pairwiseinteractiondenominators = repeat([crossedbasedenominator], ncrossedfactors)
        elseif all(f -> f == random, crossedfactortypes)
            crosseddenominators = Vector{AnovaFactor}(undef, ncrossedfactors)
            for i ∈ 1:ncrossedfactors
                otherfactors = (1:ncrossedfactors)[Not(i)]
                j = otherfactors[1]
                k = otherfactors[2]
                crosseddenominators[i] = threewayinteraction(interactionsmap[(i,j)], interactionsmap[(i,k)], interactionsmap[(1,2,3)])
            end
            pairwiseinteractiondenominators = repeat([interactionsmap[(1,2,3)]], ncrossedfactors)
        elseif count(f -> f == random, crossedfactortypes) == 1
            i = findfirst(f -> f == random, crossedfactortypes)

            crosseddenominators = Vector{AnovaFactor}(undef, ncrossedfactors)
            crosseddenominators[i] = crossedbasedenominator

            fixedindexes = (1:ncrossedfactors)[Not(i)]
            for j ∈ fixedindexes
                crosseddenominators[j] = interactionsmap[(i,j)]
            end

            fixedinteractionindex = sum(fixedindexes) - 2
            pairwiseinteractiondenominators = Vector{AnovaFactor}(undef, ncrossedfactors)
            pairwiseinteractiondenominators[fixedinteractionindex] = interactionsmap[(1,2,3)]
            pairwiseinteractiondenominators[Not(fixedinteractionindex)] .= crossedbasedenominator
        elseif count(f -> f == random, crossedfactortypes) == 2
            i = findfirst(f -> f == fixed, crossedfactortypes)
            otherfactors = (1:ncrossedfactors)[Not(i)]
            j = otherfactors[1]
            k = otherfactors[2]

            crosseddenominators = Vector{AnovaFactor}(undef, ncrossedfactors)
            crosseddenominators[i] = threewayinteraction(interactionsmap[(i,j)], interactionsmap[(i,k)], interactionsmap[(1,2,3)])
            crosseddenominators[otherfactors] .= interactionsmap[(j,k)]

            randominteractionindex = sum(otherfactors) - 2
            pairwiseinteractiondenominators = Vector{AnovaFactor}(undef, ncrossedfactors)
            pairwiseinteractiondenominators[randominteractionindex] = crossedbasedenominator
            pairwiseinteractiondenominators[Not(randominteractionindex)] .= interactionsmap[(1,2,3)]
        end

        denominators = [crosseddenominators; pairwiseinteractiondenominators; crossedbasedenominator]
    end

    # determine correct denominators for nested factors
    if nnestedfactors > 0
        nesteddenominators = Vector{AnovaFactor}(undef, nnestedfactors)
        nesteddenominators[1:end-1] = nestedfactors[2:end]
        nesteddenominators[end] = error
        append!(denominators, nesteddenominators)
    end

    denominators
end

function threewayinteraction(interaction_ab, interaction_bc, interaction_abc)
    reducedmeansquare(factor::AnovaFactor) = factor.ms ^ 2 / factor.df

    ms = interaction_ab.ms + interaction_bc.ms - interaction_abc.ms
    df = ms ^ 2 / (reducedmeansquare(interaction_ab) + reducedmeansquare(interaction_bc) + reducedmeansquare(interaction_abc))
    AnovaFactor("", ms * df, df, ms)
end

function ftest(x, y)
    f = x.ms / y.ms
    fdist = FDist(x.df, y.df)
    p = ccdf(fdist, f)
    AnovaResult(x, f, p)
end

export anova, ftest, plot

end


text:
using Requires

function __init__()
  @require Petri="4259d249-1051-49fa-8328-3f8ab9391c33" begin
    import .Petri

    # Interoperability with Petri.jl
    Petri.Model(p::AbstractPetriNet) = begin
      ts = TransitionMatrices(p)
      t_in = map(i->Dict(k=>v for (k,v) in enumerate(ts.input[i,:]) if v != 0), 1:nt(p))
      t_out = map(i->Dict(k=>v for (k,v) in enumerate(ts.output[i,:]) if v != 0), 1:nt(p))

      Δ = Dict(i=>t for (i,t) in enumerate(zip(t_in, t_out)))
      return Petri.Model(ns(p), Δ)
    end

    Petri.Model(p::Union{AbstractLabelledPetriNet, AbstractLabelledReactionNet}) = begin
      snames = [sname(p, s) for s in 1:ns(p)]
      tnames = [tname(p, t) for t in 1:nt(p)]
      ts = TransitionMatrices(p)
      t_in = map(i->LVector(;[(snames[k]=>v) for (k,v) in enumerate(ts.input[i,:]) if v != 0]...), 1:nt(p))
      t_out = map(i->LVector(;[(snames[k]=>v) for (k,v) in enumerate(ts.output[i,:]) if v != 0]...), 1:nt(p))

      Δ = LVector(;[(tnames[i]=>t) for (i,t) in enumerate(zip(t_in, t_out))]...)
      return Petri.Model(collect(values(snames)), Δ)
    end
  end
end


text:
immutable MPI_2D_Config <: MPI_areaConfig
  ##  Voodoo numbers controlling data layout.
  ##  sNx :: No. X points in sub-grid.
  ##  sNy :: No. Y points in sub-grid.
  ##  OLx :: Overlap extent in X.
  ##  OLy :: Overlat extent in Y.
  ##  nPx :: No. of processes to use in X.
  ##  nPy :: No. of processes to use in Y.
  sN :: Tuple{Int64,Int64}
  OL :: Tuple{Int64,Int64}
  nP :: Tuple{Int64,Int64}
  rank :: Int64
  size :: Int64
  function MPI_2D_Config(n :: Tuple{Int64,Int64}, l :: Tuple{Int64,Int64}, p :: Tuple{Int64,Int64})
    @assert n[1] > l[1]
    @assert n[2] > l[2]
    comm = MPI.COMM_WORLD
    new(n,l,p, MPI.Comm_rank(comm), MPI.Comm_size(comm))
  end
end

function MPI2D_InitialArray!(p :: MPI_2D_Config )
  return fill(0.0, )
end

function MPI2D_GatherAll!(p :: MPI_2D_Config, T :: Array{Float64,2}, Twhole :: Array{Float64,2})
  comm = MPI.COMM_WORLD
  (OLx, OLy) = p.OL
  (sNx, sNy) = p.sN
  (nPx, nPy) = p.nP
  rank = p.rank

  if (rank!=0)
    ##  Sending to Coupler
    sreq = MPI.Isend(T,0, rank+32,comm)
    MPI.Waitall!([sreq;])
  else
    Ti = fill(NaN, size(T))
    Twhole[1:sNx, 1:sNy] = T[OLx+(1:sNx),OLy+(1:sNy)]
    for ip = 1:(nPx*nPy-1)
      ix = mod(ip,nPx)
      iy = round(Int64, (ip-ix)/nPx)
      rreq = MPI.Irecv!(Ti, ip, ip+32, comm)
      MPI.Waitall!([rreq;])
      Twhole[ix*sNx+ (1:sNx),iy*sNy+(1:sNy)] = Ti[OLx+(1:sNx),OLy + (1:sNy)]
    end
  end
  return Twhole
end

function MPI2D_DistributeCore!(p :: MPI_2D_Config, T :: Array{Float64,2},  Twhole :: Array{Float64,2})
  comm = MPI.COMM_WORLD
  (OLx, OLy) = p.OL
  (sNx, sNy) = p.sN
  (nPx, nPy) = p.nP
  rank = p.rank

  if (rank != 0)
    ##  Receiving From Coupler
    Ti = fill(NaN, p.sN)
    csreq = MPI.Irecv!(Ti,0, rank+32,comm)
    MPI.Waitall!([csreq;])
    T[OLx+(1:sNx), OLy+(1:sNy)] = Ti
  else
    Ti = fill(NaN, size(T))
    T[(1+OLx):(OLx+sNx),(1+OLy):(OLy+sNy)] =  Twhole[1:sNx, 1:sNy]
    ## Send to Workers
    for ip = 1:(nPx*nPy-1)
      ix = mod(ip,nPx)
      iy = round(Int64, (ip-ix)/nPx)
      Ti = fill(NaN, p.sN)
      Ti[:,:] = Twhole[ix*sNx+ (1:sNx),iy*sNy+(1:sNy)]
      sreq = MPI.Isend(Ti, ip, ip+32, comm)
      MPI.Waitall!([sreq;])
    end
  end
  return Twhole
end

function MPI2D_DistributeBoundary!(p :: MPI_2D_Config, T :: Array{Float64,2},  Twhole :: Array{Float64,2})
  comm = MPI.COMM_WORLD
  #= MPI.Barrier(comm) =#
  (OLx, OLy) = p.OL
  (sNx, sNy) = p.sN
  (nPx, nPy) = p.nP
  rank = p.rank

  if (rank!=0)

    crreq = Array(MPI.Request,4)
    Twest = fill(NaN, (OLx, sNy))
    Tnorth = fill(NaN,(sNx, OLy))
    Teast = fill(NaN, (OLx, sNy))
    Tsouth = fill(NaN, (sNx, OLy))

    ##  Receiving from Coupler
    crreq[1] = MPI.Irecv!(Twest,0, 4*rank+1, comm)
    crreq[2] = MPI.Irecv!(Tnorth,0, 4*rank+2, comm)
    crreq[3] = MPI.Irecv!(Teast,0, 4*rank+3, comm)
    crreq[4] = MPI.Irecv!(Tsouth,0, 4*rank+4, comm)
    MPI.Waitall!(crreq)
    T[1:OLx,OLy+(1:sNy)] = Twest
    T[OLx+(1:sNx),sNy+OLy+(1:OLy)] = Tnorth
    T[sNx+OLx+(1:OLx),OLy+(1:sNy)] = Teast
    T[OLx+(1:sNx),1:OLy] = Tsouth
  else

    Twest = Twhole[sNx*nPx-OLx+(1:OLx),1:sNy]
    Tnorth = Twhole[1:sNx, mod(sNy+(1:OLy)-1, sNy*nPy)+1]
    Teast  = Twhole[mod(sNx+(1:OLx)-1, sNx*nPx )+1,1:sNy]
    Tsouth = Twhole[1:sNx, mod(sNy*nPy-OLy+(1:OLy)-1,sNy*nPy)+1]

    T[1:OLx,OLy+(1:sNy)] = Twest
    T[OLx+(1:sNx),sNy+OLy+(1:OLy)] = Tnorth
    T[sNx+OLx+(1:OLx),OLy+(1:sNy)] = Teast
    T[OLx+(1:sNx),1:OLy] = Tsouth

    sreq = Array(MPI.Request,4)
    for ip = 1:(nPx*nPy-1)
      ix = mod(ip,nPx)
      iy = round(Int64, (ip-ix)/nPx)

      # Sending West Boundary
      iwx = mod((ix-1), nPx)
      iwy = iy
      pendingx = (iwx+1)*sNx-OLx
      pendingy = (iwy)*sNy
      Twest = Twhole[pendingx+(1:OLx), pendingy+(1:sNy)]
      sreq[1] = MPI.Isend(Twest, ip, ip*4+1, comm)

      # Sending North Boundary
      iwx = ix
      iwy = mod(iy+1,nPy)
      pendingx = iwx*sNx
      pendingy = iwy*sNy
      Tnorth = Twhole[pendingx+(1:sNx), pendingy+(1:OLy)]
      sreq[2] = MPI.Isend(Tnorth, ip, ip*4+2, comm)

      # Sending East Boundary
      iwx = mod(ix+1,nPx)
      iwy = iy
      pendingx = iwx*sNx
      pendingy = iwy*sNy
      Twest = Twhole[pendingx+(1:OLx), pendingy+(1:sNy)]
      sreq[3] = MPI.Isend(Teast, ip, ip*4+3, comm)

      # Sending South Boundary
      iwx = ix
      iwy = mod(iy-1,nPy)
      pendingx = iwx*sNx
      pendingy = (iwy+1)*sNy-OLy
      Twest = Twhole[pendingx+(1:sNx), pendingy+(1:OLy)]
      sreq[4] = MPI.Isend(Tsouth, ip, ip*4+4, comm)

      MPI.Waitall!(sreq)

    end
  end
  return T
end

function MPI2D_GatherAll!(p :: MPI_2D_Config, T :: Array{Float64,3}, Twhole :: Array{Float64,3})
  for ii = 1:size(T,3)
    T1 = T[:,:,ii]
    Twhole1 = Twhole[:,:,ii]
    MPI2D_GatherAll!(p,T1, Twhole1)
  end
  return Twhole1
end

function MPI2D_DistributeCore!(p :: MPI_2D_Config, T :: Array{Float64,3}, Twhole :: Array{Float64,3})
  for ii = 1:size(T,3)
    T1 = T[:,:,ii]
    Twhole1 = Twhole[:,:,ii]
    MPI2D_DistributeCore!(p,T1,  Twhole1)
  end
  return Twhole1
end

function MPI2D_DistributeBoundary!(p :: MPI_2D_Config, T :: Array{Float64,3},  Twhole :: Array{Float64,3})
  for ii = 1:size(T,3)
    T1 = T[:,:,ii]
    Twhole1 = Twhole[:,:,ii]
    MPI2D_DistributeBoundary!(p,T1, Twhole1)
  end
  return Twhole1
end

export MPI_2D_Config
export MPI2D_GatherAll!
export MPI2D_DistributeCore!
export MPI2D_DistributeBoundary!


text:
module std_msgs

# importing RosTime
using ..CxxROS

using Cxx
using CxxStd: StdString

cxx"""
#include <std_msgs/String.h>
#include <std_msgs/Time.h>
"""

TVals = Union{}

for dtype in [ :String, :Time ]
    cppname = string("std_msgs::", dtype)
    cxxtval = Expr(:macrocall, Symbol("@cxxt_str"), cppname)
    symtval = Symbol(dtype, :Msg)

    @eval global const $(symtval) = $(cxxtval)
    @eval TVals = Union{TVals, $(symtval)}

    cxxcstr = Expr(:macrocall, Symbol("@icxx_str"), string(cppname, "();"))
    @eval (::Type{$(symtval)})() = $(cxxcstr)
end

const StdData = Union{String, StdString, RosTime}
const StdMsgs = TVals
#const StdMsgs = Union{StringMsg, TimeMsg}

Base.convert(::Type{S}, msg::StdMsgs) where {S<:StdData} = S(icxx"$msg.data;")

function Base.convert(::Type{T}, data::StdData) where {T<:StdMsgs}
    msg = T()
    icxx"$msg.data = $data;"
    msg
end

end #module

eval(std_msgs, :(String = StringMsg))
eval(std_msgs, :(Time	= TimeMsg))
export std_msgs

max_stars_repo_name: moyhig/CxxROS

max_issues_repo_name: moyhig/CxxROS

max_forks_repo_name: moyhig/CxxROS


text:
using ArcadeLearningEnvironment
using Test

setLoggerMode!(:error)

@testset "alien" begin
    ale = ALE_new()
    loadROM(ale, "alien")
    actionset = getLegalActionSet(ale)
    minactionset = getMinimalActionSet(ale)
    @test actionset == minactionset
    @test length(actionset) == 18
    for a in actionset
        act(ale, a)
    end
    @test false == game_over(ale)
    @test 3 == lives(ale)
    @test 18 == getFrameNumber(ale)
    @test 18 == getEpisodeFrameNumber(ale)
    @test length(getScreen(ale)) == 33600 == getScreenWidth(ale) * getScreenHeight(ale)
    getScreen(ale)
    s = getScreenRGB(ale)
    act(ale, Int32(0))
    s = getScreenGrayscale(ale)
    ALE_del(ale)
end

@testset "all roms" begin
    roms = getROMList()
    for rom in roms
        if rom in ["pacman", "surround"]
            @warn "Skipping $rom"
            continue
        end
        ale = ALE_new()
        setBool(ale, "color_averaging", true)
        setInt(ale, "frame_skip", Int32(4))
        setFloat(ale, "repeat_action_probability",
                 Float32(0.))
        @info "Testing $rom"
        loadROM(ale, rom)
        actionset = getLegalActionSet(ale)
        reset_game(ale)

        for i in 1:1000
            act(ale, rand(actionset))
            game_over(ale) && reset_game(ale)
        end

        width = getScreenWidth(ale)
        height = getScreenHeight(ale)
        getScreenGrayscale(ale)
        getScreenRGB(ale)
        getScreen(ale)

        @test true

        state_ref = cloneState(ale)
        pre_restore_state = encodeState(state_ref)
        deleteState(state_ref)
        state_ref = decodeState(pre_restore_state)
        restoreState(ale, state_ref)
        post_restore_state = ale |> cloneState |> encodeState
        @test pre_restore_state == post_restore_state
        deleteState(state_ref)

        state_ref = cloneSystemState(ale)
        pre_restore_state = encodeState(state_ref)
        state_ref = decodeState(pre_restore_state)
        restoreSystemState(ale, state_ref)
        post_restore_state = ale |> cloneSystemState |> encodeState
        @test post_restore_state == pre_restore_state
        deleteState(state_ref)

        ALE_del(ale)
    end
end


text:
module NonconvexCore

const debugging = Ref(false)
const show_residuals = Ref(false)

export  Model,
        DictModel,
        addvar!,
        getobjective,
        set_objective!,
        add_ineq_constraint!,
        add_eq_constraint!,
        add_sd_constraint!,
        decompress_symmetric,
        getmin,
        getmax,
        setmin!,
        setmax!,
        setinteger!,
        optimize,
        Workspace,
        KKTCriteria,
        IpoptCriteria,
        FunctionWrapper,
        Tolerance,
        GenericCriteria,
        KKTCriteria,
        ScaledKKTCriteria,
        IpoptCriteria

using Parameters, Zygote, ChainRulesCore, ForwardDiff, Requires
using SparseArrays, Reexport
using NamedTupleTools
@reexport using LinearAlgebra, OrderedCollections
using Reexport, Setfield
import JuMP, MathOptInterface
const MOI = MathOptInterface
using JuMP: VariableRef, is_binary, is_integer, has_lower_bound,
            has_upper_bound, lower_bound, upper_bound,
            start_value, ConstraintRef, constraint_object,
            AffExpr, objective_function, objective_sense

# General

include("utilities/params.jl")
include("functions/functions.jl")
include("functions/value_jacobian.jl")
include("functions/counting_function.jl")
include("common.jl")
include("utilities/callbacks.jl")
include("utilities/convergence.jl")

# Models

include("models/flatten.jl")
include("models/model.jl")
include("models/vec_model.jl")
include("models/dict_model.jl")
include("models/model_docs.jl")
include("models/jump.jl")
include("models/moi.jl")

end


text:
module TestArrayCommon

using Base.Test
using CategoricalArrays
using NullableArrays
using CategoricalArrays: DefaultRefType
using Compat

typealias String Compat.ASCIIString

# == currently throws an error for Nullables
(==) = isequal

# Tests of vcat of CategoricalArray and NullableCategoricalArray

# Test that mergelevels handles mutually compatible ordering
@test CategoricalArrays.mergelevels([6, 3, 4, 7], [2, 3, 5, 4], [2, 4, 8]) == ([6, 2, 3, 5, 4, 7, 8], true)
@test CategoricalArrays.mergelevels([6, 3, 4, 7], [2, 3, 6, 5, 4], [2, 4, 8]) == ([6, 3, 4, 7, 2, 5, 8], false)

for (CA, A) in ((CategoricalArray, Array), (NullableCategoricalArray, NullableArray))
    # Test that vcat of compact arrays uses a reftype that doesn't overflow
    a1 = 3:200
    a2 = 300:-1:100
    ca1 = CA(a1)
    ca2 = CA(a2)
    cca1 = compact(ca1)
    cca2 = compact(ca2)
    r = vcat(cca1, cca2)
    @test r == A(vcat(a1, a2))
    @test isa(cca1, CA{Int, 1, UInt8})
    @test isa(cca2, CA{Int, 1, UInt8})
    @test isa(r, CA{Int, 1, CategoricalArrays.DefaultRefType})
    @test isa(vcat(cca1, ca2), CA{Int, 1, CategoricalArrays.DefaultRefType})
    @test ordered(r) == false
    @test levels(r) == collect(3:300)

    # Test vcat of multidimensional arrays
    a1 = Array{Int}(2, 3, 4, 5)
    a2 = Array{Int}(3, 3, 4, 5)
    a1[1:end] = (length(a1):-1:1) + 2
    a2[1:end] = (1:length(a2)) + 10
    ca1 = CA(a1)
    ca2 = CA(a2)
    cca1 = compact(ca1)
    cca2 = compact(ca2)
    r = vcat(cca1, cca2)
    @test r == A(vcat(a1, a2))
    @test isa(r, CA{Int, 4, CategoricalArrays.DefaultRefType})
    @test ordered(r) == false
    @test levels(r) == collect(3:length(a2)+10)

    # Test concatenation of mutually compatible levels
    a1 = ["Young", "Middle"]
    a2 = ["Middle", "Old"]
    ca1 = CA(a1, ordered=true)
    ca2 = CA(a2, ordered=true)
    levels!(ca1, ["Young", "Middle"])
    levels!(ca2, ["Middle", "Old"])
    r = vcat(ca1, ca2)
    @test r == A(vcat(a1, a2))
    @test levels(r) == ["Young", "Middle", "Old"]
    @test ordered(r) == true

    # Test concatenation of conflicting ordering. This drops the ordering
    a1 = ["Old", "Young", "Young"]
    a2 = ["Old", "Young", "Middle", "Young"]
    ca1 = CA(a1, ordered=true)
    ca2 = CA(a2, ordered=true)
    levels!(ca1, ["Young", "Middle", "Old"])
    # ca2 has another order
    r = vcat(ca1, ca2)
    @test r == A(vcat(a1, a2))
    @test levels(r) == ["Young", "Middle", "Old"]
    @test ordered(r) == false
end

end


text:
module RandomMatrices
importall Distributions
using Combinatorics

import Base.rand

#If the GNU Scientific Library is present, turn on additional functionality.
_HAVE_GSL = try
  using GSL
  true
catch
   false
end

export bidrand,    #Generate random bidiagonal matrix
       tridrand,   #Generate random tridiagonal matrix
       sprand,     #Generate random sparse matrix
       eigvalrand, #Generate random set of eigenvalues
       eigvaljpdf  #Eigenvalue joint probability density

typealias Dim2 (Int, Int) #Dimensions of a rectangular matrix

#Classical Gaussian matrix ensembles
include("GaussianEnsembles.jl")

# Classical univariate distributions
include("densities/Semicircle.jl")
include("densities/TracyWidom.jl")

# Ginibre
include("Ginibre.jl")

#Generating matrices of Haar measure
include("Haar.jl")
include("HaarMeasure.jl")

#Fast histogrammer for matrix eigenvalues - hist_eig
include("FastHistogram.jl")

#Formal power series
include("FormalPowerSeries.jl")

#Statistical tests based on random matrix theory
include("StatisticalTests.jl")

#Stochastic processes
include("StochasticProcess.jl")

#Invariant ensembles
if filemode(Pkg.dir("ApproxFun")) != 0
    include("InvariantEnsembles.jl")
end

end


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.Model{Float64}())
config = MOIT.TestConfig()

bridged_mock = MOIB.Variable.SOCtoRSOC{Float64}(mock)

@testset "soc1v" begin
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/√2 + 1/2, 1/√2 - 1/2, 1/√2],
        (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[-√2]])
    MOIT.soc1vtest(bridged_mock, config)

    ceqs = MOI.get(mock, MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{Float64}, MOI.Zeros}())
    @test length(ceqs) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), ceqs[1], "ceq")

    @testset "Test mock model" begin
        var_names = ["a", "b", "c"]
        MOI.set(
            mock, MOI.VariableName(),
            MOI.get(mock, MOI.ListOfVariableIndices()), var_names)
        rsocs = MOI.get(mock, MOI.ListOfConstraintIndices{
            MOI.VectorOfVariables, MOI.RotatedSecondOrderCone}())
        @test length(rsocs) == 1
        MOI.set(mock, MOI.ConstraintName(), rsocs[1], "rsoc")

        invs2 = 1 / √2
        s = """
        variables: a, b, c
        rsoc: [a, b, c] in MathOptInterface.RotatedSecondOrderCone(3)
        ceq: [$invs2*a + $invs2*b + -1.0] in MathOptInterface.Zeros(1)
        maxobjective: $invs2*a + -$invs2*b + c
        """
        model = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model, s)
        MOIU.test_models_equal(mock, model, var_names, ["rsoc", "ceq"])
    end

    @testset "Test bridged model" begin
        var_names = ["x", "y", "z"]
        MOI.set(
            bridged_mock, MOI.VariableName(),
            MOI.get(bridged_mock, MOI.ListOfVariableIndices()), var_names)
        socs = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{
            MOI.VectorOfVariables, MOI.SecondOrderCone}())
        @test length(socs) == 1
        MOI.set(bridged_mock, MOI.ConstraintName(), socs[1], "soc")

        s = """
        variables: x, y, z
        soc: [x, y, z] in MathOptInterface.SecondOrderCone(3)
        ceq: [x + -1.0] in MathOptInterface.Zeros(1)
        maxobjective: y + z
        """
        model = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model, s)
        MOIU.test_models_equal(bridged_mock, model, var_names, ["soc", "ceq"])
    end

    @testset "Delete" begin
        xyz = MOI.get(bridged_mock, MOI.ListOfVariableIndices())

        message = string("Cannot delete variable as it is constrained with other",
                         " variables in a `MOI.VectorOfVariables`.")
        for i in eachindex(xyz)
            err = MOI.DeleteNotAllowed(xyz[i], message)
            @test_throws err MOI.delete(bridged_mock, xyz[i])
        end

        test_delete_bridged_variables(bridged_mock, xyz, MOI.SecondOrderCone, 3, (
            (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone, 0),
        ))
    end
end


text:
struct PengRobinson{S,T,A} <: CubicModel
    type::S
    tc::T
    pc::T
    ω::T
    mw::T
    vc::Union{T,Nothing}
    _a::T
    _b::T
    aij::A

    function PengRobinson(tc,pc,ω,mw,vc=nothing,aij = nothing)
        if length(tc) == 1
            type = SINGLE()
        else
            type = MULTI()
        end
        S = typeof(type)
        tc,pc,ω,mw = promote(tc,pc,ω,mw)
        T = typeof(tc)
        A = typeof(aij)
        _a = 0.45723553 .*((RGAS .*tc).^2) ./pc
        _b = 0.07779607 .*((RGAS .*tc)) ./pc
        return new{S,T,A}(type,tc,pc,ω,mw,vc,_a,_b,aij)
    end
end

volume_solver_type(model::PengRobinson) = CubicRoots()
single_sat_Approx(model::PengRobinson{SINGLE}) =PRSatApprox(model)
mol_density(model::PengRobinson{SINGLE},::CriticalPoint,unit=u"mol/(m^3)") = convert_unit(u"mol/L",unit,inv(only(model.vc)))
pressure(model::PengRobinson{SINGLE},::CriticalPoint,unit=u"Pa") = convert_unit(u"Pa",unit,only(model.pc))
temperature(model::PengRobinson{SINGLE},::CriticalPoint,unit=u"K") = convert_unit(u"K",unit,only(model.tc))
mol_volume(model::PengRobinson{SINGLE},::CriticalPoint,unit=u"m^3/mol") = convert_unit(u"m^3/mol",unit,only(model.vc))
acentric_factor(model::PengRobinson{SINGLE}) = only(model.ω)
molecular_weight(model::PengRobinson{SINGLE}) = only(model.mw)
molecular_weight(model::PengRobinson{MULTI}) = model.mw

function mol_density(model::PengRobinson{MULTI},::CriticalPoint,unit=u"mol/(m^3)")
    return convert_unit.(u"mol/L",unit,1 ./ model.vc)
end

function pressure(model::PengRobinson{MULTI},::CriticalPoint,unit=u"Pa")
    return convert_unit.(u"Pa",unit,model.pc)
end

function temperature(model::PengRobinson{MULTI},::CriticalPoint,unit=u"K")
    return convert_unit.(u"K",unit,model.tc)
end

function mol_volume(model::PengRobinson{MULTI},::CriticalPoint,unit=u"m^3/mol")
    return convert_unit.(u"m^3/mol",unit,model.vc)
end

function acentric_factor(model::PengRobinson{MULTI})
    return model.ω
end

function PengRobinson(;tc,pc,ω,mw,vc=nothing,aij=nothing)
    return PengRobinson(tc,pc,ω,mw,vc,aij)
end

function PengRobinson(model::ThermoModel)
    tc = temperature(model,CriticalPoint())
    pc = pressure(model,CriticalPoint())
    ω = acentric_factor(model)
    mw = molecular_weight(model)
    return PengRobinson(tc,pc,ω,mw)
end

function cubic_aα(model::PengRobinson,t,i,j)

    m_poly = (0.37464,1.54226,0.26992)
    _1 = one(t)
    aᵢ =model._a[i]
    tcᵢ = model.tc[i]
    sqrt_trᵢ = min(sqrt(t/tcᵢ),_1) #α(t) is one for supercritical values
    ωᵢ = model.ω[i]
    mᵢ = evalpoly(ωᵢ,m_poly)
    if i === j
        return  aᵢ* (_1+mᵢ * √(_1-sqrt_trᵢ))^2
    else
        aⱼ =model._a[j]
        tcⱼ = model.tc[j]
        sqrt_trⱼ = min(sqrt(t/tcⱼ),_1)
        ωⱼ = model.ω[j]
        mⱼ = evalpoly(ωⱼ,m_poly)
        sqrt_αᵢ= (_1+mᵢ * √(_1-sqrt_trᵢ))
        sqrt_αⱼ= (_1+mⱼ * √(_1-sqrt_trⱼ))
        return sqrt(aᵢ*aⱼ)*sqrt_αᵢ*sqrt_αⱼ
    end
end

function cubic_ab(mt::SinglePT,model::PengRobinson{SINGLE},v,t)
    a = cubic_aα(model,t,1,1)
    b = only(model._b)
    return a,b
end

function cubic_ab(mt::MultiPT,model::PengRobinson{MULTI},p,t,x)
    #two options to introduce alpha:
    #here: it will allocate, but less ops
    bi = model._b
    b = dot(bi,x)
    #here: it will not allocate, but more ops
    sss= (i,j)->cubic_aα(model,t,i,j)
    a = cubic_mixing_rule(sss, x, model.aij)
    return a,b
end

function cubic_abp(mt::SingleVT,model::PengRobinson{SINGLE},v,t)
    a,b = cubic_ab(QuickStates.pt(),model,v,t) #v is ignored
    _1 = one(b)
    denom = evalpoly(v,(-b*b,2*b,_1))
    p = RGAS*t/(v-b) - a/denom
    return a,b,p
end

function cubic_abp(mt::MultiVT,model::PengRobinson{MULTI},v,t,x)
    a,b = cubic_ab(QuickStates.ptx(),model,v,t,x) #v is ignored
    _1 = one(b)
    denom = evalpoly(v,(-b*b,2*b,_1))
    p = RGAS*t/(v-b) - a/denom
    return a,b,p
end

function fugacity_coeff_impl(mt::SingleVT,model::PengRobinson{SINGLE},v,t)
    a,b,p =  cubic_abp(mt,model,v,t)
    RTinv = 1/(RGAS*t)
    A = a*p*RTinv*RTinv
    B = b*p*RTinv
    z = p*v*RTinv
     _1 = one(z)
     logϕ = z - _1 - log(z-B) - A/z
end

const PRΔ1 = 1+√2
const PRΔ2 = 1-√2
const ΔPRΔ = 2*√2

function  αR_impl(mt::MultiVT,model::PengRobinson{MULTI},rho,t,x)
    R = RGAS
    RTinv = 1/(RGAS*t)
    v = inv(rho)
    a,b,p =  cubic_abp(mt,model,v,t,x)
    -log(1-b*rho) - a*RTinv*log((PRΔ1*b*rho+1)/(PRΔ2*b*rho+1))/(ΔPRΔ*b)
end

function cubic_poly(mt::SinglePT,model::PengRobinson{SINGLE},p,t)
    a,b = cubic_ab(QuickStates.pt(),model,p,t)
    RTinv = 1/(RGAS*t)
    A = a*p*RTinv*RTinv
    B = b*p*RTinv
    k0 = B*(B*(B+1.0)-A)
    k1 = -B*(2*B+1.0) + A
    k2 = -1.0
    k3 = 1.0
    return (k0,k1,k2,k3)
end

function cubic_poly(mt::MultiPT,model::PengRobinson{MULTI},p,t,x)
    a,b = cubic_ab(QuickStates.ptx(),model,p,t,x)
    RTinv = 1/(RGAS*t)
    A = a*p*RTinv*RTinv
    B = b*p*RTinv
    k0 = B*(B*(B+1.0)-A)
    k1 = -B*(2*B+1.0) + A
    k2 = -1.0
    k3 = 1.0
    return (k0,k1,k2,k3)
end

struct PRSatApprox{M} <: SaturationModel
    model::M
    function PRSatApprox(model::CubicModel)
        T = typeof(model)
        return new{T}(model)
    end
end

function PRSatApprox(tc,pc,mw=copy(pc),vc=nothing,ω=nothing,aij = nothing)
    model = PengRobinson(tc,pc,mw,vc,ω,aij)
    return PRSatApprox(model)
end

function PRSatApprox(;tc,pc,mw=copy(pc),vc=nothing,ω=nothing,aij=nothing)
    model =  PengRobinson(tc,pc,mw,vc,ω,aij)
    return PRSatApprox(model)
end

function pressure_impl(mt::SingleSatT,model::PRSatApprox,t)
    a = -2.605272488488440
    b = -9.017571450539830
    c = -19.896014683288000
    d = 0.579677284391001
    e = 22.501011849124900
    f = -0.041440165126830
    tc = only(model.model.tc)
    pc = only(model.model.pc)
    tr = t/tc
    atr = tc/t
    num = evalpoly(atr,(a,c,e))
    denom = evalpoly(atr,(1.0,b,d,f))
    pr = exp(num/denom)*tr
    p = pr*pc
end

initial_temperature(model::PRSatApprox,p) = critical_sat_interpolation(model,p)


text:
using Flux
using Random

"""
SequentialData
    n_data: number of data
    vstr: vector of strings
    charlist: list of char used in strings and added <SOS>, <EOS>, <PAD>
    n_seq: maximum number of sequence of string
    data: processed sequencial data
"""
mutable struct SequentialData
    n_data::Integer
    vstr::Vector{String}
    charlist::Vector{String}
    n_char::Integer
    n_seq::Integer
    data
    function SequentialData(vstr::Vector{String}, charlist::Vector{String})
        n_data = length(vstr)
        n_char = length(charlist)
        n_seq = maximum(map(length, vstr)) + 1
        vstr_token = map(str -> push!(string.(collect(str)), "<EOS>"), vstr)
        data = map(str_token->append!(str_token, collect(Iterators.repeated("<PAD>", n_seq-length(str_token)))), vstr_token)
        data_array = reduce(hcat, data)
        data_sequential = [Flux.onehotbatch(data_array[i, :], charlist) for i in 1:n_seq]
        new(n_data, vstr, charlist, n_char, n_seq, data_sequential)
    end
end

"""
make dataset chunk, SequencialData splited by batch_size
"""
function get_shuffled_traindata(data::SequentialData, batch_size::Integer)
    n_data = length(data.vstr)
    vindices = collect(Iterators.partition(randperm(n_data), batch_size))
    return [([x[:, indices] for x in data.data],) for indices in vindices]
end

"""
from One-hot matrix to index of charlist
"""
onehot_to_indices(onehot) = reduce(vcat, mapslices(x->argmax(x), onehot, dims=1))

struct Encoder
    embed
    gru
    function Encoder(features, n_char)
        embed = Flux.glorot_uniform(features, n_char) |> param
        gru = GRU(features, features)
        return new(embed, gru)
    end
end
Flux.@treelike Encoder # enable params(e::Encoder)

function (e::Encoder)(x)
    Flux.reset!(e.gru)
    vindices = map(onehot_to_indices, x)
    x = map(x1 -> e.embed[:, x1], vindices)
    hidden = e.gru.(x)[end]
    return hidden
end

struct Decoder
    embed
    gru
    dense
    n_char::Integer
    n_seq::Integer
    function Decoder(features, n_char, n_seq)
        embed = Flux.glorot_uniform(features, n_char) |> param
        gru = GRU(features, features)
        dense = Dense(features, n_char, relu)
        return new(embed, gru, dense, n_char, n_seq)
    end
end
Flux.@treelike Decoder

function (d::Decoder)(h)
    Flux.reset!(d.gru)
    _, batch_size = size(h)
    d.gru.state = h
    EOS = d.embed[:, [d.n_char-1 for _ in 1:batch_size]]
    x = d.gru(EOS)
    x = d.dense(x)
    categorical = softmax(x)
    for i in 1:d.n_seq-1
        indices = onehot_to_indices(x)
        x2 = d.embed[:, indices]
        x = d.gru(x2) |> d.dense
        categorical = hcat(categorical, softmax(x))
    end
    categorical
end

struct Seq2Seq
    encoder::Encoder
    decoder::Decoder
    n_char::Integer
    function Seq2Seq(features::Integer, data::SequentialData)
        encoder = Encoder(features, data.n_char)
        decoder = Decoder(features, data.n_char, data.n_seq)
        return new(encoder, decoder, data.n_char)
    end
end
Flux.@treelike Seq2Seq

function loss_func(m::Seq2Seq)
    loss_total = 0
    weight = ones(m.n_char)
    weight[end] = 0 # <PAD> weight to zero
    return function loss(seq_x)
        hidden_state = m.encoder(reverse(seq_x))
        predict = m.decoder(hidden_state)
        t = reduce(hcat, seq_x)
        loss_batch = Flux.crossentropy(predict, t, weight=weight)
        # global loss_total += loss_batch
        return loss_batch
    end
end

max_stars_repo_name: Soyukke/FluxModels

max_issues_repo_name: Soyukke/FluxModels

max_forks_repo_name: Soyukke/FluxModels


text:
using Test
include("../../utils/funcion-objetivo.jl")

@testset "Función objetivo" begin
    data = [-1 0; 1 1; -1 -3; 1 -2; 9 -1; 1 2; 10 20]
    class = ['-', '+',  '-',   '-',  '-', '+' , '+' ]
    F=CrearFuncionObjetivo(data, class, 0.5, 0.3)
    @test 0<= F([0.1,0.5]) <= 100
end


text:
export propagate
export fresnel_kernel, rs_kernel
export point_source_propagate
export lens_propagate, four_f_propagate, mla_propagate

"""
    fresnel_kernel(fx, fy, z, λ, n=1.0)

Calculates the Fresnel propagation kernel in Fourier space.
`fx`, `f_y` are the spatial frequency. `z` is the propagation distance.
`λ` is the wavelength and `n` the refractive index.
"""
function fresnel_kernel(fx, fy, z, λ, n=1)
    κ = calc_κ(λ, n)
    return @. exp(1im * 2 * π * κ * z) * exp(-1im * π * z / κ * (fx^2 + fy^2))
end

"""
    rs_kernel(fx, fy, z, λ, n=1.0)

Calculates the Rayleigh-Sommerfeld propagation kernel in Fourier space.
`fx`, `f_y` are the spatial frequencies. `z` is the propagation distance.
`λ` is the wavelength and `n` the refractive index.
"""
function rs_kernel(fx, fy, z, λ=λ0, n=1)
    κ = calc_κ(λ, n)
    c = complex(κ^2 - (fx^2 + fy^2))
    return exp(1im * z * 2 * π * sqrt(c))
end

"""
    propagate!(arr, L, z; kernel=rs_kernel, λ=550e-9, n=1)

Propagate an electric field `arr` with a array size (in physical dimensions like meter etc) of `L` over the distance `z`.
Per default `kernel=rs_kernel` (Rayleigh-Sommerfeld) is the propagation kernel.
`λ` is the wavelength and `n` the refractive index.
"""
function propagate!(arr, L, z; kernel=rs_kernel, λ=λ0, n=1, fft_plan=plan_fft!(arr))
    # trivial propagation
    if iszero(z)
        # return a copy instead of arr
        # therefore we can rely on, that we don't modify the initial array
        return copy(arr)
    end

    # array with the frequencies in Fourier space
    freq_x = to_gpu_or_cpu(arr, Zygote.@ignore eltype(arr).(fftfreq(size(arr)[2], size(arr)[2] / L)'))
    freq_y = to_gpu_or_cpu(arr, Zygote.@ignore eltype(arr).(fftfreq(size(arr)[1], size(arr)[1] / L)))
    arr_ft = fft_plan * arr

    κ = eltype(arr)(calc_κ(λ, n))
    c_sqrt = Zygote.@ignore (κ^2 .- (freq_x.^2 .+ freq_y.^2))
    c = Zygote.@ignore 1im * eltype(arr)(z * 2π)
    c_exp = Zygote.@ignore (exp.(c .* sqrt.(c_sqrt)))
    out_ft = arr_ft .* c_exp
    #out_ft = arr_ft .* kernel.(freq_x, freq_y, Ref(z), Ref(λ), Ref(n))
    p_inv = Zygote.@ignore inv(fft_plan)
    out = p_inv.scale .* (p_inv.p * out_ft)
    return out
end

function propagate(arr, L, z; kernel=rs_kernel, λ=λ0, n=1)
    propagate!(complex.(arr), L, z, kernel=kernel, λ=λ, n=n)
end

function propagate_tilted_plane(arr, L, z, tilting_ϕ, dim=1; kernel=rs_kernel, λ=λ0, n=1, fft_plan=plan_fft!(arr))
    Zygote.@ignore GC.gc(true)
    Δzs = z .+ fftpos(L, size(arr, 1)) .* tan(deg2rad(tilting_ϕ))

    arr2 = propagate!(arr, L, Δzs[1], fft_plan=fft_plan)
    if dim == 1
        out = arr2[1:1, :]
    else
        out = arr2[:, 1:1]
    end
    for (i, Δz) in enumerate(Δzs)
        Δz = Δzs.step
        if isone(i)
            continue
        end
        arr2 = propagate!(arr2, L, Δz, kernel=rs_kernel, λ=λ, n=n, fft_plan=fft_plan)
        if dim == 1
            out = vcat(out, arr2[i:i, :])
        else
            out = hcat(out, arr2[:, i:i])
        end
    end
    return out
end

"""
    point_source_propagate(L, size, point; λ=550e-9, n=1, dtype=ComplexF64)

Propagate a point source in a field of width `L` with array size `size`
over a distance `z` in a medium with refractive index `n`.
Point source is located at position `point(x, y, z)` and then propagated to
a location z=0. `z<0` in `point` means forward propagation.

This is based on the analytical solution in real space and not on Fourier
space propagation. The latter one suffers from artifacts while microscopic large `L`.
This function should be always preferred for point sources.
"""
function point_source_propagate(L, size, point::Point; λ=λ0, n=1, dtype=ComplexF64)
    x0, y0, z = point.x, point.y, point.z
    out = zeros(dtype, size)
    k = calc_k(λ, n)
    # if we are at the plane of the point source
    # return numerical delta
    if iszero(z)
        j = argmin(map(x -> abs(x - x0), fftpos(L, size[2])))
        i = argmin(map(y -> abs(y - y0), fftpos(L, size[1])))
        out[i, j] = 1
        return out
    end
    # calculate values
    @inbounds for (j, x) in enumerate(fftpos(L, size[2]))
        for (i, y) in enumerate(fftpos(L, size[1]))
            # z>0 means that we back propagate the PS to z=0
            r = dtype(sign(-z) * sqrt((x-x0)^2 + (y-y0)^2 + z^2))
            out[i, j] = 1/r .* exp(1im * k * r);
        end
    end
    out = out ./ out[argmax(abs2.(out))]
    return out
end

"""
    lens_propagate(arr, L, f; λ=550e-9, n=1, d=nothing)

Propagate an electrical field of width `L` at a distance `d` in
front of a lens with focal length `f` in medium with refractive index `n`
to the back focal plane of the lens.
Per default `d=nothing` meaning that we set `d=f`.
Assume same sampling and size in both dimensions.
As return parameter we get the resulting field and the new field size.
Essentially the lens (assuminig infinite large of that lens) performs
a scaled Fourier transform. Therefore we get a new field size.

 # Reference
Based on:
* "Computational Fourier Optics. A MATLAB Tutorial", D. Voelz, (2011).
* Goodman, Joseph W. Introduction to Fourier optics
"""
function lens_propagate(arr, L, f; λ=λ0, n=1, d=nothing, fft_plan=plan_fft!(complex.(arr)))
    if isnothing(d)
        d = f
    end
    κ = calc_κ(λ, n)
    # shift to center again
    if typeof(arr) <: Array{<:Real}
        arr = complex.(arr)
    end
    out_f = fftshift(fft_plan * ifftshift(arr))
    dx = L / size(arr)[2]
    dy = L / size(arr)[1]

    # lens performs scaled fourier transform. Therefore new field size
    L_new = λ * f / dx
    aux_x =  Zygote.@ignore fftpos(L_new, size(arr)[2])'
    aux_y = Zygote.@ignore fftpos(L_new, size(arr)[1])
    x = Zygote.@ignore to_gpu_or_cpu(arr, aux_x)
    y = Zygote.@ignore to_gpu_or_cpu(arr, aux_y)
    c_exp = eltype(arr)(π * κ / f * (1 - d / f))
    c = Zygote.@ignore  .- 1im ./ (λ .* f) .* cis.(c_exp .* (x.^2 .+ y.^2))
    out = out_f .* c .* eltype(out_f)(dx * dy)

    return out, L_new
end

"""
    four_f_propagate(arr, L, f1, f2, NA; fft_plan=fft_plan, d=nothing)

Propagate the electrical field `arr` (field size `L`) from front
focal plane to the back focal plane of a 4f system.
The focal length of the first lens is `f1` and the second lens `f`.
Magnification of the system is then `f2/f1`.
`d` is the distance from lens to the point of the `arr` in front
of the lens. If `d == nothing` then we propagate from object plane to
image plane.

"""
function four_f_propagate(arr, L, f1, f2, NA; d=nothing, fft_plan=plan_fft!(complex.(arr)))
	radius = NA * 2 * f1 / 2

    out = lens_propagate(arr, L, f1, fft_plan=fft_plan, d=d)
    E1 = out[1]
    L1 = out[2]
    # check that the field size is large enough that the pupil fits
    @assert radius < L1 / 2
	E1_ = circ(E1, L1, radius)

    # tuple of E2 and L2
    out2 = lens_propagate(E1_, L1, f2, fft_plan=fft_plan)
    return out2
end

"""
    mla_propagate(E, L, fmla, p_mla, z1, z2)

Propagate an electrical field `E` placed a distance `z1` in front
of a microlns array, to a distance `z2` behind the microlens array.
`L` is the size of the field and `fmla` is the focal length of the lenslets.
`p_mla` is microlenses pitch.
"""
function mla_propagate(E, L, fmla, p_mla, z1, z2)

    E_front = propagate(E, L, z1)
	mla = micro_lens_array(fmla, L, p_mla, size(E))
	E_mla = E_front .* mla
	E_behind = propagate(E_mla, L, z2);
	return E_behind, real(mla)
end

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
module GraphGrinder

using LightGraphs, MetaGraphs
using MLDataPattern, Flux, Mill

include("graph2bag.jl")
include("matchvertex2edges.jl")
include("message_passing.jl")
include("simple_explainer.jl")
include("graphutils.jl")

export graph2bag, vertices2bag, SimpleExplainer, get_vertex_labels, strip_lone_vertices, component_including_vertex

end # module


text:
##############################################################################
##
## Find common prefixes (up to lim. -1 means Inf)
## Assumes length(s1) <= length(s2)
##############################################################################

function common_prefix(s1::AbstractStringorGraphemeIterator, s2::AbstractStringorGraphemeIterator, lim::Integer = -1)
    start1 = start(s1)
    start2 = start(s2)
    l = 0
    while !done(s1, start1) && (l < lim || lim < 0)
        ch1, nextstart1 = next(s1, start1)
        ch2, nextstart2 = next(s2, start2)
        ch1 != ch2 && break
        l += 1
        start1, start2 = nextstart1, nextstart2
    end
    return l, start1, start2
end

##############################################################################
##
## Hamming
##
##############################################################################

function evaluate(dist::Hamming, s1::AbstractStringorGraphemeIterator, s2::AbstractStringorGraphemeIterator, len1::Integer, len2:: Integer)
    count = 0
    for (ch1, ch2) in zip(s1, s2)
        count += ch1 != ch2
    end
    count += len2 - len1
    return count
end

##############################################################################
##
## Levenshtein
## Source: http://blog.softwx.net/2014/12/optimizing-levenshtein-algorithm-in-c.html
##
##############################################################################

type Levenshtein <: SemiMetric end
function evaluate(dist::Levenshtein, s1::AbstractStringorGraphemeIterator, s2::AbstractStringorGraphemeIterator, len1::Integer, len2::Integer)

    # prefix common to both strings can be ignored
    k, start1, start2 = common_prefix(s1, s2)
    done(s1, start1) && return len2 - k

    # distance initialized to first row of matrix
    # => distance between "" and s2[1:i}
    v0 = Array(Int, len2 - k)
    @inbounds for i2 in 1:(len2 - k)
        v0[i2] = i2
    end
    current = zero(0)
    state1 = start1
    i1 = 0
    while !done(s1, state1)
        i1 += 1
        ch1, state1 = next(s1, state1)
        left = (i1 - 1)
        current = (i1 - 1)
        state2 = start2
        i2 = 0
        while !done(s2, state2)
            i2 += 1
            ch2, state2 = next(s2, state2)
            #  update
            above, current, left = current, left, v0[i2]
            if ch1 != ch2
                # substitution
                current = min(current + 1,
                                above + 1,
                                left + 1)
            end
            v0[i2] = current
        end
    end
    return current
end

##############################################################################
##
## Damerau Levenshtein
## Source: http://blog.softwx.net/2015/01/optimizing-damerau-levenshtein_15.html
##
##############################################################################

type DamerauLevenshtein <: SemiMetric end

function evaluate(dist::DamerauLevenshtein, s1::AbstractStringorGraphemeIterator, s2::AbstractStringorGraphemeIterator, len1::Integer, len2::Integer)

    # prefix common to both strings can be ignored
    k, start1, start2 = common_prefix(s1, s2)
    done(s1, start1) && return len2 - k

    v0 = Array(Int, len2 - k)
    @inbounds for i2 in 1:(len2 - k)
        v0[i2] = i2
    end
    v2 = Array(Int, len2 - k)

    ch1, = next(s1, start1)
    current = 0
    state1 = start1
    i1 = 0
    while !done(s1, state1)
        i1 += 1
        prevch1 = ch1
        ch1, state1 = next(s1, state1)
        ch2, = next(s2, start2)
        left = (i1 - 1)
        current = i1
        nextTransCost = 0
        state2 = start2
        i2 = 0
        while !done(s2, state2)
            i2 += 1
            prevch2 = ch2
            ch2, state2 = next(s2, state2)
            above = current
            thisTransCost = nextTransCost
            nextTransCost = v2[i2]
            # cost of diagonal (substitution)
            v2[i2] = current = left
            # left now equals current cost (which will be diagonal at next iteration)
            left = v0[i2]
            if ch1 != ch2
                # insertion
                if left < current
                    current = left
                end
                # deletion
                if above < current
                    current = above
                end
                current += 1
                if i1 != 1 && i2 != 1 && ch1 == prevch2 && prevch1 == ch2
                    thisTransCost += 1
                    if thisTransCost < current
                        current = thisTransCost
                    end
                end
            end
            v0[i2] = current
        end
    end
    return current
end

##############################################################################
##
## Jaro
##
##############################################################################

type Jaro <: SemiMetric end

function evaluate(dist::Jaro, s1::AbstractStringorGraphemeIterator, s2::AbstractStringorGraphemeIterator, len1::Integer, len2::Integer)
    # if len2 == 0, m = 0 so should be 1.0 according to wikipedia. Nope.
    len2 == 0 && return 0.0

    maxdist = max(0, div(len2, 2) - 1)
    m = 0 # matching characters
    t = 0 # half number of transpositions
    flag = fill(false, len2)
    prevpos = 0

    i1 = 0
    startstate2 = start(s2)
    starti2 = 0
    for ch1 in s1
        i1 += 1
        if starti2 < i1 - maxdist - 1
            startstate2 = nextind(s2, startstate2)
            starti2 += 1
        end
        i2 = starti2
        state2 = startstate2
        while !done(s2, state2) && i2 < i1 + maxdist
            ch2, state2 = next(s2, state2)
            i2 += 1
            if ch1 == ch2 && !flag[i2]
                m += 1
                # if match is before the index of previous match
                if i2 < prevpos
                    t += 1
                end
                prevpos = max(i2, prevpos)
                flag[i2] = true
                break
            end
        end
    end
    m == 0.0 && return 1.0
    score = (m / len1 + m / len2 + (m - t) / m) / 3.0
    return 1.0 - score
end

jaro(s1::AbstractStringorGraphemeIterator, s2::AbstractStringorGraphemeIterator) = evaluate(Jaro(), s1, s2)


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function deredd(Eby::T, by::T, m1::T, c1::T, ub::T) where {T<:AbstractFloat}
    Rm1 = -0.33
    Rc1 = 0.19
    Rub = 1.53
    Eby0 = Eby > 0 ? Eby : zero(T)
    return by - Eby0, m1 - Rm1*Eby0, c1 - Rc1*Eby0, ub - Rub*Eby0
end

"""
    deredd(Eby, by, m1, c1, ub) -> by0, m0, c0, ub0

### Purpose ###

Deredden stellar Stromgren parameters given for a value of E(b-y)

### Arguments ###

* `Eby`: color index E(b-y), scalar (E(b-y) = 0.73*E(B-V))
* `by`: b-y color (observed)
* `m1`: Stromgren line blanketing parameter (observed)
* `c1`: Stromgren Balmer discontinuity parameter (observed)
* `ub`: u-b color (observed)

All arguments can be either scalars or arrays all of the same length.

### Output ###

The 4-tuple `(by0, m0, c0, ub0)`.

* `by0`: b-y color (dereddened)
* `m0`: line blanketing index (dereddened)
* `c0`: Balmer discontinuity parameter (dereddened)
* `ub0`: u-b color (dereddened)

These are scalars or arrays of the same length as the input arguments.

### Example ###

```jldoctest
julia> using AstroLib

julia> deredd(0.5, 0.2, 1.0, 1.0, 0.1)
(-0.3, 1.165, 0.905, -0.665)
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
deredd(Eby::Real, by::Real, m1::Real, c1::Real, ub::Real) =
    deredd(promote(float(Eby), float(by), float(m1), float(c1), float(ub))...)

function deredd(Eby::AbstractArray{E}, by::AbstractArray{<:Real},
                m1::AbstractArray{<:Real}, c1::AbstractArray{<:Real},
                ub::AbstractArray{<:Real}) where {E<:Real}
    @assert length(Eby) == length(by) == length(m1) == length(c1) == length(ub)
    typeeby = float(E)
    by0 = similar(Eby, typeeby)
    m0  = similar(Eby, typeeby)
    c0  = similar(Eby, typeeby)
    ub0 = similar(Eby, typeeby)
    for i in eachindex(Eby)
        by0[i], m0[i], c0[i], ub0[i] = deredd(Eby[i], by[i], m1[i], c1[i], ub[i])
    end
    return by0, m0, c0, ub0
end


text:
"""
A simple example of creating a server with HTTP.jl. It handles creating, deleting,
updating, and retrieving Animals from a dictionary through 4 different routes
"""
using HTTP

# modified Animal struct to associate with specific user
mutable struct Animal
    id::Int
    userId::Base.UUID
    type::String
    name::String
end

# use a plain `Dict` as a "data store"
const ANIMALS = Dict{Int, Animal}()
const NEXT_ID = Ref(0)
function getNextId()
    id = NEXT_ID[]
    NEXT_ID[] += 1
    return id
end

# "service" functions to actually do the work
function createAnimal(req::HTTP.Request)
    animal = JSON2.read(IOBuffer(HTTP.payload(req)), Animal)
    animal.id = getNextId()
    ANIMALS[animal.id] = animal
    return HTTP.Response(200, JSON2.write(animal))
end

function getAnimal(req::HTTP.Request)
    animalId = HTTP.URIs.splitpath(req.target)[5] # /api/zoo/v1/animals/10, get 10
    animal = ANIMALS[parse(Int, animalId)]
    return HTTP.Response(200, JSON2.write(animal))
end

function updateAnimal(req::HTTP.Request)
    animal = JSON2.read(IOBuffer(HTTP.payload(req)), Animal)
    ANIMALS[animal.id] = animal
    return HTTP.Response(200, JSON2.write(animal))
end

function deleteAnimal(req::HTTP.Request)
    animalId = HTTP.URIs.splitpath(req.target)[5] # /api/zoo/v1/animals/10, get 10
    delete!(ANIMALS, parse(Int, animal.id))
    return HTTP.Response(200)
end

# define REST endpoints to dispatch to "service" functions
const ANIMAL_ROUTER = HTTP.Router()
HTTP.@register(ANIMAL_ROUTER, "POST", "/api/zoo/v1/animals", createAnimal)
# note the use of `*` to capture the path segment "variable" animal id
HTTP.@register(ANIMAL_ROUTER, "GET", "/api/zoo/v1/animals/*", getAnimal)
HTTP.@register(ANIMAL_ROUTER, "PUT", "/api/zoo/v1/animals", updateAnimal)
HTTP.@register(ANIMAL_ROUTER, "DELETE", "/api/zoo/v1/animals/*", deleteAnimal)

HTTP.serve(ANIMAL_ROUTER, ip"127.0.0.1", 8080)


text:
## LocalRule
export LocalRule

"""
    LocalRule()

Transition rule for Metropolis-Hastings sampling where at every step a random
site is switched to another random state.
"""
struct LocalRule <: MCMCRule
end

function propose_step!(σp::Union{AState,ADoubleState}, s::MetropolisSampler{LocalRule},
                       net::NeuralNetwork, c, rc)
    flipat = rand(c.rng, 1:nsites(c.hilb))
    @inbounds old_val, new_val = flipat!(c.rng, c.σp, c.hilb, flipat)
end

function propose_step!(σp::Union{AStateBatch,ADoubleStateBatch},
                       s::MetropolisSampler{LocalRule},
                       net::NeuralNetwork, c, rc)
    sites_to_flip = 1:nsites(c.hilb)

    for σ=states(σp)
        flipat = rand(c.rng, sites_to_flip)
        @inbounds old_val, new_val = flipat!(c.rng, σ, c.hilb, flipat)
    end
end


text:
module InteractiveChaos

using AbstractPlotting, Observables, MakieLayout

"""
    subscript(i::Int)
Transform `i` to a string that has `i` as a subscript.
"""
function subscript(i::Int)
    if i < 0
        "₋"*subscript(-i)
    elseif i == 1
        "₁"
    elseif i == 2
        "₂"
    elseif i == 3
        "₃"
    elseif i == 4
        "₄"
    elseif i == 5
        "₅"
    elseif i == 6
        "₆"
    elseif i == 7
        "₇"
    elseif i == 8
        "₈"
    elseif i == 9
        "₉"
    elseif i == 0
        "₀"
    else
        join(subscript.(digits(i)))
    end
end

const MARKER = Circle(Point2f0(0, 0), Float32(1)) # allows pixel size (zoom independent)
const DEFAULT_BG = RGBf0(0.99, 0.98, 1.0)
using AbstractPlotting: px
randomcolor(args...) = RGBAf0(0.9 .* (rand(), rand(), rand())..., 0.75)

# JULIADYNAMICS_COLORS = to_color.(("#7a60bb", "#202020", "#1ba5aa"))
JULIADYNAMICS_COLORS = to_color.(["#7d53e7", "#202020", "#17c7cd"])

function colors_from_map(cmap, α, N)
    N == 1 && return [RGBAf0(0, 0, 0, 1)]
    cs = [RGBAf0(c.r, c.g, c.b, α) for c in AbstractPlotting.to_colormap(cmap, N)]
end

function pushupdate!(o::Observable, v)
    push!(o[], v)
    o[] = o[]
    return o
end

function darken_color(c, f = 1.2)
    tc = to_color(c)
    return RGBAf0(c.r/f, c.g/f, c.b/f, c.alpha)
end

include("numericdata/plot_dataset.jl")
include("numericdata/trajectory_highlighter.jl")
include("chaos/orbitdiagram.jl")
include("chaos/poincareclick.jl")
include("billiards/defs_plotting.jl")
include("billiards/defs_animating.jl")
include("billiards/interactive_billiard.jl")

end