text:
@inline eftRound{T<:AbstractFloat}(hi::T, lo::T, ::RoundingMode{:ToZero}) =
     (((signbit(hi)==signbit(lo)) | (lo==zero(T))) ? hi : nextNearerToZero(hi))

@inline eftRound{T<:AbstractFloat}(hi::T, lo::T, ::RoundingMode{:FromZero}) =
    (((signbit(hi)!=signbit(lo)) | (lo == zero(T))) ? hi : nextAwayFromZero(hi))

@inline eftRound{T<:AbstractFloat}(hi::T, lo::T, ::RoundingMode{:Up}) =
    (((signbit(lo) | (lo == zero(T)))) ? hi : nextFloat(hi))

@inline eftRound{T<:AbstractFloat}(hi::T, lo::T, ::RoundingMode{:Down}) =
    (((signbit(lo) & (lo != zero(T)))) ? prevFloat(hi) : hi)

@inline eftRound{T<:AbstractFloat}(hi::T, lo::T, ::RoundingMode{:Nearest}) =
    (hi)

ext: jl


text:
using MagneticReadHead
using Test
using OrderedCollections

foo(a) = 1
foo(x,y) = 2
foo(x,y, zs...) = 3

@testset "subnames" begin
    name2arg =Dict([:x=>1, :y=>2])
    @test MagneticReadHead.subnames(name2arg, :x) == 1
    @test MagneticReadHead.subnames(name2arg, 15) == 15

    @test MagneticReadHead.subnames(name2arg, :(f(x)+y)) ==:(f(1)+2)

end


text:
# This file defines the SimulationReport type. This type holds the reports for
#   a particular time resolution of a particular simulation.

# The SimulationReport type does not require any other times.
requiredTypes = []

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

export SimulationReport
"""
This type defines a set of reports of a `ManpowerSimulation` for a specific time
grid.

The type contains the following fields:
* `timeGrid::Vector{Float64}`: the time grid of the report. The entries in this
vector are sorted from smallest to largest, and there are no duplicate entries.
* `activeCount::Vector{Int}`: the number of active personnel in the simulation
at each time point of the time grid.
* `fluxIn::Vector{Int}`: the number of personnel that became active between
every pair of consecutive time steps, with the final time point included.
* `fluxOut::Vector{Int}`: the number of personnel that became inactive between
every pair of consecutive time steps, with the final time point included.
* `fluxOutBreakdown::Dict{String, Vector{Int}}`: the number of personnel that
became inactive between every pair of consecutive time steps, with the final
time point included, broken down by reason for becoming inactive.
* `ageGrid::Vector{Float64}`: the grid of ages for the age distribution report.
* `ageDistribution::Array{Int, 2}`: the distribution of the age of active
personnel at every time point of the time grid, where each row corresponds to
one time point.
* `ageStats::Array{Float64, 5}`: basic statistics about the age distribution of
the active personnel at every time point of the time grid. The five columns hold
the following information, in this order: mean age, standard deviation, median
age, minimum age, maximum age.
"""
type SimulationReport

    timeGrid::Vector{Float64}
    activeCount::Vector{Int}
    fluxIn::Vector{Int}
    fluxOut::Vector{Int}
    fluxOutBreakdown::Dict{String, Vector{Int}}
    ageGrid::Vector{Float64}
    ageDist::Array{Int, 2}
    ageStats::Array{Float64, 2}

    # Default constructor
    function SimulationReport()

        newSimRep = new()
        newSimRep.timeGrid = Vector{Float64}()
        newSimRep.activeCount = Vector{Int}()
        newSimRep.fluxIn = Vector{Int}()
        newSimRep.fluxOut = Vector{Int}()
        newSimRep.fluxOutBreakdown = Dict{String, Vector{Int}}()
        newSimRep.ageGrid = Vector{Float64}()
        newSimRep.ageDist = Array{Int}( 0, 0 )
        newSimRep.ageStats = Array{Float64}( 0, 0 )
        return newSimRep

    end  # SimulationReport()

    # Constructor on grid
    function SimulationReport( timeGrid::Vector{Float64} )

        newSimRep = SimulationReport()
        setTimeGrid( newSimRep, timeGrid )
        return newSimRep

    end  # SimulationReport( timeGrid )

end  # type SimulationReport


text:
filename() = "day14.txt"

#=
Read data into tuples of commands.

Example:
julia> read_data()
4-element Array{NamedTuple,1}:
 (command = :setmask, value = "XXXXXXXXXXXXXXXXXXXXXXXXXXXXX1XXXX0X",
        setter = [2 => 0, 7 => 1])
 (command = :assign, address = "8", value = 11)
 (command = :assign, address = "7", value = 101)
 (command = :assign, address = "8", value = 0)
=#
function read_data()
    function parse_command(line)
        left, right = split(line, " = ")
        if left == "mask"
            right_reversed = reverse(right)
            bit_positions = findall(in(['0','1']), right_reversed)
            bit_setter = [p => parse(Int, right_reversed[p]) for p in bit_positions]
            (command = :setmask, value = right, setter = bit_setter)
        else
            address = match(r"mem\[(.*)\]", left).captures[1]
            (command = :assign, address = address, value = parse(Int, right))
        end
    end
    lines = readlines(filename())
    return parse_command.(lines)
end

"""
Set a specific bit at `pos` position (counting from right)
for the value.  The `bit` may be 1 or 0. Returning the result
as a decimal value.
"""
function setbit(value::Int, pos::Int, bit)
    str = bitstring(value)
    rpos = 64 - pos + 1
    return parse(Int, string(str[1:rpos-1], "$bit", str[rpos+1:end]); base = 2)
end

"""
Apply bitmask to a value. The `setter` argument contains an array of
pairs. Each pair is interpreted as pos => bit.
"""
function apply_bitmask(val::Int, setter)
    for s in setter
        val = setbit(val, first(s), last(s))
    end
    return val
end

"""
Execute the program.
- For bitmask command, save a reference of the setter.
- For assign command, take the value and apply bitmasks, then assign to memory.
"""
function part1()
    program = read_data()
    local setter
    mem = Dict()
    for instruction in program
        if instruction.command == :setmask
            setter = instruction.setter
        else
            mem[instruction.address] = apply_bitmask(instruction.value, setter)
        end
    end
    return sum(values(mem))
end

"""
Read data for part 2. This is needed because part2 interprets the instructions
differently:
- '1': set bit
- '0': ignore
- 'X': floater

The `:setmask` instruction now contains a new `floater` property, which contains
the positions where the bit floats.
"""
function read_data2()
    function parse_command(line)
        left, right = split(line, " = ")
        if left == "mask"
            right_reversed = reverse(right)
            setter = [i => c for (i,c) in enumerate(right_reversed) if c == '1']
            floater = findall(==('X'), right_reversed)
            (command = :setmask, value = right, setter = setter, floater = floater)
        else
            address = match(r"mem\[(.*)\]", left).captures[1]
            (command = :assign, address = parse(Int, address), value = parse(Int, right))
        end
    end
    lines = readlines("day14.txt")
    return parse_command.(lines)
end

# Make a new setter
# Arguments:
# - `floater`: an array of bit positions
# - `tup`: tuple of bit combinations
make_setter(floater, tup) =[floater[i] => tup[i] for i in 1:length(floater)]

function part2()
    program = read_data2()
    local setter, floater
    mem = Dict{Int,Int}()
    for instruction in program
        if instruction.command == :setmask
            setter = instruction.setter
            floater = instruction.floater
        else
            temp_address = apply_bitmask(instruction.address, setter)
            for bit_tuple in Iterators.product([[0,1] for _ in 1:length(floater)]...)
                addr = apply_bitmask(temp_address, make_setter(floater, bit_tuple))
                mem[addr] = instruction.value
            end
        end
    end
    return sum(values(mem)), mem
end

# Visualization
using Plots, Statistics

function part2_animated(board_size = 261550, max_value = 969683639)
    anim = Animation()
    marker_size = 2.0
    marker_alpha = 0.8
    marker_color = :green
    marker_shape = :rect
    scatter([1],[1], xlims = (1, board_size), ylims = (1, board_size),
        markercolor = :black,  # imaginary
        markersize = marker_size,
        markeralpha = marker_alpha,
        markershape = marker_shape,
        markerstrokewidth = 0,
        legend = nothing, ticks = nothing, showaxis = false, grid = nothing,
        background_color = :black,
        foreground_color_axis = :black, gridlinewidth = 0)
    program = read_data2()
    local setter, floater
    mem = Dict{Int,Int}()
    cnt = 0
    xs = Int[]
    ys = Int[]
    vs = BigInt[]
    for instruction in program
        if instruction.command == :setmask
            setter = instruction.setter
            floater = instruction.floater
        else
            temp_address = apply_bitmask(instruction.address, setter)
            for bit_tuple in Iterators.product([[0,1] for _ in 1:length(floater)]...)
                addr = apply_bitmask(temp_address, make_setter(floater, bit_tuple))
                mem[addr] = instruction.value
                x, y = addr ÷ board_size, addr % board_size
                push!(xs, x)
                push!(ys, y)
                push!(vs, instruction.value)
            end
        end
        cnt += 1
        if cnt % 5 == 0
            scatter!(xs, ys,
                markercolor = marker_color,
                markersize = marker_size,
                markeralpha = marker_alpha,
                markershape = marker_shape,
                markerstrokewidth = 0,
                title = "AoC Day 14 - Iteration $cnt")
            frame(anim)
            resize!(xs, 0)
            resize!(ys, 0)
        end
    end
    if length(xs) > 0
        scatter!(xs, ys,
            markercolor = marker_color,
            markersize = marker_size,
            markeralpha = marker_alpha,
            markershape = marker_shape,
            markerstrokewidth = 0,
            title = "AoC Day 14 - Iteration $cnt")
    end
    frame(anim)
    gif(anim, "day14_anim.gif", fps = 13)
end


text:
using BinaryBuilder, Pkg

name = "MKL"
version = v"2022.0.0"

# Bash recipes for building across all platforms
script = read(joinpath(@__DIR__, "script.sh"), String)
script_macos = read(joinpath(@__DIR__, "script_macos.sh"), String)

non_reg_ARGS = filter(arg -> arg != "--register", ARGS)

platform_sources = [
    (
        platform = Platform("x86_64", "linux"; libc="glibc"),
        source = ArchiveSource(
            "https://anaconda.org/intel/mkl/2022.0.1/download/linux-64/mkl-2022.0.1-intel_117.tar.bz2",
            "62a3567efd4fa413181753dfc36c15bdfc8ec500eae4e5dcb104ed0d24136cfa";
            unpack_target = "mkl-x86_64-linux-gnu"
        ),
        # We need to run autofix on Linux, because here libmkl_rt doesn't
        # have a soname, so we can't ccall it without specifying the path:
        # https://github.com/JuliaSparse/Pardiso.jl/issues/69
        autofix = true,
        script = script,
    ),
    (
        platform = Platform("i686", "linux"; libc="glibc"),
        source = ArchiveSource(
            "https://anaconda.org/intel/mkl/2022.0.1/download/linux-32/mkl-2022.0.1-intel_117.tar.bz2",
            "fd800f09432a214dfffe9b973f24fdfd374aa013d6b92fbd570c744aa7eef5b2";
            unpack_target = "mkl-i686-linux-gnu"
        ),
        autofix = true,
        script = script,
    ),
    (
        platform = Platform("x86_64", "macos"),
        source = ArchiveSource(
            "https://anaconda.org/intel/mkl/2022.0.0/download/osx-64/mkl-2022.0.0-intel_105.tar.bz2",
            "fc026812458fd2053a6327bb5a42f47f97258147c39111bd6376dfabe51d05f5";
            unpack_target = "mkl-x86_64-apple-darwin14"
        ),
        # Need to disable autofix: updating linkage of libmkl_intel_thread.dylib on
        # macOS causes runtime issues:
        # https://github.com/JuliaPackaging/Yggdrasil/issues/915.
        autofix = false,
        script = script_macos,
    ),
    (
        platform = Platform("i686", "windows"),
        source = ArchiveSource(
            "https://anaconda.org/intel/mkl/2022.0.0/download/win-32/mkl-2022.0.0-intel_115.tar.bz2",
            "045c6f3ca31eca1e07785980152ad6a74513bead1f20beab440eded9596452a9";
            unpack_target = "mkl-i686-w64-mingw32"
        ),
        autofix = false,
        script = script,
    ),
    (
        platform = Platform("x86_64", "windows"),
        source = ArchiveSource(
            "https://anaconda.org/intel/mkl/2022.0.0/download/win-64/mkl-2022.0.0-intel_115.tar.bz2",
            "ba5a1f936b4e9fec0e3b5ae3ce44b56d06a6645a7c9456aa298c7d7373fd34c2";
            unpack_target = "mkl-x86_64-w64-mingw32"
        ),
        autofix = false,
        script = script,
    ),
]

# The products that we will ensure are always built
products = [
    LibraryProduct(["libmkl_core", "mkl_core", "mkl_core.2"], :libmkl_core),
    LibraryProduct(["libmkl_rt", "mkl_rt", "mkl_rt.2"], :libmkl_rt),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("IntelOpenMP_jll"),
]

non_reg_ARGS = filter(arg -> arg != "--register", ARGS)
include("../../fancy_toys.jl")
filter!(p -> should_build_platform(triplet(first(p))), platform_sources)

for (idx, (platform, source, autofix, script)) in enumerate(platform_sources)
    # Use "--register" only on the last invocation of build_tarballs
    if idx < length(platform_sources)
        args = non_reg_ARGS
    else
        args = ARGS
    end
    build_tarballs(args, name, version, [source], script, [platform], products, dependencies; lazy_artifacts = true, autofix = autofix)
end

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
import Base: @deprecate, @deprecate_binding, depwarn

### Changed in PR#629
# - HDF5.Dataset.xfer from ::hid_t to ::HDF5.Properties
@deprecate h5d_read(dataset_id::hid_t, memtype_id::hid_t, buf::AbstractArray, xfer::Properties) h5d_read(dataset_id, memtype_id, buf, xfer.id) false
@deprecate h5d_write(dataset_id::hid_t, memtype_id::hid_t, buf::AbstractArray, xfer::Properties) h5d_write(dataset_id, memtype_id, buf, xfer.id) false
@deprecate h5d_write(dataset_id::hid_t, memtype_id::hid_t, str::AbstractString, xfer::Properties) h5d_write(dataset_id, memtype_id, str, xfer.id) false
@deprecate h5d_write(dataset_id::hid_t, memtype_id::hid_t, x::T, xfer::Properties) where {T<:Union{ScalarType, Complex{<:ScalarType}}} h5d_write(dataset_id, memtype_id, x, xfer.id) false
@deprecate h5d_write(dataset_id::hid_t, memtype_id::hid_t, strs::Array{S}, xfer::Properties) where {S<:AbstractString} h5d_write(dataset_id, memtype_id, strs, xfer.id) false
@deprecate h5d_write(dataset_id::hid_t, memtype_id::hid_t, v::VLen{T}, xfer::Properties) where {T<:Union{ScalarType,CharType}} h5d_write(dataset_id, memtype_id, v, xfer.id) false
# - p_create lost toclose argument
@deprecate p_create(class, toclose::Bool, pv...) p_create(class, pv...)

### Changed in PR#632
# - using symbols instead of strings for property keys
@deprecate setindex!(p::Properties, val, name::AbstractString) setindex!(p, val, Symbol(name))

function Base.getindex(parent::Union{File,Group}, path::AbstractString, prop1::AbstractString, val1, pv...)
    depwarn("getindex(::Union{HDF5.File, HDF5.Group}, path, props...) with string key and value argument pairs is deprecated. Use keywords instead.", :getindex)
    props = (prop1, val1, pv...)
    return getindex(parent, path; [Symbol(props[i]) => props[i+1] for i in 1:2:length(props)]...)
end
function Base.setindex!(parent::Union{File,Group}, val, path::AbstractString, prop1::AbstractString, val1, pv...)
    depwarn("setindex!(::Union{HDF5.File, HDF5.Group}, val, path, props...) with string key and value argument pairs is deprecated. Use keywords instead.", :setindex!)
    props = (prop1, val1, pv...)
    return setindex!(parent, val, path; [Symbol(props[i]) => props[i+1] for i in 1:2:length(props)]...)
end

function h5open(filename::AbstractString, mode::AbstractString, pv...; kws...)
    depwarn("h5open with string key and value argument pairs is deprecated. Use keywords instead.", :h5open)
    return h5open(filename, mode; kws..., [Symbol(pv[i]) => pv[i+1] for i in 1:2:length(pv)]...)
end
function h5write(filename, name::AbstractString, data, pv...)
    depwarn("h5write with string key and value argument pairs is deprecated. Use keywords instead.", :h5write)
    return h5write(filename, name, data; [Symbol(pv[i]) => pv[i+1] for i in 1:2:length(pv)]...)
end
function h5read(filename, name::AbstractString, pv...)
    depwarn("h5read with string key and value argument pairs is deprecated. Use keywords instead.", :h5read)
    return h5read(filename, name; [Symbol(pv[i]) => pv[i+1] for i in 1:2:length(pv)]...)
end
function h5read(filename, name::AbstractString, indices::Tuple{Vararg{Union{AbstractRange{Int},Int,Colon}}}, pv...)
    depwarn("h5read with string key and value argument pairs is deprecated. Use keywords instead.", :h5read)
    return h5read(filename, name, indices; [Symbol(pv[i]) => pv[i+1] for i in 1:2:length(pv)]...)
end
function d_create(parent::Union{File,Group}, path::AbstractString, dtype::Datatype, dspace::D, prop1::AbstractString, val1, pv...) where D <: Union{Dataspace, Dims, Tuple{Dims,Dims}}
    depwarn("d_create with string key and value argument pairs is deprecated. Use keywords instead.", :d_create)
    props = (prop1, val1, pv...)
    return d_create(parent, path, dtype, dspace; [Symbol(props[i]) => props[i+1] for i in 1:2:length(props)]...)
end
function d_create(parent::Union{File,Group}, path::AbstractString, dtype::Type, dspace, prop1::AbstractString, val1, pv...)
    depwarn("d_create with string key and value argument pairs is deprecated. Use keywords instead.", :d_create)
    props = (prop1, val1, pv...)
    return d_create(parent, path, dtype, dspace; [Symbol(props[i]) => props[i+1] for i in 1:2:length(props)]...)
end
function p_create(class, prop1::AbstractString, val1, pv...)
    depwarn("p_create with string key and value argument pairs is deprecated. Use keywords instead.", :p_create)
    props = (prop1, val1, pv...)
    return p_create(class; [Symbol(props[i]) => props[i+1] for i in 1:2:length(props)]...)
end

### Changed in PR#652
# - read takes array element type, not Array with eltype
import Base: read
@deprecate read(obj::DatasetOrAttribute, ::Type{A}, I...) where {A<:Array} read(obj, eltype(A), I...)

### Changed in PR#657
# - using keywords instead of HDF5Properties objects

# deprecation helpers to avoid having default values on the equivalent low-level
# constructors in HDF5.jl
function __d_create(parent::Union{File, Group}, path::AbstractString, dtype::Datatype,
                   dspace::Dataspace, lcpl::Properties,
                   dcpl::Properties = DEFAULT_PROPERTIES,
                   dapl::Properties = DEFAULT_PROPERTIES,
                   dxpl::Properties = DEFAULT_PROPERTIES)
    d_create(parent, path, dtype, dspace, lcpl, dcpl, dapl, dxpl)
end
# a_create doesn't take property lists, so just bind the helper name directly
const __a_create = a_create

for (fsym, ptype) in ((:d_create, Union{File, Group}),
                      (:a_create, Union{File, Object}),
                     )
    chainsym = Symbol(:__, fsym)
    depsig = "$fsym(parent::$ptype, name::AbstractString, data, plists::HDF5Properties...)"
    usesig = "$fsym(parent::$ptype, name::AbstractString, data; properties...)"
    warnstr = "`$depsig` with property lists is deprecated, use `$usesig` with keywords instead"
    @eval begin
        function ($fsym)(parent::$ptype, name::AbstractString, data,
                         prop1::Properties, plists::Properties...)
            depwarn($warnstr, $(QuoteNode(fsym)))
            dtype = datatype(data)
            dspace = dataspace(data)
            obj = try
                ($chainsym)(parent, name, dtype, dspace, prop1, plists...)
            finally
                close(dspace)
            end
            return obj, dtype
        end
    end
end
for (fsym, ptype) in ((:d_write, Union{File,Group}),
                      (:a_write, Union{File,Object}),
                     )
    crsym = Symbol(:__, replace(string(fsym), "write" => "create"))
    depsig = "$fsym(parent::$ptype, name::AbstractString, data, plists::HDF5Properties...)"
    usesig = "$fsym(parent::$ptype, name::AbstractString, data; properties...)"
    warnstr = "`$depsig` with property lists is deprecated, use `$usesig` with keywords instead"
    @eval begin
        function ($fsym)(parent::$ptype, name::AbstractString, data,
                         prop1::Properties, plists::Properties...)
            depwarn($warnstr, $(QuoteNode(fsym)))
            dtype = datatype(data)
            obj = ($crsym)(parent, name, dtype, dataspace(data), prop1, plists...)
            try
                writearray(obj, dtype.id, data)
            catch exc
                o_delete(obj)
                rethrow(exc)
            finally
                close(obj)
                close(dtype)
            end
        end
    end
end
function Base.write(parent::Union{File,Group}, name::AbstractString, data::Union{T,AbstractArray{T}},
                    prop1::Properties, plists::Properties...) where {T<:Union{ScalarType,<:AbstractString,Complex{<:ScalarType}}}
    depwarn("`write(parent::Union{HDF5.File, HDF5.Group}, name::AbstractString, data, plists::HDF5Properties...)` " *
            "with property lists is deprecated, use " *
            "`write(parent::Union{HDF5.File, HDF5.Group}, name::AbstractString, data; properties...)` " *
            "with keywords instead.", :write)
    # We avoid using the d_write method to prevent double deprecation warnings.
    dtype = datatype(data)
    obj = __d_create(parent, name, dtype, dataspace(data), prop1, plists...)
    try
        writearray(obj, dtype.id, data)
    catch exc
        o_delete(obj)
        rethrow(exc)
    finally
        close(obj)
        close(dtype)
    end
end
function Base.write(parent::Dataset, name::AbstractString, data::Union{T,AbstractArray{T}},
                    prop1::Properties, plists::Properties...) where {T<:Union{ScalarType,<:AbstractString}}
    depwarn("`write(parent::HDF5Dataset, name::AbstractString, data, plists::HDF5Properties...)` " *
            "with property lists is deprecated, use " *
            "`write(parent::HDF5Dataset, name::AbstractString, data; properties...)` " *
            "with keywords instead.", :write)
    # We avoid using the a_write method to prevent double deprecation warnings.
    dtype = datatype(data)
    obj = __a_create(parent, name, dtype, dataspace(data), prop1, plists...)
    try
        writearray(obj, dtype.id, data)
    catch exc
        o_delete(obj)
        rethrow(exc)
    finally
        close(obj)
        close(dtype)
    end
end

### Changed in PR#664
# - normalized naming of C function wrappers
@deprecate_binding h5f_get_intend h5f_get_intent false
@deprecate_binding hf5start_swmr_write h5f_start_swmr_write false
@deprecate_binding h5d_oappend h5do_append false

### Changed in PR#678
# - normalized constants names to C definitions
@deprecate_binding Haddr haddr_t false
@deprecate_binding Herr herr_t false
@deprecate_binding Hid hid_t false
@deprecate_binding Hsize hsize_t false
@deprecate_binding Hssize hssize_t false
@deprecate_binding Htri htri_t false
@deprecate_binding Hvl_t hvl_t false

### Changed in PR#688
# - normalized more C type names
@deprecate_binding H5Ginfo H5G_info_t false
@deprecate_binding H5LInfo H5L_info_t false
@deprecate_binding H5Oinfo H5O_info_t false

### Changed in PR#689
# - switch from H5Rdereference1 to H5Rdereference2
@deprecate h5r_dereference(obj_id, ref_type, ref) h5r_dereference(obj_id, H5P_DEFAULT, ref_type, ref) false

### Changed in PR#690, PR#695
# - rename exported bindings
# - remove exports in > v0.14
@deprecate_binding HDF5Attribute HDF5.Attribute
@deprecate_binding HDF5File HDF5.File
@deprecate_binding HDF5Group HDF5.Group
@deprecate_binding HDF5Dataset HDF5.Dataset
@deprecate_binding HDF5Datatype HDF5.Datatype
@deprecate_binding HDF5Dataspace HDF5.Dataspace
@deprecate_binding HDF5Object HDF5.Object
@deprecate_binding HDF5Properties HDF5.Properties
@deprecate_binding HDF5Vlen HDF5.VLen
@deprecate_binding HDF5ChunkStorage HDF5.ChunkStorage
@deprecate_binding HDF5ReferenceObj HDF5.Reference false
@deprecate_binding HDF5Opaque HDF5.Opaque false

### Changed in PR#691
# - Make Reference() construct a null reference
@deprecate_binding HDF5ReferenceObj_NULL HDF5.Reference() false

### Changed in PR#693
@deprecate_binding ScalarOrString Union{ScalarType,String} false
@deprecate_binding HDF5Scalar ScalarType false
@deprecate_binding HDF5BitsKind BitsType false

### Changed in PR#695
import Base: names
@deprecate names(x::Union{Group,File,Attributes}) keys(x) false

### Changed in PR#694
@deprecate has(parent::Union{File,Group,Dataset}, path::AbstractString) Base.haskey(parent, path)
@deprecate exists(parent::Union{File,Group,Dataset,Datatype,Attributes}, path::AbstractString) Base.haskey(parent, path)


text:
using BenchmarkTools, Random
using Mill, Flux

using Plots, LaTeXStrings; pgfplotsx()

Random.seed!(42)
params = Base.Iterators.product([100, 200, 300], [100, 200, 300]) |> collect |> vec
for (m, n) in params
    @show m, n
    # A = rand(Float32, m, n)
    A = PreImputingMatrix(rand(Float32, m, n))
    # B = rand(Float32, n, m)
    B = [rand() < 0.5 ? missing : rand() for _ in 1:n, _ in 1:m]
    @btime $A * $B
end

xticklabels = [latexstring("$a|$b") for (a,b) in params]

t1 = [22.859, 50.551, 95.680, 38.428, 82.887, 200.255, 54.975, 133.416, 232.263] ./ 1000
l1 = "matrix(m,n) * Full Matrix(n,m)"
t2 = [3.868, 17.397, 39.650, 8.290, 35.292, 78.747, 13.223, 53.783, 121.800]
l2 = "Matrix(m,n) * Half Full Matrix(n,m)"
t3 = [22.114, 59.371, 96.551, 29.781, 67.402, 208.026, 36.066, 127.561, 227.967] ./ 1000
l3 = "Preimputingmatrix(m,n) * Full Matrix(n,m)"
t4 = [55.889, 119.639, 196.984, 94.915, 204.092, 335.075, 136.566, 298.371, 487.540] ./ 1000
l4 = "PreImputingmatrix(m,n) * Half Full Matrix(n,m)"
groupedbar(xticklabels, hcat(t1, t2, t3, t4); ylabel=L"ms", xlabel="n | m", labels=[l1 l2 l3 l4], title="PreImputing (log y)", yscale=:log10)

savefig("pre_imputing.png")

Random.seed!(42)
params = Base.Iterators.product([100, 200, 300], [100, 200, 300]) |> collect |> vec
for (m, n) in params
    @show m, n
    # A = rand(Float32, m, n)
    A = PreImputingMatrix(rand(Float32, m, n))
    # B = rand(Float32, n, m)
    B = [rand() < 0.5 ? missing : rand() for _ in 1:n, _ in 1:m]
    @btime gradient((A, B) -> sum(A*B), $A, $B)
end

xticklabels = [latexstring("$a|$b") for (a,b) in params]

t1 = [901.895 / 1000, 3.553, 7.889, 1.883, 7.232, 16.162, 2.648, 10.688, 24.344]
l1 = "Matrix(m,n) * Full Matrix(n,m)"
t2 = [909.361 / 1000, 3.650, 7.865, 1.792, 7.260, 16.200, 2.782, 10.598, 24.071]
l2 = "PreImputingMatrix(m,n) * Full Matrix(n,m)"
t3 = [1.077, 4.042, 9.111, 2.135, 7.637, 17.257, 3.075, 12.057, 24.922]
l3 = "PreImputingMatrix(m,n) * Half Full Matrix(n,m)"
groupedbar(xticklabels, hcat(t1, t2, t3); ylabel=L"ms", xlabel="n | m", labels=[l1 l2 l3], title="PreImputing (grad)")

savefig("pre_imputing_grad.png")


text:
using CalibrationErrors, Distances, Distributions, StatsBase
using CalibrationErrors: Bin, bin_eval, adddata!

using Random

Random.seed!(1234)

@testset "Simple example ($nclasses classes)" for nclasses in (2, 10, 100)
    # sample predictions and labels
    nsamples = 1_000
    predictions = rand(Dirichlet(nclasses, 1), nsamples)
    labels = rand(1:nclasses, nsamples)

    # create bin with all predictions and labels
    bin = Bin(predictions, labels)

    # check statistics
    @test bin.nsamples[] == nsamples
    @test bin.prediction_sum ≈ vec(sum(predictions, dims=2))
    @test bin.label_counts == counts(labels, nclasses)

    # check distance calculations
    for distance in (TotalVariation(), Cityblock(), Euclidean(), SqEuclidean())
        @test bin_eval(bin, distance) ≈
            nsamples * evaluate(distance, vec(mean(predictions; dims=2)),
                                proportions(labels, nclasses))
    end

    # compare with adding data to empty bin
    bin2 = Bin(eltype(predictions), nclasses)
    adddata!(bin2, predictions, labels)
    @test bin2.nsamples[] == bin.nsamples[]
    @test bin2.prediction_sum == bin.prediction_sum
    @test bin2.label_counts == bin.label_counts

    bin3 = Bin(eltype(predictions), nclasses)
    for i in axes(predictions, 2)
        adddata!(bin3, view(predictions, :, i), labels[i])
    end
    @test bin3.nsamples[] == bin.nsamples[]
    @test bin3.prediction_sum == bin.prediction_sum
    @test bin3.label_counts == bin.label_counts
end


text:
# %%
using Tmatrix

# inputs
rx, rz = 1e-6, 1.3e-6
n_max = 3
k1_r = 1e7; k1_i = 0.0
k2_r = 1.5e7; k2_i = 1e3

k1 = Complex(k1_r, k1_i)
k2 = Complex(k2_r, k2_i)

BigFloat_precision_list = 2 .^ [6,7,8,9]

println("***************************************************************")
println("***************************************************************")
println("Claculating T-matrix using Complex numbers")
println("***************************************************************")
println("***************************************************************")

println("using Float64 =====================================")
println("No symmetry exploited:")
@time T_no_symmetry_float_64 = calculate_Tmatrix_for_spheroid(rx, rz, n_max, k1, k2, rotationally_symmetric=false, symmetric_about_plane_perpendicular_z=false);

# %%
println("exploiting symmetry:")
@time T_with_symmetry_float_64 = calculate_Tmatrix_for_spheroid(rx, rz, n_max, k1, k2, rotationally_symmetric=true, symmetric_about_plane_perpendicular_z=true);

println("using BigFloat =====================================")
println("exploiting symmetry:")
T_with_symmetry_big_float = []
for precision_here = BigFloat_precision_list
    println("using precision $precision_here")
    @time T = calculate_Tmatrix_for_spheroid(rx, rz, n_max, k1, k2, rotationally_symmetric=true, symmetric_about_plane_perpendicular_z=true, BigFloat_precision = precision_here);
    append!(T_with_symmetry_big_float, T)
    println()
end

"""
n_array = collect(1:10)
elapsed_time = zeros(size(n_array))
for n = 1:length(n_array)
    elapsed_time[n] = @elapsed calculate_Tmatrix_for_spheroid(rx, rz, n_max, k1, k2, rotationally_symmetric=true, symmetric_about_plane_perpendicular_z=true, BigFloat_precision=8)
end
"""
println()
println()
println()
println("***************************************************************")
println("***************************************************************")
println("Claculating T-matrix using real numbers")
println("***************************************************************")
println("***************************************************************")

println("using Float64 =====================================")
println("No symmetry exploited:")
@time T_no_symmetry_float_64_SeparateRealImag = calculate_Tmatrix_for_spheroid_SeparateRealImag(rx, rz, n_max, k1_r, k1_i, k2_r, k2_i; rotationally_symmetric=false, symmetric_about_plane_perpendicular_z=false);

println("exploiting symmetry:")
@time T_with_symmetry_float_64_SeparateRealImag = calculate_Tmatrix_for_spheroid_SeparateRealImag(rx, rz, n_max, k1_r, k1_i, k2_r, k2_i; rotationally_symmetric=true, symmetric_about_plane_perpendicular_z=true);

println("using BigFloat =====================================")
println("exploiting symmetry:")
T_with_symmetry_big_float_SeparateRealImag = []
for precision_here = BigFloat_precision_list
    println("using precision $precision_here")
    @time T = calculate_Tmatrix_for_spheroid_SeparateRealImag(rx, rz, n_max, k1_r, k1_i, k2_r, k2_i; rotationally_symmetric=true, symmetric_about_plane_perpendicular_z=true, BigFloat_precision = precision_here);
    append!(T_with_symmetry_big_float_SeparateRealImag, T)
end


text:
using GHOST
time_start = now()
setup()
setup_parallel(5)
spdxs = execute(GHOST.PARALLELENABLER.conn,
                "SELECT spdx FROM gh_2007_2020.licenses ORDER BY spdx;",
                not_null = true) |>
    (obj -> getproperty.(obj, :spdx))
for spdx in spdxs
    queries(spdx)
end
time_end = now()
canonicalize(CompoundPeriod(time_end - time_start))


text:
using NQCModels.DiabaticModels: LargeDiabaticModel

struct LargeDiabaticCalculator{T,M} <: AbstractDiabaticCalculator{T,M}
    model::M
    potential::DependentField{Hermitian{T,Matrix{T}},Matrix{T}}
    derivative::DependentField{Matrix{Hermitian{T,Matrix{T}}},Matrix{T}}
    eigen::DependentField{LinearAlgebra.Eigen{T,T,Matrix{T},Vector{T}},Matrix{T}}
    adiabatic_derivative::DependentField{Matrix{Matrix{T}},Matrix{T}}
    nonadiabatic_coupling::DependentField{Matrix{Matrix{T}},Matrix{T}}
    tmp_mat::Matrix{T}
    stats::Dict{Symbol,Int}
    function LargeDiabaticCalculator{T}(model::M, atoms::Integer) where {T,M<:Model}
        mat = NQCModels.DiabaticModels.matrix_template(model, T)
        vec = NQCModels.DiabaticModels.vector_template(model, T)

        potential = Hermitian(zero(mat))
        derivative = [Hermitian(zero(mat)) for i=1:ndofs(model), j=1:atoms]
        eigen = Eigen(zero(vec), zero(mat)+I)
        adiabatic_derivative = [zero(mat) for i=1:ndofs(model), j=1:atoms]
        nonadiabatic_coupling = [zero(mat) for i=1:ndofs(model), j=1:atoms]
        tmp_mat = zero(mat)

        position = fill(NaN, ndofs(model), atoms)

        stats = Dict{Symbol,Int}(
            :potential=>0,
            :derivative=>0,
            :eigen=>0,
            :adiabatic_derivative=>0,
            :nonadiabatic_coupling=>0
        )

        new{T,M}(
            model,
            DependentField(potential, copy(position)),
            DependentField(derivative, copy(position)),
            DependentField(eigen, copy(position)),
            DependentField(adiabatic_derivative, copy(position)),
            DependentField(nonadiabatic_coupling, copy(position)),
            tmp_mat,
            stats,
        )
    end
end

function Calculator(model::LargeDiabaticModel, atoms::Integer, T::Type=Float64)
    LargeDiabaticCalculator{T}(model, atoms)
end

function evaluate_potential!(calc::Union{LargeDiabaticCalculator,DiabaticFrictionCalculator}, r)
    calc.stats[:potential] += 1
    NQCModels.potential!(calc.model, calc.potential, r)
    return nothing
end

function evaluate_eigen!(calc::Union{LargeDiabaticCalculator,DiabaticFrictionCalculator}, r)
    calc.stats[:eigen] += 1
    potential = get_potential(calc, r)
    eig = LinearAlgebra.eigen(potential)
    correct_phase!(eig, calc.eigen.vectors)
    copyto!(calc.eigen.vectors, eig.vectors)
    copyto!(calc.eigen.values, eig.values)
end

function correct_phase!(eig::LinearAlgebra.Eigen, old_eigenvectors::AbstractMatrix)
    @views for i=1:length(eig.values)
        if LinearAlgebra.dot(eig.vectors[:,i], old_eigenvectors[:,i]) < 0
            eig.vectors[:,i] .*= -1
        end
    end
end

function evaluate_adiabatic_derivative!(calc::Union{LargeDiabaticCalculator,DiabaticFrictionCalculator}, r)
    calc.stats[:adiabatic_derivative] += 1
    eigen = get_eigen(calc, r)
    derivative = get_derivative(calc, r)
    @inbounds for i in NQCModels.mobileatoms(calc)
        for j in NQCModels.dofs(calc)
            LinearAlgebra.mul!(calc.tmp_mat, derivative[j,i], eigen.vectors)
            LinearAlgebra.mul!(calc.adiabatic_derivative[j,i], eigen.vectors', calc.tmp_mat)
        end
    end
end

function evaluate_nonadiabatic_coupling!(calc::Union{LargeDiabaticCalculator,DiabaticFrictionCalculator}, r)
    calc.stats[:nonadiabatic_coupling] += 1
    eigen = get_eigen(calc, r)
    adiabatic_derivative = get_adiabatic_derivative(calc, r)

    evaluate_inverse_difference_matrix!(calc.tmp_mat, eigen.values)

    @inbounds for i in NQCModels.mobileatoms(calc)
        for j in NQCModels.dofs(calc)
            multiply_elementwise!(calc.nonadiabatic_coupling[j,i], adiabatic_derivative[j,i], calc.tmp_mat)
        end
    end
end

function evaluate_inverse_difference_matrix!(out, eigenvalues)
    @inbounds for i in eachindex(eigenvalues)
        for j in eachindex(eigenvalues)
            out[j,i] = 1 / (eigenvalues[i] - eigenvalues[j])
        end
        out[i,i] = zero(eltype(out))
    end
end

function multiply_elementwise!(coupling::Matrix, adiabatic_derivative::Matrix, eigenvalue_difference_matrix::Matrix)
    @inbounds for I in eachindex(coupling, adiabatic_derivative, eigenvalue_difference_matrix)
        coupling[I] = adiabatic_derivative[I] * eigenvalue_difference_matrix[I]
    end
end


text:
include("plain.jl")
include("JLD.jl")


text:
@testset "539.minimum-time-difference.jl" begin
    @test find_min_difference(["23:59", "00:00"]) == 1
    @test find_min_difference(["00:00", "23:59", "00:00"]) == 0
end


text:
"""
Compute the SCV of data.
"""
function scv(data::Vector{Float64})
    return (std(data)/mean(data))^2
end

"""
Merge PH distributions into a MAPH distribution
"""
function merge_dist(probs::Vector{Float64},dist::Vector{Any})
    merged_α = probs.*[reshape(dist[i].α, model_size(dist[i])) for i = 1:length(dist)]
    # merged_α = (probs ./ sum(probs)) .* [reshape(dist[i].α, model_size(dist[i])) for i = 1:length(dist)]

    α_temp = reduce(vcat,merged_α)'
    merged_T = [dist[i].T for i = 1:length(dist)]
    # merged_T = (probs ./ sum(probs)) .* [dist[i].T for i = 1:length(dist)]
    T_temp = cat(merged_T..., dims = (1,2))
    @show dist
    # merged_T0 = -(probs./sum(probs)) .* [sum(dist[i].T,dims=2) for i = 1:length(dist)]
    merged_T0 = -1.0.*[sum(dist[i].T,dims=2) for i = 1:length(dist)]
    @show merged_T0
    T0_temp = cat(merged_T0...,dims=(1,2))
    display(T0_temp)
    display(T_temp)

    return α_temp, T_temp, T0_temp
end

function vertical_merge_matrix(M::Matrix{Float64},p::Int64)
    m, n = size(M)

    new_M = zeros((p,n))

    for i=1:(p-1)
        new_M[i,:] = M[i,:]
    end

    new_M[p,:] = sum(M[p:m,:],dims=1)

    return new_M
end

function horizontal_merge_matrix(M::Matrix{Float64},p::Int64)
    m, n = size(M)

    new_M = zeros((m,p))

    for i=1:(p-1)
        new_M[:,i] = M[:,i]
    end

    new_M[:,p] = sum(M[:,p:m],dims=2)

    return new_M
end


text:
# KnetPtr type holds a gpu allocated pointer.  We try to minimize the number of actual
# allocations, which are slow, by reusing preallocated but garbage collected pointers.

mutable struct KnetPtr
    ptr                         # actual pointer, removed type ::Cptr for serialization
    len::Int                    # size in bytes
    dev::Int                    # id of the device the pointer belongs to
    parent::KnetPtr             # used to implement shared memory pointers

    # This is the low level KnetPtr constructor, it adds the finalizer and
    # does not assign parent which is only needed for shared pointers.

    function KnetPtr(ptr::Cptr,len::Int,dev::Int)
        kp = new(ptr,len,dev)
        finalizer(freeKnetPtr, kp)
    end

    # This constructor is used to create a shared pointer.  We need to
    # keep the parent field to prevent premature gc of the parent.  The
    # child does not need a special finalizer.

    function KnetPtr(parent::KnetPtr, offs::Int, len::Int)
        if len < 0 || offs < 1 || offs+len-1 > parent.len; throw(BoundsError()); end
        new(parent.ptr+offs-1, len, parent.dev, parent)
    end

    # This one is used by serialize:
    KnetPtr(ptr::Array{UInt8},len::Int)=new(ptr,len,-1)

end

# When Julia gc reclaims a KnetPtr object, the following special finalizer does not actually
# release the memory, but inserts it back in the appropriate pool for reuse.

function freeKnetPtr(p::KnetPtr)
    if p.ptr == C_NULL || isdefined(p,:parent); return; end
    mem = KnetMems[p.dev+1]
    mem.avail += p.len
    push!(mem.pools[p.len].free, p.ptr)
    p.ptr = C_NULL # to avoid double free by gcnode then gc.
end

# We use the KnetPool type to keep track of allocated and garbage collected pointers: We
# keep one KnetPool struct per size per device.

mutable struct KnetPool
    free::Vector{Cptr}          # pointers available for reuse
    nptr::Int                   # number of allocated pointers
    KnetPool()=new(Cptr[],0)
end

const KNETMEMINIT = 1<<24       # initial gpu memory limit

# KnetMem type keeps memory information for one device.
mutable struct KnetMem
    pools::Dict{Int,KnetPool}   # pointers of a given size
    limit::Int                  # current memory limit
    bytes::Int                  # total bytes allocated (inuse + avail)
    avail::Int                  # total bytes freed and available
    KnetMem()=new(Dict{Int,KnetPool}(),KNETMEMINIT,0,0)
end

# KnetMems[dev+1] holds memory information for device dev.
KnetMems = nothing
initKnetMems() = (global KnetMems = [ KnetMem() for i in 1:gpuCount() ])

# Blocksize determines the actual allocation size given the array size in bytes, and can be
# larger than what the array needs for increased reuse.
# blocksize(n::Int)=n
# function blocksize(n::Int)
#     b = sizeof(n)<<3
#     z = leading_zeros(n-1)
#     1<<(b-z)
# end
blocksize(n::Int,b=cbrt(2))=floor(Int,b^ceil(log(b,n)))

# The following used for debugging and record every request
arraysizes = Int[]; allocs = Int[]; blocksizes = Int[]

function inclimit!(m::KnetMem, minlimit=m.limit*6÷5)
    maxlimit = gpumem()[1] - 500_000_000 # m.bytes + gpufree()
    minlimit = min(maxlimit, minlimit)
    if m.limit < minlimit
        m.limit = minlimit
        @debug "limit=$(m.limit) bytes=$(m.bytes) avail=$(m.avail) pools=$(length(m.pools)) blocks=$(sum(p->p.nptr,values(m.pools)))"
    end
end

# This the main KnetPtr constructor.  It tries to avoid actual
# allocation which is slow.  Reusing a pointer is very fast.
# Allocating a new one is about 10 μs.  gc() is about 75 ms.

function KnetPtr(arraybytes::Int)
    dev = gpu(); if dev < 0; error("KnetPtr: bad device id $dev."); end
    if KnetMems==nothing; initKnetMems(); end
    mem = KnetMems[dev+1]
    blockbytes = blocksize(arraybytes)
    @dbg (push!(arraysizes,arraybytes); push!(blocksizes,blockbytes))
    pool = get!(KnetPool,mem.pools,blockbytes)
    if !isempty(pool.free)      # 1. best case we have one available in pool
        @dbg push!(allocs, 1)
        mem.avail -= blockbytes
        return KnetPtr(pop!(pool.free),blockbytes,dev)
    end
    if mem.bytes + blockbytes <= mem.limit # 2. allocate if we are within limit
        ptr = knetMalloc(blockbytes)
        if ptr != nothing
            @dbg push!(allocs, 2)
            mem.bytes += blockbytes
            pool.nptr += 1
            return KnetPtr(ptr,blockbytes,dev)
        end
    end
    if pool.nptr > 0            # 3. try gc if we had allocated this size before
        GC.gc();  @dbg print('-')
        # inclimit!(mem, (mem.bytes-mem.avail)*3) # *4=>18s *3=>16s *5÷2=>16s *2=>22s *6÷5=>30s
        inclimit!(mem, mem.limit*6÷5)
        if !isempty(pool.free)
            @dbg push!(allocs, 3)
            mem.avail -= blockbytes
            return KnetPtr(pop!(pool.free),blockbytes,dev)
        end
    end
    # At this point we have to either free old ptrs or increase limit
    # Let's just try limit increase for now
    inclimit!(mem, (mem.bytes + blockbytes*2))
    ptr = knetMalloc(blockbytes)
    if ptr == nothing
        Knet.gc(); @dbg print('+')  # last ditch effort
        ptr = knetMalloc(blockbytes)
    end
    if ptr != nothing
        @dbg push!(allocs, 4)
        mem.bytes += blockbytes
        pool.nptr += 1
        return KnetPtr(ptr,blockbytes,dev)
    end
    error("Out of gpu memory")
end

KnetPtr(n::Integer)=KnetPtr(Int(n))

# This does the actual allocation, returns `nothing` in case of error
function knetMalloc(nbytes::Int)
    ptr = Cptr[0]
    ret = @cudart1(cudaMalloc,(Ptr{Cptr},Csize_t),ptr,nbytes)
    ret == 0 ? ptr[1] : nothing
end

"""
    Knet.gc(dev=gpu())

cudaFree all pointers allocated on device `dev` that were previously allocated and garbage
collected. Normally Knet holds on to all garbage collected pointers for reuse. Try this if
you run out of GPU memory.
"""
function gc(dev=gpu())
    if KnetMems == nothing; return; end
    mem = KnetMems[dev+1]
    GC.gc(); GC.enable(false)
    for (n,v) in mem.pools
        for p in v.free
            @cudart(cudaFree,(Cptr,),p)
        end
        v.nptr -= length(v.free)
        mem.avail -= n * length(v.free)
        mem.bytes -= n * length(v.free)
        empty!(v.free)
    end
    GC.enable(true); GC.gc()
end

function knetgc()
    @warn "knetgc is deprecated, use Knet.gc instead" maxlog=1
    Knet.gc()
end

# Some utilities
meminfo(i=gpu())=(KnetMems==nothing ? [] : [(k,v.nptr,length(v.free)) for (k,v) in KnetMems[i+1].pools])

function gpuinfo(msg="",dev=gpu();n=10)
    msg != "" && print("$msg ")
    if nvmlfound # Libdl.find_library(["libnvidia-ml"],[]) != ""
        g = nvmlDeviceGetMemoryInfo(dev)
        println((:dev,dev,:total,g[1],:free,g[2],:used,g[3]))
    else
        dev0 = gpu(); gpu(dev)
        (free,total) = cudaMemGetInfo()
        gpu(dev0)
        println((:dev,dev,:total,total,:free,free,:used,total-free))
    end
    bytes = avail = ptrs = k = 0
    for (s,u,f) in sort(meminfo(dev), by=(x->x[1]*x[2]), rev=true)
        bytes += s*u; avail += s*f; ptrs += u; k += 1
        if k <= n; println((:bytes,s*u,:size,s,:alloc,u,:avail,f)); end
    end
    if n < k; println('⋮'); end
    println((:totbytes,bytes,:avail,avail,:ptrs,ptrs))
    if KnetMems != nothing
        mem = KnetMems[dev+1]
        println((:membytes,mem.bytes,:avail,mem.avail,:limit,mem.limit,:pools,length(mem.pools)))
    end
end

function memdbg(msg="")
    if nvmlfound # Libdl.find_library(["libnvidia-ml"],[]) != ""
        m = nvmlDeviceGetMemoryInfo()
    else
        m = (0,0,0)
    end
    c = cudaMemGetInfo()
    x = [0,0]
    for (s,u,f) in meminfo()
        x[1] += s*u
        x[2] += s*f
    end
    println("""$msg:
cudaMemGetInfo: ctotal: $(c[2]) cfree: $(c[1]) ctotal-cfree: $(c[2]-c[1])
nvmlDeviceGetMemoryInfo: ntotal: $(m[1]) nfree: $(m[2]) nused: $(m[3]) nfree+used: $(m[2]+m[3])
KnetPtr: ktotal: $(x[1]) kavail: $(x[2]) ktotal-kavail: $(x[1]-x[2])
nused-ktotal: $(m[3]-x[1])
""")
end


text:
using Pkg
Pkg.instantiate()

using PackageCompiler

create_sysimage(:MusicAlbums;
    sysimage_path="MusicAlbums.so",
    precompile_execution_file="deploy/precompile.jl")

max_stars_repo_name: vanhp/MicroServiceJulia

max_issues_repo_name: vanhp/MicroServiceJulia

max_forks_repo_name: vanhp/MicroServiceJulia


text:
# Iterated extended Kalman smoothing
mutable struct IEKS{IT} <: AbstractEK
    prior::Symbol
    order::Int
    diffusionmodel::Symbol
    smooth::Bool
    initialization::IT
    linearize_at::Any
end

"""
    IEKS(; prior=:ibm, order=1, diffusionmodel=:dynamic, initialization=TaylorModeInit(), linearize_at=nothing)

**Gaussian ODE filtering with iterated extended Kalman smoothing.**

To use it, use
`solve_ieks(prob, IEKS(), args...)`
instead of
`solve(prob, IEKS(), args...)`,
since it is implemented as an outer loop around the solver.

Currently, only the integrated Brownian motion prior `:ibm` is supported.
For the diffusionmodel, chose one of
`[:dynamic, :dynamicMV, :fixed, :fixedMV, :fixedMAP]`.
Just like the [`EK1`](@ref) it requires that the Jacobian of the rhs function is available.

See also: [`EK0`](@ref), [`EK1`](@ref), [`solve_ieks`](@ref)

# References:
- F. Tronarp, S. Särkkä, and P. Hennig: **Bayesian ODE Solvers: The Maximum A Posteriori Estimate**
"""
function IEKS(;
    prior=:ibm,
    order=1,
    diffusionmodel=:dynamic,
    initialization=TaylorModeInit(),
    linearize_at=nothing,
)
    if !isnothing(linearize_at)
        @assert linearize_at isa ProbODESolution
        @assert linearize_at.alg.prior == prior
        @assert linearize_at.alg.order == order
        @assert linearize_at.alg.diffusionmodel == diffusionmodel
        @assert linearize_at.alg.initialization == initialization
        @assert linearize_at.alg.smooth == true
    end
    return IEKS(prior, order, diffusionmodel, true, initialization, linearize_at)
end

"""
    solve_ieks(prob::AbstractODEProblem, alg::IEKS, args...; iterations=10, kwargs...)

Solve method to be used with the [`IEKS`](@ref). The IEKS works essentially by solving the
ODE multiple times. `solve_ieks` therefore wraps a call to the standard `solve` method,
passing `args...` and `kwargs...`.

Currently, this method is very simplistic - it iterates for a fixed numer of times and does
not use a stopping criterion.
"""
function solve_ieks(
    prob::DiffEqBase.AbstractODEProblem,
    alg::IEKS,
    args...;
    iterations=10,
    kwargs...,
)
    sol = nothing
    for i in 1:iterations
        alg.linearize_at = sol
        sol = solve(prob, alg, args...; kwargs...)
    end
    return sol
end


text:
#=~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~=#
# Problem set 6 solutions
# Written by Tyler Ransom
# Commented by Giuseppe Grasso
#=~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~=#

using Random
using LinearAlgebra
using Statistics
using Optim
using DataFrames
using DataFramesMeta
using CSV
using HTTP
using GLM
cd("/Users/peppegrass/Documents/GitHub/fall-2020/ProblemSets/PS6-ccp/") # GG: sets working directory
pwd() ## GG: prints working directory
readdir() # GG: equivalent to -ls- to see elements of working directory

# read in function to create state transitions for dynamic model
include("create_grids.jl")

function allwrap()

    #=
    GG: In this problem set, we will repeat the estimation of the simpliﬁed version of the Rust (1987, Econometrica) bus engine replacement
    model. Rather than solve the model by backwards recursion, we will exploit the renewal property of the replacement decision and estimate
    the model using conditional choice probabilities (CCPs).
    =#

# GG: I - PRELIMINARIES: Loading data, reshaping them, estimating flexible logit model

    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # Question 1: Follow the directions from PS5 to read in the data (the second CSV ﬁle you read in as part of PS5)
    # and reshape to “long” panel format, calling your long dataset df_long.
    #:::::::::::::::::::::::::::::::::::::::::::::::::::

    url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS5-ddc/busdata.csv"
    df = CSV.read(HTTP.get(url).body)

    # create bus id variable
    df = @transform(df, bus_id = 1:size(df,1))

    #---------------------------------------------------
    # reshape from wide to long (must do this thrice be-
    # cause DataFrames.stack() requires doing it one
    # variable at a time)
    #---------------------------------------------------
    # first reshape the decision variable
    dfy = @select(df, :bus_id,:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20,:RouteUsage,:Branded,:Zst)
    dfy_long = DataFrames.stack(dfy, Not([:bus_id,:RouteUsage,:Branded,:Zst]))
    rename!(dfy_long, :value => :Y)
    dfy_long = @transform(dfy_long, time = kron(collect([1:20]...),ones(size(df,1))))
    select!(dfy_long, Not(:variable))

    # next reshape the odometer variable
    dfx = @select(df, :bus_id,:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20)
    dfx_long = DataFrames.stack(dfx, Not([:bus_id]))
    rename!(dfx_long, :value => :Odometer)
    dfx_long = @transform(dfx_long, time = kron(collect([1:20]...),ones(size(df,1))))
    select!(dfx_long, Not(:variable))

    # next reshape the odometer variable
    dfxs = @select(df, :bus_id,:Xst1,:Xst2,:Xst3,:Xst4,:Xst5,:Xst6,:Xst7,:Xst8,:Xst9,:Xst10,:Xst11,:Xst12,:Xst13,:Xst14,:Xst15,:Xst16,:Xst17,:Xst18,:Xst19,:Xst20)
    dfxs_long = DataFrames.stack(dfxs, Not([:bus_id]))
    rename!(dfxs_long, :value => :Xst)
    dfxs_long = @transform(dfxs_long, time = kron(collect([1:20]...),ones(size(df,1))))
    select!(dfxs_long, Not(:variable))

    # join reshaped df's back together
    df_long = leftjoin(dfy_long, dfx_long, on = [:bus_id,:time])
    df_long = leftjoin(df_long, dfxs_long, on = [:bus_id,:time])
    sort!(df_long,[:bus_id,:time])

    # variables in array form
    Y      = Matrix(df[:,[:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20]])
    X      = Matrix(df[:,[:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20]])
    B      = Vector(df[:,:Branded])
    Xstate = Matrix(df[:,[:Xst1,:Xst2,:Xst3,:Xst4,:Xst5,:Xst6,:Xst7,:Xst8,:Xst9,:Xst10,:Xst11,:Xst12,:Xst13,:Xst14,:Xst15,:Xst16,:Xst17,:Xst18,:Xst19,:Xst20]])
    Zstate = Vector(df[:,:Zst])
    N      = size(Zstate,1)
    T      = size(Xstate,2)

    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # Question 2: Estimate a ﬂexible logit model where the dependent variable is the replacement decision and the right
    # hand side is a fully interacted set of the following variables:
    # • Mileage
    # • Mileage^2
    # • Route Usage
    # • Route Usage^2
    # • Branded
    # • Time period
    # • Time period^2
    # Hint: “Fully interacted” means that all terms from 1st order to 7th order (e.g. Odometer^2 × RouteUsage^2 × Branded × time^2).
    # Hint: Julia’s GLM package allows you to easily accomplish this by specifying the interacted variables with asterisks in between them,
    # e.g. Odometer * RouteUsage estimates a model that includes Odometer, Route Usage and the product of the two.
    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    θ̂_flex = glm(@formula(Y ~ Odometer * Odometer * RouteUsage *RouteUsage * Branded * time * time), df_long, Binomial(), LogitLink())
    println(θ̂_flex)
    # GG: Using asterisks in GLM package in Julia it gives you the fully interacted model
    # GG: we do this to get probabilities non-parametrically

#= GG: II - DYNAMIC ESTIMATION with CCPs
We will use the ﬂexible logit parameters to generate CCPs which we can use to compute the future value term as alternative to the
backwards recursion we did in PS5.
Recall the model from PS5, where the differenced conditional value function for running the bus (relative to replacing it) was:

(1): $ v_{1t}(x_{t},b) - v_{0t}(x_{t},b) =
        θ_{0} + θ_{1}x_{1t} + θ_{2}b                    [GG: same as static so far, flow utility]
        + β ∫ V_{t+1}(x_{t+1}, b) dF(x_{t+1} | x_{t})     [GG: recursive component]
        $

and where V_{t+1} is the value function and the integral is over transitions in the mileage states x_{t}.
By exploiting the RENEWAL PROPERTY of the decision property, we can express V_{t+1} instead as v_{0t+1}−log(p_{0t+1}).
And since v_{0t+1} corresponds to the renewal action, we know that it is equivalent to [GG: something missing in assignment text; see lecture notes].
Thus, our value function formulation can be simpliﬁed to:

(2): $ v_{1t}(x_{t},b) - v_{0t}(x_{t},b) =
        θ_{0} + θ_{1}x_{1t} + θ_{2}b
        + β ∫ log(p_{0t+1})(x_{t+1}, b) dF(x_{t+1} | x_{t})     [GG: log(p_{0t+1}) instead of V_{t+1}]
        $

and by discretizing the integral, we can simplify this even further to be:

(3): $ v_{1t}(x_{t},b) - v_{0t}(x_{t},b) =
        θ_{0} + θ_{1}x_{1t} + θ_{2}b
        + β ∑_{x_{1,t+1}} log(p_{0t+1})(x_{t+1}, b) [ f_1( x_{1,t+1} | x_{1,t},x_2 ) - f_0( x_{1,t+1} | x_{1,t},x_2 ) ]
        dF(x_{t+1} | x_{t})     [GG: log(p_{0t+1}) instead of V_{t+1}]
        $

where the f_j’s are deﬁned identically as in PS5 .
=#
    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # Question 3: Estimate the θ’s using (3) and assuming a discount factor of β = 0.9.
    # I will walk you through speciﬁc steps for how to do this:
    #
    # (a) Construct the state transition matrices using the exact same code as in this step of PS5.
    zval,zbin,xval,xbin,xtran = create_grids()
    df_state = DataFrame(Odometer = kron(ones(zbin),xval), RouteUsage = kron(zval,ones(xbin)), time = zeros(size(xtran,1)), Branded = zeros(size(xtran,1)))

    # (b) Compute the future value terms for all possible states of the model. Basically, what we want is −log p 0t+1 evaluated at every possible state of the model (t,b,x 1,t ,x 2 ). The easiest way to do this is to adjust the data that we feed into a predict()
    # function using the ﬂexible logit coefﬁcients from question number 2.
    #   • First, create a data frame that has four variables:
    #       – Odometer reading (equals kron(ones(zbin),xval))
    #       – Route usage (equals kron(ones(xbin),zval))
    #       – Branded (equals 0s of the same size as Odometer and Route usage)
    #       – time (equals 0s of the same size as Branded)
    #   • Now write a function that reads in this data frame, the ﬂexible logit estimates, and the other state variables
    #     (Xstate, Zstate, xtran, etc.)
    #   • Initialize the future value array, which should be a 3-dimensional array of zeros.
    #     The size of the ﬁrst dimension should be the total number of grid points (i.e. the number of rows of xtran).
    #     The second dimension should be 2, which is the possible outcomes of :Branded. The third dimension should be T+1.
    #     Note that the number of rows of the future value array should equal the number of rows of the state data frame.
    #   • Now write two nested for loops:
    #       – Loop over t from 2 to T
    #       – Loop over the two possible brand states { 0,1 }
    #   • Inside all of the for loops, make the following calculations
    #       – Update your state data frame so that the :time variable takes on the value of t and the :Branded variable
    #         takes on the value of b
    #       – Compute p0 using the predict() function applied to your updated data frame and the ﬂexible logit estimates
    #       – Store in the FV array the value of −β log p_0 . Remember that every row of the data frame corresponds to the rows in the
    #         state transition matrix, so you can vectorize this calculation.
    #   • Now multiply the state transitions by the future value term. This requires writing another for loop that goes over the rows in
    #     the original data frame (the one that you read in at the very beginning of this problem set). In other words, loop over i and t.
    #     To get the actual rows of the state transition matrix (since we don’t need to use all possible rows), you should re-use the similar
    #     code from PS5; something like this: FVT1[i,t] = (xtran[row1,:].-xtran[row0,:])'* FV[row0:row0+xbin-1,B[i]+1,t+1]
    #     The purpose of this loop is to map the CCPs from the each-possible-state-is-a-row data frame to the actual data frame we used
    #     to estimate the ﬂexible logit in question 2.
    #   • Your function should return FVT1 in “long panel” format. I used FVT1’[:] to make this conversion, but you should double check that
    #     your i and t indexing of your original data frame matches.

    # GG: Here we're computing FV using CCPs, there's no recursion
    function fvdata(b1,β,df,xbin,zbin,Zstate,Xstate,xtran,N,T,B)
        FV1  = zeros(xbin*zbin,2,T+1)
        FVT1 = zeros(N,T)
        for t=2:T
            for s=0:1
                @with(df, :time    .= t)
                @with(df, :Branded .= s)
                p0 = 1 .- convert(Array{Float64},predict(b1,df)) # GG: the b1 inside predict() is where we'll feed the θ̂_flex we estimated
                FV1[:,s+1,t] = -β.*log.(p0)
            end
        end

        for i=1:N # GG: this is exactly as in PS5, where we multiply the state transitions by the future value term.
            row0 = (Zstate[i]-1)*xbin+1
            for t=1:T
                row1      = row0+Xstate[i,t]-1
                FVT1[i,t] = (xtran[row1,:].-xtran[row0,:])⋅FV1[row0:row0+xbin-1,B[i]+1,t+1]
            end
        end
        return FVT1'[:]
    end
    # GG: calling the function, passing the arguments to it
    @time fvt1 = fvdata(θ̂_flex,0.9,df_state,xbin,zbin,Zstate,Xstate,xtran,N,T,df.Branded)
    df_long = @transform(df_long, fv = fvt1)

    # (c) Estimate the structural parameters.
    #   • Add the output of your future value function as a new column in the original “long panel” data frame.
    #     The easiest way to do this is df long = @transform(df long, fv = fvt1)
    #   • Now use the GLM package to estimate the structural model. Make use of the “offset” function to add the future value term as another
    #     regressor whose coefﬁcient is restricted to be 1. That is: theta_hat_ccp_glm = glm(@formula(Y ~ Odometer + Branded), df_long, Binomial(), LogitLink(), offset=df_long.fv)
    θ̂_ccp_glm = glm(@formula(Y ~ Odometer + Branded), df_long, Binomial(), LogitLink(), offset=df_long.fv)
    println(θ̂_ccp_glm)

    # GG: if you're dealing with a multinomial model, you might need to not use GLM and code it up yourself

    # (d) Optionally, you can write your own function to estimate a binary logit where you restrict the offset term to have a coefﬁcient of 1.
    # (I will include this code in my solutions.)
    # # You can also do this via Optim, and you'll get the exact same answer
    # FV = reshape(df_long.fv,T,N)'
    # @views @inbounds function likebus_ccp(α,Y,X,B,FV,N,T)
        # like = 0
        # for i=1:N
            # for t=1:T
                # v1    = α[1] + α[2]*X[i,t] + α[3]*B[i] + FV[i,t]
                # dem   = 1 + exp(v1)
                # like -= ( (Y[i,t]==1)*v1 ) - log(dem) # negative of the binary logit likelihood
            # end
        # end
        # return like
    # end
    # θ_true = [2;-.15;1]
    # # evaluate objective function
    # @time likebus_ccp(θ_true,Y,X,B,FV,N,T)
    # @time likebus_ccp(θ_true,Y,X,B,FV,N,T)
    # # estimate likelihood function
    # θ̂_optim = optimize(a -> likebus_ccp(a,Y,X,B,FV,N,T), θ_true, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true))
    # θ̂_ccp = θ̂_optim.minimizer
    # println(θ̂_ccp)

    # (e) Wrap all of your code in an empty function as you’ve done with other problem sets. Prepend your wrapper function call (at the very
    # end of the script) with @time so that you can time how long everything takes. (On my machine, everything took under 20 seconds.)
    return nothing
end
@time allwrap()
    # (f) Glory in the power of CCPs!
    #:::::::::::::::::::::::::::::::::::::::::::::::::::


text:
using Test

include("enchantments.jl")

@testset "Retrieve a card from a stack" begin
    @test card([1, 3, 4, 1], 2) == 3
    @test_throws BoundsError card([1, 3, 4, 1], 0)
end

@testset "Exchange a card in the stack" begin
    @test replace_card!([1, 3, 4, 1], 2 => 6) == [1, 6, 4, 1]
    @test_throws BoundsError replace_card!([1, 3, 4, 1], 5 => 6)
end

@testset "Insert a card at the of top the stack" begin
    @test insert_card_at_top!([1, 3, 4, 1], 8) == [1, 3, 4, 1, 8]
end

@testset "Remove a card from the stack" begin
    @test remove_card!([1, 3, 4, 1], 2) == [1, 4, 1]
    @test_throws BoundsError remove_card!([], 0)
    @test_throws BoundsError remove_card!([], 1)
    @test_throws BoundsError remove_card!([1, 3, 4, 1], 0)
    @test_throws BoundsError remove_card!([1, 3, 4, 1], 5)
end

@testset "Remove the top card from the stack" begin
    @test remove_card_from_top!([1, 3, 4, 1]) == [1, 3, 4]
    @test_throws ArgumentError remove_card_from_top!([])
end

@testset "Insert a card at the bottom of the stack" begin
    @test insert_card_at_bottom!([1, 3, 4, 1], 8) == [8, 1, 3, 4, 1]
end

@testset "Remove a card from the bottom of the stack" begin
    @test remove_card_from_bottom!([1, 3, 4, 1]) == [3, 4, 1]
    @test_throws ArgumentError remove_card_from_bottom!([])
end

@testset "Check size of the stack" begin
    @test  check_stack_size([1, 3, 4, 1], 4)
    @test !check_stack_size([1, 3, 4, 1, 5], 4)
    @test !check_stack_size([1, 3, 4], 4)
    @test  check_stack_size([], 0)
end

max_stars_repo_name: exercism-bot/julia

max_issues_repo_name: exercism-bot/julia

max_forks_repo_name: exercism-bot/julia


text:
@testset "Sampling" begin
    @testset "Basic consistency $f" for f in
        (:empty, :onetask, :chain, :shortlong, :medium)
        proj = getfield(Examples, f)

        samples = rand(proj, 10_000);

        # if we just did all tasks one after the other we get absolutely worst possible case
        worsecase = sum(maximum, proj.task_durations)
        # if we did everything in parallel we get absolutely best possible case
        # but it can never be shorter than the longest minimum duration of any single task
        bestcase = maximum(minimum, proj.task_durations)
        @test all(bestcase ≤ s ≤ worsecase for s in samples)
    end

    @testset "known distribution: proj $f, task $t" for (f,t) in (
        (:shortlong, :long), (:onetask, :a)
    )        # These projects have their duration fully controlled by 1 task
        # so should use same distribution.
        proj = getfield(Examples, f)
        expected_dist = proj.task_durations[t]

        # Make sure we would *not* reject the null hypothesios that it does indeed
        # come from this distrubition
        samples = rand(proj, 10_000);
        hypo_test = HypothesisTests.OneSampleADTest(samples, expected_dist)
        @test 0.05 < pvalue(hypo_test)
    end

    @testset "empty project" begin
        samples = rand(Examples.empty, 10_000)
        @test all(iszero, samples)
    end

    @testset "Static Length" begin
        proj = Project(
            (start=0, common=1, short=2, long=3, finish=0,),
            [
                :start => :common,
                :common => :short,
                :common =>:long,
                :short => :finish,
                :long => :finish
            ]
        )
        samples = rand(proj, 10_000)
        @test all(x->x==4, samples)
    end
end

@testset "Plotting" begin
    @testset "RecipeBase" begin
        proj = Examples.medium
        result, = RecipesBase.apply_recipe(Dict{Symbol, Any}(), proj)
        xs, = result.args
        @test length(xs) == 100_000
        @test xs isa Vector{Float64}
        @test result.plotattributes[:title] == "Completion Time"
    end

    @testset "Gadfly" begin
        proj = Examples.medium
        result = density(proj)

        xs = result.data.x
        @test length(xs) == 100_000
        @test xs isa Vector{Float64}

        geom = only(result.layers).geom
        @test geom isa Gadfly.Geom.LineGeometry
        @test geom.default_statistic isa Gadfly.Stat.DensityStatistic
    end
end


text:
# # Estimate Bernoulli draws probabilility

# We estimate a simple model of ``n`` independent Bernoulli draws, with
# probability ``α``. First, we load the packages we use.

using TransformVariables, LogDensityProblems, DynamicHMC, DynamicHMC.Diagnostics
using MCMCDiagnostics
using Parameters, Statistics, Random
import ForwardDiff              # use for AD

# Then define a structure to hold the data.
# For this model, the number of draws equal to `1` is a sufficient statistic.

"""
Toy problem using a Bernoulli distribution.

We model `n` independent draws from a ``Bernoulli(α)`` distribution.
"""
struct BernoulliProblem
    "Total number of draws in the data."
    n::Int
    "Number of draws `==1` in the data"
    s::Int
end

# Then make the type callable with the parameters *as a single argument*.  We
# use decomposition in the arguments, but it could be done inside the function,
# too.

function (problem::BernoulliProblem)(θ)
    @unpack α = θ               # extract the parameters
    @unpack n, s = problem      # extract the data
    ## log likelihood: the constant log(combinations(n, s)) term
    ## has been dropped since it is irrelevant to sampling.
    s * log(α) + (n-s) * log(1-α)
end

# We should test this, also, this would be a good place to benchmark and
# optimize more complicated problems.

p = BernoulliProblem(20, 10)
p((α = 0.5, ))

# Recall that we need to
#
# 1. transform from ``ℝ`` to the valid parameter domain `(0,1)` for more efficient sampling, and
#
# 2. calculate the derivatives for this transformed mapping.
#
# The helper packages `TransformVariables` and `LogDensityProblems` take care of
# this. We use a flat prior (the default, omitted)

t = as((α = as𝕀,))
P = TransformedLogDensity(t, p)
∇P = ADgradient(:ForwardDiff, P);

# Finally, we sample from the posterior. `chain` holds the chain (positions and
# diagnostic information), while the second returned value is the tuned sampler
# which would allow continuation of sampling.

results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000)

# To get the posterior for ``α``, we need to use `get_position` and
# then transform

posterior = transform.(t, results.chain);

# Extract the parameter.

posterior_α = first.(posterior);

# check the mean

mean(posterior_α)

# check the effective sample size

ess_α = effective_sample_size(posterior_α)

# NUTS-specific statistics

summarize_tree_statistics(results.tree_statistics)


text:
include("../utils.jl")

@testset "Axes2D" begin
    p = plot([1,2],[3,4])

    a  = gca()
    io = IOBuffer()
    show(io, MIME("text/plain"), a)
    s  = String(take!(io))

    @test isapproxstr(s, """
        GLEPlot.Axes2D
            Title:         none
            N. drawings:   1
            N. objects:    0
            Math mode:     false
            Layout origin: auto
        """)
end

@testset "Figure" begin
    p  = plot([1,2],[3,4])
    f  = gcf()

    io = IOBuffer()
    show(io, MIME("text/plain"), f)
    s  = String(take!(io))

    @test isapproxstr(s, """
        GLEPlot.Figure
            Name:          default ("_fig_")
            Size:          (12.0, 9.0)
            Bg. color:     white
            N. axes:       1
            LaTeX:         false
            Transparent:   false
        """
    )
end


text:
grade(bits::Unsigned) = count_ones(bits)

bits_scalar() = unsigned(0)
bits_first_of_grade(k) = (unsigned(1) << k) - unsigned(1)

bits_basis_vector(i) = unsigned(1) << (i - 1)
# bits_has_index(bits, i) = isone(bits >> (i - 1) & 1)

function bits_to_indices(bits)
	indices = Int[]
	i = 1
	while bits > 0
		if isone(bits & 1)
			push!(indices, i)
		end
		i += 1
		bits >>= 1
	end
	indices
end

function indices_to_bits(indices)
	bits = bits_scalar()
	for i ∈ indices
		bits += bits_basis_vector(i)
	end
	bits
end

"""
Return the smallest uint larger than the one given which has
the same number of binary ones.
Algorithm is [Gosper's hack](http://graphics.stanford.edu/~seander/bithacks.html#NextBitPermutation).

```
julia> GeometricAlgebra.next_bit_permutation(0b1011) |> bitstring
"00001101"
```
"""
function next_bit_permutation(v::Unsigned)
	I = one(v)
	t = v | (v - I)
	(t + I) | (((~t & -~t) - I) >> (trailing_zeros(v) + I))
end

struct FixedGradeBits
	grade::Int
end

function Base.iterate(fgb::FixedGradeBits, bits=bits_first_of_grade(fgb.grade))
	(bits, next_bit_permutation(bits))
end

Base.IteratorSize(::FixedGradeBits) = Base.IsInfinite()

"""
	bits_of_grade(k[, n])

Generate basis blade bits of grade `k` in ascending order.
Yields all basis blades in the dimension `n`, if given, otherwise iterate indefinitely.

Example
===
```jldoctest
julia> bits_of_grade(2, 4) .|> UInt8 .|> bitstring
```
"""

bits_of_grade(k) = FixedGradeBits(k)
bits_of_grade(k, n) = Iterators.take(FixedGradeBits(k), binomial(n, k))

"""
	bits_to_indices(bits)

Return the positions of the ones in the unsigned integer `bits`.

Used to convert between representations of a unit blade.
See also [`$(repr(indices_to_bits))`](@ref).

Examples
===
```jldoctest
julia> GeometricAlgebra.bits_to_indices(0b1001101)
4-element Vector{Int64}:
 1
 3
 4
 7
```
"""
bits_to_indices

"""
	indices_to_bits(indices)

Create unsigned integer with bits at the positions given in the vector `indices`.

Used to convert between representations of a unit blade.
See also [`$(repr(bits_to_indices))`](@ref).

Examples
===
```jldoctest
julia> GeometricAlgebra.indices_to_bits([1, 2, 5]) |> UInt16 |> bitstring
"0000000000010011"
```

"""
indices_to_bits

# # this version is 10% slower
# function bits_to_linear_index(bits::Unsigned)
# 	i = 1
# 	for b ∈ FixedGradeBits(grade(bits))
# 		if b >= bits
# 			return i
# 		end
# 		i += 1
# 	end
# end
function bits_to_linear_index(bits::Unsigned)
	i = 1
	u = bits_first_of_grade(grade(bits))
	while u < bits
		u = next_bit_permutation(u)
		i += 1
	end
	i
end
function linear_index_to_bits(i, k)
	bits = bits_scalar()
	for b ∈ Iterators.take(FixedGradeBits(k), i)
		bits = b
	end
	bits
end

"""
Compute sign flips of blade product due to transposing basis vectors.
(The full sign of the product will also depend on the basis norms.)
"""
function sign_from_swaps(a::Unsigned, b::Unsigned)
	swaps = 0
	while b > 0
		a >>= 1
		swaps += count_ones(a)*(b & 1)
		b >>= 1
	end
	iseven(swaps) ? 1 : -1
end

"""
Compute the overall factor arising from the geometric product between
repeated basis vectors.
"""
function factor_from_squares(sig, squares)
	factor = 1
	bv = 1
	while squares > 0
		if isone(squares & 1)
			factor *= sig[bv]
		end
		squares >>= 1
		bv += 1
	end
	factor
end

function geometric_prod_bits(sig, a::Unsigned, b::Unsigned)
	factor = sign_from_swaps(a, b)*factor_from_squares(sig, a & b)
	bits = a ⊻ b
	factor, bits
end