text:
using FreeTypeAbstraction: FTFont, FT_Get_Char_Index, hadvance, leftinkbound, inkwidth,
                           get_extent, ascender, descender

# Taken from JuliaPlots/AbstractPlotting.jl

mutable struct TextureAtlas
    rectangle_packer::RectanglePacker
    mapping::Dict{Any,Int} # styled glyph to index in sprite_attributes
    index::Int
    data::Matrix{Float16}
    # rectangles we rendered our glyphs into in normalized uv coordinates
    uv_rectangles::Vector{Vec4f0}
end

Base.size(atlas::TextureAtlas) = size(atlas.data)

@enum GlyphResolution High Low

const TEXTURE_RESOLUTION = Ref((2048, 2048))
const CACHE_RESOLUTION_PREFIX = Ref("high")

const HIGH_PIXELSIZE = 64
const LOW_PIXELSIZE = 32

const PIXELSIZE_IN_ATLAS = Ref(HIGH_PIXELSIZE)

function size_factor()
    return DOWN_SAMPLE_HIGH / DOWN_SAMPLE_FACTOR[]
end

function set_glyph_resolution!(res::GlyphResolution)
    if res == High
        TEXTURE_RESOLUTION[] = (2048, 2048)
        CACHE_RESOLUTION_PREFIX[] = "high"
        PIXELSIZE_IN_ATLAS[] = HIGH_PIXELSIZE
    else
        TEXTURE_RESOLUTION[] = (1024, 1024)
        CACHE_RESOLUTION_PREFIX[] = "low"
        PIXELSIZE_IN_ATLAS[] = LOW_PIXELSIZE
    end
end

function TextureAtlas(initial_size = TEXTURE_RESOLUTION[])
    return TextureAtlas(RectanglePacker(Rect2D(0, 0, initial_size...)),
                        Dict{Tuple{Char,String},Int}(), 1, zeros(Float16, initial_size...),
                        Vec4f0[])
end

function get_cache_path()
    return abspath(asset_path("fonts", ".cache",
                              "texture_atlas_$(CACHE_RESOLUTION_PREFIX[]).jls"))
end

const _default_font = FTFont[]
const _alternative_fonts = FTFont[]

function default_font()
    if isempty(_default_font)
        push!(_default_font, FTFont(asset_path("fonts", "DejaVuSans.ttf")))
    end
    return _default_font[]
end

function alternative_fonts()
    if isempty(_alternative_fonts)
        append!(_alternative_fonts,
                FTFont.(asset_path.(("fonts",), searchdir(asset_path("fonts")))))
    end
    return _alternative_fonts
end

function load_ascii_chars!(atlas)
    # for c in '\u0000':'\u00ff' #make sure all ascii is mapped linearly
    #     insert_glyph!(atlas, c, default_font())
    # end
    # for c in 'a':'z'
    #     insert_glyph!(atlas, c, default_font())
    # end
    for c in 'A':'Z'
        insert_glyph!(atlas, c, default_font())
    end
    for c in '0':'9'
        insert_glyph!(atlas, c, default_font())
    end
end

function cached_load()
    if isfile(get_cache_path())
        try
            return open(get_cache_path()) do io
                dict = Serialization.deserialize(io)
                fields = map(fieldnames(TextureAtlas)) do n
                    v = dict[n]
                    return isa(v, Vector) ? copy(v) : v # otherwise there seems to be a problem with resizing
                end
                return TextureAtlas(fields...)
            end
        catch e
            @info("You can likely ignore the following warning, if you just switched Julia versions for Makie")
            @warn(e)
            rm(get_cache_path())
        end
    end
    atlas = TextureAtlas()
    @info("Makie/AbstractPlotting is caching fonts, this may take a while. Needed only on first run!")
    load_ascii_chars!(atlas)
    to_cache(atlas) # cache it
    return atlas
end

function to_cache(atlas)
    if !ispath(dirname(get_cache_path()))
        mkpath(dirname(get_cache_path()))
    end
    open(get_cache_path(), "w") do io
        dict = Dict(map(fieldnames(typeof(atlas))) do name
                        return name => getfield(atlas, name)
                    end)
        return Serialization.serialize(io, dict)
    end
end

const global_texture_atlas = Base.RefValue{TextureAtlas}()

function get_texture_atlas()
    if isassigned(global_texture_atlas) &&
       size(global_texture_atlas[]) == TEXTURE_RESOLUTION[]
        global_texture_atlas[]
    else
        global_texture_atlas[] = cached_load() # initialize only on demand
        global_texture_atlas[]
    end
end

"""
    find_font_for_char(c::Char, font::FTFont)
Finds the best font for a character from a list of fallback fonts, that get chosen
if `font` can't represent char `c`
"""
function find_font_for_char(c::Char, font::FTFont)
    FT_Get_Char_Index(font, c) != 0 && return font
    for afont in alternative_fonts()
        if FT_Get_Char_Index(afont, c) != 0
            return afont
        end
    end
    return error("Can't represent character $(c) with any fallback font nor $(font.family_name)!")
end

function glyph_index!(atlas::TextureAtlas, c::Char, font::FTFont)
    if FT_Get_Char_Index(font, c) == 0
        for afont in alternative_fonts()
            if FT_Get_Char_Index(afont, c) != 0
                font = afont
            end
        end
    end
    return insert_glyph!(atlas, c, font)
end

function glyph_uv_width!(atlas::TextureAtlas, c::Char, font::FTFont)
    return atlas.uv_rectangles[glyph_index!(atlas, c, font)]
end

function glyph_uv_width!(c::Char)
    return glyph_uv_width!(get_texture_atlas(), c, defaultfont())
end

function glyph_boundingbox(c::Char, font::FTFont, pixelsize)
    if FT_Get_Char_Index(font, c) == 0
        for afont in alternativefonts()
            if FT_Get_Char_Index(afont, c) != 0
                font = afont
                break
            end
        end
    end
    bb, ext = FreeTypeAbstraction.metrics_bb(c, font, pixelsize)
    return bb
end

function insert_glyph!(atlas::TextureAtlas, glyph::Char, font::FTFont)
    return get!(atlas.mapping, (glyph, FreeTypeAbstraction.fontname(font))) do
        downsample = 5 # render font 5x larger, and then downsample back to desired pixelsize
        pad = 8 # padd rendered font by 6 pixel in each direction
        uv_pixel = render(atlas, glyph, font, downsample, pad)
        tex_size = Vec2f0(size(atlas.data) .- 1) # starts at 1

        idx_left_bottom = minimum(uv_pixel)# 0 based!!!
        idx_right_top = maximum(uv_pixel)

        # include padding
        left_bottom_pad = idx_left_bottom .+ pad .- 1
        # -1 for indexing offset
        right_top_pad = idx_right_top .- pad

        # transform to normalized texture coordinates
        uv_left_bottom_pad = (left_bottom_pad) ./ tex_size
        uv_right_top_pad = (right_top_pad) ./ tex_size

        uv_offset_rect = Vec4f0(uv_left_bottom_pad..., uv_right_top_pad...)
        i = atlas.index
        push!(atlas.uv_rectangles, uv_offset_rect)
        atlas.index = i + 1
        return i
    end
end

"""
    sdistancefield(img, downsample, pad)
Calculates a distance fields, that is downsampled `downsample` time,
with a padding applied of `pad`.
The padding is in units after downscaling!
"""
function sdistancefield(img, downsample, pad)
    # we pad before downsampling, so we need to have `downsample` as much padding
    pad = downsample * pad
    # padd the image
    padded_size = size(img) .+ 2pad

    # for the downsampling, we need to make sure that
    # we can divide the image size by `downsample` without reminder
    dividable_size = ceil.(Int, padded_size ./ downsample) .* downsample

    in_or_out = fill(false, dividable_size)
    # the size we fill the image up to
    wend, hend = size(img) .+ pad
    in_or_out[pad+1:wend, pad+1:hend] .= img .> (0.5 * 255)

    yres, xres = dividable_size .÷ downsample
    # divide by downsample to normalize distances!
    return Float16.(sdf(in_or_out, xres, yres) ./ downsample)
end

function render(atlas::TextureAtlas, glyph::Char, font, downsample = 5, pad = 6)
    #select_font_face(cc, font)
    if glyph == '\n' # don't render  newline
        glyph = ' '
    end
    # the target pixel size of our distance field
    pixelsize = PIXELSIZE_IN_ATLAS[]
    # we render the font `downsample` sizes times bigger
    bitmap, extent = renderface(font, glyph, pixelsize * downsample)
    # Our downsampeld & padded distancefield
    sd = sdistancefield(bitmap, downsample, pad)
    rect = Area(0, 0, size(sd)...)
    uv = push!(atlas.rectangle_packer, rect) # find out where to place the rectangle
    uv == nothing &&
        error("texture atlas is too small. Resizing not implemented yet. Please file an issue at Makie if you encounter this") #TODO resize surface
    # write distancefield into texture
    atlas.data[uv.area] = sd
    # return the area we rendered into!
    return uv.area
end

function layout_text(string::AbstractString, textsize::Union{AbstractVector,Number}, font,
                     align, justification = 0.5, lineheight = 1.0)
    offset_vec = alignment2num.(align)
    ft_font = font
    rscale = textsize

    atlas = get_texture_atlas()
    pos = zero(Vec2f0)

    glyphpos = glyph_positions(string, font, textsize, offset_vec[1], offset_vec[2],
                               lineheight, justification)

    positions = Vec2f0[]
    for (i, group) in enumerate(glyphpos)
        for gp in group
            # rotate around the alignment point (this is now at [0, 0, 0])
            push!(positions, pos .+ gp) # TODO why division by 4 necessary?
        end
        # between groups, push a random point for newline, it doesn't matter
        # what it is
        if i < length(glyphpos)
            push!(positions, zero(Vec2f0))
        end
    end

    return positions
end

function glyph_positions(str::AbstractString, font, fontscale_px, halign, valign,
                         lineheight_factor, justification)
    isempty(str) && return Vec2f0[]

    linebreak_indices = (i for (i, c) in enumerate(str) if c == '\n')

    groupstarts = [1; linebreak_indices .+ 1]
    groupstops = [linebreak_indices .- 1; length(str)]

    char_groups = map(groupstarts, groupstops) do start, stop
        return str[start:stop]
    end

    extents = map(char_groups) do group
        # TODO: scale as SVector not Number
        return [get_extent(font, char) .* SVector(fontscale_px, fontscale_px)
                for char in group]
    end

    # add or subtract kernings?
    xs = map(extents) do extgroup
        return cumsum([isempty(extgroup) ? 0.0 : -leftinkbound(extgroup[1]);
                       hadvance.(extgroup[1:end-1])])
    end

    # each linewidth is the last origin plus inkwidth
    linewidths = last.(xs) .+
                 [isempty(extgroup) ? 0.0 : inkwidth(extgroup[end]) for extgroup in extents]
    maxwidth = maximum(linewidths)

    width_differences = maxwidth .- linewidths

    xs_justified = map(xs, width_differences) do xsgroup, wd
        return xsgroup .+ wd * justification
    end

    # make lineheight a multiple of the largest lineheight in each line
    lineheights = map(char_groups) do group
        # guard from empty reduction
        isempty(group) && return 0.0f0
        return Float32(font.height / font.units_per_EM * lineheight_factor * fontscale_px)
    end

    # how to define line height relative to font size?
    ys = cumsum([0; -lineheights[2:end]])

    # x alignment
    xs_aligned = [xsgroup .- halign * maxwidth for xsgroup in xs_justified]

    # y alignment
    # first_max_ascent = maximum(hbearing_ori_to_top, extents[1])
    # last_max_descent = maximum(x -> inkheight(x) - hbearing_ori_to_top(x), extents[end])

    first_line_ascender = ascender(font) * fontscale_px

    last_line_descender = descender(font) * fontscale_px

    overall_height = first_line_ascender - ys[end] - last_line_descender

    ys_aligned = ys .- first_line_ascender .+ (1 - valign) .* overall_height

    # we are still operating in freetype units, let's convert to the chosen scale by dividing with 64
    return [Vec2f0.(xsgroup, y) for (xsgroup, y) in zip(xs_aligned, ys_aligned)]
end

function text_bb(str, font, size)
    positions = layout_text(str, Point2f0(0), size, font, Vec2f0(0),
                            Quaternionf0(0, 0, 0, 1), Mat4f0(I), 0.5, 1.0)

    scale = widths.(first.(FreeTypeAbstraction.metrics_bb.(collect(str), font, size)))
    return union(FRect3D(positions), FRect3D(positions .+ to_ndim.(Point3f0, scale, 0)))
end

make_iter(x) = repeated(x)
make_iter(x::AbstractArray) = x

function get_iter(defaultfunc, dictlike, key)
    return make_iter(get(defaultfunc, dictlike, key))
end

function alignment2num(x::Symbol)
    (x == :center) && return 0.5f0
    (x in (:left, :bottom)) && return 0.0f0
    (x in (:right, :top)) && return 1.0f0
    return 0.0f0 # 0 default, or better to error?
end

@component struct FontStorage
    atlas       :: TextureAtlas
    storage_fbo :: GLA.FrameBuffer #All Glyphs should be stored in the first color attachment
end

function FontStorage()
    atlas = get_texture_atlas()
    fbo   = GLA.FrameBuffer(size(atlas.data), GLA.Texture(atlas.data; minfilter = :linear, magfilter = :linear, anisotropic = 16.0f0))
    return FontStorage(atlas, fbo)
end

@render_program TextProgram
function shaders(::Type{TextProgram})
    return [load_shader("sprites.geom"), load_shader("sprites.vert"),
            load_shader("distance_shape.frag")]
end

@vao TextVao

struct TextUploader <: System end
Overseer.requested_components(::TextUploader) = (Text, TextVao, TextProgram, FontStorage)

function Overseer.prepare(::TextUploader, dio::Diorama)
    e = Entity(dio[DioEntity], 1)
    if isempty(dio[TextProgram])
        dio[e] = TextProgram(Program(text_shaders()))
    end
    if isempty(dio[FontStorage])
        # dio[e] = FontStorage()
    end
end

function Overseer.update(::TextUploader, m::AbstractLedger)
    text = m[Text]
    vao = m[TextVao]
    prog = singleton(m, TextProgram)
    ucolor = m[UniformColor]
    spatial = m[Spatial]
    vis = m[Visible]
    atlas = get_texture_atlas()
    for e in @entities_in(text && spatial && ucolor)
        t = text[e]
        offset_width, uv_texture_bbox = to_gl_text(t, atlas)
        nsprites = length(t.str)
        if !(e ∈ vao) || nsprites > length(vao[e])
            vao[e] = TextVao(VertexArray(generate_buffers(prog.program, GLint(0);
                                                          color           = fill(ucolor[e].color, nsprites),
                                                          rotation        = fill(Vec4f0(0), nsprites),
                                                          offset_width    = offset_width,
                                                          uv_texture_bbox = uv_texture_bbox),
                                         GL_POINTS))
        else
            GLA.upload!(vao[e]; color = fill(ucolor[e].color, nsprites),
                        rotation = fill(Vec4f0(0), nsprites), offset_width = offset_width,
                        uv_texture_bbox = uv_texture_bbox)
        end
        if !(e in vis)
            vis[e] = Visible()
        end
    end
end

to_gl_text(t::Text, storage::FontStorage) = to_gl_text(t, storage.atlas)
to_gl_text(t::Text, atlas) = to_gl_text(t.str, t.font_size, t.font, t.align, atlas)

function to_gl_text(string::AbstractString, textsize, font, align::Tuple{Symbol,Symbol},
                    atlas)
    # rscale          = Float32(textsize)
    chars  = Vector{Char}(string)
    scale  = Vec2f0[]
    offset = Vec2f0[]
    for c in string
        glyph_bb, extent = FreeTypeAbstraction.metrics_bb(c, font, textsize)
        push!(scale, widths(glyph_bb))
        push!(offset, minimum(glyph_bb))
    end
    # scale           = Vec2f0.(glyph_scale!.(Ref(atlas), chars, (font,), rscale))
    positions2d = layout_text(string, textsize, font, align)

    # aoffset         = align_offset(Point2f0(0), positions2d[end], atlas, rscale, font, align)

    uv_offset_width = glyph_uv_width!.(Ref(atlas), chars, (font,))
    out_uv_offset_width = Vec4f0[]
    for uv in uv_offset_width
        push!(out_uv_offset_width, Vec4f0(uv[1], uv[2], uv[3], uv[4]))
    end
    out_pos_scale = Vec4f0[]
    for (p, sc) in zip(positions2d .+ offset, scale)
        push!(out_pos_scale, Vec4f0(p[1], p[2], sc[1], sc[2]))
    end
    return out_pos_scale, out_uv_offset_width
end

struct TextRenderer <: AbstractRenderSystem end

function Overseer.requested_components(::TextRenderer)
    return (Spatial, UniformColor, Camera3D, TextVao, Text, TextProgram, IOTarget,
            FontStorage)
end

function Overseer.update(::TextRenderer, m::AbstractLedger)
    spat           = m[Spatial]
    prog           = singleton(m, TextProgram)
    camera         = singleton(m, Camera3D)
    vao            = m[TextVao]
    iofbo          = singleton(m, Canvas)
    persp_mat      = camera.projview
    projection_mat = camera.proj
    text           = m[Text]
    vis            = m[Visible]
    wh             = size(iofbo)

    glDisable(GL_DEPTH_TEST)
    glDepthFunc(GL_ALWAYS)
    glDisableCullFace()

    glEnablei(GL_BLEND, 0)
    glDisablei(GL_BLEND, 1)
    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)

    bind(iofbo)
    draw(iofbo)

    bind(prog)
    gluniform(prog, :resolution, Vec2f0(wh))
    gluniform(prog, :projview, persp_mat)
    gluniform(prog, :projection, projection_mat)

    #Absolutely no clue why this needs to be here?...
    if isempty(m[FontStorage])
        fontstorage = FontStorage()
        m[Entity(m[DioEntity], 1)] = fontstorage
    else
        fontstorage = singleton(m, FontStorage)
    end
    glyph_fbo = fontstorage.storage_fbo

    # Fragment uniforms
    gluniform(prog, :distancefield, 0, color_attachment(glyph_fbo, 1))
    gluniform(prog, :shape, 3)
    #TODO make this changeable
    # gluniform(prog, :stroke_color, zero(Vec4f0))
    # gluniform(prog, :glow_color, zero(Vec4f0))
    gluniform(prog, :stroke_width, 0.0f0)
    gluniform(prog, :glow_width, 0.0f0)
    gluniform(prog, :billboard, true)
    gluniform(prog, :scale_primitive, true)
    for e in @entities_in(vao && spat && text)
        e_vao, e_spat, e_text = vao[e], spat[e], text[e]
        if vis[e].visible
            gluniform(prog, :model, m[ModelMat][e].modelmat)
            gluniform(prog, :origin, e_text.offset)
            bind(e_vao)
            draw(e_vao)
        end
    end
    glDisablei(GL_BLEND, 0)
    return glEnablei(GL_BLEND, 1)
end


text:
mutable struct XAxis <: Axis
    ticklabels::Vector{Label{Vector{String}}}
    ticksublabels::Vector{Label{Vector{String}}}
    frame::Line
end

XAxis( ; ticklabels, ticksublabels, frame) = XAxis(ticklabels, ticksublabels, frame)

ext: jl


text:
module QuasiArrays
using Base, LinearAlgebra, LazyArrays, ArrayLayouts, DomainSets, FillArrays, StaticArrays
import Base: getindex, size, axes, axes1, length, ==, isequal, iterate, CartesianIndices, LinearIndices,
                Indices, IndexStyle, getindex, setindex!, parent, vec, convert, similar, copy, copyto!, zero,
                map, eachindex, eltype, first, last, firstindex, lastindex, in, reshape, permutedims, all,
                isreal, iszero, isempty, empty, isapprox, fill!, getproperty, showarg
import Base: @_inline_meta, DimOrInd, OneTo, @_propagate_inbounds_meta, @_noinline_meta,
                DimsInteger, error_if_canonical_getindex, @propagate_inbounds, _return_type,
                _maybetail, tail, _getindex, _maybe_reshape, index_ndims, _unsafe_getindex,
                index_shape, to_shape, @nloops, @ncall, unalias, _unaliascopy,
                to_index, to_indices, _to_subscript_indices, _splatmap, dataids,
                compute_stride1, compute_offset1, fill_to_length
import Base: Slice, IdentityUnitRange, ScalarIndex, RangeIndex, view, viewindexing, mightalias, ensure_indexable, index_dimsum,
                check_parent_index_match, reindex, _isdisjoint, _unsafe_ind2sub,
                _ind2sub, _sub2ind, _ind2sub_recurse, _lookup, SubArray,
                parentindices, reverse, ndims, checkbounds, uncolon,
                maybeview, unsafe_view, checkindex, checkbounds_indices,
                throw_boundserror, rdims, replace_in_print_matrix, show, summary,
                hcat, vcat, hvcat
import Base: *, /, \, +, -, ^, inv
import Base: exp, log, sqrt,
          cos, sin, tan, csc, sec, cot,
          cosh, sinh, tanh, csch, sech, coth,
          acos, asin, atan, acsc, asec, acot,
          acosh, asinh, atanh, acsch, asech, acoth
import Base: Array, Matrix, Vector
import Base: union, intersect, sort, sort!
import Base: conj, real, imag
# reducedim.jl imports
import Base: prod, sum, cumsum, diff, add_sum, mul_prod, mapreduce, max, min, count, _count, any, _any, all, _all, _sum, _prod, _mapreduce, reduced_index, check_reducedims, mapfoldl_impl
import Base: BitInteger, IEEEFloat, uniontypes, _InitialValue, safe_tail, reducedim1, _simple_count

import Base: ones, zeros, one, zero, fill

import Base.Broadcast: materialize, materialize!, BroadcastStyle, AbstractArrayStyle, Style, broadcasted, Broadcasted, Unknown,
                        newindex, broadcastable, preprocess, _eachindex, _broadcast_getindex, broadcast_shape,
                        DefaultArrayStyle, axistype, throwdm, instantiate, combine_eltypes, eltypes, combine_styles

import LinearAlgebra: transpose, adjoint, checkeltype_adjoint, checkeltype_transpose, Diagonal,
                        AbstractTriangular, pinv, inv, promote_leaf_eltypes, power_by_squaring,
                        integerpow, schurpow, tr, factorize, copy_oftype, rank, dot

import ArrayLayouts: indextype, concretize, fillzeros, OnesLayout, AbstractFillLayout, FillLayout, ZerosLayout, diagonallayout, diagonaldata, diagonal
import LazyArrays: MemoryLayout, UnknownLayout, Mul, ApplyLayout, BroadcastLayout,
                    InvOrPInv, ApplyStyle, AbstractLazyLayout, LazyLayout,
                    MulStyle, MulAddStyle, LazyArrayApplyStyle, combine_mul_styles, DefaultArrayApplyStyle,
                    Applied, flatten, _flatten, arguments, _mat_mul_arguments, _vec_mul_arguments, _mul_arguments,
                    rowsupport, colsupport, tuple_type_memorylayouts, applylayout, broadcastlayout,
                    LdivStyle, most, InvLayout, PInvLayout, sub_materialize, lazymaterialize,
                    _mul, rowsupport, DiagonalLayout, adjointlayout, transposelayout, conjlayout,
                    sublayout, call, LazyArrayStyle, layout_getindex, _broadcast2broadcastarray, _applyarray_summary, _broadcastarray_summary,
                    _broadcasted_mul, simplifiable, simplify

import Base.IteratorsMD

export AbstractQuasiArray, AbstractQuasiMatrix, AbstractQuasiVector, materialize,
       QuasiArray, QuasiMatrix, QuasiVector, QuasiDiagonal, Inclusion,
       QuasiAdjoint, QuasiTranspose, ApplyQuasiArray, ApplyQuasiMatrix, ApplyQuasiVector,
       BroadcastQuasiArray, BroadcastQuasiMatrix, BroadcastQuasiVector, indextype,
       QuasiKron, quasikron, UnionVcat, SubQuasiArray, QuasiFill, QuasiZeros, QuasiOnes, QuasiEye

import Base.Broadcast: broadcast_preserving_zero_d

import Base: oneto

abstract type AbstractQuasiArray{T,N} end
AbstractQuasiVector{T} = AbstractQuasiArray{T,1}
AbstractQuasiMatrix{T} = AbstractQuasiArray{T,2}
AbstractQuasiVecOrMat{T} = Union{AbstractQuasiVector{T}, AbstractQuasiMatrix{T}}
const AbstractQuasiOrVector{T} = Union{AbstractVector{T},AbstractQuasiVector{T}}
const AbstractQuasiOrMatrix{T} = Union{AbstractMatrix{T},AbstractQuasiMatrix{T}}
const AbstractQuasiOrArray{T} = Union{AbstractArray{T},AbstractQuasiArray{T}}

cardinality(d) = length(d)
cardinality(d::UnionDomain) = +(cardinality.(d.domains)...)

size(A::AbstractQuasiArray) = map(cardinality, axes(A))
axes(A::AbstractQuasiArray) = error("Override axes for $(typeof(A))")

include("indices.jl")
include("abstractquasiarray.jl")
include("multidimensional.jl")
include("subquasiarray.jl")
include("quasireshapedarray.jl")
include("quasipermutedims.jl")
include("quasibroadcast.jl")
include("abstractquasiarraymath.jl")
include("quasireducedim.jl")

include("quasiarray.jl")
include("quasiarraymath.jl")

include("lazyquasiarrays.jl")
include("quasiconcat.jl")

copy_oftype(A::AbstractQuasiArray{T}, ::Type{T}) where {T} = copy(A)
copy_oftype(A::AbstractQuasiArray{T,N}, ::Type{S}) where {T,N,S} = convert(AbstractQuasiArray{S,N}, A)

include("matmul.jl")
include("inv.jl")
include("quasiadjtrans.jl")
include("quasidiagonal.jl")
include("quasifill.jl")
include("dense.jl")

include("quasikron.jl")

promote_leaf_eltypes(x::AbstractQuasiArray{T}) where {T<:Number} = T
promote_leaf_eltypes(x::AbstractQuasiArray) = eltype(eltype(x)) #mapreduce(promote_leaf_eltypes, promote_type, x; init=Bool)

function isapprox(x::AbstractQuasiArray, y::AbstractQuasiArray;
    atol::Real=0,
    rtol::Real=Base.rtoldefault(promote_leaf_eltypes(x),promote_leaf_eltypes(y),atol),
    nans::Bool=false, norm::Function=norm)
    d = norm(x - y)
    if isfinite(d)
        return d <= max(atol, rtol*max(norm(x), norm(y)))
    else
        # Fall back to a component-wise approximate comparison
        return all(ab -> isapprox(ab[1], ab[2]; rtol=rtol, atol=atol, nans=nans), zip(x, y))
    end
end

end


text:
# Load packages

using FSM
using PyPlot
using Distributions
using ProgressMeter
using DataAssim

# Load data

function load_data()

  # Load data

  metdata = readdlm(joinpath(dirname(@__FILE__), "..", "data", "met_CdP_all.txt"))
  valdata = readdlm(joinpath(dirname(@__FILE__), "..", "data", "obs_CdP_all.txt"))

  prec = metdata[:, 7] + metdata[:, 8]

  metdata = hcat(metdata, prec)

  valdata[valdata .== -99] = NaN

  dates_met = map(DateTime, metdata[:,1], metdata[:,2], metdata[:,3], metdata[:,4])
  dates_val = map(DateTime, valdata[:,1], valdata[:,2], valdata[:,3], valdata[:,4])

  return dates_met, metdata, dates_val, valdata

end

# Load data

dates_met, metdata, dates_val, valdata = load_data()

# Settings

npart = 100
thres_prec = 273.6890
m_prec = 0.3051
p_corr = 1.0
timestep = 1.0

ntimes = size(metdata,1)

# Initial state

md = [FsmType(1, 1, 1, 1, 0) for i in 1:npart]

q_vec = [q_noise(randn(), randn(), randn(), randn(), randn(), randn()) for i in 1:npart]

# Allocate input array

ip = Array(FsmInput, npart)

# Allocate output array

od = zeros(ntimes, npart)

hs_sim = zeros(npart)

res = zeros(ntimes, 3)

# Initilize particles

wk = ones(npart) / npart

# Loop over time

@showprogress 1 "Computing..." for itime = 1:ntimes

  # Run the model for all particles

  for ipart = 1:npart

    # Handle input data

    year  = metdata[itime, 1]
    month = metdata[itime, 2]
    day   = metdata[itime, 3]
    hour  = metdata[itime, 4]
    SW    = metdata[itime, 5]
    LW    = metdata[itime, 6]
    P     = metdata[itime, 7] + metdata[itime, 8]
    Ta    = metdata[itime, 9]
    RH    = metdata[itime, 10]
    Ua    = metdata[itime, 11]
    Ps    = metdata[itime, 12]

    SW, LW, P, Ta, RH, Ua = perturb_input(q_vec, ipart, timestep, SW, LW, P, Ta, RH, Ua)

    Sf = precip_solid(P, Ta, thres_prec, p_corr, m_prec)

    Rf = precip_liquid(P, Ta, thres_prec, p_corr, m_prec)

    # Run the model

    ip[ipart] = FsmInput(year, month, day, hour, SW, LW, Sf, Rf, Ta, RH, Ua, Ps)

    hs_sim[ipart] = run_fsm(md[ipart], ip[ipart])

  end

  # Find observation

  iobs = find(dates_val .== dates_met[itime])

  hs_obs = valdata[iobs, 5]

  # Run particle filter

  if !isempty(hs_obs) && !isnan(hs_obs[1])

    for ipart = 1:npart
      wk[ipart] = pdf(Normal(hs_obs[1], max(0.1 * hs_obs[1], 0.05)), hs_sim[ipart]) * wk[ipart]
    end

    if sum(wk) > 0.0
      wk = wk / sum(wk)
    else
      wk = ones(npart) / npart
    end

    # Perform resampling

    Neff = 1 / sum(wk.^2)

    if round(Int64, Neff) < round(Int64, npart * 0.8)

      #println("Resampled at step: $itime")

      indx = resample(wk)

      md  = [deepcopy(md[i]) for i in indx]

      q_vec = [deepcopy(q_vec[i]) for i in indx]

      wk = ones(npart) / npart

    end

  end

  # Store results

  res[itime, 1] = sum(wk .* hs_sim)
  res[itime, 2] = minimum(hs_sim)
  res[itime, 3] = maximum(hs_sim)

end

# Plot results

fig = plt[:figure](figsize = (12,7))

plt[:style][:use]("ggplot")

plt[:plot](dates_val, valdata[:,5], linewidth = 1.2, color = "k", label = "Observed", zorder = 1)
plt[:fill_between](dates_met, res[:, 3], res[:, 2], facecolor = "r", edgecolor = "r", label = "Simulated", alpha = 0.55, zorder = 2)
plt[:legend]()

#using plots
#plot(dates_met, res[:, 3], fillrange = res[:, 2], color = "red", fillalpha = 0.5, linealpha = 0.5, label = "sim")
#plot(dates_met, [res[:, 3], res[:, 2]], fillrange = res[:, 2], color = "red", linealpha = 0.5, fillalpha = 0.5, label = "", grid = false)
#plot!(dates_val, valdata[:,5], linewidth = 1.2, color = "k")
# plot([x;x[end:-1:1]],[y1;y2[end:-1:1]], color = :red, linealpha = 0, fillalpha = 0.5, seriestype = :shape, label = "a filled area")


text:
lhs(eq::Assignment) = eq.lhs
rhs(eq::Assignment) = value(eq.rhs)

op(eq::Assignment) = operation(rhs(eq))
args(eq::Assignment) = Num.(arguments(rhs(eq)))

name(var) = value(var).name
bar(var) = Variable(Symbol(var, '̄'))  # the character is the overbar symbol (on top of first `'`)

# bar(op) = Variable()

# "Construct adjoint of a given Assignment with accumulation"
# function adj(eq::Assignment)
#     vars = args(eq)
#     adjoints = Adjoint(op(eq), vars)

#     bar_lhs = bar(lhs(eq))

#     eqns = map(zip(vars, adjoints)) do (var, adj)
#         barred = bar(var)

#         Assignment(barred, barred + bar_lhs * adj)
#     end

#     return eqns
# end

"Make a new variable by adding a number to an old variable, e.g. `x` -> `x3`"
function add_number(var, num)
    return Variable(var, num)   # Variable(Symbol(var, "_", num))
end

function numbered_variable!(num_times_used, var, increment=false)
    is_new = false

    if !haskey(num_times_used, var)
        is_new = true

        num = 0
        num_times_used[var] = 0
    else

        if increment
            num_times_used[var] += 1
        end

        num = num_times_used[var]
    end

    return is_new, add_number(var, num)
end

function adj(eq::Assignment)
    vars = value.(args(eq))

    bar_lhs = bar(lhs(eq))

    adjoints = adj(op(eq), bar_lhs, vars...)

    eqns = map(zip(vars, adjoints)) do (var, adj)
        @show var, adj, typeof(var)
        #if var isa Sym{Tangent{Real}}  # for linearization pass
            barred = bar(var)

            Assignment(barred, barred + adj)
        #end
    end

    return eqns
end

# modifies the Set assigned of variables which have already been assigned
function simple_adj(eq::Assignment, assigned)

    vars = ReversePropagation.args(eq)
    bar_lhs = bar(lhs(eq))

    adjoints = adj(op(eq), bar_lhs, vars...)

    eqns = Assignment[]

    for (var, adjoint) in zip(vars, adjoints)

        !isa(var, Num) && continue

        barred = bar(var)

        if any(x -> isequal(x, barred), assigned)
            push!(eqns, Assignment(barred, barred + adjoint))

        else
            push!(eqns, Assignment(barred, adjoint))
            push!(assigned, barred)
        end

    end

    return eqns
end

function adj(num_times_used, eq::Assignment)

    # @show num_times_used
    # @show eq

    vars = args(eq)
    # adjoints = Adjoint(op(eq), vars)

    # is_new, bar_lhs = numbered_variable!(num_times_used, bar(lhs(eq)))

    bar_lhs = bar(lhs(eq))
    is_new, numbered_bar_lhs = numbered_variable!(num_times_used, bar_lhs)

    adjoints = adj(op(eq), numbered_bar_lhs, vars...)

    eqns = Assignment[]

    # @show bar_lhs
    for (var, adj) in zip(vars, adjoints)

        # @show typeof(var)
        !isa(var, Num) && continue

        barred = bar(var)

        is_new, numbered = numbered_variable!(num_times_used, barred)

        if is_new # first use
            # Assignment(numbered, numbered_bar_lhs * adj)
            push!(eqns, Assignment(numbered, adj))

        else  # re-use so generate new variable
            new_num, barred_new = numbered_variable!(num_times_used, barred, true)
            # Assignment(barred_new, numbered + numbered_bar_lhs * adj)
            push!(eqns, Assignment(barred_new, numbered + adj))

        end
    end

    return eqns
end

function reverse_pass(vars, code, final)
    num_times_used = Dict{Num, Int}()

    num, final_bar = numbered_variable!(num_times_used, bar(final))
    reverse_code = [Assignment(final_bar, 1)]

    for eq in reverse(code)
        # @show eq
        append!(reverse_code, adj(num_times_used, eq))
    end

    final_vars = [numbered_variable!(num_times_used, bar(var))[2] for var in vars]
    return reverse_code, final_vars
end

function simple_reverse_pass(vars, forward_code)

    assigned = Set()  # which variables have already been assigned

    reverse_code = reduce(vcat, simple_adj.(reverse(forward_code), Ref(assigned)))

    final_vars = bar.(vars)

    return reverse_code, final_vars, assigned

end

"""
Return code for forward and reverse pass, as a vector of Assignments.
`final` is the output variable from the forward pass.
`gradient_vars` are the output variables from the reverse pass.
"""
function gradient_code(ex, vars)
    forward_code, final = cse_equations(ex)
    # reverse_code, gradient_vars = reverse_pass(vars, forward_code, final)

    reverse_code, gradient_vars, assigned = simple_reverse_pass(vars, forward_code)

    initialization_code = [Assignment(bar(final), 1)]  # need typed 1 and 0?

    unassigned = setdiff(gradient_vars, assigned)
    append!(initialization_code, [Assignment(var, 0) for var in unassigned])

    code = forward_code ∪ initialization_code ∪ reverse_code

    return code, final, gradient_vars
end

make_tuple(args) = Expr(:tuple, args...)
make_tuple(s::Symbol) = make_tuple([s])

# toexpr(ex::Assignment) = toexpr(Equation(ex.lhs, ex.rhs))

function gradient_expr(vars, ex)
    symbolic_code, final, gradient_vars = gradient_code(vars, ex)

    code = Expr(:block, toexpr.(symbolic_code)...)

    return code, final, gradient_vars
end

function gradient(ex::Num, vars)

    code, final_var, gradient_vars = gradient_expr(ex, vars)

    input_vars = toexpr(Symbolics.MakeTuple(vars))
    final = toexpr(final_var)
    gradient = toexpr(Symbolics.MakeTuple(gradient_vars))

    full_code = quote
        ($input_vars, ) -> begin
            $code
            return $(final), $(gradient)
        end
    end

    return eval(full_code)

end

gradient(f, vars) = gradient(f(vars), vars)


text:
""""
__Predicates__ are functions that ask "yes or no" questions of their argument[s].
You can ask of a number "Is this zero?" or "Is this one?" and these predicates
(`iszero`, `isone`) will work as expected with almost all numerical types.
The built-in numerical types let you query finiteness (`isfinite`, `isinf`).
These are the predicates made available for use with DoubleFloats:

> iszero, isnonzero, isone                 #  value == 0, value != 0, value == 1
  ispositive, isnegative,                  #  value >  0, value <  0
  isnonnegative, isnonpositive,            #  value >= 0, value <= 0
  isfinite, isinf,                         #  abs(value) != Inf, abs(value) == Inf
  isposinf, isneginf,                      #  value == Inf, value == -Inf
  isnan,                                   #  value is not a number (eg 0/0)
  issubnormal,                             #  value contains a subnormal part
  isnormal,                                #  value is finite and not subnormal
  isinteger, isfractional                  #  value == round(value)
  iseven, isodd,                           #  isinteger(value/2.0), !isinteger(value/2.0)
""" predicates

"""
    isnonzero(x)

Return `!iszero(x)`
"""
isnonzero(x::T) where {T<:AbstractFloat} = !iszero(x)

"""
    ispositive(x)

Returns `true` if `!isnegative(x)` and `isnonzero(x)`.
"""
ispositive(x::T) where {T<:AbstractFloat} = !isnegative(x) && isnonzero(x)

"""
    isnonnegative(x)

Returns `true` if `!isnegative(x)`.
"""
isnonnegative(x::T) where {T<:AbstractFloat} = !isnegative(x)
"""
    isnonpositive(x)

Returns `true` if `isnegative(x)` or `iszero(x)`.
"""
isnonpositive(x::T) where {T<:AbstractFloat} = isnegative(x) || iszero(x)

"""
    isfractional(x)

Returns `true` if `abs(x) < one(x)`.
"""
isfractional(x::T) where {T<:AbstractFloat} = abs(x) < one(T)

iszero(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    iszero(HI(x)) # && iszero(LO(x))

isnonzero(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    !iszero(HI(x)) # || !iszero(LO(x))

isone(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    isone(HI(x)) && iszero(LO(x))

ispositive(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    !signbit(HI(x)) && !iszero(HI(x))

"""
    isnegative(x)

Return `true` if `signbit(HI(x))` is `true`.
"""
isnegative(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    signbit(HI(x))

isnonnegative(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    !signbit(HI(x))

isnonpositive(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    signbit(HI(x)) || iszero(HI(x))

isfinite(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    isfinite(HI(x))

isinf(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    isinf(HI(x))

"""
    isposinf(x)

Tests whether a number positive and infinite.
"""
isposinf(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    isinf(HI(x))

"""
    isneginf(x)

Tests whether a number is negative and infinite.
"""
isneginf(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    isinf(HI(x)) && signbit(HI(x))

isnan(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    isnan(HI(x))

issubnormal(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    issubnormal(LO(x)) || issubnormal(HI(x))

"""
    isnormal(x)

Tests whether a floating point number is normal.
"""
isnormal(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    isfinite(HI(x)) && !iszero(x) && (abs(x) >= floatmin(T))

@inline isinteger(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    isinteger(HI(x)) && isinteger(LO(x))

isfractional(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    !isinteger(LO(x)) || !isinteger(HI(x))

isodd(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    if isinteger(x)
       (iszero(LO(x)) && isodd(HI(x))) || isodd(LO(x))
    else
       false
    end

iseven(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    if isinteger(x)
       (iszero(LO(x)) && iseven(HI(x))) || iseven(LO(x))
    else
       false
    end

iseven(x::T) where {T<:IEEEFloat} = isinteger(x) && iseven(BigInt(x))
isodd(x::T) where {T<:IEEEFloat} = isinteger(x) && isodd(BigInt(x))


text:
include("../src/valvecalculations.jl")

@testset "Dome pressures" begin
    #testing with R-value of zero to ignore the adjustment to PTRO
    @test domepressure_downhole(600-14.65, 0, 180) ≈ 750 atol = 5 #function takes psig but test values are psia
    @test domepressure_downhole(727-14.65, 0, 120) ≈ 820 atol = 5
end

@testset "Thornhill-Craver" begin
    @test ThornhillCraver_gaspassage_simplified(900, 1100, 140, 16) ≈ 1127 atol=1 #using psig inputs
    @test ThornhillCraver_gaspassage_simplified(1200, 1100, 140, 16) == 0

    seats = [7*4, 32, 5*8, 3*16]
    Hernandez_results = [3.15, 4.11, 6.36, 9.39] .* 1000
    results = [ThornhillCraver_gaspassage(420 - pressure_atmospheric, 850 - pressure_atmospheric, 150, s, 0.7) for s in seats]
    @test all(abs.(Hernandez_results .- results) .<= Hernandez_results .* 0.02) #2% tolerance to account for changing C_ds that aren't easily available (test cases use Winkler C_ds)

    @test ThornhillCraver_gaspassage(1200, 1100, 140, 16, 0.7) == 0
end

@testset "Valve table" begin

#EHU 256H example using Weatherford method

MDs = [0,1813, 2375, 2885, 3395]
TVDs = [0,1800, 2350, 2850, 3350]
incs = [0,0,0,0,0]
id = 2.441

well = Wellbore(MDs, incs, TVDs, id)
valves = GasliftValves([1813,2375,2885,3395], [1005,990,975,960], [0.073,0.073,0.073,0.073], [16,16,16,16])

tubing_pressures = [150,837,850,840,831]
casing_pressures = 1070 .+ [0,53,70,85,100]
temps = [135,145,148,151,153]

vdata, inj_depth = valve_calcs(valves = valves, well = well, sg_gas = 0.72, tubing_pressures = tubing_pressures, casing_pressures = casing_pressures, tubing_temps = temps, casing_temps = temps)

valve_table(vdata, inj_depth) #implicit test

results = vdata[1:4, [5,13,12,4]] #PSC, PVC, PVO, PSO

expected_results =
[1050. 1103 1124 1071;
 1023 1092 1111 1042;
 996  1080 1099 1015;
 968  1068 1087 987]

@test all(abs.(expected_results .- results) .< (expected_results .* 0.01)) #1% tolerance due to using TCFs versus PVT-based dome correction, as well as rounding errors

end #testset for valve table


text:
md"""
# Search Spaces

The search space is a set of possible architectures which the search policy may use to create initial candidates or to extend existing candidates.
Search spaces are constructed from simple components which can be combined in multiple ways, giving a lot of flexibility.

### Simple Parameter Spaces

The most fundamental building block of any search space if the [`ParSpace`](@ref):
"""
@testset "ParSpace example" begin #src
Random.seed!(NaiveGAflux.rng_default, 1) #src
ps1d = ParSpace([2,4,6,10])

# Draw from the search space.
@test ps1d() == 6
@test ps1d() == 10

# Possible to supply a random number generator.
@test ps1d(MersenneTwister(0)) == 4

# ParSpaces can be of any dimension and type.
ps2d = ParSpace(["1","2","3"], ["4","5","6","7"])

@test typeof(ps1d) == ParSpace{1, Int}
@test typeof(ps2d) == ParSpace{2, String}

@test ps2d() == ("1", "4")
end #src

md"""
### Layer Search Spaces

A `LayerSpace` is a search space of configurations of a single Flux layer. It is worth noting that NaiveGAflux separates search spaces for creation of new layers from
modification of existing layers. A layer drawn from a `LayerSpace` can thus be mutated into a new layer which is not part of the search space it was drawn from.

Lets look at how to create a search space for 2D convolutions:
"""

@testset "ConvSpace example" begin #src
Random.seed!(NaiveGAflux.rng_default, 1) #src
cs = ConvSpace{2}(outsizes=4:32, activations=[relu, elu, selu], kernelsizes=3:9)

inputsize = 16
convlayer = cs(inputsize)

@test string(convlayer) == "Conv((8, 3), 16 => 22, relu, pad=(4, 3, 1, 1))"
end #src

md"""
### Architecture Search Spaces

An architecture space is a search space over vertices in a computation graph, including how they are connected. Just as with [Layer Search Spaces](@ref),
the architecture search spaces are separate from modifications of graphs, meaning that an architecture drawn from a search space may be mutated into
an architecture which is not part of the searh space it was drawn from.

Lets look at an example of a search space of large image classifiers:
"""
@testset "ArchSpace example" begin #src
Random.seed!(NaiveGAflux.rng_default, 0) #src

# VertexSpace creates a MutableVertex of layers generated by the wrapped search space.
cs = VertexSpace(ConvSpace{2}(outsizes=8:256, activations=[identity, relu, elu], kernelsizes=3:5))
bs = VertexSpace(BatchNormSpace([identity, relu]))

# Block of `Conv->BatchNorm` and `BatchNorm->Conv` respectively.
# Need to make sure there is always at least one SizeAbsorb layer to make `fork` and `res` below play nice.
csbs = ArchSpaceChain(cs ,bs)
bscs = ArchSpaceChain(bs, cs)

# Randomly generates either `Conv`, `Conv->BatchNorm` or `BatchNorm->Conv`:
cblock = ArchSpace(ParSpace1D(cs, csbs, bscs))

# Generates between 1 and 5 (independent) samples from `csbs` in a sequence.
rep = RepeatArchSpace(cblock, 1:5)

# Generates between 2 and 4 parallel paths joined by concatenation (inception like-blocks) drawn from `rep`.
fork = ForkArchSpace(rep, 2:4)

# Generates a residual connection around what is generated by `rep`.
res = ResidualArchSpace(rep)

# ... and a residual fork.
resfork = ResidualArchSpace(fork)

# Pick one of the above randomly...
repforkres = ArchSpace(ParSpace1D(rep, fork, res, resfork))

# ...1 to 3 times.
blocks = RepeatArchSpace(repforkres, 1:3)

# End each block with subsamping through maxpooling.
ms = VertexSpace(PoolSpace{2}(windowsizes=2, strides=2, poolfuns=MaxPool))
reduction = ArchSpaceChain(blocks, ms)

# And lets do 2 to 4 reductions.
featureextract = RepeatArchSpace(reduction, 2:4)

# Adds 1 to 3 dense layers as outputs.
dense = VertexSpace(DenseSpace(16:512, [relu, selu]))
drep = RepeatArchSpace(dense, 0:2)
# Last layer has fixed output size (number of labels).
dout=VertexSpace(Shielded(), DenseSpace(10, identity))
output = ArchSpaceChain(drep, dout)

# Aaaand lets glue it together: Feature extracting `Conv`/`BatchNorm` layers -> global pooling -> dense layers.
archspace = ArchSpaceChain(featureextract, GlobalPoolSpace(), output)

# Input is 3 channel image.
samplemodel(invertex=conv2dinputvertex("input", 3)) = CompGraph(invertex, archspace(invertex))

# Sample one architecture from the search space.
graph1 = samplemodel()
@test nvertices(graph1) == 79

# And one more...
graph2 = samplemodel()
@test nvertices(graph2) == 128
end #src


text:
#import ..RvSpectML:   # Already in module
import ..RvSpectMLBase: predict_intrinsic_stellar_line_width, speed_of_light_mps

"""  `generate_spectrum(line_list, inst; time, rv, ssbz, snr_per_pixel, line_width, add_noise )`
Generate a spectrum using a line_list (as DataFrame with columns lambda and weight) and theoretical instrument.
Optionally, specify times, rv's, barycentric corrections, etc.
"""
function generate_spectrum(line_list::DataFrame, inst::AnyTheoreticalInstrument; time::Real = 0.0, rv::Real = 0.0, ssbz::Real = 0.0,
               snr_per_pixel::Real = 1000.0, line_width::Real = predict_intrinsic_stellar_line_width(5780), add_noise::Bool = true ) # where { LLT<:AbstractLineList }
   @assert hasproperty(line_list,:lambda)
   @assert hasproperty(line_list,:weight)
   λ = calc_λs(inst)
   flux = ones(size(λ))
   #var = ones(size(λ))
   doppler_factor = calc_doppler_factor(rv) * (1+ssbz)  # TODO: Fix problem... Is ssbz in m/s or v/c?
   for i in 1:size(line_list,1)
      #width = doppler_factor* line_list.lambda[i]*line_width*1000/speed_of_light_mps   #TODO PUT PACK just for testing
      width = doppler_factor*line_list.lambda[i]*line_width*1000/speed_of_light_mps
      flux .*= RvSpectMLBase.absorption_line.(λ, mid=line_list.lambda[i]*doppler_factor, depth=line_list.weight[i], width=width)
   end
   if add_noise
      flux .+= randn(size(flux)).*sqrt.(flux)./snr_per_pixel
      flux[flux .< 0.0] .= 0.0
   end
   var = flux ./ snr_per_pixel^2
   metadata = Dict{Symbol,Any}( :rv_true=>rv, :snr_per_pixel=>snr_per_pixel, :line_width=>line_width, :bjd=>time, :target=>"Simulation", :ssbz=>ssbz)
   metadata[:normalization] = :continuum
   if typeof(inst) <: AbstractInstrument1D

      return Spectra1DBasic(λ,flux,var,inst, metadata)
   elseif typeof(inst) <: AbstractInstrument2D
      return Spectra2DBasic(λ,flux,var,inst, metadata)
   else
      @error("Invalid instrument type: ", typeof(inst))
   end
end

"""  `generate_spectra_timeseries(line_list, inst; time, rv, ssbz, snr_per_pixel, line_width, add_noise )`
Generate a time series of spectra using times, a line_list (as DataFrame with columns lambda and weight) and a theoretical instrument.
Optionally, specify times, rv's, barycentric corrections, etc.
"""
function generate_spectra_timeseries(times::AbstractArray, line_list::DataFrame, inst::AnyTheoreticalInstrument, rvs::AbstractArray; ssbzs::AbstractArray = zeros(length(rvs)),
               snr_per_pixel::Real = 1000.0, line_width::Real = predict_intrinsic_stellar_line_width(5780) ) # where { LLT<:AbstractLineList }
  @assert length(times) == length(rvs) == length(ssbzs)
  @assert hasproperty(line_list,:lambda)
  @assert hasproperty(line_list,:weight)
  spectra = map(i->generate_spectrum(line_list, inst, snr_per_pixel=snr_per_pixel, line_width=line_width, time=times[i], rv=rvs[i], ssbz=ssbzs[i] ), 1:length(times) )
end

""" `demo_generate_spectrum_line( inst )`
Generate spectrum with one line for testing purposes.
"""
function demo_generate_spectrum_line(inst::AnyTheoreticalInstrument)
   ll = DataFrame(:lambda=>[5500.0], :weight=>[0.6] )
   generate_spectrum(ll, inst)
end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
module RayTracing

using UnPack
using RecipesBase
using StaticArrays
using IntervalSets
using LinearAlgebra
using NearestNeighbors
using Gridap: VectorValue, num_cells, get_grid
using Gridap.ReferenceFEs: get_faces, get_node_coordinates, num_cell_dims
using Gridap.Geometry: UnstructuredDiscreteModel, UnstructuredGrid, get_grid_topology,
    get_cell_node_ids

import Base: show, in
import NearestNeighbors: KDTree
import Gridap.Geometry: num_cells
import Gridap.ReferenceFEs: num_dims, num_nodes

include("point.jl")
include("segment.jl")
include("boundary.jl")
include("track.jl")
include("mesh.jl")
include("intersection.jl")
include("azimuthal_quad.jl")
include("trackgenerator.jl")
include("plot_recipes.jl")

export BoundaryConditions, Vaccum, Reflective, Periodic
export TrackGenerator
export trace!
export segmentize!

end


text:
## make options for the stationary GP
nmin, nmax = 2, 40
pnorm = 2.
K = transD_GP.GP.OrstUhn()
demean = true
fbounds = [-0.5 2.5]
sdev_pos = [0.05abs(diff([extrema(znall)...])[1])]
sdev_prop = 0.05diff(fbounds, dims=2)[:]
xall = permutedims(collect(znall))
xbounds = permutedims([extrema(znall)...])
λ, δ = [2], 0.1
fdataname = "nutest_"
## Initialize a stationary GP using these options
Random.seed!(12)
opt = transD_GP.OptionsStat(nmin = nmin,
                        nmax = nmax,
                        xbounds = xbounds,
                        fbounds = fbounds,
                        xall = xall,
                        λ = λ,
                        δ = δ,
                        fdataname = fdataname,
                        demean = demean,
                        sdev_prop = sdev_prop,
                        sdev_pos = sdev_pos,
                        pnorm = pnorm,
                        save_freq = 50,
                        quasimultid = false,
                        K = K,
                        updatenuisances = true
                        )
## nuisance options
# we don't treat *every* geometry param as nuisance.
# this is achieved by having the fixed parameters set as nuisances with
# zero-width bounds. These values are *fractions*
optn = transD_GP.OptionsNuisance(opt;
    sdev = [
        0.0 #zTx
        0.005 #zRx
        0.005 #x_rx
        0.0 #y_rx
        0.0 #rx_roll
        0.005 #rx_pitch
        0.0 #rx_yaw
        0.0 #tx_roll
        0.0 #tx_pitch
        0.0 #tx_yaw
           ],
    bounds = [
      zTx               zTx
      zRx - 5.0         zRx + 5.0
      x_rx - 5.0        x_rx + 5.0
      y_rx              y_rx
      rx_roll           rx_roll
      rx_pitch - 2.0    rx_pitch + 2.0
      rx_yaw            rx_yaw
      tx_roll           tx_roll
      tx_pitch          tx_pitch
      tx_yaw            tx_yaw
          ],
    updatenuisances = true,
    C = [0.160048    0.0131349   -0.183764
        0.0131349   0.00132317  -0.0154768
        -0.183764   -0.0154768    0.211669]
        )


text:
using RegressionFormulae
using StatsModels
using Test

@testset "fulldummy" begin
    include("fulldummy.jl")
end

@testset "nesting" begin
    include("nesting.jl")
end

@testset "powers of terms" begin
    include("power.jl")
end


text:
function solve()
    k = parse(Int, readline())
    a, b = [parse(Int, x) for x in split(readline())]

    if (a ÷ k == b ÷ k) && a % k != 0
        "NG"
    else
        "OK"
    end
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
"""
    SpeciesDirectory

Object for maintaining lists of species and their names. It can be queried by species name with the `in` function and the same indexing notation used for `Dict`s.

Each species is assigned an index starting from 1, and optionally a unique name. Species can be added and deleted with the functions `addspecies!` and `delete!`. Note that deleting a species will cause a shift in the indices of all the species that followed it.

See also: `addpsecies!`, `delete!`, `rename!`
"""
struct SpeciesDirectory
    list::Vector{String}
    dict::Dict{String,Int}

    SpeciesDirectory() = new(String[], Dict{String,Int}())
end # struct SpeciesDirectory

"""
    SpeciesDirectory(names)

Cerate a species directory from a list of unique `names` given as a string vector.
"""
function SpeciesDirectory(namelist::Vector{String})
    propernames = filter(x -> x ≠ "", namelist)
    allunique(propernames) || throw(SpeciesNameCollision())
    dir = SpeciesDirectory()
    push!(dir.list, namelist...)
    for (i, name) ∈ enumerate(dir.list)
        name == "" && continue
        dir.dict[name] = i
    end

    return dir
end

Base.length(dir::SpeciesDirectory) = length(dir.list)

Base.getindex(dir::SpeciesDirectory, str::String) = dir.dict[str]

function Base.getindex(dir::SpeciesDirectory, names::Vector{String})
    return map(x -> getindex(dir, x), names)
end

function Base.getindex(dir::SpeciesDirectory, names::String...)
    return collect(map(x -> getindex(dir, x), names))
end

"""
    in(i, dir) -> Bool
    in(name, dir) -> Bool

Return true if there is a species with the index number `i` or with the given `name` in the species directory `dir`.
"""
function Base.in(i::Int, dir::SpeciesDirectory)::Bool
    i < 1 && throw(BoundsError(dir.list, i))
    return length(dir.list) ≥ i ? true : false
end

Base.in(str::String, dir::SpeciesDirectory) = haskey(dir.dict, str)

"""
    addpsecies!(dir [, name])

Add a new species to the species directory `dir`. `name` can be ommited to create an unnamed species.
"""
function addspecies!(dir::SpeciesDirectory, speciesname::String)
    haskey(dir.dict, speciesname) && throw(SpeciesNameCollision([speciesname]))
    push!(dir.list, speciesname)
    i = length(dir)
    if speciesname ≠ ""
        dir.dict[speciesname] = i
    end
    return dir
end # function push!

addspecies!(dir::SpeciesDirectory) = addspecies!(dir, "")

"""
    rename!(dir, oldname, newname)
    rename!(dir, ind, newname)

Change the name of a species in the SpeciesDirectory `dir` from `oldname` to `newname`.

The species index `ind` can also be given instead of `oldname`. This way it is possible to give a name to a species that currently doesn't have one.
"""
function rename!(dir::SpeciesDirectory, i::Int, newname::String)
    haskey(dir.dict, newname) && throw(SpeciesNameCollision(newname))
    oldname = dir.list[i]
    delete!(dir.dict, oldname)
    dir.dict[newname] = i
    dir.list[i] = newname

    return nothing
end # function rename!

rename!(dir::SpeciesDirectory, oldname::String, newname::String) =
    rename!(dir, dir.dict[oldname], newname)

"""
    delete!(dir, ind)
    delete!(dir, name)

Delete the species with index `ind` or the given `name` from the species directory `dir`.

Note that deleting a species will cause a shift in the indices of all the species that followed it.
"""
function Base.delete!(dir::SpeciesDirectory, ind::Int)
    for (k, v) in pairs(dir.dict)
        if v > ind
            dir.dict[k] = v - 1
        end
    end
    dir.list[ind] ≠ "" && delete!(dir.dict, dir.list[ind])
    deleteat!(dir.list, ind)

    return nothing
end # function delete!

Base.delete!(dir::SpeciesDirectory, spname::String) =
    delete!(dir, dir.dict[spname])


text:
"""
    extract(archive, directory, [files]; [verbose=false]) -> Nothing

Extracts files from a compressed tar `archive` to the specified `directory`. If `files` is
specified only the files given will be extracted. The `verbose` flag can be used to display
additional information to STDOUT.
"""
function extract(archive, directory, files=AbstractString[]; verbose::Bool=false)
    @static if Sys.iswindows()
        cmd = pipeline(`$exe7z x $archive -y -so`, `$exe7z x -si -y -ttar -o$directory $files`)
    else
        cmd = `tar xvf $archive --directory=$directory $files`
    end

    if !verbose
        cmd = pipeline(cmd, stdout=devnull, stderr=devnull)
    end

    run(cmd)
end

"""
    isarchive(path) -> Bool

Determines if the given `path` is an archive.
"""
function isarchive(path)
    @static if Sys.iswindows()
        success(`$exe7z t $path -y`)
    else
        success(`tar tf $path`)
    end
end

"""
    readarchive(archive) -> Vector{AbstractString}

Returns the file names contained in the `archive`.
"""
function readarchive(archive)
    @static if Sys.iswindows()
        files = AbstractString[]
        header = "-" ^ 24
        content = false

        cmd = pipeline(`$exe7z x $archive -y -so`, `$exe7z l -si -y -ttar`)
        output = readchomp(pipeline(cmd, stderr=devnull))
        for line in split(output, "\r\n")
            # Extract the file name from the last column in the 7-zip output table.
            # Note: We can just write `line[54:end]` on Julia 0.5 and up (JuliaLang/julia#15624)
            len = length(line)
            file = len >= 54 ? line[54:len] : ""

            if file == header
                if !content
                    content = true
                else
                    break
                end
            elseif content
                push!(files, file)
            end
        end

        return files
    else
        output = readchomp(pipeline(`tar tf $archive`, stderr=devnull))
        return split(output, '\n')
    end
end


text:
"""
    map_to_canonical(x, xmin, xmax, basis::NodalBasis{Line})

Map the given coordinate `x` in the interval [`xmin`, `xmax`] to the corresponding
canonical coordinate `ξ` in the interval [-1, 1].
"""
function map_to_canonical(x, xmin, xmax, basis::NodalBasis{Line})
    ( 2x - (xmax + xmin) ) / (xmax - xmin)
end

"""
    map_to_canonical!(ξ, x, xmin, xmax, basis::NodalBasis{Line})

Map the given coordinate `x` in the interval [`xmin`, `xmax`] to the corresponding
canonical coordinate `ξ` in the interval [-1, 1], updating `ξ`.
"""
function map_to_canonical!(ξ, x, xmin, xmax, basis::NodalBasis{Line})
    @boundscheck begin
        @assert size(ξ) == size(x)
    end
    @inbounds for i in eachindex(ξ)
        ξ[i] = map_to_canonical(x[i], xmin, xmax, basis)
    end
    ξ
end

"""
    map_from_canonical(ξ, xmin, xmax, basis::NodalBasis{Line})

Map the given canonical coordinate `ξ` in the interval [-1, 1] to the corresponding
coordinate `x` in the interval [`xmin`, `xmax`].
"""
function map_from_canonical(ξ, xmin, xmax, basis::NodalBasis{Line})
    ( (xmax + xmin) + ξ * (xmax - xmin) ) / 2
end

"""
    map_from_canonical!(x, ξ, xmin, xmax, basis::NodalBasis{Line})

Map the given canonical coordinate `ξ` in the interval [-1, 1] to the corresponding
coordinate `x` in the interval [`xmin`, `xmax`], updating `x`.
"""
function map_from_canonical!(x, ξ, xmin, xmax, basis::NodalBasis{Line})
    @boundscheck begin
        @assert size(ξ) == size(x)
    end
    @inbounds for i in eachindex(x)
        x[i] = map_from_canonical(ξ[i], xmin, xmax, basis)
    end
    x
end


text:
using IncrementalInference
using Test

@testset "test default userdata::FactorMetadata..." begin

fgt = initfg()

addVariable!(fgt, :x1, ContinuousScalar)
addFactor!(fgt, [:x1], Prior(Normal()))

addVariable!(fgt, :x2, ContinuousScalar)
addFactor!(fgt, [:x1;:x2], LinearConditional(Normal(10,1)))

fc = solverData(getFactor(fgt, :x1x2f1))
@test length(fc.fnc.cpt[1].factormetadata.variableuserdata) == 2
@test fc.fnc.cpt[1].factormetadata.solvefor == :null

end

#


text:
@testset "polynomials" begin
    @test DiscretisedFluidQueues.gauss_lobatto_points(2.0,2.5,1)≈[2.25]
    @test DiscretisedFluidQueues.gauss_lobatto_points(2.0,3.0,2)≈[2.0;3.0]
    @test DiscretisedFluidQueues.gauss_lobatto_points(2.0,4.0,3)≈[2.0;3.0;4.0]
    nnodes = 5
    nodes = DiscretisedFluidQueues.gauss_lobatto_points(-1.0,1.0,nnodes)
    for n in 1:nnodes
        en = zeros(nnodes)
        en[n] = 1.0
        @test DiscretisedFluidQueues.lagrange_polynomials(nodes,nodes[n])≈en atol=sqrt(eps())
    end
    @test DiscretisedFluidQueues.lagrange_polynomials([1.0],1.0)≈[1.0] atol=sqrt(eps())
    @test DiscretisedFluidQueues.gauss_lobatto_weights(-1.0,1.0,1)≈[2.0]
    nnodes = 3
    nodes = DiscretisedFluidQueues.gauss_lobatto_points(-1.0,1.0,nnodes)
    @test DiscretisedFluidQueues.gauss_lobatto_weights(-1.0,1.0,nnodes)≈[1.0,4.0,1.0]./3
    x4 = x->x^4
    x4_interp = DiscretisedFluidQueues.lagrange_interpolation(x4,-1.0,1.0,5)
    x = collect(-3:0.02:3)
    @test x4_interp.(x)≈x4.(x)
    x4_quad = DiscretisedFluidQueues.gauss_lobatto_quadrature(x4,-1.0,1.0,4)
    @test x4_quad≈2/5

    x4_interp_p0 = DiscretisedFluidQueues.lagrange_interpolation(x->1.0-x4(x),-1.0,1.0,1)
    x = collect(-3:0.02:3)
    @test x4_interp_p0.(x)≈ones(size(x))
    x4_quad_p0 = DiscretisedFluidQueues.gauss_lobatto_quadrature(x->1.0-x4(x),-1.0,1.0,1)
    @test x4_quad_p0≈2.0
end


text:
@testset "periodic" begin
    A = BlockTridiagonal(Vcat([[0. 1.; 0. 0.]],Fill([0. 1.; 0. 0.], ∞)),
                        Vcat([[-1. 1.; 1. 1.]], Fill([-1. 1.; 1. 1.], ∞)),
                        Vcat([[0. 0.; 1. 0.]], Fill([0. 0.; 1. 0.], ∞)))

    Q,L = ql(A);
    @test Q.factors isa InfiniteLinearAlgebra.InfBlockBandedMatrix
    Q̃,L̃ = ql(BlockBandedMatrix(A)[Block.(1:100),Block.(1:100)])

    @test Q̃.factors[1:100,1:100] ≈ Q.factors[1:100,1:100]
    @test Q̃.τ[1:100] ≈ Q.τ[1:100]
    @test L[1:100,1:100] ≈ L̃[1:100,1:100]
    @test Q[1:10,1:10] ≈ Q̃[1:10,1:10]
    @test Q[1:10,1:12]*L[1:12,1:10] ≈ A[1:10,1:10]

    # complex non-selfadjoint
    c,a,b = [0 0.5; 0 0],[0 2.0; 0.5 0],[0 0.0; 2.0 0];
    A = BlockTridiagonal(Vcat([c], Fill(c,∞)),
                    Vcat([a], Fill(a,∞)),
                    Vcat([b], Fill(b,∞))) - 5im*I
    Q,L = ql(A)
    @test Q[1:10,1:12]*L[1:12,1:10] ≈ A[1:10,1:10]

    c,a,b = [0 0.5; 0 0],[0 2.0; 0.5 0],[0 0.0; 2.0 0];
    A = BlockTridiagonal(Vcat([c], Fill(c,∞)),
                    Vcat([a], Fill(a,∞)),
                    Vcat([b], Fill(b,∞)))
    Q,L = ql(A)
    @test Q[1:10,1:12]*L[1:12,1:10] ≈ A[1:10,1:10]
    @test L[1,1] == 0 # degenerate

    Q,L = ql(A')
    @test Q[1:10,1:12]*L[1:12,1:10] ≈ A[1:10,1:10]'
    @test L[1,1]  ≠ 0 # non-degenerate
end


text:
using RDatasets, Gadfly
Gadfly.set_default_plot_size(14cm, 8cm)

D = dataset("datasets","faithful")
D.g = D.Eruptions.>3.0

coord = Coord.cartesian(ymin=35, ymax=100)

pa = plot(D, coord,
    x=:Eruptions, y=:Waiting, group=:g,
    Geom.point, Geom.ellipse
)
pb = plot(D, coord,
    x=:Eruptions, y=:Waiting, color=:g,
    Geom.point, Geom.ellipse,
    layer(Geom.ellipse(levels=[0.99]), style(line_style=[:dot])),
    style(key_position=:none), Guide.ylabel(nothing)
)
hstack(pa,pb)


text:
## ----------------- SNode ------------------
struct SNode<:AbstractNode
    id::Int
    data::Dict{Symbol,Vector{Float64}}
    index::Vector{Int}
    𝝭::Dict{Symbol,Vector{Float64}}
end

## convert
Node(ξ::SNode) = Node(ξ.id,ξ.data)
SNode(ξ::T,η::SNode) where T<:AbstractNode = SNode(ξ.id,ξ.data,η.index,η.𝝭)
function (η::SNode)(ξ::SNode)
    ξ.index[ξ.id] = η.index[η.id]
    empty!(ξ.𝝭)
    for s in keys(η.𝝭)
        ξ.𝝭[s] = η.𝝭[s]
    end
end

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator


text:
#TODO: finish
function hurst(x::Array{Float64}; n::Int=10; C::Float64=Inf)
end

#TODO: finish
function hilbert_transform(x::Array{Float64}; n::Int=10)
end


text:
module DataAndPlotting
using HDF5
#using PyPlot
#using Plots

const ni=200
const nj=300

const ki=3
const kj=3

function indextocoord(i,ni)
    @assert ni>1
    @assert 1<= i <= ni
    (i-1)/(ni-1)
    end

#Create standing wave (our data set)
export initialize
function initialize()
    arr=zeros(Float64,ni,nj)
    for j in 1:nj, i in 1:ni
        x=indextocoord(i,ni)
        y=indextocoord(j,nj)
        arr[i,j]=sin(2π*ki*x) * sin(2π*kj*y)
    end
    arr
end

#Write to file in CSV form
export output
function output(arr::Array{Float64,2}, filename::String)
    h5write(filename, "data", arr)
    end

#Read back
export readfile
function readfile(filename::String)
    data=h5read(filename, "data")
    end

#Plot
 export myplot
function myplot(arr::Array{Float64,2},filename::String)
    #pyplot()
    heatmap(arr,clim=(-1.0,+1.0),color=:viridis)
    savefig(filename)
    end

end # module

max_stars_repo_name: laurasberna/DataAndPlotting

max_issues_repo_name: laurasberna/DataAndPlotting

max_forks_repo_name: laurasberna/DataAndPlotting


text:
__precompile__()

module NamedDimensions
using FunctionalData

export NamedDims, named, array, names, at
immutable NamedDims{T,N}
    data::Array{T,N}
    names::Array{Symbol}
    function NamedDims(a, names...)
        names = collect(names)
        length(names) > ndims(a) && error("NamedDimensions: Got $(length(names)) names $names for $(ndims(a)) dims $(size(a))")
        for i = 1:len(names), j = i+1:len(names)
            names[i]==names[j] && error("NamedDimensions: Duplicate names: $names")
        end

        if length(names) < ndims(a)
            names = vcat(dimnames(names, ndims(a)-length(names)), names)
        end
        new(a, names)
    end
end

named{T}(a::AbstractArray{T,0}, named::Vector{Symbol}) = a[1]
named(a::AbstractArray, names::Vector{Symbol}) = named(a, names...)
named(a::AbstractArray, names...) = (assert(all(x->isa(x,Symbol), names)); NamedDims{eltype(a),ndims(a)}(a, names...))
named(a::NamedDims, names::Vector{Symbol}) = named(a, names...)
named(a, args...) = a

import Base: getindex, eltype, isempty
array(a) = a
array(a::NamedDims) = a.data
array(a::NamedDims, inds...) = array(named(a, inds...))
getindex(a::NamedDims, inds...) = named(a, inds...)
eltype{T,N}(a::NamedDims{T,N}) = T
isempty(a::NamedDims) = len(a) == 0

for f in [:fst, :snd, :third, :last]
    @eval import FunctionalData.$f
    @eval $f(a::NamedDims) = named(squeezer($f(a.data)), droplast(a.names))
    @eval $f(a::NamedDims, x::Symbol) = $f(named(a, x))
end

import FunctionalData.at
at(a::NamedDims, i::Int) = a[last(a.names) => i]

for f in [:take, :takelast, :drop, :droplast, :partition, :partsoflen]
    @eval import FunctionalData.$f
    @eval $f(a::NamedDims, i::Int) = named($f(a.data,i), a.names)
end

for f in [:part, :trimmedpart]
    @eval import FunctionalData.$f
    @eval $f{T}(a::NamedDims, i::DenseArray{T,2}) = named($f(a.data,i), a.names)
    @eval $f(a::NamedDims, i::UnitRange) = named($f(a.data,i), a.names)
    @eval $f(a::NamedDims, i::AbstractVector) = named($f(a.data,i), a.names)
end

import FunctionalData.concat
function concat(a::NamedDims, b::NamedDims)
    assert(a.names == b.names)
    assert(size(fst(a)) == size(fst(b)))
    named(concat(a.data, b.data), a.names)
end

import FunctionalData.flatten
function flatten{T,N}(a::Vector{NamedDims{T,N}})
    assert(all(x->x.names == fst(a).names, a))
    assert(all(x->size(x)[1:end-1] == size(fst(a))[1:end-1], a))
    named(flatten(map(a,array)), fst(a).names)
end

import FunctionalData.stack
function stack(a, name::Symbol)
    assert(all(x->x.names == fst(a).names, a))
    assert(all(x->size(x) == size(fst(a)), a))
    named(stack(map(a,array)), concat(fst(a).names, name))
end

function stack(r, a::NamedDims)
    isempty(r) && return []
    if isa(fst(r), NamedDims)
        return stack(typed(r), last(a.names))
    else
        return stack(r)
    end
end

import FunctionalData: map, map2, map3, map4, map5

function map{T,N}(a::NamedDims{T,N}, f::Function)
    r = [f(at(a,i)) for i in 1:len(a)]
    return stack(r, a)
end
mapper2(f, N, a, b, fa, fb) = stack([f(fa(at(a,i)), fb(at(b,i))) for i in 1:len(N)], N)
squeezer(a::AbstractArray) = size(a,ndims(a)) == 1 ? squeeze(a,ndims(a)) : a
squeezer(a) = a
map2(a::NamedDims, b::NamedDims, f::Function) = mapper2(f,a,a,b,id,id)
map2(a::NamedDims, b, f::Function) = mapper2(f,a,a,b,id,squeezer)
map2(a, b::NamedDims, f::Function) = mapper2(f,b,a,b,squeezer,id)

at(a::NamedDims, inds::Tuple) = at(a.data, inds)
at(a::NamedDims, inds...) = named(a, inds...)
call(a::NamedDims, inds...) = named(a, inds...)

names(a::NamedDims) = a.names

import Base.==
==(a::NamedDims, b::NamedDims) = a.names == b.names && a.data == b.data

import Base.size
size(a::NamedDims) = size(a.data)
size(a::NamedDims, dim::Int) = size(a.data, dim)
size(a::NamedDims, dim::Symbol) = size(a.data, getind(a, dim))
import Base.ndims
ndims(a::NamedDims) = ndims(a.data)
import Base.length
length(a::NamedDims) = length(a.data)
import FunctionalData.len
len(a::NamedDims) = len(a.data)

function getind(a::NamedDims, x::Symbol)
    ind = findfirst(x.==a.names)
    if ind == 0
        error("NamedDimensions: No dimension $x in $(a.names)")
    end
    ind
end

function named(a::NamedDims, inds...)
    inds = [isa(x, Pair) && x.second == :end ? Pair(x.first, size(a, x.first)) : x for x in inds]
    isname(x) = (isa(x, Pair) && isa(x.first, Symbol)) || isa(x, Symbol)
    issqueezer(x) = isa(x, Pair) && isa(x.second, Integer)
    if length(inds)==1 && !isname(inds[1])
        return a.data[inds...]
    end
    if all(x->!isa(x, Pair), inds) && any(x->!isa(x, Symbol), inds)  # pairs of arguments
        if length(inds) % 2 != 0
            error("NamedDimensions: Expected an even number of parameters, got $(inds)")
        end
        inds = @p map2 inds[1:2:end] inds[2:2:end] Pair
        return named(a, inds...)
    end

    [assert(isa(x, Pair) || isa(x, Symbol)) for x in inds]
    inds = [isa(x, Pair) ? x : Pair(x, Colon()) for x in inds]
    ordered    = @p filter inds not*issqueezer | map fst
    squeezed   = @p filter inds issqueezer | map fst
    invalid = @p concat ordered squeezed | filter not*inside a.names
    if !isempty(invalid)
        error("NamedDimensions: Invalid dimensions $invalid . Valid dimensions are $(a.names)")
    end
    implicitorder = @p filter a.names not*inside squeezed | filter not*inside ordered
    ordered = @p concat implicitorder ordered

    d = Dict(inds...)
    inds = [get(d, x, Colon())  for x in a.names]
    newnames = Symbol[]
    for i in 1:length(a.names)
        isa(inds[i], Integer) ? nothing : push!(newnames, a.names[i])
    end

    b = (VERSION < v"0.5-" ? slice : view)(a.data, inds...)
    if ordered == newnames
        return named(b, newnames)
    else
        perm = [findfirst(x .== newnames) for x in ordered]
        return named(permutedims(b, perm), newnames[perm])
    end
end

import Base.summary
function summary{T,N}(a::NamedDims{T,N})
    @p mapvec2 a.names size(a.data) ((n,s) -> "$s $n") | join " x " | concat " :: NamedDims{$T,$N}"
end

import Base.show
show(io::IO, a::NamedDims) = (println(io, summary(a)); show(io,a.data))

import FunctionalData.showinfo
function showinfo(io::IO, a::NamedDims, comment::AbstractString = "")
    print(io, comment, length(comment) > 0 ? "  --  ": "----  ")
    @p mapvec2 a.names size(a.data) ((n,s) -> "$s $n") | join " x " | println io _
    showinfo(io, a.data; showheader = false)
    a
end

for f in [:sqrt]
    @eval import Base.$f
    @eval $f(a::NamedDims, args...) = named($f(a.data, args...), a.names)
end

for f in [:power, :square]
    @eval import FunctionalData.$f
    @eval $f(a::NamedDims, args...) = named($f(a.data, args...), a.names)
end

for f in [:mean, :maximum, :minimum, :median, :std, :var, :sum]
    @eval import Base.$f
    @eval $f(a::NamedDims) = $f(a.data)
    @eval $f(a::NamedDims, dim::Int) = named(squeeze($f(a.data,dim),dim), dropat(a.names,dim))
    @eval $f(a::NamedDims, dim::Symbol) = $f(a, getind(a, dim))
end

function broadcast(f::Function, smaller::NamedDims, larger::NamedDims)
    assert(all(x->in(x, larger.names), smaller.names))
    last = @p filter larger.names not*inside smaller.names
    orderforsmaller = @p filter larger.names inside smaller.names
    a = isempty(last) ? larger : named(larger, last)
    reorderedsmaller = smaller[orderforsmaller].data
    r = @p f reorderedsmaller a.data
    named(r, concat(orderforsmaller, last))[larger.names]
end

for f in [:.+, :.-, :.*, :./, :.\, :.^]
    @eval import Base.$f
    @eval $f(a::NamedDims, b::AbstractArray) = ndims(a.data) >= ndims(b)-1 ? named($f(a.data, b), a.names) : $f(a.data, b)
    @eval $f(a::AbstractArray, b::NamedDims) = ndims(b.data) >= ndims(a)-1 ? named($f(a, b.data), b.names) : $f(a, b.data)
    @eval function $f(a::NamedDims, b::NamedDims)
        length(a.names) < length(b.names) ? broadcast($f, a, b) : broadcast($f, b, a)
    end
end

for f in [:-, :+]
    @eval import Base.$f
    @eval $f(a::NamedDims, b::AbstractArray) = named($f(a.data, b), a.names)
    @eval $f(a::AbstractArray, b::NamedDims) = named($f(a, b.data), b.names)
    @eval function $f(a::NamedDims, b::NamedDims)
        assert(a.names == b.names)
        named($f(a.data, b.data), a.names)
   end
end

chars = collect('A':'Z')
dimnames(a::NamedDims, args...) = dimnames(a.names, args...)
function dimnames(names, n = 2)
    @p filter [Symbol("dim$(chars[i])") for i in 1:n+len(names)] not*inside names | take n
end

for (f,inds) in [(:*,[1,4]), (:/, [1,3]), (:\, [2,4])]
    @eval import Base.$f
    @eval function $f(a::NamedDims, b::NamedDims)
        named($f(a.data, b.data), concat(a.names, b.names)[$inds])
    end
    @eval $f(a::NamedDims, b::AbstractArray) = named($f(a.data, b), concat(a.names, dimnames(a))[$inds])
    @eval $f(a::AbstractArray, b::NamedDims) = named($f(a, b.data), concat(dimnames(a), b.names)[$inds])
end

end # module


text:
@testset "Quantum walks simulator" begin

  n = 5
  g = complete_graph(n)

  @testset "General" begin
    qwe = QWEvolution(Szegedy(g))

    @test_throws AssertionError QWEvolution(Szegedy(g), Dict{Symbol,Any}(), check=true)
    @test_no_error QWEvolution(Szegedy(g), Dict{Symbol,Any}(), check=false)
    @test check_qwdynamics(qwe) == nothing
  end

  @testset "Quantum walks simulator for CTQW" begin
    init = fill(ComplexF64(1/sqrt(n)), n)
    ctqwAdj = CTQW(g, :adjacency)
    ctqwLap = CTQW(g, :laplacian)
    ctqwAdjDense = CTQWDense(g, :adjacency)
    ctqwLapDense = CTQWDense(g, :laplacian)
    qweAdj = QWEvolution(ctqwAdj)
    qweLap = QWEvolution(ctqwLap)
    qweAdjDense = QWEvolution(ctqwAdjDense)
    qweLapDense = QWEvolution(ctqwLapDense)

    @test_throws AssertionError execute(qweAdj, init , -1.0)

    @test check_qwdynamics(qweAdj) == nothing

    @test execute(qweAdj, init, 1) ≈ execute(qweAdj, execute(qweAdj, init, 0.5), 0.5)

    @test_throws AssertionError execute(qweLap, init , -1.0)
    @test execute(qweLap, init, 1) ≈ execute(qweLap, execute(qweLap, init, 0.5), 0.5)

    # due to sign issue the measure function is needed
    @test measure(qweAdj, execute(qweAdj, init, 1)) ≈ measure(qweLap, execute(qweLap, init, 1))
    @test measure(qweAdjDense, execute(qweAdjDense, init, 1)) ≈ measure(qweLapDense, execute(qweLapDense, init, 1))

    @test_throws ErrorException QWEvolution(CTQW(g, :notimplemented))
  end

  @testset "Quantum walks simulator for Szegedy" begin
    qwe = QWEvolution(Szegedy(g))
    init = spzeros(n*n)
    init[1] = init[n] = init[n*n] = 1/sqrt(3.)

    steps = 10

    @test check_qwdynamics(qwe) == nothing

    @test_throws AssertionError execute(qwe, init, -1)
    @test_throws MethodError execute(qwe, init, 0.5)

    @test execute(qwe, init, steps, all = true, measure = true)[:,steps+1] == execute(qwe, init, steps, all = false, measure = true)
    @test execute(qwe, init, steps, all = true, measure = false)[steps+1] == execute(qwe, init, steps, all = false, measure = false)

  end

end


text:
__precompile__(true)

module Lint

using Base.Meta
using Compat
using Compat: TypeUtils, readline
using JSON
using AutoHashEquals
using Printf

if isdefined(Base, :unwrap_unionall)
    using Base: unwrap_unionall
else
    unwrap_unionall(x) = x
end

export LintMessage, LintContext
export lintfile, lintstr, lintpkg, lintserver, @lintpragma
export iserror, iswarning, isinfo
export test_similarity_string

const SIMILARITY_THRESHOLD = 10.0

# no-op. We have to use macro inside type declaration as it disallows actual function calls
macro lintpragma(s)
end

include("exprutils.jl")
using .ExpressionUtils

include("statictype.jl")

include("linttypes.jl")
include("guesstype.jl")
include("result.jl")

# AST Linting
include("abstracteval.jl")
include("ast.jl")
include("blocks.jl")
include("controls.jl")
include("curly.jl")
include("dict.jl")
include("dynamic.jl")
include("functions.jl")
include("generator.jl")
include("include.jl")
include("knowndeprec.jl")
include("macros.jl")
include("misc.jl")
include("modules.jl")
include("pragma.jl")
include("ref.jl")
include("types.jl")
include("variables.jl")

# Command-Line Interface
include("cli.jl")

# Server
include("server.jl")

# precompile hints
include("init.jl")
include("precompile.jl")

end


text:
function big_fraction_matrix_to_real_matrix(::Type{MatrixUtils}, arg0::FieldMatrix)
    return jcall(MatrixUtils, "bigFractionMatrixToRealMatrix", Array2DRowRealMatrix, (FieldMatrix,), arg0)
end

function block_inverse(::Type{MatrixUtils}, arg0::RealMatrix, arg1::jint)
    return jcall(MatrixUtils, "blockInverse", RealMatrix, (RealMatrix, jint), arg0, arg1)
end

function check_addition_compatible(::Type{MatrixUtils}, arg0::AnyMatrix, arg1::AnyMatrix)
    return jcall(MatrixUtils, "checkAdditionCompatible", void, (AnyMatrix, AnyMatrix), arg0, arg1)
end

function check_column_index(::Type{MatrixUtils}, arg0::AnyMatrix, arg1::jint)
    return jcall(MatrixUtils, "checkColumnIndex", void, (AnyMatrix, jint), arg0, arg1)
end

function check_matrix_index(::Type{MatrixUtils}, arg0::AnyMatrix, arg1::jint, arg2::jint)
    return jcall(MatrixUtils, "checkMatrixIndex", void, (AnyMatrix, jint, jint), arg0, arg1, arg2)
end

function check_multiplication_compatible(::Type{MatrixUtils}, arg0::AnyMatrix, arg1::AnyMatrix)
    return jcall(MatrixUtils, "checkMultiplicationCompatible", void, (AnyMatrix, AnyMatrix), arg0, arg1)
end

function check_row_index(::Type{MatrixUtils}, arg0::AnyMatrix, arg1::jint)
    return jcall(MatrixUtils, "checkRowIndex", void, (AnyMatrix, jint), arg0, arg1)
end

function check_same_column_dimension(::Type{MatrixUtils}, arg0::AnyMatrix, arg1::AnyMatrix)
    return jcall(MatrixUtils, "checkSameColumnDimension", void, (AnyMatrix, AnyMatrix), arg0, arg1)
end

function check_same_row_dimension(::Type{MatrixUtils}, arg0::AnyMatrix, arg1::AnyMatrix)
    return jcall(MatrixUtils, "checkSameRowDimension", void, (AnyMatrix, AnyMatrix), arg0, arg1)
end

function check_sub_matrix_index(::Type{MatrixUtils}, arg0::AnyMatrix, arg1::Vector{jint}, arg2::Vector{jint})
    return jcall(MatrixUtils, "checkSubMatrixIndex", void, (AnyMatrix, Vector{jint}, Vector{jint}), arg0, arg1, arg2)
end

function check_sub_matrix_index(::Type{MatrixUtils}, arg0::AnyMatrix, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(MatrixUtils, "checkSubMatrixIndex", void, (AnyMatrix, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function check_subtraction_compatible(::Type{MatrixUtils}, arg0::AnyMatrix, arg1::AnyMatrix)
    return jcall(MatrixUtils, "checkSubtractionCompatible", void, (AnyMatrix, AnyMatrix), arg0, arg1)
end

function check_symmetric(::Type{MatrixUtils}, arg0::RealMatrix, arg1::jdouble)
    return jcall(MatrixUtils, "checkSymmetric", void, (RealMatrix, jdouble), arg0, arg1)
end

function create_column_field_matrix(::Type{MatrixUtils}, arg0::Vector{FieldElement})
    return jcall(MatrixUtils, "createColumnFieldMatrix", FieldMatrix, (Vector{FieldElement},), arg0)
end

function create_column_real_matrix(::Type{MatrixUtils}, arg0::Vector{jdouble})
    return jcall(MatrixUtils, "createColumnRealMatrix", RealMatrix, (Vector{jdouble},), arg0)
end

function create_field_diagonal_matrix(::Type{MatrixUtils}, arg0::Vector{FieldElement})
    return jcall(MatrixUtils, "createFieldDiagonalMatrix", FieldMatrix, (Vector{FieldElement},), arg0)
end

function create_field_identity_matrix(::Type{MatrixUtils}, arg0::Field, arg1::jint)
    return jcall(MatrixUtils, "createFieldIdentityMatrix", FieldMatrix, (Field, jint), arg0, arg1)
end

function create_field_matrix(::Type{MatrixUtils}, arg0::Field, arg1::jint, arg2::jint)
    return jcall(MatrixUtils, "createFieldMatrix", FieldMatrix, (Field, jint, jint), arg0, arg1, arg2)
end

function create_field_matrix(::Type{MatrixUtils}, arg0::Vector{Vector{FieldElement}})
    return jcall(MatrixUtils, "createFieldMatrix", FieldMatrix, (Vector{Vector{FieldElement}},), arg0)
end

function create_field_vector(::Type{MatrixUtils}, arg0::Field, arg1::jint)
    return jcall(MatrixUtils, "createFieldVector", FieldVector, (Field, jint), arg0, arg1)
end

function create_field_vector(::Type{MatrixUtils}, arg0::Vector{FieldElement})
    return jcall(MatrixUtils, "createFieldVector", FieldVector, (Vector{FieldElement},), arg0)
end

function create_real_diagonal_matrix(::Type{MatrixUtils}, arg0::Vector{jdouble})
    return jcall(MatrixUtils, "createRealDiagonalMatrix", RealMatrix, (Vector{jdouble},), arg0)
end

function create_real_identity_matrix(::Type{MatrixUtils}, arg0::jint)
    return jcall(MatrixUtils, "createRealIdentityMatrix", RealMatrix, (jint,), arg0)
end

function create_real_matrix(::Type{MatrixUtils}, arg0::Vector{Vector{jdouble}})
    return jcall(MatrixUtils, "createRealMatrix", RealMatrix, (Vector{Vector{jdouble}},), arg0)
end

function create_real_matrix(::Type{MatrixUtils}, arg0::jint, arg1::jint)
    return jcall(MatrixUtils, "createRealMatrix", RealMatrix, (jint, jint), arg0, arg1)
end

function create_real_vector(::Type{MatrixUtils}, arg0::Vector{jdouble})
    return jcall(MatrixUtils, "createRealVector", RealVector, (Vector{jdouble},), arg0)
end

function create_real_vector(::Type{MatrixUtils}, arg0::jint)
    return jcall(MatrixUtils, "createRealVector", RealVector, (jint,), arg0)
end

function create_row_field_matrix(::Type{MatrixUtils}, arg0::Vector{FieldElement})
    return jcall(MatrixUtils, "createRowFieldMatrix", FieldMatrix, (Vector{FieldElement},), arg0)
end

function create_row_real_matrix(::Type{MatrixUtils}, arg0::Vector{jdouble})
    return jcall(MatrixUtils, "createRowRealMatrix", RealMatrix, (Vector{jdouble},), arg0)
end

function deserialize_real_matrix(::Type{MatrixUtils}, arg0::Object, arg1::JString, arg2::ObjectInputStream)
    return jcall(MatrixUtils, "deserializeRealMatrix", void, (Object, JString, ObjectInputStream), arg0, arg1, arg2)
end

function deserialize_real_vector(::Type{MatrixUtils}, arg0::Object, arg1::JString, arg2::ObjectInputStream)
    return jcall(MatrixUtils, "deserializeRealVector", void, (Object, JString, ObjectInputStream), arg0, arg1, arg2)
end

function fraction_matrix_to_real_matrix(::Type{MatrixUtils}, arg0::FieldMatrix)
    return jcall(MatrixUtils, "fractionMatrixToRealMatrix", Array2DRowRealMatrix, (FieldMatrix,), arg0)
end

function inverse(::Type{MatrixUtils}, arg0::RealMatrix)
    return jcall(MatrixUtils, "inverse", RealMatrix, (RealMatrix,), arg0)
end

function inverse(::Type{MatrixUtils}, arg0::RealMatrix, arg1::jdouble)
    return jcall(MatrixUtils, "inverse", RealMatrix, (RealMatrix, jdouble), arg0, arg1)
end

function is_symmetric(::Type{MatrixUtils}, arg0::RealMatrix, arg1::jdouble)
    return jcall(MatrixUtils, "isSymmetric", jboolean, (RealMatrix, jdouble), arg0, arg1)
end

function matrix_exponential(::Type{MatrixUtils}, arg0::RealMatrix)
    return jcall(MatrixUtils, "matrixExponential", RealMatrix, (RealMatrix,), arg0)
end

function serialize_real_matrix(::Type{MatrixUtils}, arg0::RealMatrix, arg1::ObjectOutputStream)
    return jcall(MatrixUtils, "serializeRealMatrix", void, (RealMatrix, ObjectOutputStream), arg0, arg1)
end

function serialize_real_vector(::Type{MatrixUtils}, arg0::RealVector, arg1::ObjectOutputStream)
    return jcall(MatrixUtils, "serializeRealVector", void, (RealVector, ObjectOutputStream), arg0, arg1)
end

function solve_lower_triangular_system(::Type{MatrixUtils}, arg0::RealMatrix, arg1::RealVector)
    return jcall(MatrixUtils, "solveLowerTriangularSystem", void, (RealMatrix, RealVector), arg0, arg1)
end

function solve_upper_triangular_system(::Type{MatrixUtils}, arg0::RealMatrix, arg1::RealVector)
    return jcall(MatrixUtils, "solveUpperTriangularSystem", void, (RealMatrix, RealVector), arg0, arg1)
end


text:
"""
    get_payoff(game::Dict{String,<:AbstractArray}, s::Vector{Vector{T}}) where T<:Real

Produces payoff profile for given game and actions probabilities array

**Input parameters**
* `game::Dict{String,<:AbstractArray}` - array of vectors to use in outer product
* `s::Vector{Vector{T}}` - collection of actions probabilities for each player in a game
"""
function get_payoff(game::Dict{String,<:AbstractArray},
                s::Vector{Vector{T}}) where T<:Real
    if !all(@. round(sum(s),digits = 2) ==1)
        error("Provided vectors are not probabilities distribution")
    end
    return Dict(k=>sum(v.*outer(s)) for (k,v) in game)
end

"""
    get_payoff(game::Dict{String,<:AbstractArray}, s::Vector{Vector{SymPy.Sym}})

Produces payoff profile for given game and actions probabilities symbolic array

**Input parameters**
* `game::Dict{String,<:AbstractArray}` - array of vectors to use in outer product
* `s::Vector{Vector{SymPy.Sym}}` - collection of actions probabilities for each player in a game
"""
function get_payoff(game::Dict{String,<:AbstractArray},
                s::Vector{Vector{SymPy.Sym}})
    return Dict(k=>sum(v.*outer(s)) for (k,v) in game)
end

"""
    best_reply(game::Dict{String,<:Array}, s::Vector{Vector{T}} where T<:Real,
    k::Int, epsil::F=0.0 ;return_val::String = "array") where F<:Real

Returns best reply of given player to actions of his counterparts
(defined by strategy profile) in the form of array or convex hull

**Input parameters**
* `game::Dict{String,<:Array}` - dictionary of players and their payoff matrices
* `s::Vector{Vector{T}} where T<:Real` - collection of actions probabilities for each player in a game
* `k::Int` - player for which the best reply is returned
* `epsil::F = 0.0` - probability of error (if the answer is disturbed)
* `return_val::String = "array"` - type of returned value ("array" or "chull")
"""
function best_reply(game::Dict{String,<:Array},
            s::Vector{Vector{T}} where T<:Real,
            k::Int,
            epsil::F = 0.0;
            return_val::String="array") where F<:Real

    (epsil < 0 || epsil > 1) && throw(ArgumentError("Probability of error (epsil) should be in range 0-1")) && return nothing

    action_no = size(game["player1"])[k]
    s_temp = deepcopy(s)
    payoffs = Vector{Real}()
    for i in 1:length(s_temp[k])
        s_temp[k] = zeros(length(s[k]))
        s_temp[k][i] = 1
        push!(payoffs, get_payoff(game, s_temp)["player"*string(k)])
    end
    eyemat = Matrix(1I, action_no, action_no)
    (epsil != 0 && rand() < epsil) ? pos = rand([false true], action_no) : pos = (payoffs .== maximum(payoffs))
    all(pos .== 0) ? pos[rand(1:action_no)] = 1 : nothing
    val2ret = eyemat[pos, :]
    if return_val == "chull"
        polyh = polyhedron(vrep(val2ret), CDDLib.Library())
        return convexhull(polyh, polyh)
    else
        return val2ret
    end
end

"""
    is_nash_q(game::Dict{String,<:Array}, s::Vector{Vector{T}}) where T<:Real

Returns a dictionary of logical vectors indicating whether given profile is Nash equilibrium.

**Input parameters**
* `game::Dict{String,<:Array}` - dictionary of players and their payoff matrices
* `s::Vector{Vector{T}})` - collection of actions probabilities for each player
"""
function is_nash_q(game::Dict{String,<:Array},
                s::Vector{Vector{T}}) where T<:Real

    nashqs = Dict{String,Bool}()
    for k in 1:length(game)
        st = s[k]
        br = best_reply(game, s, k, return_val = "chull")
        polyh = polyhedron(vrep([st]), CDDLib.Library())
        pint = intersect(br, polyh)
        nashqs["player"*string(k)] = (npoints(pint) != 0)
    end
    nashqs["Nash"] = all(values(nashqs))
    return nashqs
end

"""
    iterate_best_reply(game::Dict{String,<:Array},
            s::Vector{Vector{T}} where T<:Real,
            epsil::F = 0.0, it_num::Int = 10)  where F<:Real

Returns an array of arrays representing each players'
game strategies (action probabilities) history in every iteration

**Input parameters**
* `game::Dict{String,<:Array}` - dictionary of players and their payoff matrices
* `s::Vector{Vector{T}} where T<:Real` - collection of actions probabilities for each player
* `epsil::F = 0.0` - probability of error (if the answer is disturbed)
* `it_num::Int = 10` - an integer specifying how many iterations should be repeated
"""
function iterate_best_reply(game::Dict{String,<:Array},
    s::Vector{Vector{T}} where T<:Real,
    epsil::F = 0.0,
    it_num::Int = 10)  where F<:Real #TODO epsil not default
    s_temp = deepcopy(s)
    s_history = [] #TODO not Any type
    push!(s_history, s_temp)
    for i in 1:it_num
        br = Array{Int64}[]
        for i in 1:length(game)
            push!(br, best_reply(game, s_temp, i, epsil))
        end
        s_temp = _select_random(game, s_temp, br)
        push!(s_history, s_temp)
    end
    println("After ", it_num, " iterations")
    all(values(is_nash_q(game, s_history[length(s_history)]))) ?
        println("Nash equilibrium found") : println("Nash equilibrium not found")
    return s_history
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libgeotiff_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libgeotiff")
JLLWrappers.@generate_main_file("libgeotiff", UUID("06c338fa-64ff-565b-ac2f-249532af990e"))
end  # module libgeotiff_jll


text:
@params mutable struct BESOResult{T}
    topology::AbstractVector{T}
    objval::T
    change::T
    converged::Bool
    fevals::Int
end

@params struct BESO{T, TO, TC} <: TopOptAlgorithm
    obj::TO
    constr::TC
    vars::Vector{T}
    topology::Vector{T}
    er::T
    maxiter::Int
    p::T
    sens::AbstractVector{T}
    old_sens::AbstractVector{T}
    obj_trace::MVector{<:Any, T}
    tol::T
    sens_tol::T
    result::BESOResult{T}
end

function BESO(obj::Objective{<:ComplianceFunction}, constr::Constraint{<:VolumeFunction}; maxiter = 200, tol = 0.0001, p = 3.0, er = 0.02, sens_tol = tol/100, k = 10)
    T = typeof(obj.comp)
    topology = zeros(T, getncells(obj.problem.ch.dh.grid))
    result = BESOResult(topology, T(NaN), T(NaN), false, 0)
    black = obj.problem.black
    white = obj.problem.white
    nvars = length(topology) - sum(black) - sum(white)
    vars = zeros(T, nvars)
    sens = zeros(T, nvars)
    old_sens = zeros(T, nvars)
    obj_trace = zeros(MVector{k, T})

    return BESO(obj, constr, vars, topology, er, maxiter, p, sens, old_sens, obj_trace, tol, sens_tol, result)
end

update_penalty!(b::BESO, p::Number) = (b.p = p)

function (b::BESO{T, TO, TC})(x0 = copy(b.obj.solver.vars)) where {TO<:Objective{<:ComplianceFunction}, TC<:Constraint{<:VolumeFunction}, T}
    @unpack sens, old_sens, er, tol, maxiter = b
    @unpack obj_trace, topology, sens_tol, vars = b
    @unpack varind, black, white = b.obj.f.problem
    @unpack total_volume, cellvolumes = b.constr.f
    V = b.constr.s
    k = length(obj_trace)

    # Initialize the topology
    for i in 1:length(topology)
        if black[i]
            topology[i] = 1
        elseif white[i]
            topology[i] = 0
        else
            topology[i] = round(x0[varind[i]])
            vars[varind[i]] = topology[i]
        end
    end

    # Calculate the current volume fraction
    true_vol = vol = dot(topology, cellvolumes) / total_volume
    # Main loop
    change = T(1)
    iter = 0
    while (change > tol || true_vol > V) && iter < maxiter
        iter += 1
        if iter > 1
            old_sens .= sens
        end
        vol = max(vol*(1-er), V)
        for j in max(2, k-iter+2):k
            obj_trace[j-1] = obj_trace[j]
        end
        setpenalty!(b.obj, b.p)
        obj_trace[k] = b.obj(vars, sens)
        rmul!(sens, -1)
        if iter > 1
            @. sens = (sens + old_sens) / 2
        end
        l1, l2 = minimum(sens), maximum(sens)
        while (l2 - l1) / l2 > sens_tol
            th = (l1 + l2) / 2
            for i in 1:length(topology)
                if !black[i] && !white[i]
                    topology[i] = T(sign(sens[varind[i]] - th) > 0)
                    vars[varind[i]] = topology[i]
                end
            end
            if dot(topology, cellvolumes) - vol * total_volume > 0
                l1 = th
            else
                l2 = th
            end
        end
        true_vol = dot(topology, cellvolumes) / total_volume
        if iter >= k
            l = sum(@view obj_trace[1:k÷2])
            h = sum(@view obj_trace[k÷2+1:k])
            change = abs(l-h)/h
        end
    end

    objval = obj_trace[k]
    converged = change <= tol
    fevals = iter
    @pack! b.result = change, objval, converged, fevals

    return b.result
end


text:
"""
    Direction

Sweep or contraction direction.
"""
@enum Direction begin
    Left
    Right
end

"""
    flip(dir::Direction)

Flip `dir` (`Left` to `Right` and back).
"""
flip(dir::Direction) = dir == Left ? Right : Left

"""
    realize(x::Number)

Make real, but carefully.
"""
realize(x::Number) = error("Unrecognized numerical type")
realize(x::Real) = x
function realize(x::Complex)
    abs(imag(x)) < 1e-14 || error("Non-zero imaginary component")
    real(x)
end

"""
    S_vn(eigvals::Vector{Float64})

Von Neumann entropy of `eigvals`.
"""
S_vn(eigvals::Vector{Float64}) = -sum(x * log(x) for x in eigvals if x > 0)

"""
    println_result(msg...)

Output a result so that it stands out.
"""
println_result(msg...) = printstyled(msg..., "\n"; color=:green, bold=true)


text:
using FeatherFiles
using DataValues
using IteratorInterfaceExtensions
using TableTraits
using Test
using Arrow
using QueryTables

@testset "FeatherFiles" begin

source = [(Name="John", Age=34., Children=2),
    (Name="Sally", Age=54., Children=1),
    (Name="Jim", Age=34., Children=0)]

output_filename = tempname() * ".feather"

source |> save(output_filename)

try
    sink = load(output_filename) |> IteratorInterfaceExtensions.getiterator |> collect

    @test source == sink

    featherfile = load(output_filename)

    @test IteratorInterfaceExtensions.isiterable(featherfile) == true
    @test TableTraits.isiterabletable(featherfile) == true
    @test TableTraits.supports_get_columns_copy_using_missing(featherfile) == true
    ff_as_cols = TableTraits.get_columns_copy_using_missing(featherfile)
    @test ff_as_cols == (Name=["John", "Sally", "Jim"], Age=[34., 54., 34.], Children=[2,1,0])
finally
    GC.gc()
    GC.gc()
    # rm(output_filename)
end

source2 = [(Name=DataValue("John"), Age=DataValue(34.), Children=DataValue{Int}()),
    (Name=DataValue("Sally"), Age=DataValue{Float64}(), Children=DataValue(1)),
    (Name=DataValue{String}(), Age=DataValue(34.), Children=DataValue(0))]

output_filename2 = tempname() * ".feather"

source2 |> save(output_filename2)

try
    sink2 = load(output_filename2) |> IteratorInterfaceExtensions.getiterator |> collect

    @test source2 == sink2

    featherfile = load(output_filename2)
    @test IteratorInterfaceExtensions.isiterable(featherfile) == true
    @test TableTraits.supports_get_columns_copy_using_missing(featherfile) == true
    ff_as_cols = TableTraits.get_columns_copy_using_missing(featherfile)
    @test isequal(ff_as_cols, (Name=["John", "Sally", missing], Age=[34., missing, 34.], Children=[missing,1,0]))
finally
    GC.gc()
    GC.gc()
    # rm(output_filename2)
end

ar = load(output_filename2)

@test sprint((stream,data)->show(stream, "text/html", data), ar) == "

Name	Age	Children

"John"	34.0	#NA

"Sally"	#NA	1

#NA	34.0	0

"
@test sprint((stream,data)->show(stream, "application/vnd.dataresource+json", data), ar) == "{\"schema\":{\"fields\":[{\"name\":\"Name\",\"type\":\"string\"},{\"name\":\"Age\",\"type\":\"number\"},{\"name\":\"Children\",\"type\":\"integer\"}]},\"data\":[{\"Name\":\"John\",\"Age\":34.0,\"Children\":null},{\"Name\":\"Sally\",\"Age\":null,\"Children\":1},{\"Name\":null,\"Age\":34.0,\"Children\":0}]}"
@test sprint(show, ar) == "3x3 Feather file\nName    │ Age  │ Children\n────────┼──────┼─────────\n\"John\"  │ 34.0 │ #NA     \n\"Sally\" │ #NA  │ 1       \n#NA     │ 34.0 │ 0       "
@test showable("text/html", ar) == true
@test showable("application/vnd.dataresource+json", ar) == true

source3 = DataTable(a=[1,2,NA], b=[3,4,5])

output_filename3 = tempname() * ".feather"

source3 |> save(output_filename3)

ar = load(output_filename3)

@test sprint(show, ar) == "3x2 Feather file\na   │ b\n────┼──\n1   │ 3\n2   │ 4\n#NA │ 5"

end

@testset "Missing Conversion" begin

v1 = FeatherFiles.DataValueArrowVector(NullablePrimitive([2.0, missing, 5.0, 7.0]))
@test getindex(v1, 3) == DataValue{Float64}(5.0)
@test getindex(v1, 2) == DataValue{Float64}()
@test size(v1) == size(v1.data)
@test IndexStyle(v1) == IndexLinear()

v2 = FeatherFiles.DataValueArrowVector(DictEncoding(["fire", "walk", "with", missing, "me"]))
@test getindex(v2, 1) == DataValue{String}("fire")
@test getindex(v2, 4) == DataValue{String}()
@test IndexStyle(v2) == IndexLinear()

v3 = FeatherFiles.MissingDataValueVector([DataValue{Int64}(), DataValue{Int64}(18), DataValue{Int64}(54)])
@test getindex(v3, 2) == 18
@test getindex(v3, 1) === missing
@test IndexStyle(v3) == IndexLinear()

end