text:
#Подсчитать число и вертикальных и горизонтальных прямолинейных перегородок (прямоугольных - нет)
include("library.jl")
function count_obstacles(r)
    back_path = BackPath(r, (Sud, West))
    print(counterPartition!(r), "\tHorizontal\n")
    BackPath(r, (Nord, West))
    print(counterPartition!(r, Sud), "\tVertical\n")
    BackPath(r, (Sud, West))
    back!(r, back_path)
end

max_stars_repo_name: listova-nastya/mirea-progs

max_issues_repo_name: listova-nastya/mirea-progs

max_forks_repo_name: listova-nastya/mirea-progs


text:
export Node, Dag, NodeType, Leaf, Inner,
       children, has_children, num_children, isleaf, isinner,
       reset_counter, foreach, foreach_reset, clear_data, filter,
       nload, nsave,
       foldup, foldup_aggregate,
       count_parents, foreach_down,
       num_nodes, num_edges, tree_num_nodes, tree_num_edges, in,
       inodes, innernodes, leafnodes, linearize,
       left_most_descendent, right_most_descendent,
       node_stats, inode_stats, leaf_stats

#####################
# types and traits
#####################

"""
A node in a directed graph. This type hierarchy is organized according
to the type of graph the node is part of.
"""
abstract type Node end

"A node in a directed acyclic graph (of which it is the root)"
abstract type Dag <: Node end

"""
A trait hierarchy denoting types of `Dag` nodes
`NodeType` defines an orthogonal type hierarchy of node types, so we can dispatch on node type regardless of the graph type.
See @ref{https://docs.julialang.org/en/v1/manual/methods/#Trait-based-dispatch-1}
"""
abstract type NodeType end

"The trait of leaf nodes (nodes without children)"
struct Leaf <: NodeType end

"The trait of inner nodes (nodes that have children)"
struct Inner <: NodeType end

#####################
# basic fields and methods
#####################

# Each `Node` is required to have fields
#  - `counter::UInt32`
#  - `data:Any`
# and a specialized method for the following functions.

"Get the node type trait of the given `Node`"
@inline NodeType(node::Node) = NodeType(typeof(node))

"Get the children of a given inner node"
function children end

#####################
# derived node functions
#####################

"Does the node have children?"
@inline has_children(n::Node)::Bool = has_children(NodeType(n), n)
@inline has_children(::Inner, n::Node)::Bool = (@assert !isempty(children(n)); true)
@inline has_children(::Leaf, n::Node)::Bool = false

"Get the number of children of a given inner node"
@inline num_children(n::Node)::Int = num_children(NodeType(n), n)
@inline num_children(::Inner, n::Node)::Int = length(children(n))
@inline num_children(::Leaf, n::Node)::Int = 0

"Is the node a leaf node?"
@inline isleaf(n::Node) = NodeType(n) isa Leaf

"Is the node an Inner node?"
@inline isinner(n::Node) = NodeType(n) isa Inner

#####################
# traversal
#####################

"Set the counter field throughout this graph (assumes no node in the graph already has the value)"
function reset_counter(node::Dag, v::Int=0)
    if node.counter != v
        node.counter = v
        if isinner(node)
            for c in children(node)
                reset_counter(c, v)
            end
        end
    end
    nothing # returning nothing helps save some allocations and time
end
import Base.foreach #extend

"Apply a function to each node in a graph, bottom up"
function foreach(f::Function, node::Dag; reset=true)
    @assert node.counter == 0 "Another algorithm is already traversing this circuit and using the `counter` field. You can use `reset_counter` to reset the counter. "
    foreach_rec(f, node)
    reset && reset_counter(node)
    nothing # returning nothing helps save some allocations and time
end

function foreach(node::Dag, f_leaf::Function, f_inner::Function)
    foreach(node) do n
        isinner(n) ? f_inner(n) : f_leaf(n)
    end
    nothing # returning nothing helps save some allocations and time
end

"Apply a function to each node in a graph, bottom up, without resetting the counter"
function foreach_rec(f::Function, node::Dag)
    if (node.counter += 1) == 1
        if isinner(node)
            for c in children(node)
                foreach_rec(f, c)
            end
        end
        f(node)
    end
    nothing # returning nothing helps save some allocations and time
end

"Apply a function to each node in a graph, bottom up, while resetting the counter"
function foreach_reset(f::Function, node::Dag)
    if node.counter != 0
        node.counter = 0
        if isinner(node)
            for c in children(node)
                foreach_reset(f, c)
            end
        end
        f(node)
    end
    nothing # returning nothing helps save some allocations and time
end

"Set all the data fields in a circuit to `nothing`"
function clear_data(node::Dag)
    foreach(x -> x.data = nothing, node)
end

# TODO: consider adding a top-down version of foreach, by either linearizing into a List,
# or by keeping a visit counter to identify processing of the last parent.

import Base.filter #extend

"""Retrieve list of nodes in graph matching predicate `p`"""
function filter(p::Function, root::Dag, ::Type{T} = Union{})::Vector where T
    results = Vector{T}()
    foreach(root) do n
        if p(n)
            if !(n isa eltype(results))
                results = collect(typejoin(eltype(results), typeof(n)), results)
            end
            push!(results, n)
        end
    end
    results
end

"Default getter to obtain data associated with a node"
@inline nload(n) = n.data

"Default setter to assign data associated with a node"
@inline nsave(n,v) = n.data = v

"""
    foldup(node::Dag,
        f_leaf::Function,
        f_inner::Function,
        ::Type{T}; nload = nload, nsave = nsave, reset=true)::T where {T}

Compute a function bottom-up on the graph.
`f_leaf` is called on leaf nodes, and `f_inner` is called on inner nodes.
Values of type `T` are passed up the circuit and given to `f_inner` as a function on the children.
"""
function foldup(node::Dag, f_leaf::Function, f_inner::Function,
               ::Type{T}; nload = nload, nsave = nsave, reset=true)::T where {T}
    @assert node.counter == 0 "Another algorithm is already traversing this circuit and using the `counter` field. You can use `reset_counter` to reset the counter. "
    v = foldup_rec(node, f_leaf, f_inner, T; nload, nsave)
    reset && reset_counter(node)
    v
end

"""
Compute a function bottom-up on the graph, without resetting the counter.
`f_leaf` is called on leaf nodes, and `f_inner` is called on inner nodes.
Values of type `T` are passed up the circuit and given to `f_inner` as a function on the children.
"""
function foldup_rec(node::Dag, f_leaf::Function, f_inner::Function, ::Type{T};
                    nload = nload, nsave = nsave)::T where {T}
    if (node.counter += 1) != 1
        return nload(node)::T
    else
        v = if isinner(node)
                callback(c) = (foldup_rec(c, f_leaf, f_inner, T; nload, nsave)::T)
                f_inner(node, callback)::T
            else
                f_leaf(node)::T
            end
        return nsave(node, v)::T
    end
end

"""
Compute a function bottom-up on the circuit.
`f_leaf` is called on leaf nodes, and `f_inner` is called on inner nodes.
Values of type `T` are passed up the circuit and given to `f_inner` in aggregate
as a vector from the children.
"""
# TODO: see whether we could standardize on `foldup` and remove this version?
function foldup_aggregate(node::Dag, f_leaf::Function, f_inner::Function, ::Type{T};
                          nload = nload, nsave = nsave, reset=true)::T where {T}
    @assert node.counter == 0 "Another algorithm is already traversing this circuit and using the `counter` field. You can use `reset_counter` to reset the counter. "
    v = foldup_aggregate_rec(node, f_leaf, f_inner, T; nload, nsave)
    reset && reset_counter(node)
    return v
end

"""
Compute a function bottom-up on the circuit, without resetting the counter.
`f_leaf` is called on leaf nodes, and `f_inner` is called on inner nodes.
Values of type `T` are passed up the circuit and given to `f_inner` in aggregate
as a vector from the children.
"""
# TODO: see whether we could standardize on `foldup` and remove this version?
function foldup_aggregate_rec(node::Dag, f_leaf::Function, f_inner::Function, ::Type{T};
                              nload = nload, nsave = nsave)::T where {T}
    if (node.counter += 1) != 1
        return nload(node)::T
    else
        v = if isinner(node)
            child_values = Vector{T}(undef, num_children(node))
            map!(c -> foldup_aggregate_rec(c, f_leaf, f_inner, T; nload, nsave)::T,
                                           child_values, children(node))
            f_inner(node, child_values)::T
        else
            f_leaf(node)::T
        end
        return nsave(node, v)::T
    end
end

"Set the `counter` field of each node to its number of parents"
function count_parents(node::Dag)
    foreach(noop,node;reset=false)
end

"Apply a function to each node in a graph, top down"
function foreach_down(f::Function, node::Dag; setcounter=true)
    setcounter && count_parents(node)
    foreach_down_rec(f, node)
    nothing
end

"Apply a function to each node in a graph, top down, without setting the counter first"
function foreach_down_rec(f::Function, n::Node)
    if ((n.counter -= 1) == 0)
        f(n)
        if isinner(n)
            for c in children(n)
                foreach_down_rec(f, c)
            end
        end
    end
    nothing
end

#####################
# methods using circuit traversal
#####################

"""
    num_nodes(node::Dag)

Count the number of nodes in the `Dag`
"""
function num_nodes(node::Dag)
    count::Int = 0
    foreach(node) do n
        count += 1
    end
    count
end

"""
    tree_num_nodes(node::Dag)::BigInt

Compute the number of nodes in of a tree-unfolding of the `Dag`.
"""
function tree_num_nodes(node::Dag)::BigInt
    @inline f_leaf(n) = one(BigInt)
    @inline f_inner(n, call) = (1 + mapreduce(c -> call(c), +, children(n)))
    foldup(node, f_leaf, f_inner, BigInt)
end

"""
    tree_num_edges(node::Dag)::BigInt

Compute the number of edges in the tree-unfolding of the `Dag`.
"""
function tree_num_edges(node::Dag)::BigInt
    @inline f_leaf(n) = zero(BigInt)
    @inline f_inner(n, call) = (num_children(n) + mapreduce(c -> call(c), +, children(n)))
    foldup(node, f_leaf, f_inner, BigInt)
end

"Number of edges in the `Dag`"
function num_edges(node::Dag)
    count::Int = 0
    foreach(node) do n
        count += num_children(n)
    end
    count
end

"Is the node contained in the `Dag`?"
function Base.in(needle::Dag, circuit::Dag)
    contained::Bool = false
    foreach(circuit) do n
        contained |= (n == needle)
    end
    contained
end

"Get the list of inner nodes in a given graph"
inodes(c::Dag) = filter(isinner, c)

"Get the list of inner nodes in a given graph"
innernodes(c::Dag) = inodes(c)

"Get the list of leaf nodes in a given graph"
leafnodes(c::Dag) = filter(isleaf, c)

"Order the `Dag`'s nodes bottom-up in a list (with optional element type)"
@inline linearize(r::Dag, ::Type{T} = Union{}) where T = filter(x -> true, r, typejoin(T,typeof(r)))

"""
Return the left-most descendent.
"""
function left_most_descendent(root::Dag)::Dag
    while isinner(root)
        root = children(root)[1]
    end
    root
end

"""
Return the right-most descendent.
"""
function right_most_descendent(root::Dag)::Dag
    while isinner(root)
        root = children(root)[end]
    end
    root
end

#####################
# debugging methods (not performance critical)
#####################

# When you suspect there is a bug but execution halts, it may be because of
# pretty printing a huge recursive graph structure.
# To safeguard against that case, we set a default show:
Base.show(io::IO, c::Node) = print(io, "$(typeof(c))($(hash(c)))")

"Give count of types and fan-ins of all nodes in the graph"
node_stats(c::Dag) = merge(leaf_stats(c), inode_stats(c))

"Give count of types and fan-ins of inner nodes in the graph"
function inode_stats(c::Dag)
    groups = groupby(e -> (typeof(e),num_children(e)), inodes(c))
    map_values(v -> length(v), groups, Int)
end

"Give count of types of leaf nodes in the graph"
function leaf_stats(c::Node)
    groups = groupby(e -> typeof(e), leafnodes(c))
    map_values(v -> length(v), groups, Int)
end

ext: jl


text:
using SpecialFunctions, StatsFuns

@testset "logmvgamma" begin
    @testset "type behavior" for eltya in (Float32, Float64)
        p = rand(1:50)
        a = rand(eltya)
        # add p since loggamma is only define for positive arguments
        @test typeof(logmvgamma(p, a + p)) == eltya
    end

    @testset "consistent with loggamma" for eltya in (Float32, Float64)
        #  Γ₁(a) = Γ(a), Γ₂(a) = √π Γ(a) Γ(a - 0.5), etc
        a = rand(eltya) + 1 # add one since loggamma is only define for positive arguments
        @test logmvgamma(1, a) ≈ loggamma(a)
        @test logmvgamma(2, a) ≈    eltya(0.5logπ) + loggamma(a) + loggamma(a - eltya(0.5))
        @test logmvgamma(3, a) ≈ eltya((3/2)*logπ) + loggamma(a) + loggamma(a - eltya(0.5)) + loggamma(a - one(a))
    end

    @testset "consistent with itself" for eltya in (Float32, Float64)
        #  Γᵢ(a) = (π^{i-1/2}) Γ(a) Γᵢ₋₁(a - 0.5)
        for p in 1:50
            a = rand(eltya) + p # add p since loggamma is only define for positive arguments
            @test logmvgamma(p, a) ≈ eltya((p/2-1/2)*logπ) + loggamma(a) + logmvgamma(p - 1, a - eltya(0.5))
        end
    end
end

@testset "logmvbeta" begin
    @testset "symmetry" for eltya in (Float32, Float64)
                        for eltyb in (Float32, Float64)
            #  Bᵢ(a, b) = Bᵢ(b, a)
            for p in 1:50
                a = rand(eltya) + p
                b = rand(eltyb) + p
                @test logmvbeta(p, a, b) ≈ logmvbeta(p, b, a)
            end
        end
    end

    @testset "consistent with logbeta" for eltya in (Float32, Float64)
                                       for eltyb in (Float32, Float64)
            #  B₁(a, b) = B(a, b)
            a = rand(eltya)
            b = rand(eltyb)
            @test logmvbeta(1, a, b) ≈ logbeta(a, b)
        end
    end

    @testset "type promotion behaves" for eltya in (Float32, Float64)
                                      for eltyb in (Float32, Float64)
            a = rand(eltya)
            b = rand(eltyb)
            T = Base.promote_eltype(eltya, eltyb)
            @test typeof(logmvbeta(1, a, b)) == T
        end
    end
end

# https://github.com/JuliaStats/StatsFuns.jl/issues/115
@testset "support of binomial distribution" begin
    @test iszero(binompdf(1, 0.5, prevfloat(1.0)))
    @test iszero(binompdf(1, 0.5, nextfloat(1.0)))
    @test binomlogpdf(1, 0.5, prevfloat(1.0)) == -Inf
    @test binomlogpdf(1, 0.5, nextfloat(1.0)) == -Inf
end

@testset "binom special cases" begin
    for (n, p, k) in ((5, 0.0, 0), (5, 1.0, 5))
        @test iszero(binomlogpdf(n, p, k))
        @test isone(binompdf(n, p, k))
    end
end


text:
using Plots, HDF5

@testset "HDF5_Plots" begin
	fname = "tmpplotsave.hdf5"
	hdf5()

	x = 1:10
	psrc=plot(x, x.*x); #Create some plot
	Plots.hdf5plot_write(psrc, fname)

	#Read back file:
	gr() #Choose some fast backend likely to work in test environment.
	pread = Plots.hdf5plot_read(fname)

	#Make sure data made it through:
	@test psrc.subplots[1].series_list[1][:x] == pread.subplots[1].series_list[1][:x]
	@test psrc.subplots[1].series_list[1][:y] == pread.subplots[1].series_list[1][:y]

	#display(pread) #Don't display. Regression env might not support
end #testset


text:
#######################################################################
# Analysis of lower bounds
#######################################################################

include("data.jl")
include("process-code.jl")
include("process-text.jl")
include("process-pkgs.jl")

max_stars_repo_name: prl-julia/julia-sub

max_issues_repo_name: prl-julia/julia-sub

max_forks_repo_name: prl-julia/julia-sub


text:
module DateTests

using Test
using Dates

include("accessors.jl")
include("adjusters.jl")
include("query.jl")
include("periods.jl")
include("ranges.jl")
include("rounding.jl")
include("types.jl")
include("io.jl")
include("arithmetic.jl")
include("conversions.jl")

end

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
"""
    using PackageCompatUI
    compat_ui()

Interactive terminal interface to the `[compat]` section of
`Project.toml` for the *currently active project*.
"""
module PackageCompatUI

export compat_ui

using Pkg
using Downloads: download
using JSON3
using Scratch
using Git
using Dates

include("menu.jl")

"""
    compat_ui()

Interactive terminal interface to the `[compat]` section of
`Project.toml` for the currently active project.

*Keyword arguments*:

* `pagesize`: Maximum number of lines used for the UI. Default is 20.
* `dates`: Whether to show registration date for package versions.
  Default is true. If disabled the UI starts faster.

*Controls*:

Move up and down with arrow keys. Arrow right and left to enter and
leave packages. Use Enter to toggle compatible versions. Quit and save
with `q`. Online help with `?`.

When toggling a compatible version, all semantically versioning
compatible versions are also set or unset.

*Colors*:

* Red: No compat information available for the package.
* Yellow: Compat declared for some versions of the package, but not the latest.
* Green: Compatible with the latest version of the package.
* Gray: Package is not registered.
"""
function compat_ui(;pagesize = 20, dates = true)
    if VERSION < v"1.6.0-DEV.1554"
        println("Julia version must be at least 1.6.0-DEV.1554.")
        return
    end
    project = Pkg.Types.read_project(Base.active_project())
    needs_compat = filter(p -> !Pkg.Types.is_stdlib(last(p)),
                          merge(project.deps, project.extras))
    selected = 1
    versions = Dict{String, Vector{String}}()
    version_dates = Dict{String, Dict{String, String}}()
    for (package, uuid) in needs_compat
        print("loading data for ", rpad(package, 60), "\r")
        versions[package], version_dates[package] = find_versions(uuid, dates)
    end
    print(rpad("\r", 80), "\r")
    versions["julia"] = get_julia_versions()
    version_dates["julia"] = Dict{String, String}()
    packages = vcat("julia", sort(collect(keys(needs_compat))))
    data = MenuData(packages, extract_compat_string(project.compat),
                    versions, version_dates)
    menu = Menu(data, pagesize)
    menu.cursor[] = selected
    request(menu, cursor = menu.cursor)
    set_project_compat!(project, data.compat)
    Pkg.Types.write_project(project, Base.active_project())
    return
end

# Compatibility layer to handle Pkg internals changes between Julia 1.6
# and Julia 1.7.
extract_compat_string(compat::Dict) = Dict(k => extract_compat_string(v)
                                           for (k, v) in compat)
extract_compat_string(compat::String) = compat
extract_compat_string(compat) = compat.str
set_project_compat!(project_compat::Dict{String, String}, new_compat) =
    merge!(empty!(project_compat), new_compat)
function set_project_compat!(project, compat)
    if project.compat isa Dict{String, String}
        project.compat = compat
    else
        for (package, compat_string) in compat
            Pkg.Operations.set_compat(project, package, compat_string)
        end
    end
end

function gitcmd(repo)
    git = Git.git(["-C", repo])
    return (x...) -> readlines(`$git $x`)
end

# Find registry version information for the package with the given
# `uuid`. This may pick up information from multiple registries.
# Look up registration date of versions if `find_dates` is true.
function find_versions(uuid, find_dates)
    # TODO: Change this to be a vector of `VersionNumber` unless that
    # causes problems later.
    versions = String[]
    dates = Dict{String, String}()
    if VERSION >= v"1.7-"
        for registry in Pkg.Registry.reachable_registries()
            if haskey(registry.pkgs, uuid)
                package = registry.pkgs[uuid]
                package_dir = package.path
                info = Pkg.Registry.registry_info(package)
                append!(versions, string.(keys(info.version_info)))
                if find_dates
                    find_version_dates!(dates, registry.name,
                                        registry.repo, package_dir)
                end
            end
        end
    else
        for registry in Pkg.Types.collect_registries()
            reg_data = Pkg.Types.read_registry(joinpath(registry.path,
                                                        "Registry.toml"))
            if haskey(reg_data["packages"], string(uuid))
                package_dir = reg_data["packages"][string(uuid)]["path"]
                versions_file = joinpath(registry.path, package_dir,
                                         "Versions.toml")
                append!(versions, keys(Pkg.TOML.parsefile(versions_file)))
                if find_dates
                    find_version_dates!(dates, reg_data["name"],
                                        reg_data["repo"], package_dir)
                end
            end
        end
    end

    return sort(versions, by = VersionNumber), dates
end

fetched_this_session = Set{String}()

function find_version_dates!(dates, registry_name, registry_repo, package_dir)
    repo = @get_scratch!(registry_name)
    git = gitcmd(repo)
    if !isdir(joinpath(repo, "refs"))
        println("\nCloning the $(registry_name) Registry. This might take some time.")
        git("clone", "--bare", registry_repo, ".")
        println()
    else
        if registry_name ∉ fetched_this_session
            git("fetch", "-q", "origin", "master:master")
            push!(fetched_this_session, registry_name)
        end
    end
    time = 0

    for blame in git("blame", "-p", joinpath(package_dir, "Versions.toml"))
        if startswith(blame, "committer-time")
            time = parse(Int, last(split(blame)))
        elseif startswith(blame, "\t")
            m = match(r"\[\"(.*)\"\]", blame)
            if !isnothing(m)
                version = only(m.captures)
                dates[version] = Dates.format(unix2datetime(time), "yyyy-mm-dd")
            end
        end
    end
end

function get_julia_versions()
    io = IOBuffer()
    download("https://julialang-s3.julialang.org/bin/versions.json", io)
    seekstart(io)
    json = JSON3.read(io)
    versions = sort(VersionNumber.(String.(keys(json))))
    filter!(v -> v >= v"1" && isempty(v.build) && isempty(v.prerelease),
            versions)
    return string.(versions)
end

function find_compatible_versions(versions, compat)
    if isempty(compat)
        return String[]
    end
    semver = Pkg.Types.semver_spec(compat)
    return filter(v -> VersionNumber(v) in semver, versions)
end

function compute_compat(version_strings, compatible_version_strings)
    versions = VersionNumber.(base_version.(version_strings))
    compatible_versions = VersionNumber.(base_version.(compatible_version_strings))
    compat = String[]

    # 0.0.x versions. Just list all of them individually.
    append!(compat, string.(filter(is00x, compatible_versions)))

    # 0.x.y versions.
    xy_versions = filter(is0xy, versions)
    for x in unique([v.minor for v in xy_versions])
        x_versions = filter(v -> v.minor == x, xy_versions)
        for (i, v) in enumerate(x_versions)
            if v in compatible_versions
                if issubset(x_versions[(i + 1):end], compatible_versions)
                    push!(compat, compact_version(v))
                    break
                end
                push!(compat, string("=", v))
            end
        end
    end

    # x.y.z versions.
    xyz_versions = filter(isxyz, versions)
    for x in unique([v.major for v in xyz_versions])
        x_versions = filter(v -> v.major == x, xyz_versions)
        for (i, v) in enumerate(x_versions)
            if v in semver_spec(compat)
                continue
            end
            if v in compatible_versions
                if issubset(x_versions[(i + 1):end], compatible_versions)
                    push!(compat, compact_version(v))
                    continue
                end
                y = v.minor
                y_versions = filter(v -> v.minor == y, x_versions[(i + 1):end])
                if issubset(y_versions, compatible_versions)
                    push!(compat, compact_version(v, "~"))
                    continue
                end
                push!(compat, string("=", v))
            end
        end
    end
    return join(compat, ", ")
end

# Remove prerelease and build numbers.
function base_version(version::String)
    v = VersionNumber(version)
    return string(VersionNumber(v.major, v.minor, v.patch, (), ()))
end

# Extension of `Pkg.Types.semver_spec` for vector of strings.
function semver_spec(s::Vector{String})
    if isempty(s)
        return Pkg.Types.VersionSpec([])
    end
    return Pkg.Types.semver_spec(join(s, ", "))
end

is00x(v::VersionNumber) = v.major == v.minor == 0
is0xy(v::VersionNumber) = v.major == 0 && v.minor > 0
isxyz(v::VersionNumber) = v.major > 0

function compact_version(v::VersionNumber, s = "")
    if v.patch > 0
        return string(s, v)
    end
    if v.minor > 0 || s == "~"
        return string(s, v.major, ".", v.minor)
    end
    return string(s, v.major)
end

end # module


text:
using Random
using Sockets
include("header.jl")
include("black_jack.jl")
include("slot_machine.jl")
include("roulette.jl")
include("withdraw.jl")
include("strings.jl")
include("cryptomat/cryptomat.jl")
include("tables.jl")
include("restaurant.jl")
include("bathroom.jl")

function gamble(p::Player)
    while true
        print_dict(p, "spacer")
        print_dict(p, "gamble_0")
        printGames(p)
        print_dict(p, "gamble_1")
        s = readline(p.socket)
        if s == ""
            print_dict(p, "repeat")
            continue
        elseif s == "black_jack"
            p.current_game = black_jack
            break
        elseif s == "slot_machine"
            p.current_game = slot_machine
            break
        elseif s == "roulette"
            p.current_game = roulette
            break
        else
            print_dict(p, "irritated")
            print_dict(p, "repeat")
            continue
        end
    end
    print_dict(p, "spacer")
    while true
        if p.current_game == slot_machine
            break
        end
        print_dict(p, "table_0")
        s = readline(p.socket)
        if s == ""
            print_dict(p, "repeat")
            continue
        elseif s == "j"
            if join_table(p)
                break
            else
                print_dict(p, "gamble_3")
                p.status = reception
                return
            end

        elseif s == "c"
            create_table(p)
            continue
        else
            print_dict(p, "irritated")
            print_dict(p, "repeat")
            continue
        end
    end

    if p.current_game == black_jack
        print_dict(p, "black_jack_welcome")
    elseif p.current_game == slot_machine
        print_dict(p, "slot_machine_welcome")
    elseif p.current_game == roulette
        print_dict(p, "roulette_welcome")
    end

    while true
        print_dict(p, "spacer")
        if p.current_game == black_jack
            play_black_jack(p)
        elseif p.current_game == slot_machine
            play_slot_machine(p)
        elseif p.current_game == roulette
            play_roulette(p)
        end

        if p.status == reception
            print_dict(p, "gamble_4")
            return
        end
        while true
            print_dict(p, "gamble_2")
            s = readline(p.socket)
            if s == ""
                print_dict(p, "repeat")
                continue
            elseif s == "n"
                print_dict(p, "gamble_3")
                p.status = reception
                return
            elseif s == "y"
                break
            else
                print_dict(p, "irritated")
                print_dict(p, "repeat")
                continue
            end
        end
    end
end

function receptionDesk(p::Player)
    print_dict(p, "spacer")
    printBalance(p)
    print_dict(p, "reception_0")
    s = readline(p.socket)
    if s == ""
        print_dict(p, "reception_1")
    elseif s == "g"
        p.status = gambling
    elseif s == "w"
        withdraw(p)
    elseif s == "b"
        print_dict(p, "reception_2")
        bathroom(p)
    elseif s == "r"
        restaurant(p)
    elseif s == "l"
        print_dict(p, "spacer")
        if p.balance < 0
            print_dict(p, "debt_1")
        else
            print_dict(p, "exit")
            close(p.socket)
        end
    end
end

function main(socket)
    p = Player(0, reception, slot_machine, socket, rand(Int), "", "", false)
    print_dict(p, "spacer")
    print_dict(p, "welcome")
    while true
        if p.status == reception
            receptionDesk(p)
        elseif p.status == gambling
            gamble(p)
        end
    end
end
####################################
####################################

function preload_aesSeed()
    f = open_file_try("assets/aes.seed", "r")
    aesSeed = read(f)
    close(f)
    return aesSeed
end

DEBUG = true

aesSeed = preload_aesSeed()
server = listen(IPv6(0),6969)
println("Waiting for people to enter the casino..")
while true
    socket = accept(server)
    println("Accepted: $(getsockname(socket))")
    @async begin
        try
            main(socket)
        catch err
            if DEBUG
                if isa(err, Base.IOError) || isa(err, ArgumentError)
                    println("connection ended with $err")
                else
                    bt = catch_backtrace()
                    println()
                    showerror(stderr, err, bt)
                end
            else
                println("connection ended with $err")
            end

        end
    end
end


text:
export estimate_parameters, uniform_parameters, estimate_parameters_em, estimate_parameters_cached!

using LogicCircuits: num_nodes
using StatsFuns: logsumexp, logaddexp, logsubexp
using CUDA
using LoopVectorization

"""
Maximum likilihood estimation of parameters given data

use_gpu: If set to `true`, use gpu learning no matter which device `data` is in.

bagging support: If `pc` is a SharedProbCircuit and data is an array of DataFrames
  with the same number of "components", learn each circuit with its corresponding
  dataset.
"""
function estimate_parameters(pc::ProbCircuit, data; pseudocount::Float64,
                             use_sample_weights::Bool = true, use_gpu::Bool = isgpu(data), entropy_reg::Float64 = 0.0)
    estimate_single_circuit_parameters(pc, data; pseudocount, use_sample_weights, use_gpu, entropy_reg)
end
function estimate_parameters(spc::SharedProbCircuit, data; pseudocount::Float64,
                             use_sample_weights::Bool = true, use_gpu::Bool = isgpu(data), entropy_reg::Float64 = 0.0)
    @assert num_components(spc) == length(data) "SharedProbCircuit and data have different number of components: $(num_components(spc)) and $(length(data)), resp."

    map(1 : num_components(spc)) do component_idx
        estimate_single_circuit_parameters(spc, data[component_idx]; pseudocount, use_sample_weights, use_gpu,
                                           component_idx, entropy_reg)
    end
end
function estimate_single_circuit_parameters(pc::ProbCircuit, data; pseudocount::Float64,
                                            use_sample_weights::Bool = true, use_gpu::Bool = isgpu(data),
                                            component_idx::Integer = 0, entropy_reg::Float64 = 0.0)
    if isweighted(data)
        # `data' is weighted according to its `weight' column
        data, weights = split_sample_weights(data)
    else
        use_sample_weights = false
    end

    @assert isbinarydata(data) || isfpdata(data) "Probabilistic circuit parameter estimation for binary/floating point data only"
    bc = BitCircuit(pc, data)
    if isgpu(data)
        use_gpu = true
    end
    params = if use_gpu
        if use_sample_weights
            estimate_parameters_gpu(to_gpu(bc), data, pseudocount; weights, entropy_reg)
        else
            estimate_parameters_gpu(to_gpu(bc), data, pseudocount; entropy_reg)
        end
    else
        if use_sample_weights
            estimate_parameters_cpu(bc, data, pseudocount; weights, entropy_reg)
        else
            estimate_parameters_cpu(bc, data, pseudocount; entropy_reg)
        end
    end

    if pc isa SharedProbCircuit
        estimate_parameters_cached!(pc, bc, params, component_idx)
    else
        estimate_parameters_cached!(pc, bc, params)
    end
    params
end

function estimate_parameters_cached!(pc::SharedProbCircuit, bc, params, component_idx; exp_update_factor = 0.0)
    log_exp_factor = log(exp_update_factor)
    log_1_exp_factor = log(1.0 - exp_update_factor)
    foreach(pc) do pn
        if is⋁gate(pn)
            if num_children(pn) == 1
                pn.log_probs[:, component_idx] .= zero(Float64)
            else
                id = (bc.node2id[pn]::⋁NodeIds).node_id
                @inbounds els_start = bc.nodes[1,id]
                @inbounds els_end = bc.nodes[2,id]
                @inbounds @views pn.log_probs[:, component_idx] .= logaddexp.(log_exp_factor .+ pn.log_probs[:, component_idx], log_1_exp_factor .+ params[els_start:els_end])
                @assert isapprox(sum(exp.(pn.log_probs[:, component_idx])), 1.0, atol=1e-3) "Parameters do not sum to one locally: $(sum(exp.(pn.log_probs))); $(pn.log_probs)"
                pn.log_probs[:, component_idx] .-= logsumexp(pn.log_probs[:, component_idx]) # normalize away any leftover error
            end
        end
    end
    nothing
end
function estimate_parameters_cached!(pc::ProbCircuit, pbc; exp_update_factor = 0.0)
    if isgpu(pbc)
        pbc = to_cpu(pbc)
    end
    bc = pbc.bitcircuit
    params = pbc.params
    estimate_parameters_cached!(pc, bc, params; exp_update_factor)
end
function estimate_parameters_cached!(pc::ProbCircuit, bc, params; exp_update_factor = 0.0)
    log_exp_factor = log(exp_update_factor)
    log_1_exp_factor = log(1.0 - exp_update_factor)

    foreach(pc) do pn
        if is⋁gate(pn)
            if num_children(pn) == 1
                pn.log_probs .= zero(Float64)
            else
                id = (bc.node2id[pn]::⋁NodeIds).node_id
                @inbounds els_start = bc.nodes[1,id]
                @inbounds els_end = bc.nodes[2,id]
                @inbounds @views pn.log_probs .= logaddexp.(log_exp_factor .+ pn.log_probs, log_1_exp_factor .+ params[els_start:els_end])
                @assert isapprox(sum(exp.(pn.log_probs)), 1.0, atol=1e-3) "Parameters do not sum to one locally: $(sum(exp.(pn.log_probs))); $(pn.log_probs)"
                pn.log_probs .-= logsumexp(pn.log_probs) # normalize away any leftover error
            end
        end
    end

    nothing
end

function estimate_parameters_cpu(bc::BitCircuit, data, pseudocount; weights = nothing, entropy_reg::Float64 = 0.0)
    # no need to synchronize, since each computation is unique to a decision node
    if weights === nothing && isbinarydata(data)
        node_counts = Vector{UInt}(undef, num_nodes(bc))
    else
        node_counts = Vector{Float64}(undef, num_nodes(bc))
    end

    # Batch the data for a cleaner parameter learning pipeline
    if !isbatched(data)
        data = [data]
        if weights !== nothing
            weights = [weights]
        end
    end

    # rescale pseudocount using the average weight of samples
    if weights !== nothing
        pseudocount = pseudocount * mapreduce(sum, +, weights) / num_examples(data)
    end

    edge_counts::Vector{Float64} = zeros(Float64, num_elements(bc))
    parent_node_counts::Vector{Float64} = zeros(Float64, num_elements(bc))

    @inline function weighted_count_ones(bits::UInt64, start_idx::Number, end_idx::Number, weights)
        count::Float64 = 0.0
        op_bits::UInt64 = bits
        @inbounds for idx = 0 : (end_idx - start_idx)
            count += weights[start_idx + idx] * (op_bits & UInt64(0x1))
            op_bits = (op_bits >> 1)
        end
        count
    end

    @inline function on_node(flows, values, dec_id, weights::Nothing)
        node_counts[dec_id] = sum(1:size(flows,1)) do i
            count_ones(flows[i, dec_id])
        end
    end
    @inline function on_node(flows, values, dec_id, weights)
        node_counts[dec_id] = sum(1:size(flows,1)) do i
            weighted_count_ones(flows[i, dec_id], i * 64 - 63, min(i * 64, length(weights)), weights)
        end
    end
    @inline function on_node(flows, values::Matrix{<:AbstractFloat}, dec_id, weights::Nothing)
        node_counts[dec_id] = sum(1:size(flows,1)) do i
            flows[i, dec_id]
        end
    end
    @inline function on_node(flows, values::Matrix{<:AbstractFloat}, dec_id, weights)
        node_counts[dec_id] = sum(1:size(flows,1)) do i
            flows[i, dec_id] * weights[i]
        end
    end

    @inline function estimate(element, decision, edge_count)
        edge_counts[element] += edge_count
        parent_node_counts[element] += node_counts[decision]
    end

    @inline function on_edge(flows, values, prime, sub, element, grandpa, single_child, weights::Nothing)
        if !single_child
            edge_count = sum(1:size(flows,1)) do i
                count_ones(values[i, prime] & values[i, sub] & flows[i, grandpa])
            end
            estimate(element, grandpa, edge_count)
        end # no need to estimate single child params, they are always prob 1
    end
    @inline function on_edge(flows, values, prime, sub, element, grandpa, single_child, weights)
        if !single_child
            edge_count = sum(1:size(flows,1)) do i
                weighted_count_ones(values[i, prime] & values[i, sub] & flows[i, grandpa],
                                    i * 64 - 63, min(i * 64, length(weights)), weights)
            end
            estimate(element, grandpa, edge_count)
        end # no need to estimate single child params, they are always prob 1
    end
    @inline function on_edge(flows, values::Matrix{<:AbstractFloat}, prime, sub, element, grandpa, single_child, weights::Nothing)
        if !single_child
            edge_count = sum(1:size(flows,1)) do i
                if values[i, grandpa] == 0.0
                    0.0
                else
                    values[i, prime] * values[i, sub] / values[i, grandpa] * flows[i, grandpa]
                end
            end
            estimate(element, grandpa, edge_count)
        end # no need to estimate single child params, they are always prob 1
    end
    @inline function on_edge(flows, values::Matrix{<:AbstractFloat}, prime, sub, element, grandpa, single_child, weights)
        if !single_child
            edge_count = sum(1:size(flows,1)) do i
                if values[i, grandpa] == 0.0
                    0.0
                else
                    values[i, prime] * values[i, sub] / values[i, grandpa] * flows[i, grandpa] * weights[i]
                end
            end
            estimate(element, grandpa, edge_count)
        end # no need to estimate single child params, they are always prob 1
    end

    v, f = nothing, nothing
    if weights != nothing
        map(zip(data, weights)) do (d, w)
            v, f = satisfies_flows(bc, d, v, f; on_node = on_node, on_edge = on_edge, weights = w)
        end
    else
        map(data) do d
            v, f = satisfies_flows(bc, d, v, f; on_node = on_node, on_edge = on_edge, weights = nothing)
        end
    end

    # Entropy regularization
    if entropy_reg > 1e-8
        total_data_counts = (weights === nothing) ? Float64(num_examples(data)) : mapreduce(sum, +, weights)
        apply_entropy_reg_cpu(bc; log_params = edge_counts, edge_counts, total_data_counts, pseudocount, entropy_reg)
    else
        # Reuse `edge_counts` to store log_params to save space and time.
        for i = 1 : num_elements(bc)
            @inbounds edge_counts[i] = log((edge_counts[i] + pseudocount / num_elements(bc.nodes, bc.elements[1, i])) / (parent_node_counts[i] + pseudocount))
        end

        edge_counts # a.k.a. log_params
    end
end

function estimate_parameters_gpu(bc::BitCircuit, data, pseudocount; weights = nothing, entropy_reg::Float64 = 0.0)
    # rescale pseudocount using the average weight of samples
    if weights !== nothing
        if isbatched(data)
            pseudocount = pseudocount * mapreduce(sum, +, weights) / num_examples(data)
        else
            pseudocount = pseudocount * sum(weights) / size(weights, 1)
        end
    end

    node_counts::CuVector{Float32} = CUDA.zeros(Float32, num_nodes(bc))
    edge_counts::CuVector{Float32} = CUDA.zeros(Float32, num_elements(bc))
    # need to manually cudaconvert closure variables
    node_counts_device = CUDA.cudaconvert(node_counts)
    edge_counts_device = CUDA.cudaconvert(edge_counts)

    @inline function on_node(flows, values, dec_id, chunk_id, flow, weight::Nothing)
        c::Float32 = CUDA.count_ones(flow) # cast for @atomic to be happy
        CUDA.@atomic node_counts_device[dec_id] += c
    end
    @inline function on_node(flows, values, dec_id, bit_idx, flow, weight::Float32)
        c::Float32 = ((flow >> bit_idx) & UInt64(0x1)) * weight # cast for @atomic to be happy
        CUDA.@atomic node_counts_device[dec_id] += c
    end
    @inline function on_node(flows, values, dec_id, bit_idx, flow::AbstractFloat, weight::Nothing)
        c::Float32 = flow # cast for @atomic to be happy
        CUDA.@atomic node_counts_device[dec_id] += c
    end
    @inline function on_node(flows, values, dec_id, bit_idx, flow::AbstractFloat, weight::Float32)
        c::Float32 = flow * weight # cast for @atomic to be happy
        CUDA.@atomic node_counts_device[dec_id] += c
    end

    @inline function on_edge(flows, values, prime, sub, element, grandpa, chunk_id, edge_flow, single_child, weight::Nothing)
        if !single_child
            c::Float32 = CUDA.count_ones(edge_flow) # cast for @atomic to be happy
            CUDA.@atomic edge_counts_device[element] += c
        end
    end
    @inline function on_edge(flows, values, prime, sub, element, grandpa, bit_idx, edge_flow, single_child, weight::Float32)
        if !single_child
            c::Float32 = ((edge_flow >> bit_idx) & UInt64(0x1)) * weight # cast for @atomic to be happy
            CUDA.@atomic edge_counts_device[element] += c
        end
    end
    @inline function on_edge(flows, values, prime, sub, element, grandpa, sample_idx,
                             edge_flow::AbstractFloat, single_child, weight::Nothing)
        if !single_child
            c::Float32 = edge_flow # cast for @atomic to be happy
            CUDA.@atomic edge_counts_device[element] += c
        end
    end
    @inline function on_edge(flows, values, prime, sub, element, grandpa, sample_idx,
                             edge_flow::AbstractFloat, single_child, weight::Float32)
        if !single_child
            c::Float32 = edge_flow * weight # cast for @atomic to be happy
            CUDA.@atomic edge_counts_device[element] += c
        end
    end

    if isbatched(data)
        v, f = nothing, nothing
        if weights != nothing
            map(zip(data, weights)) do (d, w)
                if w != nothing
                    w = to_gpu(w)
                end
                v, f = satisfies_flows(to_gpu(bc), to_gpu(d), v, f; on_node = on_node, on_edge = on_edge, weights = w)
            end
        else
            map(data) do d
                v, f = satisfies_flows(to_gpu(bc), to_gpu(d), v, f; on_node = on_node, on_edge = on_edge, weights = nothing)
            end
        end
    else
        if weights != nothing
            weights = to_gpu(weights)
        end
        v, f = satisfies_flows(to_gpu(bc), to_gpu(data); on_node = on_node, on_edge = on_edge, weights = weights)
    end

    CUDA.unsafe_free!(v) # save the GC some effort
    CUDA.unsafe_free!(f) # save the GC some effort

    # Entropy regularization
    if entropy_reg > 1e-8
        total_data_counts = (weights === nothing) ? Float64(num_examples(data)) : mapreduce(sum, +, weights)
        apply_entropy_reg_gpu(bc; log_params = edge_counts, edge_counts, total_data_counts, pseudocount, entropy_reg)
    else
        # TODO: reinstate simpler implementation once https://github.com/JuliaGPU/GPUArrays.jl/issues/313 is fixed and released
        @inbounds parents = bc.elements[1,:]
        @inbounds parent_counts = node_counts[parents]
        @inbounds parent_elcount = bc.nodes[2,parents] .- bc.nodes[1,parents] .+ 1
        params = log.((edge_counts .+ (pseudocount ./ parent_elcount))
                        ./ (parent_counts .+ pseudocount))

        to_cpu(params)
    end
end

"""
Uniform distribution
"""
function uniform_parameters(pc::ProbCircuit; perturbation::Float64 = 0.0)
    foreach(pc) do pn
        if is⋁gate(pn)
            if num_children(pn) == 1
                pn.log_probs .= 0.0
            else
                if perturbation < 1e-8
                    pn.log_probs .= log.(ones(Float64, num_children(pn)) ./ num_children(pn))
                else
                    unnormalized_probs = map(x -> 1.0 - perturbation + x * 2 * perturbation, rand(num_children(pn)))
                    pn.log_probs .= log.(unnormalized_probs ./ sum(unnormalized_probs))
                end
            end
        end
    end
end

"""
Expectation maximization parameter learning given missing data
"""
function estimate_parameters_em(pc::ProbCircuit, data; pseudocount::Float64, entropy_reg::Float64 = 0.0,
                                use_sample_weights::Bool = true, use_gpu::Bool = isgpu(data),
                                exp_update_factor::Float64 = 0.0, update_per_batch::Bool = false)
    if update_per_batch && isbatched(data)
        estimate_parameters_em_per_batch(pc, data; pseudocount, entropy_reg,
                                         use_sample_weights, use_gpu, exp_update_factor)
    else
        if isweighted(data)
            # `data' is weighted according to its `weight' column
            data, weights = split_sample_weights(data)
        else
            use_sample_weights = false
        end

        pbc = ParamBitCircuit(pc, data)
        if isgpu(data)
            use_gpu = true
        elseif use_gpu && !isgpu(data)
            data = to_gpu(data)
        end

        params = if use_gpu
            if !isgpu(data)
                data = to_gpu(data)
            end
            if use_sample_weights
                estimate_parameters_gpu(to_gpu(pbc), data, pseudocount; weights, entropy_reg)
            else
                estimate_parameters_gpu(to_gpu(pbc), data, pseudocount; entropy_reg)
            end
        else
            if use_sample_weights
                estimate_parameters_cpu(pbc, data, pseudocount; weights, entropy_reg)
            else
                estimate_parameters_cpu(pbc, data, pseudocount; entropy_reg)
            end
        end

        estimate_parameters_cached!(pc, pbc.bitcircuit, params; exp_update_factor)

        params
    end
end
function estimate_parameters_em_per_batch(pc::ProbCircuit, data; pseudocount::Float64, entropy_reg::Float64 = 0.0,
                                          use_sample_weights::Bool = true, use_gpu::Bool = isgpu(data), exp_update_factor = 0.0)
    if isgpu(data)
        use_gpu = true
    elseif use_gpu && !isgpu(data)
        data = to_gpu(data)
    end

    pbc = ParamBitCircuit(pc, data)
    if use_gpu
        pbc = to_gpu(pbc)
    end

    reuse_v, reuse_f = nothing, nothing
    reuse_counts = use_gpu ? (nothing, nothing) : (nothing, nothing, nothing, nothing, nothing)

    for idx = 1 : length(data)
        pbc, reuse_v, reuse_f, reuse_counts = estimate_parameters_em(pbc, data[idx]; pseudocount, entropy_reg,
                                                                     use_gpu, reuse_v, reuse_f, reuse_counts,
                                                                     exp_update_factor)
    end

    estimate_parameters_cached!(pc, pbc)

    pbc.params
end
function estimate_parameters_em(pbc::ParamBitCircuit, data; pseudocount::Float64, entropy_reg::Float64 = 0.0,
                                use_sample_weights::Bool = true, use_gpu::Bool = isgpu(data),
                                reuse_v = nothing, reuse_f = nothing, reuse_counts = nothing,
                                exp_update_factor = 0.0
                               )
    if isweighted(data)
        # `data' is weighted according to its `weight' column
        data, weights = split_sample_weights(data)
    else
        use_sample_weights = false
    end

    if isgpu(data)
        use_gpu = true
    elseif use_gpu && !isgpu(data)
        data = to_gpu(data)
    end

    if reuse_counts === nothing
        reuse_counts = use_gpu ? (nothing, nothing) : (nothing, nothing, nothing, nothing, nothing)
    end

    params, v, f, reuse_counts = if use_gpu
        if !isgpu(pbc)
            pbc.bitcircuit = to_gpu(pbc.bitcircuit)
            pbc.params = to_gpu(pbc.params)
        end
        if use_sample_weights
            estimate_parameters_gpu(pbc, data, pseudocount; weights, reuse = true, reuse_v, reuse_f, reuse_counts, entropy_reg)
        else
            estimate_parameters_gpu(pbc, data, pseudocount; reuse = true, reuse_v, reuse_f, reuse_counts, entropy_reg)
        end
    else
        if use_sample_weights
            estimate_parameters_cpu(pbc, data, pseudocount; weights, reuse = true, reuse_v, reuse_f, reuse_counts, entropy_reg)
        else
            estimate_parameters_cpu(pbc, data, pseudocount; reuse = true, reuse_v, reuse_f, reuse_counts, entropy_reg)
        end
    end

    # Update the parameters to `pbc`
    if use_gpu # GPU
        tempparam = Vector{Float64}(undef, length(params))
        tempparam .= to_cpu(params)
        @inbounds @views pbc.params .+= log(exp_update_factor)
        @inbounds @views params .+= log(1.0 - exp_update_factor)
        delta = @inbounds @views @. CUDA.ifelse(pbc.params == params, CUDA.zero(params), CUDA.abs(pbc.params - params))
        @inbounds @views @. pbc.params = CUDA.max(pbc.params, params) + CUDA.log1p(CUDA.exp(-delta))

        CUDA.unsafe_free!(params)
        CUDA.unsafe_free!(delta)
    else # CPU
        @inbounds @views pbc.params = logaddexp.(pbc.params .+ log(exp_update_factor), params .+ log(1.0 - exp_update_factor))
    end

    pbc, v, f, reuse_counts
end

function estimate_parameters_cpu(pbc::ParamBitCircuit, data, pseudocount; weights = nothing, reuse::Bool = false,
                                 reuse_v = nothing, reuse_f = nothing,
                                 reuse_counts = (nothing, nothing, nothing, nothing, nothing), entropy_reg::Float64 = 0.0)
    # `data` is batched?
    data_batched = isbatched(data)

    # no need to synchronize, since each computation is unique to a decision node
    node_counts::Vector{Float64} = similar!(reuse_counts[1], Vector{Float64}, num_nodes(pbc.bitcircuit))
    @inbounds @views @. @avx node_counts .= typemin(eltype(node_counts)) # log(0)

    # rescale pseudocount using the average weight of samples
    if weights !== nothing
        if data_batched
            pseudocount = pseudocount * mapreduce(sum, +, weights) / num_examples(data)
        else
            pseudocount = pseudocount * sum(weights) / size(weights, 1)
        end
    end

    bc = pbc.bitcircuit
    params = pbc.params

    edge_counts::Vector{Float64} = similar!(reuse_counts[2], Vector{Float64}, num_elements(bc))
    parent_node_counts::Vector{Float64} = similar!(reuse_counts[3], Vector{Float64}, num_elements(bc))
    @inbounds @views edge_counts[:] .= zero(Float64)
    @inbounds @views parent_node_counts[:] .= zero(Float64)

    # Buffer to save some allocations
    buffer::Vector{Float64} = similar!(reuse_counts[4], Vector{Float64},
        data_batched ? num_examples(data[1]) : num_examples(data))

    # Pre-compute log(weights) to save computation time
    if weights !== nothing
        log_weights = similar!(reuse_counts[5], Vector{Float64},
            data_batched ? num_examples(data[1]) : num_examples(data))

        if !data_batched
            @inbounds @views @avx log_weights .= log.(weights)
        end
    else
        log_weights = nothing
    end

    # For batched dataset, we want to enable 'estimate' for parent_node_counts only in the last minibatch
    estimate_flag = true

    @inline function on_node(flows, values, dec_id, weights::Nothing)
        @inbounds @views @. @avx buffer = flows[:, dec_id]
        node_counts[dec_id] = logaddexp(node_counts[dec_id], logsumexp(buffer))
    end
    @inline function on_node(flows, values, dec_id, weights)
        @inbounds @views @. @avx buffer = flows[:, dec_id] + log_weights[:]
        node_counts[dec_id] = logaddexp(node_counts[dec_id], logsumexp(buffer))
    end

    @inline function estimate(element, decision, edge_count)
        edge_counts[element] += exp(edge_count)
        if estimate_flag # For batched dataset, we only accumulate parent_node_counts after all node_counts have been cumulated
            parent_node_counts[element] += exp(node_counts[decision])
        end
    end

    @inline function on_edge(flows, values, prime, sub, element, grandpa, single_child, weights::Nothing)
        θ = eltype(flows)(params[element])
        if !single_child
            @inbounds @views @. @avx buffer = values[:, prime] + values[:, sub] - values[:, grandpa] + flows[:, grandpa] + θ
            buffer .= ifelse.(isnan.(buffer[:]), typemin(eltype(flows)), buffer)

            edge_count = logsumexp(buffer)

            estimate(element, grandpa, edge_count)
        end # no need to estimate single child params, they are always prob 1
    end
    @inline function on_edge(flows, values, prime, sub, element, grandpa, single_child, weights)
        θ = eltype(flows)(params[element])
        if !single_child
            @inbounds @views @. @avx buffer = values[:, prime] + values[:, sub] - values[:, grandpa] + flows[:, grandpa] + θ + log_weights
            buffer .= ifelse.(isnan.(buffer), typemin(eltype(flows)), buffer)

            edge_count = logsumexp(buffer)

            estimate(element, grandpa, edge_count)
        end # no need to estimate single child params, they are always prob 1
    end

    if data_batched
        if weights !== nothing
            v, f = reuse_v, reuse_f
            for idx = 1 : length(data)
                d = data[idx]
                w = weights[idx]

                # Resize buffer if the current minibatch has a different size
                if size(buffer, 1) != num_examples(d)
                    resize!(buffer, num_examples(d))
                end

                # Resize log_weights if the current minibatch has a different size
                if data_batched
                    if size(log_weights, 1) != num_examples(d)
                        resize!(log_weights, num_examples(d))
                    end
                    @inbounds @views @avx log_weights .= log.(w)
                end

                # Only accumulate parent_node_counts in the final call to marginal_flows
                if idx == length(data)
                    estimate_flag = true
                else
                    estimate_flag = false
                end

                v, f = marginal_flows(pbc, d, v, f; on_node = on_node, on_edge = on_edge, weights = w)
            end
        else
            v, f = reuse_v, reuse_f
            for idx = 1 : length(data)
                d = data[idx]

                # Resize buffer if the current minibatch has a different size
                if size(buffer, 1) != num_examples(d)
                    resize!(buffer, num_examples(d))
                end

                if idx == length(data)
                    estimate_flag = true
                else
                    estimate_flag = false
                end

                v, f = marginal_flows(pbc, d, v, f; on_node = on_node, on_edge = on_edge, weights = nothing)
            end
        end
    else
        v, f = marginal_flows(pbc, data, reuse_v, reuse_f; on_node, on_edge, weights)
    end

    if entropy_reg > 1e-8
        total_data_counts = (weights === nothing) ? Float64(num_examples(data)) : mapreduce(sum, +, weights)
        edge_counts = apply_entropy_reg_cpu(bc; log_params = edge_counts, edge_counts, total_data_counts,
                                            pseudocount, entropy_reg)
    else
        # `edge_counts` now becomes "params"
        @simd for i = 1 : num_elements(bc)
            num_els = num_elements(bc.nodes, bc.elements[1, i])
            if num_els == 1
                @inbounds edge_counts[i] = zero(eltype(edge_counts)) # log(1)
            else
                @inbounds edge_counts[i] = log((edge_counts[i] + pseudocount / num_elements(bc.nodes, bc.elements[1, i])) / (parent_node_counts[i] + pseudocount))
            end
        end
    end

    if reuse
        # Also return the allocated vars v, f, counts for future reuse
        edge_counts, v, f, (node_counts, edge_counts, parent_node_counts, buffer, log_weights)
    else
        edge_counts # a.k.a. log_probs
    end
end

function estimate_parameters_gpu(pbc::ParamBitCircuit, data, pseudocount; weights = nothing, reuse::Bool = false,
                                 reuse_v = nothing, reuse_f = nothing, reuse_counts = (nothing, nothing),
                                 entropy_reg::Float64 = 0.0)
    # `data` is batched?
    data_batched = isbatched(data)

    # rescale pseudocount using the average weight of samples
    if weights !== nothing
        if data_batched
            if isgpu(weights)
                pseudocount = pseudocount * mapreduce(sum, +, to_cpu(weights)) / num_examples(data)
            else
                pseudocount = pseudocount * mapreduce(sum, +, weights) / num_examples(data)
            end
        else
            pseudocount = pseudocount * sum(weights) / size(weights, 1)
        end
    end

    bc = pbc.bitcircuit
    node_counts::CuVector{Float64} = similar!(reuse_counts[1], CuVector{Float64}, num_nodes(bc))
    edge_counts::CuVector{Float64} = similar!(reuse_counts[2], CuVector{Float64}, num_elements(bc))
    @inbounds @views node_counts[:] .= zero(Float64)
    @inbounds @views edge_counts[:] .= zero(Float64)
    # need to manually cudaconvert closure variables
    node_counts_device = CUDA.cudaconvert(node_counts)
    edge_counts_device = CUDA.cudaconvert(edge_counts)

    @inline function on_node(flows, values, dec_id, chunk_id, flow, weight::Nothing)
        c::Float64 = exp(flow) # cast for @atomic to be happy
        CUDA.@atomic node_counts_device[dec_id] += c
    end
    @inline function on_node(flows, values, dec_id, chunk_id, flow, weight::AbstractFloat)
        c::Float64 = exp(flow) * weight # cast for @atomic to be happy
        CUDA.@atomic node_counts_device[dec_id] += c
    end

    @inline function on_edge(flows, values, prime, sub, element, grandpa, chunk_id, edge_flow, single_child, weight::Nothing)
        if !single_child
            c::Float64 = exp(edge_flow) # cast for @atomic to be happy
            CUDA.@atomic edge_counts_device[element] += c
        end
    end
    @inline function on_edge(flows, values, prime, sub, element, grandpa, chunk_id, edge_flow, single_child, weight::AbstractFloat)
        if !single_child
            c::Float64 = exp(edge_flow) * weight # cast for @atomic to be happy
            CUDA.@atomic edge_counts_device[element] += c
        end
    end

    if data_batched
        v, f = reuse_v, reuse_f
        if weights != nothing
            map(zip(data, weights)) do (d, w)
                if w != nothing
                    w = to_gpu(w)
                end
                v, f = marginal_flows(pbc, d, v, f; on_node = on_node, on_edge = on_edge, weights = w)

                nothing # Return nothing to save some time
            end
        else
            map(data) do d
                v, f = marginal_flows(pbc, d, v, f; on_node = on_node, on_edge = on_edge, weights = nothing)

                nothing # Return nothing to save some time
            end

            nothing # Return nothing to save some time
        end
    else
        if weights != nothing
            weights = to_gpu(weights)
        end

        v, f = marginal_flows(pbc, data, reuse_v, reuse_f; on_node, on_edge, weights)
    end

    if entropy_reg > 1e-8
        total_data_counts = (weights === nothing) ? Float64(num_examples(data)) : mapreduce(sum, +, weights)
        params = apply_entropy_reg_gpu(bc; log_params = edge_counts, edge_counts, total_data_counts, pseudocount, entropy_reg)
    else
        # TODO: reinstate simpler implementation once https://github.com/JuliaGPU/GPUArrays.jl/issues/313 is fixed and released
        @inbounds parents = bc.elements[1,:]
        @inbounds parent_counts = node_counts[parents]
        @inbounds parent_elcount = bc.nodes[2,parents] .- bc.nodes[1,parents] .+ 1
        params = log.((edge_counts .+ (pseudocount ./ parent_elcount))
                        ./ (parent_counts .+ pseudocount))
        params = @inbounds @views @. ifelse(parent_elcount == 1, zero(params), params)

        CUDA.unsafe_free!(parents)
        CUDA.unsafe_free!(parent_counts)
        CUDA.unsafe_free!(parent_elcount)
    end

    # Only free the memory if the reuse memory is not provided
    if !reuse
        CUDA.unsafe_free!(v) # save the GC some effort
        CUDA.unsafe_free!(f) # save the GC some effort
        CUDA.unsafe_free!(node_counts) # save the GC some effort
        CUDA.unsafe_free!(edge_counts) # save the GC some effort
    end

    if reuse
        # Also return the allocated vars v, f, counts for future reuse
        if !isgpu(params)
            to_gpu(params), v, f, (node_counts, edge_counts)
        else
            params, v, f, (node_counts, edge_counts)
        end
    else
        to_cpu(params) # a.k.a. log_probs
    end
end

"""
Add entropy regularization to a deterministic (see LogicCircuits.isdeterministic) probabilistic
circuit. `alpha` is a hyperparameter that balances the weights between the likelihood and the
entropy: \argmax_{\theta} L(\theta) = ll_mean(\theta) + alpha * entropy(\theta).
"""
function apply_entropy_reg_cpu(bc::BitCircuit; log_params::Vector{Float64}, edge_counts::Vector{Float64},
                               total_data_counts::Float64, pseudocount::Float64 = 0.1,
                               entropy_reg::Float64 = 0.0)::Vector{Float64}
    a = Vector{Float64}(undef, 0) # For reuse
    log_probs = Vector{Float64}(undef, 0) # For reuse
    p_exp_logprob = Vector{Float64}(undef, 0) # For reuse

    node_log_probs = Vector{Float64}(undef, num_nodes(bc)) # Probability of each decision/sum node
    @inbounds @views node_log_probs .= -Inf
    @inbounds node_log_probs[end] = 0.0

    for i = num_nodes(bc) : -1 : 1 # Traverse the circuit top-down
        @inbounds child_ele_start = bc.nodes[1,i]
        @inbounds child_ele_end = bc.nodes[2,i]
        num_eles = child_ele_end - child_ele_start + 1

        if num_eles > 1 # No need to do anything if the decision node has 1 child.
            resize!(a, num_eles)
            resize!(log_probs, num_eles)
            resize!(p_exp_logprob, num_eles)

            @inbounds @views a .= edge_counts[child_ele_start: child_ele_end] .+ (pseudocount / num_eles)
            @inbounds b = entropy_reg * exp(node_log_probs[i] + log(total_data_counts + pseudocount))

            @inbounds @views log_probs .= log.(edge_counts[child_ele_start: child_ele_end] .+ (pseudocount / num_eles))
            @inbounds @views log_probs .-= logsumexp(log_probs)

            for _ = 1 : 3
                y = sum(b .* log_probs .- a .* exp.(-log_probs)) / num_eles
                for _ = 1 : 4
                    @inbounds @views p_exp_logprob .= a .* exp.(-log_probs)
                    @inbounds @views log_probs .+= (p_exp_logprob .- b .* log_probs .+ y) ./ (p_exp_logprob .+ b .+ 1e-4)
                    @inbounds @views log_probs .-= logsumexp(log_probs) # Project the parameters back to its valid space
                end
            end

            # Update parameters
            @inbounds @views log_params[child_ele_start: child_ele_end] .= log_probs
        end

        # Update child nodes' log_prob
        if child_ele_start >= 1
            for child_ele_id = child_ele_start : child_ele_end
                prime_idx = bc.elements[2, child_ele_id]
                sub_idx = bc.elements[3, child_ele_id]
                node_log_probs[prime_idx] = logaddexp(node_log_probs[prime_idx], node_log_probs[i] + log_params[child_ele_id])
                node_log_probs[sub_idx] = logaddexp(node_log_probs[sub_idx], node_log_probs[i] + log_params[child_ele_id])
            end
        end
    end

    log_params
end
function apply_entropy_reg_gpu(bc::BitCircuit; log_params, edge_counts,
                               total_data_counts::Float64, pseudocount::Float64 = 0.1, entropy_reg::Float64 = 0.0)
    # TODO: if this is performance-critical, consider change to native GPU code
    log_params = to_cpu(log_params)
    edge_counts = to_cpu(edge_counts)
    bc = to_cpu(bc)

    a = Vector{Float64}(undef, 0) # For reuse
    log_probs = Vector{Float64}(undef, 0) # For reuse
    p_exp_logprob = Vector{Float64}(undef, 0) # For reuse

    node_log_probs = Vector{Float64}(undef, num_nodes(bc)) # Probability of each decision/sum node
    @inbounds @views node_log_probs .= -Inf
    @inbounds node_log_probs[end] = 0.0

    for i = num_nodes(bc) : -1 : 1 # Traverse the circuit top-down
        @inbounds child_ele_start = bc.nodes[1,i]
        @inbounds child_ele_end = bc.nodes[2,i]
        num_eles = child_ele_end - child_ele_start + 1

        if num_eles > 1 # No need to do anything if the decision node has 1 child.
            resize!(a, num_eles)
            resize!(log_probs, num_eles)
            resize!(p_exp_logprob, num_eles)

            @inbounds @views a .= edge_counts[child_ele_start: child_ele_end] .+ (pseudocount / num_eles)
            @inbounds b = entropy_reg * exp(node_log_probs[i] + log(total_data_counts + pseudocount))

            @inbounds @views log_probs .= log.(edge_counts[child_ele_start: child_ele_end] .+ (pseudocount / num_eles))
            @inbounds @views log_probs .-= logsumexp(log_probs)

            for _ = 1 : 3
                y = sum(b .* log_probs .- a .* exp.(-log_probs)) / num_eles
                for _ = 1 : 4
                    @inbounds @views p_exp_logprob .= a .* exp.(-log_probs)
                    @inbounds @views log_probs .+= (p_exp_logprob .- b .* log_probs .+ y) ./ (p_exp_logprob .+ b .+ 1e-4)
                    @inbounds @views log_probs .-= logsumexp(log_probs) # Project the parameters back to its valid space
                end
            end

            # Update parameters
            @inbounds @views log_params[child_ele_start: child_ele_end] .= log_probs
        end

        # Update child nodes' log_prob
        if child_ele_start >= 1
            for child_ele_id = child_ele_start : child_ele_end
                prime_idx = bc.elements[2, child_ele_id]
                sub_idx = bc.elements[3, child_ele_id]
                node_log_probs[prime_idx] = logaddexp(node_log_probs[prime_idx], node_log_probs[i] + log_params[child_ele_id])
                node_log_probs[sub_idx] = logaddexp(node_log_probs[sub_idx], node_log_probs[i] + log_params[child_ele_id])
            end
        end
    end

    log_params
end


text:
"""
    MetropolisHastings{D}

`MetropolisHastings` has one field, `proposal`.
`proposal` is a `Proposal`, `NamedTuple` of `Proposal`, or `Array{Proposal}` in the shape of your data.
For example, if you wanted the sampler to return a `NamedTuple` with shape

```julia
x = (a = 1.0, b=3.8)
```

The proposal would be

```julia
proposal = (a=StaticProposal(Normal(0,1)), b=StaticProposal(Normal(0,1)))
````

Other allowed proposals are

```
p1 = StaticProposal(Normal(0,1))
p2 = StaticProposal([Normal(0,1), InverseGamma(2,3)])
p3 = StaticProposal((a=Normal(0,1), b=InverseGamma(2,3)))
p4 = StaticProposal((x=1.0) -> Normal(x, 1))
```

The sampler is constructed using

```julia
spl = MetropolisHastings(proposal)
```

When using `MetropolisHastings` with the function `sample`, the following keyword
arguments are allowed:

- `init_params` defines the initial parameterization for your model. If
none is given, the initial parameters will be drawn from the sampler's proposals.
- `param_names` is a vector of strings to be assigned to parameters. This is only
used if `chain_type=Chains`.
- `chain_type` is the type of chain you would like returned to you. Supported
types are `chain_type=Chains` if `MCMCChains` is imported, or
`chain_type=StructArray` if `StructArrays` is imported.
"""
struct MetropolisHastings{D} <: MHSampler
    proposal::D
end

StaticMH(d) = MetropolisHastings(StaticProposal(d))
RWMH(d) = MetropolisHastings(RandomWalkProposal(d))

# default function without RNG
propose(spl::MetropolisHastings, args...) = propose(Random.GLOBAL_RNG, spl, args...)

# Propose from a vector of proposals
function propose(
    rng::Random.AbstractRNG,
    spl::MetropolisHastings{<:AbstractArray},
    model::DensityModel
)
    proposal = map(p -> propose(rng, p, model), spl.proposal)
    return Transition(model, proposal)
end

function propose(
    rng::Random.AbstractRNG,
    spl::MetropolisHastings{<:AbstractArray},
    model::DensityModel,
    params_prev::Transition
)
    proposal = map(spl.proposal, params_prev.params) do p, params
        propose(rng, p, model, params)
    end
    return Transition(model, proposal)
end

# Make a proposal from one Proposal struct.
function propose(
    rng::Random.AbstractRNG,
    spl::MetropolisHastings{<:Proposal},
    model::DensityModel
)
    proposal = propose(rng, spl.proposal, model)
    return Transition(model, proposal)
end

function propose(
    rng::Random.AbstractRNG,
    spl::MetropolisHastings{<:Proposal},
    model::DensityModel,
    params_prev::Transition
)
    proposal = propose(rng, spl.proposal, model, params_prev.params)
    return Transition(model, proposal)
end

# Make a proposal from a NamedTuple of Proposal.
function propose(
    rng::Random.AbstractRNG,
    spl::MetropolisHastings{<:NamedTuple},
    model::DensityModel
)
    proposal = _propose(rng, spl.proposal, model)
    return Transition(model, proposal)
end

function propose(
    rng::Random.AbstractRNG,
    spl::MetropolisHastings{<:NamedTuple},
    model::DensityModel,
    params_prev::Transition
)
    proposal = _propose(rng, spl.proposal, model, params_prev.params)
    return Transition(model, proposal)
end

@generated function _propose(
    rng::Random.AbstractRNG,
    proposal::NamedTuple{names},
    model::DensityModel
) where {names}
    isempty(names) && return :(NamedTuple())
    expr = Expr(:tuple)
    expr.args = Any[:($name = propose(rng, proposal.$name, model)) for name in names]
    return expr
end

@generated function _propose(
    rng::Random.AbstractRNG,
    proposal::NamedTuple{names},
    model::DensityModel,
    params_prev::NamedTuple
) where {names}
    isempty(names) && return :(NamedTuple())
    expr = Expr(:tuple)
    expr.args = Any[
        :($name = propose(rng, proposal.$name, model, params_prev.$name)) for name in names
    ]
    return expr
end

transition(sampler, model, params) = transition(model, params)
transition(model, params) = Transition(model, params)

# Called to update proposal when the first sample is drawn
trackstep!(proposal::Proposal, params) = nothing
trackstep!(proposal::AbstractArray, params) = foreach(trackstep!, proposal, params)
trackstep!(proposal::NamedTuple, params) = foreach(trackstep!, proposal, params)

# Called to update proposal when a new step is performed
# The last argument determines if the step is an acceptance step
trackstep!(proposal::Proposal, params,
           ::Union{Val{true}, Val{false}}) = nothing

trackstep!(proposal::AbstractArray, params, accept::Union{Val{true},Val{false}}) =
           foreach((prop, par) -> trackstep!(prop, par, accept), proposal, params)

trackstep!(proposal::NamedTuple, params, accept::Union{Val{true},Val{false}}) =
           foreach((prop, par) -> trackstep!(prop, par, accept), proposal, params)

# Define the first sampling step.
# Return a 2-tuple consisting of the initial sample and the initial state.
# In this case they are identical.
function AbstractMCMC.step(
    rng::Random.AbstractRNG,
    model::DensityModel,
    spl::MHSampler;
    init_params=nothing,
    kwargs...
)
    if init_params === nothing
        transition = propose(rng, spl, model)
    else
        transition = AdvancedMH.transition(spl, model, init_params)
    end

    trackstep!(spl.proposal, transition.params)
    return transition, transition
end

# Define the other sampling steps.
# Return a 2-tuple consisting of the next sample and the the next state.
# In this case they are identical, and either a new proposal (if accepted)
# or the previous proposal (if not accepted).
function AbstractMCMC.step(
    rng::Random.AbstractRNG,
    model::DensityModel,
    spl::MHSampler,
    params_prev::AbstractTransition;
    kwargs...
)
    # Generate a new proposal.
    params = propose(rng, spl, model, params_prev)

    # Calculate the log acceptance probability.
    logα = logdensity(model, params) - logdensity(model, params_prev) +
        logratio_proposal_density(spl, params_prev, params)

    # Decide whether to return the previous params or the new one.
    if -Random.randexp(rng) < logα
        trackstep!(spl.proposal, params.params, Val(true))
        return params, params
    else
        trackstep!(spl.proposal, params_prev.params, Val(false))
        return params_prev, params_prev
    end
end

function logratio_proposal_density(
    sampler::MetropolisHastings, params_prev::Transition, params::Transition
)
    return logratio_proposal_density(sampler.proposal, params_prev.params, params.params)
end


text:
module PseudoArcLengthContinuation
	using Plots, JLD2, Printf, Dates, LinearMaps, BlockArrays
	using Setfield: setproperties, @set
	using Parameters: @with_kw, @unpack

	include("BorderedArrays.jl")
	include("LinearSolver.jl")
	include("EigSolver.jl")
	include("LinearBorderSolver.jl")
	include("Newton.jl")
	include("Continuation.jl")
	include("Bifurcations.jl")
	include("Predictor.jl")
	include("DeflationOperator.jl")
	include("BorderedProblem.jl")
	include("Plotting.jl")
	include("Utils.jl")
	include("FoldCont.jl")
	include("HopfCont.jl")

	include("periodicorbit/PeriodicOrbitUtils.jl")
	include("periodicorbit/Shooting.jl")
	include("periodicorbit/PeriodicOrbitFD.jl")
	include("periodicorbit/FloquetQaD.jl")

	export	DefaultLS, GMRESIterativeSolvers, GMRESIterativeSolvers!, GMRESKrylovKit,
			DefaultEig, EigArpack, EigIterativeSolvers, EigKrylovKit, geteigenvector, AbstractEigenSolver
	export	BorderedProblem, JacobianBorderedProblem, LinearSolverBorderedProblem
	export	MatrixBLS, BorderingBLS, MatrixFreeBLS, BorderedArray
	export	DeflationOperator, DeflatedProblem, DeflatedLinearSolver, scalardM
	export	SecantPred, BorderedPred, NaturalPred
	export	NewtonPar, newton, newtonDeflated, newtonPALC, newtonFold, newtonHopf, newtonBordered
	export	ContinuationPar, ContResult, continuation, continuation!, continuationFold, continuationHopf, continuationPOTrap, continuationBordered
	export	PALCIterable, iterate, PALCStateVariables, solution, getu, getp
	export	FoldPoint, FoldProblemMinimallyAugmented, FoldLinearSolveMinAug, foldPoint
	export	HopfPoint, HopfProblemMinimallyAugmented, HopfLinearSolveMinAug
	export	PeriodicOrbitTrapProblem, continuationPOTrap
	export	Flow, ShootingProblem, PoincareShootingProblem, continuationPOShooting, getPeriod, AbstractShootingProblem, extractPeriodShooting
	export	FloquetQaDTrap, FloquetQaDShooting
	export	plotBranch, plotBranch!, guessFromHopf
end


text:
# Autogenerated wrapper script for unixODBC_jll for x86_64-w64-mingw32
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-2.dll")
JLLWrappers.@declare_library_product(libodbccr, "libodbccr-2.dll")
JLLWrappers.@declare_library_product(libodbcinst, "libodbcinst-2.dll")
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,
        "bin\\libodbc-2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libodbccr,
        "bin\\libodbccr-2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libodbcinst,
        "bin\\libodbcinst-2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        dltest,
        "bin\\dltest.exe",
    )

    JLLWrappers.@init_executable_product(
        isql,
        "bin\\isql.exe",
    )

    JLLWrappers.@init_executable_product(
        iusql,
        "bin\\iusql.exe",
    )

    JLLWrappers.@init_executable_product(
        odbc_config,
        "bin\\odbc_config.exe",
    )

    JLLWrappers.@init_executable_product(
        odbcinst,
        "bin\\odbcinst.exe",
    )

    JLLWrappers.@init_executable_product(
        slencheck,
        "bin\\slencheck.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
weighted sampling: takes (unnormalized) log probs and returns index
"""
function wsample_logprob(logprobs::AbstractVector{<:Real})
  logdenom = logsumexp(logprobs)
  p = exp.(logprobs .- logdenom)
  return Distributions.wsample(p)
end
function wsample_logprob(logprobs::AbstractVector{<:Real}, n::Int)
  return [wsample_logprob(logprobs) for _ in 1:n]
end

Random.rand(rng::Random.AbstractRNG, d::Array{T}) where {T <: Distribution} = rand.(d)

function update_mean!(current_mean, x, iter)
  current_mean .+= (x - current_mean) / iter
end

function update_cov!(current_cov, current_mean, x, iter)
  d = x - current_mean
  current_cov .= current_cov * (iter - 1)/iter + (d*d') * (iter - 1)/iter^2
  current_cov .= Matrix(LinearAlgebra.Symmetric(current_cov))
end

eye(T::Type, n::Integer) = Matrix{T}(LinearAlgebra.I(n))
eye(n::Integer) = eye(Float64, n)

# Quantiles for arbitrary array.
function quantiles(X, q; dims, drop=false)
  Q = mapslices(x -> quantile(x, q), X, dims=dims)
  out = drop ? dropdims(Q, dims=dims) : Q
  return out
end


text:
#=
Consider all integer combinations of ab for 2 ≤ a ≤ 5 and 2 ≤ b ≤ 5:

22=4, 23=8, 24=16, 25=32
32=9, 33=27, 34=81, 35=243
42=16, 43=64, 44=256, 45=1024
52=25, 53=125, 54=625, 55=3125
If they are then placed in numerical order, with any repeats removed, we get the following sequence of 15 distinct terms:

4, 8, 9, 16, 25, 27, 32, 64, 81, 125, 243, 256, 625, 1024, 3125

How many distinct terms are in the sequence generated by ab for 2 ≤ a ≤ 100 and 2 ≤ b ≤ 100?
=#
function calc(up)
  list = BigInt[]
  for a in 2:up, b in 2:up
    c = BigInt(a)^BigInt(b)
    push!(list, c)
  end
  length(sort(unique(list)))
end
@time println(calc(100))

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
@i function ⊖(exp)(out!::GVar, x::GVar{T}) where T
    value(out!) -= exp(value(x))
    @anc anc1 = zero(T)
    anc1 += exp(value(x))
    grad(x) += grad(out!) * anc1
    anc1 -= exp(value(x))
end

@i function ⊖(log)(out!::GVar, x::GVar{T}) where T
    value(out!) -= log(value(x))
    grad(x) += grad(out!) / value(x)
end

@i function ⊖(sin)(out!::GVar, x::GVar{T}) where T
    value(out!) -= sin(value(x))
    @anc anc1 = zero(T)
    anc1 += cos(value(x))
    grad(x) += grad(out!) * anc1
    anc1 -= cos(value(x))
end

@i function ⊖(cos)(out!::GVar, x::GVar{T}) where T
    value(out!) -= cos(value(x))
    @anc anc1 = zero(T)
    anc1 -= sin(value(x))
    grad(x) += grad(out!) * anc1
    anc1 += sin(value(x))
end

for op in [:exp, :log, :sin, :cos]
    @eval @nograd ⊖($op)(out!, x::GVar)
    @eval @nograd ⊖($op)(out!::GVar, x)
end

@i function IROT(a!::GVar, b!::GVar, θ::GVar)
    IROT(value(a!), value(b!), value(θ))
    NEG(value(θ))
    value(θ) ⊖ π/2
    ROT(grad(a!), grad(b!), value(θ))
    grad(θ) += value(a!) * grad(a!)
    grad(θ) += value(b!) * grad(b!)
    value(θ) ⊕ π/2
    NEG(value(θ))
    ROT(grad(a!), grad(b!), π/2)
end

@i function IROT(a!::GVar, b!::GVar, θ)
    IROT(value(a!), value(b!), θ)
    NEG(θ)
    θ ⊖ π/2
    ROT(grad(a!), grad(b!), θ)
    θ ⊕ π/2
    NEG(θ)
    ROT(grad(a!), grad(b!), π/2)
end

@nograd IROT(a!, b!, θ::GVar)