text:
##############################################################################################################################################
# Implements the Elo algorithm as published by Arpad Elo in the 1960s.
# Two differences from the standard implementations. first, the scores are done with natural powers rather than the common Elo style of base
# ten and a scaling of 400 points. To get back, you use.
#
# Natural = Elo_style * ln(10) / 400
#
# In addition, Elo are traditionally started at 1000, where in this case, we set zero as the default rating (so negative scores are possible.)
#
# The second difference is the implementation of time periods, which matches Elo's original approach - re-rating each time period or tournament
# rather than treating each game separately. You can recover the more common game-by-game approach by simply omitting the time period column.
# This algorithm is also much faster (fast_fit!, below)
##############################################################################################################################################

#Struct is an Elo object - in this version it's a dict and a base kfactor
struct Elo
    kfac::Real
    ratings::Dict
    default_rating::Real
end

function Elo(;kfac = 0.31, ratings = Dict{Int64, Float64}(), default_rating = 0.0)
    #Set up a new Elo struct, either using key words or defaults if not provided
    return Elo(kfac, ratings, default_rating)
end

function fit!(m::Elo, games::DataFrame; verbose = false)
    #= Wrapper function to choose the right approach to fit a dataframe of games =#
    titles = names(games)

    if !issubset(["P1", "P2", "P1_wins", "P2_wins"], titles)
        error("For rating, we need at least four columns named P1, P2, P1_wins, P2_wins")
    end

    #Depending on whether time periods are specified or not
    if !in("Period", names(games))
        verbose && println("No distinct periods found. Treating each game as a separate rating period and using fastelo algorithm.")
        fast_fit!.(Ref(m), games.P1, games.P2, games.P1_wins, games.P2_wins)
    else
        period_fit!(m, games[!, [:P1, :P2, :P1_wins, :P2_wins, :Period]])
    end
    return nothing
end

function fast_fit!(m::Elo, P1, P2, P1_wins, P2_wins)
    #=
    This function is used when matches themsleves define time periods. It
    simply updates the two players' elo ratings once per match
    Do not use when one is doing tournament or by-day elo updates where one might
    have many-one/one-many matches.
    =#
    surprise = m.kfac * (P1_wins - predict(m, P1, P2) * (P1_wins + P2_wins))
    m.ratings[P1] = get(m.ratings, P1, m.default_rating) + surprise
    m.ratings[P2] = get(m.ratings, P2, m.default_rating) - surprise
end

function period_fit!(m::Elo, original_games::DataFrame)
    #=
    Proper fit function, using defined periods to re-rate in. Takes only
    columns of a dataframe named (P1, P2, P1_wins, P2_wins, Day)
    =#
    #Duplicate whole rating dataframe with reversed results
    games = dupe_for_rating(original_games)

    #Split into days
    for day_games in groupby(games, :Period, sort = true)
        #Add elo predictions to each game in each day
        day_games = DataFrame(day_games)
        day_games[!, :Predict] = predict.(Ref(m), day_games.P1, day_games.P2) .* (day_games.P1_wins .+ day_games.P2_wins)
        #Groupby Player1 including results
        aggregated_by_P1_games = combine(groupby(day_games, :P1), :P1_wins => sum => :P1_wins, :P2_wins => sum => :P2_wins, :Predict => sum => :Predict)
        update!(m, aggregated_by_P1_games)
    end
end

function predict(m::Elo, i::Integer, j::Integer; kwargs...)
    #Returns the predicted result for any two players in the existing Elo rating dictionary, for a single game (%)
    return 1.0 / (1.0 + exp((get(m.ratings, j, m.default_rating) - get(m.ratings, i, m.default_rating))))
end

function update!(m::Elo, games)
    #Update function allowing for period in which there may be duplicate games
    for row in eachrow(games)
        if ismissing(row.P1_wins)
            return nothing
        end
        m.ratings[row.P1] = get(m.ratings, row.P1, m.default_rating) + m.kfac * (row.P1_wins - row.Predict)
    end
end

function display_ratings(elo::Elo; rating_type = "natural")
    #Show all ratings for players in DataFrame
    a = collect(keys(elo.ratings))
    b = map.(x -> elo.ratings[x], a)

    if rating_type == "elo"
        b = convert_natural_to_elo.(b)
    end

    players = DataFrame(player_code = a, player_rating = b)
    return sort!(players, :player_rating, rev = true)
end

function display_ratings(elo::Elo, players; rating_type = "natural")
    #Show ratings in Dataframe for selected players
    player_dict = Dict(players[:, 1] .=> players[:, 2])
    ratings = display_ratings(elo, rating_type = rating_type)
    ratings.player_names = map.(x -> player_dict[x], ratings.player_code)
    return ratings
end

function convert_natural_to_elo(x)
    #Convert natural ratings (used internally) to the normal Elo scale
    return x * 400 / log(10) + 1000
end


text:
using DenoiseDNA
using Documenter

DocMeta.setdocmeta!(DenoiseDNA, :DocTestSetup, :(using DenoiseDNA); recursive=true)

makedocs(;
    modules=[DenoiseDNA],
    authors="Arthur Newbury",
    repo="https://github.com/EvoArt/DenoiseDNA.jl/blob/{commit}{path}#{line}",
    sitename="DenoiseDNA.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://EvoArt.github.io/DenoiseDNA.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/EvoArt/DenoiseDNA.jl",
)

ext: jl


text:
# TODO: exports
export {{ (parser.unchecked_functions.keys() + parser.checked_functions.keys())|sort|join(', ') }}

# Macro for calling SCIP functions that return misc. types
macro scip_ccall(func, args...)
    return quote
        ccall(($func, "libscipopt"), $(args...))
    end
end

# Macro for calling SCIP functions that have checked return codes
macro scip_ccall_check(func, args...)
    return quote
        ret = ccall(($func, "libscipopt"), SCIP_RETCODE, $(args...))
        if ret != SCIP_OKAY
            error("SCIP Error: $ret")
        end
    end
end

# SCIP function wrappers: unchecked functions
{% for func_name, (ret_type, arg_types, arg_names, arg_vals) in parser.unchecked_functions.items() %}{{ func_name }}({{ arg_names|join(', ') }}) = @scip_ccall("{{ func_name }}", {% if ret_type %}{{ ret_type }}{% else %}Void{% endif %}, ({{ arg_types|join(', ') }}{% if arg_types|length == 1 %},{% endif %}){% if arg_vals %}, {{ arg_vals|join(', ') }}{% endif %})
{% endfor %}
# SCIP function wrappers: unchecked functions
{% for func_name, (arg_types, arg_names, arg_vals) in parser.checked_functions.items() %}{{ func_name }}({{ arg_names|join(', ') }}) = @scip_ccall_check("{{ func_name }}", ({{ arg_types|join(', ') }}{% if arg_types|length == 1 %},{% endif %}){% if arg_vals %}, {{ arg_vals|join(', ') }}{% endif %})
{% endfor %}


text:
"""
   invariantize_embedding(
      embedding::Array{Float64, 2};
      max_point_remove::Int = ceil(Int, size(embedding, 1)*0.05)
      )

If `remove_points = true`, iteratively remove the last point of an embedding
until it is invariant. If it is not possible to render the embedding invariant,
return an empty array. The default is to try to remove a maximum of ~5% of the
points of the original embedding before giving up. The number of points we're
allowed to remove can be set by providing the named argument `max_point_remove`.

If `remove_points = false`, incrementally move last point of the embedding
towards the origin until it lies within the convex hull of all preceding points.
"""
function invariantize_embedding(
      embedding::Array{Float64, 2};
      remove_points = false,
      max_point_remove::Int = ceil(Int, size(embedding, 1)*0.05)
      )

   if size(unique(embedding, 1)) < size(embedding)
      warn("Embedding points are not unique. Returning nothing.")
      return nothing
   end

   #=
   # Keep track of the embedding's centerpoint and the original position of the
   # last point in the embedding, so we can move the last point along a line
   # from its orinal position towards the embedding's center, until the point
   # lies inside the convex hull of the preceding points.
   =#
   embedding_center = mean(embedding, 1)
   lastpoint = embedding[end, :]
   direction = embedding_center.' - lastpoint

   pts_removed = 0
   is_invariant = false
   percent_moved = 0.0

   while !is_invariant && pts_removed <= max_point_remove
      pts_removed += 1

      if invariantset(embedding[1:(size(embedding, 1) - pts_removed), :])
         is_invariant = true
      else
         percent_moved += 1
         warn("Moved point $percent_moved % towards origin to fit inside
               convex hull of previous points")
         embedding[end, :] = lastpoint + direction*(percent_moved / 100)
      end
   end

   if is_invariant
      return embedding[1:(size(embedding, 1) - pts_removed), :]
   else
      warn("Could not make embedding invariant. Returning nothing.")
      return nothing
   end
end

max_stars_repo_name: kahaaga/InvariantDistribution

max_issues_repo_name: kahaaga/InvariantDistribution

max_forks_repo_name: kahaaga/InvariantDistribution


text:
using DataStructures, StringBuilders, Missings

const keyfilename = "keyfile.csv"

const OD = OrderedDict

categories = OD("Liability" =>
                ["Transfer", "Other"],
                "Expense" =>
                [OD("Food" =>
                    ["Cafe", "Groceries", "Restaurant"]),
                 OD("Alcohol" =>
                    ["Shop", "Bar"]),
                 OD("Health" =>
                    ["Sport", "Other"]),
                 OD("Transport" =>
                    ["Private", "Public"]),
                 OD("FixedCosts" =>
                    ["Rent", "Internet", "Electricity", "Phone", "Daycare", "CSN", "Insurance", "BankCard"]),
                 "Shopping", "Fun", "Travel", "Other"],
                "Revenue" =>
                ["ParentalSupport", "Salary", "ExternalTransfer", "Other"])

cats = String[]
sb = StringBuilder()
n = 0
for (k, vs) in categories
    append!(sb, k, "\n")
    for v in vs
        if v isa OD
            for (k2, v2s) in v
                append!(sb, "\t", k2, "\n")
                for v2 in v2s
                    n += 1
                    append!(sb, "\t\t$n) ", v2, "\n")
                    push!(cats, "$k:$k2:$v2")
                end
            end
        else
            n += 1
            append!(sb, "\t$n) ", v, "\n")
            push!(cats, "$k:$v")
        end
    end
end
msg = String(sb)
print(msg)

function ask(comment, amount, person, date)
    print_with_color(:blue, comment, ", ", amount, " SEK, ", person, ", ", date, '\n')
    run(`googler -c se --np -n 5 $comment`)
    txt = readline()
    if all(isnumber, txt)
        i = parse(Int, txt)
        if 1 ≤ i ≤ n
            return cats[i]
        end
    end
    print_with_color(:red, "must be an integer between 1 and $n\n")
    ask(comment, amount, person, date)
end

function getkey(keyfile)
    key = Dict{String, String}()
    open(keyfile, "r") do o
        for l in eachline(o)
            comment, category = split(l, ",")
            @assert !haskey(key, comment) "found duplicate comment: $comment"
            key[comment] = category
        end
    end
    key
end

function addcomment!(o, key, comment, amount, person, date)
    category = ask(comment, amount, person, date)
    println(o, comment, ",", category)
    key[comment] = category
end

function categorise!(table, folder)
    keyfile = joinpath(folder, keyfilename)
    key = getkey(keyfile)
    table[:category] = ""
    open(keyfile, "a") do o
        for row in DataFrames.eachrow(table)
            comment = row[:description]
            row[:category] = haskey(key, comment) ? key[comment] : addcomment!(o, key, comment, row[:amount], row[:person], row[:date])
        end
    end
    table
end


text:
import Reinforce: action, actions, finished, ismdp, maxsteps, reset!, reward, state, step!

function test_ep_iteration()
  @info "interface::iteration::maxsteps"
  env = FooEnv()
  π = FooPolicy()
  ep = Episode(env, π)

  # given the default of `finished` is `false`, iteration should hit `maxsteps`
  for (i, (s, a, r, s′)) ∈ enumerate(ep)
    @test s == i
    @test r == -1
    @test a ∈ [1, 2, 3]
    @test s′ == i + 1
  end

  @test ep.niter  == 3
  @test ep.total_reward == -3
end  # function test_ep_iteration

function test_ep_finished()
  @info "interface::iteration::finished"
  env = FooEnv()
  π = FooPolicy()
  ep = Episode(env, π)

  for (s, a, r, s′) ∈ ep
    nothing
  end

  # @eval finished(::FooEnv, s′) = false

  @test ep.niter == 1
  @test ep.total_reward == -1
end  # function test_ep_iteration

function test_run_episode()
  @info "interface::run_episode"

  env = FooEnv()
  π = FooPolicy()

  run_episode(env, π) do sars
    s, a, r, s′ = sars

    @test a ∈ [1, 2, 3]
    @test r == -1
    @test s + 1 == s′
  end
end  # function test_run_episode

@testset "interface" begin
  @info "interface::iteration"
  test_ep_iteration()
  begin
    @eval finished(::FooEnv, s′) = (s′ == 2)
    test_ep_finished()
    @eval finished(::FooEnv, s′) = false  # reset to default
  end

  test_run_episode()
end  # @testset "env interface"


text:
abstract type AnatomicalImage <: Imaging end

"""
"""
struct T1Weighted <: AnatomicalImage end
modalitylabel(::T1Weighted) = "T1w"

struct T2Weighted <: AnatomicalImage end
modalitylabel(::T2Weighted) = "T2w"

struct T1RhoMap <: AnatomicalImage end
modalitylabel(::T1RhoMap) = "T1rho"

struct T1Map <: AnatomicalImage end
modalitylabel(::T1Map) = "T1map"

struct T2Map <: AnatomicalImage end
modalitylabel(::T2Map) = "T2map"

struct T2Star <: AnatomicalImage end
modalitylabel(::T2Star) = "T2star"

"""
    Fluid-Attenuated Inversion Recovery (FLAIR)
"""
struct FluidAttenuatedInversionRecovery <: AnatomicalImage end
const FLAIR = FluidAttenuatedInversionRecovery
modalitylabel(::FLAIR) = "FLAIR"

"""
    Fast Imaging with Steady-state free Precision
"""
struct FastImagingwithSteadyStateFreePrecission <: FunctionalMRI end
const FISPImage = FastImagingwithSteadyStateFreePrecission

struct ProtonDensity <: AnatomicalImage end
modalitylabel(::ProtonDensity) = "PD"

struct ProtonDensityMap <: AnatomicalImage end
modalitylabel(::ProtonDensityMap) = "PDmap"

"""
    Magnetization Prepared Rapid Acquisition With Gradient Echo (MPRAGE)
"""
struct MagnetizationPreparedRapidAcquisitionWithGradientEcho <: AnatomicalImage end
const MPRAGEImage = MagnetizationPreparedRapidAcquisitionWithGradientEcho

"""
    Fast Spoiled Gradient Echo (FSPGR)
"""
struct FastSpoiledGradientEcho <: AnatomicalImage end
const FSPGRImage = FastSpoiledGradientEcho

"""
    RapidAcquisitionWithRefocusedEchoes
"""
struct RapidAcquisitionWithRefocusedEchoes <: AnatomicalImage end
const RAREImage = RapidAcquisitionWithRefocusedEchoes
const TurboSpinEcho = RapidAcquisitionWithRefocusedEchoes
const FastSpinEcho = RapidAcquisitionWithRefocusedEchoes

static magnetic field B0;

"""
    Fast low angle shot (FLASH)
"""
struct FastLowAngleShot <: AnatomicalImage end
const FLASH = FastLowAngleShot
modalitylabel(::FLASH) = "FLASH"

struct CombinedPDT2 <: AnatomicalImage end
modalitylabel(::CombinedPDT2) = "PDT2"

struct InplaneT1 <: AnatomicalImage end
modalitylabel(::InplaneT1) = "inplaneT1"

struct InplaneT2 <: AnatomicalImage end
modalitylabel(::InplaneT2) = "inplaneT2"

struct Angiography <: AnatomicalImage end
modalitylabel(::Angiography) = "angio"

struct SusceptibilityWeightedImaging <: AnatomicalImage end
const SWI = SusceptibilityWeightedImaging
modalitylabel(::SWI) = "SWImagandphase"


text:
module LocalMemoize

export @memo, @memoize, get_memos, clear_memos

const MemoizedFunctions = ObjectIdDict()

get_memos(::Type{T}) where T = haskey(MemoizedFunctions, T) ?
    MemoizedFunctions[T] :
    error("No memoized cache found")
get_memos(func) = get_memos(typeof(func))
clear_memos(::Type{T}) where T = haskey(MemoizedFunctions, T) ?
    map(MemoizedFunctions[T]) do cache_table
        for key in keys(cache_table)
            delete!(cache_table, key)
        end
        cache_table
    end :
    error("No memoized cache found")
clear_memos(func) = clear_memos(typeof(func))

macro memo(expr)
    esc(expr)
end

macro memo(args...)
    nothing
end

macro memoize(expr)
    if expr isa Expr && expr.head == :function
        assign_syntax = false
    elseif expr isa Expr && expr.head == :(=) && (expr.args[1].head == :call ||
        expr.args[1].head == :(::) && expr.args[1].args[1].head == :call)
        assign_syntax = true
    else
        error("@memoize must be used on a function declaration.")
    end
    funcdecl = expr.args[1]
    funcrett = nothing
    if funcdecl.head == :(::)
        funcrett = funcdecl.args[2]
        funcdecl = funcdecl.args[1]
    end
    funcname = funcdecl.args[1]
    funcargs = funcdecl.args[2:end]
    funcbody = expr.args[2]
    memo_argnames = Symbol[]
    memo_argtypes = Symbol[]
    argdef_list = []
    argname_list = Symbol[]
    argtype_list = Symbol[]
    check_memocall = arg -> arg isa Expr && arg.head == :macrocall &&
        arg.args[1] == Symbol("@memo") && length(arg.args) > 1
    for arg in funcargs
        if check_memocall(arg)
            tmp = filter(x -> !(x isa LineNumberNode), arg.args)
            length(tmp) != 2 && error("@memo in argument list only must have one and only one argument.")
            argdef = tmp[2]
        else
            argdef = arg
        end
        if argdef isa Symbol
            argname = argdef
            argtype = :Any
        elseif argdef isa Expr && (argdef.head == :(=) || argdef.head == :kw)
            argdef1 = argdef.args[1]
            if argdef1 isa Symbol
                argname = argdef1
                argtype = :Any
            elseif argdef1 isa Expr && argdef1.head == :(::)
                argname = argdef1.args[1]
                argtype = argdef2.args[2]
            end
        elseif argdef isa Expr && argdef.head == :(::)
            argname = argdef.args[1]
            argtype = argdef.args[2]
        elseif argdef isa Expr && argdef.head == :parameters
        else
            error("Invalid argument list")
        end
        if check_memocall(arg)
            push!(memo_argnames, argname)
            push!(memo_argtypes, argtype)
        end
        push!(argdef_list, argdef)
        push!(argname_list, argname)
        push!(argtype_list, argtype)
    end

    lines = 0
    check_memolist = args -> all(arg -> arg isa Symbol, args)
    check_assignment = arg -> arg isa Expr && arg.head == :(=) &&
        arg.args[1] isa Symbol
    memolist = Symbol[]
    for (i, expression) in enumerate(funcbody.args)
        if check_memocall(expression)
            tmp = filter(x -> !(x isa LineNumberNode), expression.args)
            if check_memolist(tmp[2:end])
                append!(memolist, tmp[2:end])
            elseif check_assignment(tmp[2])
                push!(memolist, tmp[2].args[1])
            else
                error("@memo in body must be used on an assignment or a list of variables to be memoized.")
            end
            lines = i
        end
    end
    memoize_return = lines == 0
    if memoize_return && length(memo_argnames) != 0
        error("In function memoization @memo can not be used on arguments.")
    end
    if memoize_return
        memo_argnames = argname_list
        memo_argtypes = argtype_list
    end

    getsymbol = name -> Symbol(string(gensym(), "_", funcname, "_", name))

    func_cachename = getsymbol("memocache")
    cachename = getsymbol("cache")
    index_name = getsymbol("index")

    cache_struct = getsymbol("cache_t")
    dict_type = :(Dict{Tuple{$(memo_argtypes...)}, $(cache_struct)})

    if memoize_return
        original_func_name = getsymbol("original")
        original_func = :($original_func_name($(argdef_list...))=$funcbody)
        func_main = :(
            function $funcname($(argdef_list...))
                $index_name = Base.ht_keyindex($func_cachename, ($(memo_argnames...),))
                if $index_name >= 0
                    @inbounds return $func_cachename.vals[$index_name]
                else
                    $func_cachename[($(memo_argnames...),)] = $original_func_name($(argname_list...))
                end
            end
        )
        funcrett != nothing && map([original_func, func_main]) do func
            func.args[1] = :($(func.args[1])::$funcrett)
        end
        result = quote
            $original_func
            const $cache_struct = Core.Inference.return_type($original_func_name, ($(argtype_list...),))
            const $func_cachename = $dict_type()
            $func_main
            if !haskey(LocalMemoize.MemoizedFunctions, typeof($funcname))
                LocalMemoize.MemoizedFunctions[typeof($funcname)] = Any[]
            end
            push!(LocalMemoize.MemoizedFunctions[typeof($funcname)], $func_cachename)
            $funcname
        end
    else
        gettype_func_name = getsymbol("gettype")
        func1_name = getsymbol("1")
        func2_name = getsymbol("2")

        # For getting memoized type
        gettype_func = :(function $gettype_func_name($(argdef_list...)) end)
        # For first time run
        func1 = :(function $func1_name($func_cachename::$dict_type, $(argdef_list...)) end)
        # For memoized runs
        func2 = :(function $func2_name($cachename::$cache_struct, $(argdef_list...)) end)
        gettype_func_body = gettype_func.args[end].args
        func1_body = func1.args[end].args
        func2_body = func2.args[end].args

        for (i, expression) in enumerate(funcbody.args)
            if check_memocall(expression)
                tmp = filter(x -> !(x isa LineNumberNode), expression.args)
                if !check_memolist(tmp[2:end])
                    push!(gettype_func_body, tmp[2])
                    push!(func1_body, tmp[2])
                end
            else
                push!(gettype_func_body, expression)
                push!(func1_body, expression)
            end
            if i == lines
                push!(gettype_func_body, :(
                    return ($(memolist...),)
                ))
                push!(func1_body, :(
                    $func_cachename[($(memo_argnames...),)] = $cache_struct($(memolist...))
                ))
                for memovar in memolist
                    push!(func2_body, :(
                        $memovar = getfield($cachename, $(QuoteNode(memovar)))
                    ))
                end
            elseif i > lines
                push!(func2_body, expression)
            end
        end
        func_main = :(function $funcname($(argdef_list...))
            $index_name = Base.ht_keyindex($func_cachename, ($(memo_argnames...),))
            if $index_name >= 0
                @inbounds $cachename = $func_cachename.vals[$index_name]
                $func2_name($cachename, $(argname_list...))
            else
                $func1_name($func_cachename, $(argname_list...))
            end
        end)
        funcrett != nothing && map([func1, func2, func_main]) do func
            func.args[1] = :($(func.args[1])::$funcrett)
        end
        result = quote
            let gettype_func = $gettype_func
                types = Core.Inference.return_type(gettype_func, ($(argtype_list...),))
                cache_struct_name = $(QuoteNode(cache_struct))
                struct_code = :(struct $cache_struct_name end)
                argnames = $memolist
                for (argname, typ) in zip(argnames, types.types)
                    push!(struct_code.args[end].args, :($argname::$typ))
                end
                eval(struct_code)
            end
            const $func_cachename = $dict_type()
            $func1
            $func2
            $func_main
            if !haskey(LocalMemoize.MemoizedFunctions, typeof($funcname))
                LocalMemoize.MemoizedFunctions[typeof($funcname)] = Any[]
            end
            push!(LocalMemoize.MemoizedFunctions[typeof($funcname)], $func_cachename)
            $funcname
        end
    end
    esc(result)
end

end # module


text:
struct DNAdataset{S <: Real}
    #=
    data can be
        1. simulated data: Vector{sim_dna_str_w_motif}
        2. fasta that has labeled ground truth
            (those contigugous substring with CAP latters):
            NamedTuple{(:str, :labeled_ground_truth),
                Tuple{Vector{String}, Vector{String}}}
        3. fasta file without labeled ground tuth:
                Vector{String}
    =#
    raw_data::Union{Vector{sim_dna_str_w_motif},
                Vector{String},
                Vector{public_dna_str_w_motif},
                Nothing
                }

    # motif:
    # - if it's a simulated data, then it has a motif type
    # - if data is loaded from fasta as in case 2 or 3 above,
    #   then it is a Nothing type
    motif::motif_type

    # number of data points
    N::Int

    # number of rows is the length of each transformed data_pt
    # number of columns is the total number of data entries
    # currently work for data of the same length only
    dat::Matrix{S}

    """
    Input:
        motif: a motif type that generates the data
        num_data_pts: how many data points to generate
        GPU: (true/false) use GPU to store the wyk-data
            (so can use GPU for training)
    """
    function DNAdataset{S}(motif::motif_type,
                        num_data_pts::Int
        ) where {S <: Real}
        raw_data = sample_backgound_with_motif_multiple(motif, num_data_pts);
            new(
                raw_data,
                motif,
                num_data_pts,
                sim_data_2_dummy(raw_data;F=S)
            )
    end

    """
    Input:
        filepath: fasta filepath

        has_ground_truth: (true/false) the input fasta
            file has each string contained capped-letter
            DNA-alphabets that are labeled as motifs.

        same_len: (true/false) each string in the input
            fasta file has the same length.

        GPU: (true/false) use GPU to store the wyk-data
            (so can use GPU for training)
    """
    function DNAdataset{S}(dna_read::Vector{String},
                        motif::motif_type,
                        raw_data::Vector{public_dna_str_w_motif}
        ) where {S <: Real}
        new(raw_data,
            motif,
            length(dna_read),
            DNA_strings_2_dummy_vectors_same_length(dna_read;F=S)
        );
    end
    function DNAdataset{S}(dna_read::Vector{String}) where {S <: Real}
        new(nothing,
        nothing,
        length(dna_read),
        DNA_strings_2_dummy_vectors_same_length(dna_read; F=S)
        );
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.7

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Random, StatsBase, DataFrames, Plots
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.22"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
begin
	lambda, N = 1, 10^6
	xGrid = 0:6
end

# ╔═╡ 84a73e9a-2165-11eb-3620-d7ddf9c07aff
begin
	expDist = Exponential(1/lambda)
	floorData = counts(convert.(Int,floor.(rand(expDist,N))), xGrid)/N
	geomDist = Geometric(1-MathConstants.e^-lambda)
end

# ╔═╡ 84a76f3c-2165-11eb-0eeb-1bb55507bb51
begin
	plot( xGrid, floorData,
		line=:stem, marker=:circle,
		c=:blue, ms=10, msw=0, lw=4,
		label="Floor of Exponential")
	plot!( xGrid, pdf.(geomDist,xGrid),
		line=:stem, marker=:xcross,
		c=:red, ms=6, msw=0, lw=2,
		label="Geometric", ylims=(0,1),
		xlabel="x", ylabel="Probability")
end

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.22"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═84a73e9a-2165-11eb-3620-d7ddf9c07aff
# ╠═84a76f3c-2165-11eb-0eeb-1bb55507bb51
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
module ABCDGraphGenerator

using Random
using StatsBase

include("pl_sampler.jl")
include("graph_sampler.jl")

end # module


text:
##  Copyright 2021 Myrl Marmarelis
##
##  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 StatsBase, Statistics
using StaticArrays
using PDMats
using Distributions
using LinearAlgebra
using Random
using ProgressMeter: Progress, next!

import Base.@kwdef

abstract type Simulation end

# more of a tactic than a strategy
@kwdef struct RejectionSim <: Simulation
  n_sample_attempts :: Int
end

@kwdef struct HamiltonianSim <: Simulation
  n_trajectories :: Int
  collision_limit :: Int
end

# we shall also use the prior to draw an initial estimate/guess
# parameters {N,M} exist solely to inform the methods downstream
# for multiple issues: support a block-diagonal prior?
@kwdef struct Prior{ N, M, T }
  precision_scale :: AbstractPDMat{T}
  precision_dof :: T
  mean_loc :: AbstractVector{T}
  mean_scale :: T
  dirichlet_weights :: AbstractVector{T}
end

function make_impartial_prior(N::Int, M::Int;
    precision_scale::T, precision_dof::T, mean_scale::T,
    dirichlet_weight::T) where T
  Prior{N,M,T}(
    precision_scale = PDMat(fill(precision_scale, N) |> Diagonal |> Matrix),
    mean_loc = fill(zero(T), N),
    dirichlet_weights = fill(dirichlet_weight, M);
    precision_dof, mean_scale)
end

function simulate_utilities(simulation::RejectionSim, status, voters, vote, mixture, indifference)
  n_local_failures = rejection_sample_utilities!(
    voters, vote, mixture, simulation.n_sample_attempts, indifference)
  n_failures = ( status === nothing ?
    n_local_failures : (status.n_failures + n_local_failures) )
  (; n_failures )
end

# should be more than reasonable for Julia to infer fully concrete types here
function simulate_utilities(simulation::HamiltonianSim, status, voters, vote, mixture, indifference)
  for (voter_index, voter) in enumerate(voters)
    @assert( obeys_ranking(voter.utility, voter_index, vote, indifference), # disable for high-performance sessions?
      "ranking violated: $voter, $voter_index" )
  end
  local_incompleteness = hamiltonian_sample_utilities!(
    voters, vote, mixture, simulation.n_trajectories,
    simulation.collision_limit, indifference)
  incompleteness = ( status === nothing ?
    local_incompleteness : (status.incompleteness + local_incompleteness) )
  (; incompleteness )
end

# we do marginal (univariate) moments, not covariance and higher-order relations..
# note that multivariate normal does not carry any information above second order,
# but the mixture and the predictive posterior could (and almost always do).
# what if I included rolling quantiles for each voter's utility?
function simulate(prior::Prior{N,M,T}, simulation::Simulation,
    vote::Union{MultiIssueVote, IssueVote}, n_rounds::Int; seed::Int,
    n_utility_moments::Int = 1, n_burnin::Int = 0, indifference::Bool = false,
    verbose::Bool = false, gc_interval::Int = 0, do_coincident::Bool = true
    )::NamedTuple where {N,M,T}
  L = N * N
  Random.seed!(seed)
  if vote isa IssueVote
    vote = (vote,) :: MultiIssueVote
  end # sanitize.(vote.choices, vote.n_candidates) would have been fabulous syntax
  vote = map(vote) do issue
    map(choice -> sanitize(choice, issue.n_candidates), issue)
  end
  n_total_rounds = n_rounds + n_burnin
  initial_cohorts = sample_mixture_posteriors(VoterRealization{N,M,T}[];
    prior.precision_scale, prior.precision_dof, prior.mean_loc, prior.mean_scale)
  initial_shares = rand(prior.dirichlet_weights|>Dirichlet) #fill(one(T)/M, M)
  mixture = VoterMixture{N,M,T,L}(initial_cohorts, initial_shares)
  n_voters = length(vote[1].choices)
  voters = map(1:n_voters) do voter
    VoterRealization{N,M,T}(
      rand(initial_shares|>Categorical), zero(Utility{N,T}) )
  end # the above needs to be adjusted to actually properly fit the rankings,
  # so that we're not resting on all the walls (if initialized at zero, and
  # in the Hamiltonian case) and possibly passing right through them...
  #rejection_sample_utilities!(
  #  voters, vote, mixture, Int(1e10), indifference) # make sure no failure occurs?
  voters = map(1:n_voters) do voter
    VoterRealization{N,M,T}(
      rand(initial_shares|>Categorical),
      shoehorn_ranking(@SVector(randn(T, N)),
          voter, vote, indifference) )
  end
  moments = [ zero(Utility{N,T})
    for m in 1:n_utility_moments, v in 1:n_voters ]
  coincident = ( do_coincident ?
    zeros(Int, n_voters, n_voters) : nothing )
  mixtures = VoterMixture{N,M,T}[]
  sim_status = nothing
  progress = Progress(n_total_rounds, enabled=verbose)
  for trial in 1:n_total_rounds
    sim_status = simulate_utilities(simulation, sim_status,
      voters, vote, mixture, indifference)
    if trial <= n_burnin
      status = nothing # don't count them if still in burnin stage
    else
      for v in 1:n_voters
        for m in 1:n_utility_moments
          moments[m, v] += voters[v].utility .^ m
        end
      end
      if do_coincident
        count_coincidences!(coincident, voters)
      end
    end
    shares = sample_mixture_shares(voters, prior.dirichlet_weights)
    cohorts = sample_mixture_posteriors(voters; prior.precision_scale,
      prior.precision_dof, prior.mean_loc, prior.mean_scale)
    mixture = VoterMixture{N,M,T,L}(cohorts, shares)
    push!(mixtures, mixture)
    report = () -> [(:means, mean(mixture))]
    next!(progress, showvalues=report)
    if gc_interval > 0 && trial % gc_interval == 0
      GC.gc(true) # full round to reorganize memory and avoid continuous swapping from fragmented objects in virtual memory?
    end
  end
  sampled_mixtures = @view mixtures[(1+n_burnin):end]
  moments ./= n_rounds
  (; mixtures=sampled_mixtures, coincident, moments, sim_status... )
end


text:
@inline function NNlib._safe_add!(dx::CuDeviceArray{T, 4}, value, ix, iy, c, n) where T
    @inbounds CUDA.@atomic dx[ix, iy, c, n] += value
end

function grid_sample_kernel!(n_elem, output, input, grid, padding_mode)
    index = (threadIdx().x - 1) + (blockIdx().x - 1) * blockDim().x
    if index < n_elem
        iW, iH, iC, _ = size(input)
        _, gW, gH, _ = size(grid)

        w = index % gW + 1
        h = (index ÷ gW) % gH + 1
        n = index ÷ (gW * gH) + 1
        NNlib._grid_sample_kernel!(output, input, grid, padding_mode, w, h, n, iW, iH, iC)
    end
    nothing
end

function ∇grid_sample_kernel!(n_elem, dx, dgrid, Δ, input, grid, padding_mode)
    index = (threadIdx().x - 1) + (blockIdx().x - 1) * blockDim().x
    if index < n_elem
        iW, iH, iC, _ = size(input)
        _, gW, gH, _ = size(grid)

        w = index % gW + 1
        h = (index ÷ gW) % gH + 1
        n = index ÷ (gW * gH) + 1
        NNlib._∇grid_sample_kernel!(dx, dgrid, Δ, input, grid, padding_mode, w, h, n, iW, iH, iC)
    end
    nothing
end

function NNlib.grid_sample(x::CuArray{T, 4}, grid::CuArray{V, 4}; padding_mode = :zeros) where {T, V}
    pad = Val(padding_mode)
    _, _, xC, xN = size(x)
    _, gW, gH, _ = size(grid)
    n_elem = gW * gH * xN
    y = similar(x, T, (gW, gH, xC, xN))

    kernel = @cuda launch=false grid_sample_kernel!(n_elem, y, x, grid, pad)
    config = launch_configuration(kernel.fun; max_threads=256)
    threads = min(n_elem, config.threads)
    blocks = cld(n_elem, threads)
    kernel(n_elem, y, x, grid, pad; threads=threads, blocks=blocks)
    y
end

function NNlib.∇grid_sample(Δ::CuArray{T, 4}, x::CuArray{T, 4}, grid::CuArray{V, 4}; padding_mode = :zeros) where {T, V}
    pad = Val(padding_mode)
    xN = size(x, 4)
    _, gW, gH, _ = size(grid)
    n_elem = gW * gH * xN
    dx, dgrid = CUDA.zeros(T, size(x)), similar(grid)

    kernel = @cuda launch=false ∇grid_sample_kernel!(n_elem, dx, dgrid, Δ, x, grid, pad)
    config = launch_configuration(kernel.fun; max_threads=256)
    threads = min(n_elem, config.threads)
    blocks = cld(n_elem, threads)
    kernel(n_elem, dx, dgrid, Δ, x, grid, pad; threads=threads, blocks=blocks)
    dx, dgrid
end


text:
module CLIMA

include("Utilities/ParametersType/ParametersType.jl")
include("Utilities/PlanetParameters/PlanetParameters.jl")
include("Utilities/RootSolvers/RootSolvers.jl")
include("Utilities/MoistThermodynamics/MoistThermodynamics.jl")
include("Atmos/Parameterizations/SurfaceFluxes/SurfaceFluxes.jl")
include("Atmos/Parameterizations/TurbulenceConvection/TurbulenceConvection.jl")
include("Mesh/Topologies.jl")
include("Mesh/Grids.jl")
include("Arrays/MPIStateArrays.jl")
include("Atmos/Dycore/AtmosDycore.jl")
include("DGmethods/SpaceMethods.jl")
include("DGmethods/DGBalanceLawDiscretizations.jl")
include("ODESolvers/ODESolvers.jl")
include("ODESolvers/LowStorageRungeKuttaMethod.jl")
include("ODESolvers/StrongStabilityPreservingRungeKuttaMethod.jl")
include("ODESolvers/GenericCallbacks.jl")
include("InputOutput/Vtk/Vtk.jl")
include("misc.jl")

end

max_stars_repo_name: jkozdon/CLIMA

max_issues_repo_name: jkozdon/CLIMA

max_forks_repo_name: jkozdon/CLIMA


text:
struct MaxIteration <: AbstractTermination
    K::Int

    function (::Type{MaxIteration})(K::Int)
        return new(K)
    end
end

function terminate(maxiter::MaxIteration,k::Integer,fval::AbstractFloat,x::AbstractVector,g::AbstractVector)
    return k > maxiter.K
end

struct GradientNorm <: AbstractTermination
    ϵ::Float64

    function (::Type{GradientNorm})(ϵ::Float64)
        return new(ϵ)
    end
end

function terminate(gradnorm::GradientNorm,k::Integer,fval::AbstractFloat,x::AbstractVector,g::AbstractVector)
    return k > 1 && norm(g) <= gradnorm.ϵ
end


text:
__precompile__(true)

module VSL

    export BRNGType, VSL_BRNG_MCG31, VSL_BRNG_R250, VSL_BRNG_MRG32K3A, VSL_BRNG_MCG59, VSL_BRNG_WH, VSL_BRNG_SOBOL,
           VSL_BRNG_NIEDERR, VSL_BRNG_MT19937, VSL_BRNG_MT2203, VSL_BRNG_SFMT19937, VSL_BRNG_NONDETERM, VSL_BRNG_ARS5,
           VSL_BRNG_PHILOX4X32X10

    export BasicRandomNumberGenerator

    export VSL_RNG_METHOD_STD, VSL_RNG_METHOD_BOXMULLER, VSL_RNG_METHOD_BOXMULLER2, VSL_RNG_METHOD_ICDF,
           VSL_RNG_METHOD_GNORM, VSL_RNG_METHOD_CJA, VSL_RNG_METHOD_BTPE, VSL_RNG_METHOD_H2PE, VSL_RNG_METHOD_PTPE,
           VSL_RNG_METHOD_POISNORM, VSL_RNG_METHOD_NBAR, VSL_RNG_METHOD_STD_ACCURATE,
           VSL_RNG_METHOD_BOXMULLER2_ACCURATE, VSL_RNG_METHOD_ICDF_ACCURATE, VSL_RNG_METHOD_GNORM_ACCURATE,
           VSL_RNG_METHOD_CJA_ACCURATE

    export MatrixStorageType, VSL_MATRIX_STORAGE_FULL, VSL_MATRIX_STORAGE_PACKED, VSL_MATRIX_STORAGE_DIAGONAL

    export VSLDistribution, VSLContinuousDistribution, VSLDiscreteDistribution,
           Uniform, Gaussian, GaussianMV, Exponential, Laplace, Weibull, Cauchy, Rayleigh, Lognormal,
           Gumbel, Gamma, Beta,
           UniformDiscrete, UniformBits, UniformBits32, UniformBits64, Bernoulli, Geometric, Binomial, Hypergeometric,
           Poisson, PoissonV, NegBinomial

    import MKL_jll
    const libmkl = MKL_jll.libmkl_rt

    const BUFFER_LENGTH = 16

    include("errors.jl")
    include("brngs.jl")
    include("methods.jl")
    include("common.jl")

    include("continuous_distributions.jl")
    include("discrete_distributions.jl")

end


text:
"""
HAMILTONIAN MONTE CARLO TYPE INSTANTIATION

Approximate Bayesian Inference via Hamiltonian Monte Carlo.
"""
immutable HMC
  U       ::Function
  K       ::Function
  dU      ::Function
  dK      ::Function
  init_est::Array{Float64}
  d       ::Int64
end

function mcmc(parameters::HMC;
  leapfrog_params::Dict{Symbol, Real} = Dict([:ɛ => 0.05, :τ => 20]),
  set_seed::Int64 = 123,
  r::Int64 = 1000)

  U, K, dU, dK, w, d = parameters.U, parameters.K, parameters.dU, parameters.dK, parameters.init_est, parameters.d
  ɛ, τ = leapfrog_params[:ɛ], leapfrog_params[:τ]
  H(x::AbstractArray{Float64}, p::AbstractArray{Float64}) = U(x) + K(p)

  if typeof(set_seed) == Int64
    srand(set_seed)
  end

  chain = zeros(r, d);
  chain[1, :] = w
  for i in 1:(r - 1)
    w = chain[i, :]
    p = randn(length(w))
    oldE = H(w, p)

    for j in 1:τ
      p = p - (ɛ / 2) * dU(w)
      w = w + ɛ * dK(p)
      p = p - (ɛ / 2) * dU(w)
    end

    newE = H(w, p)
    dE = newE - oldE

    if dE[1] < 0
      chain[i + 1, :] = w
    elseif rand(Uniform()) < exp(-dE)[1]
      chain[i + 1, :] = w
    else
      chain[i + 1, :] = chain[i, :]
    end
  end

  return chain
end


text:
"""
    dominates(p1::Array, p2::Array)

    Return true if each element in p1 is not less than the corresponding element in p2 and at least one element in p1 is bigger than the corresponding element in p2.

# Arguments
- `p1::Array`: Numeric array of n elements.
- `p2::Array`: Numeric array of n elements.

# Examples
```julia-repl
julia> dominates([1,2,3], [1,2,1])
true
julia> dominates([0,0,0,0], [1,0,0,0])
false
```

# References
Deb, Kalyanmoy, et al. "A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II."
International conference on parallel problem solving from nature. Springer, Berlin, Heidelberg, 2000.
"""
function dominates(p1::Array, p2::Array)::Bool
    n = length(p1)
    notworse = count(i -> p1[i] < p2[i], 1:n)
    better   = count(i -> p1[i] > p2[i], 1:n)
    return (notworse == 0) && (better > 0)
end

"""
    ndsranks(data)

    Sort multidimensional data usin non-dominated sorting algorithm.

# Arguments
- `data::DataFrame`: DataFrame of variables.

# References
Deb, Kalyanmoy, et al. "A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II."
International conference on parallel problem solving from nature. Springer, Berlin, Heidelberg, 2000.
"""
function ndsranks(data::DataFrame)::Array{Int}

    mat = convert(Matrix, data)

    return ndsranks(mat)

end

"""
    ndsranks(data)

    Sort multidimensional data usin non-dominated sorting algorithm.

# Arguments
- `data::Matrix`: n x k matrix of observations where n is number of observations and k is number of variables.

# Examples
```julia-repl

```

# References
Deb, Kalyanmoy, et al. "A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II."
International conference on parallel problem solving from nature. Springer, Berlin, Heidelberg, 2000.
"""
function ndsranks(data::Matrix)::Array{Int64}

    n, p = size(data)

    ranks = zeros(Int64, n)

    mat = convert(Matrix, data)

    @inbounds for i in 1:n
        @inbounds for j in 1:n
            if i != j
                if dominates(mat[i,:], mat[j,:])
                    ranks[i] += 1
                end
            end
        end
    end

    return ranks
end

function nds(data::DataFrame)

    ranks = ndsranks(data)

    bestIndex = sortperm(ranks) |> last

    result = NDSResult(
        ranks,
        bestIndex
    )

    return result
end

max_stars_repo_name: ersagunkuruca/JMcDM

max_issues_repo_name: ersagunkuruca/JMcDM

max_forks_repo_name: ersagunkuruca/JMcDM


text:
## module for pyqtgraph
## More tedious than needed as we list some keywords, rather than pass in through kwargs...
module PyQtGraph

using PyCall
@pyimport pyqtgraph
@pyimport numpy as np
using PySide

export GraphicsWindow, GraphicsLayoutWidget, PlotWidget, ScatterPlotWidget
export ScatterPlotItem
export plot, addPlot, addViewBox, nextRow, removeItem
export pyqtgraph

abstract GraphicsContainer <: PySide.QtWidget

import Base.getindex
getindex(x::GraphicsContainer, i::Symbol) = x.w[i]

type GraphicsWindow <: GraphicsContainer
    w::PyObject
    GraphicsWindow(args...) = new(pyqtgraph.GraphicsWindow(args...))
end

# Used for laying out GraphicsWidgets in a grid.
type GraphicsLayoutWidget <: GraphicsContainer
    w::PyObject
    GraphicsLayoutWidget(parent) = new(pyqtgraph.GraphicsLayoutWidget(PySide.project(parent)))
end

## these should be at QtWidget level
## http://www.pyqtgraph.org/documentation/graphicsItems/graphicslayout.html#pyqtgraph.GraphicsLayout.__init__
resize(widget::GraphicsWindow, width::Int, height::Int) = qinvoke(widget, :resize, width, height)
setWindowTitle(widget::GraphicsWindow, args...) = qinvoke(widget, :setWindowTitle, args...)

plot(widget::GraphicsContainer, args...; kwargs...) = GraphicsPlot(qinvoke(widget, :plot, args...; kwargs...))
addPlot(widget::GraphicsContainer; title="") = GraphicsPlot(PySide.project(widget)[:addPlot](title=title))
addViewBox(widget::GraphicsContainer) = GraphicsPlot(qinvoke(widget, :addViewBox))
## layout
nextRow(widget::GraphicsContainer, args...) = qinvoke(widget, :nextRow, args...)
removeItem(widget::GraphicsContainer, item) = qinvoke(widget, :removeItem, item)

## ??? addViewBox

## PyQtGraphicsPlotItem
## Here we make these methods so basically we have
## w, dev = PlotWidget()
## dev.plot(x=1:10, y=1:20)
type GraphicsPlot <: PySide.QtWidget
    w::PyObject
    plot
    addPoints                   # points? add to API for scatter plot
    addItem
    removeItem
    clear
    addLegend
    addLine
    showGrid
    showAxis
    hideAxis
    setLabel
    setXRange
    setYRange
    enableAutoRange
    setLogMode

    function GraphicsPlot(w)
        self = new(PySide.project(w), nothing)

        ## no kwargs with anonymous functions? Here we have to define generic...
        ## cf. http://www.pyqtgraph.org/documentation/graphicsItems/plotdataitem.html#pyqtgraph.PlotDataItem.__init__
        ## x, x,y, [x y]
        ## http://www.pyqtgraph.org/documentation/graphicsItems/scatterplotitem.html#pyqtgraph.ScatterPlotItem.setData
        ## axis is in ["‘left’, ‘bottom’, ‘right’, or ‘top’]
        plot(args...; pen="b",
             brush=nothing,
             symbol=nothing, symbolPen=nothing, symbolSize=nothing, symbolBrush=nothing,
             fillLevel=nothing) =
          w[:plot](args...,  pen=pyqtgraph.mkPen(pen),
                   brush=pyqtgraph.mkBrush(brush),
                   symbol=symbol, symbolPen=symbolPen, symbolBrush=symbolBrush,
                   fillLevel=fillLevel)
        ## XXX change name; better defaults
        addPoints(args...;  symbol="o", symbolSize=5,   symbolBrush=(255,255,255, 50), pen=nothing, symbolPen=nothing) =
          plot(args...; symbol=symbol, symbolSize=symbolSize, symbolBrush=symbolBrush, pen=pen, symbolPen=symbolPen)
        addItem(item, args...; kwargs...) =    qinvoke(w,:addItem, item, args...; kwargs...)
        removeItem(item, args...; kwargs...) = qinvoke(w,:removeItem, item, args...; kwargs...)
        addLegend(;size=nothing, offset=(30,30)) = qinvoke(w, :addLegend, size=size, offset=offset)
        addLine(;x=nothing, y=nothing, z=nothing, kwargs...) = qinvoke(w,:addLine; x=x,y=y,z=z, kwargs...)
        showGrid(; x=true, y=true, alpha=nothing)  = qinvoke(w, :showGrid, x=x, y=y, alpha=alpha)
        showAxis(where::String, show::Bool=true; kwargs...) = qinvoke(w, :showAxis, where; show=show, kwargs...)
        hideAxis(axis::String) = w[:hideAxis](axis)
        setLabel(axis::String, text=nothing, units=nothing; kwargs...) = qinvoke(w, :setLabel, axis; text=text, units=units, kwargs...)
        setLogMode(args...; x=true, y=true, kwargs...) = qinvoke(w, :setLogMode, args...; x=x, y=y, kwargs...)
        enableAutoRange(args...; kwargs...) = qinvoke(w, :enableAutoRange, args...; kwargs...)
        setXRange(args...; kwargs...) = qinvoke(w, :setXRange, args...; kwargs...)
        setYRange(args...; kwargs...) = qinvoke(w, :setYRange, args...; kwargs...)

        self.plot = plot
        self.addPoints = addPoints
        self.addItem = addItem
        self.removeItem = removeItem
        self.clear = () -> w[:clear]()
        self.addLegend = addLegend
        self.addLine = addLine
        self.showGrid = showGrid
        self.showAxis = showAxis
        self.hideAxis = hideAxis
        self.setLabel = setLabel
        self.setLogMode = setLogMode
        self.enableAutoRange = enableAutoRange
        self.setXRange = setXRange
        self.setYRange = setYRange

        self
    end
end

getindex(x::GraphicsPlot, i::Symbol) = x.w[i]

## A GraphicsViewItem with a single plot
## return (container, device)
## XXX doesn't workPlotWidget(parent) = GraphicsPlot(pyqtgraph.PlotWidget(PySide.project(parent))[:getPlotItem]())
function PlotWidget(parent)
    lyt = GraphicsLayoutWidget(parent)
    plot = addPlot(lyt)
    (lyt, plot)
end

ScatterPlotWidget(parent) = GraphicsPlot(pyqtgraph.ScatterPlotWidget(PySide.project(parent)))

type ScatterPlotItem <: PySide.QtWidget
    w
    clear
    addPoints
    setBrush
    setData
    setPen
    setSize
    setSymbol
    function ScatterPlotItem(args...; pxMode::Bool=true, symbol="o", pen=nothing, brush=nothing, size=10, data=nothing)

        self = new(pyqtgraph.ScatterPlotItem(args..., pxMode=pxMode, symbol=symbol, pen=pen, brush=brush, size=size, data=data))

        addPoints(args...;kwargs...) = qinvoke(w, :addPoints, args...; kwargs...)
        setBrush(args...;kwargs...) = qinvoke(w, :setBrush, args...; kwargs...)
        setData(args...;kwargs...) = qinvoke(w, :setData, args...; kwargs...)
        setPen(args...;kwargs...) = qinvoke(w, :setPen, args...; kwargs...)
        setSize(size; update=true, dataSet=nothing, mask=nothing) = qinvoke(w, :setSize, size, update=update, dataSet=dataSet, mask=mask)
        setSymbol(symbol; update=true, dataSet=nothing, mask=nothing) = qinvoke(w, :setSymbol, symbol, update=update, dataSet=dataSet, mask=mask)
        self.clear = () -> qinvoke(w, :clear)
        self.addPoints = addPoints
        self.setBrush = setBrush
        self.setData = setData
        self.setPen = setPen
        self.setSize = setSize
        self.setSymbol = setSymbol
        self
    end
end

## PlotCurveItem

## signals sigRegionChagned, sigRangeChanged, sigXRangeChangeed, sigYRangeChangeed, ...

## Pen http://www.pyqtgraph.org/documentation/functions.html#pyqtgraph.mkPen
# mkPen(color)
# mkPen(color, width=2)
# mkPen(cosmetic=False, width=4.5, color='r')
# mkPen({'color': "FF0", width: 2})
# mkPen(None)   # (no pen)
mkPen(args...) = pyqtgraph.mkPen(args...)

## Brush: http://www.pyqtgraph.org/documentation/functions.html#pyqtgraph.mkColor
# ‘c’	one of: r, g, b, c, m, y, k, w
# R, G, B, [A]	integers 0-255
# (R, G, B, [A])	tuple of integers 0-255
# float	greyscale, 0.0-1.0
# int	see intColor()
# (int, hues)	see intColor()
# “RGB”	hexadecimal strings; may begin with ‘#’
# “RGBA”
# “RRGGBB”
# “RRGGBBAA”
# QColor	QColor instance; makes a copy.
mkBrush(args...) = pyqtgraph.mkBrush(args...)

end


text:
""" aasd"""
function fooo(x)
    # In 0.6, this could just be a quoted expression instead of a separate file, but
    # in 0.5, `quote` doesn't handle docstrings properly.
    x+2
end


text:
export text, vline, hline, box, colorbar

# NOTE objects don't have the overwrite property. so with or witout ! they don't
# overwrite the axes they get placed on.

"""
    text(text, anchor)

### Example

    text((1, 1.5), t"x^{\\star}"; fontsize=15)
"""
function text(
            anchor::T2R;
            text::String,
            axes::Option{Axes}=nothing,
            o...
        )::PreviewFigure

    axes = check_axes(axes)
    # create object
    t = Text2D(;
            anchor = Float64.(anchor),
            text
    )
    set_properties!(t; o...)
    # add to axes
    push!(axes.objects, t)
    return PreviewFigure(gcf())
end

"""
    _line(...)

Internal function to draw a horizontal or vertical line from axis to axis.
See `vline` and `hline`.
"""
function _line(
            anchor::Float64,
            horiz::Bool;
            axes=nothing,
            o...
        )::PreviewFigure

    axes = check_axes(axes)
    # create object
    l = StraightLine2D(;
            anchor,
            horiz
    )
    set_properties!(l; o...)
    # add to axes
    push!(axes.objects, l)
    return PreviewFigure(gcf())
end

"""
    vline(anchor; opts...)

Add a vertical line at a given point from the xaxis to the x2axis.
Note that this does not adapt the xaxis limits so if the point specified is
outside of the current xaxis limits you won't see anything (until you adjust
the limits accordingly).
See also `hline`.

### Example

    vline(0; color="red")
"""
vline(anchor::Real; o...) = _line(Float64(anchor), false; o...)

"""
    hline(anchor; opts...)

Add a horizontal line at a given point from the yaxis to the y2axis.
See also `vline`.

### Example

    hline(0; color="red")

"""
hline(anchor::Real; o...) = _line(Float64(anchor), true; o...)

"""
    line(from, to; options)

Add a line from `from` (in `(x, y)` format) to `to` (same format).

### Example

    line((0, 0), (1, 1); ls="--")

Will draw a dashed line from (x=0, y=0) to (x=1, y=1).
"""
line(a::T2R, b::T2R; o...) =
    plot!([a[1], b[1]], [a[2], b[2]]; o...)

"""
    box(anchor; opts...)

Add a box at position `anchor`.

### Example

    box((0, 0), (1, 1))

Will draw a box at anchor `(0, 0)` with width `1` and height `1`.
"""
function box(
            anchor::T2R,
            size::T2R;
            axes=nothing,
            o...
        )::PreviewFigure

    axes = check_axes(axes)
    b = Box2D(;
            anchor = Float64.(anchor),
            size   = Float64.(size)
    )
    set_properties!(b; o...)
    push!(axes.objects, b)
    return PreviewFigure(gcf())
end

"""
    colorbar(zmin, zmax, cmap; o...)

Add a colorbar on a side of the axes corresponding to the span `[zmin,zmax]` with a given colormap
`cmap` (vector of colors).
"""
function colorbar(
            zmin::Real,
            zmax::Real,
            cmap::Vector{Color};
            axes::Option{Axes}=nothing,
            o...
        )::PreviewFigure

    axes = check_axes(axes)
    # check if there is an existing colorbar, if so remove it (see later)
    mask = isa.(axes.objects, Colorbar)

    # Create object
    b = Colorbar(;
            zmin = Float64(zmin),
            zmax = Float64(zmax),
            cmap
    )
    # axes adjustment
    if b.position ∈ ("left","right")  # horizontal
        axes.scale = "0.8 auto"
    else
        axes.scale = "auto 0.8"
    end
    set_properties!(b; o...)

    if any(mask)
        axes.objects[findfirst(mask)] = b
    else
        push!(axes.objects, b)
    end
    return PreviewFigure(gcf())
end

"""
    colorbar(h; o...)

Adds a colorbar to a drawing handle generated by a heatmap.

### Example

    h = heatmap(randn(20, 20))
    colorbar(h)

"""
colorbar(dh::DrawingHandle{Heatmap}; o...) =
    colorbar(dh.drawing.zmin, dh.drawing.zmax, dh.drawing.cmap; o...)


text:
using LearningHorse
using Statistics
using Test

@testset "LearningHorse.jl" begin
    include("utils.jl")
    @info "Complete the test for utils"
    test_home_dir = pwd()
    include("Preprocessing.jl")
    @info "Complete the test for Preprocessing"
    cd(test_home_dir)
    include("LossFunction.jl")
    @info "Complete the test for LossFunction"
    include("./Regression.jl")
    @info "Complete the test for Regression"
    include("Classification.jl")
    @info "Complete the test for Classification"
    include("NeuralNetwork.jl")
    @info "Complete the test for NeuralNetwork"
end


text:
using ElectronProgressBars
using Test

@testset "ElectronProgressBars.jl" begin
    # Write your own tests here.
end


text:
@testset "98.validate-binary-search-tree.jl" begin
    @test is_valid_BST(TreeNode{Int}([5, 1, 4, nothing, nothing, 3, 6])) == false
    @test is_valid_BST(TreeNode{Int}([5, 1, 6])) == true
    @test is_valid_BST(TreeNode{Int}([5, 4, 6, nothing, nothing, 3, 7])) == false
    @test is_valid_BST(TreeNode{Int}([2, 2, 2])) == false
end


text:
using RestJuliaDB
using Test
using JSON3

const predicate_test = """{
    "home":{
        "\$eq":"gore"
    },
    "age":{
        "\$gt":25,
        "\$lt":45
    },
    "\$or":[
        {
            "job":{
                "\$eq":"programmer"
            }
        }, {
            "retired":{
                "\$eq":true
            }
        }
    ]
}"""

const goodrow1 = Dict([:home => "gore", :age => 30, :retired => true, :job => "manager"])
const goodrow2 = Dict([:home => "gore", :age => 30, :retired => false, :job => "programmer"])
const badrow1 = Dict([:home => "napier", :age => 30, :retired => true, :job => "manager"]) # wrong home
const badrow2 = Dict([:home => "gore", :age => 60, :retired => true, :job => "manager"]) # wrong age
const badrow3 = Dict([:home => "gore", :age => 30, :retired => false, :job => "manager"]) # not retired and job is wrong

@testset "predicate test" begin
    local pt = JSON3.read(predicate_test, Dict)
    local pred,fields = JuliaDBQuery.predicate_for(pt)
    @test pred(goodrow1) == true
    @test pred(goodrow2) == true
    @test pred(badrow1) == false
    @test pred(badrow2) == false
    @test pred(badrow3) == false
    @test fields ∋ :home
    @test fields ∋ :age
    @test fields ∋ :retired
    @test fields ∋ :job
    @test length(fields) == 4
end

# const sa_filter = """
# {
#     "sa2_code": {
#         "\$eq": 100100
#     },
#     "hour": {
#         "\$eq": 1
#     }
# }
# """

const sa_filter = """
{
    "sa2_code": {
        "\$eq": 100100
    },
    "hour": {
        "\$eq": 12
    }
}
"""

@testset "database test" begin
    local q = JSON3.read(sa_filter, Dict)
    local data = @show JuliaDBQuery.load_table("data")
    local results = @show JuliaDBQuery.run_basic_query(data, q)
    @test length(results) != 0
end


text:
using Random
abstract type Layer end
#This type isn't Optimized.
abstract type NParam <: Layer end
#This type is Optimized.
abstract type Param <: Layer end
export Param, NParam

include("../activations.jl")
include("./Pooling.jl")
include("./Conv.jl")

"""
    NetWork(layers...)

Connect multiple layers, and build a NeuralNetwork.
NetWork also supports index.
You can also add layers later using the add_layer!() Function.

# Example
```
julia> N = NetWork(Dense(10=>5, relu), Dense(5=>1, relu))

julia> N[1]

Dense(IO:10=>5, σ:relu)
```
"""
struct NetWork
    net::Dict
    function NetWork()
        new(Dict())
    end
end

function NetWork(layers...)
    N = NetWork()
    add_layer!(N, layers...)
    return N
end

function Base.getindex(N::NetWork, i::Int)
    return N.net[i]
end

function (N::NetWork)(x; rf=false)
    if !rf
        layers = N.net
        for i in 1 : length(layers)
            x = layers[i](x)
        end
        return x
    else
        record = Dict()
        layers = N.net
        record[1] = (nothing, x)
        for i in 1 : length(layers)
            x = layers[i](x, record, i)
        end
        return record
    end
end

function Base.show(io::IO, N::NetWork)
    layers = N.net
    for i in 1 : length(layers)
        print("Layer$i : ")
        if i != length(layers)
            println(layers[i])
        else
            print(layers[i])
        end
    end
end

"""
    Dense(in=>out, σ; set_w = "Xavier", set_b = zeros)
Crate a traditinal `Dense` layer, whose forward propagation is given by:
    y = σ.(W * x .+ b)
The input of `x` should be a Vactor of length `in`, (Sorry for you couldn't learn using batch. I'll implement)

# Example
```
julia> D = Dense(5=>2, relu)
Dense(IO:5=>2, σ:relu)

julia> D(rand(Float64, 5)) |> size
(2,)
```
"""
struct Dense <: Param
    w::AbstractArray
    b::AbstractArray
    activation
    function Dense(io::Pair{<:Integer, <:Integer}, activation; set_w = "Xavier", set_b = zeros)
        new(dense_w(io..., set_w), set_b(io[2]), activation)
    end
end

function Base.show(io::IO, D::Dense)
    o, i = size(D.w)
    σ = D.activation
    print("Dense(IO:$i => $o, σ:$σ)")
end

function (layer::Dense)(x::AbstractArray)
    w, b, σ = layer.w, layer.b, layer.activation
    σ.(w * x .+ b)
end

function (layer::Dense)(x::AbstractArray, record::Dict, i)
    w, b, σ = layer.w, layer.b, layer.activation
    z = w * x
    a = σ.(z .+ b)
    record[i+1] = (z, a)
    a
end

function (layer::Dense)(Δ, z, back::Bool)
    #TODO:Corresponds to a batch of x.
    w, σ = layer.w, layer.activation
    Σ = zeros(size(w)[2])
    for i in 1 : length(Σ)
        Σ[i] += w[:, i] ⋅ Δ
    end
    println("Σ:", size(Σ))
    println("z:", size(z))
    Σ.*σ.(z, true)
end

"""
    Flatten()
This layer change the dimentions Image to Vector.

# Example
```
julia> F = Flatten()
Flatten(())

julia> F(rand(10, 10, 2, 5)) |> size
(1000, )
```
"""
mutable struct Flatten <: NParam
    csize::Tuple
    function Flatten()
        return new(())
    end
end

function (F::Flatten)(x)
    F.csize = size(x)
    return reshape(x, length(x))
end

function (F::Flatten)(x, record::Dict, i)
    F.csize = size(x)
    a = reshape(x, length(x))
    record[i+1] = (a, a)
    return a
end

(F::Flatten)(Δ, z, back::Bool) = reshape(Δ, F.csize...)

"""
    Dropout(p)
This layer dropout the input data.

# Example
julia> D = Dropout(0.25)
Dropout(0.25)

julia> D(rand(10))
10-element Array{Float64,1}:
 0.0
 0.3955865029078952
 0.8157710047424143
 1.0129613533211907
 0.8060508293474877
 1.1067504108970596
 0.1461289547292684
 0.0
 0.04581776023870532
 1.2794087133638332
"""
struct Dropout <: NParam
    p::Float64
end

function Base.show(io::IO, D::Dropout)
    print("Dropout(")
    print(string(D.p))
    print(")")
end

dropout_kernel(y::T, p, q) where {T} = (y > p) ? float(1 / q) : float(0)

function (D::Dropout)(x)
    y = rand!(similar(x))
    y = dropout_kernel.(y, D.p, 1 - D.p)
    x .* y
end

function (D::Dropout)(x::AbstractArray, record::Dict, i)
    a = D(x)
    record[i+1] = (a, a)
    a
end

(D::Dropout)(Δ, z, back::Bool) = Δ

"""
     add_layer!(model, layers...)
This function add layers to model. You can add layers when you create a NeuralNetwork with `NetWork()`, and you can also use this function to add layers later.

# Example
julia> N = NetWork()

julia> add_layer!(N, Dense(10=>5, relu), Dense(5=>1, relu))

julia> N
Layer1 : Dense(IO:10 => 5, σ:relu)
Layer2 : Dense(IO:5 => 1, σ:relu)
"""
function add_layer!(model::NetWork, obj::T) where T<: Layer
    l = length(model.net)
    model.net[l+1] = obj
end

function add_layer!(model::NetWork, layers...)
    for layer in layers
        add_layer!(model, layer)
    end
end

"""
    @epochs n ex
This macro cruns `ex` `n` times. Basically this is useful for learning NeuralNetwork.

# Example
julia> a = 1

julia> @epochs 1000 a+=1
progress:1000/1000
julia>a
1001
"""
macro epochs(n, ex)
    for i in 1 : n
        progress = "progress:"*string(i)*"/"*string(n)*"\r"
        print(progress)
        eval(ex)
    end
end


text:
"""
    log_metadata!(file::Logfile, metadata[...])

Write top level metadata to the given logfile. This is useful for
describing your whole network in whatever aspect you wish.

```julia
file = Logfile("My Network", tempdir())
# metadata can be specified using keyword arguments
log_metadata!(file, description = "foo", layers = 3, learningrate = 0.001)
# providing a Dict is also possible
log_metadata!(file, Dict("nhidden" => 200, "lambda" = 0.01))
```
"""
log_metadata!(lf::Logfile; args...) = log_metadata!(lf, args)

function log_metadata!(lf::Logfile, args)
    h5open(lf.path, "r+") do root
        for (name_symb, value) in args
            name = string(name_symb)
            a_overwrite(root, name, value)
        end
    end
    nothing
end


text:
module TestParallelRadixSort
using Test

if lowercase(get(ENV, "CI", "false")) == "true" && VERSION ≥ v"1.7-"
    # Use https://github.com/perrutquist/StaticNumbers.jl/pull/12 as otherwise
    # the test would be extremely slow.
    Pkg = Base.require(Base.PkgId(Base.UUID(0x44cfe95a1eb252eab672e2afdf69b78f), "Pkg"))
    Pkg.add(
        Pkg.PackageSpec(;
            name = "StaticNumbers",
            url = "https://github.com/tkf/StaticNumbers.jl.git",
            rev = "kwargs.data",
        ),
    )
end

@testset "$file" for file in sort([
    file for file in readdir(@__DIR__) if match(r"^test_.*\.jl$", file) !== nothing
])
    include(file)
end

end  # module


text:
# lkdv.jl
##########################
using FFTW, DSP, LinearAlgebra, Plots

# initial data
N=256;
n=128;
M=1000;
u=zeros(ComplexF64, N,M);
v=zeros(ComplexF64, N,M);
for p=33:96
    u[p,1]=1/4;
    u[p+128,1]=1/8;
end

# setting
E=zeros(ComplexF64, N, N);
for k=1:N
    E[k,k]=1;
end

t=1/100;

# sqrt(M) times Fourier matrix
omega=exp(im*2*pi/N);
F=zeros(ComplexF64, N, N);
for k=1:N
    for l=1:N
    F[k,l]=omega^(-(k-1)*(l-1));
    end
end

v[:,1]=F*u[:,1];
D1=zeros(ComplexF64, N, N)
D3=zeros(ComplexF64, N, N)
for k=1:n
    D1[k,k]=k-1;
    D1[k+n,k+n]=k-n;
    D3[k,k]=(k-1)^3;
    D3[k+n,k+n]=(k-n)^3;
end

G=F'*(im*8*pi^3*t*N^(-4))*D3*F;
A=E+G+G^2/2+G^3/6;

for l=2:M
    u[:,l]=A*u[:,l-1];
end

v=zeros(ComplexF64, N,M)
for l=1:M
    v[:,l]=u[:,M-l+1];
end

p=plot(real.(u'),
       title="Linearlized KdV flow",
       st=:heatmap,
       color=:jet,
       yaxis="time t",
       xaxis="position x")
plot!(yticks = ([0 1000;], [0 10]))
plot!(xticks = ([0 256;], [0 256]))

q=surface(real.(v),
        title="Linearlized KdV flow",
        camera=(80,60),
        color=:jet,
        xaxis="time t",
        yaxis="position x")
plot!(xticks = ([0 1000;], [10 0]))
plot!(yticks = ([0 256;], [0 256]))

plot(p, q,
     layout=(1,2),
     size=(1600,600),
     margin=Plots.Measures.Length(:mm, 10.0))
savefig("lkdv.png")


text:
@testset "4.2.1.3 (g tan)^p (a+b cos)^m" begin
    (a, b, c, d, e, f, g, m, p, x, ) = @variables a b c d e f g m p x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(g*tan(e+f*x))^p*(a+a*cos(e+f*x))^m=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(g*tan(e+f*x))^p*(a+a*cos(e+f*x))^m=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*tan(e+f*x)^p*(a+a*cos(e+f*x))^m=#

    @test_int [tan(x)^4/(a + a*cos(x)), x, 5, atanh(sin(x))/(2*a) - (sec(x)*tan(x))/(2*a) + tan(x)^3/(3*a)]
    @test_int [tan(x)^3/(a + a*cos(x)), x, 5, -(sec(x)/a) + sec(x)^2/(2*a)]
    @test_int [tan(x)^2/(a + a*cos(x)), x, 4, -(atanh(sin(x))/a) + tan(x)/a]
    @test_int [tan(x)^1/(a + a*cos(x)), x, 4, -(log(cos(x))/a) + log(1 + cos(x))/a]
    @test_int [cot(x)^1/(a + a*cos(x)), x, 5, -(atanh(cos(x))/(2*a)) + (cot(x)*csc(x))/(2*a) - csc(x)^2/(2*a)]
    @test_int [cot(x)^2/(a + a*cos(x)), x, 5, -(cot(x)^3/(3*a)) - csc(x)/a + csc(x)^3/(3*a)]
    @test_int [cot(x)^3/(a + a*cos(x)), x, 6, (3*atanh(cos(x)))/(8*a) - cot(x)^4/(4*a) - (3*cot(x)*csc(x))/(8*a) + (cot(x)^3*csc(x))/(4*a)]
    @test_int [cot(x)^4/(a + a*cos(x)), x, 6, -(cot(x)^5/(5*a)) + csc(x)/a - (2*csc(x)^3)/(3*a) + csc(x)^5/(5*a)]

    @test_int [tan(3*x)/(1 + cos(3*x))^2, x, 3, -(1/(3*(1 + cos(3*x)))) - (1/3)*log(cos(3*x)) + (1/3)*log(1 + cos(3*x))]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*tan(e+f*x)^p*(a+a*cos(e+f*x))^(m/2)=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*(g*tan(e+f*x))^(p/2)*(a+a*cos(e+f*x))^m=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*(g*tan(e+f*x))^(p/2)*(a+a*cos(e+f*x))^(m/2)=#

    #= ::Section:: =#
    #=Integrands*of*the*form*(g*tan(e+f*x))^p*(a+a*cos(e+f*x))^m*(A+B*cos(e+f*x))=#

    #= ::Section:: =#
    #=Integrands*of*the*form*(g*tan(e+f*x))^p*(a+a*cos(e+f*x))^m*(A+B*cos(e+f*x)+C*cos(e+f*x)^2)=#

    #= ::Title:: =#
    #=Integrands*of*the*form*(g*tan(e+f*x))^p*(a+b*cos(e+f*x))^m=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(g*tan(e+f*x))^p*(a+b*cos(e+f*x))^m=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*tan(e+f*x)^p*(a+b*cos(e+f*x))^m=#

    #= ::Subsubsection:: =#
    #=m>0=#

    #= ::Subsubsection::Closed:: =#
    #=m<0=#

    @test_int [tan(x)^4/(a + b*cos(x)), x, 6, (2*(a - b)^(3/2)*(a + b)^(3/2)*atan((sqrt(a - b)*tan(x/2))/sqrt(a + b)))/a^4 + (b*(3*a^2 - 2*b^2)*atanh(sin(x)))/(2*a^4) - ((4*a^2 - 3*b^2)*tan(x))/(3*a^3) - (b*sec(x)*tan(x))/(2*a^2) + (sec(x)^2*tan(x))/(3*a)]
    @test_int [tan(x)^3/(a + b*cos(x)), x, 3, ((a^2 - b^2)*log(cos(x)))/a^3 - ((a^2 - b^2)*log(a + b*cos(x)))/a^3 - (b*sec(x))/a^2 + sec(x)^2/(2*a)]
    @test_int [tan(x)^2/(a + b*cos(x)), x, 6, -((2*sqrt(a - b)*sqrt(a + b)*atan((sqrt(a - b)*tan(x/2))/sqrt(a + b)))/a^2) - (b*atanh(sin(x)))/a^2 + tan(x)/a]
    @test_int [tan(x)^1/(a + b*cos(x)), x, 4, -(log(cos(x))/a) + log(a + b*cos(x))/a]
    @test_int [cot(x)^1/(a + b*cos(x)), x, 3, log(1 - cos(x))/(2*(a + b)) + log(1 + cos(x))/(2*(a - b)) - (a*log(a + b*cos(x)))/(a^2 - b^2)]
    @test_int [cot(x)^2/(a + b*cos(x)), x, 7, -((2*a^2*atan((sqrt(a - b)*tan(x/2))/sqrt(a + b)))/((a - b)^(3/2)*(a + b)^(3/2))) - (a*cot(x))/(a^2 - b^2) + (b*csc(x))/(a^2 - b^2)]
    @test_int [cot(x)^3/(a + b*cos(x)), x, 4, -(((a - b*cos(x))*csc(x)^2)/(2*(a^2 - b^2))) - ((2*a + b)*log(1 - cos(x)))/(4*(a + b)^2) - ((2*a - b)*log(1 + cos(x)))/(4*(a - b)^2) + (a^3*log(a + b*cos(x)))/(a^2 - b^2)^2]
    @test_int [cot(x)^4/(a + b*cos(x)), x, 12, (2*a^4*atan((sqrt(a - b)*tan(x/2))/sqrt(a + b)))/((a - b)^(5/2)*(a + b)^(5/2)) + (a^3*cot(x))/(a^2 - b^2)^2 - (a*cot(x)^3)/(3*(a^2 - b^2)) - (a^2*b*csc(x))/(a^2 - b^2)^2 - (b*csc(x))/(a^2 - b^2) + (b*csc(x)^3)/(3*(a^2 - b^2))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*tan(e+f*x)^p*(a+b*cos(e+f*x))^(m/2)=#

    @test_int [cot(x)/sqrt(3 - cos(x)), x, 5, -atanh(sqrt(3 - cos(x))/2)/2 - atanh(sqrt(3 - cos(x))/sqrt(2))/sqrt(2)]

    @test_int [tan(x)*sqrt(a + b*cos(x)), x, 4, 2*sqrt(a)*atanh(sqrt(a + b*cos(x))/sqrt(a)) - 2*sqrt(a + b*cos(x))]
    @test_int [tan(x)/sqrt(a + b*cos(x)), x, 3, (2*atanh(sqrt(a + b*cos(x))/sqrt(a)))/sqrt(a)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(g*tan(e+f*x))^(p/2)*(a+b*cos(e+f*x))^m=#

    @test_int [sqrt(e*tan(c + d*x))/(a + b*cos(c + d*x)), x, 9, -((2*sqrt(2)*sqrt(cos(c + d*x))*Elliptic.Pi(-(sqrt(-a + b)/sqrt(a + b)), asin(sqrt(sin(c + d*x))/sqrt(1 + cos(c + d*x))), -1)*sqrt(e*tan(c + d*x)))/(sqrt(-a + b)*sqrt(a + b)*d*sqrt(sin(c + d*x)))) + (2*sqrt(2)*sqrt(cos(c + d*x))*Elliptic.Pi(sqrt(-a + b)/sqrt(a + b), asin(sqrt(sin(c + d*x))/sqrt(1 + cos(c + d*x))), -1)*sqrt(e*tan(c + d*x)))/(sqrt(-a + b)*sqrt(a + b)*d*sqrt(sin(c + d*x)))]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*(g*tan(e+f*x))^(p/2)*(a+b*cos(e+f*x))^(m/2)=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(g*tan(e+f*x))^p*(a+b*cos(e+f*x))^m*with*p*symbolic=#

    @test_int [(g*tan(e + f*x))^p*(a + b*cos(e + f*x))^m, x, 1, (g*cot(e + f*x))^p*(g*tan(e + f*x))^p*Unintegrable((a + b*cos(e + f*x))^m/(g*cot(e + f*x))^p, x)]

    #= [(g*tan(e + f*x))^p*(a + b*cos(e + f*x))^3, x, 0, 0]
    @test_int [(g*tan(e + f*x))^p*(a + b*cos(e + f*x))^2, x, 0, 0]
    @test_int [(g*tan(e + f*x))^p*(a + b*cos(e + f*x))^1, x, 0, 0]
    @test_int [(g*tan(e + f*x))^p/(a + b*cos(e + f*x))^1, x, 0, -((g*AppellF1(1 - p, (1 - p)/2, (1 - p)/2, 2 - p, (a + b)/(b + a*sec(e + f*x)), (-a + b)/(b + a*sec(e + f*x)))*(-((a*(1 - sec(e + f*x)))/(b + a*sec(e + f*x))))^((1 - p)/2)*((a*(1 + sec(e + f*x)))/(b + a*sec(e + f*x)))^((1 - p)/2)*(g*tan(e + f*x))^(-1 + p)*(-tan(e + f*x)^2)^((1 - p)/2 + (1/2)*(-1 + p)))/(a*f*(1 - p)))]
    @test_int [(g*tan(e + f*x))^p/(a + b*cos(e + f*x))^2, x, 0, 0]
    @test_int [(g*tan(e + f*x))^p/(a + b*cos(e + f*x))^3, x, 0, 0] =#

    #= ::Section:: =#
    #=Integrands*of*the*form*(g*tan(e+f*x))^p*(a+b*cos(e+f*x))^n*(A+B*cos(e+f*x))=#

    #= ::Section:: =#
    #=Integrands*of*the*form*(g*tan(e+f*x))^p*(a+b*cos(e+f*x))^n*(A+B*cos(e+f*x)+C*cos(e+f*x)^2)=#
    end

max_stars_repo_path: test/4 Trig functions/4.2 Cosine/4.2.1.3 (g tan)^p (a+b cos)^m.jl


text:
# Plot plasma β in a cut plane from 3D PIC outputs.
#
# Hongyang Zhou, hyzhou@umich.edu

using Batsrus, PyPlot, Printf

filename = "3d_box.out"
data = readdata(filename)

sequence = 65 # cut plane index from -
p_ = 18 # thermal pressure index

X = @view data.x[:,:,:,1]
Y = @view data.x[:,:,:,2]
Z = @view data.x[:,:,:,3]

fig, ax = plt.subplots(1, 1)
fig.set_size_inches(3.3, 6)

W  = @view data.w[:,:,:,p_]
Bx = @view data.w[:,:,:,5]
By = @view data.w[:,:,:,6]
Bz = @view data.w[:,:,:,7]

cut1 = @view X[:,sequence,:]
cut2 = @view Z[:,sequence,:]
Pcut  = @view W[:,sequence,:]
Bxcut = @view Bx[:,sequence,:]
Bycut = @view By[:,sequence,:]
Bzcut = @view Bz[:,sequence,:]
PBcut = hypot.(Bx, By, Bz)

c = ax.contourf(cut1, cut2, Pcut./PBcut)
fig.colorbar(c; ax)
ax.axis("scaled")
title(data.head[:wnames][p_])

xlabel("x"); ylabel("z")

dim = [0.125, 0.013, 0.2, 0.045]
str = @sprintf "it=%d, time=%4.2f" data.head[:it] data.head[:time]
at = matplotlib.offsetbox.AnchoredText(str,
   loc="lower left", prop=Dict("size"=>8), frameon=true,
   bbox_to_anchor=(0., 1.),
   bbox_transform=ax.transAxes)
at.patch.set_boxstyle("round,pad=0.,rounding_size=0.2")
ax.add_artist(at)


text:
# Autogenerated wrapper script for qwt_jll for x86_64-linux-gnu
export qwt

using Qt_jll
using GLEW_jll
using GLFW_jll
using GLU_jll
JLLWrappers.@generate_wrapper_header("qwt")
JLLWrappers.@declare_library_product(qwt, "libqwt.so.6.3")
function __init__()
    JLLWrappers.@generate_init_header(Qt_jll, GLEW_jll, GLFW_jll, GLU_jll)
    JLLWrappers.@init_library_product(
        qwt,
        "lib/libqwt.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ImageFiltering, Colors, ColorVectorSpace
using Base.Test

@testset "gradient" begin
    y, x = 1:5, (1:7)'
    for (img, ey, ex) in ((y .* ones(x), 1, 0),
                          (ones(y) .* x, 0, 1),
                          (map(Gray, y.*ones(Float64, x)), Gray(1), Gray(0)),
                          (map(v->RGB(v,0,0), y.*ones(Float64, x)), RGB(1,0,0), RGB(0,0,0)))
        for kernelfunc in (KernelFactors.ando3, KernelFactors.sobel, KernelFactors.prewitt,
                           KernelFactors.ando4, KernelFactors.ando5,
                           Kernel.ando3, Kernel.sobel, Kernel.prewitt,
                           Kernel.ando4, Kernel.ando5)
            gy, gx = @inferred(imgradients(img, kernelfunc, Inner()))
            for val in gy
                @test abs(val - ey) < 1e-4
            end
            for val in gx
                @test abs(val - ex) < 1e-4
            end
            gy, gx = imgradients(img, kernelfunc, Pad{:replicate}())
            @test indices(gy) == indices(gx) == indices(img)
        end
    end
    # 3d
    y, x, z = 1:5, (1:7)', reshape(1:6, 1, 1, 6)
    for (img, ey, ex, ez) in ((y .* ones(x) .* ones(z), 1, 0, 0),
                              (ones(y) .* x .* ones(z), 0, 1, 0),
                              (ones(y) .* ones(x) .* z, 0, 0, 1))
        for kernelfunc in (KernelFactors.ando3, KernelFactors.sobel, KernelFactors.prewitt,
                           Kernel.ando3, Kernel.sobel, Kernel.prewitt)
            gy, gx, gz = @inferred(imgradients(img, kernelfunc, Inner()))
            for val in gy
                @test abs(val - ey) < 1e-4
            end
            for val in gx
                @test abs(val - ex) < 1e-4
            end
            for val in gz
                @test abs(val - ez) < 1e-4
            end
        end
    end
end

nothing


text:
module StatisticalTests

    # packages
    using Printf
    using DataAPI, DataFrames
    using Distributions, Statistics, StatsBase, StatsModels, GLM
    using Reexport
    @reexport using StatsModels
    const D = Distributions
    # end packages

    # imports
    import Base: show
    # import StatsBase: PValue, StatisticalModel, fit, coeftable, tiedrank, countmap
    # import DataAPI: levels
    # import DataFrames: DataFrame
    # import GLM: LinearModel
    import SpecialFunctions: loggamma
    # end imports

    # exports
    export coeftable, fit
    export
        # functions
        t_test,
        f_test,
        ks_test,
        wilcox_test
    # end export

    # types
    const NumVector{T<:Number} = AbstractVector{T}
    const IntVector{Int} = AbstractVector{Int}
    const FloatVector{T<:AbstractFloat} = AbstractVector{T}
    const NumMatrix{T<:Number} = AbstractMatrix{T}
    const IntMatrix{Int} = AbstractMatrix{Int}
    const FloatMatrix{T<:AbstractFloat} = AbstractMatrix{T}
    abstract type TestModel <: StatisticalModel end
    abstract type TestResult end
    # end types

    # includes
    include("utilities.jl")
    include("ttest.jl")
    include("ftest.jl")
    include("kstest.jl")
    include("wilcoxtest.jl")
    # end includes

end # module


text:
# this script would install TestImages and ImageMagick packages if not installed,
# but better you install them in advance manually to avoid ZMQ timeout in LabVIEW
# during installation.

try
    using TestImages, ImageMagick
catch
    import Pkg
    Pkg.add("TestImages")
    Pkg.add("ImageMagick")
end

is2d_img(s) = ! any(occursin.(["_3d_", "series","stack"], lowercase(s)))

function get2d_images()
    rf = TestImages.remotefiles
    rf = filter(is2d_img, rf)
    return (;img_list=rf,)
end

function send_testimage(;img_name)
    img = testimage(img_name)
    img_type = string(typeof(img))
    return (;img_type, bigarrs=(;img))
end

fns = (;get2d_images, send_testimage, )


text:
# Basic properties and constants related to molecular models

module Properties

using ..Dorothy
using ..Dorothy.Utils

export
		standard_atomic_weights, covalent_radii, vertebrate_aa_frequencies,

		wrapid, unwrapids!, unwrapnames!,

		namematcher, isvsite, iswater, isprotein, isacidresidue, isbasicresidue,
		ischargedresidue, ispolarresidue, ishydrophobicresidue, isbackbone,
		ismainchain, issidechain, isnuclacid, islipid, ision, ismonatomicion,
		ispolyatomicion, isalphahelix, ishelix310, ispihelix, isturn, isstrand,
		isbridge, iscoil, isbend, ishelix, issheet, isloop

# J. Meija et al. Pure Appl. Chem., 88(3):265-91, 2016.
# doi:10.1515/pac-2015-0305
const standard_atomic_weights = Dict(
		"H"  => 1.008,
		"He" => 4.002,
		"Li" => 6.94,
		"Be" => 9.012,
		"B"  => 10.81,
		"C"  => 12.01,
		"N"  => 14.00,
		"O"  => 15.99,
		"F"  => 18.99,
		"Ne" => 20.18,
		"Na" => 22.99,
		"Mg" => 24.30,
		"Al" => 26.98,
		"Si" => 28.08,
		"P"  => 30.97,
		"S"  => 32.06,
		"Cl" => 35.45,
		"Ar" => 39.94,
		"K"  => 39.09,
		"Ca" => 40.07,
		"Sc" => 44.95,
		"Ti" => 47.86,
		"V"  => 50.94,
		"Cr" => 51.99,
		"Mn" => 54.93,
		"Fe" => 55.84,
		"Co" => 58.93,
		"Ni" => 58.69,
		"Cu" => 63.54,
		"Zn" => 65.38,
		"Ga" => 69.72,
		"Ge" => 72.63,
		"As" => 74.92,
		"Se" => 78.97,
		"Br" => 79.90,
		"Kr" => 83.79,
		"Rb" => 85.46,
		"Sr" => 87.62,
		"Y"  => 88.90,
		"Zr" => 91.22,
		"Nb" => 92.90,
		"Mo" => 95.95,
		"Ru" => 101.0,
		"Rh" => 102.9,
		"Pd" => 106.4,
		"Ag" => 107.8,
		"Cd" => 112.4,
		"In" => 114.8,
		"Sn" => 118.7,
		"Sb" => 121.7,
		"Te" => 127.6,
		"I"  => 126.9,
		"Xe" => 131.2,
		"Cs" => 132.9,
		"Ba" => 137.3,
		"La" => 138.9,
		"Ce" => 140.1,
		"Pr" => 140.9,
		"Nd" => 144.2,
		"Sm" => 150.3,
		"Eu" => 151.9,
		"Gd" => 157.2,
		"Tb" => 158.9,
		"Dy" => 162.5,
		"Ho" => 164.9,
		"Er" => 167.2,
		"Tm" => 168.9,
		"Yb" => 173.0,
		"Lu" => 174.9,
		"Hf" => 178.4,
		"Ta" => 180.9,
		"W"  => 183.8,
		"Re" => 186.2,
		"Os" => 190.2,
		"Ir" => 192.2,
		"Pt" => 195.0,
		"Au" => 196.9,
		"Hg" => 200.5,
		"Tl" => 204.3,
		"Pb" => 207.2,
		"Bi" => 208.9,
		"Th" => 232.0,
		"Pa" => 231.0,
		"U"  => 238.03)

# B. Cordero et al. Dalton Trans., (21):2832-8, 2008. doi:0.1039/b801115j
const covalent_radii = Dict(
	"H"  => 0.31,
	"He" => 0.28,
	"Li" => 1.28,
	"Be" => 0.96,
	"B"  => 0.84,
	"C"  => 0.76,
	"N"  => 0.71,
	"O"  => 0.66,
	"F"  => 0.57,
	"Ne" => 0.58,
	"Na" => 1.66,
	"Mg" => 1.41,
	"Al" => 1.21,
	"Si" => 1.11,
	"P"  => 1.07,
	"S"  => 1.05,
	"Cl" => 1.02,
	"Ar" => 1.06,
	"K"  => 2.03,
	"Ca" => 1.76,
	"Sc" => 1.70,
	"Ti" => 1.60,
	"V"  => 1.53,
	"Cr" => 1.39,
	"Mn" => 1.39,
	"Fe" => 1.32,
	"Co" => 1.26,
	"Ni" => 1.24,
	"Cu" => 1.32,
	"Zn" => 1.22,
	"Ga" => 1.22,
	"Ge" => 1.20,
	"As" => 1.19,
	"Se" => 1.20,
	"Br" => 1.20,
	"Kr" => 1.16,
	"Rb" => 2.20,
	"Sr" => 1.95,
	"Y"  => 1.90,
	"Zr" => 1.75,
	"Nb" => 1.64,
	"Mo" => 1.54,
	"Tc" => 1.47,
	"Ru" => 1.46,
	"Rh" => 1.42,
	"Pd" => 1.39,
	"Ag" => 1.45,
	"Cd" => 1.44,
	"In" => 1.42,
	"Sn" => 1.39,
	"Sb" => 1.39,
	"Te" => 1.38,
	"I"  => 1.39,
	"Xe" => 1.40,
	"Cs" => 2.44,
	"Ba" => 2.15,
	"La" => 2.07,
	"Ce" => 2.04,
	"Pr" => 2.03,
	"Nd" => 2.01,
	"Pm" => 1.99,
	"Sm" => 1.98,
	"Eu" => 1.98,
	"Gd" => 1.96,
	"Tb" => 1.94,
	"Dy" => 1.92,
	"Ho" => 1.92,
	"Er" => 1.89,
	"Tm" => 1.90,
	"Yb" => 1.87,
	"Lu" => 1.87,
	"Hf" => 1.75,
	"Ta" => 1.70,
	"W"  => 1.62,
	"Re" => 1.51,
	"Os" => 1.44,
	"Ir" => 1.41,
	"Pt" => 1.36,
	"Au" => 1.36,
	"Hg" => 1.32,
	"Tl" => 1.45,
	"Pb" => 1.46,
	"Bi" => 1.48,
	"Po" => 1.40,
	"At" => 1.50,
	"Rn" => 1.50,
	"Fr" => 2.60,
	"Ra" => 2.21,
	"Ac" => 2.15,
	"Th" => 2.06,
	"Pa" => 2.00,
	"U"  => 1.96,
	"Np" => 1.90,
	"Pu" => 1.87,
	"Am" => 1.80,
	"Cm" => 1.69)

# J.L. King, T.H. Jukes. Science, 164(3881):788-98, 1969.
# doi:10.1126/science.164.3881.788
const vertebrate_aa_frequencies = Dict(
		"ALA" => 0.74,
		"ARG" => 0.42,
		"ASN" => 0.44,
		"ASP" => 0.59,
		"CYS" => 0.33,
		"GLU" => 0.58,
		"GLN" => 0.37,
		"GLY" => 0.74,
		"HIS" => 0.29,
		"ILE" => 0.38,
		"LEU" => 0.76,
		"LYS" => 0.72,
		"MET" => 0.18,
		"PHE" => 0.40,
		"PRO" => 0.50,
		"SER" => 0.81,
		"THR" => 0.62,
		"TRP" => 0.13,
		"TYR" => 0.33,
		"VAL" => 0.68)

function wrapid(id::Integer, max::Integer)
	while id > max
		id -= max + 1
	end
	id
end

function unwrapids!(ids::AbstractVector{<:Integer}, max::Integer)
	offset = 0
	atzero = false
	for i in eachindex(ids)
		if ids[i] == 0
			if ! atzero
				offset += max + 1
			end
			atzero = true
		else
			atzero = false
		end
		ids[i] += offset
	end
	ids
end

function unwrapnames!(names::AbstractVector{<:AbstractString})
	for i in eachindex(names)
		name = names[i]
		if length(name) > 0 && isdigit(first(name))
			names[i] = chop(name, head=1, tail=0) * first(name)
		end
	end
	names
end

function namematcher(pattern::AbstractString)
	Istar = findall(==('*'), pattern)
	nstars = length(Istar)
	if nstars == 0
		==(pattern)
	elseif nstars == 1
		if length(pattern) == 1
			name -> true
		else
			if Istar[1] == 1
				suffix = chop(pattern, head=1, tail=0)
				name -> endswith(name, suffix)
			elseif Istar[1] == length(pattern)
				prefix = chop(pattern)
				name -> startswith(name, prefix)
			else
				error("wildcard * may only appear at start or end")
			end
		end
	elseif nstars == 2
		if Istar[1] != 1 || Istar[2] != length(pattern)
			error("wildcard * may only appear at start or end")
		else
			name -> occursin(name, pattern)
		end
	else
		error("wildcard * may only appear at start or end")
	end
end

function namematcher(S::AbstractVector{<:AbstractString})
	F = [namematcher(s) for s in S]
	function (name)
		for f in F
			if f(name)
				return true
			end
		end
		false
	end
end

const vsite_name_pattern = ["MC*", "MN*", "MTRP*", "MW", "LP"]

const water_resname_pattern = ["HOH", "SOL", "TIP*", "WAT*"]

const acid_protein_resname_pattern = ["ASP", "GLU"]

const basic_protein_resname_pattern = ["ARG", "LYS"]

const charged_protein_resname_pattern =
		[acid_protein_resname_pattern..., basic_protein_resname_pattern...]

const polar_protein_resname_pattern =
		["ASN", "CYS", "GLN", "HIS", "SER", "THR", "TRP", "TYR"]

const hydrophobic_protein_resname_pattern =
		["ALA", "GLY", "ILE", "LEU", "MET", "PHE", "PRO", "VAL"]

const protein_resname_pattern =
		[charged_protein_resname_pattern...,
		 polar_protein_resname_pattern...,
		 hydrophobic_protein_resname_pattern...]

const backbone_name_pattern = ["N", "CA", "C"]

const mainchain_name_pattern =
		["H1", "H2", "H3", "N", "H", "CA", "C", "O", "OC1", "OC2", "OXT"]

const nuclacid_resname_pattern = ["A", "C", "G", "U", "DA", "DC", "DG", "DT"]

const lipid_resname_pattern = ["D3PC", "DLPC", "DLPE", "DLPG", "DMPA", "DMPC",
		"DMPG", "DOPC", "DOPE", "DOPG", "DOPS", "DPPC", "DPPE", "DPPG", "DPPS",
		"DSPG", "POPC", "POPE", "POPG", "POPS"]

const monatomic_ion_resname_pattern = ["NA", "CL", "CA", "MG", "K", "RB", "CS",
		"LI", "ZN", "H", "SR", "BA", "AL", "AG", "FE", "CU", "F", "BR", "I",
		"O", "S", "N", "P"]

const polyatomic_ion_resname_pattern =
		["CN", "CO3", "NO2", "NO3", "O2", "OH", "PO3", "PO4", "SO3", "SO4"]

const ion_resname_pattern =
		[monatomic_ion_resname_pattern..., polyatomic_ion_resname_pattern...]

const helix_ss_pattern = ["H", "G", "I", "T"]
const sheet_ss_pattern = ["E", "B"]
const loop_ss_pattern = ["C", "S"]

const isvsite = namematcher(vsite_name_pattern)

const iswater = namematcher(water_resname_pattern)

const isprotein = namematcher(sort(protein_resname_pattern,
		by=(resname -> vertebrate_aa_frequencies[resname])))

const isacidresidue = namematcher(sort(acid_protein_resname_pattern,
		by=(resname -> vertebrate_aa_frequencies[resname])))

const isbasicresidue = namematcher(sort(basic_protein_resname_pattern,
		by=(resname -> vertebrate_aa_frequencies[resname])))

const ischargedresidue = namematcher(sort(charged_protein_resname_pattern,
		by=(resname -> vertebrate_aa_frequencies[resname])))

const ispolarresidue = namematcher(sort(polar_protein_resname_pattern,
		by=(resname -> vertebrate_aa_frequencies[resname])))

const ishydrophobicresidue =
		namematcher(sort(hydrophobic_protein_resname_pattern,
		by=(resname -> vertebrate_aa_frequencies[resname])))

const isbackbonename = namematcher(backbone_name_pattern)

isbackbone(name, resname) = isbackbonename(name) && isprotein(resname)

const ismainchainname = namematcher(mainchain_name_pattern)

ismainchain(name, resname) = ismainchainname(name) && isprotein(resname)

issidechain(name, resname) = !ismainchainname(name) && isprotein(resname)

const isnuclacid = namematcher(nuclacid_resname_pattern)

const islipid = namematcher(lipid_resname_pattern)

const ision = namematcher(ion_resname_pattern)

const ismonatomicion = namematcher(monatomic_ion_resname_pattern)

const ispolyatomicion = namematcher(polyatomic_ion_resname_pattern)

const ishelix = namematcher(helix_ss_pattern)

const isalphahelix = namematcher("H")

const ishelix310 = namematcher("G")

const ispihelix = namematcher("I")

const isturn = namematcher("T")

const issheet = namematcher(sheet_ss_pattern)

const isstrand = namematcher("E")

const isbridge = namematcher("B")

const isloop = namematcher(loop_ss_pattern)

const iscoil = namematcher("C")

const isbend = namematcher("S")

end # module


text:
module AugmentedGPLikelihoods

using Reexport

using ChainRulesCore: @ignore_derivatives
using Distributions
@reexport using GPLikelihoods
using GPLikelihoods: AbstractLikelihood, AbstractLink
using IrrationalConstants
using LogExpFunctions
using MeasureBase
using MeasureTheory: For
using Random: AbstractRNG, GLOBAL_RNG
using SpecialFunctions
using TupleVectors

export nlatent

export init_aux_variables, init_aux_posterior
export aux_sample, aux_sample!, aux_full_conditional
export aux_posterior, aux_posterior!
export auglik_potential,
    auglik_precision, auglik_potential_and_precision, auglik_potential_and_precision
export expected_auglik_potential,
    expected_auglik_precision, expected_auglik_potential_and_precision

export logtilt, expected_logtilt
export aux_prior
export aug_loglik, aux_kldivergence

export ScaledLogistic

export LaplaceLikelihood, NegBinomialLikelihood, StudentTLikelihood

include("api.jl")
include("generic.jl")
include("SpecialDistributions/SpecialDistributions.jl")
using .SpecialDistributions

include("likelihoods/bernoulli.jl")
include("likelihoods/laplace.jl")
include("likelihoods/negativebinomial.jl")
include("likelihoods/poisson.jl")
include("likelihoods/studentt.jl")

include("utils.jl")
include("TestUtils.jl")

end


text:
# Some useful data structures that are used in multiple tests

using DataStructures

using QXContexts

# expected results for ghz exmaple files included in examples/ghz folder
ghz_results = OrderedDict{String, ComplexF32}(
    "11001" => 0 + 0im,
    "10000" => 0 + 0im,
    "00011" => 0 + 0im,
    "00000" => 1/sqrt(2) + 0im,
    "11000" => 0 + 0im,
    "10010" => 0 + 0im,
    "10111" => 0 + 0im,
    "01010" => 0 + 0im,
    "01101" => 0 + 0im,
    "11111" => 1/sqrt(2) + 0im,
)

# sample set of commands used for testing compute graph
sample_cmds = AbstractCommand[
    LoadCommand(:t1, :data_1, [2,2,2]),
    ReshapeCommand(:t1_r, :t1, [[1],[2,3]]),
    LoadCommand(:t2, :data_2, [4,2]),
    ViewCommand(:t1_s, :t1_r, :v1, 2, 4),
    ViewCommand(:t2_s, :t2, :v1, 1, 4),
    OutputCommand(:t3, 1, 2),
    ContractCommand(:t4, [1,3], :t1_s, [1,2], :t2_s, [2,3]),
    ContractCommand(:t5, [2], :t3, [1], :t4, [1,2]),
    SaveCommand(:result, :t5)
]

# matching tensors for sample commands
sample_tensors = Dict{Symbol, AbstractArray}(
    :data_1 => rand(2,2,2),
    :data_2 => rand(4,2)
)


text:
# Constraint relating to network components or quantities not introduced by FlexPlan

## Power balance

# Power balance including candidate storage
function constraint_power_balance_acne_dcne_strg(pm::_PM.AbstractDCPModel, n::Int, i::Int, bus_arcs, bus_arcs_ne, bus_arcs_dc, bus_gens, bus_convs_ac, bus_convs_ac_ne, bus_loads, bus_shunts, bus_storage, bus_storage_ne, pd, qd, gs, bs)
    p                = _PM.var(pm, n, :p)
    pg               = _PM.var(pm, n, :pg)
    pconv_grid_ac_ne = _PM.var(pm, n, :pconv_tf_fr_ne)
    pconv_grid_ac    = _PM.var(pm, n, :pconv_tf_fr)
    pconv_ac         = _PM.var(pm, n, :pconv_ac)
    pconv_ac_ne      = _PM.var(pm, n, :pconv_ac_ne)
    p_ne             = _PM.var(pm, n, :p_ne)
    ps               = _PM.var(pm, n, :ps)
    ps_ne            = _PM.var(pm, n, :ps_ne)
    v                = 1

    JuMP.@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) + sum(pconv_grid_ac[c] for c in bus_convs_ac) + sum(pconv_grid_ac_ne[c] for c in bus_convs_ac_ne)  == sum(pg[g] for g in bus_gens) - sum(ps[s] for s in bus_storage) -sum(ps_ne[s] for s in bus_storage_ne) - sum(pd[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*v^2)
end

# Power balance (without DC equipment) including candidate storage
function constraint_power_balance_acne_strg(pm::_PM.AbstractWModels, n::Int, i::Int, bus_arcs, bus_arcs_ne, bus_gens, bus_loads, bus_shunts, bus_storage, bus_storage_ne, pd, qd, gs, bs)
    p     = _PM.var(pm, n, :p)
    q     = _PM.var(pm, n, :q)
    p_ne  = _PM.var(pm, n, :p_ne)
    q_ne  = _PM.var(pm, n, :q_ne)
    pg    = _PM.var(pm, n, :pg)
    qg    = _PM.var(pm, n, :qg)
    ps    = _PM.var(pm, n, :ps)
    qs    = _PM.var(pm, n, :qs)
    ps_ne = _PM.var(pm, n, :ps_ne)
    qs_ne = _PM.var(pm, n, :qs_ne)
    w     = _PM.var(pm, n, :w, i)

    JuMP.@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) == sum(pg[g] for g in bus_gens) - sum(ps[s] for s in bus_storage) - sum(ps_ne[s] for s in bus_storage_ne) - sum(pd[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*w)
    JuMP.@constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_ne[a] for a in bus_arcs_ne) == sum(qg[g] for g in bus_gens) - sum(qs[s] for s in bus_storage) - sum(qs_ne[s] for s in bus_storage_ne) - sum(qd[d] for d in bus_loads) + sum(bs[s] for s in bus_shunts)*w)
end

# Power balance including candidate storage & flexible demand
function constraint_power_balance_acne_dcne_flex(pm::_PM.AbstractDCPModel, n::Int, i::Int, bus_arcs, bus_arcs_ne, bus_arcs_dc, bus_gens, bus_convs_ac, bus_convs_ac_ne, bus_loads, bus_shunts, bus_storage, bus_storage_ne, gs, bs)
    p                = _PM.var(pm, n, :p)
    pg               = _PM.var(pm, n, :pg)
    pconv_grid_ac_ne = _PM.var(pm, n, :pconv_tf_fr_ne)
    pconv_grid_ac    = _PM.var(pm, n, :pconv_tf_fr)
    pconv_ac         = _PM.var(pm, n, :pconv_ac)
    pconv_ac_ne      = _PM.var(pm, n, :pconv_ac_ne)
    p_ne             = _PM.var(pm, n, :p_ne)
    ps               = _PM.var(pm, n, :ps)
    ps_ne            = _PM.var(pm, n, :ps_ne)
    pflex            = _PM.var(pm, n, :pflex)
    v                = 1

    JuMP.@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) + sum(pconv_grid_ac[c] for c in bus_convs_ac) + sum(pconv_grid_ac_ne[c] for c in bus_convs_ac_ne)  == sum(pg[g] for g in bus_gens) - sum(ps[s] for s in bus_storage) -sum(ps_ne[s] for s in bus_storage_ne) - sum(pflex[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*v^2)
end

# Power balance (without DC equipment) including candidate storage & flexible demand
function constraint_power_balance_acne_flex(pm::_PM.AbstractWModels, n::Int, i::Int, bus_arcs, bus_arcs_ne, bus_gens, bus_loads, bus_shunts, bus_storage, bus_storage_ne, gs, bs)
    p     = _PM.var(pm, n, :p)
    q     = _PM.var(pm, n, :q)
    p_ne  = _PM.var(pm, n, :p_ne)
    q_ne  = _PM.var(pm, n, :q_ne)
    pg    = _PM.var(pm, n, :pg)
    qg    = _PM.var(pm, n, :qg)
    ps    = _PM.var(pm, n, :ps)
    qs    = _PM.var(pm, n, :qs)
    ps_ne = _PM.var(pm, n, :ps_ne)
    qs_ne = _PM.var(pm, n, :qs_ne)
    pflex = _PM.var(pm, n, :pflex)
    qflex = _PM.var(pm, n, :qflex)
    w     = _PM.var(pm, n, :w, i)

    JuMP.@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) == sum(pg[g] for g in bus_gens) - sum(ps[s] for s in bus_storage) - sum(ps_ne[s] for s in bus_storage_ne) - sum(pflex[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*w)
    JuMP.@constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_ne[a] for a in bus_arcs_ne) == sum(qg[g] for g in bus_gens) - sum(qs[s] for s in bus_storage) - sum(qs_ne[s] for s in bus_storage_ne) - sum(qflex[d] for d in bus_loads) + sum(bs[s] for s in bus_shunts)*w)
end

## Candidate AC branches

# Activate a candidate AC branch depending on the investment decisions in the candidate's horizon.
function constraint_ne_branch_activation(pm::_PM.AbstractPowerModel, n::Int, i::Int, horizon::Vector{Int})
    indicator = _PM.var(pm, n, :branch_ne, i)
    investments = _PM.var.(Ref(pm), horizon, :branch_ne_investment, i)

    JuMP.@constraint(pm.model, indicator == sum(investments))
end

# Force investment decisions in two different nws to be the same.
function constraint_ne_branch_investment_same(pm::_PM.AbstractPowerModel, n_1::Int, n_2::Int, i::Int)
    investment_1 = _PM.var(pm, n_1, :branch_ne_investment, i)
    investment_2 = _PM.var(pm, n_2, :branch_ne_investment, i)

    JuMP.@constraint(pm.model, investment_1 == investment_2)
end

## Candidate DC branches

# Activate a candidate DC branch depending on the investment decisions in the candidate's horizon.
function constraint_ne_branchdc_activation(pm::_PM.AbstractPowerModel, n::Int, i::Int, horizon::Vector{Int})
    indicator = _PM.var(pm, n, :branchdc_ne, i)
    investments = _PM.var.(Ref(pm), horizon, :branchdc_ne_investment, i)

    JuMP.@constraint(pm.model, indicator == sum(investments))
end

## Candidate converters

# Activate a candidate AC/DC converter depending on the investment decisions in the candidate's horizon.
function constraint_ne_converter_activation(pm::_PM.AbstractPowerModel, n::Int, i::Int, horizon::Vector{Int})
    indicator = _PM.var(pm, n, :conv_ne, i)
    investments = _PM.var.(Ref(pm), horizon, :conv_ne_investment, i)

    JuMP.@constraint(pm.model, indicator == sum(investments))
end