text:
module MonthlyDates
    using Printf
    using Dates
    import Dates: UTInstant, value, DatePart
    using RecipesBase: RecipesBase, @recipe
    # Quarter defined in Julia 1.6
    if isdefined(Dates, :Quarter)
        import Dates: Quarter
    else
        include("Quarter.jl")
    end

    ##############################################################################
    ##
    ## MonthlyDate
    ##
    ##############################################################################
    """
    `MonthlyDate`

    `MonthlyDate` wraps a `UTInstant{Month}` and interprets it according to the proleptic
    Gregorian calendar.
    """
    struct MonthlyDate <: TimeType
    	instant::UTInstant{Month}
        MonthlyDate(v::UTInstant{Month}) = new(v)
    end
    UTm(x) = UTInstant(Month(x))

    """
    `MonthlyDate(y, [m]) -> MonthlyDate`

    Construct a `MonthlyDate` type by parts.
    """
    function MonthlyDate(y::Int, m::Int = 1)
        1 <= m <= 12 || throw(ArgumentError("Month: $m out of range (1:12)"))
        MonthlyDate(UTm(12 * (y - 1) + m))
    end
    MonthlyDate(y::Year, m::Month = Month(1)) = MonthlyDate(value(y), value(m))
    MonthlyDate(y, m = 1) = MonthlyDate(Int64(y), Int64(m))

    """
    `MonthlyDate(dt::Date) -> MonthlyDate`

    Convert a `Date` to a `MonthlyDate`
    """
    MonthlyDate(dt::TimeType) = convert(MonthlyDate, dt)

    Base.convert(::Type{MonthlyDate}, dt::Date) = MonthlyDate(yearmonth(dt)...)
    Base.convert(::Type{MonthlyDate}, dt::DateTime) = MonthlyDate(yearmonth(dt)...)
    Base.convert(::Type{MonthlyDate}, x::Month) = MonthlyDate(UTInstant(x))

    Dates.yearmonth(dt::MonthlyDate) =  (1, 1) .+ divrem(value(dt) - 1, 12)
    Base.convert(::Type{Date}, dt::MonthlyDate) = Date(yearmonth(dt)...)
    Base.convert(::Type{DateTime}, dt::MonthlyDate) = DateTime(yearmonth(dt)...)
    Base.convert(::Type{Month}, dt::MonthlyDate) = Month(value(dt))

    Base.promote_rule(::Type{MonthlyDate}, x::Type{Date}) = Date
    Base.promote_rule(::Type{MonthlyDate}, x::Type{DateTime}) = DateTime

    Dates.eps(::Type{MonthlyDate}) = Month(1)
    Dates.zero(::Type{MonthlyDate}) = Month(0)

    #accessors (only bigger periods)
    Dates.month(dt::MonthlyDate) = 1 + rem(value(dt) - 1, 12)
    Dates.year(dt::MonthlyDate) =  1 + div(value(dt) - 1, 12)
    Dates.quarterofyear(dt::MonthlyDate) = 1 + div(month(dt) - 1, 3)

    Dates.Month(dt::MonthlyDate) = Month(month(dt))
    Quarter(dt::MonthlyDate) = Quarter(quarterofyear(dt))
    Dates.Year(dt::MonthlyDate) = Year(year(dt))

    # arithmetics
    Base.:+(dt::MonthlyDate, m::Month) = MonthlyDate(UTm(value(dt) + value(m)))
    Base.:-(dt::MonthlyDate, m::Month) = MonthlyDate(UTm(value(dt) - value(m)))
    Base.:+(dt::MonthlyDate, p::Quarter) = dt + Month(p)
    Base.:-(dt::MonthlyDate, p::Quarter) = dt - Month(p)
    Base.:+(dt::MonthlyDate, p::Year) = dt + Month(p)
    Base.:-(dt::MonthlyDate, p::Year) = dt - Month(p)

    # range
    function Dates.guess(a::MonthlyDate, b::MonthlyDate, c::Month)
    	Int64(div(value(b - a), value(c)))
    end
    function Dates.guess(a::MonthlyDate, b::MonthlyDate, c::Quarter)
    	Int64(div(value(b - a), 3 * value(c)))
    end
    function Dates.guess(a::MonthlyDate, b::MonthlyDate, c::Year)
    	Int64(div(value(b - a), 12 * value(c)))
    end

    # adjusters
    Dates.firstdayofquarter(dt::MonthlyDate) = Dates.firstdayofquarter(Date(dt))
    Dates.lastdayofquarter(dt::MonthlyDate) = Dates.lastdayofquarter(Date(dt))
    Dates.firstdayofmonth(dt::MonthlyDate) = Date(dt)
    Dates.lastdayofmonth(dt::MonthlyDate) = Dates.lastdayofmonth(Date(dt))

    # parse
    function Base.parse(::Type{MonthlyDate}, s::AbstractString, df::DateFormat)
        MonthlyDate(parse(Date, s, df))
    end
    const MonthlyFormat = dateformat"yyyy-mm"
    Dates.default_format(::Type{MonthlyDate}) = MonthlyFormat

    function MonthlyDate(d::AbstractString, format::AbstractString;
        locale::Dates.Locale = Dates.ENGLISH)
        parse(MonthlyDate, d, DateFormat(format, locale))
    end

    function MonthlyDate(d::AbstractString, format::DateFormat = ISODateFormat)
        parse(MonthlyDate, d, format)
    end

    # show
    function Base.print(io::IO, dt::MonthlyDate)
        y,m = yearmonth(dt)
        yy = y < 0 ? @sprintf("%05i", y) : lpad(y, 4, "0")
        mm = lpad(m, 2, "0")
        print(io, "$(yy)-$mm")
    end

    Base.show(io::IO, ::MIME"text/plain", dt::MonthlyDate) = print(io, dt)
    if VERSION >= v"1.5-"
        Base.show(io::IO, dt::MonthlyDate) = print(io, MonthlyDate, "(\"", dt, "\")")
    else
        Base.show(io::IO, dt::MonthlyDate) = print(io, dt)
    end
    if VERSION >= v"1.4-"
        Base.typeinfo_implicit(::Type{MonthlyDate}) = true
    end

    #Plot
    # Issue is that I would like to choose ticks myself. Or at least convert depending on whether all ticks are Integer or not, which is not possible
    @recipe function f(::Type{MonthlyDate}, dt::MonthlyDate)
	   (value, dt -> string(MonthlyDate(UTm(round(dt)))))
	end

	export MonthlyDate
    ##############################################################################
    ##
    ## QuarterlyDate
    ##
    ##############################################################################
    """
    QuarterlyDate

    QuarterlyDate wraps a UTInstant{Quarter} and interprets it according to the proleptic
    Gregorian calendar.
    """
    struct QuarterlyDate <: TimeType
    	instant::UTInstant{Quarter}
        QuarterlyDate(v::UTInstant{Quarter}) = new(v)
    end
    UTQ(x) = UTInstant(Quarter(x))

    """
    `QuarterlyDate(y, [q]) -> QuarterlyDate`

    Construct a `QuaterlyDate` type by parts.
    """
    function QuarterlyDate(y::Int, q::Int = 1)
        1 <= q <= 4 || throw(ArgumentError("Quarter: $q out of range (1:4)"))
        QuarterlyDate(UTQ(4 * (y - 1) + q))
    end
    QuarterlyDate(y::Year, q::Quarter = Quarter(1)) = QuarterlyDate(value(y), value(q))
    QuarterlyDate(y, q = 1) = QuarterlyDate(Int64(y), Int64(q))

    """
    `QuarterlyDate(dt::Date) -> QuarterlyDate`

    Convert a `Date` to a `QuarterlyDate`
    """
    QuarterlyDate(dt::TimeType) = convert(QuarterlyDate, dt)
    Base.convert(::Type{QuarterlyDate}, dt::MonthlyDate) = QuarterlyDate(UTQ(((value(dt) - 1) ÷ 3 + 1)))
    Base.convert(::Type{QuarterlyDate}, dt::Date) = QuarterlyDate(year(dt), quarterofyear(dt))
    Base.convert(::Type{QuarterlyDate}, dt::DateTime) = QuarterlyDate(year(dt), quarterofyear(dt))
    Base.convert(::Type{QuarterlyDate}, x::Quarter) = QuarterlyDate(UTInstant(x))

    function Dates.yearmonth(dt::QuarterlyDate)
        y, q = divrem(value(dt) - 1, 4)
        return 1 + y, 1 + 3 * q
    end
    Base.convert(::Type{MonthlyDate}, dt::QuarterlyDate) = MonthlyDate(UTm(((value(dt) - 1) * 3 + 1)))
    Base.convert(::Type{Date}, dt::QuarterlyDate) = Date(yearmonth(dt)...)
    Base.convert(::Type{DateTime}, dt::QuarterlyDate) = DateTime(yearmonth(dt)...)
    Base.convert(::Type{Quarter}, dt::QuarterlyDate) =  Quarter(value(dt))

    Base.promote_rule(::Type{QuarterlyDate}, x::Type{MonthlyDate}) = MonthlyDate
    Base.promote_rule(::Type{QuarterlyDate}, x::Type{Date}) = Date
    Base.promote_rule(::Type{QuarterlyDate}, x::Type{DateTime}) = DateTime

    Dates.eps(::Type{QuarterlyDate}) = Quarter(1)
    Dates.zero(::Type{QuarterlyDate}) = Quarter(0)

    #accessors (only bigger periods)
    Dates.month(dt::QuarterlyDate) = 1 + 3 * rem(value(dt) - 1, 4)
    Dates.year(dt::QuarterlyDate) = 1 + div(value(dt) - 1, 4)
    Dates.quarterofyear(dt::QuarterlyDate) = 1 + rem(value(dt) - 1, 4)
    Quarter(dt::QuarterlyDate) = Quarter(quarterofyear(dt))
    Dates.Year(dt::QuarterlyDate) = Year(year(dt))

    # arithmetics
    Base.:+(dt::QuarterlyDate, q::Quarter) = QuarterlyDate(UTQ(value(dt) + value(q)))
    Base.:-(dt::QuarterlyDate, q::Quarter) = QuarterlyDate(UTQ(value(dt) - value(q)))
    Base.:+(dt::QuarterlyDate, p::Year) = dt + Quarter(p)
    Base.:-(dt::QuarterlyDate, p::Year) = dt - Quarter(p)

    # range
    function Dates.guess(a::QuarterlyDate, b::QuarterlyDate, c::Quarter)
    	Int64(div(value(b - a), value(c)))
    end
    function Dates.guess(a::QuarterlyDate, b::QuarterlyDate, c::Year)
    	Int64(div(value(b - a), 4 * value(c)))
    end

    # adjusters
    Dates.firstdayofquarter(dt::QuarterlyDate) = Date(dt)
    Dates.lastdayofquarter(dt::QuarterlyDate) = Dates.lastdayofquarter(Date(dt))

    # parse
    function Dates.tryparsenext(d::DatePart{'q'}, str, i, len)
        return Dates.tryparsenext_base10(str, i, len, Dates.min_width(d), Dates.max_width(d))
    end

    function Dates.format(io, d::DatePart{'q'}, dt)
        print(io, string(quarterofyear(dt)))
    end

    function parse_quarterly(str::AbstractString)
        if occursin('Q', str) # if the string has 'Q', try the default format "2018-Q1"
            parse_quarterly(str, QuarterlyFormat)
        else # else assume it's formatted as Date
            parse_quarterly(str, ISODateFormat)
        end
    end

    function parse_quarterly(str::AbstractString, df::T) where {T<:DateFormat}
        if Dates.DatePart{'q'} ∉ Dates._directives(T) # if
            return QuarterlyDate(parse(Date, str, df))
        else
            return parse(QuarterlyDate, str, df)
        end
    end

    # the new key 'q' to Dates.CONVERSION_SPECIFIERS not added during precompilation
    # therefore construction of QuarterlyFormat with from dateformat string will fail
    # manually construct QuarterlyFormat
    const QuarterlyFormat = DateFormat{Symbol("yyyy-Qm"), Tuple{Dates.DatePart{'y'}, Dates.Delim{String, 2}, Dates.DatePart{'q'}}}(
        (
            Dates.DatePart{'y'}(4,false), Dates.Delim{String,2}("-Q"), Dates.DatePart{'q'}(1, false)),
            Dates.ENGLISH
        )
    Dates.default_format(::Type{QuarterlyDate}) = QuarterlyFormat

    function QuarterlyDate(d::AbstractString, format::AbstractString;
        locale::Dates.Locale = Dates.ENGLISH)
        parse_quarterly(d, DateFormat(format, locale))
    end

    QuarterlyDate(d::AbstractString, format::DateFormat) = parse_quarterly(d, format)
    QuarterlyDate(d::AbstractString) = parse_quarterly(d)

    # show
    function Base.print(io::IO, dt::QuarterlyDate)
        y,m = yearmonth(dt)
        yy = y < 0 ? @sprintf("%05i", y) : lpad(y, 4, "0")
        q = (m - 1) ÷ 3 + 1
        print(io, "$(yy)-Q$q")
    end
    Base.show(io::IO, ::MIME"text/plain", dt::QuarterlyDate) = print(io, dt)
    if VERSION >= v"1.5-"
        Base.show(io::IO, dt::QuarterlyDate) = print(io, QuarterlyDate, "(\"", dt, "\")")
    else
        Base.show(io::IO, dt::QuarterlyDate) = print(io, dt)
    end

    if VERSION >= v"1.4-"
        Base.typeinfo_implicit(::Type{QuarterlyDate}) = true
    end

    #Plot
    @recipe function f(::Type{QuarterlyDate}, x::QuarterlyDate)
        (value, dt -> string(QuarterlyDate(UTQ(round(dt)))))
    end

    export QuarterlyDate

    # executed at runtime to avoid issues with precompiling dicts
    function __init__()
        Dates.CONVERSION_SPECIFIERS['q'] = Quarter
        Dates.CONVERSION_DEFAULTS[Quarter] = Int64(1)
        Dates.CONVERSION_TRANSLATIONS[QuarterlyDate] = (Year, Quarter)
    end

end


text:
export @roc, rocconvert, rocfunction

struct Kernel{F,TT}
    agent::HSAAgent
    mod::ROCModule
    fun::ROCFunction
end

# split keyword arguments to `@roc` into ones affecting the compiler,
# the execution, or both
function split_kwargs(kwargs)
    compiler_kws = [:agent, :queue]
    call_kws     = [:groupsize, :gridsize, :agent, :queue]
    compiler_kwargs = []
    call_kwargs = []
    for kwarg in kwargs
        if Meta.isexpr(kwarg, :(=))
            key,val = kwarg.args
            if isa(key, Symbol)
                if !(key in compiler_kws) && !(key in call_kws)
                    throw(ArgumentError("unknown keyword argument '$key'"))
                end
                if key in compiler_kws
                    push!(compiler_kwargs, kwarg)
                end
                if key in call_kws
                    push!(call_kwargs, kwarg)
                end
            else
                throw(ArgumentError("non-symbolic keyword '$key'"))
            end
        else
            throw(ArgumentError("non-keyword argument like option '$kwarg'"))
        end
    end

    return compiler_kwargs, call_kwargs
end

# assign arguments to variables, handle splatting
function assign_args!(code, args)
    # handle splatting
    splats = map(arg -> Meta.isexpr(arg, :(...)), args)
    args = map(args, splats) do arg, splat
        splat ? arg.args[1] : arg
    end

    # assign arguments to variables
    vars = Tuple(gensym() for arg in args)
    map(vars, args) do var,arg
        push!(code.args, :($var = $(esc(arg))))
    end

    # convert the arguments, compile the function and call the kernel
    # while keeping the original arguments alive
    var_exprs = map(vars, args, splats) do var, arg, splat
         splat ? Expr(:(...), var) : var
    end

    return vars, var_exprs
end

# fast lookup of global world age
world_age() = ccall(:jl_get_tls_world_age, UInt, ())

# slow lookup of local method age
function method_age(f, tt)::UInt
    for m in Base._methods(f, tt, 1, typemax(UInt))
        @static if VERSION >= v"1.2.0-DEV.573"
            return m[3].primary_world
        else
            return m[3].min_world
        end
    end
    throw(MethodError(f, tt))
end

macro roc(ex...)
    @assert length(ex) >= 1 "Not enough arguments provided"
    call = last(ex)
    kwargs = ex[1:end-1]
    @assert call.head == :call "Invalid kernel specification"
    f = call.args[1]
    args = call.args[2:end]
    code = quote end
    compiler_kwargs, call_kwargs = split_kwargs(kwargs)
    vars, var_exprs = assign_args!(code, args)
    push!(code.args,
        quote
            GC.@preserve $(vars...) begin
                local kernel_args = rocconvert.(($(var_exprs...),))
                local kernel_tt = Tuple{Core.Typeof.(kernel_args)...}
                local agent = get_default_agent()
                local kernel = rocfunction(agent, $(esc(f)), kernel_tt;
                                           $(map(esc, compiler_kwargs)...))
                local queue = get_default_queue(agent)
                local signal = HSASignal()
                kernel(queue, signal, kernel_args...; $(map(esc, call_kwargs)...))
                wait(signal)
            end
        end)
    return code
end

## adaptors

struct Adaptor end

# Base.RefValue isn't GPU compatible, so provide a compatible alternative
mutable struct ROCRefValue{T} <: Ref{T}
    x::T
end
Base.getindex(r::ROCRefValue) = r.x
Base.setindex!(r::ROCRefValue{T}, value::T) where T = (r.x = value;)
Adapt.adapt_structure(to::Adaptor, r::Base.RefValue) = ROCRefValue(adapt(to, r[]))

## interop with HSAArray

function Base.convert(::Type{ROCDeviceArray{T,N,AS.Global}}, a::HSAArray{T,N}) where {T,N}
    ROCDeviceArray{T,N,AS.Global}(a.size, DevicePtr{T,AS.Global}(a.handle))
end

Adapt.adapt_storage(::Adaptor, xs::HSAArray{T,N}) where {T,N} =
    convert(ROCDeviceArray{T,N,AS.Global}, xs)

# convenience function
"""
    rocconvert(x)

This function is called for every argument to be passed to a kernel, allowing
it to be converted to a GPU-friendly format. By default, the function does
nothing and returns the input object `x` as-is.

Do not add methods to this function, but instead extend the underlying
Adapt.jl package and register methods for the the `AMDGPUnative.Adaptor` type.
"""
rocconvert(arg) = adapt(Adaptor(), arg)

## APIs for manual compilation

const agecache = Dict{UInt, UInt}()
const compilecache = Dict{UInt, Kernel}()

"""
    rocfunction(f, tt=Tuple{}; kwargs...)

Low-level interface to compile a function invocation for the currently-active
GPU, returning a callable kernel object. For a higher-level interface, use
[`@roc`](@ref).

Currently, no keyword arguments are implemented.

The output of this function is automatically cached, i.e. you can simply call
`rocfunction` in a hot path without degrading performance. New code will be
generated automatically, when the function changes, or when different types or
keyword arguments are provided.
"""
@generated function rocfunction(agent::HSAAgent, f::Core.Function, tt::Type=Tuple{}; kwargs...)
    tt = Base.to_tuple_type(tt.parameters[1])
    sig = Base.signature_type(f, tt)
    t = Tuple(tt.parameters)

    precomp_key = hash(sig)  # precomputable part of the keys
    quote
        Base.@_inline_meta

        # look-up the method age
        key = hash(world_age(), $precomp_key)
        if haskey(agecache, key)
            age = agecache[key]
        else
            age = method_age(f, $t)
            agecache[key] = age
        end
        key = hash(age, key)

        # compile the function
        key = hash(kwargs, key)
        for nf in 1:nfields(f)
            # mix in the values of any captured variable
            key = hash(getfield(f, nf), key)
        end
        if !haskey(compilecache, key)
            mod, fun = compile(agent, f, tt; kwargs...)
            kernel = Kernel{f,tt}(agent, mod, fun)
            compilecache[key] = kernel
        end

        return compilecache[key]::Kernel{f,tt}
    end
end
rocfunction(f::Core.Function, tt::Type=Tuple{}; kwargs...) =
    rocfunction(get_default_agent(), f, tt; kwargs...)

@generated function (kernel::Kernel{F,TT})(queue::HSAQueue, signal::HSASignal, args...; call_kwargs...) where {F,TT}
    sig = Base.signature_type(F, TT)
    args = (:F, (:( args[$i] ) for i in 1:length(args))...)

    # filter out ghost arguments that shouldn't be passed
    to_pass = map(!isghosttype, sig.parameters)
    call_t =                  Type[x[1] for x in zip(sig.parameters,  to_pass) if x[2]]
    call_args = Union{Expr,Symbol}[x[1] for x in zip(args, to_pass)            if x[2]]

    # replace non-isbits arguments (they should be unused, or compilation would have failed)
    # alternatively, make HSARuntime allow `launch` with non-isbits arguments.
    for (i,dt) in enumerate(call_t)
        if !isbitstype(dt)
            call_t[i] = Ptr{Any}
            call_args[i] = :C_NULL
        end
    end

    # finalize types
    call_tt = Base.to_tuple_type(call_t)

    quote
        Base.@_inline_meta

        roccall(queue, signal, kernel.fun, $call_tt, $(call_args...); call_kwargs...)
    end
end

# There doesn't seem to be a way to access the documentation for the
# call-syntax, so attach it to the type
"""
    (::Kernel)(args...; kwargs...)

Low-level interface to call a compiled kernel, passing GPU-compatible
arguments in `args`.  For a higher-level interface, use [`@roc`](@ref).

The following keyword arguments are supported:
- groupsize (defaults to 1)
- gridsize (defaults to 1)
"""
Kernel

ext: jl


text:
using BenchmarkTools, BlockDiagonals, LinearAlgebra, Random, TemporalGPs

rng = MersenneTwister(123456);
P = 50;
Q = 150;

# Matrix-Matrix multiplies.

A = BlockDiagonal([randn(rng, P, P) for _ in 1:3]);
A_dense = collect(A);
At = BlockDiagonal(map(collect ∘ transpose, blocks(A)));
At_dense = collect(At);
B = randn(rng, Q, Q);
C = randn(rng, Q, Q);

@benchmark mul!($C, $A, $B, 1.0, 1.0)
@benchmark mul!($C, $A_dense, $B, 1.0, 1.0)

@benchmark mul!($C, $At', $B, 1.0, 1.0)
@benchmark mul!($C, $At_dense', $B, 1.0, 1.0)

@benchmark mul!($C, $A, $B', 1.0, 1.0)
@benchmark mul!($C, $A_dense, $B', 1.0, 1.0)

@benchmark mul!($C, $At', $B', 1.0, 1.0)
@benchmark mul!($C, $At_dense', $B', 1.0, 1.0)

@benchmark mul!($C, $B, $A, 1.0, 1.0)
@benchmark mul!($C, $B, $A_dense, 1.0, 1.0)

@benchmark mul!($C, $B, $At', 1.0, 1.0)
@benchmark mul!($C, $B, $At_dense', 1.0, 1.0)

@benchmark mul!($C, $B', $A, 1.0, 1.0)
@benchmark mul!($C, $B', $A_dense, 1.0, 1.0)

@benchmark mul!($C, $B', $At', 1.0, 1.0)
@benchmark mul!($C, $B', $At_dense', 1.0, 1.0)

# Matrix-Vector multiplies.

b = randn(rng, Q);
c = randn(rng, Q);

@benchmark mul!($c, $A, $b, 1.0, 1.0)
@benchmark mul!($c, $A_dense, $b, 1.0, 1.0)


text:
# generate examples
using Literate: Literate

EXAMPLE_DIR = joinpath(@__DIR__, "src", "literate")
GENERATED_DIR = joinpath(@__DIR__, "src", "examples")
mkpath(GENERATED_DIR)
for example in readdir(EXAMPLE_DIR)
    if endswith(example, ".jl")
        input = abspath(joinpath(EXAMPLE_DIR, example))
        script = Literate.script(input, GENERATED_DIR)
        code = strip(read(script, String))
        mdpost(str) = replace(str, "@__CODE__" => code)
        Literate.markdown(input, GENERATED_DIR; postprocess=mdpost)
        Literate.notebook(input, GENERATED_DIR; execute=true)
    elseif any(endswith.(example, [".png", ".jpg", ".gif"]))
        cp(joinpath(EXAMPLE_DIR, example), joinpath(GENERATED_DIR, example); force=true)
    else
        @warn "ignoring $example"
    end
end

# remove any .vtu files in the generated dir (should not be deployed)
cd(GENERATED_DIR) do
    foreach(file -> endswith(file, ".vtu") && rm(file), readdir())
end


text:
if is_unix()
    suffix = is_apple() ? "dylib" : "so"
    cd(joinpath(dirname(@__FILE__), "src"))
    try
        run(`make FC=ifort SUFFIX=$suffix`)
    catch
        run(`make FC=gfortran SUFFIX=$suffix`)
    end
else
    error("windows currently unsupported")
end


text:
# ------- distribution functions on the invariants --------------

function rdf(at::Atoms, rcut, transform=identity)
   if transform != identity
      return idf(2, at, rcut, transform)[1][1]
   end
   nlist = neighbourlist(at, rcut)
   return nlist.r
end

"""
`idf(N::Integer, at, rcut, transform)`

invariants distribution function => accumulates all the invariants values
arising during an N-body assembly over at.
"""
idf(N::Integer, at, rcut, transform) =
      _idf(Val(N), Val(bo2edges(N)), at, rcut, transform)

function _idf(valN::Val{N}, valM::Val{M}, at::Atoms{T}, rcut::T,
              transform) where {N, M, T}
   # compute invariants vectors to learn how many there are
   x = rand(SVector{M,T})
   I1, I2 = invariants(x)

   I1acc = [ T[] for n = 1:length(I1) ]
   I2acc = [ T[] for n = 1:length(I2) ]
   Iacc = (I1acc, I2acc)

   for (i, j, r, R) in sites(at, rcut)
      eval_site_nbody!(valN, R, rcut,
                  (Iacc, s, _1, _2, _3) -> idf_accumulator(Iacc, transform.(s)),
                  Iacc, nothing)
   end

   return I1acc, I2acc
end

function idf_accumulator(Iacc, x)
   I1, I2 = invariants(x)
   for n = 1:length(I1)
      push!(Iacc[1][n], I1[n])
   end
   for n = 1:length(I2)
      push!(Iacc[2][n], I2[n])
   end
   return Iacc
end

# # =============== Experimental:
# #   evaluate NBodyIP
#
# (V::NBodyIP)(args...) = evaluate(V, args...)
#
# evaluate(V::NBodyIP, r::Number) = evaluate(V::NBodyIP, SVector(r))
#
# evaluate(V::NBodyIP, r1::T, r2::T, r3::T) where {T <: Number} =
#       evaluate(V::NBodyIP, SVector(r1, r2, r3))
#
# function evaluate(V::NBodyIP, r::SVector{N, T}) where {N, T}
#    v = zero(T)
#    for Vn in V.components
#       if bo2edges(bodyorder(Vn)) == N
#          v += Vn(r)
#       end
#    end
#    return v
# end

# dim(V::NBPoly{N,M}) where {N, M} = M-1

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
# lplagorder.m
# Luca Brugnolini 2014
# Based on Kilian
# This program selects the lag order for local projections using AIC, SIC , AICC and HCQ

function lplagorder(z,pbar,H)

# Inputs : data, maximum lag order to consider, IRFs horizons
# Outputs: lag-orders for H horizons selected by AIC and SIC (Hx1 vectors)

(T,K) = size(z);
t     = T-pbar;
paic  = zeros(H,1);
# psic  = zeros(H,1);                       # FOR SIC
# paicc = zeros(H,1);                       # FOR HURVICH AND TSAI
# phqc  = zeros(H,1);                       # FOR HQ
Ys = z[pbar+1:T,:];
Yt = z[pbar:T-1,:];                  # RHS variable of interest: y(t)
X = ones(T-pbar,1);                  # constant term

    for j=2:pbar
       X = [X z[pbar+1-j:T-j,:]];           # other RHS variables: y(t-1)~y(t-p+1)

    end

    for j=1:H                                          # loop for horizon h of IRF
        AIC = zeros(pbar,1);                       # the vector of AIC(p)
        # SIC = zeros(pbar,1);                     # the vector of SIC(p)
        # AICC = zeros(pbar,1);                    # FOR HURVICH AND TSAI
        # HQC = zeros(pbar, 1);                    # FOR HQ
        ys = Ys[j:end,:];                          # dependent variable: yt+s
        yt = Yt[1:end-j+1,:];                      # regressor of interest: yt

    for m=1:pbar                             # loop for lag order selection
        x = X[1:end-j+1,1:(m-1)*K+1];        # other independent variables
        Mx = eye(size(x,1))-x/(x'*x)*x';     # annhilation matrix **
        B1 = (yt'*Mx*yt)\(yt'*Mx*ys);        # IRF by Local Projection **
        U = Mx*ys - Mx*yt*B1;                # residual from LPs
        SIG = U'*U/t;                        # variance of errors
        # b = t/(t-(m*K+K+1));
        AIC[m]  = log(det(SIG)) + 2*(K^2*m)/t;                        # AIC statistic
        # SIC[m]  = log(det(SIG)) + (K^2*m)*log(t)/t;                 # SIC statistic
        # AICC[m] = t*(log(det(SIG))+K) + 2*b*(K^2*m+K*(K+1)/2);      # FOR HURVICH AND TSAI
        # HQC[m]  = log(det(SIG)) + 2*log(log(t))*K^2*m/t;            # FOR HQ
    end

        paic[j] = indmin(AIC)                      # paic is the lag-order by AIC
        # psic[j]  = indmin(SIC);                  # paic is the lag-order by AIC
        # paicc[j] = indmin(AICC);                 # FOR HURVICH AND TSAI
        # phqc[j]  = indmin(HQC);                  # FOR HQ

    end

return paic

end


text:
using BenchmarkTools
using Distributions
using GLM
using HypothesisTests
using LinearAlgebra
using Plots
using Random
using RDatasets
using SpecialFunctions
using Statistics

# Ridge regression

n = 10000
m = 1000

Random.seed!(666)
X = randn(n, m)

y = 10*X[:,1] + X[:,2] + randn(n)

# Exercise

# Comparison

@btime ridge_reg(X, y, 10);

@btime ridge_reg(X, y, 10, Q, Q_inv, λ);

# Dependence on mu

μs = range(0, 1000; length=50)

ws = hcat(ridge_reg.(Ref(X), Ref(y), μs, Ref(Q), Ref(Q_inv), Ref(λ))...)

plot(μs, abs.(ws');
    label="",
    yscale=:log10,
    xlabel="mu",
    ylabel="weights: log scale",
)

# LASSO

S(x, η) = max(x-η, 0) - max(-x-η, 0)

function lasso(X, y, μ, Q, Q_inv, λ;
        max_iter = 100,
        ρ = 1e3,
        w = zeros(size(X,2)),
        u = zeros(size(X,2)),
        z = zeros(size(X,2)),
    )

    for i in 1:max_iter
        w = Q * ( (Diagonal(1 ./ (λ .+ ρ)) * ( Q_inv * (X'*y + ρ*(z-u)))))
        z = S.(w + u, μ / ρ)
        u = u + w - z
    end
    return w, u, z
end

ws = zeros(size(X,2), length(μs))

for (i, μ) in enumerate(μs)
    global w, u, z
    w, u, z = i > 1 ? lasso(X, y, μ, Q, Q_inv, λ; w, u, z) : lasso(X, y, μ, Q, Q_inv, λ)
    ws[:,i] = w
end

plot(μs, abs.(ws');
    label="",
    yscale=:log10,
    xlabel="mu",
    ylabel="weights: log scale",
)

plot(0:0.1:10, gamma;
    xlabel="x",
    ylabel="gamma(x): log scale",
    label="",
    yscale=:log10,
)

# Exercise

# Exercise

# Volume of a ball

ms = 1:15
ns = Int64.([1e1; 1e3; 1e5])

Random.seed!(666)

plt = plot(ms, m -> volume_true(m, 1);
    xlabel="dimension",
    ylabel="unit ball volume",
    legend=:topleft,
    label="True",
    line=(4,:black),
)

for n in ns
    plot!(plt, ms, m -> volume_monte_carlo.(m; n=n); label="n = $n")
end

display(plt)

# Generating from the uniform distributions

rand(Uniform(-1, 1), 10, 5)

# Sampling from distributions

d1 = Normal()
d2 = Normal(1, 1)
d3 = Normal(0, 0.01)

f1(x) = pdf(d1, x)
f2(x) = pdf(d2, x)
f3(x) = pdf(d3, x)

function plot_histogram(xs, f; kwargs...)
    plt = histogram(xs;
        label="Sampled density",
        xlabel = "x",
        ylabel = "pdf(x)",
        nbins = 85,
        normalize = :pdf,
        opacity = 0.5,
        kwargs...
    )

    plot!(plt, range(minimum(xs), maximum(xs); length=100), f;
        label="True density",
        line=(4,:black),
    )

    return plt
end

plot_histogram(rand(d1, 1000000), f1)

# Exercise

# Rejection sampling

xlims = (-10, 10)

for (f, d) in zip((f1, f2, f3), (d1, d2, d3))
    Random.seed!(666)
    xs = rejection_sampling(f, f(d.μ), xlims...)

    pl = plot_histogram(xs, f)
    display(pl)
end

# Exercise

# Hypothesis testing

Random.seed!(666)

n = 1000
xs = randn(n)

# Exercise

# One-sample t-test

OneSampleTTest(xs)

# Linear models

wages = dataset("plm", "Snmesp")
wages.W = 2. .^ (wages.W)

X = Matrix(wages[:, [:N, :Y, :I, :K, :F]])
X = hcat(ones(size(X,1)), X)
y = wages[:, :W]

w0 = (X'*X) \ (X'*y)

model = lm(@formula(W ~ 1 + N + Y + I + K + F), wages)

# Exercise

# Plot histograms

y_hat = predict(model)

plot_histogram(y_hat, x -> pdf(Normal(mean(y_hat), std(y_hat)), x))

plot_histogram(y_hat, x -> pdf(fit(Normal, y_hat), x))

test_normality = ExactOneSampleKSTest(y_hat, Normal(mean(y_hat), std(y_hat)))

# Generalized linear models

model = glm(@formula(W ~ 1 + N + Y + I + K + F), wages, InverseGaussian(), SqrtLink())

# Exercise

max_stars_repo_name: petrmvala/JuliaCourse

max_issues_repo_name: petrmvala/JuliaCourse

max_forks_repo_name: petrmvala/JuliaCourse


text:
using PyPlot
using Statistics
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCellDG
include("../../../../test/useful_routines.jl")
include("../transformation-elasticity-solver.jl")

TES = TransformationElasticitySolver

function perturbation(x, frequency, amplitude)
    return amplitude * sin.(2 * pi * frequency * x)
end

perturbed_distancefunction(x, initialposition, frequency, amplitude) =
    plane_distance_function(x, [1.0, 0.0], [initialposition, 0.0]) +
    perturbation(x[2, :], frequency, amplitude)

function perturbed_normal(y, frequency, amplitude)
    ny = 2pi * amplitude * frequency * cos.(2pi * frequency * y)
    magnitude = sqrt.(1.0 .+ ny .^ 2)
    numpts = length(y)
    nx = 1.0 ./ magnitude
    ny = ny ./ magnitude
    normals = vcat(nx', ny')
    return normals
end

function plot_levelset_normal_and_exact_normal(
    ycoords,
    levelsetnormal,
    exactnormal;
    filepath = "",
)
    fig, ax = PyPlot.subplots(2, 1, sharex = true)
    ax[1].plot(ycoords, levelsetnormal[1, :], label = "levelset")
    ax[1].plot(ycoords, exactnormal[1, :], "--", label = "exact")
    ax[1].set_title("X-component of interface normal")
    ax[1].grid()
    ax[1].legend()

    ax[2].plot(ycoords, levelsetnormal[2, :])
    ax[2].plot(ycoords, exactnormal[2, :], "--")
    ax[2].set_title("Y-component of interface normal")
    ax[2].grid()

    if length(filepath) > 0
        fig.savefig(filepath)
        return fig
    else
        return fig
    end
end

function update_cell_interface_normals_tangents_areas!(
    interfacequads,
    cellid,
    normals,
    invjac,
)
    t = CutCellDG.rotate_90(normals)
    s = CutCellDG.scale_area(t, invjac)

    idx = interfacequads.celltoquad[cellid]
    interfacequads.normals[idx] .= normals
    interfacequads.tangents[idx] .= t
    interfacequads.scaleareas[idx] .= s
end

function update_interface_normals!(interfacequads, cutmesh, exactnormalfunc)
    ncells = CutCellDG.number_of_cells(cutmesh)
    cellsign = [CutCellDG.cell_sign(cutmesh, cellid) for cellid = 1:ncells]
    cellids = findall(cellsign .== 0)
    invjac = CutCellDG.inverse_jacobian(cutmesh)

    for cellid in cellids
        cellmap = CutCellDG.cell_map(cutmesh, cellid)
        qps = cellmap(CutCellDG.points(interfacequads[1, cellid]))
        exactnormals = exactnormalfunc(qps[2, :])
        update_cell_interface_normals_tangents_areas!(
            interfacequads,
            cellid,
            exactnormals,
            invjac,
        )
    end
end

function plot_normal_stress(ycoords, parentsrr, productsrr; filepath = "")

    fig, ax = PyPlot.subplots()
    ax.plot(ycoords, parentsrr, label = "parent")
    ax.plot(ycoords, productsrr, label = "product")
    ax.grid()
    ax.set_ylabel(L"\sigma_{nn}")
    ax.set_xlabel("y")
    ax.legend()

    if length(filepath) > 0
        fig.savefig(filepath)
        return fig
    else
        return fig
    end
end

distancefunction(x) =
    perturbed_distancefunction(x, initialposition, frequency, amplitude)

initialposition = 0.5
frequency = 2.5
amplitude = 1e-3
# amplitude = 1e-10
polyorder = 2
nelmts = 17
penaltyfactor = 1e3
meshwidth = [1.0, 1.0]
numqp = required_quadrature_order(polyorder) + 2

K1, K2 = 247.0, 247.0    # GPa
mu1, mu2 = 126.0, 126.0   # GPa

lambda1 = lame_lambda(K1, mu1)
lambda2 = lame_lambda(K2, mu2)
stiffness = CutCellDG.HookeStiffness(lambda1, mu1, lambda2, mu2)
rho1 = 3.93e3           # Kg/m^3
rho2 = 3.68e3           # Kg/m^3
V01 = 1.0 / rho1
V02 = 1.0 / rho2
ΔG0Jmol = -14351.0
molarmass = 0.147
ΔG0 = ΔG0Jmol / molarmass
theta0 = -0.067

numquerypoints = 1000
querypoints = vcat(
    repeat([initialposition], numquerypoints)',
    range(0, 1, length = numquerypoints)',
)
ycoords = querypoints[2, :]

transfstress =
    CutCellDG.plane_strain_transformation_stress(lambda1, mu1, theta0)

elasticitybasis = LagrangeTensorProductBasis(2, polyorder)
levelsetbasis = HermiteTensorProductBasis(2)
quad = tensor_product_quadrature(2, 4)
dim, numpts = size(interpolation_points(levelsetbasis))
basispts = interpolation_points(elasticitybasis)

cgmesh = CutCellDG.CGMesh([0.0, 0.0], meshwidth, [nelmts, nelmts], numpts)
dgmesh = CutCellDG.DGMesh([0.0, 0.0], meshwidth, [nelmts, nelmts], basispts)

# levelset = CutCellDG.LevelSet(distancefunction, cgmesh, levelsetbasis)
levelset = CutCellDG.LevelSet(x->distancefunction(x)[1], cgmesh, levelsetbasis)

elementsize = CutCellDG.element_size(cgmesh)
minelmtsize = minimum(elementsize)
maxelmtsize = maximum(elementsize)
penalty = penaltyfactor / minelmtsize * 0.5 * (lambda1 + mu1 + lambda2 + mu2)
tol = minelmtsize^(polyorder + 1)
boundingradius = 1.5 * maxelmtsize

cutmesh = CutCellDG.CutMesh(dgmesh, levelset)
refseedpoints, seedcellids = CutCellDG.seed_zero_levelset(2, levelset, cutmesh)
spatialseedpoints =
    CutCellDG.map_to_spatial(refseedpoints, seedcellids, cutmesh)

mesh, cellquads, facequads, interfacequads =
    TES.construct_merged_mesh_and_quadratures(cutmesh, levelset, numqp)
update_interface_normals!(
    interfacequads,
    cutmesh,
    x -> perturbed_normal(x, frequency, amplitude),
)

interfacenormals = CutCellDG.collect_interface_normals(interfacequads, cutmesh)
refinterfaceqps, refcellids =
    CutCellDG.collect_interface_quadrature_points(interfacequads, +1, cutmesh)
interfaceqps = CutCellDG.map_to_spatial_on_merged_mesh(
    refinterfaceqps,
    refcellids,
    +1,
    mesh,
)
exactqpnormals = perturbed_normal(interfaceqps[2, :], frequency, amplitude)

nodaldisplacement = TES.nodal_displacement(
    mesh,
    elasticitybasis,
    cellquads,
    facequads,
    interfacequads,
    stiffness,
    theta0,
    penalty,
)

closestpoints, closestcellids = CutCellDG.closest_points_on_zero_levelset(
    querypoints,
    spatialseedpoints,
    seedcellids,
    levelset,
    tol,
    boundingradius,
)

normals = CutCellDG.collect_normals_at_spatial_points(
    closestpoints,
    closestcellids,
    levelset,
)
exactnormals = perturbed_normal(closestpoints[2, :], frequency, amplitude)

parentclosestrefpoints = CutCellDG.map_to_reference_on_merged_mesh(
    closestpoints,
    closestcellids,
    -1,
    mesh,
)
productclosestrefpoints = CutCellDG.map_to_reference_on_merged_mesh(
    closestpoints,
    closestcellids,
    +1,
    mesh,
)

################################################################################
parentstrain = CutCellDG.parent_strain(
    nodaldisplacement,
    elasticitybasis,
    parentclosestrefpoints,
    closestcellids,
    mesh,
)
parentstress = CutCellDG.parent_stress(parentstrain, stiffness)

parentradialtraction = CutCellDG.traction_force_at_points(parentstress, normals)
parentsrr = CutCellDG.traction_component(parentradialtraction, normals)

parentradialtractionalt =
    CutCellDG.traction_force_at_points(parentstress, exactnormals)
parentsrralt =
    CutCellDG.traction_component(parentradialtractionalt, exactnormals)
################################################################################
#
#
#
#
#
# ################################################################################
productstrain = CutCellDG.product_elastic_strain(
    nodaldisplacement,
    elasticitybasis,
    theta0,
    productclosestrefpoints,
    closestcellids,
    mesh,
)
productstress = CutCellDG.product_stress(productstrain, stiffness, theta0)

productradialtraction =
    CutCellDG.traction_force_at_points(productstress, normals)
productsrr = CutCellDG.traction_component(productradialtraction, normals)

productradialtractionalt =
    CutCellDG.traction_force_at_points(productstress, exactnormals)
productsrralt =
    CutCellDG.traction_component(productradialtractionalt, exactnormals)
# ################################################################################
foldername = "examples\\time-step\\perturbed-plane-interface\\analytical-normal-for-elasticity\\"

plot_normal_stress(ycoords,parentsrr,productsrr)


text:
module M2_1
function f(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    l1 = q[1]
    l2 = q[2]
    l3 = q[3]
    l1_dot = q_dot[1]
    l2_dot = q_dot[2]
    l3_dot = q_dot[3]

    2.16714051041859e-11*l1.^2 - 1.86410093239986e-11*l1.*l2 - 2.4932305151384e-11*l1.*l3 + 3.54942237399883e-14*l1 - 1.12667056568291e-11*l2.^2 + 4.12227472494603e-11*l2.*l3 - 7.5964125609303e-15*l2 - 8.05356093237194e-12*l3.^2 - 2.80736365134624e-14*l3 + 1.35847747062933e-17
end
end


text:
using SafeTestsets, Test
using Surrogates: sample, SobolSample

@safetestset "AbstractGPSurrogate" begin

using Surrogates
using SurrogatesAbstractGPs
using AbstractGPs
using Zygote

@testset "1D -> 1D" begin
    lb = 0.0
    ub = 3.0
    f = x -> log(x)*exp(x);
    x = sample(5,lb,ub,SobolSample())
    y = f.(x)
    agp1D = AbstractGPSurrogate(x,y, gp=GP(SqExponentialKernel()), Σy=0.05)
    x_new = 2.5
    y_actual = f.(x_new)
    y_predicted = agp1D([x_new])
    @test isapprox(y_predicted, y_actual, atol=0.1)
end

@testset "add points 1D" begin
    lb = 0.0
    ub = 3.0
    f = x -> x^2;
    x_points = sample(5,lb,ub,SobolSample())
    y_points = f.(x_points)
    agp1D = AbstractGPSurrogate([x_points[1]],[y_points[1]], gp=GP(SqExponentialKernel()), Σy=0.05)
    x_new = 2.5
    y_actual = f.(x_new)
    for i in 2:length(x_points)
        add_point!(agp1D, x_points[i], y_points[i]);
    end
    y_predicted = agp1D([x_new]);
    @test isapprox(y_predicted, y_actual, atol=0.1)
end

@testset "2D -> 1D" begin
    lb = [0.0; 0.0]
    ub = [2.0; 2.0]
    log_exp_f = x -> log(x[1])*exp(x[2])
    x = sample(50,lb,ub,SobolSample())
    y = log_exp_f.(x)
    agp_2D = AbstractGPSurrogate(x, y)
    x_new_2D = (2.0, 1.0)
    y_actual = log_exp_f(x_new_2D)
    y_predicted = agp_2D(x_new_2D)
    @test isapprox(y_predicted, y_actual, atol=0.1)
end

@testset "add points 2D" begin
    lb = [0.0;0.0]
    ub = [2.0;2.0]
    sphere = x -> x[1]^2 + x[2]^2
    x = sample(20, lb, ub, SobolSample())
    y = sphere.(x)
    agp_2D = AbstractGPSurrogate([x[1]],[y[1]])
    logpdf_vals = []
    push!(logpdf_vals, logpdf_surrogate(agp_2D))
    for i in 2:length(x)
        add_point!(agp_2D,x[i],y[i])
        push!(logpdf_vals, logpdf_surrogate(agp_2D))
    end
    @test first(logpdf_vals) < last(logpdf_vals) #as more points are added log marginal posterior predictive probability increases
end

@testset "check ND prediction" begin
    lb = [-1.0;-1.0;-1.0]
    ub = [1.0;1.0;1.0]
    f = x ->hypot(x...)
    x = sample(25,lb,ub,SobolSample())
    y = f.(x)
    agpND = AbstractGPSurrogate(x, y, gp = GP(SqExponentialKernel()), Σy = 0.05)
    x_new = (-0.8,0.8,0.8)
    @test agpND(x_new) ≈ f(x_new) atol=0.2
end

@testset "Optimization 1D" begin
    objective_function = x -> 2*x+1
    lb = 0.0
    ub = 6.0
    x = [2.0,4.0,6.0]
    y = [5.0,9.0,13.0]
    p = 2
    a = 2
    b = 6
    my_k_EI1 = AbstractGPSurrogate(x,y)
    surrogate_optimize(objective_function,EI(),a,b,my_k_EI1,UniformSample(),maxiters=200,num_new_samples=155)
end

@testset "Optimization ND" begin
    objective_function_ND = z -> 3*hypot(z...)+1
    x = [(1.2,3.0),(3.0,3.5),(5.2,5.7)]
    y = objective_function_ND.(x)
    theta = [2.0,2.0]
    lb = [1.0,1.0]
    ub = [6.0,6.0]
    my_k_E1N = AbstractGPSurrogate(x,y)
    surrogate_optimize(objective_function_ND,EI(),lb,ub,my_k_E1N,UniformSample())
end

@testset "check working of logpdf_surrogate 1D" begin
    lb = 0.0
    ub = 3.0
    f = x -> log(x)*exp(x);
    x = sample(5,lb,ub,SobolSample())
    y = f.(x)
    agp1D = AbstractGPSurrogate(x,y, gp=GP(SqExponentialKernel()), Σy=0.05)
    logpdf_surrogate(agp1D)
end

@testset "check working of logpdf_surrogate ND" begin
    lb = [0.0; 0.0]
    ub = [2.0; 2.0]
    f = x -> log(x[1])*exp(x[2]);
    x = sample(5,lb,ub,SobolSample())
    y = f.(x)
    agpND = AbstractGPSurrogate(x,y, gp=GP(SqExponentialKernel()), Σy=0.05)
    logpdf_surrogate(agpND)
end

lb = 0.0
ub = 3.0
n = 10
x = sample(n,lb,ub,SobolSample())
f = x -> x^2
y = f.(x)

#AbstractGP 1D
@testset "AbstractGP 1D" begin
    agp1D = AbstractGPSurrogate(x,y, gp=GP(SqExponentialKernel()), Σy=0.05)
    g = x -> agp1D'(x)
    g([2.0])
end

lb = [0.0,0.0]
ub = [10.0,10.0]
n = 5
x = sample(n,lb,ub,SobolSample())
f = x -> x[1]*x[2]
y = f.(x)

# AbstractGP ND
@testset "AbstractGPSurrogate ND" begin
    my_agp = AbstractGPSurrogate(x,y, gp=GP(SqExponentialKernel()), Σy=0.05)
    g = x ->Zygote.gradient(my_agp, x)
    #g([(2.0,5.0)])
    g((2.0,5.0))
 end

end


text:
using .JSort
import Base: read, read!, fill!, write, close

mutable struct ADCTDC
    channel::Int
    adc::Int
    tdc::Int
    ADCTDC() = new(0, 0, 0)
end

mutable struct Event
    e::Int
    de::Int
    front::Int
    back::Int
    labr::Vector{ADCTDC}
end

function Event()
    Event(0, 0, 0, 0, ADCTDC[ADCTDC() for i in 1:32])
end

abstract type Buffer end

mutable struct Buffer2D <: Buffer
    source::IO
    buffer::Matrix{UInt32}
    bp::Int
    read::Int
end
function Buffer2D(io::IO; size=2^16)
    Buffer2D(io, Matrix{UInt32}(undef, size, 2), 1, 0)
end
Buffer2D(path::AbstractString; size=2^16, rw="w") = Buffer2D(open(path, rw), size=size)

mutable struct Buffer1D <: Buffer
    source::IO
    buffer::Vector{UInt32}
    bp::Int
    read::Int
end
function Buffer1D(io::IO; size=2^16)
    Buffer1D(io, Vector{UInt32}(undef, size), 1, 0)
end
Buffer1D(path::AbstractString; size=2^16, rw="w") = Buffer1D(open(path, rw), size=size)

mutable struct Reader
    buffer::Buffer1D
    output::Matrix{Buffer2D}
    labr_size::Int
end

function Reader(parameters::Parameters; buffersize=2^15)
    output = [Buffer2D(joinpath(parameters.savepath, "edeb$(b)f$(f).bin"))
               for f in 1:8, b in 1:8]
    fsum   = [Buffer2D(joinpath(parameters.savepath, "edef$f.bin"))
              for f in 1:8]
    output = [output fsum]
    source = open(joinpath(parameters.readpath, "output.bin"), "r")
    buffer = Buffer1D(source, size=buffersize)
    Reader(buffer, output, 0)
end

function read(reader::Reader)
    event = Event()
    fillbuffer!(reader)
    read_events = 0
    while isgood(reader)
        fill!(event, reader)
        fill!(reader.output[event.front, event.back], event.e, event.de)
        # Sum over all back for a given front
        fill!(reader.output[event.front, 9], event.e, event.de)
        read_events += 1
    end
    @show read_events
    write(reader)
    close(reader)
end

function fill!(event::Event, reader::Reader)
    event.e = next!(reader)
    event.de = next!(reader)
    event.front = next!(reader)
    event.back = next!(reader)
    reader.labr_size = next!(reader)
    getbody!(reader, event.labr)
end

function fillbuffer!(reader::Reader)
    fillbuffer!(reader.buffer)
end

function getheader(reader::Reader)
    r = next!(reader.buffer, 5)
    reader.labr_size = r[5]
    r
end

function getbody!(reader::Reader, body)
    for i in 1:reader.labr_size
        body[i].channel = next!(reader)
        body[i].adc     = next!(reader)
        body[i].tdc     = next!(reader)
    end
end

function fillbuffer!(buffer::Buffer1D)
    buffer.bp = 1
    remaining = stat(buffer.source).size - buffer.read
    # Resize the buffer to the remaining space if there isn't enough
    # data left to fill it
    if remaining >= sizeof(buffer.buffer)
        read!(buffer.source, buffer.buffer)
    else
        println("Last piece")
        remaining = remaining / sizeof(eltype(buffer.buffer))|> Int
        println(remaining)
        buffer.buffer = buffer.buffer[1:remaining]
        read!(buffer.source, buffer.buffer)
    end
    buffer.read += sizeof(buffer.buffer)
end

function isgood(buffer::Buffer)::Bool
    return !eof(buffer.source) || buffer.bp < length(buffer.buffer)-1
end
isgood(r::Reader) = isgood(r.buffer)

function conditionalfill!(buffer::Buffer)
    if buffer.bp <= length(buffer.buffer)
        return
    end
    fillbuffer!(buffer)
end

function next!(buffer::Buffer)
    val = buffer.buffer[buffer.bp]
    buffer.bp += 1
    conditionalfill!(buffer)
    return val
end

function next!(buffer::Buffer, size::Int)
    [next!(buffer) for _ in 1:size]
end
next!(r::Reader) = next!(r.buffer)
next!(r::Reader, i::Int) = next!(r.buffer, i)

function fill!(buffer::Buffer, x, y)
    #println("Filling! $x, $y")
    buffer.buffer[buffer.bp, 1] = x
    buffer.buffer[buffer.bp, 2] = y
    buffer.bp += 1
    conditionalwrite(buffer)
end

function conditionalwrite(buffer::Buffer)
    if buffer.bp == size(buffer.buffer, 1)
        write(buffer)
    elseif buffer.bp > size(buffer.buffer, 1)
        throw("Corrupt buffer: index too large")
    end
end

function write(buffer::Buffer2D)
    if buffer.bp != 1
      @views write(buffer.source, transpose(buffer.buffer[1:buffer.bp, :]))
      buffer.bp = 1
    end
end

function close(buffer::Buffer)
    close(buffer.source)
end

function write(reader::Reader)
    map(write, reader.output)
end

function close(reader::Reader)
    map(close, reader.output)
    close(reader.buffer)
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
"""
Loopy belief propogation for a network.

Early stopping if change is messages < `tol` for `iters_for_convergence'
iterations. For no stopping, use tol < 0.
"""
@with_kw struct LoopyBelief <: InferenceMethod
    nsamples::Int = 500
    tol::Float64 = 1e-8
    iters_for_convergence::Int = 6
end

"""
Get the lambda-message to itself for an evidence node.
If it isn't an evidence node, this will break
"""
@inline function _evidence_lambda(nn::NodeName, evidence::Assignment, ncat::Int)
    z = zeros(ncat)
    z[evidence[nn]] = 1

    return z
end

function infer{BN<:DiscreteBayesNet}(im::LoopyBelief, inf::InferenceState{BN})

    length(inf.query) == 1 || throw(ArgumentError("There can only be one query variable"))

    nsamples, tol, iters_for_convergence = im.nsamples, im.tol, im.iters_for_convergence

    bn = inf.pgm
    nodes = names(bn)
    query = first(inf.query)
    evidence = inf.evidence
    evidence_nodes = collect(keys(evidence))

    ncat_lut = Dict(nn => ncategories(bn, nn) for nn in nodes)
    parents_lut = map(nn -> parents(bn, nn), nodes)
    children_lut = map(nn -> children(bn, nn), nodes)
    factors = map(nn -> Factor(bn, nn), nodes)

    # evidence node messages to their selves
    evidence_lambdas = map(nn -> _evidence_lambda(nn, evidence, ncat_lut[nn]), evidence_nodes)

    # the index of each node in evidence (lambda) or zero otherwise
    evidence_index = indexin(nodes, evidence_nodes)

    #=
    The messages being passed from node to node (parents or children)
        Each node has a vector containing the messages from its children
        Each message is about it, so all have the same length
    =#
    lambdas = Dict{Tuple{NodeName, NodeName}, Vector{Float64}}()

    #=
    Each node has a vector containing the messages from its parents
    Each message is about the parent, so they may have different lengths
    =#
    pis = similar(lambdas)

    #=
    Init first messages
        Instead of calculating the first πᵗ and λᵗ messages per node set all of them to 1's
        This gets rid of the edge condition for accumulating π and λ for first iteration
        I am most unsure about this...
    =#
    for (i, cpd) in enumerate(bn.cpds)

        nn = name(cpd)
        nn_ncat = ncat_lut[nn]
        nn_parents = parents_lut[i]
        nn_children = children_lut[i]

        for ch in nn_children
            # lambda from child to parents (nn)
            lambdas[(ch, nn)] = fill(1/nn_ncat, nn_ncat)

            # pi from parent (nn) to children
            #  if evidence node, set to evidence_lambda to avoid recomputing
            ev_i = evidence_index[i]
            pis[(nn, ch)] = (ev_i == 0) ? fill(1/nn_ncat, nn_ncat) :
                evidence_lambdas[ev_i]
        end
    end

    # messages are passed in parallel, so need a "new" set
    new_lambdas = similar(lambdas)
    # pi messages for evidence nodes won't be (re-)computed, so
    # just copy over *all* initial pi messages
    new_pis = deepcopy(pis)

    # begin iterating
    # use a circular buffer to keep track of the maximum change
    # in the messages per iteration
    change_index = 1
    change_per_iter = fill(Inf, iters_for_convergence)

    for iter in 1:nsamples
        max_change = -Inf

        for (i, cpd) in enumerate(bn.cpds)
            nn = name(cpd)
            nn_ncat = ncat_lut[nn]
            ϕ = factors[i]
            nn_parents = parents_lut[i]
            nn_children = children_lut[i]

            # build current lambda and pi at time t (from messages @ time t-1)

            # lambda is the product of what all its children are telling it
            #  about each instance of itself
            # nn's index in the evidence
            ev_i = evidence_index[i]
            if ev_i != 0
                # if it is evidence, then just one value will be non-zero
                # normalization will kick in and do magic ...
                lambda = evidence_lambdas[ev_i]
            else
                lambda = isempty(nn_children) ? ones(nn_ncat) :
                    reduce((x,y) -> broadcast(*, x, y), (lambdas[(ch, nn)] for ch in nn_children))
            end

            # pi is the cpd of the node weighted by the pi message that each
            #  parent sent about itself
            if isempty(nn_parents)
                pi = ϕ.potential
            else
                # multiply each probability with the pi messages from parents
                #  for that particular instantiation of the parents
                ft_temp = broadcast(*, ϕ, nn_parents,
                        [pis[(pa, nn)] for pa in nn_parents])
                pi = sum!(ft_temp, nn_parents).potential
            end

            # build lambda messages to parents
            for pa in nn_parents
                # same as above:
                # multiply each p(x|parents) by pi of all parents but one
                other_pa = setdiff(nn_parents, [pa])
                lx = isempty(other_pa) ? deepcopy(ϕ) :
                        broadcast(*, ϕ, other_pa,
                                [pis[(p, nn)] for p in other_pa])
                # sum out the other parents
                sum!(lx, other_pa)
                # weight by its current evidence lambda
                broadcast!(*, lx, nn, lambda)
                # sum out the current node
                sum!(lx, nn)
                normalize!(lx)

                new_lambdas[(nn, pa)] = lx.potential

                # find change in the message
                delta = norm(lx.potential - lambdas[(nn, pa)])
                if delta > max_change
                    max_change = delta
                end
            end

            # build pi messages to children
            if ev_i == 0
                # if an evidence node, pi will stay [0 ... 0 1 0... 0 ]
                for ch in nn_children
                    other_ch = setdiff(nn_children, [ch])
                    px = pi
                    # reduce freaks out with empty arrays
                    isempty(other_ch) ||
                        (px .*= reduce((x,y) -> broadcast(*, x, y),
                                    (lambdas[(c, nn)] for c in other_ch)))
                    normalize!(px, 1)

                    new_pis[(nn, ch)] = px

                    # find change in the message
                    delta = norm(px - pis[(nn, ch)], Inf)
                    if delta > max_change
                        max_change = delta
                    end
                end
            end
        end

        # swap 'em
        new_pis, pis = pis, new_pis
        new_lambdas, lambdas, = lambdas, new_lambdas

        change_per_iter[change_index] = max_change
        # a % b is in [0, b), so add 1
        change_index = change_index % iters_for_convergence + 1

        if maximum(change_per_iter) <= tol
            break
        end
    end

    # compute belief P(x|e)
    qi = findfirst(nodes, query)
    # lambda and pi one last time
    ϕ = factors[qi]
    nn_ncat = ncat_lut[query]
    nn_children = children_lut[qi]
    nn_parents = parents_lut[qi]

    lambda = isempty(nn_children) ? ones(nn_ncat) :
            reduce((x,y) -> broadcast(*, x, y), (lambdas[(ch, query)] for ch in nn_children))

    if isempty(nn_parents)
        pi = ϕ.potential
    else
        broadcast!(*, ϕ, nn_parents,
                [pis[(pa, query)] for pa in nn_parents])
        pi = sum!(ϕ, nn_parents).potential
    end

    ftr = Factor([query], broadcast(*, lambda, pi))
    normalize!(ftr)

    return ftr
end


text:
# ------------ Switch trace ------------ #

struct SwitchTrace{A <: Tuple, T, U <: Trace} <: Trace
    gen_fn::GenerativeFunction{T, U}
    branch::U
    retval::T
    args::A
    score::Float64
    noise::Float64
    function SwitchTrace(gen_fn::GenerativeFunction{T, U},
            branch::U,
            retval::T, args::A, score::Float64,
            noise::Float64) where {T, A, U}
        new{A, T, U}(gen_fn, branch, retval, args, score, noise)
    end
end

@inline get_choices(tr::SwitchTrace) = get_choices(tr.branch)
@inline get_retval(tr::SwitchTrace) = tr.retval
@inline get_args(tr::SwitchTrace) = tr.args
@inline get_score(tr::SwitchTrace) = tr.score
@inline get_gen_fn(tr::SwitchTrace) = tr.gen_fn
@inline Base.getindex(tr::SwitchTrace, addr) = Base.getindex(tr.branch, addr)
@inline project(tr::SwitchTrace, selection::Selection) = project(tr.branch, selection)
@inline project(tr::SwitchTrace, ::EmptySelection) = tr.noise


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

module UMFPACK

export UmfpackLU

import Base: (\), getproperty, show, size
using LinearAlgebra
import LinearAlgebra: Factorization, det, lu, ldiv!

using SparseArrays
import SparseArrays: nnz

import ..increment, ..increment!, ..decrement, ..decrement!

include("umfpack_h.jl")
struct MatrixIllConditionedException <: Exception
    msg::AbstractString
end

function umferror(status::Integer)
    if status==UMFPACK_OK
        return
    elseif status==UMFPACK_WARNING_singular_matrix
        throw(LinearAlgebra.SingularException(0))
    elseif status==UMFPACK_WARNING_determinant_underflow
        throw(MatrixIllConditionedException("the determinant is nonzero but underflowed"))
    elseif status==UMFPACK_WARNING_determinant_overflow
        throw(MatrixIllConditionedException("the determinant overflowed"))
    elseif status==UMFPACK_ERROR_out_of_memory
        throw(OutOfMemoryError())
    elseif status==UMFPACK_ERROR_invalid_Numeric_object
        throw(ArgumentError("invalid UMFPack numeric object"))
    elseif status==UMFPACK_ERROR_invalid_Symbolic_object
        throw(ArgumentError("invalid UMFPack symbolic object"))
    elseif status==UMFPACK_ERROR_argument_missing
        throw(ArgumentError("a required argument to UMFPack is missing"))
    elseif status==UMFPACK_ERROR_n_nonpositive
        throw(ArgumentError("the number of rows or columns of the matrix must be greater than zero"))
    elseif status==UMFPACK_ERROR_invalid_matrix
        throw(ArgumentError("invalid matrix"))
    elseif status==UMFPACK_ERROR_different_pattern
        throw(ArgumentError("pattern of the matrix changed"))
    elseif status==UMFPACK_ERROR_invalid_system
        throw(ArgumentError("invalid sys argument provided to UMFPack solver"))
    elseif status==UMFPACK_ERROR_invalid_permutation
        throw(ArgumentError("invalid permutation"))
    elseif status==UMFPACK_ERROR_file_IO
        throw(ErrorException("error saving / loading UMFPack decomposition"))
    elseif status==UMFPACK_ERROR_ordering_failed
        throw(ErrorException("the ordering method failed"))
    elseif status==UMFPACK_ERROR_internal_error
        throw(ErrorException("an internal error has occurred, of unknown cause"))
    else
        throw(ErrorException("unknown UMFPack error code: $status"))
    end
end

macro isok(A)
    :(umferror($(esc(A))))
end

# check the size of SuiteSparse_long
if Int(ccall((:jl_cholmod_sizeof_long,:libsuitesparse_wrapper),Csize_t,())) == 4
    const UmfpackIndexTypes = (:Int32,)
    const UMFITypes = Int32
else
    const UmfpackIndexTypes = (:Int32, :Int64)
    const UMFITypes = Union{Int32, Int64}
end

const UMFVTypes = Union{Float64,ComplexF64}

## UMFPACK

# the control and info arrays
const umf_ctrl = Vector{Float64}(undef, UMFPACK_CONTROL)
ccall((:umfpack_dl_defaults,:libumfpack), Cvoid, (Ptr{Float64},), umf_ctrl)
const umf_info = Vector{Float64}(undef, UMFPACK_INFO)

function show_umf_ctrl(level::Real = 2.0)
    old_prt::Float64 = umf_ctrl[1]
    umf_ctrl[1] = Float64(level)
    ccall((:umfpack_dl_report_control, :libumfpack), Cvoid, (Ptr{Float64},), umf_ctrl)
    umf_ctrl[1] = old_prt
end

function show_umf_info(level::Real = 2.0)
    old_prt::Float64 = umf_ctrl[1]
    umf_ctrl[1] = Float64(level)
    ccall((:umfpack_dl_report_info, :libumfpack), Cvoid,
          (Ptr{Float64}, Ptr{Float64}), umf_ctrl, umf_info)
    umf_ctrl[1] = old_prt
end

## Should this type be immutable?
mutable struct UmfpackLU{Tv<:UMFVTypes,Ti<:UMFITypes} <: Factorization{Tv}
    symbolic::Ptr{Cvoid}
    numeric::Ptr{Cvoid}
    m::Int
    n::Int
    colptr::Vector{Ti}                  # 0-based column pointers
    rowval::Vector{Ti}                  # 0-based row indices
    nzval::Vector{Tv}
    status::Int
end

Base.adjoint(F::UmfpackLU) = Adjoint(F)
Base.transpose(F::UmfpackLU) = Transpose(F)

"""
    lu(A::SparseMatrixCSC; check = true) -> F::UmfpackLU

Compute the LU factorization of a sparse matrix `A`.

For sparse `A` with real or complex element type, the return type of `F` is
`UmfpackLU{Tv, Ti}`, with `Tv` = [`Float64`](@ref) or `ComplexF64` respectively and
`Ti` is an integer type ([`Int32`](@ref) or [`Int64`](@ref)).

When `check = true`, an error is thrown if the decomposition fails.
When `check = false`, responsibility for checking the decomposition's
validity (via [`issuccess`](@ref)) lies with the user.

The individual components of the factorization `F` can be accessed by indexing:

| Component | Description                         |
|:----------|:------------------------------------|
| `L`       | `L` (lower triangular) part of `LU` |
| `U`       | `U` (upper triangular) part of `LU` |
| `p`       | right permutation `Vector`          |
| `q`       | left permutation `Vector`           |
| `Rs`      | `Vector` of scaling factors         |
| `:`       | `(L,U,p,q,Rs)` components           |

The relation between `F` and `A` is

`F.L*F.U == (F.Rs .* A)[F.p, F.q]`

`F` further supports the following functions:

- [`\\`](@ref)
- [`cond`](@ref)
- [`det`](@ref)

!!! note
    `lu(A::SparseMatrixCSC)` uses the UMFPACK library that is part of
    SuiteSparse. As this library only supports sparse matrices with [`Float64`](@ref) or
    `ComplexF64` elements, `lu` converts `A` into a copy that is of type
    `SparseMatrixCSC{Float64}` or `SparseMatrixCSC{ComplexF64}` as appropriate.
"""
function lu(S::SparseMatrixCSC{<:UMFVTypes,<:UMFITypes}; check::Bool = true)
    zerobased = S.colptr[1] == 0
    res = UmfpackLU(C_NULL, C_NULL, S.m, S.n,
                    zerobased ? copy(S.colptr) : decrement(S.colptr),
                    zerobased ? copy(S.rowval) : decrement(S.rowval),
                    copy(S.nzval), 0)
    finalizer(umfpack_free_symbolic, res)
    umfpack_numeric!(res)
    check && (issuccess(res) || throw(LinearAlgebra.SingularException(0)))
    return res
end
lu(A::SparseMatrixCSC{<:Union{Float16,Float32},Ti};
   check::Bool = true) where {Ti<:UMFITypes} =
    lu(convert(SparseMatrixCSC{Float64,Ti}, A); check = check)
lu(A::SparseMatrixCSC{<:Union{ComplexF16,ComplexF32},Ti};
   check::Bool = true) where {Ti<:UMFITypes} =
    lu(convert(SparseMatrixCSC{ComplexF64,Ti}, A); check = check)
lu(A::Union{SparseMatrixCSC{T},SparseMatrixCSC{Complex{T}}};
   check::Bool = true) where {T<:AbstractFloat} =
    throw(ArgumentError(string("matrix type ", typeof(A), "not supported. ",
    "Try lu(convert(SparseMatrixCSC{Float64/ComplexF64,Int}, A)) for ",
    "sparse floating point LU using UMFPACK or lu(Array(A)) for generic ",
    "dense LU.")))
lu(A::SparseMatrixCSC; check::Bool = true) = lu(float(A); check = check)

size(F::UmfpackLU) = (F.m, F.n)
function size(F::UmfpackLU, dim::Integer)
    if dim < 1
        throw(ArgumentError("size: dimension $dim out of range"))
    elseif dim == 1
        return Int(F.m)
    elseif dim == 2
        return Int(F.n)
    else
        return 1
    end
end

function show(io::IO, F::UmfpackLU)
    print(io, "UMFPACK LU Factorization of a $(size(F)) sparse matrix")
    F.numeric != C_NULL && print(io, '\n', F.numeric)
end

## Wrappers for UMFPACK functions

# generate the name of the C function according to the value and integer types
umf_nm(nm,Tv,Ti) = "umfpack_" * (Tv == :Float64 ? "d" : "z") * (Ti == :Int64 ? "l_" : "i_") * nm

for itype in UmfpackIndexTypes
    sym_r = umf_nm("symbolic", :Float64, itype)
    sym_c = umf_nm("symbolic", :ComplexF64, itype)
    num_r = umf_nm("numeric", :Float64, itype)
    num_c = umf_nm("numeric", :ComplexF64, itype)
    sol_r = umf_nm("solve", :Float64, itype)
    sol_c = umf_nm("solve", :ComplexF64, itype)
    det_r = umf_nm("get_determinant", :Float64, itype)
    det_z = umf_nm("get_determinant", :ComplexF64, itype)
    lunz_r = umf_nm("get_lunz", :Float64, itype)
    lunz_z = umf_nm("get_lunz", :ComplexF64, itype)
    get_num_r = umf_nm("get_numeric", :Float64, itype)
    get_num_z = umf_nm("get_numeric", :ComplexF64, itype)
    @eval begin
        function umfpack_symbolic!(U::UmfpackLU{Float64,$itype})
            if U.symbolic != C_NULL return U end
            tmp = Vector{Ptr{Cvoid}}(undef, 1)
            @isok ccall(($sym_r, :libumfpack), $itype,
                        ($itype, $itype, Ptr{$itype}, Ptr{$itype}, Ptr{Float64}, Ptr{Cvoid},
                         Ptr{Float64}, Ptr{Float64}),
                        U.m, U.n, U.colptr, U.rowval, U.nzval, tmp,
                        umf_ctrl, umf_info)
            U.symbolic = tmp[1]
            return U
        end
        function umfpack_symbolic!(U::UmfpackLU{ComplexF64,$itype})
            if U.symbolic != C_NULL return U end
            tmp = Vector{Ptr{Cvoid}}(undef, 1)
            @isok ccall(($sym_c, :libumfpack), $itype,
                        ($itype, $itype, Ptr{$itype}, Ptr{$itype}, Ptr{Float64}, Ptr{Float64}, Ptr{Cvoid},
                         Ptr{Float64}, Ptr{Float64}),
                        U.m, U.n, U.colptr, U.rowval, real(U.nzval), imag(U.nzval), tmp,
                        umf_ctrl, umf_info)
            U.symbolic = tmp[1]
            return U
        end
        function umfpack_numeric!(U::UmfpackLU{Float64,$itype})
            if U.numeric != C_NULL return U end
            if U.symbolic == C_NULL umfpack_symbolic!(U) end
            tmp = Vector{Ptr{Cvoid}}(undef, 1)
            status = ccall(($num_r, :libumfpack), $itype,
                           (Ptr{$itype}, Ptr{$itype}, Ptr{Float64}, Ptr{Cvoid}, Ptr{Cvoid},
                            Ptr{Float64}, Ptr{Float64}),
                           U.colptr, U.rowval, U.nzval, U.symbolic, tmp,
                           umf_ctrl, umf_info)
            if status != UMFPACK_WARNING_singular_matrix
                umferror(status)
            end
            U.numeric = tmp[1]
            U.status = status
            return U
        end
        function umfpack_numeric!(U::UmfpackLU{ComplexF64,$itype})
            if U.numeric != C_NULL return U end
            if U.symbolic == C_NULL umfpack_symbolic!(U) end
            tmp = Vector{Ptr{Cvoid}}(undef, 1)
            status = ccall(($num_c, :libumfpack), $itype,
                           (Ptr{$itype}, Ptr{$itype}, Ptr{Float64}, Ptr{Float64}, Ptr{Cvoid}, Ptr{Cvoid},
                            Ptr{Float64}, Ptr{Float64}),
                           U.colptr, U.rowval, real(U.nzval), imag(U.nzval), U.symbolic, tmp,
                           umf_ctrl, umf_info)
            if status != UMFPACK_WARNING_singular_matrix
                umferror(status)
            end
            U.numeric = tmp[1]
            U.status = status
            return U
        end
        function solve!(x::StridedVector{Float64}, lu::UmfpackLU{Float64,$itype}, b::StridedVector{Float64}, typ::Integer)
            if x === b
                throw(ArgumentError("output array must not be aliased with input array"))
            end
            if stride(x, 1) != 1 || stride(b, 1) != 1
                throw(ArgumentError("in and output vectors must have unit strides"))
            end
            umfpack_numeric!(lu)
            (size(b,1) == lu.m) && (size(b) == size(x)) || throw(DimensionMismatch())
            @isok ccall(($sol_r, :libumfpack), $itype,
                ($itype, Ptr{$itype}, Ptr{$itype}, Ptr{Float64},
                 Ptr{Float64}, Ptr{Float64}, Ptr{Cvoid}, Ptr{Float64},
                 Ptr{Float64}),
                typ, lu.colptr, lu.rowval, lu.nzval,
                x, b, lu.numeric, umf_ctrl,
                umf_info)
            return x
        end
        function solve!(x::StridedVector{ComplexF64}, lu::UmfpackLU{ComplexF64,$itype}, b::StridedVector{ComplexF64}, typ::Integer)
            if x === b
                throw(ArgumentError("output array must not be aliased with input array"))
            end
            if stride(x, 1) != 1 || stride(b, 1) != 1
                throw(ArgumentError("in and output vectors must have unit strides"))
            end
            umfpack_numeric!(lu)
            (size(b, 1) == lu.m) && (size(b) == size(x)) || throw(DimensionMismatch())
            n = size(b, 1)
            @isok ccall(($sol_c, :libumfpack), $itype,
                        ($itype, Ptr{$itype}, Ptr{$itype}, Ptr{Float64},
                         Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
                         Ptr{Float64}, Ptr{Cvoid}, Ptr{Float64}, Ptr{Float64}),
                        typ, lu.colptr, lu.rowval, lu.nzval,
                        C_NULL, x, C_NULL, b,
                        C_NULL, lu.numeric, umf_ctrl, umf_info)
            return x
        end
        function det(lu::UmfpackLU{Float64,$itype})
            mx = Ref{Float64}()
            @isok ccall(($det_r,:libumfpack), $itype,
                           (Ptr{Float64},Ptr{Float64},Ptr{Cvoid},Ptr{Float64}),
                           mx, C_NULL, lu.numeric, umf_info)
            mx[]
        end
        function det(lu::UmfpackLU{ComplexF64,$itype})
            mx = Ref{Float64}()
            mz = Ref{Float64}()
            @isok ccall(($det_z,:libumfpack), $itype,
                        (Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Cvoid},Ptr{Float64}),
                        mx, mz, C_NULL, lu.numeric, umf_info)
            complex(mx[], mz[])
        end
        function umf_lunz(lu::UmfpackLU{Float64,$itype})
            lnz = Ref{$itype}()
            unz = Ref{$itype}()
            n_row = Ref{$itype}()
            n_col = Ref{$itype}()
            nz_diag = Ref{$itype}()
            @isok ccall(($lunz_r,:libumfpack), $itype,
                           (Ptr{$itype},Ptr{$itype},Ptr{$itype},Ptr{$itype},Ptr{$itype},Ptr{Cvoid}),
                           lnz, unz, n_row, n_col, nz_diag, lu.numeric)
            (lnz[], unz[], n_row[], n_col[], nz_diag[])
        end
        function umf_lunz(lu::UmfpackLU{ComplexF64,$itype})
            lnz = Ref{$itype}()
            unz = Ref{$itype}()
            n_row = Ref{$itype}()
            n_col = Ref{$itype}()
            nz_diag = Ref{$itype}()
            @isok ccall(($lunz_z,:libumfpack), $itype,
                           (Ptr{$itype},Ptr{$itype},Ptr{$itype},Ptr{$itype},Ptr{$itype},Ptr{Cvoid}),
                           lnz, unz, n_row, n_col, nz_diag, lu.numeric)
            (lnz[], unz[], n_row[], n_col[], nz_diag[])
        end
        function umf_extract(lu::UmfpackLU{Float64,$itype})
            umfpack_numeric!(lu)        # ensure the numeric decomposition exists
            (lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
            Lp = Vector{$itype}(undef, n_row + 1)
            Lj = Vector{$itype}(undef, lnz) # L is returned in CSR (compressed sparse row) format
            Lx = Vector{Float64}(undef, lnz)
            Up = Vector{$itype}(undef, n_col + 1)
            Ui = Vector{$itype}(undef, unz)
            Ux = Vector{Float64}(undef, unz)
            P  = Vector{$itype}(undef, n_row)
            Q  = Vector{$itype}(undef, n_col)
            Rs = Vector{Float64}(undef, n_row)
            @isok ccall(($get_num_r,:libumfpack), $itype,
                        (Ptr{$itype},Ptr{$itype},Ptr{Float64},
                         Ptr{$itype},Ptr{$itype},Ptr{Float64},
                         Ptr{$itype},Ptr{$itype},Ptr{Cvoid},
                         Ref{$itype},Ptr{Float64},Ptr{Cvoid}),
                        Lp,Lj,Lx,
                        Up,Ui,Ux,
                        P, Q, C_NULL,
                        0, Rs, lu.numeric)
            (copy(transpose(SparseMatrixCSC(min(n_row, n_col), n_row, increment!(Lp), increment!(Lj), Lx))),
             SparseMatrixCSC(min(n_row, n_col), n_col, increment!(Up), increment!(Ui), Ux),
             increment!(P), increment!(Q), Rs)
        end
        function umf_extract(lu::UmfpackLU{ComplexF64,$itype})
            umfpack_numeric!(lu)        # ensure the numeric decomposition exists
            (lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
            Lp = Vector{$itype}(undef, n_row + 1)
            Lj = Vector{$itype}(undef, lnz) # L is returned in CSR (compressed sparse row) format
            Lx = Vector{Float64}(undef, lnz)
            Lz = Vector{Float64}(undef, lnz)
            Up = Vector{$itype}(undef, n_col + 1)
            Ui = Vector{$itype}(undef, unz)
            Ux = Vector{Float64}(undef, unz)
            Uz = Vector{Float64}(undef, unz)
            P  = Vector{$itype}(undef, n_row)
            Q  = Vector{$itype}(undef, n_col)
            Rs = Vector{Float64}(undef, n_row)
            @isok ccall(($get_num_z,:libumfpack), $itype,
                        (Ptr{$itype},Ptr{$itype},Ptr{Float64},Ptr{Float64},
                         Ptr{$itype},Ptr{$itype},Ptr{Float64},Ptr{Float64},
                         Ptr{$itype},Ptr{$itype},Ptr{Cvoid}, Ptr{Cvoid},
                         Ref{$itype},Ptr{Float64},Ptr{Cvoid}),
                        Lp,Lj,Lx,Lz,
                        Up,Ui,Ux,Uz,
                        P, Q, C_NULL, C_NULL,
                        0, Rs, lu.numeric)
            (copy(transpose(SparseMatrixCSC(min(n_row, n_col), n_row, increment!(Lp), increment!(Lj), complex.(Lx, Lz)))),
             SparseMatrixCSC(min(n_row, n_col), n_col, increment!(Up), increment!(Ui), complex.(Ux, Uz)),
             increment!(P), increment!(Q), Rs)
        end
    end
end

function nnz(lu::UmfpackLU)
    lnz, unz, = umf_lunz(lu)
    return Int(lnz + unz)
end

LinearAlgebra.issuccess(lu::UmfpackLU) = lu.status == UMFPACK_OK

### Solve with Factorization

import LinearAlgebra.ldiv!

ldiv!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} =
    ldiv!(B, lu, copy(B))
ldiv!(translu::Transpose{T,<:UmfpackLU{T}}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} =
    (lu = translu.parent; ldiv!(B, transpose(lu), copy(B)))
ldiv!(adjlu::Adjoint{T,<:UmfpackLU{T}}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} =
    (lu = adjlu.parent; ldiv!(B, adjoint(lu), copy(B)))
ldiv!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) =
    ldiv!(B, lu, copy(B))
ldiv!(translu::Transpose{Float64,<:UmfpackLU{Float64}}, B::StridedVecOrMat{<:Complex}) =
    (lu = translu.parent; ldiv!(B, transpose(lu), copy(B)))
ldiv!(adjlu::Adjoint{Float64,<:UmfpackLU{Float64}}, B::StridedVecOrMat{<:Complex}) =
    (lu = adjlu.parent; ldiv!(B, adjoint(lu), copy(B)))

ldiv!(X::StridedVecOrMat{T}, lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} =
    _Aq_ldiv_B!(X, lu, B, UMFPACK_A)
ldiv!(X::StridedVecOrMat{T}, translu::Transpose{T,<:UmfpackLU{T}}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} =
    (lu = translu.parent; _Aq_ldiv_B!(X, lu, B, UMFPACK_Aat))
ldiv!(X::StridedVecOrMat{T}, adjlu::Adjoint{T,<:UmfpackLU{T}}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} =
    (lu = adjlu.parent; _Aq_ldiv_B!(X, lu, B, UMFPACK_At))
ldiv!(X::StridedVecOrMat{Tb}, lu::UmfpackLU{Float64}, B::StridedVecOrMat{Tb}) where {Tb<:Complex} =
    _Aq_ldiv_B!(X, lu, B, UMFPACK_A)
ldiv!(X::StridedVecOrMat{Tb}, translu::Transpose{Float64,<:UmfpackLU{Float64}}, B::StridedVecOrMat{Tb}) where {Tb<:Complex} =
    (lu = translu.parent; _Aq_ldiv_B!(X, lu, B, UMFPACK_Aat))
ldiv!(X::StridedVecOrMat{Tb}, adjlu::Adjoint{Float64,<:UmfpackLU{Float64}}, B::StridedVecOrMat{Tb}) where {Tb<:Complex} =
    (lu = adjlu.parent; _Aq_ldiv_B!(X, lu, B, UMFPACK_At))

function _Aq_ldiv_B!(X::StridedVecOrMat, lu::UmfpackLU, B::StridedVecOrMat, transposeoptype)
    if size(X, 2) != size(B, 2)
        throw(DimensionMismatch("input and output arrays must have same number of columns"))
    end
    _AqldivB_kernel!(X, lu, B, transposeoptype)
    return X
end
function _AqldivB_kernel!(x::StridedVector{T}, lu::UmfpackLU{T},
                          b::StridedVector{T}, transposeoptype) where T<:UMFVTypes
    solve!(x, lu, b, transposeoptype)
end
function _AqldivB_kernel!(X::StridedMatrix{T}, lu::UmfpackLU{T},
                          B::StridedMatrix{T}, transposeoptype) where T<:UMFVTypes
    for col in 1:size(X, 2)
        solve!(view(X, :, col), lu, view(B, :, col), transposeoptype)
    end
end
function _AqldivB_kernel!(x::StridedVector{Tb}, lu::UmfpackLU{Float64},
                          b::StridedVector{Tb}, transposeoptype) where Tb<:Complex
    r, i = similar(b, Float64), similar(b, Float64)
    solve!(r, lu, Vector{Float64}(real(b)), transposeoptype)
    solve!(i, lu, Vector{Float64}(imag(b)), transposeoptype)
    map!(complex, x, r, i)
end
function _AqldivB_kernel!(X::StridedMatrix{Tb}, lu::UmfpackLU{Float64},
                          B::StridedMatrix{Tb}, transposeoptype) where Tb<:Complex
    r = similar(B, Float64, size(B, 1))
    i = similar(B, Float64, size(B, 1))
    for j in 1:size(B, 2)
        solve!(r, lu, Vector{Float64}(real(view(B, :, j))), transposeoptype)
        solve!(i, lu, Vector{Float64}(imag(view(B, :, j))), transposeoptype)
        map!(complex, view(X, :, j), r, i)
    end
end

@inline function getproperty(lu::UmfpackLU, d::Symbol)
    if d == :L || d == :U || d == :p || d == :q || d == :Rs || d == :(:)
        # Guard the call to umf_extract behaind a branch to avoid infinite recursion
        L, U, p, q, Rs = umf_extract(lu)
        if d == :L
            return L
        elseif d == :U
            return U
        elseif d == :p
            return p
        elseif d == :q
            return q
        elseif d == :Rs
            return Rs
        elseif d == :(:)
            return (L, U, p, q, Rs)
        end
    else
        getfield(lu, d)
    end
end

for Tv in (:Float64, :ComplexF64), Ti in UmfpackIndexTypes
    f = Symbol(umf_nm("free_symbolic", Tv, Ti))
    @eval begin
        function ($f)(symb::Ptr{Cvoid})
            tmp = [symb]
            ccall(($(string(f)), :libumfpack), Cvoid, (Ptr{Cvoid},), tmp)
        end

        function umfpack_free_symbolic(lu::UmfpackLU{$Tv,$Ti})
            if lu.symbolic == C_NULL return lu end
            umfpack_free_numeric(lu)
            ($f)(lu.symbolic)
            lu.symbolic = C_NULL
            return lu
        end
    end

    f = Symbol(umf_nm("free_numeric", Tv, Ti))
    @eval begin
        function ($f)(num::Ptr{Cvoid})
            tmp = [num]
            ccall(($(string(f)), :libumfpack), Cvoid, (Ptr{Cvoid},), tmp)
        end
        function umfpack_free_numeric(lu::UmfpackLU{$Tv,$Ti})
            if lu.numeric == C_NULL return lu end
            ($f)(lu.numeric)
            lu.numeric = C_NULL
            return lu
        end
    end
end

function umfpack_report_symbolic(symb::Ptr{Cvoid}, level::Real)
    old_prl::Float64 = umf_ctrl[UMFPACK_PRL]
    umf_ctrl[UMFPACK_PRL] = Float64(level)
    @isok ccall((:umfpack_dl_report_symbolic, :libumfpack), Int,
                (Ptr{Cvoid}, Ptr{Float64}), symb, umf_ctrl)
    umf_ctrl[UMFPACK_PRL] = old_prl
end

umfpack_report_symbolic(symb::Ptr{Cvoid}) = umfpack_report_symbolic(symb, 4.)

function umfpack_report_symbolic(lu::UmfpackLU, level::Real)
    umfpack_report_symbolic(umfpack_symbolic!(lu).symbolic, level)
end

umfpack_report_symbolic(lu::UmfpackLU) = umfpack_report_symbolic(lu.symbolic,4.)
function umfpack_report_numeric(num::Ptr{Cvoid}, level::Real)
    old_prl::Float64 = umf_ctrl[UMFPACK_PRL]
    umf_ctrl[UMFPACK_PRL] = Float64(level)
    @isok ccall((:umfpack_dl_report_numeric, :libumfpack), Int,
                (Ptr{Cvoid}, Ptr{Float64}), num, umf_ctrl)
    umf_ctrl[UMFPACK_PRL] = old_prl
end

umfpack_report_numeric(num::Ptr{Cvoid}) = umfpack_report_numeric(num, 4.)
function umfpack_report_numeric(lu::UmfpackLU, level::Real)
    umfpack_report_numeric(umfpack_numeric!(lu).numeric, level)
end

umfpack_report_numeric(lu::UmfpackLU) = umfpack_report_numeric(lu,4.)

end # UMFPACK module

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
"""
	$(SIGNATURES)

Macro to set common fields in structs. See
https://discourse.julialang.org/t/julia-learning-macros-metaprogramming/45753/3

# Example
```
@common_fields set1 begin
    x :: Int
    y :: Float64
end

struct Foo
    @set1
    z
end
```
"""
macro common_fields(name, definition)
    return quote
        macro $(esc(name))()
            esc($(Expr(:quote, definition)))
        end
    end
end

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


text:
type OpusArray
    encoded_stream::IO
    function OpusArray(raw_audio, fs=48000)
        buf = IOBuffer()
        Opus.save(buf, Float32[Float32(x) for x in raw_audio], fs)
        return new(buf)
    end
end

function writemime(io::IO, ::MIME"text/html", x::OpusArray)
    data = base64encode(bytestring(x.encoded_stream))
    markup = """<audio controls="controls" {autoplay}>
                <source src="data:audio/ogg;base64,$data" type="audio/ogg" />
                Your browser does not support the audio element.
                </audio>"""
    print(io, markup)
end


text:
@testset "cartpole_env" begin

    env = CartPoleEnv(; rng = MersenneTwister(123))
    RLBase.test_interfaces!(env)
    RLBase.test_runnable!(env)

    env = CartPoleEnv(;T=Float32, rng = MersenneTwister(123), thetathreshold = 90.0)
    RLBase.test_interfaces!(env)
    RLBase.test_runnable!(env)

end


text:
using LinearAlgebra: I, qr

push!(LOAD_PATH, @__DIR__)
using MultiprecisionFloats
using RKTK2
using DZMisc

setprecision(5000)
const AccurateReal = BigFloat

function log(args...)
    println(args...)
    flush(stdout)
end

function constraint_residual(x)
    [x[1]^2 + x[2]^2 + x[3]^2 - one(AccurateReal), x[1] + x[2]^2]
end

function constraint_jacobian(x)
    [dbl(x[1]) dbl(x[2]) dbl(x[3]);
     one(AccurateReal) dbl(x[2]) zero(AccurateReal)]
end

function constraint_objective(x)
    norm2(constraint_residual(x))
end

function objective(x)
    x[1] * x[3]
end

function gradient(x)
    [x[3], zero(AccurateReal), x[1]]
end

function constrained_gradient(x)
    constraint_jac = copy(transpose(constraint_jacobian(x)))
    orthonormalize_columns!(constraint_jac)
    grad = gradient(x)
    grad - constraint_jac * (transpose(constraint_jac) * grad)
end

function constrain_step(x, step)
    constraint_jac = copy(transpose(constraint_jacobian(x)))
    orthonormalize_columns!(constraint_jac)
    step - constraint_jac * (transpose(constraint_jac) * step)
end

################################################################################

function constrain(x)
    x_old, obj_old = x, constraint_objective(x)
    while true
        direction = qr(constraint_jacobian(x_old)) \ constraint_residual(x_old)
        x_new = x_old - direction
        obj_new = constraint_objective(x_new)
        if obj_new < obj_old
            x_old, obj_old = x_new, obj_new
        else
            break
        end
    end
    x_old, obj_old
end

function constrained_step_value(step_size, x, step_direction, step_norm)
    x_new, obj_new = constrain(x - (step_size / step_norm) * step_direction)
    if obj_new < AccurateReal(1e-75)
        objective(x_new)
    else
        AccurateReal(NaN)
    end
end

struct ConstrainedBFGSOptimizer
    objective_value::Ref{AccurateReal}
    last_step_size::Ref{AccurateReal}
end

const INPUT_POINT = AccurateReal[-0.5, 0.707107, -0.5]
const NUM_VARS = 3

const opt = ConstrainedBFGSOptimizer(
    objective(INPUT_POINT),
    AccurateReal(0.00001))

function main()

    x = copy(INPUT_POINT)
    inv_hess = Matrix{AccurateReal}(I, NUM_VARS, NUM_VARS)

    cons_grad = constrained_gradient(x)
    cons_grad_norm = norm(cons_grad)
    while true

        # log("Performing gradient descent step...")
        grad_step_size, obj_grad = quadratic_line_search(constrained_step_value,
            opt.objective_value[], opt.last_step_size[], x, cons_grad, cons_grad_norm)

        # log("Performing BFGS step...")
        bfgs_step = constrain_step(x, inv_hess * cons_grad)
        bfgs_step_norm = norm(bfgs_step)
        bfgs_step_size, obj_bfgs = quadratic_line_search(constrained_step_value,
            opt.objective_value[], opt.last_step_size[], x, bfgs_step, bfgs_step_norm)

        if obj_bfgs < opt.objective_value[] && obj_bfgs <= obj_grad
            x, _ = constrain(x - (bfgs_step_size / bfgs_step_norm) * bfgs_step)
            cons_grad_new = constrained_gradient(x)
            cons_grad_norm_new = norm(cons_grad_new)
            update_inverse_hessian!(inv_hess,
                -bfgs_step_size / bfgs_step_norm,
                bfgs_step,
                cons_grad_new - cons_grad,
                Vector{AccurateReal}(undef, NUM_VARS))
                opt.last_step_size[] = bfgs_step_size
                opt.objective_value[] = obj_bfgs
            cons_grad = cons_grad_new
            cons_grad_norm = cons_grad_norm_new
            log(Float64(opt.objective_value[]), " | ",
                Float64(cons_grad_norm), " | ",
                Float64(opt.last_step_size[]), " | BFGS")
        elseif obj_grad < opt.objective_value[]
            x, _ = constrain(x - (grad_step_size / cons_grad_norm) * cons_grad)
            inv_hess = Matrix{AccurateReal}(I, NUM_VARS, NUM_VARS)
            opt.last_step_size[] = grad_step_size
            opt.objective_value[] = obj_grad
            cons_grad = constrained_gradient(x)
            cons_grad_norm = norm(cons_grad)
            log(Float64(opt.objective_value[]), " | ",
                Float64(cons_grad_norm), " | ",
                Float64(opt.last_step_size[]), " | GRAD")
        else
            log("Neither BFGS step (", Float64(obj_bfgs),
                ") or gradient step (", Float64(obj_grad),
                ") are better than previous point (", Float64(opt.objective_value[]), ")")
            break
        end
    end
end

main()


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

using Markdown
using InteractiveUtils

# ╔═╡ 1d2793ad-34cb-4704-b139-f55e9af06317
using Plots

# ╔═╡ 3957e56c-2d81-4cb8-9825-963ae0dc49e4
using Dates

# ╔═╡ 1196b181-459a-47e8-b57e-f320b0471a89
begin
	dates = ["4/6/2021","4/24/2021","4/25/2021","4/26/2021","4/27/2021",
		"4/28/2021","4/30/2021","5/1/2021","5/7/2021","5/8/2021","5/22/2021",
		"5/23/2021"]
	days = [1,18,19,20,21,22,24,25,31,32,46,47]
	seconds = [44.1,22,7.83,2.12,2.04,1.11,0.938,0.831,0.631,0.535,0.366,
		0.217]
end

# ╔═╡ 9b3d18a3-a176-4402-874c-45f15d1f6e14
plot(days,seconds,xlabel="Days",ylabel="Seconds",title="Benchmark speed for Ethane/6-31G Integrals",label=false)

# ╔═╡ 59ebedbe-b0c9-439b-bec6-a0aa094e6c84
md"Here's how to do this using dates instead of days. You need to define a
date format object to parse the dates properly. If you don't specify a year,
it sets it arbitrarily to 0001, so it's less work just to add years."

# ╔═╡ 0cd27fa0-2e58-4754-973a-66d758560dbc
df = DateFormat("m/d/y")

# ╔═╡ b0c6ce0c-5909-4d08-a60a-267f1bf7a2e5
Date("4/1",df);

# ╔═╡ 0e2403f7-02f9-48a7-8954-0534fc90bcae
dates2 = Date.(dates,df);

# ╔═╡ e840276d-787b-41e2-8dff-b4f522307a0f
plot(dates2,seconds,yaxis=:log,xlabel="Dates",ylabel="Seconds",title="Benchmark speed for Ethane/6-31G Integrals",label=false)

# ╔═╡ Cell order:
# ╠═1d2793ad-34cb-4704-b139-f55e9af06317
# ╠═3957e56c-2d81-4cb8-9825-963ae0dc49e4
# ╠═1196b181-459a-47e8-b57e-f320b0471a89
# ╠═9b3d18a3-a176-4402-874c-45f15d1f6e14
# ╟─59ebedbe-b0c9-439b-bec6-a0aa094e6c84
# ╠═0cd27fa0-2e58-4754-973a-66d758560dbc
# ╠═b0c6ce0c-5909-4d08-a60a-267f1bf7a2e5
# ╠═0e2403f7-02f9-48a7-8954-0534fc90bcae
# ╠═e840276d-787b-41e2-8dff-b4f522307a0f


text:
using PGFPlotsXExamples
import PGFPlotsX

using Base.Test
using NBInclude
using LaTeXStrings
using Contour
using Colors
using DataFrames

const EXAMPLE_DIR = joinpath(@__DIR__, "..", "examples/")
const MIME_types = [PGFPlotsX.Plot, PGFPlotsX.AxisLike, PGFPlotsX.TikzPicture, PGFPlotsX.TikzDocument]

ENV["TESTING_PGFPLOTSX"] = true

function export_ans(x)
    if any(issubtype.(typeof(x), MIME_types))
        folder, file = tempdir(), tempname()
        pdffile = joinpath(folder, file) * ".pdf"
        println("saving ", typeof(x), " to ", joinpath(folder, file), ".pdf")
        PGFPlotsX.save(pdffile, x)
        rm(pdffile)
    else
        println("skipping $(typeof(x))")
    end
end

PGFPlotsX.latexengine!(PGFPlotsX.PDFLATEX)

cd(EXAMPLE_DIR) do
    for file in readdir()
        if endswith(file, ".ipynb")
            println("Testing ", file)
            nbinclude(joinpath(EXAMPLE_DIR, file); anshook = export_ans)
        end
    end
end

max_stars_repo_name: KristofferC/PGFPlotsXExamples

max_issues_repo_name: KristofferC/PGFPlotsXExamples

max_forks_repo_name: KristofferC/PGFPlotsXExamples


text:
"""
    listholidays(calendar, dt0::Dates.Date, dt1::Dates.Date) → Vector{Dates.Date}

Returns the list of holidays between `dt0` and `dt1`.
"""
function listholidays(hc::HolidayCalendar, dt0::Dates.Date, dt1::Dates.Date)
    d0 = min(dt0, dt1)
    d1 = max(dt0, dt1)
    dt_range = d0:Dates.Day(1):d1
    isholiday_vec = [ isholiday(hc, i) for i in dt_range ]
    return dt_range[isholiday_vec]
end

listholidays(calendar, dt0::Dates.Date, dt1::Dates.Date) = listholidays(convert(HolidayCalendar, calendar), dt0, dt1)

"""
    listbdays(calendar, dt0::Dates.Date, dt1::Dates.Date) → Vector{Dates.Date}

Returns the list of business days between `dt0` and `dt1`.
"""
function listbdays(hc::HolidayCalendar, dt0::Dates.Date, dt1::Dates.Date)
    d = tobday(hc, min(dt0, dt1))
    d1 = max(dt0, dt1)

    # empty result
    if d > d1
        return Vector{Dates.Date}()
    end

    n = Dates.value(d1) - Dates.value(d) + 1
    raw_vec = Vector{Dates.Date}(undef, n)
    raw_vec[1] = d
    d = advancebdays(hc, d, 1)
    i = 2
    while d <= d1
        raw_vec[i] = d
        i += 1
        d = advancebdays(hc,d,1)
    end

    return raw_vec[1:(i-1)]
end

listbdays(calendar, dt0::Dates.Date, dt1::Dates.Date) = listbdays(convert(HolidayCalendar,calendar), dt0, dt1)


text:
"""
    Formatter

A `Formatter` must implement a `format(::Formatter, ::Record)` method
which takes a `Record` and returns a `String` representation of the
log `Record`.
"""
abstract type Formatter end

const DEFAULT_FMT_STRING = "[{level} | {name}]: {msg}"
const DATE_FMT_STRING = "yyyy-mm-dd HH:MM:SS"

"""
    DefaultFormatter

The `DefaultFormatter` uses a simple format string to build
the log message. Fields from the `Record` to be used should be
wrapped curly brackets.

Ex) "[{level} | {name}]: {msg}" will print message of the form
[info | root]: my info message.
[warn | root]: my warning message.
...
"""
struct DefaultFormatter <: Formatter
    fmt_str::AbstractString
    tokens::Vector{Pair{Symbol, Bool}}
    output_tz::Dates.TimeZone

    function DefaultFormatter(
        fmt_str::AbstractString=DEFAULT_FMT_STRING,
        output_tz=localzone(),
    )
        #r"(?<={).+?(?=})
        tokens = map(eachmatch(r"({.+?})|(.+?)", fmt_str)) do m
            #println(dump(m))
            if m.captures[1] != nothing
                return Symbol(strip(m.match, ('{', '}'))) => true
            else
                return Symbol(m.match) => false
            end
        end

        new(fmt_str, tokens, output_tz)
    end
end

"""
    format(::DefaultFormatter, ::Record) -> String

Iteratively replaces entries in the format string with the appropriate fields in the
`Record`.
"""
function format(fmt::DefaultFormatter, rec::Record)
    parts = map(fmt.tokens) do token
        content = token.first
        value = content

        if token.second
            tmp_val = getproperty(rec, content)

            if content === :lookup
                name, file, line = if isa(tmp_val, StackFrame)
                    # lookup is a StackFrame
                    tmp_val.func, tmp_val.file, tmp_val.line
                else
                    "<nothing>", "", -1
                end

                value = "$(name)@$(basename(string(file))):$(line)"
            elseif content === :stacktrace
                # stacktrace is a vector of StackFrames
                str_frames = map(tmp_val) do frame
                    string(frame.func, "@", basename(string(frame.file)), ":", frame.line)
                end

                value = string(" stack:[", join(str_frames, ", "), "]")
            elseif content === :date
                value = if tmp_val isa ZonedDateTime
                    Dates.format(astimezone(tmp_val, fmt.output_tz), DATE_FMT_STRING)
                elseif tmp_val isa DateTime
                    Dates.format(tmp_val, DATE_FMT_STRING)
                else
                    tmp_val
                end
            else
                value = tmp_val
            end
        end

        return value
    end

    return string(parts...)
end

struct DictFormatter <: Formatter
    aliases::Union{Dict{Symbol, Symbol}, Nothing}
    serializer::Function
end

"""
    DictFormatter([aliases, serializer])

Formats the record to Dict that is amenable to serialization formats such as JSON and then
runs the serializer function on the produced dictionary.

# Arguments
- `aliases::Dict{Symbol, Symbol}`: Mapping where the keys represent aliases and values
represent existing record attributes to include in the dictionary (defaults to all
attributes).
- `serializer::Function`: A function that takes a Dictionary and returns a string. Defaults
to `string(dict)`.
"""
DictFormatter() = DictFormatter(nothing, string)
DictFormatter(aliases::Dict{Symbol, Symbol}) = DictFormatter(aliases, string)
DictFormatter(serializer::Function) = DictFormatter(nothing, serializer)

"""
    format(::DictFormatter, ::Record) -> Dict

Converts :date, :lookup and :stacktrace to strings and dicts respectively.
"""
function format(fmt::DictFormatter, rec::Record)
    aliases = if fmt.aliases === nothing
        Dict(key => key for key in keys(rec))
    else
        fmt.aliases
    end

    dict = Dict{Symbol, Any}()

    for (alias, key) in aliases
        tmp_val = getproperty(rec, key)

        if key === :date
            value = string(tmp_val)
        elseif key === :lookup
            value = if isa(tmp_val, StackFrame)
                Dict(
                    :name => tmp_val.func,
                    :file => basename(string(tmp_val.file)),
                    :line => tmp_val.line
                )
            else
                Dict(
                    :name => "<nothing>",
                    :file => "",
                    :line => -1
                )
            end
        elseif key === :stacktrace
            value = map(
                frame -> Dict(
                    :name => frame.func,
                    :file => basename(string(frame.file)),
                    :line => frame.line
                ),
                tmp_val
            )
        else
            value = tmp_val
        end

        dict[alias] = value
    end

    return fmt.serializer(dict)
end


text:
# Julia code
# Copyright (c) 2022 the SLICOTMath.jl developers
# Portions extracted from SLICOT-Reference distribution:
# Copyright (c) 2002-2020 NICONET e.V.

function run_mb01td(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 20
    LDA = NMAX
    LDB = NMAX
    LDWORK = NMAX-1
    A = Array{Float64,2}(undef, LDA,NMAX)
    B = Array{Float64,2}(undef, LDB,NMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
if ( N<0 || N>NMAX )
else

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       B[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
# interp call 1

    INFO = SLICOT.mb01td!(N, A, B)
    @test INFO == 0
    INFO == 0 || return

if ( INFO!=0 )
else
# interp output 1
    println(io,"B:")
    _nc = N
    _nr = N
    show(io,"text/plain",B[1:_nr,1:_nc])
    println(io,)

end # if
end # if
    close(f)
end # run_X()


text:
hx = hx - dt ( iffty(-jky ffty(ez)
	+ dt^2/24( ifftx(-kx^2 fftx ( iffty ( -jky ffty(ez)))))
        + iffty(jky^3 ffty(ez))))

hy = hy + dt ( ifftx(-jkx fftx(ez))
	+ dt^2/24 *(iffty(-ky^2 ffty(ifftx(-jkx fftx(ez))))
	+ ifftx(jkx^3 fftx(ez))))

ez = ez + dt ( ifftx(-jkx fftx(hy)) - iffty(-jky ffty(hx))
	+ dt^2/24 (
	      2 iffty( -ky^2 ffty(ifftx(-jkx fftx(hy))))
	      + ifftx( jkx^3 fftx(hy))
	      - iffty( jky^3 ffty(hx))))

ex = ex + dt ( iffty(-jky ffty(hz)
        + dt^2/24 ( ifftx(-kx^2 fftx ( iffty ( -jky ffty(hz)))))
        + iffty(jky^3 ffty(hz))))

ey = ey - dt ( ifftx(-jkx fftx(hz))
	+ dt^2/24 *(iffty(-ky^2 ffty(ifftx(-jkx fftx(hz))))
	+ ifftx(jkx^3 fftx(hz))))

hz = hz - dt (  ifftx(-jkx fftx(ey))
	      - iffty(-jky ffty(ex))
	      - dt^2/24 (
	      - iffty(-ky^2 ffty(ifftx(-jkx fftx(ey))))
	      + ifftx(-kx^2 fftx(iffty(-jky ffty(ex))))
	      + ifftx(jkx^3 fftx(ey))
              - iffty(jky^3 ffty(ex))))

omega      =   sqrt(8*pi*pi)
bz         = - cos(2*pi*x)*cos(2*pi*y)*cos(omega*t)
ex         =   cos(2*pi*x)*sin(2*pi*y)*sin(omega*t)*2*pi/ omega
ey         = - sin(2*pi*x)*cos(2*pi*y)*sin(omega*t)*2*pi/ omega
diff(ex,t) =   diff(bz,y)
diff(ey,t) = - diff(bz,x)
diff(bz,t) =   diff(ex,y) - diff(ey,x)

omega = sqrt(2*pi*pi)
ez =  cos(pi*x)*cos(pi*y)*cos(omega*t)
hx =  cos(pi*x)*sin(pi*y)*sin(omega*t) * pi / omega
hy = -sin(pi*x)*cos(pi*y)*sin(omega*t) * pi / omega

diff(hx,t) = - diff(ez,y)
diff(hy,t) =   diff(ez,x)
diff(ez,t) =   diff(hy,x) - diff(hx,y)

  d(d2hz/dy2 + d2hz/dx2)/dy
- d(d2hz/dx2 + d2hz/dy2)/dx
  d3ex/dy3 - d3dey/dx3 + d3ex/dydx2 - d3ey/dxdy2


text:
@testset "last resort" begin
    # http://www.maths.usyd.edu.au/ub/sums/puzzlehunt/2016/puzzles/A2S1_Last_Resort.pdf
    f = best_feature(["advent", "achilles", "binary", "norway", "bubbly", "yacht", "anchor"])
    @test f.description == "has 1 reverse alphabetical bigrams"
end