text:
# Use
#    @warnpcfail precompile(args...)
# if you want to be warned when a precompile directive fails
macro warnpcfail(ex::Expr)
    modl = __module__
    file = __source__.file === nothing ? "?" : String(__source__.file)
    line = __source__.line
    quote
        $(esc(ex)) || @warn """precompile directive
     $($(Expr(:quote, ex)))
 failed. Please report an issue in $($modl) (after checking for duplicates) or remove this directive.""" _file=$file _line=$line
    end
end

const __bodyfunction__ = Dict{Method,Any}()

# Find keyword "body functions" (the function that contains the body
# as written by the developer, called after all missing keyword-arguments
# have been assigned values), in a manner that doesn't depend on
# gensymmed names.
# `mnokw` is the method that gets called when you invoke it without
# supplying any keywords.
function __lookup_kwbody__(mnokw::Method)
    function getsym(arg)
        isa(arg, Symbol) && return arg
        @assert isa(arg, GlobalRef)
        return arg.name
    end

    f = get(__bodyfunction__, mnokw, nothing)
    if f === nothing
        fmod = mnokw.module
        # The lowered code for `mnokw` should look like
        #   %1 = mkw(kwvalues..., #self#, args...)
        #        return %1
        # where `mkw` is the name of the "active" keyword body-function.
        ast = Base.uncompressed_ast(mnokw)
        if isa(ast, Core.CodeInfo) && length(ast.code) >= 2
            callexpr = ast.code[end-1]
            if isa(callexpr, Expr) && callexpr.head == :call
                fsym = callexpr.args[1]
                if isa(fsym, Symbol)
                    f = getfield(fmod, fsym)
                elseif isa(fsym, GlobalRef)
                    if fsym.mod === Core && fsym.name === :_apply
                        f = getfield(mnokw.module, getsym(callexpr.args[2]))
                    elseif fsym.mod === Core && fsym.name === :_apply_iterate
                        f = getfield(mnokw.module, getsym(callexpr.args[3]))
                    else
                        f = getfield(fsym.mod, fsym.name)
                    end
                else
                    f = missing
                end
            else
                f = missing
            end
        else
            f = missing
        end
        __bodyfunction__[mnokw] = f
    end
    return f
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:order, :alias, :options),Tuple{Int64,String,RuleOptions}},Type{Rule},NonTerminal,Array{Any,1}})
    Base.precompile(Tuple{Type{PrepareAnonTerminals},Array{TerminalDef,1},Set{String},Dict{Any,Any},Int64,Nothing})
    Base.precompile(Tuple{Type{TraditionalLexer},LexerConf})
    Base.precompile(Tuple{typeof(_call_userfunc),SimplifyRule_Visitor,Tree})
    Base.precompile(Tuple{typeof(_create_unless),Array{TerminalDef,1},String})
    Base.precompile(Tuple{typeof(_find_used_symbols),Tree})
    Base.precompile(Tuple{typeof(_init_builders),Array{Rule,1},Bool,Bool,Bool})
    Base.precompile(Tuple{typeof(_transform_children),ApplyTemplates,Array{Any,1}})
    Base.precompile(Tuple{typeof(_transform_children),EBNF_to_BNF,Array{Any,1}})
    Base.precompile(Tuple{typeof(_transform_children),PrepareLiterals,Array{Any,1}})
    Base.precompile(Tuple{typeof(_transform_children),PrepareSymbols,Array{Any,1}})
    Base.precompile(Tuple{typeof(_transform_children),RuleTreeToText,Array{Any,1}})
    Base.precompile(Tuple{typeof(calculate_sets),Array{Rule,1}})
    Base.precompile(Tuple{typeof(compute_lalr1_states),LALR_Analyzer})
    Base.precompile(Tuple{typeof(compute_reads_relations),LALR_Analyzer})
    Base.precompile(Tuple{typeof(expand_rule),LALR_Analyzer,NonTerminal,Nothing})
    Base.precompile(Tuple{typeof(maybe_create_child_filter),Array{LarkSymbol,1},Bool,Bool,Nothing})
    Base.precompile(Tuple{typeof(token_func),PrepareGrammar,Val{Symbol("\"234\"")},Token})
    Base.precompile(Tuple{typeof(transformer_func),ApplyTemplates,Val{:template_usage},Meta,Array{Any,1}})
    Base.precompile(Tuple{typeof(transformer_func),EBNF_to_BNF,Val{:expr},Tree,Token})
    Base.precompile(Tuple{typeof(transformer_func),TerminalTreeToPattern,Val{:expansions},Meta,Array{Any,1}})
    Base.precompile(Tuple{typeof(traverse),Tuple{LR0ItemSet,NonTerminal},Array{Any,1},Dict{Any,Int64},Array{Any,1},DefaultDict{Any,Set},DefaultDict{Any,Set},Dict{Any,Any}})
    isdefined(Lerche, Symbol("#128#137")) && Base.precompile(Tuple{getfield(Lerche, Symbol("#128#137")),Pair{String,Rule}})
    isdefined(Lerche, Symbol("#129#138")) && Base.precompile(Tuple{getfield(Lerche, Symbol("#129#138")),String})
    isdefined(Lerche, Symbol("#141#143")) && Base.precompile(Tuple{getfield(Lerche, Symbol("#141#143")),Channel{Any}})
    isdefined(Lerche, Symbol("#153#158")) && Base.precompile(Tuple{getfield(Lerche, Symbol("#153#158")),Channel{Any}})
    isdefined(Lerche, Symbol("#2#3")) && Base.precompile(Tuple{getfield(Lerche, Symbol("#2#3")),Channel{Any}})
    isdefined(Lerche, Symbol("#52#53")) && Base.precompile(Tuple{getfield(Lerche, Symbol("#52#53")),TerminalDef})
    isdefined(Lerche, Symbol("#75#76")) && Base.precompile(Tuple{getfield(Lerche, Symbol("#75#76")),Array{Any,1}})
    isdefined(Lerche, Symbol("#81#82")) && Base.precompile(Tuple{getfield(Lerche, Symbol("#81#82")),Array{Any,1}})
    isdefined(Lerche, Symbol("#get_namespace_name#257")) && Base.precompile(Tuple{getfield(Lerche, Symbol("#get_namespace_name#257")),Token,Nothing})
    let fbody = try __lookup_kwbody__(which(create_callback, (ParseTreeBuilder,))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Nothing,typeof(create_callback),ParseTreeBuilder,))
        end
    end
    let fbody = try __lookup_kwbody__(which(init_analyser!, (LALR_Analyzer,ParserConf,))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Bool,typeof(init_analyser!),LALR_Analyzer,ParserConf,))
        end
    end
end


text:
#!/usr/bin/julia

# Copyright (C) 2021, 2022 Matheus Fernandes Bigolin <mfrdrbigolin@disroot.org>
# SPDX-License-Identifier: MIT

"Day Three, Binary Diagnostic."

include("Utils.jl")

using .Utils

"Calculate a rate with the binary `nums`, operation `op` and index `i`."
calculate_rate(nums, op, i = 1) = op(sum(nums[:,i]), size(nums, 1)//2)

"""Calculate a rating with the binary `nums`, the operation `op` and the index
`i`."""
function calculate_rating(nums, op, i = 1)
    mat_size = size(nums, 1)

    if mat_size == 1
        return nums[1,:]
    end

    rate = calculate_rate(nums, op, i)
    triaged_nums = nums[nums[:,i] .== rate,:]
    return calculate_rating(triaged_nums, op, i + 1)
end

"Given the binary `nums`, calculate the power consumption of the submarine."
function solve1(nums)
    k = size(nums, 2)

    ϵ = [calculate_rate(nums, <, i) for i in 1:k] |> bin2dec
    γ = 2^k - 1 - ϵ

    return ϵ*γ
end

"Given the binary `nums`, calculate the total rating of the submarine."
function solve2(nums)
    O₂ = calculate_rating(nums, >=) |> bin2dec
    CO₂ = calculate_rating(nums, <) |> bin2dec

    return O₂ * CO₂
end

function main()
    usage_and_exit(length(ARGS) != 1)

    input_data = readlines(ARGS[1])
    matrix_data = parse.(Int, hcat(collect.(input_data)...))'

    @assert solve1(matrix_data) == 3985686
    @assert solve2(matrix_data) == 2555739

    println(solve1(matrix_data))
    println(solve2(matrix_data))
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

ext: jl


text:
#============================ Global Parameters ===============================#
const MAX_MOVES = 256
const MAX_PLY = 128
const MATE = 32000
const MAX_QUIET_TRACK = 92
ABORT_SIGNAL = Base.Threads.Atomic{Bool}(false)

#============================== Search Parameters =============================#

const Q_FUTILITY_MARGIN = 193
const RAZOR_DEPTH = 1
const RAZOR_MARGIN = 325
const BETA_PRUNE_DEPTH = 8
const BETA_PRUNE_MARGIN = 172
const SEE_PRUNE_DEPTH = 8
const SEE_QUIET_MARGIN = -190
const SEE_NOISY_MARGIN = -25
const FUTILITY_PRUNE_DEPTH = 8
const FUTILITY_MARGIN = 392
const FUTILITY_MARGIN_NOHIST = 300
const FUTILITY_LIMIT = @SVector [12000, 6000]
const COUNTER_PRUNE_DEPTH = @SVector [4, 3]
const COUNTER_PRUNE_LIMIT = @SVector [0, -500]
const FOLLOW_PRUNE_DEPTH = @SVector [4, 3]
const FOLLOW_PRUNE_LIMIT = @SVector [-1000, -3000]
const WINDOW_DEPTH = 5
const LATE_MOVE_COUNT = @SVector [SVector{10}([0, 2, 4,  7, 11, 16, 22, 29, 37, 46]), SVector{10}([0, 4, 7, 12, 20, 30, 42, 56, 73, 92])]
const LATE_MOVE_PRUNE_DEPTH = 9
const LMRTABLE = init_reduction_table()
const PROBCUT_DEPTH = 5
const PROBCUT_MARGIN = 190

#============================ Piece square tables =============================#

const PAWN_PSQT = SVector{8}([
    SVector{4}([makescore(  0,  0), makescore(  0,  0), makescore(  0,  0), makescore( 0,   0)]),
    SVector{4}([makescore(  -5, -6), makescore( 48, 15), makescore( 18, 27), makescore( 5,  3)]),
    SVector{4}([makescore(-24,-15), makescore( -3, 3), makescore( -19,-11), makescore( -27, -9)]),
    SVector{4}([makescore(-31, 6), makescore(-14, 8), makescore(  -19, 9), makescore( 10, -33)]),
    SVector{4}([makescore( -16, 23), makescore(  14, -2), makescore(0,-11), makescore( -10,-24)]),
    SVector{4}([makescore( 4, 14), makescore(15, 23), makescore(13, -15), makescore( -30, -37)]),
    SVector{4}([makescore( 33,  28), makescore( -6, 25), makescore( 1, 12), makescore(24, -47)]),
    SVector{4}([makescore(  0,  0), makescore(  0,  0), makescore(  0,  0), makescore( 0,   0)])])

const KNIGHT_PSQT = SVector{8}([
    SVector{4}([makescore(-8, -37), makescore(-42, -59), makescore(-96, -14), makescore(-80, 32)]),
    SVector{4}([makescore( -28, -42), makescore(-108, -4), makescore(-1, -56), makescore(-26, 6)]),
    SVector{4}([makescore( -33, 7), makescore(49, -40), makescore( 16,  -15), makescore( -14, 23)]),
    SVector{4}([makescore( -9, 16), makescore( 29, -10), makescore( 66, 14), makescore( 32,  38)]),
    SVector{4}([makescore( 27, 5), makescore( 54, 19), makescore( 13, 49), makescore( 72,  43)]),
    SVector{4}([makescore( 32, -39), makescore( 74, 1), makescore( 94, 2), makescore( 130,  -18)]),
    SVector{4}([makescore( -84, -33), makescore(-32, 10), makescore( 86, -24), makescore( -3, 2)]),
    SVector{4}([makescore(-175,-117), makescore(29, -25), makescore(-46, 9), makescore( 40, 1)])])

const BISHOP_PSQT = SVector{8}([
    SVector{4}([makescore(-103, -41), makescore( -82, 59), makescore(-49,-15), makescore(-53,22)]),
    SVector{4}([makescore(-33, -62), makescore( 3, -58), makescore( 12,-23), makescore(  -9,  6)]),
    SVector{4}([makescore(-8, -16), makescore( -22, 7), makescore( -1,  -48), makescore( 11, 5)]),
    SVector{4}([makescore( 32, 8), makescore( 88,  -83), makescore( 10,  10), makescore( 46, -17)]),
    SVector{4}([makescore( 29, 21), makescore( 21,  1), makescore( 23, 4), makescore( 66, -32)]),
    SVector{4}([makescore( 21, 44), makescore(  60,  -21), makescore( 23,  -20), makescore( 17, -15)]),
    SVector{4}([makescore(-47, -50), makescore(-42, -37), makescore(  -14,  20), makescore( -13,  36)]),
    SVector{4}([makescore(-32, -80), makescore( -10, 31), makescore(-47,-16), makescore(-3,-28)])])

const ROOK_PSQT = SVector{8}([
    SVector{4}([makescore(-25,-17), makescore(-47,12), makescore(35,-20), makescore(31,-37)]),
    SVector{4}([makescore(-93,1), makescore(-58,-4), makescore(-9, -27), makescore( -42, -21)]),
    SVector{4}([makescore(-56, -15), makescore(45,-15), makescore( -44, -22), makescore( -16, -39)]),
    SVector{4}([makescore(-26, -5), makescore( -9, 19), makescore( -14, 2), makescore(16,  -29)]),
    SVector{4}([makescore(-41, 20), makescore(-12, -7), makescore( 43, 21), makescore( -15, 13)]),
    SVector{4}([makescore(-19, 8), makescore( 85, -4), makescore( -4, -1), makescore(-23, 20)]),
    SVector{4}([makescore(  57,  -14), makescore( 20, 17), makescore( 38, 1), makescore(-31, 25)]),
    SVector{4}([makescore( 38, 23), makescore(111, 18), makescore( -13, 48), makescore(-93, 66)])])

const QUEEN_PSQT = SVector{8}([
    SVector{4}([makescore( 0,-70), makescore(-5,-60), makescore(-5,-45), makescore( 5,-25)]),
    SVector{4}([makescore( 0,-55), makescore( 5,-30), makescore(10,-20), makescore(10, -5)]),
    SVector{4}([makescore( 0,-40), makescore( 5,-20), makescore(10,-10), makescore(10,  0)]),
    SVector{4}([makescore( 5,-25), makescore( 5,  0), makescore(10, 10), makescore(10, 25)]),
    SVector{4}([makescore( 0,-30), makescore(15, -5), makescore(10, 10), makescore( 5, 20)]),
    SVector{4}([makescore(-5,-40), makescore(10,-20), makescore( 5,-10), makescore(10,  0)]),
    SVector{4}([makescore(-5,-50), makescore( 5,-25), makescore(10,-25), makescore(10,-10)]),
    SVector{4}([makescore( 0,-75), makescore( 0,-50), makescore( 0,-40), makescore( 0,-35)])])

const KING_PSQT = SVector{8}([
    SVector{4}([makescore(216,  0), makescore(264, 36), makescore(216, 68), makescore(160, 60)]),
    SVector{4}([makescore(224, 44), makescore(240, 80), makescore(188,108), makescore(144,108)]),
    SVector{4}([makescore(156, 72), makescore(208,104), makescore(136,136), makescore( 96,140)]),
    SVector{4}([makescore(132, 80), makescore(152,124), makescore(112,136), makescore( 80,136)]),
    SVector{4}([makescore(124, 76), makescore(144,132), makescore( 84,160), makescore( 56,160)]),
    SVector{4}([makescore(100, 72), makescore(116,136), makescore( 64,148), makescore( 24,152)]),
    SVector{4}([makescore(72,  40), makescore(96,  96), makescore( 52, 92), makescore( 28,104)]),
    SVector{4}([makescore(48,   8), makescore(72,  48), makescore( 36, 60), makescore(  0, 64)])])

const PSQT = @SVector [PAWN_PSQT, KNIGHT_PSQT, BISHOP_PSQT, ROOK_PSQT, QUEEN_PSQT, KING_PSQT]

#=============================== PIECE VALUES =================================#

const PVALS = @SVector [makescore(211, 245), makescore(949, 863), makescore(953, 962), makescore(1220, 1490), makescore(2711, 2840), makescore(15000, 15000)]
const PVALS_MG = SVector{6}(scoreMG.(PVALS))
const PVALS_EG = SVector{6}(scoreEG.(PVALS))

#=============================== Tempo Bonus ==================================#

const TEMPO_BONUS = 22

#============================= Rook Evaluation ================================#

const ROOK_OPEN_FILE_BONUS = makescore(114, 10)
const ROOK_SEMIOPEN_FILE_BONUS = makescore(42, 5)
const ROOK_KING_FILE_BONUS = 10
const ROOK_ON_QUEEN_FILE = makescore(20, 26)

#============================= Pawn Evaluation ================================#

const PAWN_SHIELD_BONUS = makescore(13, 10)
const DOUBLE_PAWN_PENALTY = makescore(34, 84)
const ISOLATED_PAWN_PENALTY = makescore(17, 34)
const PAWN_DEFEND_PAWN_BONUS = 10
const WEAK_PAWN_PENALTY = 25
const PASS_PAWN_THREAT = SVector{7}([makescore(0, 0), makescore(6, 23), makescore(7, 37), makescore(3, 1), makescore(16, 29), makescore(142, 149), makescore(237, 216)])
const CONNECTED_PAWN_PSQT = SVector{8}([
    SVector{4}([makescore(   0,  0), makescore(   0,  0), makescore(   0,  0), makescore(   0,  0)]),
    SVector{4}([makescore(   0,-10), makescore(  10,  0), makescore(   5,  0), makescore(   5, 15)]),
    SVector{4}([makescore(  15,  0), makescore(   30, 0), makescore(  20, 10), makescore(  25, 15)]),
    SVector{4}([makescore(  10,  0), makescore(  25,  5), makescore(  10, 10), makescore(  15, 20)]),
    SVector{4}([makescore(  15,  8), makescore(  20, 15), makescore(  25, 20), makescore(  30, 20)]),
    SVector{4}([makescore(  60, 25), makescore(  51, 50), makescore(  70, 55), makescore(  85, 60)]),
    SVector{4}([makescore( 110,  0), makescore( 205, 10), makescore( 230, 30), makescore( 240, 50)]),
    SVector{4}([makescore(   0,  0), makescore(   0,  0), makescore(   0,  0), makescore(   0,  0)])])
const WEAK_UNOPPOSED = makescore(22, 52)
const WEAK_LEVER = makescore(31, 34)
const BACKWARD_PAWN = makescore(37, 11)

#=========================== Knight Evaluation ================================#

const KNIGHT_TRAP_PENALTY = makescore(50, 50)
const KNIGHT_RAMMED_BONUS = makescore(6, 14)
const KNIGHT_OUTPOST_BONUS = makescore(20, 61)
const KNIGHT_POTENTIAL_OUTPOST_BONUS = makescore(8, 29)

#=========================== Bishop Evaluation ================================#

const BISHOP_TRAP_PENALTY = makescore(80, 80)
const BISHOP_COLOR_PENALTY = makescore(3, 14)
const BISHOP_RAMMED_COLOR_PENALTY = makescore(3, 31)
const BISHOP_PAIR_BONUS = makescore(55, 77)
const BISHOP_OUTPOST_BONUS = makescore(50, 29)
const BISHOP_CENTRAL_CONTROL = makescore(60, 33)

#============================= King Evaluation ================================#

const CASTLE_OPTION_BONUS = makescore(10, 0)
const KING_PAWN_SHIELD_BONUS = 12
const KING_FLANK_ATTACK = makescore(11, 3)
const PAWNLESS_FLANK = makescore(19, 111)

const KING_SHELTER_OFFFILE = SVector{8}([
SVector{8}([makescore(-11, 1), makescore(12, -25), makescore(18, -9), makescore(9, 2), makescore(4, 2), makescore(0, 0), makescore(-4, -31), makescore(-45, 15)]),
SVector{8}([makescore(14, -9), makescore(18, -18), makescore(0, -6), makescore(-14, 2), makescore(-27, 12), makescore(-63, 60), makescore(81, 72), makescore(-24, 0)]),
SVector{8}([makescore(30, -5), makescore(11, -9), makescore(-27, 5), makescore(-11, -9), makescore(-18, -5), makescore(-11, 0), makescore(0, 58), makescore(-12, -3)]),
SVector{8}([makescore(2, 9), makescore(18, -11), makescore(2, -11), makescore(12, -21), makescore(21, -34), makescore(-54, 3), makescore(-124, 45), makescore(3, -8)]),
SVector{8}([makescore(-15, 5), makescore(2, -5), makescore(-25, 0), makescore(-18, 3), makescore(-18, -7), makescore(-38, -2), makescore(28, -17), makescore(-7, -3)]),
SVector{8}([makescore(40, -18), makescore(18, -17), makescore(-20, 0), makescore(-12, -18), makescore(3, -23), makescore(14, -20), makescore(36, -27), makescore(-23, 0)]),
SVector{8}([makescore(21, -13), makescore(-2, -16), makescore(-22, -3), makescore(-18, -9), makescore(-27, -8), makescore(-34, 27), makescore(0, 38), makescore(-9, 0)]),
SVector{8}([makescore(-9, -12), makescore(4, -18), makescore(5, 0), makescore(-3, 9), makescore(-11, 12), makescore(-9, 32), makescore(-171, 77), makescore(-17, 11)])
])

const KING_SHELTER_ONFILE = SVector{8}([
SVector{8}([makescore(0, 0), makescore(-13, -25), makescore(2, -20), makescore(-40, 14), makescore(-24, 0), makescore(1, 40), makescore(-169, -10), makescore(-48, 5)]),
SVector{8}([makescore(0, 0), makescore(22, -23), makescore(5, -9), makescore(-20, -3), makescore(-3, -14), makescore(24, 64), makescore(-186, -5), makescore(-40, 1)]),
SVector{8}([makescore(0, 0), makescore(30, -13), makescore(-4, -9), makescore(5, -19), makescore(13, -9), makescore(-89, 45), makescore(-86, -75), makescore(-20, -7)]),
SVector{8}([makescore(0, 0), makescore(-5, 7), makescore(-4, 0), makescore(-19, 0), makescore(-29, 0), makescore(-100, 30), makescore(5, -43), makescore(-24, -9)]),
SVector{8}([makescore(0, 0), makescore(10, 0), makescore(10, -7), makescore(12, -13), makescore(12, -28), makescore(-59, 13), makescore(-106, -63), makescore(-3, -9)]),
SVector{8}([makescore(0, 0), makescore(4, -10), makescore(-20, 0), makescore(-29, -8), makescore(15, -26), makescore(-40, 1), makescore(53, 36), makescore(-20, -8)]),
SVector{8}([makescore(0, 0), makescore(20, -17), makescore(10, -15), makescore(-10, -9), makescore(-29, 7), makescore(-10, 13), makescore(-58, -50), makescore(-33, 10)]),
SVector{8}([makescore(0, 0), makescore(10, -40), makescore(17, -29), makescore(-20, -6), makescore(-19, 16), makescore(-7, 20), makescore(-230, -57), makescore(-24, 0)])
])

const KING_STORM_UNBLOCKED = SVector{4}([
SVector{8}([makescore(-7, 26), makescore(116, -10), makescore(-28, 23), makescore(-22, 5), makescore(-17, 0), makescore(-10, -7), makescore(-19, 2), makescore(-25, -5)]),
SVector{8}([makescore(-6, 46), makescore(54, 9), makescore(-22, 21), makescore(-8, 8), makescore(-7, 2), makescore(2, -7), makescore(-4, 0), makescore(-14, 0)]),
SVector{8}([makescore(5, 35), makescore(14, 20), makescore(-26, 18), makescore(-14, 8), makescore(0, 0), makescore(4, 0), makescore(8, -10), makescore(0, -4)]),
SVector{8}([makescore(-5, 22), makescore(13, 18), makescore(-19, 8), makescore(-17, 0), makescore(-16, 0), makescore(4, -14), makescore(0, -10), makescore(-16, 0)])
])

const KING_STORM_BLOCKED = SVector{4}([
SVector{8}([makescore(0, 0), makescore(-18, -19), makescore(-20, -4), makescore(15, -20), makescore(7, -10), makescore(0, -20), makescore(-6, -4), makescore(14, 26)]),
SVector{8}([makescore(0, 0), makescore(-19, -37), makescore(-6, -11), makescore(32, -14), makescore(-4, -5), makescore(10, -26), makescore(-10, -13), makescore(-20, 0)]),
SVector{8}([makescore(0, 0), makescore(-31, -52), makescore(-30, -9), makescore(8, -14), makescore(0, -5), makescore(-11, -17), makescore(-16, -18), makescore(-14, 3)]),
SVector{8}([makescore(0, 0), makescore(-5, -20), makescore(-19, -20), makescore(-14, -7), makescore(-7, -10), makescore(1, -28), makescore(70, -14), makescore(10, 18)])
])

#============================ Queen Evaluation ================================#

#============================ Other Evaluation ================================#

const SPACE_BONUS = 4
const PIN_BONUS = 21
const HANGING_BONUS = makescore(6, 19)
const THREAT_BY_MINOR = @SVector [makescore(16, 24), makescore(60, 86), makescore(42, 74), makescore(131, 7), makescore(89, 123)]
const THREAT_BY_ROOK = @SVector [makescore(10, 31), makescore(58, 36), makescore(44, 88), makescore(1, 6), makescore(120, 45)]
const LAZY_THRESH = 2000
const MINOR_KING_PROTECTION = makescore(8, 6)
const RESTRICTION_BONUS = makescore(19, 2)
const KNIGHT_SAFE_CHECK = 182
const QUEEN_SAFE_CHECK  = 148
const BISHOP_SAFE_CHECK = 87
const ROOK_SAFE_CHECK   = 137
const KD_QUEEN  = 44
const KD_ROOK   = 9
const KD_BISHOP = 9
const KD_KNIGHT = 1
const THREAT_BY_PAWN = makescore(129, 28)
const THREAT_BY_PUSH = makescore(62, 43)
const THREAT_BY_KING = makescore(8, 26)

const KING_FLANK_DEFEND = makescore(10, 1)
const KING_BOX_WEAK = makescore(7, 16)
const UNSAFE_CHECK = makescore(18, 0)
const KNIGHT_ON_QUEEN = makescore(15, 10)

#========================= Mobility Evaluation ================================#

const KNIGHT_MOBILITY = @SVector [makescore(-60,-80), makescore(-50,-55), makescore(-10,-30), makescore( -5,-15),
    makescore(  5,  10), makescore( 15, 15), makescore( 20, 25), makescore( 30, 30), makescore( 35, 35)]

const BISHOP_MOBILITY = @SVector [makescore(-50,-60), makescore(-20,-25), makescore( 15, -5), makescore( 25, 15), makescore( 40, 25),
    makescore( 50, 40), makescore( 55, 55), makescore( 65, 60), makescore( 65, 65), makescore( 70, 75), makescore( 80, 80),
    makescore( 80, 85), makescore( 90, 90), makescore( 100, 100)]

const ROOK_MOBILITY = @SVector [makescore(-60,-75), makescore(-30,-20), makescore(-15, 30), makescore(-10, 55), makescore( -5, 70), makescore( 0, 80),
      makescore(  10,110), makescore( 15,120), makescore( 30,130), makescore( 30,140), makescore( 30,155), makescore( 40,165),
      makescore( 45,165), makescore( 50,170), makescore( 60,170)]

const QUEEN_MOBILITY = @SVector [makescore(-40,-35), makescore(-20,-15), makescore(  5,  10), makescore(  5, 20), makescore( 15, 35), makescore( 20, 55),
      makescore( 30, 60), makescore( 40, 75), makescore( 45, 80), makescore( 50, 90), makescore( 55, 95), makescore( 60,105),
      makescore( 60,115), makescore( 65,120), makescore( 70,125), makescore( 70,125), makescore( 70,130), makescore( 75,135),
      makescore( 80,140), makescore( 90,145), makescore( 90,150), makescore( 100,165), makescore(100,170), makescore(100,175),
      makescore(105,185), makescore(110,190), makescore(115,205), makescore(115,210)]

#============================ SCALING FACTORS =================================#

const SCALE_OCB_BISHOPS = 64
const SCALE_OCB_ONE_KNIGHT = 106
const SCALE_OCB_ONE_ROOK = 96
const SCALE_DRAW = 0
const SCALE_NORMAL = 128

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


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

using CASIM

render_option = Not_Render()
sim_horizon = 5000
num_run = 10
init_num_ac = 100
# spawn_rates = [1000., 2000., 4000., 6000.] ./ 3600
spawn_rates = [250., 500., 1000., 1500., 2000., 4000.] ./ 3600
# spawn_rates = 4000. / 3600

"""
	VICAS + Multi
  spawn_rates
"""
exp = "SpawnRates_VICAS_10km_init_ac_" * string(init_num_ac) * "/" * VICAS_policyname # to make sub-folder

if !ispath("../simulation_results/SpawnRates_VICAS_10km_init_ac_" * string(init_num_ac))
	mkdir("../simulation_results/SpawnRates_VICAS_10km_init_ac_" * string(init_num_ac))
end
cas = VICAS()
cas_ctrl = Uniform_CAS(cas)
sensor = Multi_Threat_Tracker()
sensor_ctrl = Uniform_Sensor(sensor)

# for run in 1:num_run
for run in 7:num_run
	for spawn_rate in spawn_rates
		scheme = "Multi_sr_" *
			string(Int(spawn_rate * 3600)) * "_run_" * string(run)
		airspace_ctrl = Constant_Spawn_Rate_Controller(init_num_ac, spawn_rate)
		simulate(exp, scheme, sim_horizon, airspace_ctrl, cas_ctrl, sensor_ctrl, render_option)
	end
end

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
module Checks
    export Action, Image, Annotation, Output

    struct Action
        label::String
        description::String
        identifier::String
    end

    struct Image
        alt::String
        image_url::String
        caption::Union{String, Nothing}
    end

    struct Annotation
        filename::String
        blob_href::String
        start_line::Int
        end_line::Int
        warning_level::String
        message::String
        title::Union{String, Nothing}
        raw_details::String # Documented as Union{String, Nothing}, but that errors
    end

    struct Output
        title::String
        summary::String
        text::String # Documented as Union{String, Nothing}, but that errors
        annotations::Vector{Annotation}
        images::Vector{Image}
    end

end
using .Checks

struct CheckRun <: GitHubType
    id::Union{Int, Nothing}
    head_sha::Union{String, Nothing}
    external_id::Union{String, Nothing}
    status::Union{String, Nothing}
    conclusion::Union{String, Nothing}
    started_at::Union{DateTime, Nothing}
    completed_at::Union{DateTime, Nothing}
    app::Union{App, Nothing}
    pull_requests::Union{Vector{PullRequest}, Nothing}
end
CheckRun(data::Dict) = json2github(CheckRun, data)
namefield(cr::CheckRun) = cr.id

@api_default function create_check_run(api::GitHubAPI, repo::Repo; headers = Dict(), kwargs...)
    headers["Accept"] = "application/vnd.github.antiope-preview+json"
    result = gh_post_json(api, "/repos/$(name(repo))/check-runs"; headers=headers, kwargs...)
    return CheckRun(result)
end

@api_default function update_check_run(api::GitHubAPI, repo::Repo, id::Int; headers = Dict(), kwargs...)
    headers["Accept"] = "application/vnd.github.antiope-preview+json"
    result = gh_patch_json(api, "/repos/$(name(repo))/check-runs/$(id)"; headers=headers, kwargs...)
    return CheckRun(result)
end


text:
function interpolate(obj::BivariateGridInterpolator, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::Vector{Vector{jdouble}})
    return jcall(obj, "interpolate", BivariateFunction, (Vector{jdouble}, Vector{jdouble}, Vector{Vector{jdouble}}), arg0, arg1, arg2)
end


text:
# Autogenerated wrapper script for unixODBC_jll for armv6l-linux-gnueabihf
export dltest, isql, iusql, libodbc, libodbccr, libodbcinst, odbc_config, odbcinst, slencheck

using Libiconv_jll
JLLWrappers.@generate_wrapper_header("unixODBC")
JLLWrappers.@declare_library_product(libodbc, "libodbc.so.2")
JLLWrappers.@declare_library_product(libodbccr, "libodbccr.so.2")
JLLWrappers.@declare_library_product(libodbcinst, "libodbcinst.so.2")
JLLWrappers.@declare_executable_product(dltest)
JLLWrappers.@declare_executable_product(isql)
JLLWrappers.@declare_executable_product(iusql)
JLLWrappers.@declare_executable_product(odbc_config)
JLLWrappers.@declare_executable_product(odbcinst)
JLLWrappers.@declare_executable_product(slencheck)
function __init__()
    JLLWrappers.@generate_init_header(Libiconv_jll)
    JLLWrappers.@init_library_product(
        libodbc,
        "lib/libodbc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libodbccr,
        "lib/libodbccr.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libodbcinst,
        "lib/libodbcinst.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        dltest,
        "bin/dltest",
    )

    JLLWrappers.@init_executable_product(
        isql,
        "bin/isql",
    )

    JLLWrappers.@init_executable_product(
        iusql,
        "bin/iusql",
    )

    JLLWrappers.@init_executable_product(
        odbc_config,
        "bin/odbc_config",
    )

    JLLWrappers.@init_executable_product(
        odbcinst,
        "bin/odbcinst",
    )

    JLLWrappers.@init_executable_product(
        slencheck,
        "bin/slencheck",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module OptimizeQCQP

using LinearAlgebra
using SparseArrays
using Arpack

export minimize, Problem

include("linalg.jl")

function minimize(Q0::AbstractMatrix{T}, g0::AbstractVector{T},
                  ::Nothing, ::Nothing, Q1::Union{AbstractMatrix{T},UniformScaling},
                  ::Nothing, ::Nothing, cu::Real; tol=T(1e-4)) where T
    # Solve the problem:
    #    minimize x'*Q0*x/2 + g0'*x
    #    subject to x'*Q1*x/2 <= cu

    # Pick a safe value of λ that makes Q0 + λ*Q1 positive-definite.
    # Keep track of the smallest known such value in λp
    λ, C = posdefλ(Q0, Q1)
    λp = λ

    ## Construct an initial estimate
    x = C \ (-g0)
    xQx = x'*(Q1*x)
    if xQx > 2*cu
        # Since λ only guarantees positivity and does not even try to satisfy the cu
        # constraint, it's possible that this λ is far too small. Try to estimate the right
        # value in terms of a correction Δλ where
        #          x = (C + Δλ*Q1) \ (-g)
        # and then picking Δλ so that x satisfies the constraint.
        # We can expand that inverse in two ways:
        # (1) If C dominates Δλ*Q1, use the expansion inv(C + Δλ*Q1) ≈ inv(C) - Δλ*inv(C)*Q1*inv(C)
        Δx = - (C \ (Q1 * x))  # coefficient of Δλ in correction to x
        # Solve for Δλ from (x + Δλ*Δx)'*Q1*(x + Δλ*Δx) == 2*cu
        a = Δx'*(Q1*Δx)       # coefficient of Δλ^2
        b = 2*x'*(Q1*Δx)      # coefficient of Δλ (guaranteed negative)
        c = xQx - 2*cu        # const (coefficient of Δλ^0)
        s2 = b^2 - 4*a*c
        if s2 >= 0
            ΔλC = 2*c/(-b + sqrt(s2))   # choose the smaller of the two roots
        else
            # it doesn't completely reach the criterion, so choose the minimum
            # (this is likely a sign that solution (2) will be selected)
            ΔλC = -b/(2*a)
        end
        @assert ΔλC >= 0
        xC = x + ΔλC*Δx   # this is our estimate of the solution, given this Δλ
        # (2) If Δλ*Q1 dominates C, then just use x ≈ (1/Δλ) * (Q1 \ (-g))  (dropping C altogether)
        v1 = Q1 \ (-g0)
        Δλ1 = sqrt(-(g0'*v1)/(2*cu))   # initial alternative estimate of Δλ, discounting higher-order Δλ
        x1 = v1 ./ Δλ1         # the alternate estimate of x
        # Pick between them by assessing accuracy in matching
        #     (C + Δλ*Q1) * x = -g0
        errC = Q0*xC .+ (λ + ΔλC).*(Q1*xC) .+ g0
        err1 = Q0*x1 .+ (λ + Δλ1).*(Q1*x1) .+ g0
        absprod(x, y) = abs(x)*abs(y)
        errmagC = sum(absprod(x,y) for (x,y) in zip(errC, xC))
        errmag1 = sum(absprod(x,y) for (x,y) in zip(err1, x1))
        # λ = (s2 < 0 || errmagC > errmag1) ? λ + Δλ1 : λ + ΔλC
        # Use interpolation, weighted by the error
        λ += (errmag1*ΔλC + errmagC*Δλ1)/(errmag1 + errmagC)
        C = cholesky(Q0 + λ*Q1)
        x = C \ (-g0)
        xQx = x'*(Q1*x)
    end

    ## Newton's method iterative improvement
    # Algorithm 4.3 in Nodecal & Wright, generalized for Q1 != I
    # Also written to make it globally convergent
    iter = 0
    while iter < 10
        λp = min(λ, λp)
        xQxold = xQx
        q = Q1*x
        dxdλ = - (C \ q)
        Δ = cu - xQx/2
        Δλ = Δ/cu * (xQx / (2*(q'*dxdλ)))
        λold = λ
        # Global convergence iteration
        itergc = 0
        while itergc < 10
            λ = max(λold + Δλ, zero(T))
            # posdef iteration
            C = cholesky(Q0 + λ*Q1; check=false)
            iterc = 0
            while !issuccess(C) && iterc < 10
                iterc += 1
                λ = (λ + iterc*λp)/(iterc+1)   # move towards a known-safe value
                C = cholesky(Q0 + λ*Q1; check=false)
            end
            if !issuccess(C)
                λ = (λold + λp)/2
                C = cholesky(Q0 + λ*Q1)
            end
            x = C \ (-g0)
            xQx = x'*(Q1*x)
            (xQx <= 2*cu || abs(xQx - 2*cu) < abs(xQxold - 2*cu)) && break
            itergc += 1
            Δλ /= itergc+1
        end
        abs(xQx - xQxold) < tol*(xQx + xQxold) && break
        iter += 1
    end
    return x, λ
end

# Junk from attempt to incorporate higher-order correction
        # (2) If Δλ*Q1 dominates C, use the expansion
        #        inv(C + Δλ*Q1) ≈ inv(Δλ*Q1) - inv(Δλ*Q1)*C*inv(Δλ*Q1)
        #                       = (1/Δλ - λ/Δλ^2)*inv(Q1) - (1/Δλ^2)*inv(Q1)*Q0*inv(Q1)
        # v2 = Q1 \ (Q0 * v1)
        # v1Qv1, v2Qv2, v1Qv2 = v1'*Q1*v1, v2'*Q1*v2, v1'*Q1*v2
        # let λ=λ, v1Qv1=v1Qv1, v2Qv2=v2Qv2, v1Qv2=v1Qv2, cu=cu
        #     froot(Δλ) = (1/Δλ^2) * ((1-λ/Δλ)^2*v1Qv1 - 2*(1-λ/Δλ)/Δλ*v1Qv2 + 1/Δλ^2*v2Qv2) - 2*cu
        #     Δλ1 = bisectroot(froot, zero(Δλ1), Δλ1)
        # end
        # x1 = (1/Δλ1-λ/Δλ1^2).*v1 - v2./Δλ1^2         # the alternate estimate of x

# function bisectroot(f, a, b; tol=1e-3, maxiter=10, maxseek=20)
#     fa, fb = f(a), f(b)
#     f0 = min(abs(fa), abs(fb))
#     fa == 0 && return a
#     fb == 0 && return b
#     iter = 0
#     while fa*fb > 0 && iter < maxseek
#         a, fa = b, fb
#         b *= 2
#         fb = f(b)
#         fb == 0 && return b
#         iter += 1
#     end
#     iter = 0
#     while iter < maxiter
#         c = (a+b)/2
#         fc = f(c)
#         (fc == 0 || (abs(fc - fa) + abs(fc - fb)) < tol*f0) && return c
#         if sign(fc) == sign(fa)
#             a, fa = c, fc
#         else
#             b, fb = c, fc
#         end
#         iter += 1
#     end
#     return b
# end

# Solve the problem:
#    (approximately) minimize x'*Q0*x/2 + g0'*x
#    subject to cl[i] <= x'*Qi*x/2 + gi'*x <= cu[i]       for i = 1,...,m
#                l[k] <= x[k] <= u[k]                     for k = 1,...,n
#
# Note that equality constraints are handled by setting `cl[i] = cu[i]` and/or
# `l[k] = u[k]`.
#
# For positive-semidefinite Q0 and Qi, the problem is convex, otherwise it is NP-hard
# (hence the "approximate" minimization).

struct Problem{T<:Real,Q0T<:AbstractMatrix{T},g0T<:AbstractVector{T}}
    n::Int
    l::Vector{T}
    u::Vector{T}
    cl::Vector{T}
    cu::Vector{T}
    Q0::Q0T
    g0::g0T
    Qidense::Vector{Matrix{T}}
    Qisparse::Vector{SparseMatrixCSC{T,Int}}
    gidense::Vector{Vector{T}}
    gisparse::Vector{SparseVector{T,Int}}
    issparseQi::Vector{Bool}
    issparsegi::Vector{Bool}
    isposdefQi::Vector{Bool}
end

function Problem(Q0::AbstractMatrix{T}, g0::AbstractVector{T},
                 lower=fill(-T(Inf), length(g0)), upper=fill(T(Inf), length(g0))) where T
    n = length(g0)
    indices(g0) == (Base.OneTo(n),) || throw(DimensionMismatch("g0 indices must start with 1, got $(indices(g0))"))
    indices(Q0) == (Base.OneTo(n), Base.OneTo(n))  || throw(DimensionMismatch("Q0 indices must match g0, got $(indices(g0)) and $(indices(Q0))"))
    indices(lower) == (Base.OneTo(n),) || throw(DimensionMismatch("lower's indices must match g0, got $(indices(lower)) and $(indices(g0))"))
    indices(upper) == (Base.OneTo(n),) || throw(DimensionMismatch("upper's indices must match g0, got $(indices(upper)) and $(indices(g0))"))
    Qidense = Matrix{T}[]
    Qisparse = SparseMatrixCSC{T,Int}[]
    gidense = Vector{T}[]
    gisparse = SparseVector{T,Int}[]
    cl, cu = T[], T[]
    issparseQi = Bool[]
    issparsegi = Bool[]
    isposdefQi = Bool[]
    return Problem{T,typeof(Q0),typeof(g0)}(n, lower, upper, cl, cu,
        Q0, g0, Qidense, Qisparse, gidense, gisparse,
        issparseQi, issparsegi, isposdefQi)
end

iszerom(Q::Matrix)          = false
iszerom(Q::SparseMatrixCSC) = isempty(Q.nzval)

zerom(::Type{T}, n) where T = sparse(Int[], Int[], T[], n, n)
zerov(::Type{T}, n) where T = sparse(Int[], T[], n)

function addconstraint!(P::Problem{T}, Qi::AbstractMatrix, gi::AbstractVector, cl::Real, cu::Real;
                        isposdef::Bool=false) where T
    n = P.n
    indices(gi) == (Base.OneTo(n),) || throw(DimensionMismatch("gi indices must match 1:$n, got $(indices(gi))"))
    indices(Qi) == (Base.OneTo(n), Base.OneTo(n))  || throw(DimensionMismatch("Qi indices must match (1:$n, 1:$n), got $(indices(Qi))"))
    Qsp = issparse(Qi)
    gsp = issparse(gi)
    Qsp ? push!(P.Qisparse, Qi) : push!(P.Qidense, Qi)
    gsp ? push!(P.gisparse, gi) : push!(P.gidense, gi)
    push!(P.cl, cl)
    push!(P.cu, cu)
    push!(P.issparseQi, Qsp)
    push!(P.issparsegi, gsp)
    push!(P.isposdefQi, isposdef)
    return P
end
function addconstraint!(P::Problem{T}, Qi::Real, gi::Real, cl::Real, cu::Real; kwargs...) where T
    @assert Qi == 0 && gi == 0
    n = P.n
    return addconstraint!(P, zerom(T, n), zerov(T, n), cl, cu; kwargs...)
end
function addconstraint!(P::Problem{T}, Qi::Real, gi::AbstractVector, cl::Real, cu::Real; kwargs...) where T
    @assert Qi == 0
    n = P.n
    return addconstraint!(P, zerom(T, n), gi, cl, cu; kwargs...)
end
function addconstraint!(P::Problem{T}, Qi::AbstractMatrix, gi::Real, cl::Real, cu::Real; kwargs...) where T
    @assert gi == 0
    n = P.n
    return addconstraint!(P, Qi, zerov(T, n), cl, cu; kwargs...)
end

function minimize(P::Problem; usefactorization::Bool=!issparse(P.Q0))
    if usefactorization
        return minimize_factorization(P)
    else
        error("not yet implemented")
    end
end

function minimize_factorization(P::Problem{T}, x::AbstractVector, μ) where T
    m, n = length(P.issparseQi), P.n
    indices(x) == (Base.OneTo(n),) || throw(DimensionMismatch("x must have indices 1:$n, got $(linearindices(x))"))
    λx, λc, sx, sc, constrval = initialize_auxvars(P, x)
    error("unfinished")
    # if isempty(λx) && m == 1
    #     # No bounds constraints
end

function initialize_auxvars(P::Problem{T}, x) where T
    m, n = length(P.issparseQi), P.n
    # How many inequality constraints are "real"? (x[k] <= Inf does not count as "real")
    nx = count(isfinite, P.l)  + count(isfinite, P.u)
    nc = count(isfinite, P.cl) + count(isfinite, P.cu)
    # Lagrange multipliers
    λx = zeros(T, nx)
    λc = zeros(T, nc)
    # Slack variables
    sx = zeros(T, nx)
    sc = zeros(T, nc)
    # Temporaries
    constrval = fill(T(NaN), m)
    # Initialize the slacks
    si = 0
    for k = 1:n
        if isfinite(P.l[k])
            sx[si+=1] = min(zero(T), x[k] - P.l[k])
        end
    end
    for k = 1:n
        if isfinite(P.u[k])
            sx[si+=1] = min(zero(T), P.u[k] - x[k])
        end
    end
    si = 0
    for k = 1:m
        if isfinite(P.cl[k])
            sc[si+=1] = min(zero(T), evalconstr!(constrval, P, k, x) - P.cl[k])
        end
    end
    for k = 1:m
        if isfinite(P.cu[k])
            sc[si+=1] = min(zero(T), P.cu[k] - evalconstr!(constrval, P, k, x))
        end
    end
    return λx, λc, sx, sc, constrval
end

end # module


text:
using JuLIP
using JuLIP.Testing
using JuLIP: sigvol_d
using Test
using LinearAlgebra

##

h2("Testing `minimise!` with equilibration with LJ calculator to lattice")
calc = lennardjones(r0=rnn(:Al))
at = bulk(:Al, cubic=true) * 10
X0 = positions(at) |> mat
at = rattle!(at, 0.02)
set_calculator!(at, calc)
x = dofs(at)
println(@test (energy(at) == energy(at, x) == energy(calc, at)
                          == energy(calc, at, x) )
               )
minimise!(at, precond=:id, verbose=2)
X1 = positions(at) |> mat
X0 .-= X0[:, 1]
X1 .-= X1[:, 1]
F = X1 / X0
println("check that the optimiser really converged to a lattice")
@show norm(F'*F - I, Inf)
@show norm(F*X0 - X1, Inf)
@test norm(F*X0 - X1, Inf) < 1e-4

##

h2("same test but large and with Exp preconditioner")
at = bulk(:Al, cubic=true) * (20,20,2)
at = rattle!(at, 0.02)
set_calculator!(at, calc)
minimise!(at, precond = :exp, method = :lbfgs,
          robust_energy_difference = true, verbose=2
          )

##

h2("Variable Cell Test")
calc = lennardjones(r0=rnn(:Al))
at = set_pbc!(bulk(:Al, cubic=true), true)
set_calculator!(at, calc)
variablecell!(at)
x = dofs(at)
println(@test (energy(at) == energy(at, x) == energy(calc, at)
                          == energy(calc, at, x) )
               )
minimise!(at, verbose = 2)

##

h2("FF preconditioner for StillingerWeber")
at = bulk(:Si, cubic=true) * (10,10,2)
at = set_pbc!(at, true)
at = rattle!(at, 0.02)
set_calculator!(at, StillingerWeber())
P = FF(at, StillingerWeber())
minimise!(at, precond = P, method = :lbfgs,
          robust_energy_difference = true, verbose=2)

##

h2("FF preconditioner for EAM")
at = bulk(:W, cubic=true) * (10,10,2)
at = set_pbc!(at, true)
at = rattle!(at, 0.02)
X0 = positions(at)
set_positions!(at, X0)
set_calculator!(at, eam_W)
P = FF(at, eam_W)
minimise!(at, precond = P, method = :lbfgs, robust_energy_difference = true, verbose=2)

## steepest descent
set_positions!(at, X0)
set_calculator!(at, eam_W)
P = FF(at, eam_W)
minimise!(at, precond = P, method = :sd, robust_energy_difference = true, verbose=2)

##
h2("Optimise again with some different stabilisation options")
set_positions!(at, X0)
set_calculator!(at, eam_W)
P = FF(at, eam_W, stab=0.1, innerstab=0.2)
minimise!(at, precond = P, method = :lbfgs, robust_energy_difference = true, verbose=2)

##
h2("for comparison now with Exp")
set_positions!(at, X0)
minimise!(at, precond = :exp, method = :lbfgs, robust_energy_difference = true, verbose=2)

h2("Test optimisation with VariableCell")
# start with a clean `at`
at = bulk(:Al) * 2   # cubic=true,
apply_defm!(at, I + 0.02 * rand(3,3))
calc = lennardjones(r0=rnn(:Al))
set_calculator!(at, calc)
variablecell!(at)
println(@test JuLIP.Testing.fdtest(calc, at, verbose=true, rattle=0.1))

h2("For the initial state, stress/virial is far from 0:")
@show norm(virial(at), Inf)
JuLIP.Solve.minimise!(at, verbose=2)
println("After optimisation, stress/virial should be 0:")
@show norm(virial(at), Inf)
@test norm(virial(at), Inf) < 1e-4

h2("Check sigvol derivative")
println(@test fdtest(c -> JuLIP.sigvol(reshape(c, (3,3))),
                     c -> JuLIP.sigvol_d(reshape(c, (3,3)))[:],
                     rand(3,3)))

# TODO: revive this after introducing external potentials
# h2("And now with pressure . . .")
# set_constraint!(at, VariableCell(at, pressure=10.0123))
# JuLIP.Testing.fdtest(calc, at, verbose=true, rattle=0.02)
# at = bulk(:Al) * 2
# set_calculator!(at, calc)
# set_constraint!(at, VariableCell(at, pressure=0.01))
# JuLIP.Solve.minimise!(at, verbose = 2)
# @show norm(virial(at), Inf)
# @show norm(JuLIP.gradient(at), Inf)
# println(@test norm(JuLIP.gradient(at), Inf) < 1e-4)
# @info "note it is correct that virial is O(1) since we applied pressure"


text:
abstract type Recurrence end
abstract type LinearRecurrence <: Recurrence end

struct HyperRecurrence{T} <: LinearRecurrence
    func::Symbol
    arg::PolyElem{T}
    coeffs::Vector{<:FracElem{<:PolyElem{T}}}
end

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

coeffs(r::HyperRecurrence) = r.coeffs
order(r::HyperRecurrence) = length(r.coeffs) - 1

function Base.show(io::IO, r::HyperRecurrence)
    res = join(["$(c) * $(r.func)($(r.arg + (i - 1)))" for (i, c) in enumerate(r.coeffs)], " + ")
    print(io, "$(res) = 0")
end

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

function closedform(r::HyperRecurrence{T}; init=nothing) where {T}
    sz = order(r)
    vars = [initvar(string(r.func), i) for i in 0:sz - 1]
    R, b = PolynomialRing(base_ring(r.arg), vars)
    F = FractionField(R)

    sequences = petkovsek(r.coeffs, r.arg)
    @debug "After Petkovsek" sequences init

    A = [F(s(i)) for i in 0:sz-1, s in sequences]
    b = [F(x) for x in b]
    mA = MatrixSpace(F, size(A)...)(A)
    mB = MatrixSpace(F, length(b), 1)(b)
    mX = Nemo.solve(mA, mB)

    if isnothing(init)
        cfs = Array(mX)
    else
        ini = Array(init)
        cfs = []
        for x in Array(mX)
            n, d = numerator(x), denominator(x)
            n = isone(n) ? one(first(ini)) : sum(i*c for (i, c) in zip(ini, Nemo.coeffs(n)))
            d = isone(d) ? one(first(ini)) : sum(i*c for (i, c) in zip(ini, Nemo.coeffs(d)))
            push!(cfs, n//d)
        end
    end

    F = isnothing(init) ? F : base_ring(base_ring(parent(first(cfs))))
    sequences = [change_coeff_field(F, s) for s in sequences]
    S = parent(first(sequences))
    @debug "Coefficient field" F S

    cform = sum(S(c)*s for (c, s) in zip(cfs, sequences))

    @debug "Closed form" sequences mA mB mX cform
    return cform
end


text:
```Write a function on_all that applies a function to every element of a list.
Use it to print the first twenty perfect squares. The perfect squares can be
found by multiplying each natural number with itself. The first few perfect
squares are 1*1= 1, 2*2=4, 3*3=9, 4*4=16. Twelve for example is not a perfect
square because there is no natural number m so that m*m=12. (This question is
tricky if your programming language makes it difficult to pass functions as arguments.)```

function on_all(x::Int64)
   for i in x
      if isinteger(sqrt(i)) == true
         return i
      end
   end
end

function find_PS()
   c = 0
   cou = 1
   while c <= 200
      if isa(on_all(cou), Int64)
         println( on_all(cou))
         c += 1
      end
   cou += 1
   end
end

@time find_PS()

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
using HawkesSimulator ; global const H = HawkesSimulator
using Test
using Distributions,Statistics,StatsBase,LinearAlgebra
using QuadGK
using Random ; Random.seed!(0)

@testset "kernels" begin
  # integral is 1
  mykernels = let tau = rand(Uniform(0.3,0.8))
    k1 = H.KernelExp(tau)
    k2 = H.KernelStep(15.0)
    tau = rand(Uniform(0.3,0.8))
    taud = rand(Uniform(0.1,0.4))
    k3 = H.KernelAlphaDelay(tau,taud)
    (k1,k2,k3)
  end

  for ker in mykernels
    integr =  quadgk(t->H.interaction_kernel(t,ker),-3.0,20.0)[1]
    @test isapprox(integr,1;atol=0.001)
  end

  # to do - Fourier!!!

end

@testset "Flush spiketrain" begin
  hnew ,spktnew = H._flush_train!([1.,2.,3.],[ 4., 5., 12. , 13. ],3.0+3.0)
  @test all(hnew .== [1.,2,3,4,5.])
  @test all(spktnew .== [12.,13.])
end

@testset "2D Vs linear model" begin
  weights = [0.31 -0.3
            0.9  -0.15]
  inputs = [5.0 , 5.0]
  rates_analytic  = inv(I-weights)*inputs

  tau =  0.5
  mykernel = H.KernelExp(tau) #considering single population
  nneus = length(inputs)
  popstate = H.PopulationState(mykernel,nneus)
  network = H.RecurrentNetwork(popstate,weights,inputs)

  function simulate!(network,num_spikes)
    t_now = 0.0
    H.reset!(network) # clear spike trains etc
    for _ in 1:num_spikes
      t_now = H.dynamics_step!(t_now,network)
      H.flush_trains!(popstate,10.0;Tflush=5.0)
    end
    H.flush_trains!(popstate)
    return t_now
  end
  n_spikes = 200_000
  Tmax = simulate!(network,n_spikes)
  rates_ntw = H.numerical_rates(network.populations[1])
  @test all(isapprox.(rates_analytic,rates_ntw;rtol=0.2))

  # do it again, checking the one-population optimized function
  function simulate_onepop!(network,num_spikes)
    t_now = 0.0
    H.reset!(network) # clear spike trains etc
    for _ in 1:num_spikes
      t_now = H.dynamics_step_singlepopulation!(t_now,network)
      H.flush_trains!(popstate,10.0;Tflush=5.0)
    end
    H.flush_trains!(popstate)
    return t_now
  end
  n_spikes = 200_000
  Tmax = simulate!(network,n_spikes)
  rates_ntw = H.numerical_rates(network.populations[1])
  @test all(isapprox.(rates_analytic,rates_ntw;rtol=0.2))

end

@testset "Warmup phase" begin
  N1 = 33
  N2 = 12
  weights1 = fill(Inf,N1,N1)
  inputs1 = fill(Inf,N1)
  weights2 = fill(Inf,N2,N2)
  inputs2 = fill(Inf,N2)

  tau =  123.5
  mykernel1 = H.KernelExp(tau) #considering single population
  mykernel2 = H.KernelExp(rand()*tau) #considering single population

  popstate1 = H.PopulationState(mykernel1,N1)
  popstate2 = H.PopulationState(mykernel2,N2)

  population1 = H.PopulationHawkes(popstate1,inputs1,
    (H.ConnectionWeights(weights1),popstate1) )
  population2 = H.PopulationHawkes(popstate2,inputs2,
    (H.ConnectionWeights(weights2),popstate2) )

  # multi-pop constructor
  network = H.RecurrentNetwork(population1,population2)

  ##

  target_rates1 = rand(Uniform(5.,12.),N1)
  target_rates2 = rand(Uniform(12.,44.),N2)
  target_rates = [target_rates1,target_rates2]

  Ttot = 500.0

  H.do_warmup!(Ttot,network,target_rates)
  H.flush_trains!(network)

  @test all( isapprox.(H.numerical_rates(population1),target_rates1;rtol=0.2) )
  @test all( isapprox.(H.numerical_rates(population2),target_rates2;rtol=0.2) )
end

@testset "Poisson input units" begin
  N1 = 13
  N2 = 17
  target_rates1 = rand(Uniform(5.,12.),N1)
  target_rates2 = rand(Uniform(12.,44.),N2)

  population1 = let spikegen = H.SGPoisson(target_rates1)
    H.PopulationInput(spikegen,N1)
  end
  population2 = let spikegen = H.SGPoisson(target_rates2)
    H.PopulationInput(spikegen,N2)
  end

  network = H.RecurrentNetwork(population1,population2)
  function simulate!(network,num_spikes)
    t_now = 0.0
    H.reset!(network) # clear spike trains etc
    for _ in 1:num_spikes
      t_now = H.dynamics_step!(t_now,network)
    end
    H.flush_trains!(network)
    return t_now
  end

  n_spikes = 200_000
  Tmax = simulate!(network,n_spikes)
  rates_ntw1 = H.numerical_rates(network.populations[1])
  rates_ntw2 = H.numerical_rates(network.populations[2])

  @test all(isapprox.(target_rates1,rates_ntw1;rtol=0.1))
  @test all(isapprox.(target_rates2,rates_ntw2;rtol=0.1))
end

@testset "Input units, time-varying function" begin
  Nunits = 200
  function rate_fun(t::Float64,::Integer)
    return 10.0 + 6.0*cos(t*(2π)*0.5)
  end
  function rate_fun_upper(::Float64,::Integer)
    return 10.0 + 6.0
  end
  #plot(t->rate_fun(t,0),range(0,10;length=300);leg=false)
  population1 = let spikegen = H.SGPoissonFunction(rate_fun,rate_fun_upper)
    H.PopulationInput(spikegen,Nunits)
  end
  network = H.RecurrentNetwork(population1)
  function simulate!(network,num_spikes)
    t_now = 0.0
    H.reset!(network) # clear spike trains etc
    for _ in 1:num_spikes
      t_now = H.dynamics_step_singlepopulation!(t_now,network)
    end
    H.flush_trains!(network)
    return t_now
  end
  n_spikes = 50_000
  Tmax = simulate!(network,n_spikes)
  # now, compute rates!
  tmid,rats = H.instantaneous_rates(collect(1:Nunits),0.1,population1;Tend=Tmax)
  @test all(isapprox.(rats,rate_fun.(tmid,1);rtol=0.33))
end

@testset "Input units, given spike train" begin
  # post-pre spikes paradigm
  # (nor really needed here, but useful for plasticity)
  function ps_trains(rate::R,Δt_ro::R,Ttot::R;
      tstart::R = 0.05) where R
    post = collect(range(tstart,Ttot; step=inv(rate)))
    pre = post .- Δt_ro
    return [pre,post]
  end
  function get_test_pop(rate,nreps,Δt_ro,connection_test)
    Ttot = nreps/rate
    prepostspikes = ps_trains(rate,Δt_ro,Ttot)
    gen = H.SGTrains(prepostspikes)
    state = H.PopulationState(H.InputUnit(gen),2)
    return H.PopulationInputTestWeights(state,connection_test)
  end
  connection_test = H.ConnectionWeights(fill(0.0,2,2))
  nreps = 100
  population = get_test_pop(0.5,nreps,2E-3,connection_test)
  population.state.unittype.spike_generator.trains
  network = H.RecurrentNetwork(population)
  function simulate!(network,num_spikes)
    t_now = 0.0
    H.reset!(network) # clear spike trains etc
    for _ in 1:num_spikes
      t_now = H.dynamics_step_singlepopulation!(t_now,network)
    end
    H.flush_trains!(network)
    return t_now
  end
  n_spikes = nreps*2-3
  Tmax = simulate!(network,n_spikes)
  trains_out1 =  population.state.trains_history[1]
  target_trains1 = population.state.unittype.spike_generator.trains[1]
  trains_out2 =  population.state.trains_history[2]
  target_trains2 = population.state.unittype.spike_generator.trains[2]
  ntest = nreps-10
  @test all(trains_out1[1:ntest] .== target_trains1[1:ntest] )
  @test all(trains_out2[1:ntest] .== target_trains2[1:ntest] )
end

@testset "special case of exponential kernel" begin
  nneus = 1
  tauker = 12.34E-2
  trace_ker = H.Trace(tauker,nneus,H.ForDynamics())
  H.update_for_dynamics!(trace_ker,1)
  get_trace(t) = H.trace_proposal!([NaN,],t,trace_ker)[1]
  _area = quadgk(get_trace,0,2.0)[1]
  @test isapprox(_area,1.0;atol=1E-5)

  # like 2D test above, but for ExpKernel types

  nneus = 2
  tauker = 0.5
  trace_ker = H.Trace(tauker,nneus,H.ForDynamics())
  trace_useless = H.Trace(123.0,nneus,H.ForPlasticity())

  popstate = H.PopulationStateExpKernel(nneus,trace_ker,trace_useless)
  myweights = [0.31 -0.3
              0.9  -0.15]
  myinputs = [5.0 , 5.0]
  rates_analytic  = inv(I-myweights)*myinputs

  connection = H.ConnectionExpKernel(myweights,trace_ker)
  population = H.PopulationExpKernel(popstate,connection,myinputs)

  n_spikes = 10_000
  recorder = H.RecFullTrain(n_spikes,1)
  network = H.RecurrentNetworkExpKernel(population,recorder)

  function simulate!(network,num_spikes)
    t_now = 0.0
    H.reset!(network) # clear spike trains etc
    for _ in 1:num_spikes
      t_now = H.dynamics_step!(t_now,network)
    end
    return t_now
  end
  function simulate2!(network,num_spikes)
    t_now = 0.0
    H.reset!(network) # clear spike trains etc
    for _ in 1:num_spikes
      t_now = H.dynamics_step_singlepopulation!(t_now,network)
    end
    return t_now
  end
  Tmax = simulate!(network,n_spikes)
  rates_ntw1 = H.numerical_rates(recorder,nneus,Tmax)
  rates_ntw2 = H.numerical_rates(recorder,nneus,Tmax)
  @test all(isapprox.(rates_analytic,rates_ntw1;rtol=0.2))
  @test all(isapprox.(rates_analytic,rates_ntw2;rtol=0.2))
end

##


text:
module ErrorfreeArithmetic

import Base: (+),(-),(*),(/),fma,sqrt

export HasFMAonChip,
       eftIncr, eftDecr, eftSqr, eftRecip,   # single
       accInv, accInvForSign, accSqrt, accSqrtForSign,
       eftAdd, eftSub, eftMul, accDiv,       # double
       eftAddGTE, eftAddLTE, eftSubGTE, accDivForSign,
       eftAddAs2, eftMulAs4, eftMulAs2,      # triple (also eftAdd, eftAddGTE, eftMul)
       eftFMA, eftFMA, eftFMAas2, eftFMSas2,
       eftAddAs3, eftAddGTEas3,              # quadruple (also eftAdd, eftAddGTE)
       eftIncr_inline, eftDecr_inline,       # inline versions
       eftSqr_inline,
       eftAdd_inline, eftAddAs2_inline, eftAddGTE_inline,
       eftSub_inline, eftSubGTE_inline,
       eftMul_inline, eftMulAs2_inline

# all functions use type StdFloat
if !isdefined(:StdFloat)
  if isdefined(Main,:StdFloat)
    import Main:StdFloat
  else
    typealias StdFloat Union{Float64,Float32,Float16}
  end
end

include("HasFMAonChip.jl")  # const HasFMAonChip::Bool

include("single.jl")
include("double.jl")
include("triple.jl")
include("quadruple.jl")

include("renormalize.jl")
include("inline.jl")
include("doubdoub.jl")

end # module


text:
# Suppress warnings during testing.
# required for reducing logging during tests
using Memento
setlevel!(getlogger(PowerModels), "error")

# is this the best way to find a file in a package?
base_dir = string(dirname(dirname(pathof(PowerSystems))))
case5_data = PM.parse_file(joinpath(base_dir, "data/matpower/case5.m"))
case5_data = PM.replicate(case5_data, 2)

case5_dc_data = PM.parse_file(joinpath(base_dir, "data/matpower/case5_dc.m"))
case5_dc_data = PM.replicate(case5_dc_data, 2)

# TODO: currently JuMP.num_variables is the best we can do to introspect the JuMP model.
#  Ideally this would also test the number of constraints generated

@testset "PowerModels Model Build" begin
    pm = PowerSimulations.build_nip_model(case5_data, DCPPowerModel)
    @test JuMP.num_variables(pm.model) == 34
    pm = PowerSimulations.build_nip_model(case5_data, PM.ACPPowerModel)
    @test JuMP.num_variables(pm.model) == 96
    pm = PowerSimulations.build_nip_model(case5_data, PM.SOCWRPowerModel)
    @test JuMP.num_variables(pm.model) == 110
end

@testset "PM with type extensions" begin
    pm = PowerSimulations.build_nip_model(case5_data, DCPPowerModel)
    JuMP.num_variables(pm.model) == 34
    pm = PowerSimulations.build_nip_model(case5_data, ACPPowerModel)
    JuMP.num_variables(pm.model) == 96
    pm = PowerSimulations.build_nip_model(case5_dc_data, DCPPowerModel)
    JuMP.num_variables(pm.model) == 36
    pm = PowerSimulations.build_nip_model(case5_dc_data, DCPPowerModel)
    JuMP.num_variables(pm.model) == 48
end
#=
@testset "PM integration into PS" begin
    PM_dict5 = PowerSimulations.pass_to_pm(c_sys5, 24)
    PM_object5 = PowerSimulations.build_nip_model(PM_dict5, DCAngleModel);
    @test JuMP.num_variables(PM_object5.model) == 384

    PM_dict_rts = PowerSimulations.pass_to_pm(c_rts, 24)
    PM_object_rts = PowerSimulations.build_nip_model(PM_dict_rts, DCAngleModel);
    @test JuMP.num_variables(PM_object_rts.model) == 6432
end
=#


text:
function is_jll_name(name::AbstractString)::Bool
    return endswith(name, "_jll")
end

function _get_all_dependencies_nonrecursive(working_directory::AbstractString, pkg, version)
    all_dependencies = String[]
    package_relpath = get_package_relpath_in_registry(;
        package_name=pkg, registry_path=working_directory
    )
    deps = Pkg.TOML.parsefile(joinpath(working_directory, package_relpath, "Deps.toml"))
    for version_range in keys(deps)
        if version in Pkg.Types.VersionRange(version_range)
            for name in keys(deps[version_range])
                push!(all_dependencies, name)
            end
        end
    end
    unique!(all_dependencies)
    return all_dependencies
end

const guideline_allowed_jll_nonrecursive_dependencies = Guideline(;
    info="If this is a JLL package, only deps are Pkg, Libdl, and other JLL packages",
    docs=nothing,
    check=data -> meets_allowed_jll_nonrecursive_dependencies(
        data.registry_head, data.pkg, data.version
    ),
)

function meets_allowed_jll_nonrecursive_dependencies(
    working_directory::AbstractString, pkg, version
)
    # If you are a JLL package, you are only allowed to have five kinds of dependencies:
    # 1. Pkg
    # 2. Libdl
    # 3. Artifacts
    # 4. JLLWrappers
    # 5. LazyArtifacts
    # 6. TOML
    # 8. MPIPreferences
    # 7. other JLL packages
    all_dependencies = _get_all_dependencies_nonrecursive(working_directory, pkg, version)
    allowed_dependencies = ("Pkg", "Libdl", "Artifacts", "JLLWrappers", "LazyArtifacts", "TOML", "MPIPreferences")
    for dep in all_dependencies
        if dep ∉ allowed_dependencies && !is_jll_name(dep)
            return false,
            "JLL packages are only allowed to depend on $(join(allowed_dependencies, ", ")) and other JLL packages"
        end
    end
    return true, ""
end


text:
module TestMachines

using MLJBase
using Test
using Statistics
using ..Models

@load DecisionTreeRegressor

N=50
X = (a=rand(N), b=rand(N), c=rand(N))
y = 2*X.a - X.c + 0.05*rand(N)

train, test = partition(eachindex(y), 0.7);

tree = DecisionTreeRegressor(max_depth=5)

t = Machine(tree, X, y)
@test_logs (:info, r"Training") fit!(t)
@test_logs (:info, r"Training") fit!(t, rows=train)
@test_logs (:info, r"Not retraining") fit!(t, rows=train)
@test_logs (:info, r"Training") fit!(t)
t.model.max_depth = 1
# MLJBase.recursive_setproperty!(t, :(model.max_depth),  1)
@test_logs (:info, r"Updating") fit!(t)

predict(t, selectrows(X,test))
@test rms(predict(t, selectrows(X, test)), y[test]) < std(y)

mach = machine(ConstantRegressor(), X, y)
@test_logs (:info, r"Training") fit!(mach)
yhat = predict_mean(mach, X)

n = nrows(X)
rms(yhat, y) ≈ std(y)*sqrt(1 - 1/n)

# test an unsupervised univariate case:
mach = machine(UnivariateStandardizer(), float.(1:5))
@test_logs (:info, r"Training") fit!(mach)
@test isempty(params(mach))

# test a frozen NodalMachine
stand = machine(Standardizer(), source((x1=rand(10),)))
freeze!(stand)
@test_logs (:warn, r"not trained as it is frozen\.$") fit!(stand)

@testset "warnings" begin
    @test_logs((:info, r"does not support"),
               machine(tree, X, y, rand(N)))
    @test_throws DimensionMismatch machine(tree, X, y[1:end-1])
    @test_logs((:info, r"does not support"),
               @test_throws DimensionMismatch machine(tree, X, y, rand(N-1)))
    @test_logs((:info, r"does not support"),
               @test_throws ArgumentError machine(tree, X, y, fill('a', N)))
    @test_logs((:warn, r"The scitype of `y`"),
               machine(tree, X, categorical(1:N)))
    @test_logs((:warn, r"The scitype of `X`"),
               machine(tree, (x=categorical(1:N),), y))
end

@testset "weights" begin
    yraw = ["Perry", "Antonia", "Perry", "Skater"]
    X = (x=rand(4),)
    y = categorical(yraw)
    w = [2, 3, 2, 5]

    # without weights:
    mach = machine(ConstantClassifier(), X, y)
    fit!(mach)
    d1 = predict(mach)[1]
    d2 = MLJBase.UnivariateFinite([y[1], y[2], y[4]], [0.5, 0.25, 0.25])
    @test all([pdf(d1, c) ≈ pdf(d2, c) for c in MLJBase.classes(d1)])

    # with weights:
    mach = machine(ConstantClassifier(), X, y, w)
    fit!(mach)
    d1 = predict(mach)[1]
    d2 = MLJBase.UnivariateFinite([y[1], y[2], y[4]], [1/3, 1/4, 5/12])
    @test all([pdf(d1, c) ≈ pdf(d2, c) for c in MLJBase.classes(d1)])
end

end # module
true


text:
export @vary

macro vary(args...)
    params = Dict{Symbol,Dict{Symbol,Float64}}()
    for ex in args[1:end-1]
        if !(ex.head in (:comparison, :(=), :kw)
             || ex.head == :call && ex.args[1] in (:(<=), :(>=)))
            throw(ArgumentError("Expression '$ex' is neither a comparison nor an assignment."))
        end
        if ex.head == :call
            op, sym, value = ex.args
            if op in (:(<=), :(<))
                key = :upper
            elseif op in (:(>=), :(>))
                key = :lower
            end
        elseif length(ex.args) == 3
            println(ex.head)
            println(ex.args)
            sym, op, value = ex.args
            if op in (:(<=), :(<))
                key = :upper
            elseif op in (:(>=), :(>))
                key = :lower
            end
        else
            sym, value = ex.args
            key = :initial
        end
        if !isa(value, Number)
            v = eval(value)
        else
            v = value
        end
        d = Dict(key=>v)
        if !(sym in keys(params))
            merge!(params, Dict(sym=>d))
        else
            merge!(params[sym], d)
        end
    end
    typ = args[end]
    if typ.head != :call
        throw(ArgumentError("The last expression must be a type instantiation."))
    end
    replace_parameters!(typ, params)
    return :($(esc(typ)))
end

function replace_parameters!(expr, params)
    for (i, ex) in enumerate(expr.args)
        if isa(ex, Expr)
            replace_parameters!(expr.args[i], params)
        elseif isa(ex, Symbol) && ex in keys(params)
            if expr.head != :kw || (expr.head == :kw && i > 1)
                upper = get(params[ex], :upper, Inf)
                lower = get(params[ex], :lower, -Inf)
                initial = get(params[ex], :initial, (upper-lower)/2)
                expr.args[i] = Parameter(initial, lower, upper)
            end
        end
    end
end


text:
"""
    get_pis(sequences)
"""

function get_pis(sequences)
    c= countmap((collect(string(sequences...))))
    return (c['A'], c['C'], c['G'], c['T']) ./ sum([c['A'], c['C'], c['G'], c['T']]);
end

"""
    get_transition_mat(seq1, seq2; include_dash_seqs = true)
"""
function get_transition_mat(seq1, seq2; include_dash_seqs = true)
    a1, a2 = kmer_seeded_align(seq1, seq2)

    indices = ['A', 'C', 'G', 'T']
    matrix = zeros(4,4)

    for i in 1:length(a1)
        if a1[i] == '-' || a2[i] == '-'
            continue
        else
            matrix[findfirst(indices, a1[i]), findfirst(indices, a2[i])] += 1
        end
    end
    return matrix
end

"""
    estimate_distance(seq1, seq2)
"""
function estimate_distance(seq1, seq2)
    matrix = get_transition_mat(seq1, seq2)
    pi_a, pi_c, pi_g, pi_t = get_pis([seq1, seq2])
    pi_y = pi_t + pi_c
    pi_r = pi_a + pi_g

    s_1 = (matrix[2,4] + matrix[4,2])/sum(matrix)
    s_2 = (matrix[1,3] + matrix[3,1])/sum(matrix)
    v = (matrix[1,2] + matrix[2,1] + matrix[1,4] + matrix[4,1] + matrix[2,3] + matrix[3,2] + matrix[3,4] + matrix[4,3])/sum(matrix)

    a_1 = - log(1 - (pi_y * s_1)/(2 * pi_t * pi_c) - (v/(2 * pi_y) ))
    a_2 = - log(1 - (pi_r * s_2)/(2 * pi_a * pi_g) - (v/(2 * pi_r) ))
    b = - log(1 - (v/(2 * pi_y * pi_r)))

    d = ((2 * pi_t * pi_c)/pi_y)*(a_1 - (pi_r * b)) + ((2 * pi_a * pi_g)/pi_r)*(a_2 - (pi_y * b)) + (2 * pi_y * pi_r * b)

    return d
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Random

@testset "constructors" begin
    v = [0x61,0x62,0x63,0x21]
    @test String(v) == "abc!" && isempty(v)
    @test String("abc!") == "abc!"
    @test String(0x61:0x63) == "abc"

    # Check that resizing empty source vector does not corrupt string
    b = IOBuffer()
    write(b, "ab")
    x = take!(b)
    s = String(x)
    resize!(x, 0)
    empty!(x) # Another method which must be tested
    @test s == "ab"
    resize!(x, 1)
    @test s == "ab"

    @test isempty(string())
    @test eltype(GenericString) == Char
    @test firstindex("abc") == 1
    @test cmp("ab","abc") == -1
    @test typemin(String) === ""
    @test typemin("abc") === ""
    @test "abc" === "abc"
    @test "ab"  !== "abc"
    @test string("ab", 'c') === "abc"
    @test string() === ""
    codegen_egal_of_strings(x, y) = (x===y, x!==y)
    @test codegen_egal_of_strings(string("ab", 'c'), "abc") === (true, false)
    let strs = ["", "a", "a b c", "до свидания"]
        for x in strs, y in strs
            @test (x === y) == (objectid(x) == objectid(y))
        end
    end
end

@testset "{starts,ends}with" begin
    @test startswith("abcd", 'a')
    @test startswith("abcd", "a")
    @test startswith("abcd", "ab")
    @test !startswith("ab", "abcd")
    @test !startswith("abcd", "bc")
    @test endswith("abcd", 'd')
    @test endswith("abcd", "d")
    @test endswith("abcd", "cd")
    @test !endswith("abcd", "dc")
    @test !endswith("cd", "abcd")
    @test startswith("ab\0cd", "ab\0c")
    @test !startswith("ab\0cd", "ab\0d")
end

@test filter(x -> x ∈ ['f', 'o'], "foobar") == "foo"

@testset "string iteration, and issue #1454" begin
    str = "é"
    str_a = vcat(str...)
    @test length(str_a)==1
    @test str_a[1] == str[1]

    str = "s\u2200"
    @test str[1:end] == str
end

@testset "sizeof" begin
    @test sizeof("abc") == 3
    @test sizeof("\u2222") == 3
end

# issue #3597
@test string(GenericString("Test")[1:1], "X") == "TX"

@testset "parsing Int types" begin
    let b, n
    for T = (UInt8,Int8,UInt16,Int16,UInt32,Int32,UInt64,Int64,UInt128,Int128,BigInt),
            b = 2:62,
            _ = 1:10
        n = (T != BigInt) ? rand(T) : BigInt(rand(Int128))
        @test parse(T, string(n, base = b),  base = b) == n
    end
    end
end

@testset "issue #6027 - make symbol with invalid char" begin
    sym = Symbol(Char(0xdcdb))
    @test string(sym) == string(Char(0xdcdb))
    @test String(sym) == string(Char(0xdcdb))
    @test Meta.lower(Main, sym) === sym
    res = string(Meta.parse(string(Char(0xdcdb)," = 1"),1,raise=false)[1])
    @test res == """\$(Expr(:error, "invalid character \\\"\\udcdb\\\"\"))"""
end

@testset "Symbol and gensym" begin
    @test Symbol("asdf") === :asdf
    @test Symbol(:abc,"def",'g',"hi",0) === :abcdefghi0
    @test :a < :b
    @test startswith(string(gensym("asdf")),"##asdf#")
    @test gensym("asdf") != gensym("asdf")
    @test gensym() != gensym()
    @test startswith(string(gensym()),"##")
    @test_throws ArgumentError Symbol("ab\0")
    @test_throws ArgumentError gensym("ab\0")
end
@testset "issue #6949" begin
    f = IOBuffer()
    x = split("1 2 3")
    local nb = 0
    for c in x
        nb += write(f, c)
    end
    @test nb == 3
    @test String(take!(f)) == "123"
end

@testset "issue #7248" begin
    @test_throws BoundsError length("hello", 1, -1)
    @test_throws BoundsError prevind("hello", 0, 1)
    @test_throws BoundsError length("hellø", 1, -1)
    @test_throws BoundsError prevind("hellø", 0, 1)
    @test_throws BoundsError length("hello", 1, 10)
    @test nextind("hello", 0, 10) == 10
    @test_throws BoundsError length("hellø", 1, 10) == 9
    @test nextind("hellø", 0, 10) == 11
    @test_throws BoundsError checkbounds("hello", 0)
    @test_throws BoundsError checkbounds("hello", 6)
    @test_throws BoundsError checkbounds("hello", 0:3)
    @test_throws BoundsError checkbounds("hello", 4:6)
    @test_throws BoundsError checkbounds("hello", [0:3;])
    @test_throws BoundsError checkbounds("hello", [4:6;])
    @test checkbounds("hello", 1) === nothing
    @test checkbounds("hello", 5) === nothing
    @test checkbounds("hello", 1:3) === nothing
    @test checkbounds("hello", 3:5) === nothing
    @test checkbounds("hello", [1:3;]) === nothing
    @test checkbounds("hello", [3:5;]) === nothing
    @test checkbounds(Bool, "hello", 0) === false
    @test checkbounds(Bool, "hello", 1) === true
    @test checkbounds(Bool, "hello", 5) === true
    @test checkbounds(Bool, "hello", 6) === false
    @test checkbounds(Bool, "hello", 0:5) === false
    @test checkbounds(Bool, "hello", 1:6) === false
    @test checkbounds(Bool, "hello", 1:5) === true
    @test checkbounds(Bool, "hello", [0:5;]) === false
    @test checkbounds(Bool, "hello", [1:6;]) === false
    @test checkbounds(Bool, "hello", [1:5;]) === true
end

@testset "issue #15624 (indexing with out of bounds empty range)" begin
    @test ""[10:9] == ""
    @test "hello"[10:9] == ""
    @test "hellø"[10:9] == ""
    @test SubString("hello", 1, 5)[10:9] == ""
    @test SubString("hello", 1, 0)[10:9] == ""
    @test SubString("hellø", 1, 5)[10:9] == ""
    @test SubString("hellø", 1, 0)[10:9] == ""
    @test SubString("", 1, 0)[10:9] == ""
    @test_throws BoundsError SubString("", 1, 6)
    @test_throws BoundsError SubString("", 1, 1)
end

@testset "issue #22500 (using `get()` to index strings with default returns)" begin
    utf8_str = "我很喜欢Julia"

    # Test that we can index in at valid locations
    @test get(utf8_str, 1, 'X') == '我'
    @test get(utf8_str, 13, 'X') == 'J'

    # Test that obviously incorrect locations return the default
    @test get(utf8_str, -1, 'X') == 'X'
    @test get(utf8_str, 1000, 'X') == 'X'

    # Test that indexing into the middle of a character throws
    @test_throws StringIndexError get(utf8_str, 2, 'X')
end

#=
# issue #7764
let
    srep = repeat("Σβ",2)
    s="Σβ"
    ss=SubString(s,1,lastindex(s))

    @test repeat(ss,2) == "ΣβΣβ"

    @test lastindex(srep) == 7

    @test iterate(srep, 3) == ('β',5)
    @test iterate(srep, 7) == ('β',9)

    @test srep[7] == 'β'
    @test_throws BoundsError srep[8]
end
=#

# This caused JuliaLang/JSON.jl#82
@test first('\x00':'\x7f') === '\x00'
@test last('\x00':'\x7f') === '\x7f'

# make sure substrings do not accept code unit if it is not start of codepoint
let s = "x\u0302"
    @test s[1:2] == s
    @test_throws BoundsError s[0:3]
    @test_throws BoundsError s[1:4]
    @test_throws StringIndexError s[1:3]
end

@testset "issue #9781" begin
    # parse(Float64, SubString) wasn't tolerant of trailing whitespace, which was different
    # to "normal" strings. This also checks we aren't being too tolerant and allowing
    # any arbitrary trailing characters.
    @test parse(Float64,"1\n") == 1.0
    @test [parse(Float64,x) for x in split("0,1\n",",")][2] == 1.0
    @test_throws ArgumentError parse(Float64,split("0,1 X\n",",")[2])
    @test parse(Float32,"1\n") == 1.0
    @test [parse(Float32,x) for x in split("0,1\n",",")][2] == 1.0
    @test_throws ArgumentError parse(Float32,split("0,1 X\n",",")[2])
end
# test AbstractString functions at beginning of string.jl
struct tstStringType <: AbstractString
    data::Array{UInt8,1}
end
@testset "AbstractString functions" begin
    tstr = tstStringType(unsafe_wrap(Vector{UInt8},"12"))
    @test_throws MethodError ncodeunits(tstr)
    @test_throws MethodError codeunit(tstr)
    @test_throws MethodError codeunit(tstr, 1)
    @test_throws MethodError codeunit(tstr, true)
    @test_throws MethodError isvalid(tstr, 1)
    @test_throws MethodError isvalid(tstr, true)
    @test_throws MethodError iterate(tstr, 1)
    @test_throws MethodError iterate(tstr, true)
    @test_throws MethodError lastindex(tstr)

    gstr = GenericString("12")
    @test string(gstr) isa GenericString

    @test Array{UInt8}(gstr) == [49, 50]
    @test Array{Char,1}(gstr) == ['1', '2']

    @test gstr[1] == '1'
    @test gstr[1:1] == "1"
    @test gstr[[1]] == "1"

    @test s"∀∃"[big(1)] == '∀'
    @test_throws StringIndexError GenericString("∀∃")[Int8(2)]
    @test_throws BoundsError GenericString("∀∃")[UInt16(10)]

    @test first(eachindex("foobar")) === 1
    @test first(eachindex("")) === 1
    @test last(eachindex("foobar")) === lastindex("foobar")
    @test iterate(eachindex("foobar"),7) === nothing
    @test Int == eltype(Base.EachStringIndex) ==
                 eltype(Base.EachStringIndex{String}) ==
                 eltype(Base.EachStringIndex{GenericString}) ==
                 eltype(eachindex("foobar")) == eltype(eachindex(gstr))
    @test map(uppercase, "foó") == "FOÓ"
    @test nextind("fóobar", 0, 3) == 4

    @test Symbol(gstr) == Symbol("12")

    @test sizeof(gstr) == 2
    @test ncodeunits(gstr) == 2
    @test length(gstr) == 2
    @test length(GenericString("")) == 0

    @test nextind(1:1, 1) == 2
    @test nextind([1], 1) == 2

    @test length(gstr, 1, 2) == 2

    # no string promotion
    let svec = [s"12", GenericString("12"), SubString("123", 1, 2)]
        @test all(x -> x == "12", svec)
        @test svec isa Vector{AbstractString}
    end
end

@testset "issue #10307" begin
    @test typeof(map(x -> parse(Int16, x), AbstractString[])) == Vector{Int16}

    for T in [Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128]
        for i in [typemax(T), typemin(T)]
            s = "$i"
            @test tryparse(T, s) == i
        end
    end

    for T in [Int8, Int16, Int32, Int64, Int128]
        for i in [typemax(T), typemin(T)]
            f = "$(i)0"
            @test tryparse(T, f) === nothing
        end
    end
end

@testset "issue #11142" begin
    s = "abcdefghij"
    sp = pointer(s)
    @test unsafe_string(sp) == s
    @test unsafe_string(sp,5) == "abcde"
    @test typeof(unsafe_string(sp)) == String
    s = "abcde\uff\u2000\U1f596"
    sp = pointer(s)
    @test unsafe_string(sp) == s
    @test unsafe_string(sp,5) == "abcde"
    @test typeof(unsafe_string(sp)) == String

    @test tryparse(BigInt, "1234567890") == BigInt(1234567890)
    @test tryparse(BigInt, "1234567890-") === nothing

    @test tryparse(Float64, "64") == 64.0
    @test tryparse(Float64, "64o") === nothing
    @test tryparse(Float32, "32") == 32.0f0
    @test tryparse(Float32, "32o") === nothing
end

import Unicode

@testset "issue #10994: handle embedded NUL chars for string parsing" begin
    for T in [BigInt, Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128]
        @test_throws ArgumentError parse(T, "1\0")
    end
    for T in [BigInt, Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128, Float64, Float32]
        @test tryparse(T, "1\0") === nothing
    end
    let s = Unicode.normalize("tést",:NFKC)
        @test unsafe_string(Base.unsafe_convert(Cstring, Base.cconvert(Cstring, s))) == s
        @test unsafe_string(Base.unsafe_convert(Cstring, Symbol(s))) == s
    end
    @test_throws ArgumentError Base.unsafe_convert(Cstring, Base.cconvert(Cstring, "ba\0d"))

    cstrdup(s) = @static Sys.iswindows() ? ccall(:_strdup, Cstring, (Cstring,), s) : ccall(:strdup, Cstring, (Cstring,), s)
    let p = cstrdup("hello")
        @test unsafe_string(p) == "hello"
        Libc.free(p)
    end
end

@testset "iteration" begin
    @test [c for c in "ḟøøƀäṙ"] == ['ḟ', 'ø', 'ø', 'ƀ', 'ä', 'ṙ']
    @test [i for i in eachindex("ḟøøƀäṙ")] == [1, 4, 6, 8, 10, 12]
    @test [x for x in enumerate("ḟøøƀäṙ")] == [(1, 'ḟ'), (2, 'ø'), (3, 'ø'), (4, 'ƀ'), (5, 'ä'), (6, 'ṙ')]
end
@testset "isvalid edge conditions" begin
    for (val, pass) in (
            (0, true), (0xd7ff, true),
            (0xd800, false), (0xdfff, false),
            (0xe000, true), (0xffff, true),
            (0x10000, true), (0x10ffff, true),
            (0x110000, false)
        )
        @test isvalid(Char, val) == pass
    end
    for (val, pass) in (
            ("\x00", true),
            ("\x7f", true),
            ("\x80", false),
            ("\xbf", false),
            ("\xc0", false),
            ("\xff", false),
            ("\xc0\x80", false),
            ("\xc1\x80", false),
            ("\xc2\x80", true),
            ("\xc2\xc0", false),
            ("\xed\x9f\xbf", true),
            ("\xed\xa0\x80", false),
            ("\xed\xbf\xbf", false),
            ("\xee\x80\x80", true),
            ("\xef\xbf\xbf", true),
            ("\xf0\x90\x80\x80", true),
            ("\xf4\x8f\xbf\xbf", true),
            ("\xf4\x90\x80\x80", false),
            ("\xf5\x80\x80\x80", false),
            ("\ud800\udc00", false),
            ("\udbff\udfff", false),
            ("\ud800\u0100", false),
            ("\udc00\u0100", false),
            ("\udc00\ud800", false),
        )
        @test isvalid(String, val) == pass == isvalid(String(val))
        @test isvalid(val[1]) == pass
    end

    # Issue #11203
    @test isvalid(String, UInt8[]) == true == isvalid("")

    # Check UTF-8 characters
    # Check ASCII range (true),
    # then single continuation bytes and lead bytes with no following continuation bytes (false)
    for (rng,flg) in ((0:0x7f, true), (0x80:0xff, false))
        for byt in rng
            @test isvalid(String, UInt8[byt]) == flg
        end
    end
    # Check overlong lead bytes for 2-character sequences (false)
    for byt = 0xc0:0xc1
        @test isvalid(String, UInt8[byt,0x80]) == false
    end
    # Check valid lead-in to two-byte sequences (true)
    for byt = 0xc2:0xdf
        for (rng,flg) in ((0x00:0x7f, false), (0x80:0xbf, true), (0xc0:0xff, false))
            for cont in rng
                @test isvalid(String, UInt8[byt, cont]) == flg
            end
        end
    end
    # Check three-byte sequences
    for r1 in (0xe0:0xec, 0xee:0xef)
        for byt = r1
            # Check for short sequence
            @test isvalid(String, UInt8[byt]) == false
            for (rng,flg) in ((0x00:0x7f, false), (0x80:0xbf, true), (0xc0:0xff, false))
                for cont in rng
                    @test isvalid(String, UInt8[byt, cont]) == false
                    @test isvalid(String, UInt8[byt, cont, 0x80]) == flg
                end
            end
        end
    end
    # Check hangul characters (0xd000-0xd7ff) hangul
    # Check for short sequence, or start of surrogate pair
    for (rng,flg) in ((0x00:0x7f, false), (0x80:0x9f, true), (0xa0:0xff, false))
        for cont in rng
            @test isvalid(String, UInt8[0xed, cont]) == false
            @test isvalid(String, UInt8[0xed, cont, 0x80]) == flg
        end
    end
    # Check valid four-byte sequences
    for byt = 0xf0:0xf4
        if (byt == 0xf0)
            r0 = ((0x00:0x8f, false), (0x90:0xbf, true), (0xc0:0xff, false))
        elseif byt == 0xf4
            r0 = ((0x00:0x7f, false), (0x80:0x8f, true), (0x90:0xff, false))
        else
            r0 = ((0x00:0x7f, false), (0x80:0xbf, true), (0xc0:0xff, false))
        end
        for (rng,flg) in r0
            for cont in rng
                @test isvalid(String, UInt8[byt, cont]) == false
                @test isvalid(String, UInt8[byt, cont, 0x80]) == false
                @test isvalid(String, UInt8[byt, cont, 0x80, 0x80]) == flg
            end
        end
    end
    # Check five-byte sequences, should be invalid
    for byt = 0xf8:0xfb
        @test isvalid(String, UInt8[byt, 0x80, 0x80, 0x80, 0x80]) == false
    end
    # Check six-byte sequences, should be invalid
    for byt = 0xfc:0xfd
        @test isvalid(String, UInt8[byt, 0x80, 0x80, 0x80, 0x80, 0x80]) == false
    end
    # Check seven-byte sequences, should be invalid
    @test isvalid(String, UInt8[0xfe, 0x80, 0x80, 0x80, 0x80, 0x80]) == false

    # invalid Chars
    @test  isvalid('a')
    @test  isvalid('柒')
    @test !isvalid("\xff"[1])
    @test !isvalid("\xc0\x80"[1])
    @test !isvalid("\xf0\x80\x80\x80"[1])
    @test !isvalid('\ud800')
    @test  isvalid('\ud7ff')
    @test !isvalid('\udfff')
    @test  isvalid('\ue000')
end

@testset "NULL pointers are handled consistently by String" begin
    @test_throws ArgumentError unsafe_string(Ptr{UInt8}(0))
    @test_throws ArgumentError unsafe_string(Ptr{UInt8}(0), 10)
end
@testset "ascii for ASCII strings and non-ASCII strings" begin
    @test ascii("Hello, world") == "Hello, world"
    @test typeof(ascii("Hello, world")) == String
    @test ascii(GenericString("Hello, world")) == "Hello, world"
    @test typeof(ascii(GenericString("Hello, world"))) == String
    @test_throws ArgumentError ascii("Hello, ∀")
    @test_throws ArgumentError ascii(GenericString("Hello, ∀"))
end
@testset "issue #17271: lastindex() doesn't throw an error even with invalid strings" begin
    @test lastindex("\x90") == 1
    @test lastindex("\xce") == 1
end
# issue #17624, missing getindex method for String
@test "abc"[:] == "abc"

@testset "issue #18280: next/nextind must return past String's underlying data" begin
    for s in ("Hello", "Σ", "こんにちは", "😊😁")
        local s
        @test iterate(s, lastindex(s))[2] > sizeof(s)
        @test nextind(s, lastindex(s)) > sizeof(s)
    end
end
# Test cmp with AbstractStrings that don't index the same as UTF-8, which would include
# (LegacyString.)UTF16String and (LegacyString.)UTF32String, among others.

mutable struct CharStr <: AbstractString
    chars::Vector{Char}
    CharStr(x) = new(collect(x))
end
Base.iterate(x::CharStr, i::Integer=1) = iterate(x.chars, i)
Base.lastindex(x::CharStr) = lastindex(x.chars)
@testset "cmp without UTF-8 indexing" begin
    # Simple case, with just ANSI Latin 1 characters
    @test "áB" != CharStr("áá") # returns false with bug
    @test cmp("áB", CharStr("áá")) == -1 # returns 0 with bug

    # Case with Unicode characters
    @test cmp("\U1f596\U1f596", CharStr("\U1f596")) == 1   # Gives BoundsError with bug
    @test cmp(CharStr("\U1f596"), "\U1f596\U1f596") == -1
end

@testset "repeat" begin
    @inferred repeat(GenericString("x"), 1)
    @test repeat("xx",3) == repeat("x",6) == repeat('x',6) == repeat(GenericString("x"), 6) == "xxxxxx"
    @test repeat("αα",3) == repeat("α",6) == repeat('α',6) == repeat(GenericString("α"), 6) == "αααααα"
    @test repeat("x",1) == repeat('x',1) == "x"^1 == 'x'^1 == GenericString("x")^1 == "x"
    @test repeat("x",0) == repeat('x',0) == "x"^0 == 'x'^0 == GenericString("x")^0 == ""

    for S in ["xxx", "ååå", "∀∀∀", "🍕🍕🍕"]
        c = S[1]
        s = string(c)
        @test_throws ArgumentError repeat(c, -1)
        @test_throws ArgumentError repeat(s, -1)
        @test_throws ArgumentError repeat(S, -1)
        @test repeat(c, 0) == ""
        @test repeat(s, 0) == ""
        @test repeat(S, 0) == ""
        @test repeat(c, 1) == s
        @test repeat(s, 1) == s
        @test repeat(S, 1) == S
        @test repeat(c, 3) == S
        @test repeat(s, 3) == S
        @test repeat(S, 3) == S*S*S
    end
end
@testset "issue #12495: check that logical indexing attempt raises ArgumentError" begin
    @test_throws ArgumentError "abc"[[true, false, true]]
    @test_throws ArgumentError "abc"[BitArray([true, false, true])]
end

@testset "concatenation" begin
    @test "ab" * "cd" == "abcd"
    @test 'a' * "bc" == "abc"
    @test "ab" * 'c' == "abc"
    @test 'a' * 'b' == "ab"
    @test 'a' * "b" * 'c' == "abc"
    @test "a" * 'b' * 'c' == "abc"
end

@testset "unrecognized escapes in string/char literals" begin
    @test_throws Meta.ParseError Meta.parse("\"\\.\"")
    @test_throws Meta.ParseError Meta.parse("\'\\.\'")
end

@testset "thisind" begin
    let strs = Any["∀α>β:α+1>β", s"∀α>β:α+1>β",
                   SubString("123∀α>β:α+1>β123", 4, 18),
                   SubString(s"123∀α>β:α+1>β123", 4, 18)]
        for s in strs
            @test_throws BoundsError thisind(s, -2)
            @test_throws BoundsError thisind(s, -1)
            @test thisind(s, 0) == 0
            @test thisind(s, 1) == 1
            @test thisind(s, 2) == 1
            @test thisind(s, 3) == 1
            @test thisind(s, 4) == 4
            @test thisind(s, 5) == 4
            @test thisind(s, 6) == 6
            @test thisind(s, 15) == 15
            @test thisind(s, 16) == 15
            @test thisind(s, 17) == 17
            @test_throws BoundsError thisind(s, 18)
            @test_throws BoundsError thisind(s, 19)
        end
    end

    let strs = Any["", s"", SubString("123", 2, 1), SubString(s"123", 2, 1)]
        for s in strs
            @test_throws BoundsError thisind(s, -1)
            @test thisind(s, 0) == 0
            @test thisind(s, 1) == 1
            @test_throws BoundsError thisind(s, 2)
        end
    end
end

@testset "prevind and nextind" begin
    for s in Any["∀α>β:α+1>β", GenericString("∀α>β:α+1>β")]
        @test_throws BoundsError prevind(s, 0)
        @test_throws BoundsError prevind(s, 0, 0)
        @test_throws BoundsError prevind(s, 0, 1)
        @test prevind(s, 1) == 0
        @test prevind(s, 1, 1) == 0
        @test prevind(s, 1, 0) == 1
        @test prevind(s, 2) == 1
        @test prevind(s, 2, 1) == 1
        @test prevind(s, 4) == 1
        @test prevind(s, 4, 1) == 1
        @test prevind(s, 5) == 4
        @test prevind(s, 5, 1) == 4
        @test prevind(s, 5, 2) == 1
        @test prevind(s, 5, 3) == 0
        @test prevind(s, 15) == 14
        @test prevind(s, 15, 1) == 14
        @test prevind(s, 15, 2) == 13
        @test prevind(s, 15, 3) == 12
        @test prevind(s, 15, 4) == 10
        @test prevind(s, 15, 10) == 0
        @test prevind(s, 15, 9) == 1
        @test prevind(s, 16) == 15
        @test prevind(s, 16, 1) == 15
        @test prevind(s, 16, 2) == 14
        @test prevind(s, 17) == 15
        @test prevind(s, 17, 1) == 15
        @test prevind(s, 17, 2) == 14
        @test_throws BoundsError prevind(s, 18)
        @test_throws BoundsError prevind(s, 18, 0)
        @test_throws BoundsError prevind(s, 18, 1)

        @test_throws BoundsError nextind(s, -1)
        @test_throws BoundsError nextind(s, -1, 0)
        @test_throws BoundsError nextind(s, -1, 1)
        @test nextind(s, 0, 2) == 4
        @test nextind(s, 0, 20) == 26
        @test nextind(s, 0, 10) == 15
        @test nextind(s, 1) == 4
        @test nextind(s, 1, 1) == 4
        @test nextind(s, 1, 2) == 6
        @test nextind(s, 1, 9) == 15
        @test nextind(s, 1, 10) == 17
        @test nextind(s, 2) == 4
        @test nextind(s, 2, 1) == 4
        @test nextind(s, 3) == 4
        @test nextind(s, 3, 1) == 4
        @test nextind(s, 4) == 6
        @test nextind(s, 4, 1) == 6
        @test nextind(s, 14) == 15
        @test nextind(s, 14, 1) == 15
        @test nextind(s, 15) == 17
        @test nextind(s, 15, 1) == 17
        @test nextind(s, 15, 2) == 18
        @test nextind(s, 16) == 17
        @test nextind(s, 16, 1) == 17
        @test nextind(s, 16, 2) == 18
        @test nextind(s, 16, 3) == 19
        @test_throws BoundsError nextind(s, 17)
        @test_throws BoundsError nextind(s, 17, 0)
        @test_throws BoundsError nextind(s, 17, 1)

        for x in 0:ncodeunits(s)+1
            n = p = x
            for j in 1:40
                if 1 ≤ p
                    p = prevind(s, p)
                    @test prevind(s, x, j) == p
                end
                if n ≤ ncodeunits(s)
                    n = nextind(s, n)
                    @test nextind(s, x, j) == n
                end
            end
        end
    end
end

@testset "first and last" begin
    s = "∀ϵ≠0: ϵ²>0"
    @test_throws ArgumentError first(s, -1)
    @test first(s, 0) == ""
    @test first(s, 1) == "∀"
    @test first(s, 2) == "∀ϵ"
    @test first(s, 3) == "∀ϵ≠"
    @test first(s, 4) == "∀ϵ≠0"
    @test first(s, length(s)) == s
    @test first(s, length(s)+1) == s
    @test_throws ArgumentError last(s, -1)
    @test last(s, 0) == ""
    @test last(s, 1) == "0"
    @test last(s, 2) == ">0"
    @test last(s, 3) == "²>0"
    @test last(s, 4) == "ϵ²>0"
    @test last(s, length(s)) == s
    @test last(s, length(s)+1) == s
end

@testset "invalid code point" begin
    s = String([0x61, 0xba, 0x41])
    @test !isvalid(s)
    @test s[2] == reinterpret(Char, UInt32(0xba) << 24)
end

@testset "ncodeunits" begin
    for (s, n) in [""     => 0, "a"   => 1, "abc"  => 3,
                   "α"    => 2, "abγ" => 4, "∀"    => 3,
                   "∀x∃y" => 8, "🍕"  => 4, "🍕∀" => 7]
        @test ncodeunits(s) == n
        @test ncodeunits(GenericString(s)) == n
    end
end

@testset "0-step nextind and prevind" begin
    for T in [String, SubString, Base.SubstitutionString, GenericString]
        e = convert(T, "")
        @test nextind(e, 0, 0) == 0
        @test_throws BoundsError nextind(e, 1, 0)
        @test_throws BoundsError prevind(e, 0, 0)
        @test prevind(e, 1, 0) == 1

        s = convert(T, "∀x∃")
        @test nextind(s, 0, 0) == 0
        @test nextind(s, 1, 0) == 1
        @test_throws StringIndexError nextind(s, 2, 0)
        @test_throws StringIndexError nextind(s, 3, 0)
        @test nextind(s, 4, 0) == 4
        @test nextind(s, 5, 0) == 5
        @test_throws StringIndexError nextind(s, 6, 0)
        @test_throws StringIndexError nextind(s, 7, 0)
        @test_throws BoundsError nextind(s, 8, 0)

        @test_throws BoundsError prevind(s, 0, 0)
        @test prevind(s, 1, 0) == 1
        @test_throws StringIndexError prevind(s, 2, 0)
        @test_throws StringIndexError prevind(s, 3, 0)
        @test prevind(s, 4, 0) == 4
        @test prevind(s, 5, 0) == 5
        @test_throws StringIndexError prevind(s, 6, 0)
        @test_throws StringIndexError prevind(s, 7, 0)
        @test prevind(s, 8, 0) == 8
    end
end

@test unsafe_wrap(Vector{UInt8},"\xcc\xdd\xee\xff\x80") == [0xcc,0xdd,0xee,0xff,0x80]

@test iterate("a", 1)[2] == 2
@test nextind("a", 1) == 2
@test iterate("az", 1)[2] == 2
@test nextind("az", 1) == 2
@test iterate("a\xb1", 1)[2] == 2
@test nextind("a\xb1", 1) == 2
@test iterate("a\xb1z", 1)[2] == 2
@test nextind("a\xb1z", 1) == 2
@test iterate("a\xb1\x83", 1)[2] == 2
@test nextind("a\xb1\x83", 1) == 2
@test iterate("a\xb1\x83\x84", 1)[2] == 2
@test nextind("a\xb1\x83\x84", 1) == 2
@test iterate("a\xb1\x83\x84z", 1)[2] == 2
@test nextind("a\xb1\x83\x84z", 1) == 2

@test iterate("\x81", 1)[2] == 2
@test nextind("\x81", 1) == 2
@test iterate("\x81z", 1)[2] == 2
@test nextind("\x81z", 1) == 2
@test iterate("\x81\xb1", 1)[2] == 2
@test nextind("\x81\xb1", 1) == 2
@test iterate("\x81\xb1z", 1)[2] == 2
@test nextind("\x81\xb1z", 1) == 2
@test iterate("\x81\xb1\x83", 1)[2] == 2
@test nextind("\x81\xb1\x83", 1) == 2
@test iterate("\x81\xb1\x83\x84", 1)[2] == 2
@test nextind("\x81\xb1\x83\x84", 1) == 2
@test iterate("\x81\xb1\x83\x84z", 1)[2] == 2
@test nextind("\x81\xb1\x83\x84z", 1) == 2

@test iterate("\xce", 1)[2] == 2
@test nextind("\xce", 1) == 2
@test iterate("\xcez", 1)[2] == 2
@test nextind("\xcez", 1) == 2
@test iterate("\xce\xb1", 1)[2] == 3
@test nextind("\xce\xb1", 1) == 3
@test iterate("\xce\xb1z", 1)[2] == 3
@test nextind("\xce\xb1z", 1) == 3
@test iterate("\xce\xb1\x83", 1)[2] == 3
@test nextind("\xce\xb1\x83", 1) == 3
@test iterate("\xce\xb1\x83\x84", 1)[2] == 3
@test nextind("\xce\xb1\x83\x84", 1) == 3
@test iterate("\xce\xb1\x83\x84z", 1)[2] == 3
@test nextind("\xce\xb1\x83\x84z", 1) == 3

@test iterate("\xe2", 1)[2] == 2
@test nextind("\xe2", 1) == 2
@test iterate("\xe2z", 1)[2] == 2
@test nextind("\xe2z", 1) == 2
@test iterate("\xe2\x88", 1)[2] == 3
@test nextind("\xe2\x88", 1) == 3
@test iterate("\xe2\x88z", 1)[2] == 3
@test nextind("\xe2\x88z", 1) == 3
@test iterate("\xe2\x88\x83", 1)[2] == 4
@test nextind("\xe2\x88\x83", 1) == 4
@test iterate("\xe2\x88\x83z", 1)[2] == 4
@test nextind("\xe2\x88\x83z", 1) == 4
@test iterate("\xe2\x88\x83\x84", 1)[2] == 4
@test nextind("\xe2\x88\x83\x84", 1) == 4
@test iterate("\xe2\x88\x83\x84z", 1)[2] == 4
@test nextind("\xe2\x88\x83\x84z", 1) == 4

@test iterate("\xf0", 1)[2] == 2
@test nextind("\xf0", 1) == 2
@test iterate("\xf0z", 1)[2] == 2
@test nextind("\xf0z", 1) == 2
@test iterate("\xf0\x9f", 1)[2] == 3
@test nextind("\xf0\x9f", 1) == 3
@test iterate("\xf0\x9fz", 1)[2] == 3
@test nextind("\xf0\x9fz", 1) == 3
@test iterate("\xf0\x9f\x98", 1)[2] == 4
@test nextind("\xf0\x9f\x98", 1) == 4
@test iterate("\xf0\x9f\x98z", 1)[2] == 4
@test nextind("\xf0\x9f\x98z", 1) == 4
@test iterate("\xf0\x9f\x98\x84", 1)[2] == 5
@test nextind("\xf0\x9f\x98\x84", 1) == 5
@test iterate("\xf0\x9f\x98\x84z", 1)[2] == 5
@test nextind("\xf0\x9f\x98\x84z", 1) == 5

@test iterate("\xf8", 1)[2] == 2
@test nextind("\xf8", 1) == 2
@test iterate("\xf8z", 1)[2] == 2
@test nextind("\xf8z", 1) == 2
@test iterate("\xf8\x9f", 1)[2] == 2
@test nextind("\xf8\x9f", 1) == 2
@test iterate("\xf8\x9fz", 1)[2] == 2
@test nextind("\xf8\x9fz", 1) == 2
@test iterate("\xf8\x9f\x98", 1)[2] == 2
@test nextind("\xf8\x9f\x98", 1) == 2
@test iterate("\xf8\x9f\x98z", 1)[2] == 2
@test nextind("\xf8\x9f\x98z", 1) == 2
@test iterate("\xf8\x9f\x98\x84", 1)[2] == 2
@test nextind("\xf8\x9f\x98\x84", 1) == 2
@test iterate("\xf8\x9f\x98\x84z", 1)[2] == 2
@test nextind("\xf8\x9f\x98\x84z", 1) == 2

# codeunit vectors

let s = "∀x∃y", u = codeunits(s)
    @test u isa Base.CodeUnits{UInt8,String}
    @test length(u) == ncodeunits(s) == 8
    @test sizeof(u) == sizeof(s)
    @test eltype(u) === UInt8
    @test size(u) == (length(u),)
    @test strides(u) == (1,)
    @test u[1] == 0xe2
    @test u[2] == 0x88
    @test u[8] == 0x79
    @test_throws ErrorException (u[1] = 0x00)
    @test collect(u) == b"∀x∃y"
end

# issue #24388
let v = unsafe_wrap(Vector{UInt8}, "abc")
    s = String(v)
    @test_throws BoundsError v[1]
    push!(v, UInt8('x'))
    @test s == "abc"
end

# PR #25535
let v = [0x40,0x41,0x42]
    @test String(view(v, 2:3)) == "AB"
end

# make sure length for identical String and AbstractString return the same value, PR #25533
let rng = MersenneTwister(1), strs = ["∀εa∀aε"*String(rand(rng, UInt8, 100))*"∀εa∀aε",
                                   String(rand(rng, UInt8, 200))]
    for s in strs, i in 1:ncodeunits(s)+1, j in 0:ncodeunits(s)
            @test length(s,i,j) == length(GenericString(s),i,j)
    end
    for i in 0:10, j in 1:100,
        s in [randstring(rng, i), randstring(rng, "∀∃α1", i), String(rand(rng, UInt8, i))]
        @test length(s) == length(GenericString(s))
    end
end

# conversion of SubString to the same type, issue #25525
let x = SubString("ab", 1, 1)
    y = convert(SubString{String}, x)
    @test y === x
    chop("ab") === chop.(["ab"])[1]
end

max_stars_repo_name: davidanthoff/julia

max_issues_repo_name: davidanthoff/julia

max_forks_repo_name: davidanthoff/julia


text:
"""
    DocumenterGHActionContainer

[Repository](https://github.com/Nosferican/DocumenterGHActionContainer.jl)
"""
module DocumenterGHActionContainer
    for (root, dirs, files) ∈ walkdir(joinpath(@__DIR__))
        foreach(file -> include(joinpath(root, file)), filter!(file -> occursin(r"^\d{2}_\w+\.jl$", file), files))
    end
end


text:
@inline function ExaTronKernels.dgpstep(n::Int,x::CuDeviceArray{Float64,1},xl::CuDeviceArray{Float64,1},
                         xu::CuDeviceArray{Float64,1},alpha,w::CuDeviceArray{Float64,1},
                         s::CuDeviceArray{Float64,1})
    tx = threadIdx().x

    if tx <= n
        @inbounds begin
            # It might be better to process this using just a single thread,
            # rather than diverging between multiple threads.

            if x[tx] + alpha*w[tx] < xl[tx]
                s[tx] = xl[tx] - x[tx]
            elseif x[tx] + alpha*w[tx] > xu[tx]
                s[tx] = xu[tx] - x[tx]
            else
                s[tx] = alpha*w[tx]
            end
        end
    end
    CUDA.sync_threads()

    return
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

@testset "Repeat" begin
    function hex3RGB()
        clusterelements = SVector((0,0),(-1,0),(-1,1))
        colors = [colorant"red",colorant"green",colorant"blue"]
        names = ["R","G","B"]
        eltlattice = Repeat.HexBasis1()
        clusterbasis = Repeat.LatticeBasis(( -1,2),(2,-1))
        lattice = Repeat.LatticeCluster(clusterbasis,eltlattice,clusterelements)
        properties =  DataFrame(Color = colors, Name = names)
        return Repeat.ClusterWithProperties(lattice,properties)
    end

    """ Create a LatticeCluser with three elements at (0,0),(-1,0),(-1,1) coordinates in the HexBasis1 lattice"""
    function hex3cluster()
        clusterelts = SVector((0,0),(-1,0),(-1,1))
        eltlattice = Repeat.HexBasis1()
        clusterbasis = Repeat.LatticeBasis(( -1,2),(2,-1))
        return Repeat.LatticeCluster(clusterbasis,eltlattice,clusterelts)
    end

    @test [-1 2;2 -1] == Repeat.basismatrix(Repeat.clusterbasis(hex3RGB()))

    function basistest(a::Repeat.AbstractLatticeCluster)
        return Repeat.clusterbasis(a)
    end

    @test basistest(hex3cluster()) == basistest(hex3RGB())
end