text:
const strip_patterns                = @compat(UInt32(0))
const strip_corrupt_utf8            = @compat(UInt32(0x1)) << 0
const strip_case                    = @compat(UInt32(0x1)) << 1
const stem_words                    = @compat(UInt32(0x1)) << 2
const tag_part_of_speech            = @compat(UInt32(0x1)) << 3

const strip_whitespace              = @compat(UInt32(0x1)) << 5
const strip_punctuation             = @compat(UInt32(0x1)) << 6
const strip_numbers                 = @compat(UInt32(0x1)) << 7
const strip_non_letters             = @compat(UInt32(0x1)) << 8

const strip_indefinite_articles     = @compat(UInt32(0x1)) << 9
const strip_definite_articles       = @compat(UInt32(0x1)) << 10
const strip_articles                = (strip_indefinite_articles |
                                       strip_definite_articles)

const strip_prepositions            = @compat(UInt32(0x1)) << 13
const strip_pronouns                = @compat(UInt32(0x1)) << 14

const strip_stopwords               = @compat(UInt32(0x1)) << 16
const strip_sparse_terms            = @compat(UInt32(0x1)) << 17
const strip_frequent_terms          = @compat(UInt32(0x1)) << 18

const strip_html_tags               = @compat(UInt32(0x1)) << 20

const alpha_sparse      = 0.05
const alpha_frequent    = 0.95

const regex_cache = Dict{AbstractString, Regex}()
function mk_regex(regex_string)
    d = haskey(regex_cache, regex_string) ?
            regex_cache[regex_string] :
            (regex_cache[regex_string] = Regex(regex_string))
    (length(regex_cache) > 50) && empty!(regex_cache)
    d
end

##############################################################################
#
# Remove corrupt UTF8 characters
#
##############################################################################
function remove_corrupt_utf8(s::AbstractString)
    r = zeros(Char, endof(s)+1)
    i = 0
    for chr in s
        i += 1
        r[i] = (chr != 0xfffd) ? chr : ' '
    end
    return Compat.UTF8String(CharString(r))
end

remove_corrupt_utf8!(d::FileDocument) = error("FileDocument cannot be modified")

function remove_corrupt_utf8!(d::StringDocument)
    d.text = remove_corrupt_utf8(d.text)
    nothing
end

function remove_corrupt_utf8!(d::TokenDocument)
    for i in 1:length(d.tokens)
        d.tokens[i] = remove_corrupt_utf8(d.tokens[i])
    end
end

function remove_corrupt_utf8!(d::NGramDocument)
    new_ngrams = Dict{AbstractString, Int}()
    for token in keys(d.ngrams)
        new_token = remove_corrupt_utf8(token)
        if haskey(new_ngrams, new_token)
            new_ngrams[new_token] = new_ngrams[new_token] + 1
        else
            new_ngrams[new_token] = 1
        end
    end
    d.ngrams = new_ngrams
end

function remove_corrupt_utf8!(crps::Corpus)
    for doc in crps
        remove_corrupt_utf8!(doc)
    end
end

##############################################################################
#
# Conversion to lowercase
#
##############################################################################

remove_case{T <: AbstractString}(s::T) = lowercase(s)

remove_case!(d::FileDocument) = error("FileDocument cannot be modified")

function remove_case!(d::StringDocument)
    d.text = remove_case(d.text)
    nothing
end

function remove_case!(d::TokenDocument)
    for i in 1:length(d.tokens)
        d.tokens[i] = remove_case(d.tokens[i])
    end
end

function remove_case!(d::NGramDocument)
    new_ngrams = Dict{AbstractString, Int}()
    for token in keys(d.ngrams)
        new_token = remove_case(token)
        if haskey(new_ngrams, new_token)
            new_ngrams[new_token] = new_ngrams[new_token] + 1
        else
            new_ngrams[new_token] = 1
        end
    end
    d.ngrams = new_ngrams
end

function remove_case!(crps::Corpus)
    for doc in crps
        remove_case!(doc)
    end
end

##############################################################################
#
# Stripping HTML tags
#
##############################################################################
const script_tags = Regex("<script\\b[^>]*>([\\s\\S]*?)</script>")
const html_tags = Regex("<[^>]*>")

function remove_html_tags(s::AbstractString)
    s = remove_patterns(s, script_tags)
    remove_patterns(s, html_tags)
end

function remove_html_tags!(d::AbstractDocument)
    error("HTML tags can be removed only from a StringDocument")
end

function remove_html_tags!(d::StringDocument)
    d.text = remove_html_tags(d.text)
    nothing
end

function remove_html_tags!(crps::Corpus)
    for doc in crps
        remove_html_tags!(doc)
    end
end

##############################################################################
#
# Remove specified words
#
##############################################################################
function remove_words!{T <: AbstractString}(entity::(@compat Union{AbstractDocument,Corpus}),
                                    words::Vector{T})
    skipwords = Set{AbstractString}()
    union!(skipwords, words)
    prepare!(entity, strip_patterns, skip_words = skipwords)
end

function remove_whitespace!(entity::(@compat Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_whitespace! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_whitespace)
end
function remove_punctuation!(entity::(@compat Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_punctuation! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_punctuation)
end
function remove_nonletters!(entity::(@compat Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_nonletters! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_non_letters)
end
function remove_numbers!(entity::(@compat Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_numbers! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_numbers)
end
function remove_articles!(entity::(@compat Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_articles! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_articles)
end
function remove_indefinite_articles!(entity::(@compat Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_indefinite_articles! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_indefinite_articles)
end
function remove_definite_articles!(entity::(@compat Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_definite_articles! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_definite_articles)
end
function remove_prepositions!(entity::(@compat Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_prepositions! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_prepositions)
end
function remove_pronouns!(entity::(@compat Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_pronouns! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_pronouns)
end
function remove_stop_words!(entity::(@compat Union{AbstractDocument,Corpus}))
    Base.warn_once("remove_stop_words! is deprecated, Use prepare! instead.")
    prepare!(entity, strip_stopwords)
end

##############################################################################
#
# Part-of-Speech tagging
#
##############################################################################

tag_pos!(entity) = error("Not yet implemented")

##############################################################################
#
# Drop terms based on frequency
#
##############################################################################

function sparse_terms(crps::Corpus, alpha::Real = alpha_sparse)
    update_lexicon!(crps)
    update_inverse_index!(crps)
    res = Array(UTF8String, 0)
    ndocs = length(crps.documents)
    for term in keys(crps.lexicon)
        f = length(crps.inverse_index[term]) / ndocs
        if f <= alpha
            push!(res, Compat.UTF8String(term))
        end
    end
    return res
end

function frequent_terms(crps::Corpus, alpha::Real = alpha_frequent)
    update_lexicon!(crps)
    update_inverse_index!(crps)
    res = Array(UTF8String, 0)
    ndocs = length(crps.documents)
    for term in keys(crps.lexicon)
        f = length(crps.inverse_index[term]) / ndocs
        if f >= alpha
            push!(res, Compat.UTF8String(term))
        end
    end
    return res
end

# Sparse terms occur in less than x percent of all documents
remove_sparse_terms!(crps::Corpus, alpha::Real = alpha_sparse) = remove_words!(crps, sparse_terms(crps, alpha))

# Frequent terms occur in more than x percent of all documents
remove_frequent_terms!(crps::Corpus, alpha::Real = alpha_frequent) = remove_words!(crps, frequent_terms(crps, alpha))

##############################################################################
#
# Remove parts from document based on flags or regular expressions
#
##############################################################################

function prepare!(crps::Corpus, flags::UInt32; skip_patterns = Set{AbstractString}(), skip_words = Set{AbstractString}())
    ((flags & strip_sparse_terms) > 0) && union!(skip_words, sparse_terms(crps))
    ((flags & strip_frequent_terms) > 0) && union!(skip_words, frequent_terms(crps))

    ((flags & strip_corrupt_utf8) > 0) && remove_corrupt_utf8!(crps)
    ((flags & strip_case) > 0) && remove_case!(crps)
    ((flags & strip_html_tags) > 0) && remove_html_tags!(crps)

    lang = language(crps.documents[1])   # assuming all documents are of the same language - practically true
    r = _build_regex(lang, flags, skip_patterns, skip_words)
    !isempty(r.pattern) && remove_patterns!(crps, r)

    ((flags & stem_words) > 0) && stem!(crps)
    ((flags & tag_part_of_speech) > 0) && tag_pos!(crps)
    nothing
end

function prepare!(d::AbstractDocument, flags::UInt32; skip_patterns = Set{AbstractString}(), skip_words = Set{AbstractString}())
    ((flags & strip_corrupt_utf8) > 0) && remove_corrupt_utf8!(d)
    ((flags & strip_case) > 0) && remove_case!(d)
    ((flags & strip_html_tags) > 0) && remove_html_tags!(d)

    r = _build_regex(language(d), flags, skip_patterns, skip_words)
    !isempty(r.pattern) && remove_patterns!(d, r)

    ((flags & stem_words) > 0) && stem!(d)
    ((flags & tag_part_of_speech) > 0) && tag_pos!(d)
    nothing
end

function remove_patterns(s::AbstractString, rex::Regex)
    iob = IOBuffer()
    ibegin = 1
    for m in matchall(rex, s)
        len = m.offset-ibegin+1
        if len > 0
            Base.write_sub(iob, s.data, ibegin, len)
            write(iob, ' ')
        end
        ibegin = m.endof+m.offset+1
    end
    len = length(s.data) - ibegin + 1
    (len > 0) && Base.write_sub(iob, s.data, ibegin, len)
    takebuf_string(iob)
end

function remove_patterns{T <: Compat.String}(s::SubString{T}, rex::Regex)
    iob = IOBuffer()
    ioffset = s.offset
    data = s.string.data
    ibegin = 1
    for m in matchall(rex, s)
        len = m.offset-ibegin+1
        if len > 0
            Base.write_sub(iob, data, ibegin+ioffset, len)
            write(iob, ' ')
        end
        ibegin = m.endof+m.offset+1
    end
    len = s.endof - ibegin + 1
    (len > 0) && Base.write_sub(iob, data, ibegin+ioffset, len)
    takebuf_string(iob)
end

remove_patterns!(d::FileDocument, rex::Regex) = error("FileDocument cannot be modified")

function remove_patterns!(d::StringDocument, rex::Regex)
    d.text = remove_patterns(d.text, rex)
    nothing
end

function remove_patterns!(d::TokenDocument, rex::Regex)
    for i in 1:length(d.tokens)
        d.tokens[i] = remove_patterns(d.tokens[i], rex)
    end
end

function remove_patterns!(d::NGramDocument, rex::Regex)
    new_ngrams = Dict{AbstractString, Int}()
    for token in keys(d.ngrams)
        new_token = remove_patterns(token, rex)
        if haskey(new_ngrams, new_token)
            new_ngrams[new_token] = new_ngrams[new_token] + 1
        else
            new_ngrams[new_token] = 1
        end
    end
    d.ngrams = new_ngrams
    nothing
end

function remove_patterns!(crps::Corpus, rex::Regex)
    for doc in crps
        remove_patterns!(doc, rex)
    end
end

##
# internal helper methods

_build_regex(lang, flags::UInt32) = _build_regex(lang, flags, Set{AbstractString}(), Set{AbstractString}())
_build_regex{T <: AbstractString}(lang, flags::UInt32, patterns::Set{T}, words::Set{T}) = _combine_regex(_build_regex_patterns(lang, flags, patterns, words))

function _combine_regex{T <: AbstractString}(regex_parts::Set{T})
    l = length(regex_parts)
    (0 == l) && return r""
    (1 == l) && return mk_regex(pop!(regex_parts))

    iob = IOBuffer()
    write(iob, "($(pop!(regex_parts)))")
    for part in regex_parts
        write(iob, "|($part)")
    end
    mk_regex(takebuf_string(iob))
end

function _build_regex_patterns{T <: AbstractString}(lang, flags::UInt32, patterns::Set{T}, words::Set{T})
    ((flags & strip_whitespace) > 0) && push!(patterns, "\\s+")
    if (flags & strip_non_letters) > 0
        push!(patterns, "[^a-zA-Z\\s]")
    else
        ((flags & strip_punctuation) > 0) && push!(patterns, "[,;:.!?()]+")
        ((flags & strip_numbers) > 0) && push!(patterns, "\\d+")
    end
    if (flags & strip_articles) > 0
        union!(words, articles(lang))
    else
        ((flags & strip_indefinite_articles) > 0) && union!(words, indefinite_articles(lang))
        ((flags & strip_definite_articles) > 0) && union!(words, definite_articles(lang))
    end
    ((flags & strip_prepositions) > 0) && union!(words, prepositions(lang))
    ((flags & strip_pronouns) > 0) && union!(words, pronouns(lang))
    ((flags & strip_stopwords) > 0) && union!(words, stopwords(lang))

    words_pattern = _build_words_pattern(collect(words))
    !isempty(words_pattern) && push!(patterns, words_pattern)
    patterns
end

function _build_words_pattern{T <: AbstractString}(words::Vector{T})
    isempty(words) && return ""

    iob = IOBuffer()
    write(iob, "\\b(")
    write(iob, words[1])
    l = length(words)
    for idx in 2:l
        write(iob, '|')
        write(iob, words[idx])
    end
    write(iob, ")\\b")
    takebuf_string(iob)
end


text:
using GHOST
using GHOST: @unpack, @everywhere, READY, remotecall, load!
setup()
time_start = now()
println(time_start)
setup_parallel()

@unpack conn, schema, pat = GHOST.PARALLELENABLER
data = execute(conn,
               "SELECT branch FROM $(schema).repos WHERE status = 'Init' ORDER BY commits;",
               not_null = true) |>
    (obj -> getproperty.(obj, :branch))
data = [ @view(data[i:min(i + 99, lastindex(data))]) for i in 1:100:length(data) ]
function magic(nodes::AbstractVector)
    query = "fragment X on Ref { id target { ... on Commit { history(until: \$until) { totalCount } } } } query Commits(\$x:[ID!]!, \$until:GitTimestamp!){nodes(ids:\$x){...X}}"
    vars = Dict("x" => nodes, "until" => "2020-01-01T00:00:00Z")
    result = graphql(query, vars = vars)
    json = JSON3.read(result.Data)
    try
        if :data ∈ propertynames(json)
            x = DataFrame((;x.id, x.target.history.totalCount) for x in (x for x in json.data.nodes if !isnothing(x)))
            load!(x, GHOST.PARALLELENABLER.conn, "INSERT INTO gh_2007_2020.repos_chk VALUES(\$1,\$2) ON CONFLICT DO NOTHING;")
        end
        if :errors ∈ propertynames(json)
            y = [ (branch = SubString(x.message, 52, length(x.message) - 2), commits = 0) for x in values(json.errors) ]
            load!(y, GHOST.PARALLELENABLER.conn, "INSERT INTO gh_2007_2020.repos_chk VALUES(\$1,\$2) ON CONFLICT DO NOTHING;")
        end
    catch err
        println("Error in: $nodes")
    end
    nothing
end
magic(data[2])
@sync @distributed for idx in eachindex(data)
    magic(data[idx])
    println(idx)
end

ext: jl


text:
using LinearAlgebra

function E(r :: Vector{Float64}, t :: Float64) :: Vector{Float64}
    [0., 0., .01]
end

function B(r :: Vector{Float64}, t :: Float64) :: Vector{Float64}
    [0., 0., 1.]
end

function boris_step!(r :: Vector{Float64}, v :: Vector{Float64},
                     qm :: Float64, dt :: Float64, t :: Float64) :: Tuple{Vector, Vector}
    dtqm = dt * qm
    r += .5 * dt .* v
    B_rt = B(r, t)
    E_rt = E(r, t)
    p = .5 * dtqm .* B_rt
    a_sq = .25 * dtqm^2 * dot(B_rt, B_rt)
    v += .5 * dtqm .* E_rt
    v_prime = v .+ cross(v, p)
    v += 2. .* cross(v_prime, p) / (1 + a_sq)
    v += .5 * dtqm .* E_rt
    r += .5 * dt .* v
    return r, v
end

qm = 1.
dt = .25
steps = 100
r = zeros(steps, 3)
v = zeros(steps, 3)
r[1,:] = [0. 0. 0.]
v[1,:] = [.5 0. 0.]

function test!()
    t = 0.
    for i=2:steps
        t += .5 * dt
        next = boris_step!(r[i-1,:], v[i-1,:], qm, dt, t)
        r[i,:] = next[1]
        v[i,:] = next[2]
        t += .5 * dt
    end
    return r, v
end


text:
"""
Creates an expression with the same head and typ
"""
function similar_expr(x::Expr, args)
    expr = Expr(x.head)
    expr.typ = x.typ
    expr.args = args
    expr
end
similar_expr(x::Expr) = similar_expr(x, [])

function typed_expr(T, head, args...)
    expr = Expr(head, args...)
    expr.typ = T
    expr
end
function filter_expr(keep, ast)
    replace_or_drop(x-> (false, x), x-> !keep(x), ast)
end
function replace_expr(f, ast)
    replace_or_drop(f, x-> false, ast) # never drop/filter
end
function replace_or_drop(f, drop, ast::Vector, result = [])
    for elem in ast
        replace_or_drop(f, drop, elem, result)
    end
    result
end

function replace_or_drop(f, drop, ast, result = [])
    drop(ast) && return result
    replace, replacement = f(ast)
    if replace
        if isa(replacement, Tuple)
            if !isempty(replacement)
                push!(result, replacement...)
            end
        else
            push!(result, replacement)
        end
    else
        expr = if isa(ast, Expr)
            nexpr = similar_expr(ast)
            replace_or_drop(f, drop, ast.args, nexpr.args)
            nexpr
        else
            ast
        end
        push!(result, expr)
    end
    result
end

isa_applytype(x) = false
# function isa_applytype(x::Expr)
#     x.head == :. && x.args[1] == :Core && x.args[2] == QuoteNode(:apply_type)
# end
function isa_applytype(x::GlobalRef)
    x.mod == Core && x.name == :apply_type
end
function isa_applytype(x::Expr)
    x.head == :call || return false
    isa_applytype(x.args[1])
end

applytype_args(x::Expr) = x.args[2:end]

function applytype_type(x::Expr, args)
    Targs = x.args[2:end]
    Expr(:curly, Targs...)
end
function applytype_type(x::GlobalRef, args::Vector)
    Expr(:curly,  args...)
end

_normalize_ast(value) = true, value
function _normalize_ast(qn::QuoteNode)
    return true, qn.value
end

function resolve_typ(x::Expr)
    x.typ.parameters[1]
end
resolve_typ(x::GlobalRef) = eval(x)

function _normalize_ast(expr::Expr)
    if expr.head == :invoke
        lam = expr.args[1] # Ignore lambda for now
        res = similar_expr(expr, map(normalize_ast, view(expr.args, 2:length(expr.args))))
        return true, res
    # elseif expr.head == :new
    #     T = resolve_typ(expr.args[1])
    #     expr = Expr(:call, T, map(normalize_ast, expr.args[2:end])...)
    #     expr.typ = T
    #     return true, expr
    # elseif expr.head == :static_parameter# TODO do something reasonable with static and meta
    #     # TODO, can other static parameters beside literal values escape with code_typed, optimization = false?
    #     return true, expr.args[1]
    elseif expr.head == :meta
        return true, nothing
    elseif expr.head == :call
        f = expr.args[1]
        if f == GlobalRef(Core, :apply_type)
            T = unspecialized_type(expr.typ)
            return true, T
        end
        return true, similar_expr(expr, map(normalize_ast, expr.args))
    end
    return false, expr
end
function normalize_ast(expr)
    ast = Sugar.replace_expr(_normalize_ast, expr)
    Sugar.replace_or_drop(x-> (false, x), x-> x == nothing, ast)
end

function remove_inline(list)
    stack = []; result = []
    for elem in list
        if isa(elem, Expr) && elem.head == :meta
            if elem.args[1] == :push_loc
                inexpr = Expr(:inlined_func, elem.args[2:end]...)
                if !isempty(stack)
                    push!(last(stack).args, inexpr) # add to parent
                end
                push!(stack, inexpr)
            elseif elem.args[1] == :pop_loc
                if isempty(stack)
                    error("found pop_loc, but there was no push_loc")
                end
                expr = pop!(stack)
                if isempty(stack)
                    push!(result, expr)
                end
            end
        elseif !isempty(stack) # we're in the middle of an inline expr
            push!(last(stack).args, elem)
        else
            push!(result, elem)
        end
    end
    result
end


text:
abstract type Model{T} <: AbstractDataset{T}
end

users(model::Model) = model.users
items(model::Model) = model.items

checkuser(model::Model, user::Int) = (user >= 0 && user <= users(model)) ? user : error("invalid user id ($user)")
checkitem(model::Model, item::Int) = (item >= 0 && item <= items(model)) ? item : error("invalid item id ($item)")

predict(model::Model, user::Int, item::Int) where T = missing

train!(model::Model, dataset::Dataset) = nothing

function Base.getindex(model::Model{T}, user::Int, item::Int) where T
    checkuser(model, user)
    checkitem(model, item)

    value = predict(model, user, item)

    if !ismissing(value)
        return PredictRating(correct(value, model.preference), model.preference)
    end

    return MissingRating{T}()
end

function Base.getindex(model::Model{T}, dataset::Dataset{T}) where T
    predicts = Vector{UserPreference{T}}(undef, length(dataset))

    for i = 1:length(dataset)
        (u, v, r) = dataset[i]
        predicts[i] = UserPreference{T}(u, v, model[u, v])
    end

    return predicts
end


text:
# Grid test
using CommonLH, Test

function linear_test_grid(n :: Integer)
    return LinearGrid(-3.0, 4.0, n);
end

function power_spaced_test_grid(n :: Integer)
    return PowerSpacedGrid(-3.0, 4.0, n, 0.5)
end

function lin_spaced_test_grid(n :: Integer)
    return LinSpacedGrid(-3.0, 4.0, n, -0.2)
end

function general_test(g :: AbstractGrid)
    @testset "General" begin
        @test length(g) > 1
        @test minimum(g) < maximum(g)
        gridV = grid(g);
        @test length(gridV) == length(g)
        @test all(diff(gridV) .> 0)
        @test gridV[end] ≈ maximum(g)
        @test gridV[1] ≈ minimum(g)
        @test g[1] ≈ gridV[1]
        @test g[end] ≈ gridV[end]
        @test all(intervals(g) .≈ diff(gridV))
        @test sum(intervals(g)) ≈ maximum(g) - minimum(g)
        for (ig, gPoint) in enumerate(g)
            @test gPoint ≈ gridV[ig]
        end

        for ig = 1 : length(g)
            @test g[ig] ≈ gridV[ig]
        end

    end
end

function linear_test(g :: LinearGrid)
    @testset "Linear" begin
        stepV = intervals(g);
        @test all(stepV .≈ stepV[1])
    end
end

@testset "Grids" begin
    n = 5;
    for g in [linear_test_grid(n), power_spaced_test_grid(n),
        lin_spaced_test_grid(n)]

        general_test(g)
    end

    linear_test(linear_test_grid(n));
end

# ---------------


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function julia_logo_test(fname)
    Drawing(1000, 1000, fname)
    origin()
    background("white")
    for (pos, n) in Tiler(1000, 1000, 2, 2)
        gsave()
        translate(pos - Point(150, 100))
        if n == 1
            julialogo()
            fillpath()
        elseif n == 2
            randomhue()
            setline(0.3)
            julialogo(action=:stroke)
        elseif n == 3
            sethue("orange")
            julialogo(color=false)
        elseif n == 4
            julialogo(action=:clip)
            setopacity(0.6)
            for i in 1:400
                randomhue()
                gsave()
                box(Point(rand(0:250), rand(0:250)), 350, 5, :fill)
                grestore()
            end
            clipreset()
        end
        grestore()
    end
    @test finish() == true
end

fname="julia-logo-drawing.pdf"
julia_logo_test(fname)
println("...finished test: output in $(fname)")


text:
using Distributions
using LinearAlgebra
#the following are the initializer that can be used
function glorot_uniform(fan_in, fan_out)
    sd = sqrt(6.0 / fan_in + fan_out)
    m = rand(Uniform(-sd, sd), fan_in,fan_out)
    return m
end

function glorot_normal(fan_in, fan_out)
    sd = sqrt(2.0 / fan_in + fan_out)
    m = rand(Uniform(0, sd), fan_in,fan_out)
    return m
end

#The “xavier” weight initialization was found to have problems when used
#to initialize networks that use the rectified linear (ReLU) activation function.
#Added he
function he_normal(fan_in, fan_out)
    sd = sqrt(2.0 / fan_in)
    m = rand(Uniform(0, (1000*sd)), fan_in,fan_out)
end


text:
# ---
# title: 966. Vowel Spellchecker
# id: problem966
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table, String
# link: <https://leetcode.com/problems/vowel-spellchecker/description/>
# hidden: true
# ---
#
# Given a `wordlist`, we want to implement a spellchecker that converts a query
# word into a correct word.
#
# For a given `query` word, the spell checker handles two categories of spelling
# mistakes:
#
#   * Capitalization: If the query matches a word in the wordlist ( **case-insensitive** ), then the query word is returned with the same case as the case in the wordlist.
#     * Example: `wordlist = ["yellow"]`, `query = "YellOw"`: `correct = "yellow"`
#     * Example: `wordlist = ["Yellow"]`, `query = "yellow"`: `correct = "Yellow"`
#     * Example: `wordlist = ["yellow"]`, `query = "yellow"`: `correct = "yellow"`
#   * Vowel Errors: If after replacing the vowels ('a', 'e', 'i', 'o', 'u') of the query word with any vowel individually, it matches a word in the wordlist ( **case-insensitive** ), then the query word is returned with the same case as the match in the wordlist.
#     * Example: `wordlist = ["YellOw"]`, `query = "yollow"`: `correct = "YellOw"`
#     * Example: `wordlist = ["YellOw"]`, `query = "yeellow"`: `correct = ""` (no match)
#     * Example: `wordlist = ["YellOw"]`, `query = "yllw"`: `correct = ""` (no match)
#
# In addition, the spell checker operates under the following precedence rules:
#
#   * When the query exactly matches a word in the wordlist ( **case-sensitive** ), you should return the same word back.
#   * When the query matches a word up to capitlization, you should return the first such match in the wordlist.
#   * When the query matches a word up to vowel errors, you should return the first such match in the wordlist.
#   * If the query has no matches in the wordlist, you should return the empty string.
#
# Given some `queries`, return a list of words `answer`, where `answer[i]` is
# the correct word for `query = queries[i]`.
#
#
#
# **Example 1:**
#
#
#
#     Input: wordlist = ["KiTe","kite","hare","Hare"], queries = ["kite","Kite","KiTe","Hare","HARE","Hear","hear","keti","keet","keto"]
#     Output: ["kite","KiTe","KiTe","Hare","hare","","","KiTe","","KiTe"]
#
#
#
# **Note:**
#
#   * `1 <= wordlist.length <= 5000`
#   * `1 <= queries.length <= 5000`
#   * `1 <= wordlist[i].length <= 7`
#   * `1 <= queries[i].length <= 7`
#   * All strings in `wordlist` and `queries` consist only of **english**  letters.
#
#
## @lc code=start
using LeetCode

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


text:
n = 10 # rand(1:10)
A = ones(Float64, n, n)
@test matrixdepot("circul", ones(Float64,n)) ≈ A

x = [1:n;];
A = matrixdepot("circul", x)
# test symmetry
@test A + A' == (A + A')'
B = matrixdepot("circul", n)
C = matrixdepot("circul", [1:n;], n)
@test B ≈ C

println("'circul' passed test..." )


text:
"""

```
Trajectory
```

Trajectory data type.
Default reading with the Chemfiles infrastructure, except for DCD and PDB trajectory
files, if the "PDBTraj" option is provided.
See memory issue (https://github.com/chemfiles/Chemfiles.jl/issues/44)

"""
abstract type Trajectory end

function Trajectory(
    filename::String,
    solute::Selection,
    solvent::Selection;
    format::String = "",
    chemfiles = false,
)
    if !chemfiles
        if format == "dcd" || FileOperations.file_extension(filename) == "dcd"
            NamdDCD(filename, solute, solvent)
        elseif format == "PDBTraj"
            PDBTraj(filename, solute, solvent)
        else
            ChemFile(filename, solute, solvent, format = format)
        end
    else
        ChemFile(filename, solute, solvent, format = format)
    end
end

# If only one selection is provided, assume that the solute and the
# solvent are the same

function Trajectory(
    filename::String,
    solvent::Selection;
    format::String = "",
    chemfiles = false,
)
    Trajectory(filename, solvent, solvent, format = format, chemfiles = chemfiles)
end


text:
# This is part of Quiz.jl. License is MIT.

"""
`struct qProb <: qSing`\n
`Quiz.jl`'s single problem container.
"""
struct qProb <: qSing
    val::Vector{qElem}
end

"""
`(x::qProb)()`\n
`qProb` functor that returns its value as a `Vector{qElem}`.
"""
(x::qProb)() = x.val

the(x::qProb) = x()

export qProb


text:
using CSV, DataFrames
deaths = CSV.read(normpath("data", "2018-10-17-4Amod.csv"))
ages = DataFrame(age = collect(0:20),
	agest = [0;0;1;collect(5:5:90)],
	age2 =[0;1;fill(2,3);fill(3,6);fill(4,4);fill(5,2);fill(6,4)])
ages2 = by(ages[:,[:agest, :age2]], :age2,
df -> DataFrame(agest=minimum(df[:agest])))

pop = DataFrame()
pop1 = CSV.read(normpath("data", "mfsv17.csv"))
ages1 = DataFrame(age1 = collect(0:101),
	age = [0;1;fill(2,4);map((x)->fld(x,5)+2, collect(5:89));fill(20,11)])
pop1_age = join(pop1, ages1, on = :age1)

function popsexadd(sex)
	pop1_sex = pop1_age[pop1_age[:sex].==sex, [:n, :age]]
	pop_sex = by(pop1_sex, :age, df -> DataFrame(sex = sex, n = sum(df[:n])))
end

pop = vcat(popsexadd("Kv"), popsexadd("M"))

function df_reg(sex, caexpr, ageaggr = false)
	if caexpr == "pop"
		df_sub = pop[pop[:sex].==sex, [:age, :n]]
		sort!(df_sub, (:age))
	else
		df_sub = deaths[((map((x)->occursin(Regex(caexpr), x),
		convert(Vector{String}, deaths[:icd10]))) .& (deaths[:sex].==sex)),
		[:age, :n]]
	end
	df_aggr = by(df_sub, :age, df -> DataFrame(n=sum(df[:n])))
	df_aggr_age = join(df_aggr, ages, on = :age)
	if ageaggr
		df_aggr_age2 = by(df_aggr_age[:, [:n, :age2]], :age2,
		df -> DataFrame(n=sum(df[:n])))
		return join(df_aggr_age2, ages2, on = :age2)
	else
		return df_aggr_age
	end
end

function propframe(sex, ca1, ca2 = "A00-Y98", ageaggr = false)
	ca1frame = df_reg(sex, ca1, ageaggr)
	ca2frame = df_reg(sex, ca2, ageaggr)
	propframe = DataFrame(age = ca1frame[1], agest = ca1frame[:agest],
			      ca1n = ca1frame[:n],
			      ca2n = ca2frame[:n],
			      prop = ca1frame[:n]./ca2frame[:n])
	return Dict(:sex => sex, :ca1 => ca1, :ca2 => ca2, :frame => propframe)
end

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
### Gradient-based optimization tests
# The following NLopt gradient -based algorithms are compared

grad_algorithms_autodiff = [
    # good
    #:LD_MMA, # Method of Moving Asymptotes #tmp removed from tests
    (algorithm = :LD_SLSQP, skip = []), # Sequential Least-Squares Quadratic Programming
    (algorithm = :LD_CCSAQ, skip = []), # Conservative convex separable approximation
    # errors
    (algorithm = :LD_LBFGS, skip = [:f_3p_1im, :f_4p_2im, :f_4p_3im, :f_5p_3im, :f_2p]), # Low-storage BFGS
    (algorithm = :LD_TNEWTON_PRECOND_RESTART, skip = [:f_3p_1im, :f_4p_2im, :f_4p_3im, :f_5p_3im, :f_2p]), # Preconditioned truncated Newton
    (algorithm = :LD_TNEWTON_PRECOND, skip = [:f_3p_1im, :f_4p_3im, :f_5p_3im]), # Same without restarting
    (algorithm = :LD_TNEWTON_RESTART, skip = [:f_3p_1im, :f_4p_2im, :f_4p_3im, :f_2p]), # Same without preconditioning
    (algorithm = :LD_TNEWTON, skip = [:f_3p_1im, :f_4p_3im]), # Same without restarting or preconditioning
    (algorithm = :LD_VAR2, skip = [:f_3p_1im, :f_4p_2im, :f_4p_3im, :f_5p_3im, :f_2p]), # Shifted limited-memory variable-metric (rank 2)
    (algorithm = :LD_VAR1, skip = [:f_3p_1im, :f_4p_2im, :f_4p_3im, :f_5p_3im, :f_2p])  # Shifted limited-memory variable-metric (rank 1)
]

[test_alg(alg; bounds=(-1e10,1e10), loss_grad=:AUTODIFF) for alg in grad_algorithms_autodiff]

grad_algorithms_finite = [
    # good
    #:LD_MMA, # Method of Moving Asymptotes #tmp removed from tests
    (algorithm = :LD_SLSQP, skip = []), # Sequential Least-Squares Quadratic Programming
    (algorithm = :LD_CCSAQ, skip = [:f_5p_3im]), # Conservative convex separable approximation
    # errors
    (algorithm = :LD_LBFGS, skip = [:f_1p, :f_3p_1im, :f_4p_3im, :f_5p_3im, :f_4p_2im, :f_2p_1im, :f_2p]), # Low-storage BFGS
    (algorithm = :LD_TNEWTON_PRECOND_RESTART, skip = [:f_3p_1im, :f_3p_1im_dep, :f_4p_3im, :f_5p_3im]), # Preconditioned truncated Newton
    (algorithm = :LD_TNEWTON_PRECOND, skip = [:f_3p_1im, :f_4p_3im, :f_5p_3im]), # Same without restarting
    (algorithm = :LD_TNEWTON_RESTART, skip = [:f_3p_1im, :f_4p_3im, :f_5p_3im]), # Same without preconditioning
    (algorithm = :LD_TNEWTON, skip = [:f_3p_1im, :f_4p_3im, :f_5p_3im]), # Same without restarting or preconditioning
    (algorithm = :LD_VAR2, skip = [:f_3p_1im, :f_4p_2im, :f_4p_3im, :f_5p_3im, :f_2p]), # Shifted limited-memory variable-metric (rank 2)
    (algorithm = :LD_VAR1, skip = [:f_3p_1im, :f_4p_2im, :f_4p_3im, :f_5p_3im, :f_2p])  # Shifted limited-memory variable-metric (rank 1)
]

[test_alg(alg; bounds=(-1e10,1e10), loss_grad=:FINITE) for alg in grad_algorithms_finite]


text:
module SVGMesh

using Mesh
using SVG

function meshXY(fid, n::Net, styles::Dict)
	for f in n.faces
		polyline(fid, n.vertices[[f.ab.from, f.ab.to]], styles[f.ab.style])
		polyline(fid, n.vertices[[f.bc.from, f.bc.to]], styles[f.bc.style])
		polyline(fid, n.vertices[[f.ca.from, f.ca.to]], styles[f.ca.style])
	end
end

###
end

max_stars_repo_name: lawless-m/jsvg


text:
export AbstractInnerProduct, EuclideanIP, WeightedIP, euclidean_distance, euclidean_norm,
       init_weight!, change_weight!, norm2, inner, distance

abstract type AbstractInnerProduct end

struct EuclideanIP <: AbstractInnerProduct end

"""
    WeightedIP(weight)

Represents a weighted euclidean inner product.
"""
struct WeightedIP{T<:Real} <: AbstractInnerProduct
    weight::Vector{T} # The applied weights are 1 / w[i]
end
WeightedIP(x::AbstractVector) = WeightedIP(zeros(real(eltype(x)), length(x)))

"""
    distance(u, v, inner_product=EuclideanIP)

Compute the distance ||u-v|| with the norm induced by `inner_product`.
If `inner_product === nothing` the euclidean scalar product is used.
"""
function distance(x::AbstractVector, y::AbstractVector, ::EuclideanIP=EuclideanIP())
    @boundscheck length(x) == length(y)
    n = length(x)
    @inbounds d = abs2(x[1] - y[1])
    @inbounds for i in 2:n
        @fastmath d += abs2(x[i] - y[i])
    end
    sqrt(d)
end
function distance(x::AbstractVector, y::AbstractVector, in::WeightedIP)
    @boundscheck length(in.weight) == length(x) == length(y)
    w = in.weight
    @inbounds d = abs2(x[1] - y[1]) / (w[1]^2)
    for i in 2:length(x)
        @inbounds d += abs2(x[i] - y[i]) / (w[i]^2)
    end
    sqrt(d)
end

function norm2(x::AbstractVector, ::EuclideanIP=EuclideanIP())
    out = zero(real(eltype(x)))
    for i in eachindex(x)
        @inbounds out += abs2(x[i])
    end
    out
end
function norm2(x::AbstractVector, in::WeightedIP)
    @boundscheck length(in.weight) == length(x)
    out = zero(real(eltype(x)))
    w = in.weight
    for i in eachindex(x)
        @inbounds out += abs2(x[i]) / (w[i]^2)
    end
    out
end
LinearAlgebra.norm(x::AbstractVector, ip::AbstractInnerProduct) = sqrt(norm2(x, ip))

(ip::EuclideanIP)(x::AbstractVector) = LinearAlgebra.norm(x, ip)
(ip::WeightedIP)(x::AbstractVector) = LinearAlgebra.norm(x, ip)


text:
#
# diff.jl -
#
# Implement finite differences operators.
#
#-------------------------------------------------------------------------------
#
# This file is part of LazyAlgebra (https://github.com/emmt/LazyAlgebra.jl)
# released under the MIT "Expat" license.
#
# Copyright (c) 2017-2021 Éric Thiébaut.
#

module FiniteDifferences

export Diff

using MayOptimize
using LazyAlgebra
using LazyAlgebra.Foundations
import LazyAlgebra: apply!, vcreate, identical

using Base: @propagate_inbounds
import Base: show

const ArrayAxis = AbstractUnitRange{Int}
const ArrayAxes{N} = NTuple{N,ArrayAxis}

"""
    limits(r) -> (first(r), last(r))

yields the first and last value of the unit-range `r`.

"""
limits(r::AbstractUnitRange) = (first(r), last(r))

"""
    Diff([opt::MayOptimize.Vectorize,] n=1, dims=:)

yields a linear mapping that computes a finite difference approximation of the
`n`-order derivative along the dimension(s) specified by `dims`.  Arguments
`dims` is an integer, a tuple or a vector of integers specifying along which
dimension(s) to apply the operator or `:` to specify all dimensions.  If
multiple dimensions are specified, the result is as if the operator is applied
separately on the specified dimension(s).

Optional argument `opt` is the optimization level and may be specified as the
first or last argument.  By default, `opt` is assumed to be `Vectorize`,
however depending on the dimensions of the array, the dimensions of interest
and on the machine, setting `opt` to `InBounds` may be more efficient.

If `dims` is a scalar, the result, say `y`, of applying the finite difference
operator to an array, say `x`, has the same axes as `x`.  Otherwise and even
though `x` has a single dimension or `dims` is a 1-tuple, `y` has one more
dimension than `x`, the last dimension of `y` is used to store the finite
differences along each dimensions specified by `dims` and the leading
dimensions of `y` are the same as the dimensions of `x`.

More specifically, the operator created by `Diff` implements **forward finite
differences** with **flat boundary conditions**, that is to say extrapolated
entries are assumed equal to the nearest entry.

"""
struct Diff{L,D,O<:OptimLevel} <: LinearMapping end
# L = level of differentiation
# D = list of dimensions along which compute the differences
# O = optimization level

# Constructors.
function Diff(n::Integer = 1,
              dims::Union{Colon,Integer,Tuple{Vararg{Integer}},
                          AbstractVector{<:Integer}}=Colon(),
              opt::Type{<:OptimLevel} = Vectorize)
    return Diff{to_int(n), to_dims(dims), opt}()
end

function Diff(opt::Type{<:OptimLevel}, n::Integer = 1,
              dims::Union{Colon,Integer,Tuple{Vararg{Integer}},
                          AbstractVector{<:Integer}}=Colon())
    return Diff{to_int(n), to_dims(dims), opt}()
end

function Diff(n::Integer, opt::Type{<:OptimLevel})
    return Diff{to_int(n), Colon, opt}()
end

# Make a finite difference operator callable.
@callable Diff

# Two finite difference operators are identical if they have the same level of
# differentiation and list of dimensions along which compute the differences.
# Their optimization levels may be different.
identical(::Diff{L,D}, ::Diff{L,D}) where {L,D} = true

# Print operator in such a way that is similar to how the operator would be
# created in Julia.
show(io::IO, ::Diff{L,D,Opt}) where {L,D,Opt} =
    print(io, "Diff(", L, ',', (D === Colon ? ":" : D),',',
          (Opt === Debug ? "Debug" :
           Opt === InBounds ? "InBounds" :
           Opt === Vectorize ? "Vectorize" : Opt), ')')

"""
    differentiation_order(A)

yields the differentiation order of finite difference operator `A` (argument
can also be a type).

"""
differentiation_order(::Type{<:Diff{L,D,Opt}}) where {L,D,Opt} = L

"""
    dimensions_of_interest(A)

yields the list of dimensions of interest of finite difference operator `A`
(argument can also be a type).

"""
dimensions_of_interest(::Type{<:Diff{L,D,Opt}}) where {L,D,Opt} = D

"""
    optimization_level(A)

yields the optimization level for applying finite difference operator `A`
(argument can also be a type).

"""
optimization_level(::Type{<:Diff{L,D,Opt}}) where {L,D,Opt} = Opt

for f in (:differentiation_order,
          :dimensions_of_interest,
          :optimization_level)
    @eval begin
        $f(A::Diff) = $f(typeof(A))
        $f(A::Gram{<:Diff}) = $f(typeof(A))
        $f(::Type{<:Gram{T}}) where {T<:Diff} = $f(T)
    end
end

# Convert argument to `Int`.
to_int(x::Int) = x
to_int(x::Integer) = Int(x)

# Convert argument to the type parameter which specifies the list of dimensions
# of interest.
to_dims(::Colon) = Colon
to_dims(x::Int) = x
to_dims(x::Integer) = to_int(x)
to_dims(x::Tuple{Vararg{Int}}) = x
to_dims(x::Tuple{Vararg{Integer}}) = map(to_int, x)
to_dims(x::AbstractVector{<:Integer}) = to_dims((x...,))

# Drop list of dimensions from type to avoid unecessary specializations.
anydims(::Diff{L,D,P}) where {L,D,P} = Diff{L,Any,P}()
anydims(::Gram{Diff{L,D,P}}) where {L,D,P} = gram(Diff{L,Any,P}())

# Applying a separable operator is split in several stages:
#
# 1. Check arguments (so that avoiding bound checking should be safe) and deal
#    with the trivial cases α = 0 or no dimension of interest to apply the
#    operation (to simplify subsequent stages).
#
# 2. If α is non-zero, dispatch on dimension(s) along which to apply the
#    operation and on the specific values of the multipliers α and β.
#
# The second stage may be split in several sub-stages.

# Declare all possible signatures (not using unions) to avoid ambiguities.
for (P,A) in ((:Direct,  :Diff),
              (:Adjoint, :Diff),
              (:Direct,  :(Gram{<:Diff})))
    @eval function apply!(α::Number,
                          P::Type{$P},
                          A::$A,
                          x::AbstractArray,
                          scratch::Bool,
                          β::Number,
                          y::AbstractArray)
        inds, ndims = check_arguments(P, A, x, y)
        if α == 0 || ndims < 1
            # Get rid of this stupid case!
            vscale!(y, β)
        else
            # Call unsafe_apply! to dispatch on the dimensions of interest and on
            # the values of the multipliers.
            unsafe_apply!(α, P, A, x, β, y, inds)
        end
        return y
    end
end

# FIXME: This should not be necessary.
function apply!(α::Number,
                ::Type{<:Adjoint},
                A::Gram{<:Diff},
                x::AbstractArray,
                scratch::Bool,
                β::Number,
                y::AbstractArray)
    apply!(α, Direct, A, x, scratch, β, y)
end

function vcreate(::Type{Direct},
                 A::Diff{L,D,P},
                 x::AbstractArray{T,N},
                 scratch::Bool) where {L,D,P,T,N}
    if D === Colon
        return Array{T}(undef, size(x)..., N)
    elseif isa(D, Tuple{Vararg{Int}})
        return Array{T}(undef, size(x)..., length(D))
    elseif isa(D, Int)
        # if L === 1 && scratch && isa(x, Array)
        #    # First order finite difference along a single dimension.
        #    # Operation could be done in-place but we must preserve
        #    # type-stability.
        #    return x
        #else
        #    return Array{T}(undef, size(x))
        #end
        return Array{T}(undef, size(x))
    else
        error("invalid list of dimensions")
    end
end

function vcreate(::Type{Adjoint},
                 A::Diff{L,D,P},
                 x::AbstractArray{T,N},
                 scratch::Bool) where {L,D,P,T,N}
    # Checking the validity of the argument dimensions is done by applying the
    # opererator.  In-place operation never possible, so ignore the scratch
    # flag.
    if D === Colon || isa(D, Tuple{Vararg{Int}})
        return Array{T}(undef, size(x)[1:N-1])
    elseif isa(D, Int)
        return Array{T}(undef, size(x))
    else
        error("invalid list of dimensions")
    end
end

#------------------------------------------------------------------------------
# CHECKING OF ARGUMENTS

"""
    check_arguments(P, A, x, y) -> inds, ndims

checks that arguments `x` and `y` are valid for applying `P(A)`, with `A` a
separable operator, to `x` and store the result in `y`.  The result is a
2-tuple, `inds` is the axes that the arguments have in common and `ndims` is
the number of dimensions of interest.

If this function returns normally, the caller may safely assume that index
bound checking is not needed; hence, this function must throw an exception if
the dimensions/indices of `x` and `y` are not compatible or if the dimensions
of interest in `A` are out of range.  This function may also throw an exception
if the element types of `x` and `y` are not compatible.

This method must be specialized for the different types of separable operators.

"""
function check_arguments(P::Type{<:Union{Direct,Adjoint}},
                         A::Union{Diff{L,D},Gram{<:Diff{L,D}}},
                         x::AbstractArray,
                         y::AbstractArray) where {L,D}
    inds = check_axes(P, A, axes(x), axes(y))
    ndims = check_dimensions_of_interest(D, length(inds))
    return inds, ndims
end

function check_axes(P::Type{<:Union{Direct,Adjoint}},
                    A::Diff{L,D},
                    xinds::ArrayAxes,
                    yinds::ArrayAxes) where {L,D}
    if D === Colon || isa(D, Dims)
        if P === Direct
            length(yinds) == length(xinds) + 1 ||
                throw_dimension_mismatch("output array must have one more dimension than input array")
            N = (D === Colon ? length(xinds) : length(D))
            yinds[end] == 1:N ||
                throw_dimension_mismatch("last axis of output array must be 1:", N)
            yinds[1:end-1] == xinds ||
                throw_dimension_mismatch("leading axes must be identical")
            return xinds
        else
            length(yinds) == length(xinds) - 1 ||
                throw_dimension_mismatch("output array must have one less dimension than input array")
            N = (D === Colon ? length(yinds) : length(D))
            xinds[end] == 1:N ||
                throw_dimension_mismatch("last axis of input array must be 1:", N)
            xinds[1:end-1] == yinds ||
                throw_dimension_mismatch("leading axes must be identical")
            return yinds
        end
    elseif isa(D, Int)
        xinds == yinds || throw_dimension_mismatch("array axes must be identical")
        return xinds
    else
        throw(ArgumentError("invalid dimensions of interest"))
    end
end

function check_axes(P::Type{<:Operations},
                    A::Gram{<:Diff},
                    xinds::ArrayAxes,
                    yinds::ArrayAxes)
    xinds == yinds || throw_dimension_mismatch("array axes must be identical")
    return xinds
end

check_dimensions_of_interest(::Type{Colon}, ndims::Int) = ndims

check_dimensions_of_interest(dim::Int, ndims::Int) = begin
    1 ≤ dim ≤ ndims ||
        throw_dimension_mismatch("out of range dimension ", dim,
                                 "for ", ndims,"-dimensional arrays")
    return 1
end

check_dimensions_of_interest(dims::Dims{N}, ndims::Int) where {N} = begin
    for dim in dims
        1 ≤ dim ≤ ndims ||
            throw_dimension_mismatch("out of range dimension ", dim,
                                     "for ", ndims,"-dimensional arrays")
    end
    return N
end

throw_dimension_mismatch(str::String) = throw(DimensionMismatch(str))

@noinline throw_dimension_mismatch(args...) =
    throw_dimension_mismatch(string(args...))

#------------------------------------------------------------------------------

# Apply the operation along all dimensions of interest but one dimension at a
# time and knowing that α is not zero.
@generated function unsafe_apply!(α::Number,
                                  ::Type{P},
                                  A::Diff{L,D},
                                  x::AbstractArray,
                                  β::Number,
                                  y::AbstractArray,
                                  inds::ArrayAxes{N}) where {L,D,N,
                                                             P<:Union{Direct,
                                                                      Adjoint}}
    # Allocate empty vector of statements.
    exprs = Expr[]

    # Discard type parameter specifying the dimensions of interest to avoid
    # specialization on this parameter.
    push!(exprs, :(B = anydims(A)))

    # Dispatch on dimensions of interest.
    if isa(D, Int)
        # Arrays x and y have the same dimensions.
        push!(exprs, :(unsafe_apply!(α, P, B, x, β, y,
                                     inds[1:$(D-1)],
                                     inds[$D],
                                     inds[$(D+1):$N],
                                     CartesianIndex())))
    elseif D === Colon || isa(D, Dims)
        # One of x or y (depending on whether the direct or the adjoint
        # operator is applied) has an extra leading dimension used to store the
        # result computed along a given dimension.
        keep_beta = true # initially scale y by β
        dims = (D === Colon ? (1:N) : D)
        for l in 1:length(dims)
            d = dims[l]
            push!(exprs, :(unsafe_apply!(α, P, B, x,
                                         $(keep_beta ? :β : 1), y,
                                         inds[1:$(d-1)],
                                         inds[$d],
                                         inds[$(d+1):$N],
                                         CartesianIndex($l))))
            keep_beta = (P === Direct && A <: Diff)
        end
    else
        # This should never happen.
        return quote
            error("invalid list of dimensions of interest")
        end
    end

    return quote
        $(Expr(:meta, :inline))
        $(exprs...)
        nothing
    end
end

@generated function unsafe_apply!(α::Number,
                                  ::Type{P},
                                  A::Gram{<:Diff{L,D}},
                                  x::AbstractArray,
                                  β::Number,
                                  y::AbstractArray,
                                  inds::ArrayAxes{N}) where {L,D,N,
                                                             P<:Direct}
    # Allocate empty vector of statements.
    exprs = Expr[]

    # Discard type parameter specifying the dimensions of interest to avoid
    # specialization on this parameter.
    push!(exprs, :(B = anydims(A)))

    # Dispatch on dimensions of interest.  Arrays x and y have the same
    # dimensions and there is no last index `l` to specify.
    if isa(D, Int)
        push!(exprs, :(unsafe_apply!(α, P, B, x, β, y,
                                     inds[1:$(D-1)],
                                     inds[$D],
                                     inds[$(D+1):$N])))
    elseif D === Colon || isa(D, Dims)
        # β is set to 1 after first dimension of interest.
        dims = (D === Colon ? (1:N) : D)
        for l in 1:length(dims)
            d = dims[l]
            push!(exprs, :(unsafe_apply!(α, P, B, x,
                                         $(l == 1 ? :β : 1), y,
                                         inds[1:$(d-1)],
                                         inds[$d],
                                         inds[$(d+1):$N])))
        end
    else
        # This should never happen.
        return quote
            error("invalid list of dimensions of interest")
        end
    end

    return quote
        $(Expr(:meta, :inline))
        $(exprs...)
        nothing
    end
end

# Dispatch on multipliers values (α is not zero).
function unsafe_apply!(alpha::Number,
                       P::Type{<:Operations},
                       A::Union{Diff{L,Any,Opt},
                                Gram{Diff{L,Any,Opt}}},
                       x::AbstractArray,
                       beta::Number,
                       y::AbstractArray,
                       I::ArrayAxes,
                       J::ArrayAxis,
                       K::ArrayAxes,
                       l::CartesianIndex) where {L,Opt}
    if alpha == 1
        if beta == 0
            unsafe_apply!(axpby_yields_x,     1, P, A, x, 0, y, I, J, K, l)
        elseif beta == 1
            unsafe_apply!(axpby_yields_xpy,   1, P, A, x, 1, y, I, J, K, l)
        else
            β = promote_multiplier(beta, y)
            unsafe_apply!(axpby_yields_xpby,  1, P, A, x, β, y, I, J, K, l)
        end
    else
        α = promote_multiplier(alpha, y)
        if beta == 0
            unsafe_apply!(axpby_yields_ax,    α, P, A, x, 0, y, I, J, K, l)
        elseif beta == 1
            unsafe_apply!(axpby_yields_axpy,  α, P, A, x, 1, y, I, J, K, l)
        else
            β = promote_multiplier(beta, y)
            unsafe_apply!(axpby_yields_axpby, α, P, A, x, β, y, I, J, K, l)
        end
    end
    nothing
end

# Dispatch on multipliers values (α is not zero) for Gram compositions of a
# finite difference operator.
function unsafe_apply!(alpha::Number,
                       P::Type{<:Operations},
                       A::Gram{<:Diff},
                       x::AbstractArray,
                       beta::Number,
                       y::AbstractArray,
                       I::ArrayAxes,
                       J::ArrayAxis,
                       K::ArrayAxes)
    if alpha == 1
        if beta == 0
            unsafe_apply!(axpby_yields_x,     1, P, A, x, 0, y, I, J, K)
        elseif beta == 1
            unsafe_apply!(axpby_yields_xpy,   1, P, A, x, 1, y, I, J, K)
        else
            β = promote_multiplier(beta, y)
            unsafe_apply!(axpby_yields_xpby,  1, P, A, x, β, y, I, J, K)
        end
    else
        α = promote_multiplier(alpha, y)
        if beta == 0
            unsafe_apply!(axpby_yields_ax,    α, P, A, x, 0, y, I, J, K)
        elseif beta == 1
            unsafe_apply!(axpby_yields_axpy,  α, P, A, x, 1, y, I, J, K)
        else
            β = promote_multiplier(beta, y)
            unsafe_apply!(axpby_yields_axpby, α, P, A, x, β, y, I, J, K)
        end
    end
    nothing
end

#------------------------------------------------------------------------------
#
# The operator D implementing 1st order forward finite difference with flat
# boundary conditions and its adjoint D' are given by:
#
#     D = [ -1   1   0   0
#            0  -1   1   0
#            0   0  -1   1
#            0   0   0   0];
#
#     D' = [ -1   0   0   0
#             1  -1   0   0
#             0   1  -1   0
#             0   0   1   0];
#
# The row (for D) and column (for D') of zeros are to preserve the size.  This
# is needed for multi-dimensional arrays when derivatives along each dimension
# are stored into a single array.
#
# Apply 1st order finite differences along 1st dimension:
#
function unsafe_apply!(f::Function,
                       α::Number,
                       ::Type{Direct},
                       A::Diff{1,Any,Opt},
                       x::AbstractArray,
                       β::Number,
                       y::AbstractArray,
                       I::Tuple{},
                       J::ArrayAxis,
                       K::ArrayAxes,
                       l::CartesianIndex) where {Opt}
    T = real(eltype(x))
    jmin, jmax = limits(J)
    if jmin ≤ jmax
        @maybe_inbounds Opt for k in CartesianIndices(K)
            @maybe_vectorized Opt for j in jmin:jmax-1
                z = x[j+1,k] - x[j,k]
                y[j,k,l] = f(α, z, β, y[j,k,l])
            end
            let j = jmax, z = zero(T)
                y[j,k,l] = f(α, z, β, y[j,k,l])
            end
        end
    end
    nothing
end
#
# Apply 1st order finite differences along 2nd and subsequent dimensions:
#
function unsafe_apply!(f::Function,
                       α::Number,
                       ::Type{Direct},
                       A::Diff{1,Any,Opt},
                       x::AbstractArray,
                       β::Number,
                       y::AbstractArray,
                       I::ArrayAxes,
                       J::ArrayAxis,
                       K::ArrayAxes,
                       l::CartesianIndex) where {Opt}
    T = real(eltype(x))
    jmin, jmax = limits(J)
    if jmin ≤ jmax
        @maybe_inbounds Opt for k in CartesianIndices(K)
            for j in jmin:jmax-1
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    z = x[i,j+1,k] - x[i,j,k]
                    y[i,j,k,l] = f(α, z, β, y[i,j,k,l])
                end
            end
            let j = jmax, z = zero(T)
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k,l] = f(α, z, β, y[i,j,k,l])
                end
            end
        end
    end
    nothing
end
#
# Apply adjoint of 1st order finite differences along 1st dimension:
#
function unsafe_apply!(f::Function,
                       α::Number,
                       ::Type{Adjoint},
                       A::Diff{1,Any,Opt},
                       x::AbstractArray,
                       β::Number,
                       y::AbstractArray,
                       I::Tuple{},
                       J::ArrayAxis,
                       K::ArrayAxes,
                       l::CartesianIndex) where {Opt}
    T = real(eltype(x))
    jmin, jmax = limits(J)
    if jmin < jmax
        @maybe_inbounds Opt for k in CartesianIndices(K)
            let j = jmin
                z = -x[j,k,l]
                y[j,k] = f(α, z, β, y[j,k])
            end
            @maybe_vectorized Opt for j in jmin+1:jmax-1
                z = x[j-1,k,l] - x[j,k,l]
                y[j,k] = f(α, z, β, y[j,k])
            end
            let j = jmax
                z = x[j-1,k,l]
                y[j,k] = f(α, z, β, y[j,k])
            end
        end
    elseif jmin == jmax && β != 1
        let j = jmin, z = zero(T)
            @maybe_vectorized Opt for k in CartesianIndices(K)
                y[j,k] = f(α, z, β, y[j,k])
            end
        end
    end
    nothing
end
#
# Apply adjoint of 1st order finite differences along 2nd and subsequent
# dimensions:
#
function unsafe_apply!(f::Function,
                       α::Number,
                       ::Type{Adjoint},
                       A::Diff{1,Any,Opt},
                       x::AbstractArray,
                       β::Number,
                       y::AbstractArray,
                       I::ArrayAxes,
                       J::ArrayAxis,
                       K::ArrayAxes,
                       l::CartesianIndex) where {Opt}
    T = real(eltype(x))
    jmin, jmax = limits(J)
    if jmin < jmax
        @maybe_inbounds Opt for k in CartesianIndices(K)
            let j = jmin
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    z = -x[i,j,k,l]
                    y[i,j,k] = f(α, z, β, y[i,j,k])
                end
            end
            for j in jmin+1:jmax-1
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    z = x[i,j-1,k,l] - x[i,j,k,l]
                    y[i,j,k] = f(α, z, β, y[i,j,k])
                end
            end
            let j = jmax
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    z = x[i,j-1,k,l]
                    y[i,j,k] = f(α, z, β, y[i,j,k])
                end
            end
        end
    elseif jmin == jmax && β != 1
        let j = jmin, z = zero(T)
            @maybe_inbounds Opt for k in CartesianIndices(K)
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k] = f(α, z, β, y[i,j,k])
                end
            end
        end
    end
    nothing
end
#
# The Gram composition D'*D of the 1st order forward finite differences D with
# flat boundary conditions writes:
#
#     D'*D = [  1  -1   0   0   0
#              -1   2  -1   0   0
#               0  -1   2  -1   0
#               0   0  -1   2  -1
#               0   0   0  -1   1 ]
#
# Apply D'*D along 1st dimension:
#
function unsafe_apply!(f::Function,
                       α::Number,
                       ::Type{<:Union{Direct,Adjoint}},
                       A::Gram{Diff{1,Any,Opt}},
                       x::AbstractArray,
                       β::Number,
                       y::AbstractArray,
                       I::Tuple{},
                       J::ArrayAxis,
                       K::ArrayAxes) where {Opt}
    T = real(eltype(x))
    jmin, jmax = limits(J)
    if jmin < jmax
        @maybe_inbounds Opt for k in CartesianIndices(K)
            let j = jmin
                z = x[j,k] - x[j+1,k]
                y[j,k] = f(α, z, β, y[j,k])
            end
            @maybe_vectorized Opt for j in jmin+1:jmax-1
                z = T(2)*x[j,k] - (x[j-1,k] + x[j+1,k])
                y[j,k] = f(α, z, β, y[j,k])
            end
            let j = jmax
                z = x[j,k] - x[j-1,k]
                y[j,k] = f(α, z, β, y[j,k])
            end
        end
    elseif jmin == jmax && β != 1
        let j = jmin, z = zero(T)
            @maybe_vectorized Opt for k in CartesianIndices(K)
                y[j,k] = f(α, z, β, y[j,k])
            end
        end
    end
    nothing
end
#
# Apply  D'*D along 2nd and subsequent dimensions:
#
function unsafe_apply!(f::Function,
                       α::Number,
                       ::Type{<:Union{Direct,Adjoint}},
                       A::Gram{Diff{1,Any,Opt}},
                       x::AbstractArray,
                       β::Number,
                       y::AbstractArray,
                       I::ArrayAxes,
                       J::ArrayAxis,
                       K::ArrayAxes) where {Opt}
    T = real(eltype(x))
    jmin, jmax = limits(J)
    if jmin < jmax
        @maybe_inbounds Opt for k in CartesianIndices(K)
            let j = jmin
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    z = x[i,j,k] - x[i,j+1,k]
                    y[i,j,k] = f(α, z, β, y[i,j,k])
                end
            end
            for j in jmin+1:jmax-1
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    z = T(2)*x[i,j,k] - (x[i,j-1,k] + x[i,j+1,k])
                    y[i,j,k] = f(α, z, β, y[i,j,k])
                end
            end
            let j = jmax
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    z = x[i,j,k] - x[i,j-1,k]
                    y[i,j,k] = f(α, z, β, y[i,j,k])
                end
            end
        end
    elseif jmin == jmax && β != 1
        let j = jmin, z = zero(T)
            @maybe_inbounds Opt for k in CartesianIndices(K)
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k] = f(α, z, β, y[i,j,k])
                end
            end
        end
    end
    nothing
end

#------------------------------------------------------------------------------
#
# 2nd order finite differences with flat boundary conditions are computed by:
#
#     D = [-1   1   0   0   0   0
#           1  -2   1   0   0   0
#           0   1  -2   1   0   0
#           0   0   1  -2   1   0
#           0   0   0   1  -2   1
#           0   0   0   0   1  -1]
#
# Remarks:
#
#  - Applying this operator on a single dimension is self-adjoint.
#
#  - For a single dimension, this operator is the opposite of the Gram
#    composition of 1st order finite differences (backward or forward).
#
# Apply 2nd order finite differences along 1st dimension:
#
function unsafe_apply!(f::Function,
                       α::Number,
                       ::Type{Direct},
                       A::Diff{2,Any,Opt},
                       x::AbstractArray,
                       β::Number,
                       y::AbstractArray,
                       I::Tuple{},
                       J::ArrayAxis,
                       K::ArrayAxes,
                       l::CartesianIndex) where {Opt}
    T = real(eltype(x))
    jmin, jmax = limits(J)
    if jmin < jmax
        @maybe_inbounds Opt for k in CartesianIndices(K)
            let j = jmin
                z = x[j+1,k] - x[j,k]
                y[j,k,l] = f(α, z, β, y[j,k,l])
            end
            @maybe_vectorized Opt for j in jmin+1:jmax-1
                z = x[j-1,k] + x[j+1,k] - T(2)*x[j,k]
                y[j,k,l] = f(α, z, β, y[j,k,l])
            end
            let j = jmax
                z = x[j-1,k] - x[j,k]
                y[j,k,l] = f(α, z, β, y[j,k,l])
            end
        end
    elseif jmin == jmax && β != 1
        let j = jmin, z = zero(T)
            @maybe_vectorized Opt for k in CartesianIndices(K)
                y[j,k,l] = f(α, z, β, y[j,k,l])
            end
        end
    end
    nothing
end
#
# Apply 2nd order finite differences along 2nd and subsequent dimensions:
#
function unsafe_apply!(f::Function,
                       α::Number,
                       ::Type{Direct},
                       A::Diff{2,Any,Opt},
                       x::AbstractArray,
                       β::Number,
                       y::AbstractArray,
                       I::ArrayAxes,
                       J::ArrayAxis,
                       K::ArrayAxes,
                       l::CartesianIndex) where {Opt}
    T = real(eltype(x))
    jmin, jmax = limits(J)
    if jmin < jmax
        @maybe_inbounds Opt for k in CartesianIndices(K)
            let j = jmin
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    z = x[i,j+1,k] - x[i,j,k]
                    y[i,j,k,l] = f(α, z, β, y[i,j,k,l])
                end
            end
            for j in jmin+1:jmax-1
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    # Other possibility:
                    # z = (x[i,j-1,k] - x[i,j,k]) + (x[i,j+1,k] - x[i,j,k])
                    z = x[i,j-1,k] + x[i,j+1,k] - T(2)*x[i,j,k]
                    y[i,j,k,l] = f(α, z, β, y[i,j,k,l])
                end
            end
            let j = jmax
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    z = x[i,j-1,k] - x[i,j,k]
                    y[i,j,k,l] = f(α, z, β, y[i,j,k,l])
                end
            end
        end
    elseif jmin == jmax && β != 1
        let j = jmin, z = zero(T)
            @maybe_inbounds Opt for k in CartesianIndices(K)
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k,l] = f(α, z, β, y[i,j,k,l])
                end
            end
        end
    end
    nothing
end
#
# Apply adjoint of 2nd order finite differences along 1st dimension:
#
function unsafe_apply!(f::Function,
                       α::Number,
                       ::Type{Adjoint},
                       A::Diff{2,Any,Opt},
                       x::AbstractArray,
                       β::Number,
                       y::AbstractArray,
                       I::Tuple{},
                       J::ArrayAxis,
                       K::ArrayAxes,
                       l::CartesianIndex) where {Opt}
    T = real(eltype(x))
    jmin, jmax = limits(J)
    if jmin < jmax
        @maybe_inbounds Opt for k in CartesianIndices(K)
            let j = jmin
                z = x[j+1,k,l] - x[j,k,l]
                y[j,k] = f(α, z, β, y[j,k])
            end
            @maybe_vectorized Opt for j in jmin+1:jmax-1
                # Other possibility:
                # z = (x[j-1,k,l] - x[j,k,l]) + (x[j+1,k,l] - x[j,k,l])
                z = x[j-1,k,l] + x[j+1,k,l] - T(2)*x[j,k,l]
                y[j,k] = f(α, z, β, y[j,k])
            end
            let j = jmax
                z = x[j-1,k,l] - x[j,k,l]
                y[j,k] = f(α, z, β, y[j,k])
            end
        end
    elseif jmin == jmax && β != 1
        let j = jmin, z = zero(T)
            @maybe_vectorized Opt for k in CartesianIndices(K)
                y[j,k] = f(α, z, β, y[j,k])
            end
        end
    end
    nothing
end
#
# Apply 2nd order finite differences along 2nd and subsequent dimensions:
#
function unsafe_apply!(f::Function,
                       α::Number,
                       ::Type{Adjoint},
                       A::Diff{2,Any,Opt},
                       x::AbstractArray,
                       β::Number,
                       y::AbstractArray,
                       I::ArrayAxes,
                       J::ArrayAxis,
                       K::ArrayAxes,
                       l::CartesianIndex) where {Opt}
    T = real(eltype(x))
    jmin, jmax = limits(J)
    if jmin < jmax
        @maybe_inbounds Opt for k in CartesianIndices(K)
            let j = jmin
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    z = x[i,j+1,k,l] - x[i,j,k,l]
                    y[i,j,k] = f(α, z, β, y[i,j,k])
                end
            end
            for j in jmin+1:jmax-1
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    z = x[i,j-1,k,l] + x[i,j+1,k,l] - T(2)*x[i,j,k,l]
                    y[i,j,k] = f(α, z, β, y[i,j,k])
                end
            end
            let j = jmax
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    z = x[i,j-1,k,l] - x[i,j,k,l]
                    y[i,j,k] = f(α, z, β, y[i,j,k])
                end
            end
        end
    elseif jmin == jmax && β != 1
        let j = jmin, z = zero(T)
            @maybe_inbounds Opt for k in CartesianIndices(K)
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k] = f(α, z, β, y[i,j,k])
                end
            end
        end
    end
    nothing
end
#
# The Gram composition of 2nd order finite differences writes:
#
#    D'*D = [ 2  -3   1   0   0   0             (1)
#            -3   6  -4   1   0   0             (2)
#             1  -4   6  -4   1   0             (3)
#             0   1  -4   6  -4   1             (3)
#             0   0   1  -4   6  -3             (4)
#             0   0   0   1  -3   2]            (5)
#
# The above is for len ≥ 4, with len is the length of the dimension of
# interest, omitting the Eq. (5) for len = 4 and repeating Eq. (5) as necessary
# for the central rows for n ≥ 5.  For len = 3:
#
#    D'*D = [ 2  -3   1                         (1)
#            -3   6  -3                         (6)
#             1  -3   2]                        (5)
#
# For len = 2:
#
#    D'*D = [ 2  -2                             (7)
#            -2   2]                            (8)
#
# For len = 1, D = 0 and D'*D = 0 (the null 1×1 operator).
#
# Methods to apply the rows of D'D ():
#
# - Eq. (1), first row when len ≥ 3:
#
@inline @propagate_inbounds D2tD2_1(x::AbstractArray, j::Int, k) = begin
    T = real(eltype(x))
    T(2)*x[j,k] - T(3)*x[j+1,k] + x[j+2,k]
end
@inline @propagate_inbounds D2tD2_1(x::AbstractArray, i, j::Int, k) = begin
    T = real(eltype(x))
    T(2)*x[i,j,k] - T(3)*x[i,j+1,k] + x[i,j+2,k]
end
#
# - Eq. (2), second row when len ≥ 4:
#
@inline @propagate_inbounds D2tD2_2(x::AbstractArray, j::Int, k) = begin
    T = real(eltype(x))
    T(6)*x[j,k] - T(3)*x[j-1,k] - T(4)*x[j+1,k] + x[j+2,k]
end
@inline @propagate_inbounds D2tD2_2(x::AbstractArray, i, j::Int, k) = begin
    T = real(eltype(x))
    T(6)*x[i,j,k] - T(3)*x[i,j-1,k] - T(4)*x[i,j+1,k] + x[i,j+2,k]
end
#
# - Eq. (3), central rows when len ≥ 5:
#
@inline @propagate_inbounds D2tD2_3(x::AbstractArray, j::Int, k) = begin
    T = real(eltype(x))
    (x[j-2,k] + x[j+2,k]) + T(6)*x[j,k] - T(4)*(x[j-1,k] + x[j+1,k])
end
@inline @propagate_inbounds D2tD2_3(x::AbstractArray, i, j::Int, k) = begin
    T = real(eltype(x))
    (x[i,j-2,k] + x[i,j+2,k]) + T(6)*x[i,j,k] - T(4)*(x[i,j-1,k] + x[i,j+1,k])
end
#
# - Eq. (4), before last row when len ≥ 4:
#
@inline @propagate_inbounds D2tD2_4(x::AbstractArray, j::Int, k) = begin
    T = real(eltype(x))
    T(6)*x[j,k] - T(3)*x[j+1,k] - T(4)*x[j-1,k] + x[j-2,k]
end
@inline @propagate_inbounds D2tD2_4(x::AbstractArray, i, j::Int, k) = begin
    T = real(eltype(x))
    T(6)*x[i,j,k] - T(3)*x[i,j+1,k] - T(4)*x[i,j-1,k] + x[i,j-2,k]
end
#
# - Eq. (5), last row when len ≥ 3:
#
@inline @propagate_inbounds D2tD2_5(x::AbstractArray, j::Int, k) = begin
    T = real(eltype(x))
    T(2)*x[j,k] - T(3)*x[j-1,k] + x[j-2,k]
end
@inline @propagate_inbounds D2tD2_5(x::AbstractArray, i, j::Int, k) = begin
    T = real(eltype(x))
    T(2)*x[i,j,k] - T(3)*x[i,j-1,k] + x[i,j-2,k]
end
#
# - Eq. (6), central row when len = 3:
#
@inline @propagate_inbounds D2tD2_6(x::AbstractArray, j::Int, k) = begin
    T = real(eltype(x))
    T(6)*x[j,k] - T(3)*(x[j-1,k] + x[j+1,k])
end
@inline @propagate_inbounds D2tD2_6(x::AbstractArray, i, j::Int, k) = begin
    T = real(eltype(x))
    T(6)*x[i,j,k] - T(3)*(x[i,j-1,k] + x[i,j+1,k])
end
#
# - Eq. (7), first row when len = 2:
#
@inline @propagate_inbounds D2tD2_7(x::AbstractArray, j::Int, k) = begin
    z = x[j,k] - x[j+1,k]
    return z + z
end
@inline @propagate_inbounds D2tD2_7(x::AbstractArray, i, j::Int, k) = begin
    z = x[i,j,k] - x[i,j+1,k]
    return z + z
end
#
# - Eq. (8), last row when len = 2:
#
@inline @propagate_inbounds D2tD2_8(x::AbstractArray, j::Int, k) = begin
    z = x[j,k] - x[j-1,k]
    return z + z
end
@inline @propagate_inbounds D2tD2_8(x::AbstractArray, i, j::Int, k) = begin
    z = x[i,j,k] - x[i,j-1,k]
    return z + z
end
#
# Apply Gram composition of 2nd order finite differences along 1st dimension:
#
function unsafe_apply!(f::Function,
                       α::Number,
                       ::Type{<:Union{Direct,Adjoint}},
                       A::Gram{Diff{2,Any,Opt}},
                       x::AbstractArray,
                       β::Number,
                       y::AbstractArray,
                       I::Tuple{},
                       J::ArrayAxis,
                       K::ArrayAxes) where {Opt}
    T = real(eltype(x))
    jmin, jmax = limits(J)
    len = length(J)
    if len ≥ 5
        @maybe_inbounds Opt for k in CartesianIndices(K)
            let j = jmin
                y[j,k] = f(α, D2tD2_1(x,j,k), β, y[j,k])
            end
            let j = jmin+1
                y[j,k] = f(α, D2tD2_2(x,j,k), β, y[j,k])
            end
            @maybe_vectorized Opt for j in jmin+2:jmax-2
                y[j,k] = f(α, D2tD2_3(x,j,k), β, y[j,k])
            end
            let j = jmax-1
                y[j,k] = f(α, D2tD2_4(x,j,k), β, y[j,k])
            end
            let j = jmax
                y[j,k] = f(α, D2tD2_5(x,j,k), β, y[j,k])
            end
        end
    elseif len == 4
        @maybe_vectorized Opt for k in CartesianIndices(K)
            let j = jmin
                y[j,k] = f(α, D2tD2_1(x,j,k), β, y[j,k])
            end
            let j = jmin+1
                y[j,k] = f(α, D2tD2_2(x,j,k), β, y[j,k])
            end
            let j = jmax-1
                y[j,k] = f(α, D2tD2_4(x,j,k), β, y[j,k])
            end
            let j = jmax
                y[j,k] = f(α, D2tD2_5(x,j,k), β, y[j,k])
            end
        end
    elseif len == 3
        @maybe_vectorized Opt for k in CartesianIndices(K)
            let j = jmin
                y[j,k] = f(α, D2tD2_1(x,j,k), β, y[j,k])
            end
            let j = jmin+1
                y[j,k] = f(α, D2tD2_6(x,j,k), β, y[j,k])
            end
            let j = jmax
                y[j,k] = f(α, D2tD2_5(x,j,k), β, y[j,k])
            end
        end
    elseif len == 2
        @maybe_vectorized Opt for k in CartesianIndices(K)
            let j = jmin
                y[j,k] = f(α, D2tD2_7(x,j,k), β, y[j,k])
            end
            let j = jmax
                y[j,k] = f(α, D2tD2_8(x,j,k), β, y[j,k])
            end
        end
    elseif len == 1 && β != 1
        let j = jmin, z = zero(T)
            @maybe_vectorized Opt for k in CartesianIndices(K)
                y[j,k] = f(α, z, β, y[j,k])
            end
        end
    end
    nothing
end
#
# Apply Gram composition of 2nd order finite differences along 2nd and
# subsequent dimensions:
#
function unsafe_apply!(f::Function,
                       α::Number,
                       ::Type{<:Union{Direct,Adjoint}},
                       A::Gram{Diff{2,Any,Opt}},
                       x::AbstractArray,
                       β::Number,
                       y::AbstractArray,
                       I::ArrayAxes,
                       J::ArrayAxis,
                       K::ArrayAxes) where {Opt}
    T = real(eltype(x))
    jmin, jmax = limits(J)
    len = length(J)
    if len ≥ 5
        @maybe_inbounds Opt for k in CartesianIndices(K)
            let j = jmin
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k] = f(α, D2tD2_1(x,i,j,k), β, y[i,j,k])
                end
            end
            let j = jmin+1
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k] = f(α, D2tD2_2(x,i,j,k), β, y[i,j,k])
                end
            end
            for j in jmin+2:jmax-2
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k] = f(α, D2tD2_3(x,i,j,k), β, y[i,j,k])
                end
            end
            let j = jmax-1
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k] = f(α, D2tD2_4(x,i,j,k), β, y[i,j,k])
                end
            end
            let j = jmax
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k] = f(α, D2tD2_5(x,i,j,k), β, y[i,j,k])
                end
            end
        end
    elseif len == 4
        @maybe_inbounds Opt for k in CartesianIndices(K)
            let j = jmin
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k] = f(α, D2tD2_1(x,i,j,k), β, y[i,j,k])
                end
            end
            let j = jmin+1
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k] = f(α, D2tD2_2(x,i,j,k), β, y[i,j,k])
                end
            end
            let j = jmax-1
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k] = f(α, D2tD2_4(x,i,j,k), β, y[i,j,k])
                end
            end
            let j = jmax
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k] = f(α, D2tD2_5(x,i,j,k), β, y[i,j,k])
                end
            end
        end
    elseif len == 3
        @maybe_inbounds Opt for k in CartesianIndices(K)
            let j = jmin
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k] = f(α, D2tD2_1(x,i,j,k), β, y[i,j,k])
                end
            end
            let j = jmin+1
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k] = f(α, D2tD2_6(x,i,j,k), β, y[i,j,k])
                end
            end
            let j = jmax
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k] = f(α, D2tD2_5(x,i,j,k), β, y[i,j,k])
                end
            end
        end
    elseif len == 2
        @maybe_inbounds Opt for k in CartesianIndices(K)
            let j = jmin
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k] = f(α, D2tD2_7(x,i,j,k), β, y[i,j,k])
                end
            end
            let j = jmax
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k] = f(α, D2tD2_8(x,i,j,k), β, y[i,j,k])
                end
            end
        end
    elseif len == 1 && β != 1
        let j = jmin, z = zero(T)
            @maybe_inbounds Opt for k in CartesianIndices(K)
                @maybe_vectorized Opt for i in CartesianIndices(I)
                    y[i,j,k] = f(α, z, β, y[i,j,k])
                end
            end
        end
    end
    nothing
end

end # module


text:
using LinearAlgebra
using SDPA

using Test

include("example1.jl")
include("example2.jl")

include("MOI_wrapper.jl")
include("MPB_wrapper.jl")


text:
# Original Java version Copyright 2005 Google Inc.
# Julia translation Copyright 2022 Ted Dunning
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using S2
using LinearAlgebra
using Test

# Basics

# Check that the conversion between Pi radians and 180 degrees is exact.
@test Float64(π) == radians(S1Angle(π))
@test 180.0 == degrees(S1Angle(π))
@test Float64(π) == radians(degrees(180))
@test 180.0 == degrees(degrees(180))

@test 90.0 == degrees(radians(π / 2))
@test 90.0 == degrees(S1Angle(π / 2))

@test -90.0 == degrees(S1Angle(-π / 2))
@test Float64(-π / 4) == radians(degrees(-45))

@test 30.0 ≈ degrees(degrees(10) + degrees(20))
@test -10 ≈ degrees(degrees(10) - degrees(20))
@test 20 ≈ degrees(2 * degrees(10))
@test 5 ≈ degrees(degrees(10) / 2.0)
@test 1 ≈ cos(S1Angle(0))
@test 1 ≈ sin(degrees(90))
@test 1 ≈ tan(degrees(45))

@test 100.0 * π ≈ distance(S1Angle(π), 100.0)
@test 50.0 * π ≈ distance(S1Angle(π / 2), 100.0)
@test 25.0 * π ≈ distance(radians(π / 4), 100.0)

@test 0.0 ≈ radians(normalizeRange(degrees(360.0)))
@test -90.0 ≈ degrees(normalizeRange(degrees(-90.0)))
@test 180.0 ≈ degrees(normalizeRange(degrees(-180.0)))
@test 90.0 ≈ degrees(normalizeRange(degrees(90.0)))
@test 180.0 ≈ degrees(normalizeRange(degrees(180.0)))
@test -90.0 ≈ degrees(normalizeRange(degrees(270.0)))
@test 180.0 ≈ degrees(normalizeRange(degrees(540.0)))
@test 90.0 ≈ degrees(normalizeRange(degrees(-270.0)))

@test π ≈ radians(normalizeRange(radians(π)))

# -π maps to π.
@test π ≈ radians(normalizeRange(radians(-π)))

let a = degrees(90.0)
    @test a === normalizeRange(a)
end

@test 1.5 == radians(S1Angle(0.5) + S1Angle(0.75) + S1Angle(0.25))

@test 180 == degrees(normalizeRange(degrees(90) + degrees(45) + degrees(30) + degrees(15)))

@test 0.75 * π ≈ radians(S1Angle(π / 2) + degrees(45))


text:
using AmbientForcing
using PowerDynamics, Distributions
using OrdinaryDiffEq
using Test

@testset "PD_example" begin
    line_list = []
    append!(line_list, [StaticLine(from = 1, to = 2, Y = -1im / 0.02 + 4)])
    append!(line_list, [StaticLine(from = 1, to = 3, Y = -1im / 0.02 + 4)])

    node_list = []
    append!(node_list, [SlackAlgebraic(U = 1, Y_n = 0)])
    append!(node_list, [FourthOrderEq(H = 3.318, P = -0.6337, D = 0.1, Ω = 50, E_f = 0.5, T_d_dash = 0.1,  T_q_dash = 8.690,X_d_dash = 0.111,  X_q_dash = 0.103, X_d = 0.1, X_q = 0.6)])
    append!(node_list, [PQAlgebraic(P = -0.6337, Q = 0.0)]) # This is our constraint

    pg = PowerGrid(node_list, line_list)
    rpg = rhs(pg)

    g = constraint_equations(rpg)
    op = find_operationpoint(pg)
    @test isapprox(sum(g(op.vec)), 0.0, atol = 1e-8)

    Frand = random_force(rpg, [0.0,1], Uniform)
    z_new_all = ambient_forcing(rpg, op.vec, 2.0, Frand) # Our new valid inital condition
    @test isapprox(sum(g(z_new_all)), 0.0, atol = 1e-8)

    τ = 2.0 # the integration time
    idx = idx_exclusive(rpg, ["u_r_3", "u_i_3"])
    # genrate a vector only with non-vanishing componats the voltage at node 2
    Frand = random_force(rpg, [0.0, 1.0], Uniform, idx)
    z_new_node_3 = ambient_forcing(rpg, op.vec, τ, Frand)
    @test isapprox(sum(g(z_new_node_3)), 0.0, atol = 1e-8)

    idx = idx_exclusive(rpg, ["u_r_2", "u_i_2", "θ_2", "ω_2"])
    dist_vec = [[0,1], [0, 2], [0, 2π] ./τ, [-5, 5] ./τ]
    Frand = Frand = random_force(rpg, dist_vec, Uniform, idx)
    z_new_node_2 = ambient_forcing(rpg, op.vec, τ, Frand)
    @test isapprox(sum(g(z_new_node_2)), 0.0, atol = 1e-8)
end

@testset "OrdinaryDiffEq_example" begin
    # adopting the Robertson Example from the DifferentialEquation docs
    # https://diffeq.sciml.ai/stable/tutorials/advanced_ode_example/#Handling-Mass-Matrices
    function rober(du,u,p,t)
        y₁,y₂,y₃ = u
        du[1] = -0.04 * y₁ + 1e4 * y₂ * y₃
        du[2] =  0.04 * y₁ - 1e4 * y₂ * y₃ - 3e7 * y₂^2
        du[3] =  y₁ + y₂ + y₃ - 1
    nothing
  end

    M = [1. 0  0
         0  1. 0
        0  0  0];

    ode_rober = ODEFunction(rober, mass_matrix = M)
    u0 = [1.0,0.0,0.0]; # an inital condition which fulfills the constraint

    g_rober = constraint_equations(ode_rober)
    @test isapprox(sum(g_rober(u0)), 0.0, atol = 1e-8)

    # Randomly perturbing all varibales
    Frand = random_force(ode_rober, [0.0,1], Uniform)
    z_new_all = ambient_forcing(ode_rober, u0, 2.0, Frand)
    @test isapprox(sum(g_rober(z_new_all)), 0.0, atol = 1e-8)

    # only perturbing the second variable y₂
    h = [0, 1 ,0]
    z_new_2 = ambient_forcing(ode_rober, u0, 2.0, h)
    @test isapprox(sum(g_rober(z_new_2)), 0.0, atol = 1e-8)
end


text:
using Libdl, Printf

push!(Libdl.DL_LOAD_PATH, "/tmpbig/thofmann/juliagap/libgap/.libs/")
include("/tmpbig/thofmann/juliagap/LibGAP.jl/src/initialization.jl")
libgap.initialize( [ ""
                   , "-l", "/tmpbig/thofmann/juliagap/libgap/"
                   , "-T", "-r", "-A", "-q"
                   , "-m", "512m" ], ["\0"])

using Hecke

include("/users/cip/users/thofmann/CarloLibGAPTest/towards_cohomology.jl")
include("/tmpbig/thofmann/v0.6/Hecke/examples/FieldEnumeration.jl")
include("/tmpbig/thofmann/v0.6/Hecke/examples/fields.jl")

sprint_formatted(fmt, args...) = @eval @sprintf($fmt, $(args...))

grp_no = Int(eval(Meta.parse(ARGS[1])))

n, i = small_solvable_groups[grp_no]

file = sprint_formatted("%0$(length(string(length(small_solvable_groups))))d", grp_no) * "-$n-$i."* "log"

println("========================================")
println("Group: $n $i")
println("========================================")

global dbound = 0

if length(ARGS) == 2
  dbound = fmpz(eval(Meta.parse(ARGS[2])))
else
  if mod(n, 2) == 0
    dbound = Hecke.lower_discriminant_bound(n, div(n, 2))
  else
    dbound = Hecke.lower_discriminant_bound(n, 0)
  end
end

dbound = 10 * dbound

println("========================================")
println("Discriminant bound: $dbound")
println("========================================")

set_verbose_level(:FieldsNonFancy, 1)

l = []

global complex_fields = []
global real_fields = []
global all = []

global found_complex = false
global found_real = false

if mod(n, 2) != 0
  found_complex = true
end

global only_real = false

while true
  global complex_fields
  global real_fields
  global found_complex
  global found_real
  global only_real
  global dbound

  @show only_real

  global l = fields(n, i, dbound, only_real = only_real)

  for v in l
    x = v[1]
    if !found_complex
      if signature(x)[2] != 0
        push!(complex_fields, v)
      end
    end

    if !found_real
      if signature(x)[2] == 0
        push!(real_fields, v)
      end
    end
  end

  if !found_complex && length(complex_fields) > 0
    append!(all, complex_fields)
    found_complex = true
    println("== FOUND SOMETHING COMPLEX ==")
    empty!(complex_fields)
  end

  if !found_real && length(real_fields) > 0
    append!(all, real_fields)
    found_real = true
    println("== FOUND SOMETHING REAL ==")
    empty!(real_fields)
  end

  if found_complex
    if !found_real
      only_real = true
      println("== I AM IN ONLY_REAL = TRUE MODE ==")
    else
      println("== EXITING ==")
      break
    end
  end

  dbound = 100 * dbound
  println("========================================")
  println("Increasing discriminant bound: $dbound")
  println("========================================")
  println("only_real = $only_real")
  println("found_complex = $found_complex")
  println("found_real = $found_real")
end

l = all

ll = map(v -> v[1], l)
ffields = [ (x, discriminant(maximal_order(x))) for x in ll ]
sort!(ffields, lt = (x, y) -> abs(x[2]) <= abs(y[2]))

@show length(ffields)

_write_fields(ffields, file)


text:
reload("OBC")
using HDF5, JLD
using PyPlot
using Distributions

kmax = 1
D = 2
iters = 10000

loc = "/home/bana/largeresearch/rnaseq.jld"
if !isfile(loc)
    loc = "/home/bana/largeresearch/seq-data/rnaseq.jld"
end
@load loc brca_norm brca_norm_sampids brca_norm_geneids
@load loc luad_norm luad_norm_sampids luad_norm_geneids

blas_set_num_threads(1)

#pick two random genes
#data0 = brca_norm[[5,4,3],[1:40]]'
#data1 = luad_norm[[5,4,3],[1:40]]'
data0 = brca_norm[[5,3],[1:40]]'
data1 = luad_norm[[5,3],[1:40]]'
#data0 = convert(Matrix{Int}, round(brca_norm[:,[2,3]]))
#data1 = convert(Matrix{Int}, round(luad_norm[:,[2,3]]))

#println(brca_norm_geneids[[2,3]])
#println(luad_norm_geneids[[2,3]])

#plot(data0[:,1], data0[:,2], "go", label="brca")
#plot(data1[:,1], data1[:,2], "ro", label="luad")

#assert(false)

d = 100.0
######################################################################
# Class 0
######################################################################
cov = eye(D,D)
#cov = [0.5 -0.2; -0.2 1.0]
mu = ones(D,kmax)
#mu = [1.6 -1.0; 1.6 -1.0]
#mu = [0.2 0.2]'
prior = MPM.MPMPrior(D=D, kmax=kmax, kappa=10., S=eye(2)*10)
start = MPM.MPMParams(mu, #mu :: Matrix{Float64}
    cov, #sigma :: Matrix{Float64}
    ones(kmax)/kmax, #w :: Vector{Float64}
    clamp(log(data0'/d),-8.0,Inf), #lam :: Matrix{Float64}
    1) #k :: Int

pmoves = MPM.MPMPropMoves()
pmoves.lammove = 0.005
pmoves.mumove = 0.6
pmoves.priorkappa = 100.0

obj_a = MPM.MPMSampler(prior, data0, deepcopy(start), pmoves, d)
obj_a.usepriors = true

mymh_a = MPM.MHRecord(obj_a,burn=1000,thin=50)
sample(mymh_a,iters)

######################################################################
# Class 1
######################################################################

start = MPM.MPMParams(mu, #mu :: Matrix{Float64}
    cov, #sigma :: Matrix{Float64}
    ones(kmax)/kmax, #w :: Vector{Float64}
    clamp(log(data1'/d),-8.0,Inf), #lam :: Matrix{Float64}
    1) #k :: Int
obj_b = MPM.MPMSampler(prior, data1, deepcopy(start), pmoves, d)
mymh_b = MPM.MHRecord(obj_b,burn=1000,thin=50)
sample(mymh_b,iters)

######################################################################
# BEE Computation
######################################################################

#@time be1 = MPM.bee_e_data(vcat(data0,data1), mymh_a.db, mymh_b.db, 20)
#@time be2 = MPM.bee_e_data(vcat(data1,data1), mymh_a.db, mymh_b.db, 20)
@time bec = MPM.bee_e_cube(vcat(data1,data1), mymh_a.db, mymh_b.db, 20)

#@show be1
#@show be2
@show bec

assert(false)
######################################################################
# Plotting
######################################################################

include("src/plot_utils.jl")
mins, maxs, (lens, steps, grid) = MPM.get_grid(vcat(data0, data1), factor=2, N=100)
#n1, n2, gext, grid = gen_grid([0,10.0,0,10])

#gavg,lams = MPM.calc_g(grid, [start], 1000)
#imshow(reshape(gavg,N,N), extent=gext, aspect="equal", origin="lower")
#colorbar()
#plot(data0[:,1], data0[:,2], "g.", alpha=0.8)
#plot(10*exp(lams[1,:]), 10*exp(lams[2,:]), "k.")

#close("all")

ga = MPM.calc_g(grid, mymh_a.db, 20, dmean=d)
#imshow(reshape(ga,n1,n2)', extent=gext, aspect="equal", origin="lower")
#colorbar()
#contour(reshape(ga,n1,n2)', extent=gext, aspect="equal", origin="lower")
#plot(data0[:,1], data0[:,2], "g.", alpha=0.8)

#figure()
gb = MPM.calc_g(grid, mymh_b.db, 20, dmean=d)
#imshow(reshape(gb,n1,n2)', extent=gext, aspect="equal", origin="lower")
#colorbar()
#contour(reshape(gb,n1,n2)', extent=gext, aspect="equal", origin="lower")
#plot(data1[:,1], data1[:,2], "m.", alpha=0.8)

figure()
gavg = exp(min(ga,gb))
#gavg = ga - gb
#gavg = ga
gext = [mins[1], maxs[1], mins[2], maxs[2]]
imshow(reshape(gavg,lens[1],lens[2])', extent=gext, aspect="equal", origin="lower", interpolation="nearest")
colorbar()
#contour(reshape(gavg,lens[1],lens[2])', (0.0,), extent=gext, aspect="equal", origin="lower")
plot(data0[:,1], data0[:,2], "g.", alpha=0.8)
plot(data1[:,1], data1[:,2], "r.", alpha=0.8)

#println("$(MPM.error_points(mymh_a.db, mymh_b.db, [tst_data; tst_data1], [zeros(size(tst_data,1)), ones(size(tst_data1,2))] ))")

#g() = plot(tst_data[:,1]+rand(size(tst_data,1)), tst_data[:,2]+rand(size(tst_data,1)), "g.", alpha=0.3)
#fa() = plot_traces(mymh_a.db, [:mu,:sigma,:lam,:energy])
#fb() = plot_traces(mymh_b.db, [:mu,:sigma,:lam,:energy])