text:
module BenchQueueHotPotato

using BenchmarkTools
using ConcurrentCollections
using ConcurrentCollections.Implementations: atomic_modifyfield!
using Formatting: format
using ..Utils: maptasks

const HOTPOTATO = false
const NOTPOTATO = true

const Stat = typeof((npush = 0, npop = 0, t1 = time_ns(), t2 = time_ns()))

Base.@kwdef struct HotPotatoResult
    duration::Union{Float64,Nothing}
    nrepeat::Int
    stats::Vector{Stat}
end

function unfair_sleep(seconds::Real)
    t0 = time_ns()
    ns = seconds * 1e9
    while time_ns() - t0 < ns
        GC.safepoint()
        # yield()
    end
end

@noinline function hotpotato!(
    q;
    ntasks::Integer = Threads.nthreads(),
    duration::Union{Real,Nothing} = nothing,
    nrepeat::Integer = duration === nothing ? 2^15 : typemax(Int64),
    delay::Real = 0,  # not using 0.001 mentioned in Izraelevitz & Scott (2017)
    nfills::Integer = 0,
    ppush::Real = 0.5,
)
    ntasks < 1 && error("require positive `ntasks`: got $ntasks")
    local tasks
    push!(q, HOTPOTATO)
    t0 = time_ns()
    tasks = maptasks(1:ntasks) do itask
        # Core.print("$itask-th task started\n")
        local t1 = time_ns()
        local npush = 0
        local npop = 0
        for irepeat in 1:nrepeat
            # ccall(:jl_breakpoint, Cvoid, (Any,), (; irepeat, itask))
            if rand() < ppush
                push!(q, NOTPOTATO)
                npush += 1
            else
                y = popfirst!(q)
                npop += 1
                if y == HOTPOTATO
                    for _ in 1:nfills
                        push!(q, NOTPOTATO)
                    end
                    unfair_sleep(delay)
                    # sleep(delay)
                    push!(q, y)
                    npush += 1
                end
            end
            if duration !== nothing
                if (time_ns() - t0) / 1e9 > duration
                    break
                end
            end
        end
        return (; npush, npop, t1, t2 = time_ns())::Stat
    end
    return HotPotatoResult(; nrepeat, duration, stats = map(fetch, tasks))
end

mutable struct AtomicRef{T}
    @atomic x::T
end

const FAIStat = typeof((; nfai = 0, t1 = time_ns(), t2 = time_ns()))

Base.@kwdef struct FAIResult
    duration::Union{Float64,Nothing}
    nrepeat::Int
    refvalue::Int
    stats::Vector{FAIStat}
end

function fai_stats(;
    ntasks::Integer = Threads.nthreads(),
    duration::Union{Real,Nothing} = nothing,
    nrepeat::Integer = duration === nothing ? 2^15 : typemax(Int64),
    impl::Val = Val(:threads),
)
    ntasks < 1 && error("require positive `ntasks`: got $ntasks")
    local tasks
    ref = if impl === Val(:threads)
        Threads.Atomic{Int32}(0)
    else
        AtomicRef{Int32}(0)
    end
    t0 = time_ns()
    tasks = maptasks(1:ntasks) do itask
        local t1 = time_ns()
        local nfai = 0
        for irepeat in 1:nrepeat
            if impl === Val(:threads)
                Threads.atomic_add!(ref, Int32(1))
            elseif impl === Val(:unsafe)
                atomic_modifyfield!(ref, Val(:x), +, Int32(1))
            else
                @atomic ref.x += true
            end
            nfai += 1
            if duration !== nothing
                if (time_ns() - t0) / 1e9 > duration
                    break
                end
            end
        end
        return (; nfai, t1, t2 = time_ns())::FAIStat
    end
    if impl === Val(:threads)
        refvalue = ref[]
    else
        refvalue = @atomic ref.x
    end
    return FAIResult(; nrepeat, duration, refvalue, stats = map(fetch, tasks))
end

function setup(; kwargs...)
    suite = BenchmarkGroup()
    suite["channel"] = @benchmarkable hotpotato!(Channel{Bool}(Inf); $kwargs...)
    suite["dlcrq"] =
        @benchmarkable hotpotato!(DualLinkedConcurrentRingQueue{Bool}(); $kwargs...)
    suite["dlq"] = @benchmarkable hotpotato!(DualLinkedQueue{Bool}(); $kwargs...)
    return suite
end

function clear() end

function summarize(result::HotPotatoResult)
    npush = sum(stat.npush for stat in result.stats)
    npop = sum(stat.npop for stat in result.stats)
    t1min = minimum(stat.t1 for stat in result.stats)
    t2max = maximum(stat.t2 for stat in result.stats)
    duration = (t2max - t1min) / 1e9
    throughput = (npush + npop) / duration
    return (; npush, npop, duration, throughput)
end

function Base.show(io::IO, ::MIME"text/plain", result::HotPotatoResult)
    (; npush, npop, duration, throughput) = summarize(result)
    println(io, "Hot potato benchmark with ", length(result.stats), " tasks")
    println(
        io,
        format(npush; commas = true),
        " pushes and ",
        format(npop; commas = true),
        " pops in ",
        duration,
        " seconds",
    )
    print(
        io,
        "throughput = ",
        format(floor(Int, throughput); commas = true),
        " ops/seconds",
    )
end

function summarize(result::FAIResult)
    nfai = sum(stat.nfai for stat in result.stats)
    t1min = minimum(stat.t1 for stat in result.stats)
    t2max = maximum(stat.t2 for stat in result.stats)
    duration = (t2max - t1min) / 1e9
    throughput = nfai / duration
    return (; nfai, duration, throughput)
end

function Base.show(io::IO, ::MIME"text/plain", result::FAIResult)
    (; nfai, duration, throughput) = summarize(result)
    println(io, "FAI benchmark with ", length(result.stats), " tasks")
    println(io, format(nfai; commas = true), " FAIs in ", duration, " seconds")
    print(
        io,
        "throughput = ",
        format(floor(Int, throughput); commas = true),
        " ops/seconds",
    )
end

function timings(result::HotPotatoResult)
    t1min = minimum(stat.t1 for stat in result.stats)
    table =
        Iterators.map(enumerate(result.stats)) do (taskid, stat)
            (; t1, t2) = stat
            return (
                (; t = (t1 - t1min) / 1e9, event = :begin, taskid),
                (; t = (t2 - t1min) / 1e9, event = :end, taskid),
            )
        end |>
        Iterators.flatten |>
        collect
    sort!(table; by = row -> row.t)
    return table
end

end  # module

ext: jl


text:
__precompile__()
module FirstOrderLogic

using AutoHashEquals
using Match

# Declarative and supportive code.
include("types.jl")
include("helpers.jl")
include("testing_helpers.jl")

# Workhorse code.
include("parser.jl")
include("renamed_quantifiers_form.jl")
include("quantified_variables_form.jl")
include("prenex_normal_form.jl")
include("skolem_normal_form.jl")
include("conjunctive_normal_form.jl")
include("removed_quantifier_form.jl")
include("resolution.jl")

end


text:
using TransformVariables, LogDensityProblems, DynamicHMC,
    Distributions, Statistics, StatsFuns, ForwardDiff

using Random

export dynamicHMC

# import Flux
function dynamicHMC(m :: JointDistribution, _data, method=logpdf, N=1000::Int)
    ℓ(pars) = logpdf(m, merge(pars, _data), method)

    t = xform(m,_data)
    P = TransformedLogDensity(t, ℓ)
    ∇P = ADgradient(Val(:ForwardDiff), P)
    results = mcmc_with_warmup(MersenneTwister(), ∇P, N; reporter=DynamicHMC.NoProgressReport());
    samples = TransformVariables.transform.(parent(∇P).transformation, results.chain)
end

function dynamicHMC(m :: JointDistribution, _data, ::Val{Inf})
    ℓ(pars) = logpdf(m, merge(pars, _data))

    t = xform(m,_data)
    P = TransformedLogDensity(t, ℓ)
    ∇P = ADgradient(Val(:ForwardDiff), P)

    # initialization
    rng = MersenneTwister()
    results = DynamicHMC.mcmc_keep_warmup(rng, ∇P, 0; reporter = NoProgressReport())
    steps = DynamicHMC.mcmc_steps(results.sampling_logdensity, results.final_warmup_state)

    (results, steps)
end

using ResumableFunctions

export stream
@resumable function stream(f::typeof(dynamicHMC), m :: JointDistribution, _data::NamedTuple)
    t = xform(m, _data)
    (results, steps) = dynamicHMC(m, _data, Val(Inf))
    Q = results.final_warmup_state.Q
    while true
        Q, tree_stats = DynamicHMC.mcmc_next_step(steps, Q)
        @yield (merge(t(Q.q), (_ℓ = Q.ℓq,)), tree_stats)
    end
end


text:
import Base: push!, empty!, isempty, first, last, length, size, append!, pop!, convert, pushfirst!, getindex, setindex!

# Activate bonus
const enable_bonus_tests = true

mutable struct CircularBuffer{T} <: AbstractVector{T}
    data::Array{T, 1}
    first::Integer
    length::Integer

    function CircularBuffer{T}(capacity::Integer) where {T}
        new{T}(Array{T, 1}(undef, capacity), 1, 0)
    end
end

CircularBuffer(capacity::Integer) = CircularBuffer{Any}(capacity)

capacity(cb::CircularBuffer) = length(cb.data)

isempty(cb::CircularBuffer) = iszero(cb.length)
isfull(cb::CircularBuffer) = cb.length == capacity(cb)

_move_pos(cb::CircularBuffer, i) = mod1(cb.first + i, capacity(cb))
_last_pos(cb::CircularBuffer) = _move_pos(cb, cb.length - 1)

function push!(cb::CircularBuffer{T}, item::T; overwrite::Bool=false) where {T}
    if isfull(cb)
        if !overwrite
            throw(BoundsError())
        else
            cb.first = _move_pos(cb, 1)
        end
    else
        cb.length += 1
    end

    pos = _last_pos(cb)
    cb.data[pos] = item
    return cb
end

function first(cb::CircularBuffer)
    if isempty(cb)
        throw(BoundsError())
    end
    return cb.data[cb.first]
end

function popfirst!(cb::CircularBuffer)
    value = first(cb)
    cb.first = _move_pos(cb, 1)
    cb.length -= 1
    return value
end

function last(cb::CircularBuffer)
    if isempty(cb)
        throw(BoundsError())
    end
    return cb.data[_last_pos(cb)]
end

function pop!(cb::CircularBuffer)
    value = last(cb)
    cb.length -= 1
    return value
end

function empty!(cb::CircularBuffer{T}) where {T}
    cb.data = Array{T, 1}(undef, capacity(cb))
    cb.length = 0
    return cb
end

length(cb::CircularBuffer) = cb.length

size(cb::CircularBuffer) = (length(cb),)

function append!(cb::CircularBuffer, items::AbstractArray; overwrite::Bool=false)
    for el in items
        push!(cb, el, overwrite=overwrite)
    end
    return cb
end

function pushfirst!(cb::CircularBuffer{T}, item::T; overwrite::Bool=false) where {T}
    if isfull(cb)
        if !overwrite
            throw(BoundsError())
        end
    else
        cb.length += 1
    end

    cb.first = cb.first == 1 ? capacity(cb) : cb.first - 1
    cb.data[cb.first] = item
    return cb
end

function getindex(cb::CircularBuffer, i::Integer)
    checkbounds(cb, i)
    return cb.data[_move_pos(cb, i-1)]
end

function setindex!(cb::CircularBuffer{T}, val::T, i::Integer) where {T}
    checkbounds(cb, i)
    cb.data[_move_pos(cb, i-1)] = val
    return cb
end

function convert(::Type{Array}, cb::CircularBuffer{T}) where {T}
    return T[item for item in cb]
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
using JLD2
include("../models/bpModel.jl")

setRNGs(12345)

@load "data/bpData.jld" ys θ
model = makeBPModel(θ, ys)

N = 400 * model.n

rpio = RPIO{model.Tx}(N, model.n, 1)
@test rpEstimateDiscrete(model, rpio) ≈ -5415.0 atol=1.0
rpio = RPIO{model.Tx}(N, model.n, 4)
@test rpEstimateDiscrete(model, rpio) ≈ -5415.0 atol=1.0


text:
module BayesianAssignment

import CSV
import CairoMakie
import Makie

using Reexport: @reexport
@reexport begin
using AlgebraOfGraphics
using Books
using DataFrames
using Dates
using Distributions
using Random
using StableRNGs
end

export build

include("data.jl")
export dataset, first_and_last_few_rows, generate_samples

include("plots.jl")

function build()
    println("Generating content")
    fail_on_error = true
    dataset()
    Books.gen(; fail_on_error)
    extra_head = """
        <script src="https://cdn.usefathom.com/script.js" data-site="WYKYRAWQ" defer></script>
        """
    Books.build_all(; extra_head, fail_on_error)
end

end # module


text:
using Assignment01
using Test

@testset "Assignment01" begin

@testset "Question 1" begin
    @test question1 isa AbstractString
    @test isdefined(Assignment01, :my_name)
    @test isdefined(Assignment01, :my_age)
    @test !occursin("kevin", lowercase(Assignment01.my_name))
    @test Assignment01.my_age < 30
    include("../src/run.jl")
    @test my_name == Assignment01.my_name
    @test my_age == Assignment01.my_age
    @test isnothing(main())
    @test lowercase(question1) == "yes"
end

@testset "Question 2" begin
    @test question2 isa AbstractString

    scmd = string.(split(question2))
    @test first(scmd) == "ls"

    push!(scmd, normpath(joinpath(@__DIR__, "..")))
    dircontents = split(read(Cmd(scmd), String), '\n')
    @test ".git" in dircontents
    @test length(dircontents) > 10
end

@testset "Question 3" begin
    rm("q3", force=true)
    @test question3 isa AbstractString
    scmd = string.(split(question3))
    @test first(scmd) == "mkdir"
    @test length(scmd) == 2
    run(Cmd(scmd))
    @test isdir("q3")
    rm("q3", force=true)
end

@testset "Question 4" begin
    @test question4 isa Real
    @test question4 == 3
end

@testset "Question 5" begin
    @test question5 isa Float64
    @test 0 < question5 < 1
    @test question5 ≈ 7//10
end

end # testset Assignment01


text:
"""
Jarque–Bera Statistic as a normality measure for optimization
"""
function normality(x::AbstractVector{T}) where T<:Real

    P = Base.promote_op(/, T, T)

    # moments at origin
    n = P(length(x))
    m1 = mo2 = mo3 = mo4 = P(0)

    @inbounds for xi in x
        m1 += x1 = xi / n
        mo2 += x2 = x1 * x1 * n
        mo3 += x3 = x1 * x2 * n
        mo4 += x4 = x1 * x3 * n
    end

    # moments
    m2 = mo2 - m1^2
    m3 = mo3 - P(3)*m1*mo2 + P(2)*m1^3
    m4 = mo4 - P(4)*m1*mo3 + P(6)*m1^2*mo2 - P(3)*m1^2*m1^2

    S = m3 / (m2 * sqrt(m2))
    K = m4 / m2^2

    n/P(6) * (S^2 + P(1/4) * (K-P(3))^2)

end


text:
# ## benchmarking
using BenchmarkTools
using InteractiveUtils

cost!(prob.solver_data, prob.problem, mode=:nominal)
@benchmark cost!($prob.solver_data, $prob.problem, mode=:nominal)
@code_warntype cost!(prob.solver_data, prob.problem, mode=:nominal)

gradients!(prob.problem, mode=:nominal)
@benchmark gradients!($prob.problem, mode=:nominal)
@code_warntype gradients!(prob.problem, mode=:nominal)

mode = :nominal
IterativeLQR.backward_pass!(prob.policy, prob.problem, mode=mode)
@benchmark IterativeLQR.backward_pass!($prob.policy, $prob.problem, mode=:nominal)# setup=(mode=:nominal)
@benchmark IterativeLQR.backward_pass!($prob.policy, $prob.problem, mode=:cand)# setup=(mode=:nominal)
@code_warntype IterativeLQR.backward_pass!(prob.policy, prob.problem, mode=mode)

lagrangian_gradient!(prob.solver_data, prob.policy, prob.problem)
@benchmark lagrangian_gradient!($prob.solver_data, $prob.policy, $prob.problem)
@code_warntype lagrangian_gradient!(prob.solver_data, prob.policy, prob.problem)

IterativeLQR.trajectory_sensitivities(prob.problem, prob.policy, prob.solver_data)
@benchmark IterativeLQR.trajectory_sensitivities($prob.problem, $prob.policy, $prob.solver_data)
@code_warntype IterativeLQR.trajectory_sensitivities(prob.problem, prob.policy, prob.solver_data)

rollout!(prob.policy, prob.problem, step_size=prob.solver_data.step_size[1])
@benchmark rollout!($prob.policy, $prob.problem, step_size=$prob.solver_data.step_size[1])
@code_warntype rollout!(prob.policy, prob.problem, step_size=prob.solver_data.step_size[1])

augmented_lagrangian_update!(prob.problem.objective)
@benchmark augmented_lagrangian_update!($prob.problem.objective)
@code_warntype augmented_lagrangian_update!(prob.problem.objective)

constraint_violation(prob.problem.objective.constraint_data)
@benchmark constraint_violation($prob.problem.objective.constraint_data)
@code_warntype constraint_violation(prob.problem.objective.constraint_data)

initialize_controls!(prob, ū)
@benchmark initialize_controls!($prob, $ū)
@code_warntype initialize_controls!(prob, ū)

initialize_states!(prob, x̄)
@benchmark initialize_states!($prob, $x̄)
@code_warntype initialize_states!(prob, x̄)

function _forward_pass!(prob::Solver, x, u)
    initialize_controls!(prob, u)
    initialize_states!(prob, x)

    forward_pass!(prob.policy, prob.problem, prob.solver_data,
        min_step_size = 1.0e-5,
        line_search = :armijo)
end

_forward_pass!(prob, x̄, ū)
@benchmark _forward_pass!($prob, $x̄, $ū)
@code_warntype _forward_pass!(prob, x̄, ū)

A = Diagonal(ones(5))
a = rand(5)
b = ones(Int, 5)
using BenchmarkTools
@benchmark $A .= $(Diagonal(b .* a))

t = 1
a = prob.problem.objective.constraint_penalty[t]
b = prob.problem.objective.a[t]
c = prob.problem.objective.constraint_penalty_matrix[t]
@benchmark $c .= $(Diagonal(a .* b))


text:
rx = r"\s*\S+\s*"

words(text) = map(
    m -> (m.offset, m.offset + textlen(m.match), m.match, textlen(m.match)),
    eachmatch(rx, text),
)

function characters(word)
    chars = collect(word)
    widths = textwidth.(chars)
    return map(i -> (chars[i], widths[i]), 1:length(chars))
end

function style_at_each_line(text)
    lines = split(text, "\n")
    for (i, line) in enumerate(lines)
        for tag in eachmatch(OPEN_TAG_REGEX, line)
            markup = tag.match[2:(end - 1)]
            close_rx = r"(?<!\{)\{(?!\{)\/" * markup * r"\}"
            close_tag = match(close_rx, text)

            if isnothing(close_tag) && i < length(lines)
                lines[i + 1] = "{$markup}" * lines[i + 1]
            end
        end
    end
    return join(lines, '\n')
end

function split_tags_into_words(text)
    tags = map(m -> m.match[2:(end - 1)], eachmatch(OPEN_TAG_REGEX, text))

    for markup in tags
        tag = match(Regex("\\{$markup\\}"), text)
        isnothing(tag) && continue

        close_rx = r"(?<!\{)\{(?!\{)\/" * markup * r"\}"
        close_tag = match(close_rx, text)

        isnothing(close_tag) && continue

        tag_words = map(
            w -> replace(replace(w, "{" => ""), "}" => ""),
            ([w[3] for w in words(tag.match)]),
        )
        if length(tag_words) > 1
            openers = join(map(w -> "{" * rstrip(w) * "}", tag_words))
            closers = join(map(w -> "{/" * rstrip(w) * "}", tag_words))

            try
                text =
                    text[1:(tag.offset - 1)] *
                    openers *
                    text[(tag.offset + textwidth(markup) + 2):(close_tag.offset - 1)] *
                    closers *
                    text[(close_tag.offset + textwidth(markup) + 3):end]
            catch
                text =
                    text[1:prevind(text, tag.offset - 1)] *
                    openers *
                    text[nextind(text, tag.offset + textwidth(markup) + 2):prevind(
                        text,
                        close_tag.offset - 1,
                    )] *
                    closers *
                    text[prevind(text, close_tag.offset + textwidth(markup) + 3):end]
            end
        end
    end
    return text
end

"""
    reshape_text(text::AbstractString, width::Int)

Reshape a text to have a given width.

Insert newline characters in a string so that each line is within the given width.
"""
function reshape_text(text::AbstractString, width::Int)
    occursin('\n', text) && (return do_by_line(ln -> reshape_text(ln, width), text))

    textlen(text) ≤ width && return text
    text = split_tags_into_words(text)
    position = 0
    cuts = []
    for (start, _, word, word_length) in words(text)
        if position + word_length > width
            len_t = length(text[1:start])
            # we need to cut the text
            if word_length > width
                # the word is longer than the line
                word_chars = characters(word)
                position > 0 && push!(cuts, len_t + length(cuts))
                _, position = first(word_chars)  # width of first char

                for (i, (char, w)) in enumerate(word_chars[2:end])
                    if position + w > width
                        push!(cuts, len_t + length(cuts) + i)
                        position = w
                    else
                        position += w
                    end
                end

            elseif position > 0 && start > 1
                push!(cuts, len_t + length(cuts))
                position = word_length
            end
        else
            position += word_length
        end
    end

    chars = collect(text)
    for cut in cuts
        insert!(chars, cut, '\n')
    end

    out = apply_style(style_at_each_line(String(chars)))
    # do_by_line(ln -> ln * "\e[0m", out)
end


text:
using LLLplus, Documenter, MacroTools

# make.jl still doesn't seem to be running doctests.  Or at least it isn't
# running the doctests I expect it to run
setup = quote
    using LLLplus
end
DocMeta.setdocmeta!(LLLplus, :DocTestSetup, setup; recursive = true)

makedocs(
         sitename="LLLplus.jl",
         pages = [
                  "Home" => "index.md",
                  "Functions" => "functions.md"],
         format = Documenter.HTML(prettyurls = haskey(ENV, "CI")),
         strict = true,
)

deploydocs(
    repo   = "github.com/christianpeel/LLLplus.jl.git"
)


text:
module ReinforcementLearningBase

const RLBase = ReinforcementLearningBase
export RLBase

using Random
using Reexport

@reexport using CommonRLSpaces

include("inline_export.jl")
include("interface.jl")
include("CommonRLInterface.jl")
include("base.jl")

end # module


text:
export cgls

include("cgls.jl")


text:
module ClassicalDicke
export q_of_ϵ,q_sign,minimum_ϵ_for,Point,Pointθϕ,Pointθφ,discriminant_of_q_solution,
       energy_shell_volume,  maximum_P_for_ϵ,minimum_nonnegative_Q_for_ϵ,maximum_Q_for_ϵ,
       minimum_energy,minimum_energy_point, normal_frequency,phase_space_dist_squared,
       ClassicalDickeSystem,hamiltonian,DickeSystem
    using ..ClassicalSystems
    using ..PhaseSpaces
    using ParameterizedFunctions
    using QuadGK

    dicke_step = @ode_def begin
        dQ = integate_backwards*(P*ω₀-γ*P*q*Q/sqrt(4-P^2-Q^2)) # =dH/dP
        dq = integate_backwards*p*ω # =dH/dp
        dP = integate_backwards*(γ*q*Q^2/sqrt(4-P^2-Q^2)-γ*q*sqrt(4-P^2-Q^2)-Q*ω₀) # =-dH/dQ
        dp = integate_backwards*(-γ*Q*sqrt(4-P^2-Q^2)-q*ω) #-dH/dq

    end ω₀ ω γ integate_backwards
    function dicke_out_of_bounds(u,_,t)
        return 4-u[3]^2-u[1]^2<=0
    end
    varnames=[:Q,:q,:P,:p]
    """
    ```julia
    abstract type DickeSystem <: ClassicalSystems.ClassicalSystem
    ```
    Abstract subtype of [`ClassicalSystems.ClassicalSystem`](@ref) which represents the
    classical Dicke model.
    To generate it, use the concrete subtypes [`ClassicalDickeSystem`](@ref) or [`DickeBCE.QuantumDickeSystem`](@ref Main.DickeBCE.QuantumDickeSystem).
    """
    abstract type DickeSystem <: ClassicalSystems.ClassicalSystem end

    """
    ```julia
    struct ClassicalDickeSystem <: DickeSystem
    ```
    Subtype of [`ClassicalSystems.ClassicalSystem`](@ref) which represents the
    classical Dicke model with the given parameters ``ω_0``, ``ω``, and ``γ``. See Eq. (5) of Ref. [Pilatowsky2021](@cite).
    To generate this struct, use the constructor
    ```julia
    ClassicalDickeSystem(;ω₀::Real,ω::Real,γ::Real)
    ```
    For example, `system = ClassicalDickeSystem(ω₀=1, ω=1, γ=1)`.

    This struct may be passed to all functions in this module that require an instance of `DickeSystem`,
    as well as functions in other modules that require [`ClassicalSystems.ClassicalSystem`](@ref),
    such as  [`ClassicalSystems.integrate`](@ref ClassicalSystems.integrate(::ClassicalSystems.ClassicalSystem,::AbstractVector{<:Real},::Real)).
    """
    struct ClassicalDickeSystem <: DickeSystem
        parameters::Vector{Float64}
        ClassicalDickeSystem(;ω₀::Real,ω::Real,γ::Real)= new(Float64.([ω₀,ω,γ]))

    end
    function ClassicalSystems.parameters(system::ClassicalDickeSystem)
        return system.parameters
    end
    function ClassicalSystems.step(system::DickeSystem)
        return dicke_step
    end
    function ClassicalSystems.out_of_bounds(system::DickeSystem)
        return dicke_out_of_bounds
    end
    function ClassicalSystems.varnames(system::DickeSystem)
        return varnames
    end
    function Base.show(io::IO, cds::ClassicalDickeSystem)
        ω₀,ω,γ=ClassicalSystems.parameters(cds)
        print(io,"ClassicalDickeSystem(ω₀ = $ω₀, ω = $ω, γ = $γ)")
    end
    """
    ```julia
    function hamiltonian(system::DickeSystem)
    ```
    Returns a classical Hamiltonian function `h(x)` where `x=[Q,q,P,p]`, which is given by
    Eq. (5) of Ref. [Pilatowsky2021](@cite).

    # Arguments
    - `system` should be a subtype of [`DickeSystem`](@ref ClassicalDicke.DickeSystem).

    """
    function hamiltonian(system::DickeSystem)
        ω₀,ω,γ=ClassicalSystems.parameters(system)
        function H(u)
            Q,q,P,p=u
            r=1-(Q^2+P^2)/4
            (ω₀/2)*(Q^2+P^2)+(ω/2)*(q^2+p^2)+2*γ*q*Q*sqrt(r)-ω₀
        end
        return H
    end

    """
    ```julia
    function energy_shell_volume(system::DickeSystem,ϵ::Real)
    ```
    Returns the volume of the classical energy shell in the phase space, that is,
    ```math
    \\mathcal{V}(\\mathcal{M}_ϵ) = \\int_\\mathcal{M} \\text{d}\\mathbf{x} \\, \\delta(h_\\text{cl}(\\mathbf{x})- \\epsilon).
    ```

    This is computed using a modified version of
    Eq. (A2) in [Pilatowsky2021Identification](@cite), (the integral is
    performed in the variable ``Q`` instead of ``y``). That equation was originally derived in Ref. [Bastarrachea2014](@cite), with
    different units for ``\\epsilon``.

    The integration is performed using [`QuadGK`](https://github.com/JuliaMath/QuadGK.jl).
    # Arguments
    - `system` should be a subtype of [`DickeSystem`](@ref ClassicalDicke.DickeSystem).

    - `ϵ` is the energy.
    """
    function energy_shell_volume(system::DickeSystem,ϵ::Real)
        ω₀,ω,γ=ClassicalSystems.parameters(system)
        if ϵ <= minimum_energy(system)
            return 0.0
        end
        if ϵ>=ω₀
            return 8*pi^2/ω
        end

        #This function gives the height of P for a given Q, we then integrate
        #it over the allowed positive values for Q (and multiply by 2*pi*4/ω).
        Pmax(Q)=sqrt((ω*(2*ϵ + 2*ω₀) - Q^2*((-4 + Q^2)*γ^2 + ω*ω₀))/(Q^2*γ^2 + ω*ω₀))
        return quadgk(Pmax,minimum_nonnegative_Q_for_ϵ(system,ϵ),maximum_Q_for_ϵ(system,ϵ))[1]*2*pi*4/ω
    end

    """
    ```julia
    function normal_frequency(system::DickeSystem,
        sgn::Union{typeof(-),typeof(+)} = +)
    ```
    Returns the ground-state normal frequency, that is, ``\\Omega_{\\epsilon_{\\text{GS}}}^{A,B}`` at the bottom of page 3 of Ref. [Pilatowsky2021](@cite).

    # Arguments
    - `system` should be a subtype of [`DickeSystem`](@ref ClassicalDicke.DickeSystem).

    - `sgn` is `-` for ``\\Omega^A`` and `+` for ``\\Omega^B``. Defaults to `+`.

    Note: This function currently only works for the supperadiant phase.
    """
    function normal_frequency(system::DickeSystem,sgn::Union{typeof(-),typeof(+)}=+)
        ω₀,ω,γ=ClassicalSystems.parameters(system)
        γc=sqrt(ω₀*ω)/2
        if γc>γ
            error("No methods for the normal phase (please add it!)")
            return
        end
        return sqrt((1/(2*γc^4))*sgn(ω^2*γc^4 + ω₀^2*γ^4, sqrt((ω₀^2*γ^4 -ω^2*γc^4)^2 + 4*γc^8*ω^2*ω₀^2)))
    end

    """
    ```julia
    function minimum_energy_point(system::DickeSystem,
        Qsign::Union{typeof(-),typeof(+)} = +)
    ```
    Returns the ground-state coordinate, that is,  `(0, 0, 0, 0)` for the normal
    phase and ``\\mathbf{x}_\\text{GS}`` below Eq. (7) of Ref. [Pilatowsky2021](@cite) for the
    superradiant phase.

    # Arguments
    - `system` should be a subtype of [`DickeSystem`](@ref ClassicalDicke.DickeSystem).

    - `Qsign` toggles the sign of the ``Q`` coordinate in the superradiant phase,
      that is, `+` for ``\\mathbf{x}_\\text{GS}`` and  `-` for ``\\widetilde{\\mathbf{x}}_\\text{GS}`` . Defaults to `+`.
    """
    function minimum_energy_point(system::DickeSystem,signoQ::Union{typeof(-),typeof(+)}=+)
        local signoq
        ω₀,ω,γ=ClassicalSystems.parameters(system)
        γc=sqrt(ω₀*ω)/2
        if γc >= γ
            return Point(Q=0,P=0,p=0,q=0)
        end
        if signoQ==+
            signoq=-
        else
            signoq=+
        end
        return Point(Q=signoQ(sqrt(2-ω*ω₀/(2*γ^2))),q=signoq(sqrt(4*γ^2/ω^2 - ω₀^2/(4*γ^2))),P=0,p=0)
    end

    """
    ```julia
    function minimum_energy(system::DickeSystem)
    ```
    Returns the ground-state energy, as given by multiplying Eq. (15) of Ref. [Bastarrachea2014](@cite) by
    ``\\omega_0`` (because in that work they use ``\\epsilon =E/(\\omega_0 j)``, and in this module we
    use ``\\epsilon = E/j``).

    # Arguments
    - `system` should be a subtype of [`DickeSystem`](@ref ClassicalDicke.DickeSystem).

    """
    function minimum_energy(system::DickeSystem)
        ω₀,ω,γ = ClassicalSystems.parameters(system)
        γc=sqrt(ω₀*ω)/2

        if γ ≤ γc
            return -ω₀
        else
            return - ω₀*( (γc/γ)^2 + (γ/γc)^2 )/2
        end
    end

    """
    ```julia
    function discriminant_of_q_solution(system::DickeSystem;
        Q::Real,
        P::Real,
        p::Real,
        ϵ::Real)
    ```
    Returns the discriminant of the second degree equation in ``q`` given by
    ```math
        h_\\text{cl}(Q,q,P,p)=\\epsilon,
    ```
    where ``h_\\text{cl}`` is given by Eq. (5) of Ref. [Pilatowsky2021](@cite).

    # Arguments
    - `system` should be a subtype of [`DickeSystem`](@ref ClassicalDicke.DickeSystem).

    # Keyword arguments
    - `Q`, `P`, `p`, and `ϵ` are the values of ``Q``, ``P``, ``p``, and ``\\epsilon``, respectively.
    """
    function discriminant_of_q_solution(system::DickeSystem; Q::Real,P::Real,p::Real,ϵ::Real)
        ω₀,ω,γ=ClassicalSystems.parameters(system)
        r=1-(Q^2+P^2)/4
        return 4*γ^2*Q^2*r - p^2*ω^2 - ω*ω₀*(P^2 + Q^2 - 2) + 2*ω*ϵ
    end

    """
    ```julia
    function q_of_ϵ(system::DickeSystem;
        Q::Real,
        P::Real,
        p::Real,
        ϵ::Real,
        sgn::Union{typeof(-),typeof(+)}=+,
        returnNaNonError::Bool=true)
    ```
    Returns the solutions ``q_\\pm`` of the second degree equation in ``q`` given by
    ```math
        h_\\text{cl}(Q,q,P,p)=\\epsilon,
    ```
    where ``h_\\text{cl}`` is given by Eq. (5) of Ref. [Pilatowsky2021](@cite).

    # Arguments
    - `system` should be a subtype of [`DickeSystem`](@ref ClassicalDicke.DickeSystem).

    # Keyword arguments
    - `Q`, `P`, `p`, and `ϵ` are values of ``Q``, ``P``, ``p``, and ``\\epsilon``, respectively.
    - `sgn` is `+` for ``q_+`` and `-` for ``q_-``
    - If `returnNaNonError` is `true`, then `NaN` is returned if there are no solutions. If it is `false`, and error is raised.
    """
    function q_of_ϵ(system::DickeSystem;
        Q::Real,
        P::Real,
        p::Real,
        ϵ::Real,
        sgn::Union{typeof(-),typeof(+)}=+,
        returnNaNonError::Bool=true)
        Δ=discriminant_of_q_solution(system;Q=Q,P=P,p=p,ϵ=ϵ)
        if Δ<0
            if returnNaNonError
                return NaN
            else
                ϵmin= minimum_ϵ_for(system;Q=Q,P=P,p=p)
                error("The minimum energy for the given parameters is $ϵmin. There is no q such that (Q=$Q,q,P=$P,p=$p) has ϵ=$ϵ.")
            end
        end
        ω₀,ω,γ=ClassicalSystems.parameters(system)
        return sgn(-2*γ*Q*sqrt(1-(Q^2+P^2)/4),sqrt(Δ))/ω
    end
    """
    ```julia
    function q_sign(system::DickeSystem,
        x::AbstractVector{<:Real},
        ϵ::Real=hamiltonian(system)(x))
    ```
    Returns the sign of the root of the second degree equation in ``q`` given by
    ```math
        h_\\text{cl}(Q,q,P,p)=\\epsilon.
    ```
    That is, this function returns `+` if `q=x[2] ≈ q_of_ϵ(system;Q=x[1],P=x[3],p=x[4],ϵ=ϵ,sgn=+)`
    and returns `-` if `q=x[2] ≈ q_of_ϵ(system;Q=x[1],P=x[3],p=x[4],ϵ=ϵ,sgn=-)`.
    # Arguments
    - `system` should be a subtype of [`DickeSystem`](@ref ClassicalDicke.DickeSystem).

    - `x` is a vector in the form `[Q,q,P,p]`.
    - `ϵ` should be the energy of `x`. If this is not passed, it is computed using [`hamiltonian`](@ref ClassicalDicke.hamiltonian)`(system)(x)`.
    """
    function q_sign(system::DickeSystem,
        x::AbstractVector{<:Real},
        ϵ::Real=hamiltonian(system)(x))
        Q,q,P,p = x
        q₊ = q_of_ϵ(system; Q=Q, P=P, p=p, ϵ=ϵ, sgn=+)
        q₋ = q_of_ϵ(system; Q=Q, P=P, p=p, ϵ=ϵ, sgn=-)
        if abs(q - q₊) < abs(q - q₋)
          return +
        else
          return -
        end
    end
    """
    ```julia
    function minimum_ϵ_for(system::DickeSystem;
        Q::Union{Real,Nothing}=nothing,
        q::Union{Real,Nothing}=nothing,
        P::Union{Real,Nothing}=nothing,
        p::Union{Real,Nothing}=nothing)
    ```
    Returns the minimum energy ``\\epsilon`` when constraining the system to three fixed values of the coordinates ``Q``, ``q``, ``P``, ``p``.

    # Arguments
    - `system` should be a subtype of [`DickeSystem`](@ref ClassicalDicke.DickeSystem).

    # Keyword arguments
    - You may pass either ``(Q,q,P)`` or ``(q,P,p)``. The other combinanations are not implemented.

    This function can be especially useful to draw contours of the available phase space (see [Drawing contours of the available phase space](@ref))
    """
    function minimum_ϵ_for(system::DickeSystem;
        Q::Union{Real,Nothing}=nothing,
        q::Union{Real,Nothing}=nothing,
        P::Union{Real,Nothing}=nothing,
        p::Union{Real,Nothing}=nothing)
        if count(x->x==nothing,[q,Q,p,P])!=1
            error("You have to give three of Q,q,P,p")
        end
        ω₀,ω,γ=ClassicalSystems.parameters(system)
        if q==nothing
            r=1-(Q^2+P^2)/4
            return -(4*γ^2*Q^2*r - p^2*ω^2 - ω*ω₀*(P^2 + Q^2 - 2))/(2*ω)
        end
        if Q==nothing
            return (2*p^2*ω + 2*q^2*ω + P^2*ω₀- sqrt((-4 + P^2)^2* (4*q^2*γ^2 + ω₀^2)))/4
        end
        error("Sorry, that combinanation of variables is not implemented")
    end
    """
    ```julia
    function maximum_P_for_ϵ(system::DickeSystem,ϵ::Real)
    ```
    Computes the maximum value of the parameter ``P`` accessible to the system at energy ``\\epsilon``.

    # Arguments
    - `system` should be a subtype of [`DickeSystem`](@ref ClassicalDicke.DickeSystem).

    - `ϵ` is the energy.
    """
    function maximum_P_for_ϵ(system::DickeSystem,ϵ::Real)
        ω₀,ω,γ=ClassicalSystems.parameters(system)
        if ϵ>=ω₀
            return 2.0
        elseif ω*ω₀>sqrt(2)*γ*sqrt(ω*(ω₀-ϵ))
            return sqrt(2*(ω₀+ϵ)/ω₀)
        else
            return sqrt(clamp((2*√(2*γ^6*ω*(ω₀-ϵ)))/(-γ^4) + ω*ω₀/(γ^2)+4,0,4))
        end
    end
    """
    ```julia
    function maximum_Q_for_ϵ(system::DickeSystem,ϵ::Real)
    ```
    Computes the maximum value of the parameter ``Q`` accessible to the system at energy ``\\epsilon``.

    # Arguments
    - `system` should be a subtype of [`DickeSystem`](@ref ClassicalDicke.DickeSystem).

    - `ϵ` is the energy.
    """
    function maximum_Q_for_ϵ(system::DickeSystem,ϵ::Real)
        ω₀,ω,γ=ClassicalSystems.parameters(system)
        if ϵ>=ω₀
            return 2.0
        end
         return sqrt(clamp((4γ^2-ω*ω₀)/(2*γ^2) + (1/2)*sqrt(max(0,(16*γ^4+8γ^2*ϵ*ω+ω^2*ω₀^2)/γ^4)),0,4))
    end
    """
    ```julia
    function minimum_nonnegative_Q_for_ϵ(system::DickeSystem,ϵ::Real)
    ```
    Computes the minimum nonnegative value of the parameter ``Q`` accessible to the system at energy ``\\epsilon``.

    # Arguments
    - `system` should be a subtype of [`DickeSystem`](@ref ClassicalDicke.DickeSystem).

    - `ϵ` is the energy.
    """
    function minimum_nonnegative_Q_for_ϵ(system::DickeSystem,ϵ::Real)
        ω₀,ω,γ=ClassicalSystems.parameters(system)
        if ϵ >=- ω₀
            return 0.0
        end
         return sqrt(clamp((4γ^2-ω*ω₀)/(2*γ^2) - (1/2)*sqrt(max(0,(16*γ^4+8γ^2*ϵ*ω+ω^2*ω₀^2)/γ^4)),0,4))
    end

    """
    ```julia
    function Point(;Q::Real,q::Real,P::Real,p::Real)
    ```
    Returns the array `[Q,q,P,p]`. This ordering is used throughout this package.
    """
    Point(;Q::Real,q::Real,P::Real,p::Real) = Float64[Q,q,P,p]

    """
    ```julia
    function Point(system::DickeSystem;
        Q::Real,
        P::Real,
        p::Real,
        ϵ::Real,
        sgn::Union{typeof(-),typeof(+)} = +)
    ```
    Returns a list `[Q,q,P,p]`, where `q` is calculated with [`ClassicalDicke.q_of_ϵ`](@ref).
    If there are no solutions for ``q``, an error is raised.
    """
    Point(system::DickeSystem;
        Q::Real,
        P::Real,
        p::Real,
        ϵ::Real,
        sgn::Union{typeof(-),typeof(+)} = +) = Point(q=q_of_ϵ(system,Q=Q,P=P,p=p,ϵ=ϵ,sgn=sgn,returnNaNonError=false),
                                                        P=P, p=p, Q=Q)

    """
    ```julia
    function Pointθϕ(;θ::Real,ϕ::Real,q::Real,p::Real)
    ```
    Returns a list `[Q,q,P,p]`, where `Q` and `P` are calculated from `θ` and `ϕ` using [`PhaseSpaces.Q_of_θϕ`](@ref) and  [`PhaseSpaces.P_of_θϕ`](@ref).
    """
    Pointθϕ(;θ::Real,ϕ::Real,q::Real,p::Real)=Point(Q=Q_of_θϕ(θ,ϕ),q=q,P=P_of_θϕ(θ,ϕ),p=p)
    Pointθφ(;θ::Real,φ::Real,q::Real,p::Real)=Pointθϕ(θ=θ,ϕ=φ,q=q,p=p)

    """
    ```julia
    function phase_space_dist_squared(x::AbstractVector{<:Real},y::AbstractVector{<:Real})
    ```
    Returns the phase-space distance ``d_{\\mathcal{M}}(\\mathbf{x},\\mathbf{y})`` (See App. C of Ref [Pilatowsky2021](@cite)),
    where `x` and `y` are vectors in the form `[Q,q,P,p]`.
    """
    function phase_space_dist_squared(x::AbstractVector{<:Real},y::AbstractVector{<:Real})
        Q1,q1,P1,p1=x
        Q2,q2,P2,p2=y
        return (q1-q2)^2 + (p1-p2)^2 + arc_between_QP(Q1,P1,Q2,P2)^2
    end
    # """
    # ```julia
    # function WignerHWxSU2_fixed_ϵ(system::DickeSystem;Q,P,p,ϵ,j,signoq=+)
    # ```
    # Returns an instance of [`TWA.PhaseSpaceDistribution`](@ref) that samples points from the classical energy shell
    # using a [Random Walk Metropolis-Hastings algorithm](https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm) implemented in [Mamba](https://mambajl.readthedocs.io/en/latest/samplers/rwm.html?highlight=RWMVariate).
    # """
    #  function WignerHWxSU2_fixed_ϵ(system::DickeSystem;Q,P,p,ϵ,j,signoq=+)
    #         Δ(Q,P)=discriminant_of_q_solution(system;Q=Q,P=P,p=p,ϵ=ϵ)
    #         x₀=Point(system,Q=Q,P=P,p=p,ϵ=ϵ,sgn=signoq);
    #         dst=TWA.coherent_Wigner_HWxSU2(x₀,j);
    #         function probability_density(u)
    #             Q,q,P,p=u
    #             return dst.probability_density(Point(system,Q=Q,P=P,p=p,ϵ=ϵ,sgn=signoq))/sqrt(Δ(Q,P))
    #         end
    #         function logf(u)
    #             θ,ϕ,p=u
    #             Q=PhaseSpaces.Q_of_θϕ(θ,ϕ)
    #             P=PhaseSpaces.P_of_θϕ(θ,ϕ)
    #             try
    #                 return log(probability_density([Q,NaN,P,p]))#*sin(θ)?
    #             catch
    #                 return -Inf
    #             end
    #           end
    #
    #         m = RWMVariate([PhaseSpaces.θ_of_QP(Q,P),PhaseSpaces.ϕ_of_QP(Q,P),p], [1/sqrt(j),1/sqrt(j),1/sqrt(j)],logf,proposal = Normal)
    #         function sample()
    #                 θ,ϕ,p=sample!(m)
    #                 return Point(system,Q=PhaseSpaces.Q_of_θϕ(θ,ϕ),P=PhaseSpaces.P_of_θϕ(θ,ϕ),p=p,ϵ=ϵ,sgn=signoq)
    #         end
    #         for i in 1:10000
    #             sample!(m) #bake
    #         end
    #         TWA.PhaseSpaceDistribution(probability_density,sample)
    #     end

        """
        ```julia
        function classical_path_random_sampler(system::DickeSystem; ϵ::Real, dt::Real=3)
        ```
        This function returns a function `sample()`, which produces random points within the classical
        energy shell at energy `ϵ`. The function `sample()` returns points from a fixed chaotic
        trajectory picked at random separated by a fixed time interval  `dt`. If the classical dynamics are ergodic,
        **which only happens in chaotic regions,** this produces random points in the energy shell.
        # Arguments
        - `system` should be a subtype of [`DickeSystem`](@ref ClassicalDicke.DickeSystem).

        # Keyword arguments
        - `ϵ` is the energy of the energy shell from where to sample.
        - `dt` is the fixed time interval that separates the points that are returned by `sample()`.
        """
        function classical_path_random_sampler(system::DickeSystem;ϵ::Real,dt::Real=3)
             ω₀,ω,γ=ClassicalSystems.parameters(system)
             s=sqrt(16*γ^4 + 8*γ^2*ϵ*ω + ω^2*ω₀^2)/(2*γ^2)
             a=2- ω*ω₀/(2*γ^2)
             mx=sqrt(a+s)
             if a<s
                 mn=0
             else
                 mn=sqrt(a-s)
             end
             if mn==mx
                 Q=mn #estamos en la energia del estado base
             else
                 Q=Distributions.rand(Distributions.Uniform(mn,mx))
             end
             Q=Q*rand((-1,1))
             u0=ClassicalDicke.Point(system,Q=Q,P=0,p=0,ϵ=ϵ)
             function sample()
                 u0=ClassicalSystems.integrate(system,u0,dt; save_everystep=false,show_progress=false).u[end]
                 if rand((true,false))
                     u0[1]=-u0[1]
                     u0[2]=-u0[2]
                 end
                 return u0
             end
             for i in 1:10
                 sample()
             end
             return sample
         end

end


text:
using HorizonSideRobots

function line(r,i)
    while !isborder(r,HorizonSide(i % 4)) && !isborder(r,HorizonSide((i+1) % 4))
        move!(r,HorizonSide(i % 4))
        move!(r,HorizonSide((i+1) % 4))
        putmarker!(r)
    end
    while ismarker(r)
        move!(r,HorizonSide((i+2) % 4))
        move!(r,HorizonSide((i+3) % 4))
    end
end

function direction(r)
    for i in 0:3
        line(r,i)
    end
    putmarker!(r)
end

direction(r)

max_stars_repo_name: Cheburkova/mirea-progs

max_issues_repo_name: Cheburkova/mirea-progs

max_forks_repo_name: Cheburkova/mirea-progs


text:
function adjacent(game,unit::Unit)
	hadj=adjacent(unit.loc)
	adj=Unit[]
	for h in hadj
		u=game.map[h]
		if isa(u,Unit)
			push!(adj,u)
		end
	end
	return adj
end
function resetmap(game)
	for loc in game.board.grid #keys(game.map)? No, offgrid locations are safe
		game.map[loc]=0
	end
	drawboard(game)
	reveal(game.board.c)
end

function getcellgroup(game,unit::Unit)
	white=(1,1,1)
	cellgroup=[unit]
	if unit.color==white
		adju=adjacent(game,unit)
		for u in adju
			g=getcellgroup(game,u) #can this be optimized when there are several white units?
			for gu in g
				if !in(gu,cellgroup)
					push!(cellgroup,gu)
				end
			end
		end
		return cellgroup
	end
	temp=[unit]
	while !isempty(temp)
		temp2=Unit[]
		for t in temp
			for h in adjacent(t.loc)
				if in(h,keys(game.map))
					tu=game.map[h]
					if !in(tu,cellgroup) && !in(tu,temp) && !in(tu,temp2) && isa(tu,Unit) && (tu.color==unit.color || tu.color==white)
						push!(temp2,tu)
					end
				end
			end
		end
		for t2 in temp2
			push!(cellgroup,t2)
		end
		temp=temp2
	end
	return cellgroup
end

function spreadlife!(game,unit,lifemap)
	white=(1,1,1)
	hex=unit.loc
	if !in(hex,keys(lifemap))
		return lifemap
	end
	lifemap[hex].+=unit.baselife.*unit.color
	temp=Dict(hex=>unit.baselife)
	checked=[hex]
	for rad in 1:unit.ir
		temp2=Dict()
		for t in temp
			for h in adjacent(t[1],1,unit.groundlevel)
				if !in(h,keys(temp)) && !in(h,keys(temp2)) && !in(h,checked) && in(h,game.board.grid)
					lif=(unit.baselife/6/rad).*unit.color
					if game.map[h]==0 || unit.passover
						temp2[h]=lif
					elseif unit.passoverself && (game.map[h].color==unit.color || game.map[h].color==white)
						temp2[h]=lif
					elseif unit.inclusive && game.map[h]!=0 #&& !in(h,keys(temp2))
						lifemap[h].+=lif
					end
				end
				push!(checked,h)
			end
		end
		for (h2,i2) in temp2
			lifemap[h2].+=i2
		end
		temp=temp2
	end
	return lifemap
end
function spreadlife(game,unit)
	lifemap=Dict()
	for loc in game.board.grid
		lifemap[loc]=[0.0,0,0]
	end
	return spreadlife!(game,unit,lifemap)
end
function unitslive(game,color)
	lifemap=Dict()
	for loc in game.board.grid
		lifemap[loc]=[0.0,0,0]	#this is sometimes redundant and sometimes needed...
	end
	for unit in game.units
		if (color==(1,1,1) || unit.color==color || unit.color==(1,1,1)) #white, your superpowers are being removed, enjoy them while you can
			unit.live!(game,unit,lifemap)
		end
	end
	return lifemap
end
function unitslive(game)
	lifemap=Dict()
	for loc in game.board.grid
		lifemap[loc]=[0.0,0,0]
	end
	for unit in game.units
		unit.live!(game,unit,lifemap)
	end
	return lifemap
end

function placeseq!(game)
	for entry in game.sequence
		if isa(entry,Unit)
			placeunit!(game,entry)
		elseif entry==:harvest
			harvest!(game)
		elseif entry[1]==:expand
			expandboard!(game,entry[2]...,false)
		end
	end
end

function getgroup(game,unit)
	lifemap=unitslive(game,unit.color)
	connectedunits=[unit]
	spawns=Unit[]
	body=getcellgroup(game,unit)
	for unit in body
		if unit.canspawn
			push!(spawns,unit)
		end
	end
	ulocs=[unit.loc]
	temp=[unit.loc]
	while !isempty(temp)
		temp2=Tuple[]
		for t in temp
			for h in adjacent(t)
				if !in(h,ulocs) && !in(h,temp) && !in(h,temp2) && in(h,game.board.grid) && sum(lifemap[h])>0
					push!(temp2,h)
					u=game.map[h]
					if isa(u,Unit) && !in(u,connectedunits)
						push!(connectedunits,u)
					end
				end
			end
		end
		for t2 in temp2
			push!(ulocs,t2)
		end
		temp=temp2
	end
	return newgroup(spawns,body,connectedunits)
end
function getgroups(game)
	groups=Group[]
	for spawn in game.spawns
		g=getgroup(game,spawn)
		if !hasgroup(groups,g)
			push!(groups,g)
		end
	end
	return groups
end
function sync!(game::Game)
	game.lifemap=unitslive(game)
	game.groups=getgroups(game)
	game.points=checkharvest(game,false)
	GAccessor.text(game.gui[:scorelabel],pointslabel(game))
	GAccessor.text(game.gui[:newslabel],infolabel(game))
end
function placeable(game,unit::Unit)
	if haskey(game.map,unit.loc) && game.map[unit.loc]!=0
		return false
	elseif !in(unit.loc[3],unit.pl)
		return false
	end
	return true
end
function placeunit!(game,unit)
	if placeable(game,unit)
		game.map[unit.loc]=unit
		if !in(unit,game.sequence)
			push!(game.sequence,unit)
		end
		if !in(unit,game.units)
			push!(game.units,unit)
		end
		if unit.canspawn && !in(unit,game.spawns)
			push!(game.spawns,unit)
		end
	end
	return "<3"
end
function findunit(unit::Unit,units)
	i=0
	for u in units
		i+=1
		if unit.loc==u.loc && unit.color==u.color && unit.name==u.name
			return i
		end
	end
	return 0
end
function removeunit!(game,unit::Unit)
	game.map[unit.loc]=0
	push!(game.sequence,(:delete,unit))
	i=findunit(unit,game.units)
	if i>0
		deleteat!(game.units,i)
	end
	if unit.canspawn
		i=findunit(unit,game.spawns)
		deleteat!(game.spawns,i)
	end
	return "<3"
end

function newledger(game)
	ledger=Dict()
	for loc in game.board.grid
		ledger[loc]=[[0.0,0,0],0.0,0,0,0]
	end
	return ledger
end
function getpoints!(game,unit,loc,distance,ledger,partial=1) #remake cleaner. Meaning, more humanreadable? ll lci llm l lci lif. OMG this place is bugprone, readability helps.
	llm=game.lifemap[loc] #local lifemap
	ll=ledger[loc]
	lif=distance==0 ? unit.baselife : (unit.baselife/6/distance)
	lif=lif*partial
	ncol=numcolors(llm)
	l=llm.-ll[5] #local life - lightharvest
	ncoll=numcolors(l)
	points=[[0.0,0,0],0.0,0,0,0]
	if ncol==3 && ncoll==3
		points[5]=min(min(l...),lif)
		#l.-=points[5]
		ll.+=points
		lif-=points[5]
		if numcolors(l)==1
			l[1]=0;l[2]=0;l[3]=0
		end
		l=llm.-ll[5]
	end
	if ncol==1 && game.map[loc]==0 && sum(ll[1])<1
		ci=l[1]==0 ? (l[2]==0 ? 3 : 2) : 1
		if l[ci]>1;l[ci]=1;end
		h=min(l[ci],lif,1)
		points[1]=points[1].+h.*unit.color
		#l[ci]-=h
		ll[1]=ll[1]+points[1]
#if unit.color==(0,0,1);println(loc,ll,l);end #useful for debugging
	else
		#l=l-ll[2:4] #something is wrong here... Aha! Green doesn't harvest green, green harvests blue!
		l[1]-=ll[4]
		l[2]-=ll[2]
		l[3]-=ll[3]
		lifc=lif.*unit.color #nooo? Yes?
		for c in 1:3
			if lifc[c]>0
				points[c+1]+=min(l[c],l[c%3+1],lifc[c]) #no! Maybe...
				#game.lifemap[loc][c%3+1]-=points[c+1] #problem? *trollface* Solved!
				ll[c+1]+=points[c+1]
			end
		end
	end
	return points
end
function unitharvest(game,unit,ledger,partial=1)
	points=[[0.0,0,0],0.0,0,0,0]
	if unit.harvested>=1
		return points
	end

	white=(1,1,1) #remove?
	points.+=getpoints!(game,unit,unit.loc,0,ledger)
	temp=[unit.loc]
	checked=[unit.loc]
	for rad in 1:unit.ir
		temp2=[]
		for t in temp
			for h in adjacent(t,1,unit.groundlevel)
				if !in(h,temp) &&!in(h,checked) && in(h,game.board.grid)
					if game.map[h]==0 || unit.passover #why not merge these clauses?
						p=getpoints!(game,unit,h,rad,ledger,partial)
						points.+=p
						push!(temp2,h)
					elseif unit.passoverself && (game.map[h].color==unit.color || game.map[h].color==white)
						points.+=getpoints!(game,unit,h,rad,ledger) #problem! Not anymore! Why are the comments still here? Why not? I asked first! They liven up the code. OK then
						push!(temp2,h)
					elseif unit.inclusive && game.map[h]!=0
						points.+=getpoints!(game,unit,h,rad,ledger) #*seriousface* solving by deprecating game.lifemap. No! Long live the lifemap+ledger union
					end
				end
				push!(checked,h)
			end
		end
		temp=temp2
	end
	return points
end

function checkharvest(game,unit::Unit,ledger,partial=1)	#| ^
	return unit.harvest(game,unit,ledger,partial)	#L_|
end
function checkharvest(game,group::Group,ledger)
	points=[[0.0,0,0],0.0,0,0,0]
	for unit in group.units
		p=checkharvest(game,unit,ledger)
		if !in(unit,group.body)
			points.+=p./3
		else
			points.+=p
		end
	end
	return points
end
function checkharvest(game::Game,sync::Bool=true)
	#game=deepcopy(game) #this is silly, remove ! from unit.harvest! Done, still updates lifemap thou. Not anymore, now puts harvest data in a ledger
	if sync
		sync!(game) #this shouldn't count as modifying the state of the game since if they arent correct the state is incorrect and should always be updated. Maybe we don't have to sync thou since it is always synced upon new move... Robustness over velocity!
	end
	ledger=newledger(game)
	points=[[0.0,0,0],0.0,0,0,0]
	for group in game.groups
		#points+=checkharvest(game,group,ledger)
# !
		for unit in group.body
			points+=checkharvest(game,unit,ledger)
			unit.harvested=true
		end
	end
	for group in game.groups
# !
		for unit in group.units
			if float(unit.harvested)<1
				partial=min(1/3,1-unit.harvested)
				points+=checkharvest(game,unit,ledger,partial)
				unit.harvested+=partial
			end
		end
		group.harvested=true
	end
	for group in game.groups
		group.harvested=false
		for unit in group.units
			unit.harvested=false
		end
	end
#may need to do this to avoid units in multiple groups to multiharvest, but can't do it now because it should harvest to the group where it is in the body, so have to harvest with bodies first. Or let it harvest partially
	return points
end

function bonds(game)
	nc=0
	#checked=[]
	for (loc,unit) in game.map #rewrite with game.units?
	#	push!(checked,loc) #this didnt work...
		if isa(unit,Unit)
			for c in adjacent(unit.loc)
				if in(c,keys(game.map))
					ac=game.map[c]
					if isa(ac,Unit) && ac.color!=unit.color
						nc+=1
					end
				end
			end
		end
	end
	return Int(nc/2)
end
function pointslabel(game,sync::Bool=true)
	points=game.points #checkharvest(game,sync)
	bp=round.(points[1];digits=1)
	points[1]=0
	points=round.(points;digits=1)
	return "Points!\nLite:\nRed\t$(bp[1])\nGreen\t$(bp[2])\nBlue\t$(bp[3])\nLife:\nRed\t$(points[2])\nGreen\t$(points[3])\nBlue\t$(points[4])\nLight:\t$(points[5])"
end
function infolabel(game)
	rgb=[0.0,0,0]
	for unit in game.units
		rgb.+=unit.color
	end
	rgb=round.(rgb;digits=1)
	return "Information!\nUnits: $(length(game.units))\nRed: $(rgb[1])\nGreen: $(rgb[2])\nBlue: $(rgb[3])\nBonds: $(bonds(game))"
end
function undo!(game) #wont undo captures? Maybe when reloading sequence. There aren't captures anymore. Wont undo board expansions. Maybe it should? Not much use now... Easy right, just pop the seq and reload
	removeunit!(game.units[end])
	if !game.colock
		game.colind-=1
		if game.colind<1
			game.colind=game.colmax
		end
	end
	return hex
end

function drawboard(game,ctx,w,h)
	game.board.size=game.window[2]/(game.board.shells*game.board.sizemod)
	size=game.board.size
	rectangle(ctx, 0, 0, w, h)
	set_source_rgb(ctx, game.board.bgcolor...)
	fill(ctx)
	set_source_rgb(ctx, game.color...)
	arc(ctx, size, size, 3size, 0, 2pi)
	fill(ctx)
	offx=game.board.offsetx+game.board.panx
	offy=game.board.offsety+game.board.pany
	for (lo,lif) in game.lifemap
		if sum(lif)>0 && lo[3]==2
			offset=(offx,offy)
			#if lo[3]==1
			#	offset=offset.+(-cos(pi/6)*size,sin(pi/6)*size)
			#elseif lo[3]==3
			#	offset=offset.+(-cos(pi/6)*size,-sin(pi/6)*size)
			#end
			plo=hex_to_pixel(lo[1],lo[2],size)
			ploc=(plo[1]+offset[1]+w/2,plo[2]+offset[2]+h/2)
			rad=size*0.866/2
			col=lif./(lif.+3)
			set_source_rgb(ctx, col...)
			arc(ctx,ploc[1],ploc[2],rad*2, 0, 2pi)
			fill(ctx)
		end
	end
	set_source_rgb(ctx, game.board.gridcolor...)
	for loc in game.board.grid
		if loc[3]==2
			x,y=hex_to_pixel(loc[1],loc[2],size)
			hexlines(ctx,x+w/2+offx,y+h/2+offy,size)
		end
	end
	for unit in game.units
		offset=(offx,offy)
		if unit.loc[3]==1
			offset=offset.+(-cos(pi/6)*size,sin(pi/6)*size)
		elseif unit.loc[3]==3
			offset=offset.+(-cos(pi/6)*size,-sin(pi/6)*size)
		end
		loc=hex_to_pixel(unit.loc[1],unit.loc[2],size)
		floc=(loc[1]+offset[1]+w/2,loc[2]+offset[2]+h/2)
		rad=size*0.866/2
		#unit border:
		set_source_rgb(ctx,game.board.gridcolor...)
		arc(ctx, floc[1],floc[2],rad+1, 0, 2pi)
		stroke(ctx)
		set_source_rgb(ctx,unit.color...)
		arc(ctx,floc[1],floc[2],rad, 0, 2pi) #why isn't the circle radius the distance between locs? Whyyyy whyyyy someone pleeease fiiix
		fill(ctx)
		if !isempty(unit.graphic)
			set_source_rgb(ctx,game.board.bgcolor...)
			points=Point[]
			for p in unit.graphic
				push!(points,Point(floc[1]+rad*p[1],floc[2]+rad*p[2]))
			end
			polygon(ctx,points)
			fill(ctx)
		end
	end
	#showall(game.board.win) #should probably look up the difference between all these revealing methods
	reveal(game.board.c)
end
function drawboard(game::Game)
	ctx=getgc(game.board.c)
	h=height(game.board.c)
	w=width(game.board.c)
	drawboard(game,ctx,w,h)
end

function expandboard!(game::Game,shells::Integer=6,initlocs=[(6,6,2)],reveal=true)
	patch=makegrid(shells,initlocs)
	for loc in patch
		if !in(loc,keys(game.map))
			game.map[loc]=0
			push!(game.board.grid,loc)
		end
	end
	push!(game.sequence,(:expand,[shells,initlocs]))
	if reveal #without this there is sometimes some severe error (when loadicing)
		drawboard(game)
	end
	return "<3"
end
function zoom(game,factor)
	game.board.sizemod*=factor
	drawboard(game)
end
function pass!(game,nomax::Bool=false,reverse::Bool=false)
	max=game.colmax
	if nomax
		max=length(game.colors)
	end
	if reverse
		game.colind=game.colind-1
		if game.colind==0
			game.colind=max
		end
	else
		game.colind=game.colind%max+1
	end
	game.color=game.colors[game.colind]
	drawboard(game)
end

function center(game,hex)
	#game.board.sizemod=Gtk.G_.value(game.gui[:zadj])/10
	game.board.size=game.window[2]/(game.board.shells*game.board.sizemod)
	loc=hex_to_pixel(hex[1],hex[2],game.board.size)
	game.board.offsetx=-loc[1]+Gtk.G_.value(game.gui[:xadj])
	game.board.offsety=-loc[2]+Gtk.G_.value(game.gui[:yadj])
	drawboard(game)
	return true
end


text:
module GraphQLGen

using RBNF
using Expronicon

include("types.jl")
include("parser.jl")
include("dag.jl")
include("codegen.jl")

include("ExprPrettify/expr_prettify.jl")
include("print.jl")

include("generate.jl")

"""
    parse(schema::AbstractString)

Parse a GraphQL schema at top-level to AST.
"""
function parse(schema::AbstractString)
    tokens = RBNF.runlexer(GQL, schema)
    ast, ctx = RBNF.runparser(document, tokens)
    ctx.tokens.current > ctx.tokens.length ||
        throw(Meta.ParseError("invalid syntax in GraphQL program"))
    return ast
end

end


text:
# EXCLUDE FROM TESTING

using CuArrays, CUDAnative, CUDAdrv, CUDAapi

CUDAnative.initialize()
const dev = device()
const cap = capability(dev)

using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10
BenchmarkTools.DEFAULT_PARAMETERS.gcsample = true

using SpecialFunctions

## scalar CPU version

@inline cndf2(in::Float32) = 0.5f0 + 0.5f0 * erf(0.707106781f0 * in)

function blackscholes_cpu(sptprice::Float32, strike::Float32, rate::Float32,
                      volatility::Float32, time::Float32)
    logterm = log10(sptprice / strike)
    powterm = .5f0 * volatility * volatility
    den = volatility * sqrt(time)
    d1 = (((rate + powterm) * time) + logterm) / den
    d2 = d1 - den
    NofXd1 = cndf2(d1)
    NofXd2 = cndf2(d2)
    futureValue = strike * exp(-rate * time)
    c1 = futureValue * NofXd2
    call = sptprice * NofXd1 - c1
    return call - futureValue + sptprice
end

## vectorized CPU version

@inline cndf2(in::AbstractArray{Float32}) = 0.5f0 .+ 0.5f0 .* erf.(0.707106781f0 .* in)

function blackscholes_cpu(sptprice::AbstractArray{Float32},
                      strike::AbstractArray{Float32},
                      rate::AbstractArray{Float32},
                      volatility::AbstractArray{Float32},
                      time::AbstractArray{Float32})
    logterm = log10.(sptprice ./ strike)
    powterm = .5f0 .* volatility .* volatility
    den = volatility .* sqrt.(time)
    d1 = (((rate .+ powterm) .* time) .+ logterm) ./ den
    d2 = d1 .- den
    NofXd1 = cndf2(d1)
    NofXd2 = cndf2(d2)
    futureValue = strike .* exp.(- rate .* time)
    c1 = futureValue .* NofXd2
    call = sptprice .* NofXd1 .- c1
    return call .- futureValue .+ sptprice
end

## native CUDA version

@inline cndf2_cuda(in::Float32) = 0.5f0 + 0.5f0 * CUDAnative.erf(0.707106781f0 * in)

function blackscholes_kernel(sptprice::AbstractArray{Float32},
                             strike::AbstractArray{Float32},
                             rate::AbstractArray{Float32},
                             volatility::AbstractArray{Float32},
                             time::AbstractArray{Float32},
                             out::AbstractArray{Float32})
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x

    if i <= size(sptprice, 1)
        logterm = CUDAnative.log10(sptprice[i] / strike[i])
        powterm = 0.5f0 * volatility[i] * volatility[i]
        den = volatility[i] * CUDAnative.sqrt(time[i])
        d1 = (((rate[i] + powterm) * time[i]) + logterm) / den
        d2 = d1 - den
        NofXd1 = cndf2_cuda(d1)
        NofXd2 = cndf2_cuda(d2)
        futureValue = strike[i] * CUDAnative.exp(-rate[i] * time[i])
        c1 = futureValue * NofXd2
        call = sptprice[i] * NofXd1 - c1
        out[i] = call - futureValue + sptprice[i]
    end

    return
end

## scalar CuArrays version

function blackscholes_gpu(sptprice::Float32, strike::Float32, rate::Float32,
                          volatility::Float32, time::Float32)
    logterm = CUDAnative.log10(sptprice / strike)
    powterm = .5f0 * volatility * volatility
    den = volatility * CUDAnative.sqrt(time)
    d1 = (((rate + powterm) * time) + logterm) / den
    d2 = d1 - den
    NofXd1 = cndf2_cuda(d1)
    NofXd2 = cndf2_cuda(d2)
    futureValue = strike * CUDAnative.exp(-rate * time)
    c1 = futureValue * NofXd2
    call = sptprice * NofXd1 - c1
    return call - futureValue + sptprice
end

## vectorized CuArrays version

@inline cndf2_cuarr(in::AbstractArray{Float32}) = 0.5f0 .+ 0.5f0 .* CUDAnative.erf.(0.707106781f0 .* in)

function blackscholes_gpu(sptprice::AbstractArray{Float32},
                          strike::AbstractArray{Float32},
                          rate::AbstractArray{Float32},
                          volatility::AbstractArray{Float32},
                          time::AbstractArray{Float32})
    logterm = CUDAnative.log10.(sptprice ./ strike)
    powterm = .5f0 .* volatility .* volatility
    den = volatility .* CUDAnative.sqrt.(time)
    d1 = (((rate .+ powterm) .* time) .+ logterm) ./ den
    d2 = d1 .- den
    NofXd1 = cndf2_cuarr(d1)
    NofXd2 = cndf2_cuarr(d2)
    futureValue = strike .* CUDAnative.exp.(- rate .* time)
    c1 = futureValue .* NofXd2
    call = sptprice .* NofXd1 .- c1
    return call .- futureValue .+ sptprice
end

## non-native CUDA C version

const cuda_source = "$(tempname()).cu"
const cuda_ptx = "$(tempname()).ptx"

open(cuda_source, "w") do io
    print(io, """
        extern "C" __global__ void blackscholes_kernel(const float *sptprice,
                                                       const float *strike,
                                                       const float *rate,
                                                       const float *volatility,
                                                       const float *time,
                                                       float *out,
                                                       size_t n)
        {
            int i = blockIdx.x * blockDim.x + threadIdx.x;
            if (i < n) {
                float logterm = log10(sptprice[i] / strike[i]);
                float powterm = 0.5 * volatility[i] * volatility[i];
                float den = volatility[i] * sqrt(time[i]);
                float d1 = (((rate[i] + powterm) * time[i]) + logterm) / den;
                float d2 = d1 - den;
                float NofXd1 = 0.5 + 0.5 * erf(0.707106781 * d1);
                float NofXd2 = 0.5 + 0.5 * erf(0.707106781 * d2);
                float futureValue = strike[i] * exp(-rate[i] * time[i]);
                float c1 = futureValue * NofXd2;
                float call = sptprice[i] * NofXd1 - c1;
                out[i] = call - futureValue + sptprice[i];
            }
        }
    """)
end

toolkit = CUDAapi.find_toolkit()
nvcc = CUDAapi.find_cuda_binary("nvcc", toolkit)
toolchain = CUDAapi.find_toolchain(toolkit)
flags = `-ccbin=$(toolchain.host_compiler) -arch=sm_$(cap.major)$(cap.minor)`
run(`$nvcc $flags -ptx -o $cuda_ptx $cuda_source`)

const cuda_module = CuModuleFile(cuda_ptx)
const cuda_function = CuFunction(cuda_module, "blackscholes_kernel")

## main

function checksum(reference, result)
    reference_sum = sum(reference)
    result_sum = sum(result)
    diff = abs(1-reference_sum/result_sum)
    if diff>0.01
        warn("checksum failed: $result_sum instead of $reference_sum (relative difference: $diff)")
        println(stacktrace())
    end
end

function main(iterations)
    sptprice   = Float32[ 42.0 for i = 1:iterations ]
    strike     = Float32[ 40.0 + (i / iterations) for i = 1:iterations ]
    rate       = Float32[ 0.5 for i = 1:iterations ]
    volatility = Float32[ 0.2 for i = 1:iterations ]
    time       = Float32[ 0.5 for i = 1:iterations ]

    timings = Dict()

    reference = blackscholes_cpu.(sptprice, strike, rate, volatility, time)

    let benchmark = @benchmarkable begin
                out = blackscholes_cpu.($sptprice, $strike, $rate,
                                        $volatility, $time)
            end setup=(
                out = zeros($sptprice)
            ) teardown=(
                checksum($reference, out)
            )
        timings["Single-threaded (scalar)"] = run(benchmark)
    end

    let benchmark = @benchmarkable begin
                out = blackscholes_cpu($sptprice, $strike, $rate,
                                       $volatility, $time)
            end setup=(
                out = zeros($sptprice)
            ) teardown=(
                checksum($reference, out)
            )
        timings["Single-threaded (vectorized)"] = run(benchmark)
    end

    let benchmark = @benchmarkable begin
                cudacall(cuda_function,
                         Tuple{Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat},
                               Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Csize_t},
                         sptprice_dev, strike_dev, rate_dev, volatility_dev,
                         time_dev, out, n; blocks=grid, threads=block)
                synchronize()
            end setup=(
                sptprice_dev = CUDAdrv.CuArray($sptprice);
                strike_dev = CUDAdrv.CuArray($strike);
                rate_dev = CUDAdrv.CuArray($rate);
                volatility_dev = CUDAdrv.CuArray($volatility);
                time_dev = CUDAdrv.CuArray($time);
                out = CUDAdrv.CuArray{Float32}(size($strike));

                n = size($sptprice, 1);
                block = min(n, 1024);
                grid = ceil(Integer, n/block)
            ) teardown=(
                checksum($reference, Array(out))
            )
        timings["CUDA C (kernel)"] = run(benchmark)
    end

    let benchmark = @benchmarkable begin
                @cuda blocks=grid threads=block blackscholes_kernel(sptprice_dev, strike_dev, rate_dev,
                                                                    volatility_dev, time_dev, out)
                synchronize()
            end setup=(
                sptprice_dev = CUDAdrv.CuArray($sptprice);
                strike_dev = CUDAdrv.CuArray($strike);
                rate_dev = CUDAdrv.CuArray($rate);
                volatility_dev = CUDAdrv.CuArray($volatility);
                time_dev = CUDAdrv.CuArray($time);
                out = CUDAdrv.CuArray{Float32}(size($strike));

                n = size($sptprice, 1);
                block = min(n, 1024);
                grid = ceil(Integer, n/block)
            ) teardown=(
                checksum($reference, Array(out))
            )
        timings["CUDAnative.jl (kernel)"] = run(benchmark)
    end

    let benchmark = @benchmarkable begin
                out = blackscholes_gpu.(sptprice_dev, strike_dev, rate_dev,
                                        volatility_dev, time_dev)
                synchronize()
            end setup=(
                sptprice_dev = CuArrays.CuArray($sptprice);
                strike_dev = CuArrays.CuArray($strike);
                rate_dev = CuArrays.CuArray($rate);
                volatility_dev = CuArrays.CuArray($volatility);
                time_dev = CuArrays.CuArray($time);
                out = zeros($sptprice)
            ) teardown=(
                checksum($reference, Array(out))
            )
        timings["CuArrays.jl (scalar)"] = run(benchmark)
    end

    let benchmark = @benchmarkable begin
                out = blackscholes_gpu(sptprice_dev, strike_dev, rate_dev,
                                       volatility_dev, time_dev)
                synchronize()
            end setup=(
                sptprice_dev = CuArrays.CuArray($sptprice);
                strike_dev = CuArrays.CuArray($strike);
                rate_dev = CuArrays.CuArray($rate);
                volatility_dev = CuArrays.CuArray($volatility);
                time_dev = CuArrays.CuArray($time);
                out = zeros($sptprice)
            ) teardown=(
                checksum($reference, Array(out))
            )
        timings["CuArrays.jl (vectorized)"] = run(benchmark)
    end

    return timings
end

function driver()
    iterations = 10^7
    timings = main(iterations)

    println()
    println("Timings:")
    for (test, trials) in timings
        println("* $test: ", BenchmarkTools.prettytime(time(trials)))
    end

    println()
    println("Rates:")
    for (test, trials) in timings
        println("* $test: ", 1e9*iterations/time(trials), " ops/sec")
    end
end

driver()

rm(cuda_source)
rm(cuda_ptx)


text:
module LoggerTests

using Logging
using QXContexts.Logger
using Test
using Dates
using MPI

@testset "Logger Tests" begin
    @testset "QXLogger" begin
        @testset "INFO test" begin
            io = IOBuffer()
            global_logger(QXLogger(io; show_info_source=true))
            @info "info_test"

            log = split(String(take!(io)), "\n")[1:end-1]
            df = DateFormat("[yyyy/mm/dd-HH:MM:SS.sss]");

            for l in log
                log_elem = split(l, " ")
                em = collect(eachmatch(r"\[(.*?)\]", log_elem[1])) # capture values in []
                @test DateTime(em[1].match, df) !== nothing
                @test match(r"[rank=/\d+/]", em[2].match) !== nothing
                @test em[3].match == "[host=$(gethostname())]"
                @test log_elem[2] == "INFO"
                @test log_elem[3] == "info_test"
            end
        end

        @testset "WARN test" begin
            io = IOBuffer()
            global_logger(QXLogger(io; show_info_source=true))

            @warn "warn_test"; line_num = @__LINE__

            log = split(String(take!(io)), "\n")[1:end-1]
            df = DateFormat("[yyyy/mm/dd-HH:MM:SS.sss]");

            for l in log
                log_elem = split(l, " ")
                em = collect(eachmatch(r"\[(.*?)\]", log_elem[1])) # capture values in []
                @test DateTime(em[1].match, df) !== nothing
                @test match(r"[rank=/\d+/]", em[2].match) !== nothing
                @test em[3].match == "[host=$(gethostname())]"
                @test log_elem[2] == "WARN"
                @test log_elem[3] == "warn_test"
                @test log_elem[4] == "-@->"
                file_line = splitext(log_elem[5])
                @test file_line[1] == splitext(@__FILE__)[1]
                @test split(file_line[2], ":")[2] == string(line_num)
            end
        end

        @testset "ERROR test" begin
            io = IOBuffer()
            global_logger(QXLogger(io; show_info_source=true))

            line_num = @__LINE__; @error "error_test"

            log = split(String(take!(io)), "\n")[1:end-1]
            df = DateFormat("[yyyy/mm/dd-HH:MM:SS.sss]");

            for l in log
                log_elem = split(l, " ")
                em = collect(eachmatch(r"\[(.*?)\]", log_elem[1])) # capture values in []
                @test DateTime(em[1].match, df) !== nothing
                @test match(r"[rank=/\d+/]", em[2].match) !== nothing
                @test em[3].match == "[host=$(gethostname())]"
                @test log_elem[2] == "ERROR"
                @test log_elem[3] == "error_test"
                @test log_elem[4] == "-@->"
                file_line = splitext(log_elem[5])
                @test file_line[1] == splitext(@__FILE__)[1]
                @test split(file_line[2], ":")[2] == string(line_num)
            end
        end
    end

    @testset "QXLoggerMPIPerRank" begin

        @testset "WARN test" begin

            mktempdir() do path
                # Fail to create logger if  MPI not initialised
                if !MPI.Initialized()
                    @test_throws """MPI is required for this logger. Pleasure ensure MPI is initialised. Use `QXLogger` for non-distributed logging""" QXLoggerMPIPerRank()
                    MPI.Init()
                end

                global_logger(QXLoggerMPIPerRank(; show_info_source=true, path=path))

                line_num = @__LINE__; @warn "warn_test"

                df = DateFormat("[yyyy/mm/dd-HH:MM:SS.sss]");
                @test isdir(joinpath(path, "QXContexts_io_" * string(global_logger().session_id)))

                log = readlines(joinpath(path, "QXContexts_io_" * string(global_logger().session_id), "rank_0.log"))

                for l in log
                    log_elem = split(l, " ")
                    em = collect(eachmatch(r"\[(.*?)\]", log_elem[1])) # capture values in []
                    @test DateTime(em[1].match, df) !== nothing
                    @test match(r"[rank=/\d+/]", em[2].match) !== nothing
                    @test em[3].match == "[host=$(gethostname())]"
                    @test log_elem[2] == "WARN"
                    @test log_elem[3] == "warn_test"
                    @test log_elem[4] == "-@->"
                    file_line = splitext(log_elem[5])
                    @test file_line[1] == splitext(@__FILE__)[1]
                    @test split(file_line[2], ":")[2] == string(line_num)
                end
            end
        end

    end

end

end


text:
module ReconstructionEvaluations

using MATLAB
using HDF5
using StatsBase
using PyPlot
using Graphics
using FastConv

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

export
    load_edges,
    edges_to_syn_dicts,
    segm_overlap, overlap_in_chunks, create_segID_map,
    map_synapses,
    get_indexed_seg_IDs,
    build_count_table,
    compute_nri,
    merge_columns,
    split_column,
    remove_synapse,
    add_synapse,
    merge_sensitivity,
    split_sensitivity,
    NRI,
    #synapse recovery estimates
    find_trunk_mapping, synapse_recovery,
    Vec3, Point3, BoundingCube,
    # limits in world coordinates
    is_inside, xmin, xmax, ymin, ymax, zmin, zmax, center,
    xrange, yrange, zrange, shift, deform,
    width, height, depth,
    # clustering methods
    read_sparse, write_sparse,
    load_seg_sizes, classify_pre_post,
    hist_seg_sizes

include("utils/include.jl")

include("eval/include.jl")

include("analysis/include.jl")
using .SBMs

include("CronUtils/include.jl")

include("experiments/experiments.jl")

end


text:
"""
    process_streaks(df::AbstractDataFrame; photometry = false)

Take as input a processed pokes datafame `df` of 1 single session and returns a dataframe were each row are information about a poking streak.
A poking streak ends when the animal changes side. For photometry data a set of other columns regarding indexes for traces allignment are computed.
"""

function process_streaks(df::AbstractDataFrame; photometry = false)
    vals_per_streak = [:Stim, :StimFreq, :Wall, :Protocol, :Block, :StreakInBlock, :Side, :ReverseStreak]
    vals_per_day = [:Box,:StimDay];
    booleans=[:Reward,:Stim,:Wall,:Correct,:StimDay]#columns to convert to Bool
    for x in booleans
        if x in propertynames(df)
            df[:,x] = eltype(df[:,x]) == Bool ? df[:,x] : occursin.("true",df[:,x],)
        end
    end
    streak_table = combine(groupby(df, :Streak)) do dd
        lr_idx = findlast(dd[!,:Reward] .== 1).== nothing ? 0 : findlast(dd[!,:Reward] .== 1)
        dt = DataFrame(
        Num_pokes = size(dd,1),
        Num_Rewards = length(findall(dd[!,:Reward].==1)),
        Start_Reward = dd[1,:Reward],
        Last_Reward = lr_idx,
        Last_Reward_Time = lr_idx == 0 ? dd[1,:PokeOut] -  dd[1,:PokeIn] : dd[lr_idx,:PokeOut] -  dd[1,:PokeIn],
        Prev_Reward = findlast(dd[!,:Reward] .== 1).== nothing ? 0 : findprev(dd[!,:Reward] .==1, findlast(dd[!,:Reward] .==1)-1),
        Trial_Duration = (dd[end,:PokeOut]-dd[1,:PokeIn]),
        Poking_time = sum(dd.PokeDur),
        Start = (dd[1,:PokeIn]),
        Stop = (dd[end,:PokeOut]),
        PreInterpoke = size(dd,1) > 1 ? maximum(skipmissing(dd[!,:PreInterpoke])) : missing,
        PostInterpoke = size(dd,1) > 1 ? maximum(skipmissing(dd[!,:PostInterpoke])) : missing,
        CorrectStart = dd[1,:Correct],
        CorrectLeave = !dd[end,:Correct]
        )

        for v in vals_per_streak
            if v in propertynames(dd)
                dt[!,v] .= dd[1, v]
            end
        end

        return dt
    end
    streak_table[!,:Prev_Reward] = [x .== nothing ? 0 : x for x in streak_table[:,:Prev_Reward]]
    streak_table[!,:AfterLast] = streak_table[!,:Num_pokes] .- streak_table[!,:Last_Reward];
    streak_table[!,:AfterLast_Duration] = streak_table[!,:Trial_Duration] .- streak_table[!,:Last_Reward_Time];
    streak_table[!,:BeforeLast] = streak_table[!,:Last_Reward] .- streak_table[!,:Prev_Reward].-1;
    prov = lead(streak_table[!,:Start],default = 0.0) .- streak_table[!,:Stop];
    streak_table[!,:Travel_to] = [x.< 0 ? 0 : x for x in prov]

    for v in vals_per_day
        if v in propertynames(df)
            streak_table[!,v] .= df[1, v]
        end
    end

    # for photometry datas a set of indexes are extracted to allign the data in the pokes dataframe
    # This indexes are to be reported in the streaks dataframe as well
    if photometry
        frames = combine(groupby(df, :Streak)) do dd
            dd = DataFrame(
            In = dd[1,:In],
            Out = dd[end,:Out],
            LR_In = findlast(dd[:Reward])==0 ? NaN : dd[findlast(dd[:Reward]),:In],
            LR_Out = findlast(dd[:Reward])==0 ? NaN : dd[findlast(dd[:Reward]),:Out]
            )
        end
        streak_table[:In] = frames[:In]
        streak_table[:Out] = frames[:Out]
        streak_table[:LR_In] = frames[:LR_In]
        streak_table[:LR_Out] = frames[:LR_Out]
    end

    return streak_table
end


text:
"""
An interval type tailored specifically to store intervals of
indices of an NDSparse object. Some of the operations on this
like `in` or `<` may be controversial for a generic Interval type.
"""
struct Interval{T}
    first::T
    last::T
end

# desired properties:
Base.eltype(int::Interval{T}) where {T} = T
Base.first(int::Interval) = int.first
Base.last(int::Interval) = int.last
Base.isempty(int::Interval) = first(int) > last(int)
Base.isless(x::Interval, y::Interval) = x.last < y.first
Base.in(x, int::Interval) = first(int) <= x <= last(int)
Base.in(x::AbstractRange, int::Interval) = hasoverlap(Interval(first(x),last(x)), int)
Base.in(x::AbstractArray, int::Interval) = any(a in int for a in x)
Base.in(x::Colon, int::Interval) = true
Base.in(int::Interval, x::Union{AbstractArray, Colon}) = x in int

function hasoverlap(i1::Interval, i2::Interval)
    (isempty(i1) || isempty(i2)) && return false

    (first(i2) <= last(i1) && first(i1) <= last(i2)) ||
    (first(i1) <= last(i2) && first(i2) <= last(i1))
end

boxintervals(i) = map(Interval, first(i), last(i))

function boxhasoverlap(a,b)
    all(map(hasoverlap, boxintervals(a), boxintervals(b)))
end

function boxmerge(a, b)
    c = map(merge, boxintervals(a), boxintervals(b))
    Interval(map(first, c), map(last, c))
end

function Base.intersect(i1::Interval, i2::Interval)
    Interval(max(first(i1), first(i2)), min(last(i1), last(i2)))
end
function Base.merge(i1::Interval, i2::Interval)
    Interval(min(first(i1), first(i2)), max(last(i1), last(i2)))
end

function _map(f, i::Interval)
    fst, lst = f(first(i)), f(last(i))
    if lst < fst
        throw(ArgumentError(
            "map on $(typeof(i)) is only allowed on monotonically increasing functions"
           )
        )
    end
    Interval(fst, lst)
end

Base.show(io::IO, i::Interval) = (show(io, i.first); print(io, ".."); show(io, i.last))


text:
# SPD-License-Identifier: MIT

using ArgParse
using BayesianModels
using CUDA
using HDF5
using JLD2
using LinearAlgebra
using Random
using TOML

CUDA.allowscalar(false)

include("utils.jl")

const DEFAULTS_MODEL = Dict(
    "ncomponents" => 1,
    "covtype" => "full",
    "priortype" => "categorical",
    "noise_init" => 0.1,
    "pstrength" => 1
)

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--checkpoint-rate", "-c"
            default = 100
            arg_type = Int
            help = "save a checkpoint of the model every C updates"
        "--learning-rate", "-l"
            default = 0.1
            arg_type = Float64
            help = "constant learning rate of the training"
        "--mini-batch-size", "-m"
            default = 1
            arg_type = Int
            help = "number of utterances per mini-batch"
        "--single-precision", "-S"
            action = :store_true
            help = "use float32 data type for computations"
        "--start-from", "-s"
            default = ""
            help = "start training from the given checkpoint model"
        "--start-update", "-k"
            default = 1
            arg_type = Int
            help = "start training the given update number"
        "--steps", "-n"
            default = 1
            arg_type = Int
            help = "number of training steps"
        "--use-gpu", "-g"
            action = :store_true
            help = "use a GPU for training the model"
        "config"
            required = true
            help = "GMM configuration file"
        "stats"
            required = true
            help = "statistics of the data to initialize the model"
        "features"
            required = true
            help = "HDF5 archive containing the features"
        "outdir"
            required = true
            help = "directory"
    end
    s.description = """
    Train a Gaussian Mixture Model on a set of utterances. The list
    of utterances to use is read from the standard input. Upon
    completion, The final model will be stored in outdir/final.jld2.
    """
    parse_args(s)
end

function main(args)
    config = TOML.parsefile(args["config"])
    outdir = args["outdir"]
    T = args["single-precision"] ? Float32 : Float64

    utts = collect(eachline(stdin))

    model = nothing
    modelconf = nothing
    if isempty(args["start-from"])
        modelconf = merge(DEFAULTS_MODEL, config["model"])
        f = h5open(args["stats"], "r")
        μ = read(f["mean"])
        Σ = read(f["cov"])
        close(f)
        model = makemodel(T, modelconf, μ, Σ)
    else
        dict = load(args["start-from"])
        model, modelconf = loadmodel(args["start-from"])
    end

    if args["use-gpu"]
        init_gpu()
        model |> gpu!
    end

    checkrate = args["checkpoint-rate"]
    lrate = args["learning-rate"]
    steps = args["steps"]
    mb_size = args["mini-batch-size"]
    inferenceconf = Dict(
        "learning-rate" => lrate,
        "steps" => steps,
        "mb_size" => mb_size
    )
    params = filter(isbayesianparam, getparams(model))
    cache = Dict()
    h5open(args["features"], "r") do ffea
        for step in args["start-update"]:steps
            # Clear cache for gradient computation
            empty!(cache)

            # Randomize the utterance list
            shuffle!(utts)[1:mb_size]

            # Load the minibatch data
            X = hcat([read(ffea[utts[n]]) for n in 1:mb_size]...)
            X = convert(Array{T}, X)
            if args["use-gpu"] X = X |> CuArray end

            # One gradient step
            scale = length(utts) / mb_size
            𝓛, llh, KL = elbo(model, X; cache = cache, stats_scale = scale,
                              detailed = true)
            ∇ξ = ∇elbo(model, cache, params)
            gradstep(∇ξ, lrate = lrate)

            # Log the status of the training in JSON format
            norm = scale * size(X,2)
            norm𝓛 = 𝓛 / norm
            normllh = llh / norm
            normKL = KL / norm
            println("{\"step\": $step, \"nsteps\":$steps, \"elbo\": $norm𝓛, \"llh\":$normllh, \"kl\":$normKL}")

            # Checkpoint
            if step % checkrate == 0
                model |> cpu!
                savemodel(outdir, model, modelconf, inferenceconf, step)
                if args["use-gpu"] model |> gpu!  end
                params = filter(isbayesianparam, getparams(model))
            end
        end
    end

    model |> cpu!
    savemodel(outdir, model, modelconf, inferenceconf)
end

args = parse_commandline()
main(args)

max_stars_repo_name: MartinKocour/SpeechLab

max_issues_repo_name: MartinKocour/SpeechLab

max_forks_repo_name: MartinKocour/SpeechLab


text:
test = """
    Player 1:
    9
    2
    6
    3
    1

    Player 2:
    5
    8
    4
    7
    10
    """

function parse_input(input::AbstractString)::Vector{Vector}
    decks = split.(split(input, "\n\n"), '\n'; keepempty = false)
    @assert all(startswith.(first.(decks), "Player"))
    return [parse.(Int, deck[2:end]) for deck in decks]
end

function game(combat::Function, decks)::Int
    decks = deepcopy(decks)
    winner = combat(decks)
    winning = decks[winner]
    return sum(winning .* reverse(keys(winning)))
end

function combat!(decks)::Int
    while !any(isempty.(decks))
        simple_step!(decks)
    end
    return findfirst(!isempty, decks)
end

function simple_step!(decks)
    play = popfirst!.(decks)
    winner = argmax(play)
    append!(decks[winner], play[[winner, 3 - winner]])
    return decks
end

function recursive_combat!(decks)::Int
    deal_hashes = Set()
    while !any(isempty.(decks))
        h = hash(decks)
        h ∈ deal_hashes && return 1
        push!(deal_hashes, h)
        recursive_step!(decks)
    end
    return findfirst(!isempty, decks)
end

function recursive_step!(decks)
    play = popfirst!.(decks)
    if all(@. play <= length(decks))
        subdecks = getindex.(decks, range.(1, play; step = 1))
        winner = recursive_combat!(subdecks)
    else
        winner = argmax(play)
    end
    append!(decks[winner], play[[winner, 3 - winner]])
    return decks
end

q1(decks) = game(combat!, decks)
q2(decks) = game(recursive_combat!, decks)

let decks = parse_input(test)
    @assert q1(decks) == 306
    @assert q2(decks) == 291
end

decks = parse_input(read("day22.in", String))
println("Q1: ", q1(decks))
println("Q2: ", q2(decks))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
#using TensorToolbox, Test, LinearAlgebra

T1=dimtree([8,9,5,6,7])
T2=dimtree([2,4,6,8,9])
T3=dimtree(7)

X=float(reshape(collect(1:32),(2,2,2,2,2)))
println("\n\nTesting truncation of tensor X of size ",size(X)," to htensor.")
H1=htrunc(X,T1);
err=norm(full(H1)-X)
println("norm(full(H)-X) = ",err)
@test err ≈ 0 atol=1e-10

H2=htrunc(X,T2)
err=norm(full(H2)-X)
println("norm(full(H)-X) = ",err)
@test err ≈ 0 atol=1e-10

Y=float(reshape(collect(1:2^7),(2,2,2,2,2,2,2)));
H3=htrunc(Y,T3);
@test norm(full(H3)-Y) ≈ 0 atol=1e-10

X=rand(5,4,3,2,4,2);
H=htrunc(X,dimtree(ndims(X)),maxrank=20);
@test norm(full(H)-X) ≈ 0 atol=1e-10

println("\n\n...Testing size of htensor H : ", size(H))
@test size(H) == size(X)
println("\n...Testing ndims of htensor H : ", ndims(H))
@test ndims(H) == ndims(X)

A=MatrixCell(undef,ndims(H))
for n=1:ndims(H)
  A[n]=rand(rand(1:10),size(H,n))
end
err=norm(full(ttm(H,A)) - ttm(full(H),A))
println("\n...Testing ttm for ttensor H and array of matrices A : norm(full(ttm(H,A))-ttm(full(H),A)) = ", err)
@test err ≈ 0 atol=1e-10

println("\n\n...Testing norm of htensor H.")
err=abs(norm(H) - norm(full(H)))
println("|norm(H) - norm(full(H))| = ", err)
@test err ≈ 0 atol=1e-10

println("\n\n...Testing scalar multiplication 3*H.")
err=norm(full(3*H) - 3*full(H))
println("norm(full(3*H) - 3*full(H)) = ", err )
@test err ≈ 0 atol=1e-10

println("\n\nTesting functions plus and minus for htensors H1 and H2.")
T=dimtree(3)
H1=htrunc(rand(5,4,3),T)
H2=htrunc(rand(5,4,3),T)
H=H1+H2
F=full(H1)+full(H2);
err=norm(full(H) - F)
println("norm(full(H1+H2) - (full(H1)+full(H2))) = ",err )
@test err ≈ 0 atol=1e-10
H=H1-H2
F=full(H1)-full(H2);
err=norm(full(H) - F)
println("norm(full(H1+H2) - (full(H1)+full(H2))) = ",err )
@test err ≈ 0 atol=1e-10

println("\n\n...Testing inner product.")
H=innerprod(H1,H2)
err=abs(H - innerprod(full(H1),full(H2)))
println("|innerprod(H1,H2) - innerprod(full(H1),full(H2))| = ",err )
@test err ≈ 0 atol=1e-10

println("\n\nTesting function dropdims for squeezed tensor X and squeezed htensor Y.")
X=rand(1,2,4,3,2)
H=htrunc(X)
Y=dropdims(H)
err = norm(full(Y)-dropdims(X,dims=1))
println("norm(full(Y)-dropdims(X,dims=1)) = ",err)
@test err ≈ 0 atol=1e-10
X=rand(2,1,4,3,2)
H=htrunc(X)
Y=dropdims(H)
err = norm(full(Y)-dropdims(X,dims=2))
println("norm(full(Y)-dropdims(X,dims=2)) = ",err)
@test err ≈ 0 atol=1e-10
X=rand(2,4,1,3,2)
H=htrunc(X)
Y=dropdims(H)
err = norm(full(Y)-dropdims(X,dims=3))
println("norm(full(Y)-dropdims(X,dims=3)) = ",err)
@test err ≈ 0 atol=1e-10
X=rand(2,4,3,1,2)
H=htrunc(X)
Y=dropdims(H)
err = norm(full(Y)-dropdims(X,dims=4))
println("norm(full(Y)-dropdims(X,dims=4)) = ",err)
@test err ≈ 0 atol=1e-10
X=rand(2,4,3,2,1)
H=htrunc(X)
Y=dropdims(H)
err = norm(full(Y)-dropdims(X,dims=5))
println("norm(full(Y)-dropdims(X,dims=5)) = ",err)
@test err ≈ 0 atol=1e-10

H=randhtensor([4,5,2])
println("\n\n...Testing reorthogonalization for random htensor H of size ", size(H),".")
println("Has orthogonal factor matrices: ", H.isorth)
println("After reorthogonalization: ", reorth!(H).isorth)
@test H.isorth

H=htrunc(rand(5,5,5),dimtree(3),maxrank=2);
println("Testing hrank: ",hrank(H))
@test hrank(H)==[1, 2, 2, 2, 2]


text:
"""
Convert number to opposite sign
"""
opposite(x::T) where {T<:Real} = -1*x


text:
using AdvancedHMC.Adaptation: WelfordVar, NaiveVar, WelfordCov, NaiveCov, add_sample!, get_var, get_cov, reset!
using Test, LinearAlgebra, Distributions

# Check that the estimated variance is approximately correct.
@testset "Online v.s. naive v.s. true var/cov estimation" begin
    D = 10
    n_samples = 100_000

    var_welford = WelfordVar(D)
    var_naive = NaiveVar()
    var_estimators = [var_welford, var_naive]
    cov_welford = WelfordCov(D)
    cov_naive = NaiveCov()
    cov_estimators = [cov_welford, cov_naive]
    estimators = [var_estimators..., cov_estimators...]

    for dist in [MvNormal(zeros(D), ones(D)), Dirichlet(D, 1)]
        for _ = 1:n_samples
            s = rand(dist)
            for estimator in estimators
                add_sample!(estimator, s)
            end
        end

        @test get_var(var_welford) ≈ get_var(var_naive) atol=0.1D
        for estimator in var_estimators
            @test get_var(estimator) ≈ var(dist) atol=0.1D
        end

        @test get_cov(cov_welford) ≈ get_cov(cov_naive) atol=0.1D^2
        for estimator in cov_estimators
            @test get_cov(estimator) ≈ cov(dist) atol=0.1D^2
        end

        for estimator in estimators
            reset!(estimator)
        end
    end
end


text:
using Revise
using VCFTools
using MendelImpute
using GeneticVariation
using Random
using StatsBase
using CodecZlib
using ProgressMeter
using JLD2, FileIO, JLSO
using BenchmarkTools
using GroupSlices
using TimerOutputs
using LinearAlgebra

BLAS.set_num_threads(1)

Threads.nthreads()

Random.seed!(2020)
width   = 64
tgtfile = "target.chr20.typedOnly.maf0.01.masked.vcf.gz"
reffile = "ref.chr20.w$width.maf0.01.excludeTarget.jlso"
outfile = "mendel.chr20.imputed.target.vcf.gz"
@time ph = phase(tgtfile, reffile, outfile = outfile, width = width,
    dynamic_programming = false);

# import imputed result and compare with true
X_mendel = convert_gt(Float32, outfile);
X_complete = convert_gt(Float32, "target.chr20.full.vcf.gz");
n, p = size(X_mendel);
error_rate = sum(X_mendel .!= X_complete) / n / p

max_stars_repo_name: Hua-Zhou/MendelImpute