text:
# Julia wrapper for header: xkbcommon.h
# Automatically generated using Clang.jl

function xkb_keysym_get_name(keysym, buffer, size)
    ccall((:xkb_keysym_get_name, libxkbcommon), Cint, (xkb_keysym_t, Cstring, Cint), keysym, buffer, size)
end

function xkb_keysym_from_name(name, flags)
    ccall((:xkb_keysym_from_name, libxkbcommon), xkb_keysym_t, (Cstring, xkb_keysym_flags), name, flags)
end

function xkb_keysym_to_utf8(keysym, buffer, size)
    ccall((:xkb_keysym_to_utf8, libxkbcommon), Cint, (xkb_keysym_t, Cstring, Cint), keysym, buffer, size)
end

function xkb_keysym_to_utf32(keysym)
    ccall((:xkb_keysym_to_utf32, libxkbcommon), UInt32, (xkb_keysym_t,), keysym)
end

function xkb_keysym_to_upper(ks)
    ccall((:xkb_keysym_to_upper, libxkbcommon), xkb_keysym_t, (xkb_keysym_t,), ks)
end

function xkb_keysym_to_lower(ks)
    ccall((:xkb_keysym_to_lower, libxkbcommon), xkb_keysym_t, (xkb_keysym_t,), ks)
end

function xkb_context_new(flags)
    ccall((:xkb_context_new, libxkbcommon), Ptr{xkb_context}, (xkb_context_flags,), flags)
end

function xkb_context_ref(context)
    ccall((:xkb_context_ref, libxkbcommon), Ptr{xkb_context}, (Ptr{xkb_context},), context)
end

function xkb_context_unref(context)
    ccall((:xkb_context_unref, libxkbcommon), Cvoid, (Ptr{xkb_context},), context)
end

function xkb_context_set_user_data(context, user_data)
    ccall((:xkb_context_set_user_data, libxkbcommon), Cvoid, (Ptr{xkb_context}, Ptr{Cvoid}), context, user_data)
end

function xkb_context_get_user_data(context)
    ccall((:xkb_context_get_user_data, libxkbcommon), Ptr{Cvoid}, (Ptr{xkb_context},), context)
end

function xkb_context_include_path_append(context, path)
    ccall((:xkb_context_include_path_append, libxkbcommon), Cint, (Ptr{xkb_context}, Cstring), context, path)
end

function xkb_context_include_path_append_default(context)
    ccall((:xkb_context_include_path_append_default, libxkbcommon), Cint, (Ptr{xkb_context},), context)
end

function xkb_context_include_path_reset_defaults(context)
    ccall((:xkb_context_include_path_reset_defaults, libxkbcommon), Cint, (Ptr{xkb_context},), context)
end

function xkb_context_include_path_clear(context)
    ccall((:xkb_context_include_path_clear, libxkbcommon), Cvoid, (Ptr{xkb_context},), context)
end

function xkb_context_num_include_paths(context)
    ccall((:xkb_context_num_include_paths, libxkbcommon), UInt32, (Ptr{xkb_context},), context)
end

function xkb_context_include_path_get(context, index)
    ccall((:xkb_context_include_path_get, libxkbcommon), Cstring, (Ptr{xkb_context}, UInt32), context, index)
end

function xkb_context_set_log_level(context, level)
    ccall((:xkb_context_set_log_level, libxkbcommon), Cvoid, (Ptr{xkb_context}, xkb_log_level), context, level)
end

function xkb_context_get_log_level(context)
    ccall((:xkb_context_get_log_level, libxkbcommon), xkb_log_level, (Ptr{xkb_context},), context)
end

function xkb_context_set_log_verbosity(context, verbosity)
    ccall((:xkb_context_set_log_verbosity, libxkbcommon), Cvoid, (Ptr{xkb_context}, Cint), context, verbosity)
end

function xkb_context_get_log_verbosity(context)
    ccall((:xkb_context_get_log_verbosity, libxkbcommon), Cint, (Ptr{xkb_context},), context)
end

function xkb_context_set_log_fn(context, log_fn)
    ccall((:xkb_context_set_log_fn, libxkbcommon), Cvoid, (Ptr{xkb_context}, Ptr{Cvoid}), context, log_fn)
end

function xkb_keymap_new_from_names(context, names, flags)
    ccall((:xkb_keymap_new_from_names, libxkbcommon), Ptr{xkb_keymap}, (Ptr{xkb_context}, Ptr{xkb_rule_names}, xkb_keymap_compile_flags), context, names, flags)
end

function xkb_keymap_new_from_file(context, file, format, flags)
    ccall((:xkb_keymap_new_from_file, libxkbcommon), Ptr{xkb_keymap}, (Ptr{xkb_context}, Ptr{FILE}, xkb_keymap_format, xkb_keymap_compile_flags), context, file, format, flags)
end

function xkb_keymap_new_from_string(context, string, format, flags)
    ccall((:xkb_keymap_new_from_string, libxkbcommon), Ptr{xkb_keymap}, (Ptr{xkb_context}, Cstring, xkb_keymap_format, xkb_keymap_compile_flags), context, string, format, flags)
end

function xkb_keymap_new_from_buffer(context, buffer, length, format, flags)
    ccall((:xkb_keymap_new_from_buffer, libxkbcommon), Ptr{xkb_keymap}, (Ptr{xkb_context}, Cstring, Cint, xkb_keymap_format, xkb_keymap_compile_flags), context, buffer, length, format, flags)
end

function xkb_keymap_ref(keymap)
    ccall((:xkb_keymap_ref, libxkbcommon), Ptr{xkb_keymap}, (Ptr{xkb_keymap},), keymap)
end

function xkb_keymap_unref(keymap)
    ccall((:xkb_keymap_unref, libxkbcommon), Cvoid, (Ptr{xkb_keymap},), keymap)
end

function xkb_keymap_get_as_string(keymap, format)
    ccall((:xkb_keymap_get_as_string, libxkbcommon), Cstring, (Ptr{xkb_keymap}, xkb_keymap_format), keymap, format)
end

function xkb_keymap_min_keycode(keymap)
    ccall((:xkb_keymap_min_keycode, libxkbcommon), xkb_keycode_t, (Ptr{xkb_keymap},), keymap)
end

function xkb_keymap_max_keycode(keymap)
    ccall((:xkb_keymap_max_keycode, libxkbcommon), xkb_keycode_t, (Ptr{xkb_keymap},), keymap)
end

function xkb_keymap_key_for_each(keymap, iter, data)
    ccall((:xkb_keymap_key_for_each, libxkbcommon), Cvoid, (Ptr{xkb_keymap}, xkb_keymap_key_iter_t, Ptr{Cvoid}), keymap, iter, data)
end

function xkb_keymap_key_get_name(keymap, key)
    ccall((:xkb_keymap_key_get_name, libxkbcommon), Cstring, (Ptr{xkb_keymap}, xkb_keycode_t), keymap, key)
end

function xkb_keymap_key_by_name(keymap, name)
    ccall((:xkb_keymap_key_by_name, libxkbcommon), xkb_keycode_t, (Ptr{xkb_keymap}, Cstring), keymap, name)
end

function xkb_keymap_num_mods(keymap)
    ccall((:xkb_keymap_num_mods, libxkbcommon), xkb_mod_index_t, (Ptr{xkb_keymap},), keymap)
end

function xkb_keymap_mod_get_name(keymap, idx)
    ccall((:xkb_keymap_mod_get_name, libxkbcommon), Cstring, (Ptr{xkb_keymap}, xkb_mod_index_t), keymap, idx)
end

function xkb_keymap_mod_get_index(keymap, name)
    ccall((:xkb_keymap_mod_get_index, libxkbcommon), xkb_mod_index_t, (Ptr{xkb_keymap}, Cstring), keymap, name)
end

function xkb_keymap_num_layouts(keymap)
    ccall((:xkb_keymap_num_layouts, libxkbcommon), xkb_layout_index_t, (Ptr{xkb_keymap},), keymap)
end

function xkb_keymap_layout_get_name(keymap, idx)
    ccall((:xkb_keymap_layout_get_name, libxkbcommon), Cstring, (Ptr{xkb_keymap}, xkb_layout_index_t), keymap, idx)
end

function xkb_keymap_layout_get_index(keymap, name)
    ccall((:xkb_keymap_layout_get_index, libxkbcommon), xkb_layout_index_t, (Ptr{xkb_keymap}, Cstring), keymap, name)
end

function xkb_keymap_num_leds(keymap)
    ccall((:xkb_keymap_num_leds, libxkbcommon), xkb_led_index_t, (Ptr{xkb_keymap},), keymap)
end

function xkb_keymap_led_get_name(keymap, idx)
    ccall((:xkb_keymap_led_get_name, libxkbcommon), Cstring, (Ptr{xkb_keymap}, xkb_led_index_t), keymap, idx)
end

function xkb_keymap_led_get_index(keymap, name)
    ccall((:xkb_keymap_led_get_index, libxkbcommon), xkb_led_index_t, (Ptr{xkb_keymap}, Cstring), keymap, name)
end

function xkb_keymap_num_layouts_for_key(keymap, key)
    ccall((:xkb_keymap_num_layouts_for_key, libxkbcommon), xkb_layout_index_t, (Ptr{xkb_keymap}, xkb_keycode_t), keymap, key)
end

function xkb_keymap_num_levels_for_key(keymap, key, layout)
    ccall((:xkb_keymap_num_levels_for_key, libxkbcommon), xkb_level_index_t, (Ptr{xkb_keymap}, xkb_keycode_t, xkb_layout_index_t), keymap, key, layout)
end

function xkb_keymap_key_get_syms_by_level(keymap, key, layout, level, syms_out)
    ccall((:xkb_keymap_key_get_syms_by_level, libxkbcommon), Cint, (Ptr{xkb_keymap}, xkb_keycode_t, xkb_layout_index_t, xkb_level_index_t, Ptr{Ptr{xkb_keysym_t}}), keymap, key, layout, level, syms_out)
end

function xkb_keymap_key_repeats(keymap, key)
    ccall((:xkb_keymap_key_repeats, libxkbcommon), Cint, (Ptr{xkb_keymap}, xkb_keycode_t), keymap, key)
end

function xkb_state_new(keymap)
    ccall((:xkb_state_new, libxkbcommon), Ptr{xkb_state}, (Ptr{xkb_keymap},), keymap)
end

function xkb_state_ref(state)
    ccall((:xkb_state_ref, libxkbcommon), Ptr{xkb_state}, (Ptr{xkb_state},), state)
end

function xkb_state_unref(state)
    ccall((:xkb_state_unref, libxkbcommon), Cvoid, (Ptr{xkb_state},), state)
end

function xkb_state_get_keymap(state)
    ccall((:xkb_state_get_keymap, libxkbcommon), Ptr{xkb_keymap}, (Ptr{xkb_state},), state)
end

function xkb_state_update_key(state, key, direction)
    ccall((:xkb_state_update_key, libxkbcommon), xkb_state_component, (Ptr{xkb_state}, xkb_keycode_t, xkb_key_direction), state, key, direction)
end

function xkb_state_update_mask(state, depressed_mods, latched_mods, locked_mods, depressed_layout, latched_layout, locked_layout)
    ccall((:xkb_state_update_mask, libxkbcommon), xkb_state_component, (Ptr{xkb_state}, xkb_mod_mask_t, xkb_mod_mask_t, xkb_mod_mask_t, xkb_layout_index_t, xkb_layout_index_t, xkb_layout_index_t), state, depressed_mods, latched_mods, locked_mods, depressed_layout, latched_layout, locked_layout)
end

function xkb_state_key_get_syms(state, key, syms_out)
    ccall((:xkb_state_key_get_syms, libxkbcommon), Cint, (Ptr{xkb_state}, xkb_keycode_t, Ptr{Ptr{xkb_keysym_t}}), state, key, syms_out)
end

function xkb_state_key_get_utf8(state, key, buffer, size)
    ccall((:xkb_state_key_get_utf8, libxkbcommon), Cint, (Ptr{xkb_state}, xkb_keycode_t, Cstring, Cint), state, key, buffer, size)
end

function xkb_state_key_get_utf32(state, key)
    ccall((:xkb_state_key_get_utf32, libxkbcommon), UInt32, (Ptr{xkb_state}, xkb_keycode_t), state, key)
end

function xkb_state_key_get_one_sym(state, key)
    ccall((:xkb_state_key_get_one_sym, libxkbcommon), xkb_keysym_t, (Ptr{xkb_state}, xkb_keycode_t), state, key)
end

function xkb_state_key_get_layout(state, key)
    ccall((:xkb_state_key_get_layout, libxkbcommon), xkb_layout_index_t, (Ptr{xkb_state}, xkb_keycode_t), state, key)
end

function xkb_state_key_get_level(state, key, layout)
    ccall((:xkb_state_key_get_level, libxkbcommon), xkb_level_index_t, (Ptr{xkb_state}, xkb_keycode_t, xkb_layout_index_t), state, key, layout)
end

function xkb_state_serialize_mods(state, components)
    ccall((:xkb_state_serialize_mods, libxkbcommon), xkb_mod_mask_t, (Ptr{xkb_state}, xkb_state_component), state, components)
end

function xkb_state_serialize_layout(state, components)
    ccall((:xkb_state_serialize_layout, libxkbcommon), xkb_layout_index_t, (Ptr{xkb_state}, xkb_state_component), state, components)
end

function xkb_state_mod_name_is_active(state, name, type)
    ccall((:xkb_state_mod_name_is_active, libxkbcommon), Cint, (Ptr{xkb_state}, Cstring, xkb_state_component), state, name, type)
end

function xkb_state_mod_index_is_active(state, idx, type)
    ccall((:xkb_state_mod_index_is_active, libxkbcommon), Cint, (Ptr{xkb_state}, xkb_mod_index_t, xkb_state_component), state, idx, type)
end

function xkb_state_key_get_consumed_mods2(state, key, mode)
    ccall((:xkb_state_key_get_consumed_mods2, libxkbcommon), xkb_mod_mask_t, (Ptr{xkb_state}, xkb_keycode_t, xkb_consumed_mode), state, key, mode)
end

function xkb_state_key_get_consumed_mods(state, key)
    ccall((:xkb_state_key_get_consumed_mods, libxkbcommon), xkb_mod_mask_t, (Ptr{xkb_state}, xkb_keycode_t), state, key)
end

function xkb_state_mod_index_is_consumed2(state, key, idx, mode)
    ccall((:xkb_state_mod_index_is_consumed2, libxkbcommon), Cint, (Ptr{xkb_state}, xkb_keycode_t, xkb_mod_index_t, xkb_consumed_mode), state, key, idx, mode)
end

function xkb_state_mod_index_is_consumed(state, key, idx)
    ccall((:xkb_state_mod_index_is_consumed, libxkbcommon), Cint, (Ptr{xkb_state}, xkb_keycode_t, xkb_mod_index_t), state, key, idx)
end

function xkb_state_mod_mask_remove_consumed(state, key, mask)
    ccall((:xkb_state_mod_mask_remove_consumed, libxkbcommon), xkb_mod_mask_t, (Ptr{xkb_state}, xkb_keycode_t, xkb_mod_mask_t), state, key, mask)
end

function xkb_state_layout_name_is_active(state, name, type)
    ccall((:xkb_state_layout_name_is_active, libxkbcommon), Cint, (Ptr{xkb_state}, Cstring, xkb_state_component), state, name, type)
end

function xkb_state_layout_index_is_active(state, idx, type)
    ccall((:xkb_state_layout_index_is_active, libxkbcommon), Cint, (Ptr{xkb_state}, xkb_layout_index_t, xkb_state_component), state, idx, type)
end

function xkb_state_led_name_is_active(state, name)
    ccall((:xkb_state_led_name_is_active, libxkbcommon), Cint, (Ptr{xkb_state}, Cstring), state, name)
end

function xkb_state_led_index_is_active(state, idx)
    ccall((:xkb_state_led_index_is_active, libxkbcommon), Cint, (Ptr{xkb_state}, xkb_led_index_t), state, idx)
end
# Julia wrapper for header: xkbcommon-x11.h
# Automatically generated using Clang.jl

function xkb_x11_setup_xkb_extension(connection, major_xkb_version, minor_xkb_version, flags, major_xkb_version_out, minor_xkb_version_out, base_event_out, base_error_out)
    ccall((:xkb_x11_setup_xkb_extension, libxkbcommon_x11), Cint, (Ptr{xcb_connection_t}, UInt16, UInt16, xkb_x11_setup_xkb_extension_flags, Ptr{UInt16}, Ptr{UInt16}, Ptr{UInt8}, Ptr{UInt8}), connection, major_xkb_version, minor_xkb_version, flags, major_xkb_version_out, minor_xkb_version_out, base_event_out, base_error_out)
end

function xkb_x11_get_core_keyboard_device_id(connection)
    ccall((:xkb_x11_get_core_keyboard_device_id, libxkbcommon_x11), Int32, (Ptr{xcb_connection_t},), connection)
end

function xkb_x11_keymap_new_from_device(context, connection, device_id, flags)
    ccall((:xkb_x11_keymap_new_from_device, libxkbcommon_x11), Ptr{xkb_keymap}, (Ptr{xkb_context}, Ptr{xcb_connection_t}, Int32, xkb_keymap_compile_flags), context, connection, device_id, flags)
end

function xkb_x11_state_new_from_device(keymap, connection, device_id)
    ccall((:xkb_x11_state_new_from_device, libxkbcommon_x11), Ptr{xkb_state}, (Ptr{xkb_keymap}, Ptr{xcb_connection_t}, Int32), keymap, connection, device_id)
end


text:
# save to disk/read from disk
printstyled("  Native file I/O\n", color=:light_green)
savfile1 = "test.dat"
fs = 100.0
nc = 10
nx = 2^12

D = randSeisData(nc, nx=nx, s=1.0)
t = mk_t(nx, D.t[1][1,2])
for i in 1:D.n
  D.fs[i] = fs
  D.t[i] = copy(t)
end
S = convert(NodalData, D)
S.ox = (rand()-0.5) * 360.0
S.oy = (rand()-0.5) * 90.0
S.oz = rand()*1000.0
S.info["foo"] = "bar"
pop_nodal_dict!(S.info)

printstyled("    NodalChannel\n", color=:light_green)
C = getindex(S, 1)
wseis(savfile1, C)
R = rseis(savfile1)[1]
@test R == C

printstyled("    NodalData\n", color=:light_green)
wseis(savfile1, S)
R = rseis(savfile1)[1]
@test R == S
rm(savfile1)

ext: jl


text:
push!(LOAD_PATH, pwd())
using Distributions
using PyPlot
using BVB

max_stars_repo_name: richardkwo/boosting-VB

max_issues_repo_name: richardkwo/boosting-VB

max_forks_repo_name: richardkwo/boosting-VB


text:
"""
`AbstractPolygon{N}`

Supertype for `N`-vertices polygon.
"""
abstract type AbstractPolygon{N} <: AbstractVector{N} end

Base.haslength(::Type{AbstractPolygon}) = false
Base.haslength(::Type{<:AbstractPolygon{N}}) where N = true
Base.length(::Type{<:AbstractPolygon{N}}) where N = N

Base.haslength(::AbstractPolygon{N}) where N = true
Base.length(::AbstractPolygon{N}) where N = N
Base.size(::AbstractPolygon{N}) where N = (N,)
Base.eltype(::AbstractPolygon) = UInt32
Base.IndexStyle(::Type{<:AbstractPolygon}) = Base.IndexLinear()

"""
`Polygon{N} <: AbstractPolygon{N}`

Fixedsize polygon with `N`-vertices representing offset into zeroindexed vertex array.
"""
struct Polygon{N} <: AbstractPolygon{N}
  data :: NTuple{N,UInt32}
end

Polygon{L}(array::Array) where L = Polygon{L}((array...,))
Base.getindex(polygon::Polygon, i) = polygon.data[i]
function Base.read(io::IO, ::Type{Polygon{N}}) where N
  data = read(io, N * sizeof(UInt32))
  Base.unsafe_load(Base.unsafe_convert(Ptr{Polygon{N}}, data))
end

"""
`Triangle <: AbstractPolygon{3}`

Alias for `Polygon{3}`
"""
const Triangle = Polygon{3}
Triangle(a,b,c) = Polygon{3}((a,b,c))

edges(poly::Polygon{N}) where N = [(poly[i], poly[i%N+1]) for i in 1:N]


text:
println("\n\n\nStarting runtests.jl $(join(ARGS, " ")) ...")
using FactCheck, KShiftsClustering

getdata(n) = 10*rand(1, n) + 0.5
facts("all") do
    centers = kshifts(getdata(1_000_000), 10)
    @fact round(Int, sort(vec(centers)))  -->  collect(1:10)

    centers = kshifts(getdata(1000), 10)
    for i = 1:100
        kshifts!(centers, getdata(1000))
    end
    @fact round(Int, sort(vec(centers)))  -->  collect(1:10)

    centers, ids = kshiftmedoids(rand(10, 100000), 10)

    centers = kshifts(rand(2,1000), 10)
end


text:
module FirthRegression
using Reexport, Compat
@reexport using GLM
using GLM: wrkwt!, installbeta!
export FirthGLM, penalized_deviance, anova

"""
Second derivative of mu with respect to eta.
"""
function mueta2 end
mueta2(::LogitLink, η) = (e = exp(η); (e*(1-e))/(1+e)^3)
mueta2(::LogLink, η) = exp(η)
# Need to figure out how to compute the deviance for non-canonical links if we want to use them
# mueta2(::CauchitLink, η) = -2*η/(pi*(1+η^2)^2)
# mueta2(::ProbitLink, η) = -exp(-η^2/2)*η/sqrt2π
# mueta2(::InverseLink, η) = 2/η^3

@compat typealias FirthGlmResp{T,D<:Union{Binomial,Poisson},L} GlmResp{T,D,L}

type FirthGLM{G<:FirthGlmResp,L<:LinPred} <: GLM.AbstractGLM
    rr::G
    pp::L
    fit::Bool
end

type FirthGLMTest{G<:FirthGlmResp,L<:LinPred} <: GLM.AbstractGLM
    rr::G
    pp::L
    pp_full::L
    fit::Bool
end

"""
Update just the factorization of the X'WX, without computing new betas.
"""
function updatefact end
function updatefact!{T}(p::GLM.DensePredChol{T}, wt::Vector{T})
    scr = scale!(p.scratch, wt, p.X)
    cholfact!(At_mul_B!(GLM.cholfactors(p.chol), scr, p.X), :U)
end
@eval function updatefact!{T}(p::GLM.SparsePredChol{T}, wt::Vector{T})
    scr = scale!(p.scratch, wt, p.X)
    XtX = p.Xt*scr
    p.chol = $(if VERSION >= v"0.4.0-dev+3307"
        :(cholfact(Symmetric{eltype(XtX),typeof(XtX)}(XtX, 'U')))
    else
        :(cholfact(scale!(XtX + XtX', convert(eltype(XtX), 1/2))))
    end)
end

updatefact!(m::FirthGLM) = updatefact!(m.pp, wrkwt!(m.rr))
function updatefact!(m::FirthGLMTest)
    wt = wrkwt!(m.rr)
    updatefact!(m.pp, wt)
    updatefact!(m.pp_full, wt)
end

"""
Compute delbeta, without recomputing the factorization of X'WX.
"""
function updatedelbeta end
function updatedelbeta!{T}(p::GLM.DensePredChol{T}, r::Vector{T})
    A_ldiv_B!(p.chol, At_mul_B!(p.delbeta, p.scratch, r))
    p
end
function updatedelbeta!{T}(p::GLM.SparsePredChol{T}, r::Vector{T})
    p.delbeta = p.chol\Ac_mul_B!(p.delbeta, p.scratch, r)
end

"""
Compute the diagonal of the hat matrix.
"""
function hatdiag end
hatdiag(p::GLM.DensePredChol, wt::Vector) = vec(sumabs2(p.X/p.chol[:U], 2)).*wt
hatdiag(p::GLM.SparsePredChol, wt::Vector) = vec(sumabs2(p.chol[:L]\p.Xt[p.chol[:p], :], 1)).*wt

hatdiag(m::FirthGLM) = hatdiag(m.pp, m.rr.wrkwts)
hatdiag(m::FirthGLMTest) = hatdiag(m.pp_full, m.rr.wrkwts)

"""
Adjust working residuals using the given hat matrix diagonal `h`.
"""
function adjustwrkresid!(m::@compat(Union{FirthGLM,FirthGLMTest}), h::Vector)
    r = m.rr
    link = r.l
    eta = r.eta
    mueta = r.mueta
    var = r.var
    wrkresid = r.wrkresid
    @inbounds @simd for i = 1:length(eta)
        # See Kosmidis and Firth, 2009, 4⋅3
        wrkresid[i] += h[i]*mueta2(link, eta[i])/(2*mueta[i]^3/var[i])
    end
    r
end

"""
Compute the penalized deviance. This requires that updatefact! has been
called since the working weights were last updated.
"""
function penalized_deviance end
penalized_deviance(m::FirthGLM) = deviance(m.rr) - logdet(m.pp.chol)
penalized_deviance(m::FirthGLMTest) = deviance(m.rr) - logdet(m.pp_full.chol)

function GLM._fit!(m::@compat(Union{FirthGLM,FirthGLMTest}), verbose::Bool, maxIter::Integer, minStepFac::Real,
                  convTol::Real, start)
    m.fit && return m
    maxIter >= 1 || error("maxIter must be positive")
    0 < minStepFac < 1 || error("minStepFac must be in (0, 1)")

    cvg = false; p = m.pp; r = m.rr
    lp = r.mu
    if start != nothing
        copy!(p.beta0, start)
        fill!(p.delbeta, 0)
    else
        delbeta!(p, wrkresp(r), wrkwt!(r))
        installbeta!(p)
    end
    devold = Inf
    for i=1:maxIter
        f = 1.0
        local dev
        linpred!(lp, p)
        try
            updatemu!(r, lp)
            updatefact!(m)
            dev = penalized_deviance(m)
        catch e
            isa(e, DomainError) ? (dev = Inf) : rethrow(e)
        end
        dev == NaN && (dev = Inf)

        while dev > devold
            f /= 2.; f > minStepFac || error("step-halving failed at beta0 = $(p.beta0)")
            verbose && println("step-halving $f: $dev")
            try
                updatemu!(r, linpred(p, f))
                updatefact!(m)
                dev = penalized_deviance(m)
            catch e
                isa(e, DomainError) ? (dev = Inf) : rethrow(e)
            end
            dev == NaN && (dev = Inf)
        end

        installbeta!(p, f)

        crit = (devold - dev)/dev
        verbose && println("$i: $dev, $crit")
        if crit < convTol; cvg = true; break end

        devold = dev
        adjustwrkresid!(m, hatdiag(m))
        updatedelbeta!(p, r.wrkresid)
    end
    cvg || error("failure to converge in $maxIter iterations")
    m.fit = true
    m
end

cols(X::AbstractMatrix) = [sub(X, :, i) for i = 1:size(X, 2)]
cols(X::SparseMatrixCSC) = [X[:, i] for i = 1:size(X, 2)]

immutable ANOVAResult
    chisq::Float64
    df::Int
    p::Float64
end

"""
Performs analysis of deviance for two nested design matrices with
a Firth penalized approach. The penalty for both models is determined
based on the Fisher information matrix for the larger model.
"""
function anova(m1::FirthGLM, m2::FirthGLM; fitargs...)
    m1.rr.y == m2.rr.y || throw(ArgumentError("model responses (y) do not match"))
    (m1.rr.d == m2.rr.d && m1.rr.l == m2.rr.l) ||
        throw(ArgumentError("model distributions or link functions do not match"))
    m1.rr.wts == m2.rr.wts || throw(ArgumentError("model weights do not match"))
    m1.rr.offset == m2.rr.offset || throw(ArgumentError("model offsets do not match"))

    if size(m1.pp.X, 2) == size(m2.pp.X, 2)
        throw(ArgumentError("one design matrix must be nested inside the other"))
    elseif size(m1.pp.X, 2) > size(m2.pp.X, 2)
        full, reduced = m1, m2
    else
        full, reduced = m2, m1
    end

    # Figure out which columns are nested
    fullcols = cols(full.pp.X)
    reduced_colidx = findin(fullcols, cols(reduced.pp.X))
    if length(reduced_colidx) != size(reduced.pp.X, 2)
        throw(ArgumentError("one design matrix must be nested inside the other"))
    end

    # Fit the full model if necessary
    if !full.fit
        fit!(full; fitargs...)
    end
    fulldev = penalized_deviance(full)

    # Fit the reduced model
    if !reduced.fit
        # If the second model has not been fit, we can reuse its
        # response and linear predictor
        rr = reduced.rr
        pp = reduced.pp
        startcoef = coef(full)[reduced_colidx]
    else
        rr = typeof(reduced.rr)(reduced.rr.y, reduced.rr.d, reduced.rr.l,
                                copy(reduced.rr.eta), similar(reduced.rr.y),
                                reduced.rr.offset, reduced.rr.wts)
        pp = GLM.cholpred(reduced.pp.X)
        startcoef = coef(reduced)
    end
    reducedmodel = FirthGLMTest(rr, pp, GLM.cholpred(full.pp.X), false)
    fit!(reducedmodel; start=startcoef, fitargs...)
    reduceddev = penalized_deviance(reducedmodel)

    chisq = reduceddev - fulldev
    df = size(full.pp.X, 2) - length(reduced_colidx)
    ANOVAResult(chisq, df, ccdf(Chisq(df), chisq))
end

end # module


text:
module PETScX

using MPI

MPI.Initialized() || MPI.Init()

using Libdl

include("const.jl")
include("startup.jl")
include("lib.jl")
include("init.jl")
include("viewer.jl")
include("options.jl")
include("vec.jl")
include("mat.jl")
# include("matshell.jl")
# include("dm.jl")
# include("dmda.jl")
# include("ksp.jl")
# include("ref.jl")
# include("pc.jl")
# include("snes.jl")
# include("sys.jl")

end


text:
#
#Simple example of a Poisson problem with a variable coefficient:
#
#           y
#             |
#             |           u_y = 0
#          Ly ----------------------------------
#             |                                |
#             |                                |
#             |                                |
#       u = 0 |   -(au_x)_x - (au_y)_y = f     | u = 0
#             |                                |
#             |                                |
#             |                                |
#           -----------------------------------|-----  x
#                         u_y = 0               Lx
#
#Here,
#
#      u = sin(ωx x) * cos(ωy y)
#      a = exp(x-y)
#      f = ω exp(x-y) ( 2ω sin(ω x) cos(ω y) - cos ω(x-y) )
#      ωx = mπ / Lx
#      ωy = nπ / Ly
#

using FinElt
using FinElt.PlanarPoisson

include("params.jl")
const ωx = pi / Lx
const ωy = pi / Ly

function exact_u(x)
    return sin(ωx*x[1]) * cos(ωy*x[2])
end

function a(x)
    return exp(x[1]-x[2])
end

function f(x)
    return exp(x[1]-x[2]) * (
             ωx * ( ωx * sin(ωx*x[1]) - cos(ωx*x[1]) ) * cos(ωy*x[2])
           + ωy * ( ωy * cos(ωy*x[2]) - sin(ωy*x[2]) ) * sin(ωx*x[1])
           )
end

ell_f(v, z) = P1SourceTimesFunc!(v, z, f)
lin_functionals = [ ("Omega", ell_f) ]

maxerr = zeros(refinements+1)
@printf("%10s  %12s  %8s  %8s\n\n",
        "N", "max error", "rate", "seconds")
for k = 0:refinements
    start = time()
    mesh = read_msh_file("rect$k.msh")
    vp = VariationalProblem(mesh, ["Left", "Right"])
    add_bilin_form!(vp, "Omega", grad_dot_grad!, a)
    add_lin_functnl!(vp, "Omega", source_times_func!, f)
    # Implicit zero boundary conditions.
    A, b = assembled_linear_system(vp)
    ufree = A \ b
    uh = complete_soln(ufree, vp)
    u = get_nodal_vals(exact_u, mesh)
    finish = time()
    maxerr[k+1] = maximum(abs.(uh-u))
    N = length(ufree)
    if k == 0
        @printf("%10d  %12.4e\n", N, maxerr[k+1])
    else
        rate = log2(maxerr[k]/maxerr[k+1])
        elapsed = finish - start
        @printf("%10d  %12.4e  %8.4f  %8.4f\n",
                N, maxerr[k+1], rate, elapsed)
    end
end


text:
Tables.isrowtable(::Type{<:DESolution}) = true
Tables.columns(x::DESolution) = Tables.columntable(Tables.rows(x))

struct DESolutionRows{T, U}
    names::Vector{Symbol}
    types::Vector{Type}
    lookup::Dict{Symbol, Int}
    t::T
    u::U
end

DESolutionRows(names, types, t, u) = DESolutionRows(names, types, Dict(nm => i for (i, nm) in enumerate(names)), t, u)

Base.length(x::DESolutionRows) = length(x.u)
Base.eltype(::Type{DESolutionRows{T, U}}) where {T, U} = DESolutionRow{eltype(T), eltype(U)}
Base.iterate(x::DESolutionRows, st=1) = st > length(x) ? nothing : (DESolutionRow(x.names, x.lookup, x.t[st], x.u[st]), st + 1)

function Tables.rows(sol::DESolution)
    VT = eltype(sol.u)
    if VT <: AbstractArray
        N = length(sol.u[1])
        names = [:timestamp, (has_syms(sol.prob.f) ? (sol.prob.f.syms[i] for i = 1:N) : (Symbol("value", i) for i = 1:N))...]
        types = Type[eltype(sol.t), (eltype(sol.u[1]) for i = 1:N)...]
    else
        names = [:timestamp, has_syms(sol.prob.f) ? sol.prob.f.syms[1] : :value]
        types = Type[eltype(sol.t), VT]
    end
    return DESolutionRows(names, types, sol.t, sol.u)
end

Tables.schema(x::DESolutionRows) = Tables.Schema(x.names, x.types)

struct DESolutionRow{T, U} <: Tables.AbstractRow
    names::Vector{Symbol}
    lookup::Dict{Symbol, Int}
    t::T
    u::U
end

Tables.columnnames(x::DESolutionRow) = getfield(x, :names)
Tables.getcolumn(x::DESolutionRow, i::Int) = i == 1 ? getfield(x, :t) : getfield(x, :u)[i - 1]
Tables.getcolumn(x::DESolutionRow, nm::Symbol) = nm === :timestamp ? getfield(x, :t) : getfield(x, :u)[getfield(x, :lookup)[nm] - 1]

IteratorInterfaceExtensions.isiterable(sol::DESolution) = true
IteratorInterfaceExtensions.getiterator(sol::DESolution) =
    Tables.datavaluerows(Tables.rows(sol))
#TableTraits.isiterabletable(sol::DESolution) = true
Tables.istable(::EnsembleSolution) = false


text:
import Base.BLAS.axpy!
using Base.Distributed

function inplace_train_vectors!(vm::VectorModel, doc::DenseArray{Tw},
		window_length::Int,
		start_lr::Float64, total_words::Float64, words_read::DenseArray{Int64},
		total_ll::DenseArray{Float64}; batch::Int=10000,
		context_cut::Bool = true, sense_treshold::Float64=1e-32)

	N = length(doc)
	in_grad = zeros(Tsf, M(vm), T(vm))
	out_grad = zeros(Tsf, M(vm))

	z = zeros(T(vm))
	senses = 0.
	max_senses = 0.

	tic()
	for i in 1:N
		x = doc[i]
		lr1 = max(start_lr * (1 - words_read[1] / (total_words+1)), start_lr * 1e-4)
		lr2 = lr1

		random_reduce = context_cut ? rand(1:window_length-1) : 0
		window = window_length - random_reduce

		z[:] = 0.

		n_senses = var_init_z!(vm, x, z)
		senses += n_senses
		max_senses = max(max_senses, n_senses)
		for j in max(1, i - window):min(N, i + window) #
			if i == j continue end
			var_update_z!(vm, x, doc[j], z)
		end

		exp_normalize!(z)

		for j in max(1, i - window):min(N, i + window)
			if i == j continue end
			y = doc[j]

			ll = in_place_update!(vm, x, y, z, lr1, in_grad, out_grad, sense_treshold)

			total_ll[2] += 1
			total_ll[1] += (ll - total_ll[1]) / total_ll[2]

		end

		words_read[1] += 1

		#variational update for q(pi_v)
		var_update_counts!(vm, x, z, lr2)

		if i % batch == 0
			time_per_kword = batch / toq() / 1000
			@printf("%.2f%% %.4f %.4f %.4f %.2f/%.2f %.2f kwords/sec\n",
					words_read[1] / (total_words / 100),
					total_ll[1], lr1, lr2, senses / i, max_senses, time_per_kword)
			tic()
		end

		if words_read[1] > total_words break end
	end
	toq()
end

function in_place_update!{Tw <: Integer}(vm::VectorModel,
		x::Tw, y::Tw, z::DenseArray{Float64}, lr::Float64,
		in_grad::DenseArray{Tsf, 2}, out_grad::DenseArray{Tsf}, sense_treshold::Float64)

	return ccall((:inplace_update, "superlib"), Float32,
		(Ptr{Float32}, Ptr{Float32},
			Int, Int, Ptr{Float64},
			Int,
			Ptr{Int32}, Ptr{Int8}, Int64,
			Ptr{Float32}, Ptr{Float32},
			Float32, Float32),
		sdata(vm.In), sdata(vm.Out),
			M(vm), T(vm), z,
			x,
			view(vm.path, :, y), view(vm.code, :, y), size(vm.code, 1),
			in_grad, out_grad,
			lr, sense_treshold)
end

function var_init_z!(vm::VectorModel, x::Integer, z::DenseArray{Float64})
	return expected_logpi!(z, vm, x)
end

function var_update_z!{Tw <: Integer}(vm::VectorModel,
		x::Tw, y::Tw, z::DenseArray{Float64}, num_meanings::Int=T(vm))
	ccall((:update_z, "superlib"), Void,
		(Ptr{Float32}, Ptr{Float32},
			Int, Int, Ptr{Float64}, Int,
			Ptr{Int32}, Ptr{Int8}, Int64),
		vm.In, vm.Out,
			M(vm), num_meanings, z, x,
			view(vm.path, :, y), view(vm.code, :, y), size(vm.path, 1))
end

function var_update_counts!(vm::VectorModel, x::Integer,
		local_counts::DenseArray{Float64}, lr::Float64)
	counts = view(vm.counts, :, x)
	for k in 1:T(vm)
		counts[k] += lr * (local_counts[k] * vm.frequencies[x] - counts[k])
	end
end

function inplace_train_vectors!(vm::VectorModel, dict::Dictionary, path::AbstractString,
		window_length::Int; batch::Int = 64000, start_lr::Float64 = 0.025,
		log_path::Union{AbstractString, Void} = nothing, threshold::Float64 = Inf,
		context_cut::Bool = true, epochs::Int = 1, init_count::Float64=-1, sense_treshold::Float64=1e-32)
	for w in 1:V(vm)
		vm.counts[1, w] = init_count > 0 ? init_count : vm.frequencies[w]
	end

	nbytes = filesize(path)
	train_words = Float64(sum(vm.frequencies)) * epochs

	words_read = shared_zeros(Int64, (1,))
	total_ll = shared_zeros(Float64, (2,))

	function do_work(id::Int)
		file = open(path)

		bytes_per_worker = convert(Int, floor(nbytes / nworkers()))

		start_pos = bytes_per_worker * (id-1)
		end_pos = start_pos+bytes_per_worker

		seek(file, start_pos)
		align(file)
		buffer = zeros(Int32, batch)
		while words_read[1] < train_words
			doc = read_words(file, start_pos, end_pos, dict, buffer,
				vm.frequencies, threshold, words_read, train_words)

			println("$(length(doc)) words read, $(position(file))/$end_pos")
			if length(doc) == 0
				break
			end

			inplace_train_vectors!(vm, doc, window_length,
				start_lr, train_words, words_read, total_ll;
				context_cut = context_cut, sense_treshold = sense_treshold)
		end

		close(file)
	end

	refs = Array{Future}(nworkers())

	for i in 1:nworkers()
		refs[i] = remotecall(do_work, i+1,  i)
	end

	for i in 1:nworkers()
		fetch(refs[i])
	end

	println("Learning complete $(words_read[1]) / $train_words")

	return words_read[1]
end


text:
""
function run_hopf_iv(file, model_type::Type, optimizer; kwargs...)
    return _PMs.run_model(file, model_type, optimizer, build_hopf_iv; ref_extensions=[ref_add_xfmr!],  solution_processors=[ _HPM.sol_data_model!], multinetwork=true, kwargs...)
end

""
function build_hopf_iv(pm::AbstractPowerModel)
    bounded= true

    for (n, network) in _PMs.nws(pm)
        _PMs.variable_bus_voltage(pm, nw=n, bounded=bounded)
        _PMs.variable_bus_voltage_magnitude_sqr(pm, nw=n, bounded=bounded)
        variable_transformer_voltage(pm, nw=n, bounded=false)

        _PMs.variable_branch_current(pm, nw=n, bounded=bounded)
        _PMs.variable_dcline_current(pm, nw=n)
        variable_transformer_current(pm, nw=n, bounded=false)

        variable_load_current(pm, nw=n, bounded=bounded)
        variable_gen_current(pm, nw=n, bounded=bounded)

    end

    for (n, network) in _PMs.nws(pm)
        for i in _PMs.ids(pm, :ref_buses, nw=n)
            constraint_ref_bus(pm, i, nw=n)
        end

        for i in _PMs.ids(pm, :bus, nw=n)
            constraint_current_balance(pm, i, nw=n)
            constraint_vm_auxiliary_variable(pm, i, nw=n)
            constraint_voltage_harmonics_relative_magnitude(pm, i, nw=n)
        end

        for g in _PMs.ids(pm, :gen, nw=n)
            _PMs.constraint_gen_active_bounds(pm, g, nw=n)
            _PMs.constraint_gen_reactive_bounds(pm, g, nw=n)
        end

        for i in _PMs.ids(pm, :load, nw=n)
            constraint_load_power(pm, i, nw=n)
        end

        for b in _PMs.ids(pm, :branch, nw=n)
            _PMs.constraint_current_from(pm, b, nw=n)
            _PMs.constraint_current_to(pm, b, nw=n)

            _PMs.constraint_voltage_drop(pm, b, nw=n)

            _PMs.constraint_current_limit(pm, b, nw=n)
        end

        for t in _PMs.ids(pm, :xfmr, nw=n)
            constraint_transformer_core_excitation(pm, t, nw=n)
            constraint_transformer_core_voltage_drop(pm, t, nw=n)
            constraint_transformer_core_voltage_balance(pm, t, nw=n)
            constraint_transformer_core_current_balance(pm, t, nw=n)

            constraint_transformer_winding_config(pm, t, nw=n)
            constraint_transformer_winding_current_balance(pm, t, nw=n)
        end

        # for d in _PMs.ids(pm, :dcline, nw=n)
        #     _PMs.constraint_dcline_power_losses(pm, d, nw=n)
        # end
    end

    #constraints across harmonics
    fundamental = 1
    for i in _PMs.ids(pm, :bus, nw=fundamental)
        constraint_voltage_magnitude_rms(pm, i)
        constraint_voltage_thd(pm, i, fundamental=fundamental)
    end

    for b in _PMs.ids(pm, :branch, nw=fundamental)
        constraint_current_limit_rms(pm, b)
    end

    for g in _PMs.ids(pm, :gen, nw=fundamental)
        constraint_active_filter(pm, g, fundamental=fundamental)
    end

    # _PMs.objective_min_fuel_and_flow_cost(pm)
    # objective_current_distortion_minimization(pm)
    objective_voltage_distortion_minimization(pm)
end

function ref_add_xfmr!(ref::Dict{Symbol,<:Any}, data::Dict{String,<:Any}) ## data not actually needed!
    _PMs.apply_pm!(_ref_add_xfmr!, ref, data)
end

function _ref_add_xfmr!(ref::Dict{Symbol,<:Any}, data::Dict{String,<:Any})
    if !haskey(ref, :xfmr)
        # error(_LOGGER, "required xfmr data not found")
        ref[:xfmr] = Dict()
        ref[:xfmr_arcs_from] = Dict()
        ref[:xfmr_arcs_to] = Dict()
        ref[:xfmr_arcs] = Dict()
        ref[:bus_arcs_xfmr] = Dict((i, []) for (i,bus) in ref[:bus])

    else
        ref[:xfmr] = Dict(x for x in ref[:xfmr] if  x.second["f_bus"] in keys(ref[:bus]) &&
                                                    x.second["t_bus"] in keys(ref[:bus])
            )

        ref[:xfmr_arcs_from] = [(t,xfmr["f_bus"],xfmr["t_bus"]) for (t,xfmr) in ref[:xfmr]]
        ref[:xfmr_arcs_to]   = [(t,xfmr["t_bus"],xfmr["f_bus"]) for (t,xfmr) in ref[:xfmr]]

        ref[:xfmr_arcs] = [ref[:xfmr_arcs_from]; ref[:xfmr_arcs_to]]

        bus_arcs_xfmr = Dict((i, []) for (i,bus) in ref[:bus])
        for (t,i,j) in ref[:xfmr_arcs]
            push!(bus_arcs_xfmr[i], (t,i,j))
        end
        ref[:bus_arcs_xfmr] = bus_arcs_xfmr
    end
end


text:
@testset "kernelpdmat" begin
    rng = MersenneTwister(123456)
    A = rand(rng, 10, 5)
    vecA = (RowVecs(A), ColVecs(A))
    a = rand(rng, 10)
    k = SqExponentialKernel()
    for obsdim in [1, 2]
        @test all(
            Matrix(kernelpdmat(k, A; obsdim=obsdim)) .≈
            Matrix(PDMat(kernelmatrix(k, A; obsdim=obsdim))),
        )
        @test kernelpdmat(k, vecA[obsdim]) == kernelpdmat(k, A; obsdim=obsdim)
        # @test_throws ErrorException kernelpdmat(k,ones(100,100),obsdim=obsdim)
    end
    @test @test_deprecated(kernelpdmat(k, A)) == kernelpdmat(k, ColVecs(A))
end


text:
# Toy linear ODE with 2 DOFs
# u_1_t - a * u_1 = 0
# u_2_t - b * u_1 - c * u_2 = 0

import GridapODEs.ODETools: ODEOperator
import GridapODEs.ODETools: AffineODEOperator
import GridapODEs.ODETools: ConstantODEOperator
import GridapODEs.ODETools: allocate_cache
import GridapODEs.ODETools: update_cache!
import GridapODEs.ODETools: allocate_residual
import GridapODEs.ODETools: jacobian!
import GridapODEs.ODETools: jacobian_t!
import GridapODEs.ODETools: jacobian_and_jacobian_t!
import GridapODEs.ODETools: allocate_jacobian
import GridapODEs.ODETools: residual!

struct ODEOperatorMock{T<:Real,C} <: ODEOperator{C}
  a::T
  b::T
  c::T
end

function residual!(r::AbstractVector,op::ODEOperatorMock,t::Real,u::AbstractVector,u_t::AbstractVector,ode_cache)
  r .= 0
  r[1] = u_t[1] - op.a * u[1]
  r[2] = u_t[2] - op.b * u[1] - op.c * u[2]
  r
end

function allocate_residual(op::ODEOperatorMock,u::AbstractVector,cache)
  zeros(2)
end

function jacobian!(J::AbstractMatrix,op::ODEOperatorMock,t::Real,u::AbstractVector,u_t::AbstractVector,ode_cache)
  J[1,1] += -op.a
  J[2,1] += -op.b
  J[2,2] += -op.c
  J
end

function jacobian_t!(J::AbstractMatrix,op::ODEOperatorMock,t::Real,u::AbstractVector,u_t::AbstractVector,du_t_u::Real,ode_cache)
  J[1,1] += 1.0*du_t_u
  J[2,2] += 1.0*du_t_u
  J
end

function jacobian_and_jacobian_t!(J::AbstractMatrix,op::ODEOperatorMock,t::Real,u::AbstractVector,u_t::AbstractVector,du_t_u::Real,ode_cache)
  jacobian!(J,op,t,u,u_t,ode_cache)
  jacobian_t!(J,op,t,u,u_t,du_t_u,ode_cache)
  J
end

function allocate_jacobian(op::ODEOperatorMock,u::AbstractVector,cache)
  zeros(2,2)
end

allocate_cache(op::ODEOperatorMock) = nothing
update_cache!(cache,op::ODEOperatorMock,t::Real) = nothing


text:
# separable.jl : separable Gaussian random fields

## SeparableCovarianceFunction ##
struct SeparableCovarianceFunction{d,V}
    cov::V
end

"""
	SeparableCovarianceFunction(cov...)

	Create a separable covariance function in `length(cov)` dimensions for the covariance structures `cov`. Usefull for defining anisotropic covariance functions, or if the usual KL expansion is too expensive.

# Examples
```
julia> e = Exponential(0.1)
exponential (λ=0.1, σ=1.0, p=2.0)

julia> m = Matern(0.01,1.0)
Matérn (λ=0.01, ν=1.0, σ=1.0, p=2.0)

julia> c = SeparableCovarianceFunction([e,m])
2d separable covariance function [ exponential (λ=0.1, σ=1.0, p=2.0), Matérn (λ=0.01, ν=1.0, σ=1.0, p=2.0) ]

```
See also: [`CovarianceFunction`](@ref), [`KarhunenLoeve`](@ref)
"""
SeparableCovarianceFunction(cov::Vector{T}) where {T<:CovarianceStructure} = SeparableCovarianceFunction{length(cov),Vector{T}}(cov)
SeparableCovarianceFunction(cov::T...) where {T<:CovarianceStructure}  = SeparableCovarianceFunction([cov...])

function GaussianRandomField(mean::Array{T} where {T<:Real},cov::SeparableCovarianceFunction{d,T} where {T},kl::KarhunenLoeve{n},pts::V...;kwargs...) where {d,n,V<:AbstractVector}
    all(size(mean).==length.(pts)) || throw(DimensionMismatch("size of the mean does not correspond to the dimension of the points"))
    length(pts) == d || throw(DimensionMismatch("number of point ranges must be equal to the dimension of the covariance function"))

    # generate the 1d grf's
    data = SpectralData[]
    for i in 1:length(cov.cov)
        cov_ = CovarianceFunction(1,cov.cov[i])
		if ( typeof(cov.cov[i]) <: Exponential && cov.cov[i].p == 1.)
			push!(data,compute_analytic(cov_,n,pts[i]))
		else
        	push!(data,_GaussianRandomField(mean,cov_,kl,pts[i];kwargs...).data)
		end
    end
    eigenval = [data[i].eigenval for i in 1:d]

    # determine n-d eigenvalues
    p = Base.product([data[i].eigenval for i in 1:d]...)
    m = map(prod,p)
    idx = sortperm(m[:],rev=true)
    pidx = Base.product(broadcast(:,1,length.(eigenval))...)

	GaussianRandomField{typeof(cov),typeof(kl),typeof(pts)}(mean,cov,pts,(collect(pidx)[idx],data))
end

# zero-mean GRF
GaussianRandomField(cov::SeparableCovarianceFunction{d,T} where {T},kl,pts::V...;kwargs...) where {d,V<:AbstractVector} = GaussianRandomField(zeros(eltype(pts[1]),length.(pts)...),cov,kl,pts...;kwargs...)

# constant mean GRF
GaussianRandomField(mean::R where {R<:Real},cov::SeparableCovarianceFunction{d,T} where {T},kl,pts::V...;kwargs...) where {d,V<:AbstractVector} = GaussianRandomField(mean*ones(eltype(pts[1]),length.(pts)...),cov,kl,pts...;kwargs...)

# return number of dimension
ndims(::SeparableCovarianceFunction{d}) where {d} = d

# sample function
function sample(grf::GaussianRandomField{C,KarhunenLoeve{n}} where {C<:SeparableCovarianceFunction}; xi::Vector{T} = randn(n)) where {n,T<:Real}
    length(xi) == n || throw(DimensionMismatch("length of random points vector must be equal to $(n)"))
    (order,data) = grf.data
	x = zeros(T,prod(size(grf.mean)))
	d = length(grf.cov.cov)
	for i in 1:n
        ev = prod([data[j].eigenval[order[i][j]] for j = 1:d])
        ef = kron([data[j].eigenfunc[:,order[i][j]] for j = 1:d]...)
		x += xi[i]*ev.*ef
	end
	grf.mean + prod([grf.cov.cov[i].σ for i = 1:d])*reshape(x, size(grf.mean))
end

function sample(grf::GaussianRandomField{C,KarhunenLoeve{n}} where {C<:SeparableCovarianceFunction{1}}; xi::Vector{T} = randn(n)) where {n,T<:Real}
    length(xi) == n || throw(DimensionMismatch("length of random points vector must be equal to $(n)"))
    (order,data) = grf.data
	x = data[1].eigenfunc*(data[1].eigenval.*xi)
	grf.mean + grf.cov.cov[1].σ*x
end

function show(io::IO, s::SeparableCovarianceFunction{d}) where {d}
	str = "$(d)d separable covariance function [ $(s.cov[1])"
	for i in 2:length(s.cov)
		str *= ", $(s.cov[i])"
    end
	str *= " ]"
	print(io, str)
end


text:
using PackageCompiler
using Setfield
using Pkg
using Random

const trace_dir = abspath(@__DIR__,"../traces/")
const trace_file = Vector{UInt8}()

function trace()
    global trace_dir,trace_file
    !isdir(trace_dir) && mkdir(trace_dir)

    empty!(trace_file)
    for c in joinpath(trace_dir,"trace_$(rand(UInt32)).jl")
        push!(trace_file,UInt8(c))
    end
    push!(trace_file,C_NULL)

    opts = Base.JLOptions()
    opts = @set opts.trace_compile = pointer(trace_file)
    opts = @set opts.compile_enabled = 2
    unsafe_store!(Base.cglobal(:jl_options,Base.JLOptions),opts)
    return nothing
end

function untrace()
    global trace_dir,cstr
    opts = Base.JLOptions()
    opts = @set opts.trace_compile = C_NULL
    opts = @set opts.compile_enabled = 1
    unsafe_store!(Base.cglobal(:jl_options,Base.JLOptions),opts)
    return nothing
end

function grep_possible_packages!(package_set::Set{Symbol},expr::Expr)
    (length(expr.args) == 0) && return
    s = expr.args[1]
    if expr.head == :.
        if typeof(s) === Symbol && s !== :.
            push!(package_set,s)
        end
    end
    foreach(expr.args) do a
        if typeof(a) == Expr
            grep_possible_packages!(package_set,a)
        end
    end
    nothing
end

semantic_fix(line) = begin
    line = replace(line,"(\"" => "(raw\"")
    line = replace(line,"DatePart{Char(" => "DatePart{(")
    line = replace(line,r"(?<!#)#s(?=\d)" => "T")
    line = replace(line,")()" => ")")
end

try_parse_line(line,linenum = 0,filename ="") = begin
    try
        expr = Meta.parse(line, raise = true)
        if expr.head != :incomplete
            return expr
        else
            @warn "skipping line with parsing error" linenum filename
            return nothing
        end
    catch e
        @warn "skipping line with parsing error" linenum filename
        return nothing
    end
end

brute_build_local() = brute_build_julia(replace = false)
export brute_build_local

function brute_build_julia(;clear_traces = true, replace = true)
    !isdir(trace_dir) && begin
        @info "no trace files found"
        return
    end
    if Sys.isunix()
        ENV["JULIA_CC"] = "/usr/bin/gcc"
    end
    blacklist = push!(Fezzik.blacklist(),"Main","##benchmark#","###compiledcall","VSCode")
    statements = Set{String}()
    packages = Set{Symbol}()
    for fname in readdir(trace_dir)
        counter = 0
        fname = abspath(trace_dir,fname)
        for line in eachline(fname)
            try
                counter += 1
                any(occursin.(blacklist,[line])) && continue
                # line = semantic_fix(line);
                expr = try_parse_line(line,counter,fname);
                isnothing(expr) && continue
                grep_possible_packages!(packages,expr)
                push!(statements,line);
            catch e
                statement = line
                @info "Failed to parse statement" statement
            end
        end
    end

    out_file = abspath(@__DIR__,"../","precomp.jl")
    @info "generating precompile"

    for pkg in packages
        @eval import $pkg
    end
    @info "dry run3"
    statements = filter(shuffle(collect(statements))) do st
        try
            res = eval(Meta.parse(st))
            if !res
                @warn st
            end
            return res
        catch err
            @show err
            return false
        end
    end

    open(out_file, "w") do io
        for pkg in packages
            println(io, "import $(pkg)")
        end
        for line in statements
            println(io, "$line")
        end
    end
    @info "used $(length(statements)) precompile statements"
    @show out_file
    compile_incremental(packages,out_file,replace)
    @info "DONE!!"

    clear_traces && isdir(trace_dir) && for f in readdir(trace_dir)
        file = joinpath(trace_dir,f)
        try
            rm(file;force = true)
        catch err
            @warn "failed to remove trace file" err file
        end
    end

    exit()
    nothing
end
export brute_build_julia

import Libdl
sysimage_size() = stat(unsafe_string(Base.JLOptions().image_file)).size/1024/1024
export sysimage_size

revert(;force = false) = begin
    PackageCompiler.restore_default_sysimage()
end

compile_incremental(packages,file,replace) = begin
    packages = filter(collect(packages)) do mod
        mod == :Base ? false : mod == :Core ? false : true
    end

    missing_packages = packages_not_in_project(project_ctx(),packages)
    try Pkg.add(missing_packages) catch err
        unresolved = map(x->String(x[1]),eachmatch(r"[*] (\w+)",err.msg))
        @show unresolved
        packages = setdiff(packages,Symbol.(unresolved))
        missing_packages = setdiff(missing_packages,unresolved)
        !isempty(missing_packages)  && Pkg.add(missing_packages)
    end

    @show Base.JLOptions().image_file |> unsafe_string |> deepcopy
    if replace
        PackageCompiler.create_sysimage(packages; precompile_statements_file = file , replace_default = true)
    else
        PackageCompiler.create_sysimage(packages; precompile_execution_file = file , sysimage_path = pwd() * "/JuliaSysimage." * PackageCompiler.Libdl.dlext,
        base_sysimage = Base.JLOptions().image_file |> unsafe_string |> deepcopy , include_transitive_dependencies = false, sysimage_build_args = `-O2 --check-bounds=no`)
    end
end

function packages_not_in_project(ctx, packages)
    packages_in_project = collect(keys(ctx.env.project.deps))
    if ctx.env.pkg !== nothing
        push!(packages_in_project, ctx.env.pkg.name)
    end
    setdiff(string.(packages), packages_in_project)
end

function project_ctx(project = dirname(Base.active_project()))
    @show project
    project_toml_path = Pkg.Types.projectfile_path(project; strict=true)
    if project_toml_path === nothing
        error("could not find project at $(repr(project))")
    end
    return Pkg.Types.Context(env=Pkg.Types.EnvCache(project_toml_path))
end


text:
# ---
# title: 200. Number of Islands
# id: problem200
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search, Breadth-first Search, Union Find
# link: <https://leetcode.com/problems/number-of-islands/description/>
# hidden: true
# ---
#
# Given an `m x n` 2d `grid` map of `'1'`s (land) and `'0'`s (water), return
# _the number of islands_.
#
# An **island** is surrounded by water and is formed by connecting adjacent
# lands horizontally or vertically. You may assume all four edges of the grid
# are all surrounded by water.
#
#
#
# **Example 1:**
#
#
#
#     Input: grid = [
#       ["1","1","1","1","0"],
#       ["1","1","0","1","0"],
#       ["1","1","0","0","0"],
#       ["0","0","0","0","0"]
#     ]
#     Output: 1
#
#
# **Example 2:**
#
#
#
#     Input: grid = [
#       ["1","1","0","0","0"],
#       ["1","1","0","0","0"],
#       ["0","0","1","0","0"],
#       ["0","0","0","1","1"]
#     ]
#     Output: 3
#
#
#
#
# **Constraints:**
#
#   * `m == grid.length`
#   * `n == grid[i].length`
#   * `1 <= m, n <= 300`
#   * `grid[i][j]` is `'0'` or `'1'`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
push!(LOAD_PATH,"../src/")
using FunctionOperators

## I "stole" that script from CMBLensing.jl (marius311)
# https://github.com/marius311/CMBLensing.jl/tree/gh-pages

using Documenter
# ensure in right directory
cd(dirname(@__FILE__))

# convert Jupyter notebooks to markdown and do some simple preprocessing before
# feeding into Documenter.jl

function convert_to_markdown(file)
    run(`jupyter nbconvert ../examples/$file --to markdown --template src/documenter.tpl --output-dir src-staging`)
    return "src-staging/$(replace(file, "ipynb"=>"md"))"
end

function convert_equations_and_search_lines(file)
    contents = read(file, String)
    contents = replace(contents, r"\$\$(.*?)\$\$"s => s"""```math
    \g<1>
    ```""")
    contents = replace(contents, r"\* \$(.*?)\$" => s"* ``\g<1>``") # starting a line with inline math screws up tex2jax for some reason
    contents = replace(contents, r"search.*$" => "")
    write(file, contents)
    return file
end

rm("src-staging", force=true, recursive=true)
mkdir("src-staging")
for file in readdir("../examples")
    if endswith(file, "ipynb")
        file |> convert_to_markdown |> convert_equations_and_search_lines
    elseif !startswith(file, ".")
        cp("src/$file", "src-staging/$file")
    end
end
cp("src/reference.md","src-staging/reference.md")
cp("../README.md","src-staging/index.md")
cp("src/assets/", "src-staging/assets")

# highlight output cells (i.e. anything withouout a language specified) white

@eval Documenter.Writers.HTMLWriter function mdconvert(c::Markdown.Code, parent::MDBlockContext; kwargs...)
    @tags pre code
    language = isempty(c.language) ? "none" : c.language
    pre[".language-$(language)"](code[".language-$(language)"](c.code))
end

# adds the MyBinder button on each page that is a notebook

#=@eval Documenter.Writers.HTMLWriter function render_article(ctx, navnode)
    @tags article header footer nav ul li hr span a img

    header_links = map(Documents.navpath(navnode)) do nn
        title = mdconvert(pagetitle(ctx, nn); droplinks=true)
        nn.page === nothing ? li(title) : li(a[:href => navhref(ctx, nn, navnode)](title))
    end

    topnav = nav(ul(header_links))

    pageurl = get(getpage(ctx, navnode).globals.meta, :EditURL, getpage(ctx, navnode).source)
    nbpath = foldl(replace,["src-staging"=>"../examples",".md"=>".ipynb"], init=pageurl)
    if isfile(nbpath)
        url = "https://mybinder.org/v2/gh/hakkelt/FunctionOperators.jl/master?urlpath=lab/tree/$(basename(nbpath))"
        push!(topnav.nodes, a[".edit-page", :href => url](img[:src => "https://mybinder.org/badge_logo.svg"]()))
    end

    art_header = header(topnav, hr(), render_topbar(ctx, navnode))

    # build the footer with nav links
    art_footer = footer(hr())
    if navnode.prev !== nothing
        direction = span[".direction"]("Previous")
        title = span[".title"](mdconvert(pagetitle(ctx, navnode.prev); droplinks=true))
        link = a[".previous", :href => navhref(ctx, navnode.prev, navnode)](direction, title)
        push!(art_footer.nodes, link)
    end

    if navnode.next !== nothing
        direction = span[".direction"]("Next")
        title = span[".title"](mdconvert(pagetitle(ctx, navnode.next); droplinks=true))
        link = a[".next", :href => navhref(ctx, navnode.next, navnode)](direction, title)
        push!(art_footer.nodes, link)
    end

    pagenodes = domify(ctx, navnode)
    article["#docs"](art_header, pagenodes, art_footer)
end=#

makedocs(
    sitename="FunctionOperators.jl",
    source = "src-staging",
    pages = [
        "index.md",
        "Tutorial.md",
        "reference.md"
    ],
    format = Documenter.HTML(assets = ["assets/style.css", "assets/favicon.ico"])
)

deploydocs(
    repo = "github.com/hakkelt/FunctionOperators.jl.git",
)


text:
"""
    RegressionKalmanState{Fl <: AbstractFloat}

Similar to the univariate Kalman filter, but exploits the structure of the regression system
matrices - in particular, the fact that the variance of the state (`P`) is always zero, as
well as `R` and `Q`.
"""
mutable struct RegressionKalmanState{Fl<:AbstractFloat}
    v::Fl
    F::Fl
    a::Vector{Fl}
    llk::Fl

    function RegressionKalmanState(a1::Vector{Fl}) where Fl
        return new{Fl}(zero(Fl), zero(Fl), a1, zero(Fl))
    end
end

function save_kalman_state_in_filter_output!(
    filter_output::FilterOutput{Fl}, kalman_state::RegressionKalmanState{Fl}
) where Fl
    filter_output.a[1] = copy(kalman_state.a)
    filter_output.P[1] = copy(fill(zero(Fl), 1, 1))
    for t in 1:length(filter_output.v)
        filter_output.v[t] = copy(fill(kalman_state.v, 1))
        filter_output.F[t] = copy(fill(kalman_state.F, 1, 1))
        filter_output.a[t + 1] = copy(kalman_state.a)
        filter_output.att[t] = copy(kalman_state.a)
        filter_output.P[t + 1] = copy(fill(zero(Fl), 1, 1))
        filter_output.Ptt[t] = copy(fill(zero(Fl), 1, 1))
        # There is no diffuse part in UnivariateKalmanState
        filter_output.Pinf[t] = copy(fill(zero(Fl), size(filter_output.Ptt[t])))
    end
    return filter_output
end

mutable struct RegressionKalmanFilter{Fl<:AbstractFloat} <: KalmanFilter
    a1::Vector{Fl}
    kalman_state::RegressionKalmanState

    function RegressionKalmanFilter(a1::Vector{Fl}) where Fl
        kalman_state = RegressionKalmanState(a1)
        return new{Fl}(a1, kalman_state)
    end
end

function reset_filter!(filter::RegressionKalmanFilter{Fl}) where Fl
    copyto!(filter.kalman_state.a, filter.a1)
    filter.kalman_state.llk = zero(Fl)
    return filter
end

function regression_update_v!(
    kalman_state::RegressionKalmanState{Fl}, y::Fl, Z::Vector{Fl}
) where Fl
    kalman_state.v = y - dot(Z, kalman_state.a)
    return kalman_state
end

function update_llk!(kalman_state::RegressionKalmanState{Fl}) where Fl
    kalman_state.llk -= (
        HALF_LOG_2_PI + 0.5 * (log(kalman_state.F) + kalman_state.v^2 / kalman_state.F)
    )
    return kalman_state
end

function regression_update_kalman_state!(kalman_state, y, Z)
    if isnan(y)
        kalman_state.v = NaN
    else
        regression_update_v!(kalman_state, y, Z)
        update_llk!(kalman_state)
    end
    return kalman_state
end

function optim_kalman_filter(
    sys::LinearUnivariateTimeVariant{Fl}, filter::RegressionKalmanFilter{Fl}
) where Fl
    return filter_recursions!(filter.kalman_state, sys)
end

function kalman_filter!(
    filter_output::FilterOutput, sys::StateSpaceSystem, filter::RegressionKalmanFilter{Fl}
) where Fl
    filter_recursions!(filter_output, filter.kalman_state, sys)
    return filter_output
end

function filter_recursions!(
    kalman_state::RegressionKalmanState{Fl}, sys::LinearUnivariateTimeVariant{Fl}
) where Fl
    kalman_state.F = sys.H[1]
    @inbounds for t in eachindex(sys.y)
        regression_update_kalman_state!(kalman_state, sys.y[t], sys.Z[t])
    end
    return kalman_state.llk
end

function filter_recursions!(
    filter_output::FilterOutput,
    kalman_state::RegressionKalmanState{Fl},
    sys::LinearUnivariateTimeVariant,
) where Fl
    kalman_state.F = sys.H[1]
    @inbounds for t in eachindex(sys.y)
        regression_update_kalman_state!(kalman_state, sys.y[t], sys.Z[t])
    end
    save_kalman_state_in_filter_output!(filter_output, kalman_state)
    return filter_output
end


text:
# function _precompile_()
#     ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
#     precompile(Tuple{typeof(SubMatrixSelectionSVD.α), Base.UnitRange{Int64}, Base.LinAlg.Symmetric{Float64, Array{Float64, 2}}})
#     precompile(Tuple{getfield(SubMatrixSelectionSVD, Symbol("##projectionscorefiltered#10")), Int64, Array{Float64, 1}, typeof(identity), Array{Float64, 2}, Base.UnitRange{Int64}, Array{Float64, 1}})
#     # precompile(Tuple{getfield(SubMatrixSelectionSVD, Symbol("##smssvd#19")), Int64, Int64, typeof(identity), Array{Float64, 2}, Int64, Array{Float64, 1}})
#     precompile(Tuple{getfield(SubMatrixSelectionSVD, Symbol("#kw##projectionscorefiltered")), Array{Any, 1}, typeof(SubMatrixSelectionSVD.projectionscorefiltered), Array{Float64, 2}, Base.UnitRange{Int64}, Array{Float64, 1}})
#     precompile(Tuple{typeof(SubMatrixSelectionSVD._αfiltered), Array{Float64, 2}, Base.UnitRange{Int64}, Array{Float64, 1}, Array{Float64, 1}})
#     precompile(Tuple{typeof(SubMatrixSelectionSVD.smssvd), Array{Float64, 2}, Int64, Array{Float64, 1}})
#     precompile(Tuple{typeof(SubMatrixSelectionSVD.smssvd), Array{Float64, 2}, Int64})
#     precompile(Tuple{typeof(SubMatrixSelectionSVD.smssvd), Array{Float64, 2}, Array{Int64, 1}})
# end


text:
# Infinite stream processing version 2, don't give <unk> as an input to model, which can handle infinite amount of tokenized-raw text data
const SOS = "<s>"
const EOS = "</s>"
const UNK = "<unk>"

# character constants
const PAD = Char(0x0011)
const SOW = Char(0x0012)
const EOW = Char(0x0013)

function create_vocab(vocabfile::AbstractString)
    ercount = 0
    result = Dict{AbstractString, Int}(SOS=>1, EOS=>2, UNK=>3)
    open(vocabfile) do f
        for line in eachline(f)
            words= split(line)
            try
                @assert(length(words) == 1, "The vocabulary file seems broken")
            catch e
                ercount += 1
                if (length(words) == 0)
                    if (ercount == 5)
                        info("There is empty line in vocabulary file")
                        ercount = 0
                    end
                    continue
                else
                    warn("Something unexpected happen in vocabfile")
                end
            end
            word = words[1]
            if length(word) > 65
                isinteractive() && warn("Too long words in file $vocabfile")
                word = word[1:65]
            end
            get!(result, word, 1+length(result))
        end
    end
    return result
end

"""
Reads from infinite stream modifies the s paramater in the light of sequence lengths,
it modifies the s argument and such that it contains what i called information tuple:
(outid, inid) -> outid is the the output vocabulary to use in softmax, inid is the input vocabulary to use in char level lstm
"""
function readstream!(f::IO,
                     s::Dict{Int64, Array{Any, 1}},
                     vocab::Dict{AbstractString, Int64},
                     realvocab::Dict{AbstractString, Int64};
                     maxlines=100, llimit=3, ulimit=300)

    k = 0
    while (k != maxlines)
        eof(f) && return false
        words = split(readline(f))
        wlen = length(words)
        (wlen < llimit) || (wlen > ulimit-2) && continue
        seq = Array(Tuple{Int32, Int32}, wlen+2)
        seq[1] = (vocab[SOS], vocab[SOS])
        for i=1:wlen
            word = words[i]
            index = get(vocab, word, vocab[UNK])
            realindex = get(realvocab, word, realvocab[UNK]) # no unking activity desired
            seq[i+1] = (index, realindex)
        end
        seq[end] = (vocab[EOS], vocab[EOS])
        skey = wlen+2
        (!haskey(s, skey)) && (s[skey] = Any[])
        push!(s[skey], seq)
        k += 1
    end
end

function mbatch(sequences::Array{Any, 1}, batchsize::Int)
    seqlen = length(sequences[1])
    data = Array(Any, seqlen)
    for cursor=1:seqlen
        d = Array(Tuple{Int32, Int32}, batchsize)
        for i=1:batchsize
            d[i] = sequences[i][cursor]
        end
        data[cursor] = d
    end
    return data
end

function nextbatch(f::IO,
                   sdict::Dict{Int64, Array{Any, 1}},
                   vocab::Dict{AbstractString, Int64},
                   realvocab::Dict{AbstractString, Int64},
                   batchsize;
                   o...)

    slens = collect(filter(x->length(sdict[x])>=batchsize, keys(sdict)))
    if length(slens) < 10
        readstream!(f, sdict, vocab, realvocab; o...)
        slens = collect(filter(x->length(sdict[x])>=batchsize, keys(sdict)))
        if length(slens) == 0
            return nothing
        end
    end
    slen = rand(slens)
    sequence = sdict[slen][1:batchsize]
    deleteat!(sdict[slen], 1:batchsize)
    return mbatch(sequence, batchsize)
end

"""
bytpe: out -> give "unk"ed version of the word, give normal version otherwise.
"""
function ibuild_sentence(i2w::Array{AbstractString, 1}, sequence::Array{Any, 1}, kth::Int; btype=:out, verbose=true)
    sentence = Any[]
    for i=1:length(sequence)
        z = sequence[i][kth]
        index = ((btype == :out) ? z[1] : z[2])
        push!(sentence, i2w[index])
    end
    if verbose
        for item in sentence; print("$item ");end;
    else
        return sentence
    end
end

function create_chvocab(word_vocab::Dict{AbstractString, Int64})
    res = Dict{Char, Int}(PAD=>1, SOW=>2, EOW=>3)
    for k in keys(word_vocab)
        for ch in k
            if ch == PAD || ch == SOW || ch == EOW
                warn("$ch is used in vocabulary")
            end
            (ch == ' ') && continue
            get!(res, ch, 1+length(res))
        end
    end
    return res
end

function create_chvocab(f::AbstractString, ulimit::Int=20000)
    res = Dict{Char, Int}(PAD=>1, SOW=>2, EOW=>3)
    counter = 0
    open(f) do file
        for line in eachline(file)
            counter += 1; (counter > ulimit) && break
            word = line
            # word = split(line)[2] no need because not to use count file
            for c in word
                if c == PAD || c == SOW || c == EOW
                    warn("$c is used in vocabulary")
                end
                (c == ' ') && continue
                get!(res, c, 1+length(res))
            end
        end
    end
    return res
end

longest_word{T}(word_vocab::Dict{T, Int}) = findmax(map(length, keys(word_vocab)))[1]

function ibuild_wordflup(i2c::Array{Char, 1}, data::Array{Any, 1}, kth::Int; verbose=false)
    word = Any[]
    for i=1:length(data)
        z = data[i][kth]
        push!(word, i2c[z])
    end
    if verbose
        for item in word; print("$item");end;println()
    else
        return word
    end
end

"""
Find the maximum length word in vocabulary, padds each word based on that length
"""
function create_conv_vocab(word_vocab_all::Dict{AbstractString, Int})
    maxlen = longest_word(word_vocab_all)
    result = Dict{AbstractString, Int}()
    for word in keys(word_vocab_all)
        word2 = string(SOW) * word * string(EOW)
        new_word = string(word2, string(PAD)^(maxlen-length(word)))
        result[new_word] = word_vocab_all[word]
    end
    return result
end

function cbatch4conv(wids::Array{Tuple{Int32, Int32}, 1}, i2w_all_conv::Array{AbstractString, 1}, ch::Dict{Char, Int}, atype::DataType)
    batchsize = length(wids)
    words = mapreduce(x->i2w_all_conv[x[2]], string, "", wids) # concatanete each word
    maxlen = length(words)

    data = zeros(Int32, maxlen)
    mask = zeros(Int32, maxlen, 1)
    for (i,c) in enumerate(words)
        if c == PAD
            data[i] = ch[PAD]
        else
            data[i] = get(ch, c, ch[PAD])
            mask[i, 1] = 1
        end
    end
    mcon = convert(atype, mask)
    return (data, mcon)
end

function calculate_coverage(vocabfile::AbstractString, vocabsize::Int)
    total = covered = 0.0
    open(vocabfile) do file
        counter = 0
        for line in eachline(file)
            counter += 1
            occurence = parse(Float64, split(line)[1])
            (counter <= vocabsize) && (covered += occurence)
            total += occurence
        end
    end
    rate = (covered / total) * 100
    println("Vocabsize $vocabsize | Coverage rate $rate")
    return rate
end

word_arrange(pad_critic::Int, x::AbstractString) = string(string(PAD)^(pad_critic - length(x) -2), string(SOW), x, string(EOW))

function charlup3(wids::Array{Tuple{Int32,Int32},1}, i2w_all::Array{AbstractString, 1}, ch::Dict{Char, Int})
    words = map(x->i2w_all[x[2]], wids)
    pad_critic = findmax(map(length, words))[1] + 2 # + 2 for start and and end characters
    batchsize = length(words)

    data = [ zeros(Int32, batchsize) for i=1:pad_critic ]
    masks = [ ones(Float32, batchsize, 1) for i=1:pad_critic ]
    words = map(x->word_arrange(pad_critic, x), words)

    @inbounds for i=1:batchsize
        word = words[i]
        for (cursor, char) in enumerate(word)
            (char == PAD) && (masks[cursor][i] = 0.0)
            data[cursor][i] = get(ch, char, ch[PAD])
        end
    end
    return data, masks
end

## deprecated code
# """
# wids : batch of words in (outid, inid) -> outid is the softmax layer id, inid is the char level id,
# charlup only returns the input to char level lstms
# """
# function charlup(wids::Array{Tuple{Int32,Int32},1}, i2w_all::Array{AbstractString, 1}, ch::Dict{Char, Int})
#     words = map(x->i2w_all[x[2]], wids)
#     critic = findmax(map(length, words))[1]

#     batchsize = length(words)
#     data = Array(Any, critic+2)
#     data[1] = fill!(zeros(Int32, batchsize), ch[SOW])
#     masks = Array(Any, critic+2)
#     masks[1] = ones(Int32, batchsize, 1)
#     for cursor=1:critic+1 # to pad EOW to the end of the word
#         d = Array(Int32, batchsize)
#         mask = ones(Float32, batchsize, 1)
#         @inbounds for i=1:batchsize
#             word = words[i]
#             if length(word) < critic
#                 if length(word) >= cursor
#                     try
#                         d[i] = get(ch, word[cursor], ch[PAD]) #ch[word[cursor]], there may be unk characters
#                     catch
#                         d[i] = ch[PAD]
#                     end
#                 elseif length(word)+1 == cursor
#                     d[i] = ch[EOW]
#                 else
#                     d[i] = ch[PAD]
#                     mask[i] = 0
#                 end
#             else
#                 if cursor>critic
#                     d[i] = ch[EOW]
#                 else
#                     try
#                         d[i] = get(ch, word[cursor], ch[PAD]) # unking operation
#                     catch
#                         d[i] = ch[PAD]
#                     end
#                 end
#             end
#         end
#         data[cursor + 1] = d
#         masks[cursor + 1] = mask
#     end
#     return data, masks
# end

# function charlup2(wids::Array{Tuple{Int32,Int32},1}, i2w_all::Array{AbstractString, 1}, ch::Dict{Char, Int})
#     words = map(x->i2w_all[x[2]], wids)
#     pad_critic = findmax(map(length, words))[1] + 2 # + 2 for start and and end characters

#     batchsize = length(words)
#     data = Array(Any, pad_critic)
#     masks = Array(Any, pad_critic)

#     for cursor=1:pad_critic-1
#         d = Array(Int32, batchsize)
#         mask = ones(Float32, batchsize, 1)
#         @inbounds for i=1:batchsize
#             word = words[i]
#             tostart = pad_critic - (length(word)+2) + 1
#             if cursor == tostart
#                 d[i] = ch[SOW]
#             elseif cursor > tostart
#                 try
#                    ch_place = cursor - tostart # it shows which chracter of word in
#                    d[i] = ch[word[ch_place]]
#                 catch
#                     d[i] = ch[PAD]
#                     mask[i] = 0.0
#                 end
#             else
#                 d[i] = ch[PAD]
#                 mask[i] = 0.0
#             end
#         end
#         data[cursor] = d
#         masks[cursor] = mask
#     end
#     data[end] = fill!(zeros(Int32, batchsize), ch[EOW])
#     masks[end] = ones(Float32, batchsize, 1)
#     return data, masks
# end

max_stars_repo_name: kirnap/KU-dependencyParser

max_issues_repo_name: kirnap/KU-dependencyParser

max_forks_repo_name: kirnap/KU-dependencyParser


text:
using TimeZones.TZData: TimeOffset, value, hour, minute, second, hourminutesecond

# Time seconds constructor
t = TimeOffset(5025)
@test hour(t) == 1
@test minute(t) == 23
@test second(t) == 45
@test value(t) == 5025
@test hourminutesecond(t) == (1, 23, 45)

t = TimeOffset(-5025)
@test hour(t) == -1
@test minute(t) == -23
@test second(t) == -45
@test value(t) == -5025
@test hourminutesecond(t) == (-1, -23, -45)

t = TimeOffset(0,61,61)
@test t == TimeOffset(1,2,1)
@test value(t) == 3721

t = TimeOffset(1,-23,-45)
@test t == TimeOffset(0,36,15)
@test value(t) == 2175

# TimeOffset String constructor
@test TimeOffset("1") == TimeOffset(1,0,0)  # See Pacific/Apia rules for an example.
@test TimeOffset("1:23") == TimeOffset(1,23,0)
@test TimeOffset("1:23:45") == TimeOffset(1,23,45)
@test TimeOffset("-1") == TimeOffset(-1,0,0)
@test TimeOffset("-1:23") == TimeOffset(-1,-23,0)
@test TimeOffset("-1:23:45") == TimeOffset(-1,-23,-45)
@test_throws ArgumentError TimeOffset("1:-23:45")
@test_throws ArgumentError TimeOffset("1:23:-45")
@test_throws ArgumentError TimeOffset("1:23:45:67")

@test string(TimeOffset(1,23,45)) == "01:23:45"
@test string(TimeOffset(-1,-23,-45)) == "-01:23:45"
@test string(TimeOffset(0,0,0)) == "00:00:00"
@test string(TimeOffset(0,-1,0)) == "-00:01:00"
@test string(TimeOffset(0,0,-1)) == "-00:00:01"
@test string(TimeOffset(24,0,0)) == "24:00:00"

# Math
@test TimeOffset(1,23,45) + TimeOffset(-1,-23,-45) == TimeOffset(0)
@test TimeOffset(1,23,45) - TimeOffset(1,23,45) == TimeOffset(0)

# TimeOffset show function
t = TimeOffset(1,23,45)
@test sprint(print, t) == "01:23:45"
@test sprint(show, t) == "01:23:45"


text:
## https://github.com/SciML/DifferentialEquations.jl/issues/849
function ArrayInterface.lu_instance(x::ComponentMatrix)
    T = eltype(x)
    noUnitT = typeof(zero(T))
    luT = LinearAlgebra.lutype(noUnitT)
    ipiv = Vector{LinearAlgebra.BlasInt}(undef, 0)
    info = zero(LinearAlgebra.BlasInt)
    return LU{luT}(similar(x), ipiv, info)
end

# Helpers for dealing with adjoints and such
_first_axis(x::AbstractComponentVecOrMat) = getaxes(x)[1]

_second_axis(x::AbstractMatrix) = FlatAxis()
_second_axis(x::ComponentMatrix) = getaxes(x)[2]
_second_axis(x::AdjOrTransComponentVecOrMat) = getaxes(x)[2]

_out_axes(::typeof(*), a, b::AbstractVector) = (_first_axis(a), )
_out_axes(::typeof(*), a, b::AbstractMatrix) = (_first_axis(a), _second_axis(b))
_out_axes(::typeof(\), a, b::AbstractVector) = (_second_axis(a), )
_out_axes(::typeof(\), a, b::AbstractMatrix) = (_second_axis(a), _second_axis(b))
_out_axes(::typeof(/), a::AbstractMatrix, b) = (_first_axis(a), _first_axis(b))

# Arithmetic
for op in [:*, :\, :/]
    @eval begin
        function Base.$op(A::AbstractComponentVecOrMat, B::AbstractComponentVecOrMat)
            C = $op(getdata(A), getdata(B))
            ax = _out_axes($op, A, B)
            return ComponentArray(C, ax)
        end
    end
    for (adj, Adj) in zip([:adjoint, :transpose], [:Adjoint, :Transpose])
        @eval begin
            function Base.$op(aᵀ::$Adj{T,<:ComponentVector}, B::AbstractComponentMatrix) where {T}
                cᵀ = $op(getdata(aᵀ), getdata(B))
                ax2 = _out_axes($op, aᵀ, B)[2]
                return $adj(ComponentArray(cᵀ', ax2))
            end
            function Base.$op(A::$Adj{T,<:CV}, B::CV) where {T<:Real, CV<:ComponentVector{T}}
                return $op(getdata(A), getdata(B))
            end
        end
    end
end

max_stars_repo_name: jonniedie/ComponentArrays

max_issues_repo_name: jonniedie/ComponentArrays

max_forks_repo_name: jonniedie/ComponentArrays


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/UMAT.jl/blob/master/LICENSE

"""
!cmname    - DPC_LH1  = Drucker-Prager Linear hardening
!cmname    - DPC_EXP1 = Drucker-Prager Exponential hardening

! ------------------- ABAQUS:
!
  PROPS(1)  = 200.0d0   ! Sigma_t
  PROPS(2)  = 200.0d0   ! Sigma_c
  PROPS(3)  =  100.00d2   ! K_Zero
  PROPS(4)  =   100.0d0   ! C_const
  PROPS(5)  =   1.0     ! N_Power
  PROPS(6)  =   5.0d0    ! Eta
  PROPS(7)  =   200.0d0   ! A_Const
  PROPS(8)  =   0.0d0   ! B_Const
  PROPS(9)  =   2.0D5   ! E
  PROPS(10) =   0.3d0   ! NU
  PROPS(11) =   0.0d0   ! DT (jos todellinen DT=0)
  PROPS(12) = 1.0       ! Print All to DAT
  PROPS(13) = 1.e-13      ! FTOLER
  PROPS(14) = 20.0      ! MAXITER
  PROPS(15) = 1.0       ! ELNUM
  PROPS(16) = 1.0       ! IPNUM
  PROPS(17) = 2.0       ! MTYP
  PROPS(18) = 1.0       ! HARDTYP
  DTIME = .00011      ! Aika-askel DT
  !DTIME = 1.0        ! Aika-askel DT

  cmname  = 'DPC_LH1'
  !cmname  = 'DPO_EXP1'
  ndi = 3
  nshr = 3
  Stress = 0.0d0
  ddsdde = 0.0d0
  dstran = 0.0d0
  stran  = 0.0d0
  Stran_Old = 0.0D0
  dstran(3) =  5.0d-5                 ! Venymainkrementti (vakio)
  dstran(2) = -PROPS(10) * dstran(3)  ! Venymainkrementti (vakio)
  dstran(1) = -PROPS(10) * dstran(3)  ! Venymainkrementti (vakio)
"""
function DruckerPragerMaterial(;Sigma_t = 200.,
                               Sigma_c = 200.,
                               K_Zero = 100.,
                               C_const = 100.,
                               N_Power = 1.,
                               Eta = 5.,
                               A_Const = 200.,
                               B_Const = 0.,
                               E = 2.e5,
                               NU = 0.3,
                               DT = 0.,
                               PRNTYP = false, # remember to convert float
                               FTOLER = 1.e-13,
                               MAXITER = 20,   # remember to convert float
                               ELNUM = 1,      # remember to convert float
                               IPNUM = 1,      # remember to convert float
                               MTYP = 2,       # remember to convert float
                               HARDTYP = 1,    # remember to convert float
                               DPC="DPC_LH1")

    umat_other = UmatOtherState(NTENS=6, NSTATV=1, NPROPS=18, CMNAME=DPC)
    umat_params = [Sigma_t, Sigma_c, K_Zero, C_const, N_Power, Eta, A_Const,
                   B_Const, E, NU, DT, float(PRNTYP), FTOLER, float(MAXITER),
                   float(ELNUM), float(IPNUM), float(MTYP), float(HARDTYP)]

    NTENS=6
    NSTATV=1
    NPROPS=18
    return UmatMaterial(NTENS=NTENS,
                        NSTATV=NSTATV,
                        NPROPS=NPROPS,
                        #drivers=UmatDriverState(NTENS=NTENS),
                        #ddrivers=UmatDriverState(NTENS=NTENS),
                        parameters=UmatParameterState(NPROPS=NPROPS, PROPS=umat_params),
                        #variables=UmatVariableState(NTENS=NTENS, NSTATV=NSTATV),
                        #variables_new=UmatVariableState(NTENS=NTENS, NSTATV=NSTATV),
                        #dparameters=UmatParameterState(NPROPS=NPROPS),
                        umat_other=umat_other,
                        lib_path=joinpath(lib_dir,"libdrucker_prager_plasticity." * dlext))
end


text:
using GLFW

struct KeyEvent
    key::GLFW.Key
    scancode::Int32
    action::GLFW.Action
    mods::Int32
end

abstract type AbstractContext end
struct VulkanContext <: AbstractContext end
struct OpenGLContext <: AbstractContext end

shift = 1
ctrl = 2
alt = 4

function ismod(name::AbstractString, event::KeyEvent)
    if name == "alt"
        return div(event.mods, 4) == 1
    elseif name == "ctrl"
        return div(mod(event.mods, 4), 2) == 1
    elseif name == "shift"
        return mod(mod(event.mods, 4), 2) == 1
    end
end

function iskey(name::AbstractString, event::KeyEvent)
    if occursin("+", name)
        if occursin("ctrl", name)
            consumed_name = replace(name, "ctrl+" => "")
            return ismod("ctrl", event)  && iskey(consumed_name, event)
        elseif occursin("shift", name)
            consumed_name = replace(name, "shift+" => "")
            return ismod("shift", event)  && iskey(consumed_name, event)
        elseif occursin("alt", name)
            consumed_name = replace(name, "alt+" => "")
            return ismod("alt", event)  && iskey(consumed_name, event)
        end
    else
        return getfield(GLFW, Symbol(uppercase("key_$name"))) == event.key
    end
end

function iskey(name::AbstractString, action::AbstractString, event::KeyEvent)
    if action == "released"
        action = "release"
    elseif action == "pressed"
        action = "press"
    end
    return iskey(name, event) && isaction(action, event)
end

function isaction(name::AbstractString, event::KeyEvent)
    getfield(GLFW, Symbol(uppercase(name))) == event.action
end

function wrap_key_callback(f::Function)
    function callback(window::GLFW.Window, key::GLFW.Key, scancode::Int32, action::GLFW.Action, mods::Int32)
        event = KeyEvent(key, scancode, action, mods)
        f(window, event)
    end
    return callback
end

function window_loop(context::OpenGLContext, renderer::Function)
    function loop(window::GLFW.Window)
        renderer()
        GLFW.SwapBuffers(window)
        GLFW.PollEvents()
    end
end

function window_loop(context::VulkanContext, renderer::Function)
    function loop(window::GLFW.Window)
        renderer()
        GLFW.PollEvents()
    end
end

function provide_window_hints(context::VulkanContext)
    GLFW.WindowHint(GLFW.CLIENT_API, GLFW.NO_API)
end

function provide_window_hints(context::OpenGLContext)
end

function create_window(window_size::Tuple{Integer,Integer}, name, context::OpenGLContext)
    provide_window_hints(context)
    window = GLFW.CreateWindow(window_size..., name)
    GLFW.MakeContextCurrent(window)
    return window
end

function create_window(window_size::Tuple{Integer,Integer}, name, context::VulkanContext)
    provide_window_hints(context)
    window = GLFW.CreateWindow(window_size..., name)
    return window
end

function run_window(window, loop; key_callback = nothing)
    if !isnothing(key_callback)
        callback = wrap_key_callback(key_callback)
        GLFW.SetKeyCallback(window, callback)
    end
    while !GLFW.WindowShouldClose(window)
        loop(window)
    end
    println("Good bye!")
    GLFW.DestroyWindow(window)
end

function run_window(window_size::Tuple{Integer,Integer}; key_callback = nothing, renderer = ()->nothing, name::AbstractString = "Julia Game Engine", context::AbstractContext = OpenGLContext())
    window = create_window(window_size, name, context)
    loop = window_loop(context, renderer)
    run_window(window, loop, key_callback = key_callback)
end


text:
#transfer
transfer_left(v::MPSBondTensor, A::GenericMPSTensor, Ab::GenericMPSTensor=A) =
    _transpose_as(_transpose_front(Ab)' * _transpose_front(_transpose_front(v)*_transpose_tail(A)), v)
transfer_right(v::MPSBondTensor, A::GenericMPSTensor, Ab::GenericMPSTensor=A) =
    _transpose_as(_transpose_tail(_transpose_front(A)*_transpose_front(v)) * _transpose_tail(Ab)', v)

#transfer for 2 mpo tensors
transfer_left(v::MPSBondTensor,A::MPOTensor,B::MPOTensor) = @tensor t[-1;-2] := v[1,2]*A[2,3,-2,4]*conj(B[1,3,-1,4])
transfer_right(v::MPSBondTensor,A::MPOTensor,B::MPOTensor) = @tensor t[-1;-2] := A[-1,3,1,4]*conj(B[-2,3,2,4])*v[1,2]

#transfer, but there are utility legs in the middle that are passed through
transfer_left(v::AbstractTensorMap{S,N1,N2},A::GenericMPSTensor{S,N},Ab::GenericMPSTensor{S,N}=A) where {S,N1,N2,N} =
    _permute_as(Ab'*permute(_permute_front(v)*_permute_tail(A),tuple(1,ntuple(x->N1+N2-1+x,N-1)...),tuple(ntuple(x->x+1,N1+N2-2)...,N1+N2+N-1)),v)

transfer_right(v::AbstractTensorMap{S,N1,N2},A::GenericMPSTensor{S,N},Ab::GenericMPSTensor{S,N}=A) where {S,N1,N2,N} =
    _permute_as(permute(A*_permute_tail(v),tuple(1,ntuple(x->N+x,N1+N2-2)...),tuple(ntuple(x->x+1,N-1)...,N1+N2+N-1))*_permute_tail(Ab)',v)

#mpo transfer
transfer_left(v::MPSTensor,O::MPOTensor,A::MPSTensor,Ab::MPSTensor) = @tensor v[-1 -2;-3] := v[4,2,1]*A[1,3,-3]*O[2,5,-2,3]*conj(Ab[4,5,-1])
transfer_right(v::MPSTensor,O::MPOTensor,A::MPSTensor,Ab::MPSTensor) = @tensor v[-1 -2;-3] := A[-1,4,5]*O[-2,2,3,4]*conj(Ab[-3,2,1])*v[5,3,1]

#utility, allowing transfering with arrays
function transfer_left(v,A::AbstractArray,Ab::AbstractArray=A;rvec=nothing,lvec=nothing)
    for (a,b) in zip(A,Ab)
        v = transfer_left(v,a,b)
    end

    if rvec != nothing && lvec != nothing
        if v isa MPSBondTensor #normal transfer
            @tensor v[-1;-2]-=rvec[1,2]*v[2,1]*lvec[-1,-2]
        elseif v isa MPSTensor #utiity leg in the middle
            @tensor v[-1 -2;-3]-=rvec[1,2]*v[2,-2,1]*lvec[-1,-3]
        else #what have you just given me?
            @assert false
        end
    end

    return v
end
function transfer_right(v,A::AbstractArray,Ab::AbstractArray=A;rvec=nothing,lvec=nothing)
    for (a,b) in Iterators.reverse(zip(A,Ab))
        v = transfer_right(v,a,b)
    end

    if rvec != nothing && lvec != nothing
        if v isa MPSBondTensor #normal transfer
            @tensor v[-1;-2]-=lvec[1,2]*v[2,1]*rvec[-1,-2]
        elseif v isa MPSTensor #utiity leg in the middle
            @tensor v[-1 -2;-3]-=lvec[1,2]*v[2,-2,1]*rvec[-1,-3]
        else
            @assert false
        end
    end

    return v
end
transfer_left(v,O::AbstractArray,A::AbstractArray,Ab::AbstractArray) = reduce((v,x)->transfer_left(v,x[1],x[2],x[3]),zip(O,A,Ab),init=v)
transfer_right(v,O::AbstractArray,A::AbstractArray,Ab::AbstractArray) = reduce((v,x)->transfer_right(v,x[1],x[2],x[3]),Iterators.reverse(zip(O,A,Ab)),init=v)


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    Map <: AbstractDict{MOI.VariableIndex,AbstractBridge}

Mapping between bridged constraints and the bridge that bridged the constraint.
"""
struct Map <: AbstractDict{MOI.ConstraintIndex,AbstractBridge}
    # Constraint Index of bridged constraint -> Bridge.
    # It is set to `nothing` when the constraint is deleted.
    bridges::Vector{Union{Nothing,AbstractBridge}}
    # Constraint Index of bridged constraint -> Constraint type.
    constraint_types::Vector{Tuple{Type,Type}}
    # The order of the keys is used in `keys_of_type` which is used by
    # `ListOfConstraintIndices`. Therefore they need to be in the order
    # of creation so we need `OrderedDict` and not `Dict`.
    # For `VariableIndex` constraints: (variable, set type) -> bridge
    single_variable_constraints::OrderedDict{Tuple{Int64,Type},AbstractBridge}

    function Map()
        return new(
            Union{Nothing,AbstractBridge}[],
            Tuple{Type,Type}[],
            OrderedDict{Tuple{Int64,Type},AbstractBridge}(),
        )
    end
end

# Implementation of `AbstractDict` interface.

function Base.isempty(map::Map)
    return all(bridge -> bridge === nothing, map.bridges) &&
           isempty(map.single_variable_constraints)
end

function Base.empty!(map::Map)
    empty!(map.bridges)
    empty!(map.constraint_types)
    empty!(map.single_variable_constraints)
    return map
end

function Base.length(map::Map)
    return count(bridge -> bridge !== nothing, map.bridges) +
           length(map.single_variable_constraints)
end

_index(ci::MOI.ConstraintIndex) = ci.value

_index(ci::MOI.ConstraintIndex{MOI.VectorOfVariables}) = -ci.value

function Base.haskey(map::Map, ci::MOI.ConstraintIndex{F,S}) where {F,S}
    return 1 <= _index(ci) <= length(map.bridges) &&
           map.bridges[_index(ci)] !== nothing &&
           (F, S) == map.constraint_types[_index(ci)]
end

Base.getindex(map::Map, ci::MOI.ConstraintIndex) = map.bridges[_index(ci)]

function Base.haskey(
    map::Map,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,S},
) where {S}
    return haskey(map.single_variable_constraints, (ci.value, S))
end

function Base.getindex(
    map::Map,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,S},
) where {S}
    return map.single_variable_constraints[(ci.value, S)]
end

function Base.delete!(map::Map, ci::MOI.ConstraintIndex)
    map.bridges[_index(ci)] = nothing
    return map
end

function Base.delete!(
    map::Map,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,S},
) where {S}
    delete!(map.single_variable_constraints, (ci.value, S))
    return map
end

function Base.values(map::Map)
    return Base.Iterators.flatten((
        # See comment in `values(::Variable.Map)`.
        Base.Iterators.Filter(bridge -> bridge !== nothing, map.bridges),
        values(map.single_variable_constraints),
    ))
end

# Implementation of iterate: it should combine non-`VariableIndex` constraints
# and `VariableIndex` constraints.
function _iterate_sv(
    map::Map,
    elem_state = iterate(map.single_variable_constraints),
)
    if elem_state === nothing
        return nothing
    end
    (i, S), bridge = elem_state[1]
    ci = MOI.ConstraintIndex{MOI.VariableIndex,S}(i)
    return ci => bridge, (2, elem_state[2])
end

_index(index, F, S) = MOI.ConstraintIndex{F,S}(index)

function _index(index, F::Type{MOI.VectorOfVariables}, S)
    return MOI.ConstraintIndex{F,S}(-index)
end

function _iterate(map::Map, state = 1)
    while state ≤ length(map.bridges) && map.bridges[state] === nothing
        state += 1
    end
    if state > length(map.bridges)
        return _iterate_sv(map)
    end
    F, S = map.constraint_types[state]
    return _index(state, F, S) => map.bridges[state], (1, state + 1)
end

Base.iterate(map::Map) = _iterate(map)

function Base.iterate(map::Map, state)
    if state[1] == 1
        return _iterate(map, state[2])
    end
    return _iterate_sv(map, iterate(map.single_variable_constraints, state[2]))
end

# Custom interface for information needed by `AbstractBridgeOptimizer`s that is
# not part of the `AbstractDict` interface.

"""
    number_of_type(map::Map, C::Type{<:MOI.ConstraintIndex})

Return the number of keys of type `C` in `map`.
"""
function number_of_type(map::Map, ::Type{MOI.ConstraintIndex{F,S}}) where {F,S}
    return count(i -> haskey(map, _index(i, F, S)), eachindex(map.bridges))
end

function number_of_type(
    map::Map,
    ::Type{MOI.ConstraintIndex{MOI.VariableIndex,S}},
) where {S}
    ret = 0
    for (_, key) in keys(map.single_variable_constraints)
        ret += Int(key == S)
    end
    return ret
end

"""
    keys_of_type(map::Map, C::Type{<:MOI.ConstraintIndex})

Return a list of all the keys of type `C` in `map` in order order in which they
were created with `add_key_for_bridge`.
"""
function keys_of_type(map::Map, C::Type{MOI.ConstraintIndex{F,S}}) where {F,S}
    return Base.Iterators.Filter(
        ci -> haskey(map, ci),
        MOI.Utilities.LazyMap{C}(i -> _index(i, F, S), eachindex(map.bridges)),
    )
end

function keys_of_type(
    map::Map,
    C::Type{MOI.ConstraintIndex{MOI.VariableIndex,S}},
) where {S}
    return MOI.Utilities.LazyMap{C}(
        key -> C(key[1]),
        Base.Iterators.Filter(
            key -> key[2] == S,
            keys(map.single_variable_constraints),
        ),
    )
end

"""
    list_of_key_types(map::Map)

Return a list of all the different concrete type of keys in `map`.
"""
function list_of_key_types(map::Map)
    list = Set{Tuple{Type,Type}}()
    for i in eachindex(map.bridges)
        if map.bridges[i] !== nothing
            push!(list, map.constraint_types[i])
        end
    end
    for (_, key) in keys(map.single_variable_constraints)
        push!(list, (MOI.VariableIndex, key))
    end
    return list
end

"""
    variable_constraints(map::Map, vi::MOI.VariableIndex)

Return the list of all keys corresponding to [`MathOptInterface.VariableIndex`](@ref)
constraints on the variable `vi`.
"""
function variable_constraints(map::Map, vi::MOI.VariableIndex)
    cis = MOI.ConstraintIndex{MOI.VariableIndex}[]
    for key in keys(map.single_variable_constraints)
        if key[1] == vi.value
            push!(cis, MOI.ConstraintIndex{MOI.VariableIndex,key[2]}(vi.value))
        end
    end
    return cis
end

"""
    vector_of_variables_constraints(map::Map)

Return the list of all keys that correspond to
[`MathOptInterface.VectorOfVariables`](@ref) constraints.
"""
function vector_of_variables_constraints(map::Map)
    return MOI.Utilities.LazyMap{MOI.ConstraintIndex{MOI.VectorOfVariables}}(
        i -> MOI.ConstraintIndex{map.constraint_types[i]...}(-i),
        Base.Iterators.Filter(
            i ->
                map.bridges[i] !== nothing &&
                    map.constraint_types[i][1] == MOI.VectorOfVariables,
            eachindex(map.bridges),
        ),
    )
end

"""
    has_bridges(map::Map)::Bool

Return a `Bool` indicating whether any bridge was added yet.

Note that it returns `false` even if all bridges were deleted while `isempty`
would return `true`.

It is computed in `O(1)` while `isempty` needs `O(n)` hence it is used by
[`MOI.Bridges.AbstractBridgeOptimizer`](@ref) to shortcut operations in case
variable bridges are not used.
"""
function has_bridges(map::Map)
    return !isempty(map.bridges) || !isempty(map.single_variable_constraints)
end

"""
    add_key_for_bridge(
        map::Map,
        bridge::AbstractBridge,
        func::MOI.AbstractFunction,
        set::MOI.AbstractSet,
    )

Return a new constraint index `ci` and store the mapping `ci => bridge`.
"""
function add_key_for_bridge end

function add_key_for_bridge(
    map::Map,
    bridge::AbstractBridge,
    ::F,
    ::S,
) where {F<:MOI.AbstractFunction,S<:MOI.AbstractSet}
    push!(map.bridges, bridge)
    push!(map.constraint_types, (F, S))
    return _index(length(map.bridges), F, S)
end

function add_key_for_bridge(
    map::Map,
    bridge::AbstractBridge,
    func::MOI.VariableIndex,
    ::S,
) where {S<:MOI.AbstractScalarSet}
    map.single_variable_constraints[(func.value, S)] = bridge
    return MOI.ConstraintIndex{MOI.VariableIndex,S}(func.value)
end

"""
    EmptyMap <: AbstractDict{MOI.ConstraintIndex, AbstractBridge}

Empty version of [`Map`](@ref). It is used by
[`MOI.Bridges.Variable.SingleBridgeOptimizer`](@ref) as it does
not bridge any constraint.
"""
struct EmptyMap <: AbstractDict{MOI.VariableIndex,AbstractBridge} end

Base.isempty(::EmptyMap) = true

Base.empty!(::EmptyMap) = nothing

Base.keys(::EmptyMap) = MOI.Utilities.EmptyVector{MOI.VariableIndex}()

Base.values(::EmptyMap) = MOI.Utilities.EmptyVector{AbstractBridge}()

has_bridges(::EmptyMap) = false

number_of_type(::EmptyMap, ::Type{<:MOI.ConstraintIndex}) = 0


text:
"""
    UnitSimplexOracle(right_side)

Represents the scaled unit simplex:
```
C = {x ∈ R^n_+, ∑x ≤ right_side}
```
"""
struct UnitSimplexOracle{T} <: LinearMinimizationOracle
    right_side::T
end

UnitSimplexOracle{T}() where {T} = UnitSimplexOracle{T}(one(T))

UnitSimplexOracle(rhs::Integer) = UnitSimplexOracle{Rational{BigInt}}(rhs)

"""
LMO for scaled unit simplex:
`∑ x_i = τ`
Returns either vector of zeros or vector with one active value equal to RHS if
there exists an improving direction.
"""
function compute_extreme_point(lmo::UnitSimplexOracle{T}, direction; v = nothing, kwargs...) where {T}
    idx = argmin(direction)
    if direction[idx] < 0
        return ScaledHotVector(lmo.right_side, idx, length(direction))
    end
    return ScaledHotVector(zero(T), idx, length(direction))
end

function convert_mathopt(
    lmo::UnitSimplexOracle{T},
    optimizer::OT;
    dimension::Integer,
    kwargs...,
) where {T,OT}
    MOI.empty!(optimizer)
    τ = lmo.right_side
    n = dimension
    (x, _) = MOI.add_constrained_variables(optimizer, [MOI.Interval(0.0, 1.0) for _ in 1:n])
    MOI.add_constraint(
        optimizer,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(n), x), 0.0),
        MOI.LessThan(τ),
    )
    return MathOptLMO(optimizer)
end

"""
Dual costs for a given primal solution to form a primal dual pair
for scaled unit simplex.
Returns two vectors. The first one is the dual costs associated with the constraints
and the second is the reduced costs for the variables.
"""
function compute_dual_solution(::UnitSimplexOracle{T}, direction, primalSolution) where {T}
    idx = argmax(primalSolution)
    critical = min(direction[idx], 0)
    lambda = [critical]
    mu = direction .- lambda
    return lambda, mu
end

"""
    ProbabilitySimplexOracle(right_side)

Represents the scaled probability simplex:
```
C = {x ∈ R^n_+, ∑x = right_side}
```
"""
struct ProbabilitySimplexOracle{T} <: LinearMinimizationOracle
    right_side::T
end

ProbabilitySimplexOracle{T}() where {T} = ProbabilitySimplexOracle{T}(one(T))

ProbabilitySimplexOracle(rhs::Integer) = ProbabilitySimplexOracle{Float64}(rhs)

"""
LMO for scaled probability simplex.
Returns a vector with one active value equal to RHS in the
most improving (or least degrading) direction.
"""
function compute_extreme_point(lmo::ProbabilitySimplexOracle{T}, direction; v = nothing, kwargs...) where {T}
    idx = argmin(direction)
    return ScaledHotVector(lmo.right_side, idx, length(direction))
end

function convert_mathopt(
    lmo::ProbabilitySimplexOracle{T},
    optimizer::OT;
    dimension::Integer,
    kwargs...,
) where {T,OT}
    MOI.empty!(optimizer)
    τ = lmo.right_side
    n = dimension
    (x, _) = MOI.add_constrained_variables(optimizer, [MOI.Interval(0.0, 1.0) for _ in 1:n])
    MOI.add_constraint(
        optimizer,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(n), x), 0.0),
        MOI.EqualTo(τ),
    )
    return MathOptLMO(optimizer)
end

"""
Dual costs for a given primal solution to form a primal dual pair
for scaled probability simplex.
Returns two vectors. The first one is the dual costs associated with the constraints
and the second is the reduced costs for the variables.
"""
function compute_dual_solution(
    ::ProbabilitySimplexOracle{T},
    direction,
    primal_solution;
    kwargs...,
) where {T}
    idx = argmax(primal_solution)
    lambda = [direction[idx]]
    mu = direction .- lambda
    return lambda, mu
end