text:
const tf = TensorFlow

function current_device(sess)
    dvs = tf_graph(sess).op_context.devices
    isempty(dvs) ? nothing : dvs[end]
end

function run_external_optimizer(snp)
    p = open(`/home/keno/tensorflow/bazel-bin/tensorflow/compiler/xla/tools/external_optimizer`; read=true, write=true)
    buf = IOBuffer()
    writeproto(buf, snp)
    write(p, take!(buf))
    close(p.in)
    buf = IOBuffer(read(p))
    readproto(buf, HloSnapshot())
end

# Manually implement the xrt_compile op to plaster over 1.12 vs nightly
# differences. Nightly adds a second output.
function make_xrt_compile(graph)
    alloc = placeholder(String)
    desc = tf.NodeDescription(graph, "XRTCompile", tf.get_name("XRTCompile"))
    tf.add_input(desc, alloc)
    cc = tf.Tensor(tf.Operation(desc), 1)
    cc, alloc
end

mutable struct XRTCompilation
    sess
    device::Union{Nothing, TensorFlow.Device}
    # XLA Program Shape
    shape::ProgramShape
    # Julia Return Type
    rt::Type
    h::Int64
    global compile
    function compile(sess, comp::XLAComputation, rt::Type)
        buf = IOBuffer()
        # Enable this to proxy through an external XLA optimizer, to
        # try out optimizations before they're deployed to TPUs.
        #comp.hlo_snapshot = run_external_optimizer(comp.hlo_snapshot)
        writeproto(buf, comp)
        op, alloc = make_xrt_compile(tf_graph(sess))
        res = new(sess,
            current_device(sess),
            comp.hlo_snapshot.hlo.hlo_module.host_program_shape, rt,
            run(tf_session(sess), op, Dict(alloc => String(take!(buf)))))
        finalizer(close, res)
        res
    end
end
function Base.show(io::IO, comp::XRTCompilation)
    print(io, "XRTCompilation handle")
end

function Base.setproperty!(c::XRTCompilation, args...)
    error("XRTCompilation may not be modified")
end

function Base.close(c::XRTCompilation)
    if tf_session(c.sess).ptr == C_NULL
        # The session was closed so our resources are already freed
        Core.setfield!(c, :h, -1)
        return
    end
    f() = as_default(tf_graph(c.sess)) do
        try
            run(tf_session(c.sess), TensorFlow.Ops.xrt_release_compilation_handle(c.h))
        catch err
            # This runs as a finalizer. The error gets printed using the
            # C printer rather than the julia one. Transform the error
            # message to make sure it's readable
            isa(err, TensorFlow.TFException) || rethrow(err)
            error(string("TensorFlow error: ", string(err.status)))
        end
    end
    if c.device === nothing
        f()
    else
        with_device(f, c.device)
    end
    Core.setfield!(c, :h, -1)
end

function make_default_execution_config()
    return XRTExecutionConfig(rng_seed=rand(UInt32) | 0x01)
end

function make_run_op(com::XRTCompilation, inputs...; device=com.device, config=make_default_execution_config())
    iob = PipeBuffer();
    writeproto(iob, config)
    str = String(take!(iob))
    op = as_default(tf_graph(com.sess)) do
        handles = map(x->gethandle!(com.sess, x), Iterators.filter(x->representable(typeof(x)), inputs))
        # Make sure everything is on the same device
        @assert all(input->input.device == device, handles)
        str = TensorFlow.constant(str)
        _op() = TensorFlow.Ops.xrt_execute(com.h, str, collect(map(x->x.h, handles)))
        return device === nothing ? _op() : with_device(_op, device)
    end
end

mutable struct XRTAllocation
    sess
    device::Union{Nothing, TensorFlow.Device}
    h::Int64
    function XRTAllocation(sess, literal::LiteralProto)
        buf = IOBuffer()
        device = current_device(sess)
        device_ordinal = 0
        if device !== nothing && device.parts[end].kind[2] == "TPU"
            device_ordinal = device.parts[end].index - 1
        end
        writeproto(buf, XLAAllocation(
            device_ordinal = device_ordinal,
            value = literal))
        local alloc
        op = as_default(tf_graph(sess)) do
            alloc = placeholder(String)
            TensorFlow.Ops.xrt_allocate(alloc)
        end
        h = run(sess, op, Dict(alloc=>String(take!(buf))))
        #ccall(:jl_, Cvoid, (Any,), "Got allocation $(h)")
        res = new(sess,
            current_device(sess),
            h)
        finalizer(close, res)
        res
    end

    global _run
    function _run(com::XRTCompilation, inputs...; config=make_default_execution_config())
        op = make_run_op(com, inputs...; config=config)
        h = run(com.sess, op, Dict())
        #ccall(:jl_, Cvoid, (Any,), "Got allocation $(h)")
        res = @GC.preserve com inputs new(com.sess, com.device, h)
        finalizer(close, res)
        res
    end

    global run
    function run(com::XRTCompilation, inputs...; device=com.device, config=make_default_execution_config())
        if typeof(device) <: DeviceIndex
            device = tf_tpu_device(com.sess, device)
        end
        op = make_run_op(com, inputs...; device=device, config=config)
        h = @GC.preserve com inputs run(com.sess, op, Dict(); async=true)
        #ccall(:jl_, Cvoid, (Any,), "Got allocation $(h)")
        res = new(com.sess, com.device, h)
        finalizer(close, res)

        return if com.rt <: XRTArray
            T = convert(Type, XlaType(com.shape.result.element_type))
            dims = (com.shape.result.dimensions...,)
            XRTArray{T, dims, length(dims)}(res)
        else
            XRTRemoteStruct{com.rt}(res)
        end
    end

    global run_on_devices
    function run_on_devices(compilation_handle, args...; devices=all_tpu_devices(compilation_handle.sess))
		sess = compilation_handle.sess

		# Helper functions to prepare an argument for sending to the TPU
		prepare_arg(x::Number) = XRTArray(sess, x)
		prepare_arg(x::AbstractArray) = XRTArray(sess, x)
		prepare_arg(x::Tuple) = prepare_arg.(x)
		prepare_arg(x) = XRTRemoteStruct(sess, x)

		# For each TPU device
		worker_ops = Any[]
		all_remote_handles = Any[compilation_handle]
		for device in devices
			# Prepare arguments
			compilation_args = with_device(tf_tpu_device(sess, device)) do
				# Drop `Val`'s because we don't need to pass those in
				return prepare_arg.(filter(x -> !isa(x, Val), collect(args)))
			end
			push!(all_remote_handles, compilation_args)

			# Place an XRT execute op on the device
			push!(worker_ops, make_xrt_execute_on(device, compilation_handle, compilation_args...))
		end

        # Next, launch these worker ops
        task = @async begin
			rets = @GC.preserve all_remote_handles run(sess, worker_ops; async=true)

			# Wrap rets as XRTAllocations, then as XRTRemoteStructs
			rets = [new(sess, tf_tpu_device(sess, devices[idx]), rets[idx]) for idx in 1:length(devices)]
			return XRTRemoteStruct{compilation_handle.rt}.(rets)
		end

        return task
    end
end
function Base.setproperty!(c::XRTAllocation, args...)
    error("XRTAllocation may not be modified")
end

const tf = TensorFlow
global counter = 0
function Base.close(c::XRTAllocation)
     if tf_session(c.sess).ptr == C_NULL
        # The session was closed so our resources are already freed
        Core.setfield!(c, :h, -1)
        return
    end
    #ccall(:jl_, Cvoid, (Any,), "Start Releasing allocation $(c.h)")
    global counter
    counter += 1
    desc = tf.NodeDescription(tf_graph(c.sess), "Const", "ReleaseAllocationHandle$(counter)/Const")
    desc["value"] = TensorFlow.RawTensor(c.h)
    desc["dtype"] = typeof(c.h)
    cc = tf.Tensor(tf.Operation(desc))
    desc2 = tf.NodeDescription(tf_graph(c.sess), "XRTReleaseAllocationHandle", "ReleaseAllocationHandle$(counter)")
    tf.add_input(desc2, cc)
    c.device !== nothing && tf.set_device(desc2, c.device)
    #ccall(:jl_, Cvoid, (Any,), "Releasing allocation $(c.h)")
    try
        run(c.sess, tf.Tensor(tf.Operation(desc2)))
    catch err
        isa(err, TensorFlow.TFException) || rethrow(err)
        error(string("TensorFlow error: ", string(err.status)))
    end
    #ccall(:jl_, Cvoid, (Any,), "Done releasing allocation $(c.h)")
    Core.setfield!(c, :h, -1)
end

mutable struct XRTStorage{T, N}
    localstorage::Union{Array{T, N}, Nothing}
    remotestorage::Union{XRTAllocation, Nothing}
    function XRTStorage{T, N}(localstorage::Union{Array{T, N}, Nothing},
                              remotestorage::Union{XRTAllocation, Nothing}) where {T, N}
        @assert localstorage !== nothing || remotestorage !== nothing
        new{T, N}(localstorage, remotestorage)
    end
    function XRTStorage{T, N}(a::Array{T, N}) where {T, N}
        new{T, N}(a, nothing)
    end
    function XRTStorage{T, N}(a::XRTAllocation) where {T, N}
        new{T, N}(nothing, a)
    end
end

struct XRTArray{T, Dims, N} <: AbstractArray{T, N}
    storage::XRTStorage{T, N}
    @noinline function XRTArray{T, Dims, N}(a::Array{T, N}) where {T, Dims, N}
        @assert size(a) == Dims
        @assert length(Dims) == N
        new{T, Dims, N}(XRTStorage{T, N}(a))
    end
    @noinline function XRTArray{T, (), 0}(a::T) where {T<:XLAScalar}
        XRTArray{T, (), 0}(fill(a))
    end
    @noinline function XRTArray{T, (), 0}(a::T) where {T<:XRTArray}
        # Can share storage
        storage = a.storage.remotestorage !== nothing ?
            XRTStorage{T, 0}(a.storage.remotestorage) :
            XRTStorage{T, 0}(fill(a))
        new{T, (), 0}(storage)
    end
    @noinline function XRTArray{T, Dims, N}(h::XRTAllocation) where {T, Dims, N}
        @assert length(Dims) == N
        new{T, Dims, N}(XRTStorage{T, N}(h))
    end
end

function read_literal(rs::XRTAllocation)
    sess, dev, h = rs.sess, rs.device, rs.h
    op() = as_default(tf_graph(sess)) do
        run(sess, TensorFlow.Ops.xrt_read_literal(h))::String
    end
    literal = dev === nothing ? op() : with_device(op, dev)
    readproto(IOBuffer(literal), LiteralProto())
end

function Base.convert(::Type{Array}, A::XRTArray)
    if A.storage.localstorage !== nothing
        return copy(A.storage.localstorage)
    end
    rs = A.storage.remotestorage
    lit = read_literal(rs)
    A.storage.localstorage = convert(Array, lit)
    return copy(A.storage.localstorage)
end

@noinline function Base.convert(::Type{T}, A::XRTArray{T, (), 0}) where {T}
    convert(Array, A)[]::T
end

function gethandle!(sess, A::XRTArray)
    if A.storage.remotestorage !== nothing
        @assert A.storage.remotestorage.sess === sess
        return A.storage.remotestorage
    end
    A.storage.remotestorage = XRTAllocation(sess, convert(LiteralProto, A.storage.localstorage))
    A.storage.remotestorage
end
function gethandle!(sess, x::XRTAllocation)
    @assert sess == x.sess
    return x
end

const XRTMatrix{T, Dims} = XRTArray{T, Dims, 2} where {T, Dims}
const XRTVector{T, Dims} = XRTArray{T, Dims, 1} where {T, Dims}
XRTArray(A::AbstractArray) = XRTArray(collect(A)::Array)
function XRTArray(a::Array{T}) where {T}
    XRTArray{T, size(a), ndims(a)}(a)
end
function XRTArray(a::XLAScalar)
    XRTArray{typeof(a), (), 0}(a)
end
function XRTArray(a::Tuple{XLAScalar})
    XRTArray{typeof(a[1]), (), 0}(a[1])
end
function XRTArray(a::Tuple{XRTArray})
    XRTArray{typeof(a[1]), (), 0}(a[1])
end
XRTArray(sess, A::AbstractArray) = XRTArray(sess, collect(A)::Array)
XRTArray(sess, a::XLAScalar) = XRTArray(sess, fill(a))
function XRTArray(sess, a::Array{T}) where {T}
    ret = XRTArray{T, size(a), ndims(a)}(a)
    gethandle!(sess, ret)
    ret
end

Base.eltype(A::Type{<:XRTArray{T}}) where {T} = T
Base.size(A::Type{<:XRTArray{T, Dims}} where T) where {Dims} = Dims
Base.size(A::Type{<:XRTArray{T, Dims}} where T, i) where {Dims} = i <= length(Dims) ? Dims[i] : 1
@inline function Base.axes(A::Type{<:XRTArray{<:Any, Dims}}, d) where {Dims}
    d <= length(Dims) ? axes(A)[d] : Base.OneTo(1)
end
Base.length(A::Type{<:XRTArray}) = prod(size(A))

Base.eltype(A::XRTArray{T}) where {T} = T
Base.size(A::XRTArray{T, Dims}) where {T, Dims} = Dims
Base.size(A::XRTArray{T, Dims}, i) where {T, Dims} = i <= length(Dims) ? Dims[i] : 1
@inline function Base.axes(A::XRTArray{<:Any, Dims}, d) where {Dims}
    d <= length(Dims) ? axes(A)[d] : Base.OneTo(1)
end

import .xla: Shape
function Shape(::Type{XRTArray{T, Dims, N}} where N) where {T, Dims}
    Shape(T, Dims)
end
Shape(T::Type, SHP::Tuple) = Shape(convert(XlaType, T), SHP)
function Shape(::Type{<:XRTArray{T, SHP1}}, SHP2::Tuple) where {T, SHP1}
    # Map arrays of arrays to higher rank arrays
    Shape(T, (SHP1..., SHP2...))
end
function Shape(XLAT::XlaType, SHP::Tuple)
    perm = sortperm(collect(SHP); rev=true)
    Shape(
        element_type = XLAT.which,
        dimensions = Int64[SHP...],
        layout = Layout(
            format = Format.DENSE,
            minor_to_major = perm .- 1,
            max_sparse_elements = 0
        )
    )
end

Base.convert(::Type{Shape}, AT::Type{XRTArray{T, Dims, N}} where N) where {T, Dims} =
    Shape(AT)

shape(A::XRTArray) = Shape(typeof(A))

Base.isempty(A::XRTArray{T, Dims}) where {T, Dims} = prod(Dims) == 0

Base.print_array(io::IO, A::XRTArray) = Base.print_array(io, convert(Array, A))
Base.show_vector(io::IO, A::XRTArray, opn='[', cls=']') =
    Base.show_vector(io, convert(Array, A), opn, cls)
Base._show_nonempty(io::IO, A::XRTArray, prefix::String) =
    Base._show_nonempty(io, convert(Array, A), prefix)
Base._show_nonempty(io::IO, A::XRTMatrix, prefix::String) =
    Base._show_nonempty(io, convert(Array, A), prefix)

# References a remote struct. We currently have a tension between
# local structs of remote data and a remote reference to a struct.
# Ideally we'd unify them eventually.
struct XRTRemoteStruct{T}
    storage::XRTAllocation
end
function gethandle!(sess, A::XRTRemoteStruct)
    @assert sess == A.storage.sess
    A.storage
end
Base.show(io::IO, x::XRTRemoteStruct{T}) where {T} =
    println("Remote $(T)")

struct_to_literal(A::XRTArray) = convert(LiteralProto, convert(Array, A))
struct_to_literal(x::XLA.XLAScalar) = struct_to_literal(XRTArray(x))
function struct_to_literal(x)
    T = typeof(x)
    LiteralProto(
        shape = dtype_to_shape(T),
        tuple_literals = collect(struct_to_literal(getfield(x, i)) for i = 1:fieldcount(T) if representable(fieldtype(T, i)))
    )
end

function literal_to_struct(::Type{<:XRTArray}, x::LiteralProto)
    # TODO: We could instead request the remote to create this for
    # us. Otherwise we're doing a round trip here. Good enough for now.
    XRTArray(convert(Array, x))
end

function literal_to_struct(::Type{<:XLA.XLAScalar}, x::LiteralProto)
    # TODO: We could instead request the remote to create this for
    # us. Otherwise we're doing a round trip here. Good enough for now.
    convert(Array, x)[]
end

function literal_to_struct(T::Type, x::LiteralProto)
    reconstructed_fields = Any[]
    cur_proto_field = 1
    for i = 1:fieldcount(T)
        fT = fieldtype(T, i)
        if representable(fT)
            push!(reconstructed_fields,
                literal_to_struct(fT, x.tuple_literals[cur_proto_field]))
            cur_proto_field += 1
        else
            push!(reconstructed_fields, fT.instance)
        end
    end
    eval(Expr(:new, T, reconstructed_fields...))
end

function XRTRemoteStruct{T}(sess, x::T) where {T}
    XRTRemoteStruct{T}(XRTAllocation(sess, struct_to_literal(x)))
end
XRTRemoteStruct(sess, x::T) where {T} = XRTRemoteStruct{T}(sess, x)

function Base.convert(::Type{T}, strct::XRTRemoteStruct{T}) where {T}
    proto = read_literal(strct.storage)
    literal_to_struct(T, proto)
end

Base.fetch(strct::XRTRemoteStruct{T}) where {T} = convert(T, strct)


text:
#   This file is part of UnitSystems.jl
#   It is licensed under the MIT license
#   UnitSystems Copyright (C) 2020 Michael Reed
#       _           _                         _
#      | |         | |                       | |
#   ___| |__   __ _| | ___ __ __ ___   ____ _| | __ _
#  / __| '_ \ / _` | |/ / '__/ _` \ \ / / _` | |/ _` |
# | (__| | | | (_| |   <| | | (_| |\ V / (_| | | (_| |
#  \___|_| |_|\__,_|_|\_\_|  \__,_| \_/ \__,_|_|\__,_|
#
#   https://github.com/chakravala
#   https://crucialflow.com

# angle

@pure radian(U::UnitSystem) = angle(one(U),U,Metric)
@pure steradian(U::UnitSystem) = solidangle(one(U),U,Metric)
@pure degree(U::UnitSystem) = angle(turn(U)/two(U)^3/three(U)^2/five(U),U,Metric)
@pure gradian(U::UnitSystem) = angle(turn(U)/two(U)^4/five(U)^2,U,Metric)
@pure arcminute(U::UnitSystem) = degree(U)/two(U)^2/three(U)/five(U)
@pure arcsecond(U::UnitSystem) = arcminute(U)/two(U)^2/three(U)/five(U)

@pure apm(U::UnitSystem) = one(U)/minute(U)
@pure rpm(U::UnitSystem) = turn(U)/minute(U)

# time

@pure second(U::UnitSystem) = time(one(U),U,Metric)
@pure minute(U::UnitSystem) = two(U)^2*three(U)*five(U)*second(U)
@pure hour(U::UnitSystem) = two(U)^2*three(U)*five(U)*minute(U)
@pure day(U::UnitSystem) = two(U)^3*three(U)*hour(U)
@pure year(U::UnitSystem) = aⱼ*day(U)
@pure radarmile(U::UnitSystem) = two(U)*nauticalmile(U)/lightspeed(U)

# length

@pure meter(U::UnitSystem) = length(one(U),U,Metric)
@pure earthmeter(U::UnitSystem) = length(one(U),U,Meridian)
#@pure navigationmeter(U::UnitSystem) = greatcircle(U)/two(U)^9/five(U)^7
@pure angstrom(U::UnitSystem) = hecto(U)*pico(U)*meter(U)
@pure foot(U::UnitSystem) = length(one(U),U,English)
@pure inch(U::UnitSystem) = length(one(U),U,IPS)
#@pure rackunit(U::UnitSystem) = foot(U)*seven(U)/two(U)^4/three(U)
@pure yard(U::UnitSystem) = three(U)*foot(U)
@pure surveyfoot(U::UnitSystem) = length(one(U),U,Survey)
@pure statutemile(U::UnitSystem) = length(two(U)^5*three(U)*five(U)*eleven(U),U,Survey)
@pure earthradius(U::UnitSystem) = sqrt(earthmass(U)*gravitation(U)/gforce(U))
#@pure navigationradius(U::UnitSystem) = length(sqrt(earthmass(Metric)*gravitation(Metric)/(gforce(Metric) - gravitation(Metric)*(lunarmass(Metric)/lunardistance(Metric)^2+solarmass(Metric)/astronomicalunit(Metric)^2))),U,Metric)
@pure greatcircle(U::UnitSystem) = normal(turn(U))*earthradius(U)
#@pure greatcircle(U::UnitSystem) = normal(turn(U))*navigationradius(U)
@pure nauticalmile(U::UnitSystem) = length(one(U),U,Nautical)
#@pure navigationmile(U::UnitSystem) = greatcircle(U)/two(U)^5/three(U)^3/five(U)^2
@pure astronomicalunit(U::UnitSystem) = length(𝟏,U,IAU)
@pure lunardistance(U::UnitSystem) = length(𝟏,U,IAUE)
@pure jupiterdistance(U::UnitSystem) = length(𝟏,U,IAUJ)
@pure mile(U::UnitSystem) = length(two(U)^5*three(U)*five(U)*eleven(U),U,English)
@pure admiraltymile(U::UnitSystem) = length(two(U)^6*five(U)*nineteen(U),U,English)
@pure meridianmile(U::UnitSystem) = length(two(U)^4*five(U)^5/three(U)^3,U,Metric)
@pure lightyear(U::UnitSystem) = year(U)*lightspeed(U)
@pure parsec(U::UnitSystem) = astronomicalunit(U)*two(U)^7*three(U)^4*five(U)^3/turn(U)

# area

@pure barn(U::UnitSystem) = area((two(U)*five(U))^-28,U,Metric)
@pure hectare(U::UnitSystem) = area(hecto(U)*hecto(U),U,Metric)
@pure acre(U::UnitSystem) = area(two(U)^-7/five(U),U,MPH)
@pure surveyacre(U::UnitSystem) = area(two(U)^3*three(U)^2*five(U)*eleven(U)^2,U,Survey)
#@pure township(U::UnitSystem) = two(U)^9*three(U)^2*five(U)*surveyacre(U)
#@pure footballfield(U::UnitSystem) = area(two(U)^8*three(U)^2*five(U)^2,U,English)

# volume

@pure gallon(U::UnitSystem) = volume(three(U)*seven(U)*eleven(U),U,IPS)
@pure liter(U::UnitSystem) = volume(inv((two(U)*five(U))^3),U,Metric)
@pure quart(U::UnitSystem) = gallon(U)/two(U)^2
@pure pint(U::UnitSystem) = quart(U)/two(U)
@pure cup(U::UnitSystem) = pint(U)/two(U)
@pure fluidounce(U::UnitSystem) = cup(U)/two(U)^3
@pure teaspoon(U::UnitSystem) = five(U)*milli(U)*liter(U)
@pure tablespoon(U::UnitSystem) = three(U)*teaspoon(U)
#@pure oilbarrel(U::UnitSystem) = two(U)*three(U)*seven(U)*gallon(U)

# speed

@pure bubnoff(U::UnitSystem) = meter(U)/year(U)
@pure ips(U::UnitSystem) = inch(U)/second(U)
@pure fps(U::UnitSystem) = foot(U)/second(U)
@pure fpm(U::UnitSystem) = foot(U)/minute(U)
@pure ms(U::UnitSystem) = meter(U)/second(U)
@pure kmh(U::UnitSystem) = kilo(U)*meter(U)/hour(U)
@pure mph(U::UnitSystem) = mile(U)/hour(U)
@pure knot(U::UnitSystem) = nauticalmile(U)/hour(U)
@pure mps(U::UnitSystem) = mile(U)/second(U)

# mass

@pure grain(U::UnitSystem) = milli(U)*pound(U)/seven(U)
@pure gram(U::UnitSystem) = mass(one(U),U,Gauss)
@pure earthgram(U::UnitSystem) = mass(milli(U),U,Meridian)
@pure kilogram(U::UnitSystem) = mass(one(U),U,Metric)
@pure tonne(U::UnitSystem) = mass(kilo(U),U,Metric)
@pure ton(U::UnitSystem) = mass(two(U)*kilo(U),U,English)
@pure pound(U::UnitSystem) = mass(one(U),U,English)
@pure ounce(U::UnitSystem) = mass(two(U)^-4,U,English)
@pure slug(U::UnitSystem) = mass(one(U),U,British)
@pure slinch(U::UnitSystem) = mass(one(U),U,IPS)
@pure hyl(U::UnitSystem) = mass(one(U),U,GravitationalMetric)

# force

@pure dyne(U::UnitSystem) = force(one(U),U,Gauss)
@pure newton(U::UnitSystem) = force(one(U),U,Metric)
@pure poundal(U::UnitSystem) = force(one(U),U,FPS)
@pure poundforce(U::UnitSystem) = force(one(U),U,English)
@pure kilopond(U::UnitSystem) = force(one(U),U,MetricEngineering)

# pressure

@pure pascal(U::UnitSystem) = pressure(one(U),U,Metric)
@pure bar(U::UnitSystem) = pressure(hecto(U)*kilo(U),U,Metric)
@pure barye(U::UnitSystem) = pressure(one(U),U,Gauss)
@pure psi(U::UnitSystem) = pressure(one(U),U,IPS)
@pure technicalatmosphere(U::UnitSystem) = kilopond(U)/(centi(U)*meter(U))^2
@pure atmosphere(U::UnitSystem) = pressure(atm,U,Metric)
@pure inchmercury(U::UnitSystem) = pressure(inHg,U,Metric)
@pure torr(U::UnitSystem) = pressure(atm/(two(U)^3*five(U)*nineteen(U)),U,Metric)

# energy

@pure electronvolt(U::UnitSystem) = elementarycharge(U)*electricpotential(one(U),U,SI2019)
@pure erg(U::UnitSystem) = energy(one(U),U,Gauss)
@pure joule(U::UnitSystem) = energy(one(U),U,Metric)
@pure footpound(U::UnitSystem) = poundforce(U)*foot(U)
@pure meancalorie(U::UnitSystem) = energy(two(U)^2*five(U)*three(U)^2/fourtythree(U),U,InternationalMean)
@pure kilocalorie(U::UnitSystem) = energy(two(U)^5*five(U)^4*three(U)^2/fourtythree(U),U,International)
@pure calorie(U::UnitSystem) = kilocalorie(U)/(two(U)*five(U))^3
@pure earthcalorie(U::UnitSystem) = molaramount(temperature(calorie(U),Metric,Meridian),Metric,Meridian)
@pure thermalunit(U::UnitSystem) = mass(temperature(kilocalorie(U),Metric,English),Metric,English)
@pure tontnt(U::UnitSystem) = giga(U)*calorie(U)
@pure gasgallon(U::UnitSystem) = two(U)*three(U)*nineteen(U)*kilo(U)*thermalunit(U)

# power

@pure watt(U::UnitSystem) = power(one(U),U,Metric)
@pure tonsrefrigeration(U::UnitSystem) = frequency(two(U)*five(U)/three(U),U,Metric)*thermalunit(U)
@pure boilerhorsepower(U::UnitSystem) = frequency(Constant(1339)/(two(U)^4*three(U)^2),U,Metric)*thermalunit(U)
# thermalconductivity_water(British) ≈ 0.5778
@pure thermalconductivity_water(U::UnitSystem) = thermalconductivity((two(U)^2*three(U)*five(U))^2/thermalunit(U),U,Metric)
@pure horsepower(U::UnitSystem) = power(two(U)*five(U)^2*eleven(U),U,British)
@pure horsepowerwatt(U::UnitSystem) = power(two(U)^4*three(U)^3/five(U)*normal(twopi(U)),U,British)
@pure horsepowermetric(U::UnitSystem) = power(three(U)*five(U)^2,U,GravitationalMetric)
@pure electricalhorsepower(U::UnitSystem) = power(Constant(746),U,Metric)

# electromagnetic

@pure coulomb(U::UnitSystem) = charge(one(U),U,Metric)
@pure ampere(U::UnitSystem) = current(one(U),U,Metric)
@pure volt(U::UnitSystem) = electricpotential(one(U),U,Metric)
@pure henry(U::UnitSystem) = inductance(one(U),U,Metric)
@pure ohm(U::UnitSystem) = resistance(one(U),U,Metric)
@pure siemens(U::UnitSystem) = conductance(one(U),U,Metric)
@pure farad(U::UnitSystem) = capacitance(one(U),U,Metric)
@pure weber(U::UnitSystem) = magneticflux(one(U),U,Metric)
@pure tesla(U::UnitSystem) = magneticfluxdensity(one(U),U,Metric)
@pure abcoulomb(U::UnitSystem) = charge(one(U),U,EMU)
@pure abampere(U::UnitSystem) = current(one(U),U,EMU)
@pure abvolt(U::UnitSystem) = electricpotential(one(U),U,EMU)
@pure abhenry(U::UnitSystem) = inductance(one(U),U,EMU)
@pure abohm(U::UnitSystem) = resistance(one(U),U,EMU)
@pure abmho(U::UnitSystem) = conductance(one(U),U,EMU)
@pure abfarad(U::UnitSystem) = capacitance(one(U),U,EMU)
@pure maxwell(U::UnitSystem) = magneticflux(one(U),U,EMU)
@pure gauss(U::UnitSystem) = magneticfluxdensity(one(U),U,EMU)
@pure oersted(U::UnitSystem) = magneticfield(one(U),U,EMU)
@pure gilbert(U::UnitSystem) = abampere(U)/two(U)/turn(U)
@pure statcoulomb(U::UnitSystem) = charge(one(U),U,ESU)
@pure statampere(U::UnitSystem) = current(one(U),U,ESU)
@pure statvolt(U::UnitSystem) = electricpotential(one(U),U,ESU)
@pure stathenry(U::UnitSystem) = inductance(one(U),U,ESU)
@pure statohm(U::UnitSystem) = resistance(one(U),U,ESU)
@pure statmho(U::UnitSystem) = conductance(one(U),U,ESU)
@pure statfarad(U::UnitSystem) = capacitance(one(U),U,ESU)
@pure statweber(U::UnitSystem) = magneticflux(one(U),U,ESU)
@pure stattesla(U::UnitSystem) = magneticfluxdensity(one(U),U,ESU)
@pure earthcoulomb(U::UnitSystem) = charge(one(U),U,Meridian)

# temperature

#@pure freezing(U::UnitSystem) = temperature(T₀-milli,U,Metric)
@pure boiling(U::UnitSystem) = temperature(T₀+Constant(99.9839),U,Metric)
@pure sealevel(U::UnitSystem) = temperature(T₀+𝟑*𝟓,U,Metric)
@pure kelvin(U::UnitSystem) = temperature(one(U),U,Metric)
@pure rankine(U::UnitSystem) = temperature(one(U),U,English)
@pure celsius(U::UnitSystem) = temperature(T₀,U,Metric)
@pure fahrenheit(U::UnitSystem) = temperature(Constant(459.67),U,English)
#@pure delisle(U::UnitSystem) = temperature(two(U)/three(U),U,Metric)
#@pure reaumur(U::UnitSystem) = temperature(five(U)/two(U)^2,U,Metric)

# mole

@pure mole(U::UnitSystem) = molaramount(one(U),U,Metric)
@pure earthmole(U::UnitSystem) = molaramount(one(U),U,Meridian)
@pure poundmole(U::UnitSystem) = molaramount(one(U),U,English)
@pure slugmole(U::UnitSystem) = molaramount(one(U),U,British)
@pure slinchmole(U::UnitSystem) = molaramount(one(U),U,IPS)

# photometric

@pure lumen(U::UnitSystem) = luminousflux(one(U),U,Metric)
@pure candela(U::UnitSystem) = luminousintensity(one(U),U,Metric)
@pure lux(U::UnitSystem) = illuminance(one(U),U,Metric)
@pure footcandle(U::UnitSystem) = illuminance(one(U),U,English)
@pure phot(U::UnitSystem) = illuminance(one(U),U,Gauss)
@pure nit(U::UnitSystem) = luminance(one(U),U,Metric)
@pure apostilb(U::UnitSystem) = luminance(two(U)/turn(U),U,Metric)
@pure stilb(U::UnitSystem) = luminance(one(U),U,Gauss)
@pure lambert(U::UnitSystem) = luminance(two(U)/turn(U),U,Gauss)
@pure footlambert(U::UnitSystem) = luminance(two(U)/turn(U),U,English)
@pure bril(U::UnitSystem) = centi(U)*nano(U)*lambert(U)

#const neper = Metric(one(U),log(𝟙))
#const bel = Metric(one(U),log10(𝟙))
#const decibel = Metric(one(U),dB(𝟙))
@pure hertz(U::UnitSystem) = one(U)/second(U)
@pure kayser(U::UnitSystem) = wavenumber(one(U),U,Gauss)
@pure diopter(U::UnitSystem) = wavenumber(one(U),U,Metric)
@pure gforce(U::UnitSystem) = specificforce(one(U),U,English)
@pure galileo(U::UnitSystem) = specificforce(one(U),U,Gauss)
@pure eotvos(U::UnitSystem) = specificforce(nano(U),U,Gauss)/length(one(U),U,Gauss)
@pure darcy(U::UnitSystem) = area(milli(U)/normal(atmosphere(Metric)),U,Gauss)
@pure poise(U::UnitSystem) = viscosity(one(U),U,Gauss)
@pure reyn(U::UnitSystem) = viscosity(one(U),U,IPS)
@pure stokes(U::UnitSystem) = diffusivity(one(U),U,Gauss)
@pure rayl(U::UnitSystem) = specificimpedance(one(U),U,Metric)
@pure katal(U::UnitSystem) = catalysis(one(U),U,Metric)
@pure mpge(U::UnitSystem) = mile(U)/gasgallon(U)
@pure langley(U::UnitSystem) = calorie(U)/(centi(U)*meter(U))^2
@pure jansky(U::UnitSystem) = fluence((Constant(1.0)*deci(U))^26,U,Metric)
@pure solarflux(U::UnitSystem) = hecto(U)^2*jansky(U)
@pure curie(U::UnitSystem) = Constant(37)*giga(U)*hertz(U)
@pure sievert(U::UnitSystem) = energy(one(U),U,Metric)/mass(U,Metric)
@pure rem(U::UnitSystem) = centi(U)*sievert(U)
@pure roentgen(U::UnitSystem) = chargedensity(one(U),U,ESU)/density(Constant(1.293),U,Metric)

ext: jl


text:
mutable struct BibliographyNode
	pages::Vector{String}
	bibfile::String
    bibdatabase::Dict{String, Any}
    database::Dict{String, Tuple{Any, Any}}
    style::Union{Symbol, Nothing}
    build::String
    source::String
    keys::Vector

    function BibliographyNode(;
    		Pages = [],
    		BibFile = "src/assets/bibliography.bib",
            Style = nothing,
            build,
            source,
            others...
        )

        bibdatabase = read_bib_file(BibFile)
        database = Dict{String, Tuple{Any, Any}}()

        new(Pages, BibFile, bibdatabase, database, Style, build, source, [])
    end
end

struct CitationNode
    key::String
    alttext::Union{String, Nothing}
    bibnode::Union{BibliographyNode, Nothing}
end


text:
module TestExpressionsAliases

import Test: @testset, @test
import Volcanito: get_alias

test_cases = (
    (
        input = :x,
        output = (:x, :x),
    ),
    (
        input = :(x = a + sin(b)),
        output = (:x, :(a + sin(b))),
    ),
    (
        input = :(x = 1),
        output = (:x, 1),
    ),
    (
        input = :(x = 1 + y),
        output = (:x, :(1 + y)),
    ),
    (
        input = :(1 + y),
        output = (Symbol("1 + y"), :(1 + y)),
    ),
)

@testset "get_alias" begin
    for case in test_cases
        @test get_alias(case.input) == case.output
    end
end

end


text:
# Input to Configuration

function load_lammps(file_path :: String; atom_names = nothing, radii = [3.5], weights = [1.0], boundary_type = ["p", "p", "p"], units = "lj")
    lines = readlines(file_path)
    dict = Dict([("atoms", 0), ("bonds", 0), ("angles", 0), ("dihedrals", 0), ("impropers", 0),
            ("atom types", 0), ("bond types", 0), ("angle types", 0), ("dihedral types", 0), ("improper types", 0)])

    for key = keys(dict)
        line_no = findall(occursin.(key, lines))
        if isempty(line_no)
            continue
        else
            dict[key] = parse(Int, split(lines[line_no[1]], " ")[1])
        end
    end

    Atoms = Vector{Atom}(undef, dict["atoms"])
    Angles = Vector{Angle}(undef, dict["angles"])
    Bonds = Vector{Bond}(undef, dict["bonds"])
    Impropers = Vector{Improper}(undef, dict["impropers"])
    Dihedrals = Vector{Dihedral}(undef, dict["dihedrals"])

    num_atoms           = dict["atoms"]
    num_atom_types      = dict["atom types"]
    num_bond_types      = dict["bond types"]
    num_angle_types     = dict["angle types"]
    num_dihedral_types  = dict["dihedral types"]
    num_improper_types  = dict["improper types"]

    line_no = findall(occursin.("xlo", lines))[1]
    xbounds_vec = split(lines[line_no], " "); line_no += 1
    ybounds_vec = split(lines[line_no], " ");        line_no += 1
    zbounds_vec = split(lines[line_no], " ");      line_no += 1
    x_bounds            = [parse(Float64, xbounds_vec[1]), parse(Float64, xbounds_vec[2])]
    y_bounds            = [parse(Float64, ybounds_vec[1]), parse(Float64, ybounds_vec[2])]
    z_bounds            = [parse(Float64, zbounds_vec[1]), parse(Float64, zbounds_vec[2])]

    domain = Domain( SVector([x_bounds, y_bounds, z_bounds]), SVector(boundary_type) )

    if length(radii) == num_atom_types
        radii = radii
    else
        radii = radii[1] .* ones(num_atom_types)
    end
    if length(weights) == num_atom_types
        weights = weights
    else
        weights = weights .* ones(num_atom_types)
    end

    # Read Masses
    masses = zeros(num_atom_types)
    line_num_masses = findall(occursin.("Masses", lines))[1] + 1
    for i = 1:num_atom_types
        masses[i] = parse(Float64, split(lines[line_num_masses + i])[2])
    end

    # Read Per Atom data

    # Parse atom names
    if atom_names == nothing
        atom_names = [Symbol(i) for i = 1:num_atom_types]
    end

    # Prep for Atom positions
    line_num_atoms = findall(occursin.("Atoms", lines))[1] + 1

    # Prep for Velocities
    ind_velocities = findall(occursin.("Velocities", lines))
    if ~isempty(ind_velocities)
        line_num_velocity = ind_velocities[1] + 1
    end

    # Iterate through each atom
    for i = 1:num_atoms
        info = split(lines[line_num_atoms+i])
        position = SA_F64[parse(Float64, info[3]), parse(Float64, info[4]), parse(Float64, info[5])]

        ind_velocities = findall(occursin.("Velocities", lines))
        if ~isempty(ind_velocities)
            info_vel = split(lines[line_num_velocity+i])
            velocity = SA_F64[parse(Float64, info_vel[2]), parse(Float64, info_vel[3]), parse(Float64, info_vel[4])]
        else
            velocity = SVector(zeros(3))
        end
        Atoms[parse(Int, info[1])] = Atom(masses[parse(Int, info[2])],
                                          position,
                                          velocity,
                                          atom_names[parse(Int, info[2])]
                                )
    end

    c = Configuration(SVector(Atoms), num_atom_types,
                    SVector(Angles), num_angle_types,
                    SVector(Bonds), num_bond_types,
                    SVector(Impropers), num_improper_types,
                    SVector(Dihedrals), num_dihedral_types,
                    SVector(atom_names), SVector(masses),
                    SVector(radii), SVector(weights),
                    domain, units)

    return c
end

function load_lammps(file_paths :: Vector{String}; atom_names = nothing, radii = [3.5], weights = [1.0], boundary_type = ["p", "p", "p"])
    n = length(file_paths)

    r = Vector{Configuration}(undef, n)
    for i = 1:n
        r[i] = load_lammps(file_paths[i]; atom_names = nothing, radii = [3.5], weights = [1.0], boundary_type = ["p", "p", "p"])
    end

    return r
end

# SAVE
function save_lammps_data(c::Configuration; file = "./DATA")
    open(file, "w") do f
        write(f, "LAMMPS Description \n")
        write(f, "\n")

        num_atoms       = length(c.Atoms)
        num_bonds       = length(c.Bonds)
        num_angles      = length(c.Angles)
        num_dihedrals   = length(c.Dihredrals)
        num_impropers   = length(c.Impropers)
        #Specify holistic information
        for (d, text) in zip([num_atoms, num_bonds, num_angles, num_dihedrals, num_impropers],
                             ["atoms", "bonds", "angles", "dihedrals", "impropers"])
            write(f, "$d $text\n")
        end

        write(f, "\n")

        #Specify type information
        for (d, text) in zip([c.num_atom_types, c.num_bond_types, c.num_angle_types, c.num_dihedral_types, c.num_improper_types],
                            ["atom types", "bond types", "angle types", "dihedral types", "improper types"])
            if d == 0
                continue
            else
                write(f, "$d $text\n")
            end
        end

        write(f, "\n")

        #Specify bounds
        for (bound, text) in zip(c.domain.bounds, ["xlo xhi", "ylo yhi", "zlo zhi"])
            lo = bound[1]
            hi = bound[2]
            write(f, "$lo $hi $text\n")
        end

        write(f, "\n")

        #Specify Masses
        write(f, "Masses\n")
        write(f, "\n")
        for i = 1:c.num_atom_types
            mass = c.masses[i]
            write(f, "$i $mass\n")
        end

        write(f, "\n")
        #Specify Positions
        write(f, "Atoms\n")
        write(f, "\n")
        for i = 1:num_atoms
            pos = c.Atoms[i].Position
            atom_id = findall(c.atom_names .== pos.type)[1]
            x, y, z = pos.x, pos.y, pos.z
            write(f, "$i $atom_id $x $y $z\n")
        end

    end

end


text:
using ITensors
using ITensors.Ops
using Zygote

s = siteinds("S=1/2", 4)

function f1(x)
  y = ITensor(Op("Ry", 1; θ=x), s)
  return y[1, 1]
end

@show x = 2.0

@show f1(x)
@show f1'(x)

function f2(x)
  y = exp(ITensor(Op("Ry", 1; θ=x), s))
  return y[1, 1]
end

@show x = 2.0

@show f2(x)
@show f2'(x)

function f3(x)
  y = Op("Ry", 1; θ=x) + Op("Ry", 1; θ=x)
  return y[1].params.θ
end

@show x = 2.0

@show f3(x)
@show f3'(x)

function f4(x)
  y = ITensor(Op("Ry", 1; θ=x) + Op("Ry", 1; θ=x), s)
  return y[1, 1]
end

@show x = 2.0

@show f4(x)
@show f4'(x)

function f5(x)
  y = exp(ITensor(Op("Ry", 1; θ=x) + Op("Ry", 1; θ=x), s))
  return y[1, 1]
end

@show x = 2.0

@show f5(x)
@show f5'(x)

function f6(x)
  y = ITensor(exp(Op("Ry", 1; θ=x) + Op("Ry", 1; θ=x)), s)
  return y[1, 1]
end

@show x = 2.0

@show f6(x)
@show f6'(x)

function f7(x)
  y = ITensor(2 * Op("Ry", 1; θ=x), s)
  return y[1, 1]
end

@show x = 2.0

@show f7(x)
@show f7'(x)

function f8(x)
  y = ITensor(2 * (Op("Ry", 1; θ=x) + Op("Ry", 1; θ=x)), s)
  return y[1, 1]
end

@show x = 2.0

@show f8(x)
@show f8'(x)

function f9(x)
  y = ITensor(Op("Ry", 1; θ=x) * Op("Ry", 2; θ=x), s)
  return y[1, 1]
end

@show x = 2.0

@show f9(x)
@show f9'(x)

function f10(x)
  y = ITensor(exp(-x * Op("X", 1) * Op("X", 2)), s)
  return norm(y)
end

@show x = 2.0

@show f10(x)
@show f10'(x)

V = randomITensor(s[1], s[2])

function f11(x)
  y = exp(-x * Op("X", 1) * Op("X", 2))
  y *= exp(-x * Op("X", 1) * Op("X", 2))
  U = Prod{ITensor}(y, s)
  return norm(U(V))
end

@show x = 2.0

@show f11(x)
@show f11'(x)

function f12(x)
  y = exp(-x * (Op("X", 1) + Op("Z", 1) + Op("Z", 1)); alg=Trotter{1}(1))
  U = Prod{ITensor}(y, s)
  return norm(U(V))
end

@show x = 2.0

@show f12(x)
@show f12'(x)

## ## XXX: Error in vcat!
## function f13(x)
##   y = -x * (Op("X", 1) * Op("X", 2) + Op("Z", 1) * Op("Z", 2))
##   U = ITensor(y, s)
##   return norm(U * V)
## end
##
## @show x = 2.0
##
## @show f13(x)
## @show f13'(x)
##
## ## XXX: Error in vcat!
## function f14(x)
##   y = exp(-x * (Op("X", 1) * Op("X", 2) + Op("Z", 1) * Op("Z", 2)); alg=Trotter{1}(1))
##   U = ITensor(y, s)
##   return norm(U * V)
## end
##
## @show x = 2.0
##
## @show f14(x)
## @show f14'(x)


text:
using SimpleTraits.BaseTraits

@test istrait(IsAnything{Any})
@test istrait(IsAnything{Union{}})
@test istrait(IsAnything{Int})

@test !istrait(IsNothing{Any})
@test !istrait(IsNothing{Union{}})
@test !istrait(IsNothing{Int})

@test istrait(IsBits{Int})
@test !istrait(IsBits{Vector{Int}})

@test istrait(IsImmutable{Float64})
@test !istrait(IsImmutable{Vector{Int}})

@test istrait(IsCallable{Function})

a = collect(1:5)
b = view(a, 2:3)
c = view(a, 1:2:5)
@test istrait(IsContiguous{typeof(b)})
@test !istrait(IsContiguous{typeof(c)})

@test istrait(IsIndexLinear{Vector})
@test !istrait(IsIndexLinear{AbstractArray})

struct T9867 end
@test istrait(IsCallable{T9867})

if VERSION<v"1-"
    @test_throws ErrorException istrait(IsIterator{Base.UnitRange})
else
    @test istrait(IsIterator{Base.UnitRange})
end


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

# Conversion/Promotion

"""
    Date(dt::DateTime) -> Date

Converts a `DateTime` to a `Date`. The hour, minute, second, and millisecond parts of
the `DateTime` are truncated, so only the year, month and day parts are used in
construction.
"""
Date(dt::TimeType) = convert(Date,dt)

"""
    DateTime(dt::Date) -> DateTime

Converts a `Date` to a `DateTime`. The hour, minute, second, and millisecond parts of
the new `DateTime` are assumed to be zero.
"""
DateTime(dt::TimeType) = convert(DateTime,dt)

Base.convert(::Type{DateTime},dt::Date) = DateTime(UTM(value(dt)*86400000))
Base.convert(::Type{Date},dt::DateTime) = Date(UTD(days(dt)))
"""
    convert{T<:Real}(::Type{T}, dt::DateTime) -> T
Converts a DateTime value `dt` to a number of type `T`. The returned value corresponds to the number of Rata Die milliseconds since epoch.
See `convert(DateTime, x::Real)` for inverse.
"""
Base.convert{R<:Real}(::Type{R},x::DateTime) = convert(R,value(x))
"""
    convert{T<:Real}(::Type{T}, dt::Date) -> T
Converts a Date value `dt` to a number of type `T`. The returned value corresponds to the number of Rata Die days since epoch.
See `convert(Date, x::Real)` for inverse.
"""
Base.convert{R<:Real}(::Type{R},x::Date)     = convert(R,value(x))
"""
    convert{T<:Real}(::Type{DateTime}, x::T) -> DateTime
Converts a number of type `T` to a DateTime. `x` should be the number of Rata Die milliseconds since epoch.
See `convert(Int64,dt::DateTime)` for inverse.
"""
Base.convert{R<:Real}(::Type{DateTime}, x::R) = DateTime(UTM(x))
"""
    convert{T<:Real}(::Type{Date}, x::T) -> Date
Converts a number of type `T` to a Date. `x` should be the number of Rata Die days since epoch.
See `convert(Int64,dt::Date)` for inverse.
"""
Base.convert{R<:Real}(::Type{Date}, x::R) = Date(UTD(x))

### External Conversions
const UNIXEPOCH = value(DateTime(1970)) #Rata Die milliseconds for 1970-01-01T00:00:00

"""
    unix2datetime(x) -> DateTime

Takes the number of seconds since unix epoch `1970-01-01T00:00:00` and converts to the
corresponding `DateTime`.
"""
function unix2datetime(x)
    rata = UNIXEPOCH + round(Int64, Int64(1000) * x)
    return DateTime(UTM(rata))
end
"""
    datetime2unix(dt::DateTime) -> Float64

Takes the given `DateTime` and returns the number of seconds
since the unix epoch `1970-01-01T00:00:00` as a `Float64`.
"""
datetime2unix(dt::DateTime) = (value(dt) - UNIXEPOCH)/1000.0

"""
    now() -> DateTime

Returns a `DateTime` corresponding to the user's system time including the system timezone
locale.
"""
function now()
    tv = Libc.TimeVal()
    tm = Libc.TmStruct(tv.sec)
    return DateTime(tm.year+1900,tm.month+1,tm.mday,tm.hour,tm.min,tm.sec,div(tv.usec,1000))
end

"""
    today() -> Date

Returns the date portion of `now()`.
"""
today() = Date(now())

"""
    now(::Type{UTC}) -> DateTime

Returns a `DateTime` corresponding to the user's system time as UTC/GMT.
"""
now(::Type{UTC}) = unix2datetime(time())

"""
    rata2datetime(days) -> DateTime

Takes the number of Rata Die days since epoch `0000-12-31T00:00:00` and returns the
corresponding `DateTime`.
"""
rata2datetime(days) = DateTime(yearmonthday(days)...)

"""
    datetime2rata(dt::TimeType) -> Int64

Returns the number of Rata Die days since epoch from the given `Date` or `DateTime`.
"""
datetime2rata(dt::TimeType) = days(dt)

# Julian conversions
const JULIANEPOCH = value(DateTime(-4713,11,24,12))

"""
    julian2datetime(julian_days) -> DateTime

Takes the number of Julian calendar days since epoch `-4713-11-24T12:00:00` and returns the
corresponding `DateTime`.
"""
function julian2datetime(f)
    rata = JULIANEPOCH + round(Int64, Int64(86400000) * f)
    return DateTime(UTM(rata))
end

"""
    datetime2julian(dt::DateTime) -> Float64

Takes the given `DateTime` and returns the number of Julian calendar days since the julian
epoch `-4713-11-24T12:00:00` as a `Float64`.
"""
datetime2julian(dt::DateTime) = (value(dt) - JULIANEPOCH)/86400000.0

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using CirculantCov
import CirculantCov as CC
using Test

@testset "φ_grid, θ_grid" begin

    θspan = (2.3,2.7)
    N     = 250
    type  = :equicosθ
    θ, θ∂ = CC.θ_grid(; θspan, N, type)
    @test all(θspan[1] .<  θ  .<  θspan[2])
    @test all(θspan[1] .<= θ∂ .<= θspan[2])
    @test length(θ)  == N
    @test length(θ∂) == N+1

    θspan = (2.3,2.7)
    N     = 250
    type  = :equiθ
    θ, θ∂ = CC.θ_grid(; θspan, N, type)
    @test all(θspan[1] .<  θ  .<  θspan[2])
    @test all(θspan[1] .<= θ∂ .<= θspan[2])
    @test length(θ)  == N
    @test length(θ∂) == N+1

end

@testset "counterclock_Δφ and fullcircle" begin

    φstart, φstop = 0, 1.0        ; @test CC.counterclock_Δφ(φstart, φstop) ≈ φstop - φstart
    φstart, φstop = 0, π          ; @test CC.counterclock_Δφ(φstart, φstop) ≈ φstop - φstart
    φstart, φstop = 0, 3π/2       ; @test CC.counterclock_Δφ(φstart, φstop) ≈ φstop - φstart
    φstart, φstop = 0, 2pi - .001 ; @test CC.counterclock_Δφ(φstart, φstop) ≈ φstop - φstart

    φstart, φstop = 2π - 0.5, 1.5 ; @test CC.counterclock_Δφ(φstart, φstop) ≈ 2.0

    φstart, φstop =  1.0        , 0 ; @test isapprox(CC.counterclock_Δφ(φstart, φstop), 2pi - 1.0, rtol=1e-5)
    φstart, φstop =  π          , 0 ; @test isapprox(CC.counterclock_Δφ(φstart, φstop), π, rtol=1e-5)
    φstart, φstop =  3π/2       , 0 ; @test isapprox(CC.counterclock_Δφ(φstart, φstop), π/2, rtol=1e-5)
    φstart, φstop =  2pi - .001 , 0 ; @test isapprox(CC.counterclock_Δφ(φstart, φstop), .001, rtol=1e-5)

    @test CC.fullcircle(2pi * (0:9) / 10)[2]   == 1
    @test CC.fullcircle(pi/2 * (0:9) / 10)[2]  == 4
    @test CC.fullcircle(2pi/3 * (0:9) / 10)[2] == 3
    @test all(CC.fullcircle(2pi * (0:9) / 10)[1]  .≈ 2pi * (0:9) / 10)
    @test all(CC.fullcircle(pi/2 * (0:9) / 10)[1]  .≈ 2pi * (0:39) / 40)
    @test all(CC.fullcircle(2pi/3 * (0:9) / 10)[1]  .≈ 2pi * (0:29) / 30)

end

# @testset "CirculantCov: periodize" begin

    #= 3 ways to test out periodizing

    n = 1024*3    # even
    #n = (768-1)*3 # odd
    f = 2π*(0:n-1)/n .|> cos .|> cφ -> exp(-(1-cφ)*100) * (1 - ((1-cφ)*100)^2/2) .|> complex

    freq_mult=3
    nfm = n÷freq_mult
    @assert nfm == n//freq_mult

    # -----
    f′ = sum( circshift(f, k*nfm) for k=0:freq_mult-1)
    f  .|> real  |> plot
    f′ .|> real |> plot
    f′[1:nfm] .|> real |> plot

    # -----
    fk = zeros(eltype(f), size(f))
    fk[1:freq_mult:end] = fft(f)[1:freq_mult:end]
    f′′ = freq_mult.* ifft(fk)

    f .|> real |> plot
    f′′ .|> real |> plot
    f′′[1:nfm] .|> real |> plot

    # -----
    f′′′ = ifft(fft(f)[1:freq_mult:end])
    f    .|> real |> plot
    f′′′ .|> real |> plot
    =#

# end


text:
using IBMQJulia
using SimpleMock
using Yao, YaoBlocks
using Test

# todo add tests wrt standard circuits(qft etc) to confirm that qobj is correct
# verify the results
qc = chain(1, put(1=>I2))
test_token = "8e87a83bbe4f5ad0aa953094fb8df853b07b2a86dadf010261eef2a65cd524df29e5bc38bbe2f8064155226bdfb6dcaecc8e3a6a029e402c8c7389ec0cef3574"
user = authenticate(test_token)

mock(readline => Mock(() -> "1")) do _
    global reg = createreg(user)
end
job = apply!(reg, [qc])
sleep(10)
stat = status(job)
if stat == "COMPLETED"
    res = getresult(job)
end

@testset "IBMQJulia.jl" begin
    @testset "API" begin
        include("testapi.jl")
    end

    @testset "Qobj" begin
        include("testqobj.jl")
    end

    @testset "Single Quibit Unitary Gates" begin
        @test isunitary(U1(0.5))
        @test isunitary(U2(0.5, 0.6))
        @test isunitary(U3(0.5, 0.6, 0.4))
        u1 = U1(π/2)
        u2 = U2(π/2, π/6)
        u3 = U3(π/2, π/6, π/4)

        for fs in [u1, u2, u3]
            @test eval(YaoBlocks.parse_ex(dump_gate(fs), 1)) == fs
            @test Yao.iparams_eltype(fs) == Float64
        end
        @test Yao.getiparams(u1) == (π/2)
        @test Yao.getiparams(u2) == (π/2, π/6)
        @test Yao.getiparams(u3) == (π/2, π/6, π/4)

    end
end


text:
"""
Shows the differentation backends available for use with `DerivableFunctions.jl`.
"""
diff_backends() = vcat([:Symbolic, :ForwardDiff], AddedBackEnds(Val(length(methods(AddedBackEnds)) > 1)))
AddedBackEnds(::Val{false}) = Symbol[]

GetDeriv(ADmode::Symbol, args...; kwargs...) = GetDeriv(Val(ADmode), args...; kwargs...)
GetGrad(ADmode::Symbol, args...; kwargs...) = GetGrad(Val(ADmode), args...; kwargs...)
GetJac(ADmode::Symbol, args...; kwargs...) = GetJac(Val(ADmode), args...; kwargs...)
GetHess(ADmode::Symbol, args...; kwargs...) = GetHess(Val(ADmode), args...; kwargs...)
GetDoubleJac(ADmode::Symbol, args...; kwargs...) = GetDoubleJac(Val(ADmode), args...; kwargs...)
GetMatrixJac(ADmode::Symbol, args...; kwargs...) = GetMatrixJac(Val(ADmode), args...; kwargs...)

# ForwardDiff as standard
GetDeriv(F::Function, args...; kwargs...) = GetDeriv(Val(:ForwardDiff), F, args...; kwargs...)
GetGrad(F::Function, args...; kwargs...) = GetGrad(Val(:ForwardDiff), F, args...; kwargs...)
GetJac(F::Function, args...; kwargs...) = GetJac(Val(:ForwardDiff), F, args...; kwargs...)
GetHess(F::Function, args...; kwargs...) = GetHess(Val(:ForwardDiff), F, args...; kwargs...)
GetDoubleJac(F::Function, args...; kwargs...) = GetDoubleJac(Val(:ForwardDiff), F, args...; kwargs...)
GetMatrixJac(F::Function, args...; kwargs...) = GetMatrixJac(Val(:ForwardDiff), F, args...; kwargs...)

"""
    GetDeriv(ADmode::Val, F::Function; kwargs...) -> Function
Returns a function which computes the scalar derivative of `F` out-of-place via a backend specified by `ADmode`.

Example:
```julia
Derivative = GetDeriv(Val(:ForwardDiff), x->exp(-x^2))
Derivative(5.0)
```
For available backends, see `diff_backends()`.
"""
function GetDeriv(ADmode::Val, F::Function, args...; kwargs...)
    Dv = _GetDeriv(ADmode; kwargs...)
    EvaluateDeriv(X::Number, args...) = Dv(F, X, args...)
    EvaluateDeriv(X::Num) = _GetDerivPass(F, X)
end
function GetDeriv(ADmode::Val{:Symbolic}, F::Function, args...; verbose::Bool=true, kwargs...)
    M = try GetSymbolicDerivative(F, 1, :derivative; kwargs...)   catch;  nothing  end
    if isnothing(M)
        verbose && (@warn "Unable to compute symbolic derivative of $F, falling back to ForwardDiff.")
        GetDeriv(Val(:ForwardDiff), F)
    else M end
end
GetDeriv(ADmode::Val; Kwargs...) = (Dv=_GetDeriv(ADmode; Kwargs...);    EvaluateDerivative(F::Function, X, args...; kwargs...) = Dv(F, X, args...; kwargs...))
GetDeriv(ADmode::Val, F::DFunction, args...; Kwargs...) = EvaldF(F)

"""
    GetGrad(ADmode::Val, F::Function; kwargs...) -> Function
Returns a function which computes the gradient of `F` out-of-place via a backend specified by `ADmode`.

Example:
```julia
Gradient = GetGrad(Val(:ForwardDiff), x->x[1]^2 - x[2]^3)
Gradient(rand(2))
```
For available backends, see `diff_backends()`.
"""
function GetGrad(ADmode::Val, F::Function, args...; kwargs...)
    Dv = _GetGrad(ADmode; kwargs...)
    EvaluateGradient(X::AbstractVector{<:Number}, args...) = Dv(F, X, args...)
    EvaluateGradient(X::AbstractVector{<:Num}) = _GetGradPass(F, X)
end
function GetGrad(ADmode::Val{:Symbolic}, F::Function, m::Int=GetArgLength(F), args...; verbose::Bool=true, kwargs...)
    M = try GetSymbolicDerivative(F, m, :gradient; kwargs...)   catch;  nothing  end
    if isnothing(M)
        verbose && (@warn "Unable to compute symbolic derivative of $F, falling back to ForwardDiff.")
        GetGrad(Val(:ForwardDiff), F, m, args...)
    else M end
end
GetGrad(ADmode::Val; Kwargs...) = (Dv=_GetGrad(ADmode; Kwargs...);  EvaluateGradient(F::Function, X, args...; kwargs...) = Dv(F, X, args...; kwargs...))
GetGrad(ADmode::Val, F::DFunction, args...; Kwargs...) = EvaldF(F)

"""
    GetJac(ADmode::Val, F::Function; kwargs...) -> Function
Returns a function which computes the Jacobian of `F` out-of-place via a backend specified by `ADmode`.

Example:
```julia
Jacobian = GetJac(Val(:ForwardDiff), x->[x[1]^2, -x[2]^3, x[1]*x[2]])
Jacobian(rand(2))
```
For available backends, see `diff_backends()`.
"""
function GetJac(ADmode::Val, F::Function, args...; kwargs...)
    Dv = _GetJac(ADmode; kwargs...)
    EvaluateJacobian(X::AbstractVector{<:Number}, args...) = Dv(F, X, args...)
    EvaluateJacobian(X::AbstractVector{<:Num}) = _GetJacPass(F, X)
end
function GetJac(ADmode::Val{:Symbolic}, F::Function, m::Int=GetArgLength(F), args...; verbose::Bool=true, kwargs...)
    M = try GetSymbolicDerivative(F, m, :jacobian; kwargs...)   catch;  nothing  end
    if isnothing(M)
        verbose && (@warn "Unable to compute symbolic derivative of $F, falling back to ForwardDiff.")
        GetJac(Val(:ForwardDiff), F, m, args...)
    else M end
end
GetJac(ADmode::Val; Kwargs...) = (Dv=_GetJac(ADmode; Kwargs...);  EvaluateJacobian(F::Function, X, args...; kwargs...) = Dv(F, X, args...; kwargs...))
GetJac(ADmode::Val, F::DFunction, args...; Kwargs...) = EvaldF(F)

"""
    GetHess(ADmode::Val, F::Function; kwargs...) -> Function
Returns a function which computes the Hessian of `F` out-of-place via a backend specified by `ADmode`.

Example:
```julia
Hessian = GetHess(Val(:ForwardDiff), x->x[1]^2 -x[2]^3 + x[1]*x[2])
Hessian(rand(2))
```
For available backends, see `diff_backends()`.
"""
function GetHess(ADmode::Val, F::Function, args...; kwargs...)
    Dv = _GetHess(ADmode; kwargs...)
    EvaluateHess(X::AbstractVector{<:Number}, args...) = Dv(F, X, args...)
    EvaluateHess(X::AbstractVector{<:Num}) = _GetHessPass(F, X)
end
function GetHess(ADmode::Val{:Symbolic}, F::Function, m::Int=GetArgLength(F), args...; verbose::Bool=true, kwargs...)
    M = try GetSymbolicDerivative(F, m, :hessian; kwargs...)   catch;  nothing  end
    if isnothing(M)
        verbose && (@warn "Unable to compute symbolic derivative of $F, falling back to ForwardDiff.")
        GetHess(Val(:ForwardDiff), F, m, args...)
    else M end
end
GetHess(ADmode::Val; Kwargs...) = (Dv=_GetHess(ADmode; Kwargs...);  EvaluateHessian(F::Function, X, args...; kwargs...) = Dv(F, X, args...; kwargs...))
GetHess(ADmode::Val, F::DFunction, args...; Kwargs...) = EvalddF(F)

_MakeTuple(Tup::Int) = (Tup,);    _MakeTuple(Tup::Tuple) = Tup
function _SizeTuple(F::Function, m::Int)
    T = try size(F(rand(m))) catch; size(F(rand())) end
    _MakeTuple(T)
end
"""
    GetMatrixJac(ADmode::Val, F::Function; kwargs...) -> Function
Returns a function which computes the Jacobian of an array-valued function `F` out-of-place via a backend specified by `ADmode`.

Example:
```julia
Jacobian = GetMatrixJac(Val(:ForwardDiff), x->[x[1]^2 x[2]^3; x[1]*x[2] 2])
Jacobian(rand(2))
```
For available backends, see `diff_backends()`.
"""
function GetMatrixJac(ADmode::Val, F::Function, m::Int=GetArgLength(F), f::Tuple=_SizeTuple(F,m), args...; kwargs...)
    Dv = _GetJac(ADmode; kwargs...)
    EvaluateMatrixJacobian(X::AbstractVector{<:Number}, args...) = reshape(Dv(vec∘F, X, args...), f..., m)
    EvaluateMatrixJacobian(X::Number) = reshape(Dv(vec∘F∘(z::AbstractVector->z[1]), [X]), f..., m)
    EvaluateMatrixJacobian(X::Union{<:Num,<:AbstractVector{<:Num}}) = _GetMatrixJacPass(F, X)
end
function GetMatrixJac(ADmode::Val{:Symbolic}, F::Function, m::Int=GetArgLength(F), f::Tuple=_SizeTuple(F,m), args...; verbose::Bool=true, kwargs...)
    M = try GetSymbolicDerivative(F, m, :matrixjacobian; kwargs...)   catch;  nothing  end
    if isnothing(M)
        verbose && (@warn "Unable to compute symbolic derivative of $F, falling back to ForwardDiff.")
        GetMatrixJac(Val(:ForwardDiff), F, m, f, args...)
    else M end
end
GetMatrixJac(ADmode::Val; Kwargs...) = (Dv=_GetMatrixJac(ADmode; Kwargs...);    EvaluateMatrixJacobian(F::Function, X, args...; kwargs...) = Dv(F, X, args...; kwargs...))
GetMatrixJac(ADmode::Val, F::DFunction, args...; Kwargs...) = EvaldF(F)

# For emergencies: needs an extra evaluation of function to determine length(Func(p))
function _GetMatrixJac(ADmode::Val; kwargs...)
    Dv = _GetJac(ADmode; kwargs...)
    Functor(Func::Function, X::AbstractVector{<:Number}, args...) = reshape(Dv(vec∘Func, X, args...), size(Func(X))..., length(X))
    Functor(Func::Function, X::Number, args...) = reshape(Dv(vec∘Func∘(z::AbstractVector->z[1]), [X], args...), size(Func(X))..., 1)
end

"""
    GetDoubleJac(ADmode::Val, F::Function; kwargs...) -> Function
Returns a function which computes the Jacobian of the Jacobian for a vector-valued function `F` out-of-place via a backend specified by `ADmode`.

Example:
```julia
DoubleJacobian = GetDoubleJac(Val(:ForwardDiff), x->[x[1]^2, x[1]*x[2]^3])
DoubleJacobian(rand(2))
```
For available backends, see `diff_backends()`.
"""
function GetDoubleJac(ADmode::Val, F::Function, m::Int=GetArgLength(F), f::Int=length(F(rand(m))), args...; kwargs...)
    Dv = _GetJac(ADmode; kwargs...)
    if f == 1
        EvaluateDoubleJac(X::AbstractVector{<:Number}) = reshape(Dv(vec∘(z->Dv(F,z)), X), m, m)
        EvaluateDoubleJac(X::AbstractVector{<:Num}) = _GetDoubleJacPass(F, X)
    else
        EvaluateDoubleJacobian(X::AbstractVector{<:Number}) = reshape(Dv(vec∘(z->Dv(F,z)), X), f, m, m)
        EvaluateDoubleJacobian(X::AbstractVector{<:Num}) = _GetDoubleJacPass(F, X)
    end
end
function GetDoubleJac(ADmode::Val{:Symbolic}, F::Function, m::Int=GetArgLength(F), f::Int=length(F(rand(m))), args...; verbose::Bool=true, kwargs...)
    M = try GetSymbolicDerivative(F, m, :doublejacobian; kwargs...)   catch;  nothing  end
    if isnothing(M)
        verbose && (@warn "Unable to compute symbolic derivative of $F, falling back to ForwardDiff.")
        GetDoubleJac(Val(:ForwardDiff), F, m, f, args...)
    else M end
end
GetDoubleJac(ADmode::Val; Kwargs...) = (Dv=_GetDoubleJac(ADmode; Kwargs...);    EvaluateDoubleJacobian(F::Function, X, args...; kwargs...) = Dv(F, X, args...; kwargs...))
GetDoubleJac(ADmode::Val, F::DFunction, args...; Kwargs...) = EvalddF(F)

# For emergencies: needs an extra evaluation of function to determine length(Func(p))
function _GetDoubleJac(ADmode::Val; kwargs...)
    Dv = _GetJac(ADmode; kwargs...)
    Functor(Func::Function, X) = reshape(Dv(vec∘(z->Dv(Func,z)), X), length(Func(X)), length(X), length(X))
end

# Fall back to ForwardDiff as standard
_GetDeriv(ADmode::Val{true}; kwargs...) = _GetDeriv(Val(:ForwardDiff); kwargs...)
_GetGrad(ADmode::Val{true}; kwargs...) = _GetGrad(Val(:ForwardDiff); kwargs...)
_GetJac(ADmode::Val{true}; kwargs...) = _GetJac(Val(:ForwardDiff); kwargs...)
_GetHess(ADmode::Val{true}; kwargs...) = _GetHess(Val(:ForwardDiff); kwargs...)
_GetMatrixJac(ADmode::Val{true}; kwargs...) = _GetMatrixJac(Val(:ForwardDiff); kwargs...)
_GetDoubleJac(ADmode::Val{true}; kwargs...) = _GetDoubleJac(Val(:ForwardDiff); kwargs...)

# User has passed either Num or Vector{Num} to function, try to perfom symbolic passthrough
_GetDerivPass(F::Function, X) = SymbolicPassthrough(F(X), X, :derivative)
_GetGradPass(F::Function, X) = SymbolicPassthrough(F(X), X, :gradient)
_GetJacPass(F::Function, X) = SymbolicPassthrough(F(X), X, :jacobian)
_GetHessPass(F::Function, X) = SymbolicPassthrough(F(X), X, :hessian)
_GetMatrixJacPass(F::Function, X) = SymbolicPassthrough(F(X), X, :matrixjacobian)
_GetDoubleJacPass(F::Function, X) = SymbolicPassthrough(F(X), X, :doublejacobian)

## Implement this for new backends
_GetDeriv(ADmode::Val{:ForwardDiff}; kwargs...) = ForwardDiff.derivative
_GetGrad(ADmode::Val{:ForwardDiff}; kwargs...) = ForwardDiff.gradient
_GetJac(ADmode::Val{:ForwardDiff}; kwargs...) = ForwardDiff.jacobian
_GetHess(ADmode::Val{:ForwardDiff}; kwargs...) = ForwardDiff.hessian
##

# Error messages for unloaded backends
_GetDeriv(ADmode::Val{T}; kwargs...) where T = throw("Backend $T does not exist or is not currently loaded.")
_GetGrad(ADmode::Val{T}; kwargs...) where T = throw("Backend $T does not exist or is not currently loaded.")
_GetJac(ADmode::Val{T}; kwargs...) where T = throw("Backend $T does not exist or is not currently loaded.")
_GetHess(ADmode::Val{T}; kwargs...) where T = throw("Backend $T does not exist or is not currently loaded.")
# _GetMatrixJac and _GetDoubleJac already implemented in terms of _GetJac

GetGrad!(ADmode::Symbol, args...; kwargs...) = GetGrad!(Val(ADmode), args...; kwargs...)
GetJac!(ADmode::Symbol, args...; kwargs...) = GetJac!(Val(ADmode), args...; kwargs...)
GetHess!(ADmode::Symbol, args...; kwargs...) = GetHess!(Val(ADmode), args...; kwargs...)
GetMatrixJac!(ADmode::Symbol, args...; kwargs...) = GetMatrixJac!(Val(ADmode), args...; kwargs...)

# No Passthrough for these
GetGrad!(ADmode::Val, args...; kwargs...) = _GetGrad!(ADmode, args...; kwargs...)
GetJac!(ADmode::Val, args...; kwargs...) = _GetJac!(ADmode, args...; kwargs...)
GetHess!(ADmode::Val, args...; kwargs...) = _GetHess!(ADmode, args...; kwargs...)
GetMatrixJac!(ADmode::Val, args...; kwargs...) = _GetMatrixJac!(ADmode, args...; kwargs...)

# ForwardDiff as standard
GetGrad!(F::Function, args...; kwargs...) = GetGrad!(Val(:ForwardDiff), F, args...; kwargs...)
GetJac!(F::Function, args...; kwargs...) = GetJac!(Val(:ForwardDiff), F, args...; kwargs...)
GetHess!(F::Function, args...; kwargs...) = GetHess!(Val(:ForwardDiff), F, args...; kwargs...)
GetMatrixJac!(F::Function, args...; kwargs...) = GetMatrixJac!(Val(:ForwardDiff), F, args...; kwargs...)

# Evaluation of differentation operations into pre-specified arrays for functions which are themselves out-of-place
"""
    GetGrad!(ADmode::Val, F::Function; kwargs...) -> Function
Returns a function which computes gradients in-place via a backend specified by `ADmode`.
The function returned by `GetGrad!` has argument structure `(Y::AbstractVector, X::AbstractVector)` where the gradient of `F` evaluated at `X` is saved into `Y`.

Example:
```julia
Gradient! = GetGrad!(Val(:ForwardDiff), x->x[1]^2 - x[2]^3)
Y = Vector{Float64}(undef, 2)
Gradient!(Y, rand(2))
```
For available backends, see `diff_backends()`.
"""
function GetGrad!(ADmode::Val, F::Function; kwargs...)
    Dv = _GetGrad!(ADmode; kwargs...)
    EvaluateGradient!(Y::AbstractVector{<:Number}, X::AbstractVector{<:Number}, args...) = Dv(Y, F, X, args...)
    EvaluateGradient!(Y::AbstractVector{<:Num}, X::AbstractVector{<:Num}) = _GetGradPass!(Y, F, X)
end

"""
    GetJac!(ADmode::Val, F::Function; kwargs...) -> Function
Returns a function which computes Jacobians in-place via a backend specified by `ADmode`.
The function returned by `GetJac!` has argument structure `(Y::AbstractMatrix, X::AbstractVector)` where the Jacobian of `F` evaluated at `X` is saved into `Y`.

Example:
```julia
Jacobian! = GetJac!(Val(:ForwardDiff), x->[x[1]^2, -x[2]^3, x[1]*x[2]])
Y = Matrix{Float64}(undef, 3, 2)
Jacobian!(Y, rand(2))
```
For available backends, see `diff_backends()`.
"""
function GetJac!(ADmode::Val, F::Function; kwargs...)
    Dv = _GetJac!(ADmode; kwargs...)
    EvaluateJacobian!(Y::AbstractMatrix{<:Number}, X::AbstractVector{<:Number}, args...) = Dv(Y, F, X, args...)
    EvaluateJacobian!(Y::AbstractMatrix{<:Num}, X::AbstractVector{<:Num}) = _GetJacPass!(Y, F, X)
end

"""
    GetHess!(ADmode::Val, F::Function; kwargs...) -> Function
Returns a function which computes Hessians in-place via a backend specified by `ADmode`.
The function returned by `GetHess!` has argument structure `(Y::AbstractMatrix, X::AbstractVector)` where the Hessian of `F` evaluated at `X` is saved into `Y`.

Example:
```julia
Hessian! = GetHess!(Val(:ForwardDiff), x->x[1]^2 -x[2]^3 + x[1]*x[2])
Y = Matrix{Float64}(undef, 2, 2)
Hessian!(Y, rand(2))
```
For available backends, see `diff_backends()`.
"""
function GetHess!(ADmode::Val, F::Function; kwargs...)
    Dv = _GetHess!(ADmode; kwargs...)
    EvaluateHess!(Y::AbstractMatrix{<:Number}, X::AbstractVector{<:Number}, args...) = Dv(Y, F, X, args...)
    EvaluateHess!(Y::AbstractMatrix{<:Num}, X::AbstractVector{<:Num}) = _GetHessPass!(Y, F, X)
end

"""
    GetMatrixJac!(ADmode::Val, F::Function; kwargs...) -> Function
Returns a function which computes Jacobians in-place for array-valued functions via a backend specified by `ADmode`.
The function returned by `GetMatrixJac!` has argument structure `(Y::AbstractArray, X::AbstractVector)` where the Jacobian of `F` evaluated at `X` is saved into `Y`.

Example:
```julia
Jacobian! = GetMatrixJac!(Val(:ForwardDiff), x->[x[1]^2 x[2]^3; x[1]*x[2] 2])
Y = Array{Float64}(undef, 2, 2, 2)
Jacobian!(Y, rand(2))
```
For available backends, see `diff_backends()`.
"""
function GetMatrixJac!(ADmode::Val, F::Function; kwargs...)
    Dv = _GetMatrixJac!(ADmode; kwargs...)
    EvaluateMatrixJacobian(Y::AbstractArray{<:Number}, X::AbstractVector{<:Number}, args...) = Dv(Y, F, X, args...)
    EvaluateMatrixJacobian(Y::AbstractArray{<:Num}, X::AbstractVector{<:Num}) = _GetMatrixJacPass!(Y, F, X)
end

# Need to extend this to functions F which are themselves also in-place
_GetGradPass!(Y, F::Function, X) = copyto!(Y, SymbolicPassthrough(F(X), X, :gradient))
_GetJacPass!(Y, F::Function, X) = copyto!(Y, SymbolicPassthrough(F(X), X, :jacobian))
_GetHessPass!(Y, F::Function, X) = copyto!(Y, SymbolicPassthrough(F(X), X, :hessian))
_GetMatrixJacPass!(Y, F::Function, X) = copyto!(Y, SymbolicPassthrough(F(X), X, :matrixjacobian))

_GetGrad!(ADmode::Val{T}; kwargs...) where T = throw("Backend $T does not exist or is not currently loaded.")
_GetJac!(ADmode::Val{T}; kwargs...) where T = throw("Backend $T does not exist or is not currently loaded.")
_GetHess!(ADmode::Val{T}; kwargs...) where T = throw("Backend $T does not exist or is not currently loaded.")
_GetMatrixJac!(ADmode::Val{T}; kwargs...) where T = throw("Backend $T does not exist or is not currently loaded.")

# Fall back to ForwardDiff as standard
_GetGrad!(ADmode::Val{true}; kwargs...) = _GetGrad!(Val(:ForwardDiff); kwargs...)
_GetJac!(ADmode::Val{true}; kwargs...) = _GetJac!(Val(:ForwardDiff); kwargs...)
_GetHess!(ADmode::Val{true}; kwargs...) = _GetHess!(Val(:ForwardDiff); kwargs...)
_GetMatrixJac!(ADmode::Val{true}; kwargs...) = _GetMatrixJac!(Val(:ForwardDiff); kwargs...)

## Implement this for new backends
_GetGrad!(ADmode::Val{:ForwardDiff}; kwargs...) = ForwardDiff.gradient!
_GetJac!(ADmode::Val{:ForwardDiff}; kwargs...) = ForwardDiff.jacobian!
_GetHess!(ADmode::Val{:ForwardDiff}; kwargs...) = ForwardDiff.hessian!
_GetMatrixJac!(ADmode::Val{:ForwardDiff}; kwargs...) = _GetJac!(ADmode; kwargs...) # DELIBERATE!!!! _GetJac!() recognizes output format from given Array
##

_ConsistencyCheck(Fexpr, var, deriv::Symbol) = _ConsistencyCheck(Fexpr, var, Val(deriv))
_ConsistencyCheck(Fexpr::AbstractVector{<:Num}, var::AbstractVector{<:Num}, ::Union{Val{:jacobian},Val{:doublejacobian}}) = nothing
_ConsistencyCheck(Fexpr::Num, var::AbstractVector{<:Num}, ::Union{Val{:gradient},Val{:hessian}}) = nothing
_ConsistencyCheck(Fexpr::Num, var::Num, ::Val{:derivative}) = nothing
_ConsistencyCheck(Fexpr::AbstractArray{<:Num}, var::Union{<:Num,<:AbstractVector{<:Num}}, ::Val{:matrixjacobian}) = nothing
function _ConsistencyCheck(Fexpr, var, deriv::Val{T}) where T
    if T ∉ [:derivative, :gradient, :jacobian, :hessian, :doublejacobian, :matrixjacobian]
        throw("Invalid deriv type: $T.")
    else
        throw("Got :$T but Fexpr=$(typeof(Fexpr)) and arg=$(typeof(var)).")
    end
end

"""
Executes symbolic derivative as specified by `deriv::Symbol`.
"""
function SymbolicPassthrough(Fexpr::Union{<:AbstractArray{<:Num},<:Num}, var::Union{<:AbstractVector{<:Num},<:Num}, deriv::Symbol=:jacobian; simplify::Bool=true)
    _ConsistencyCheck(Fexpr, var, deriv)

    SymbolicDoubleJacobian(V::AbstractVector{<:Num}, z::Num; simplify::Bool=true) = SymbolicDoubleJacobian(V, [z]; simplify=simplify)
    SymbolicDoubleJacobian(V::AbstractVector{<:Num}, z::AbstractVector{<:Num}; simplify::Bool=true) = SymbolicMatrixJacobian(Symbolics.jacobian(V,z),z; simplify=simplify)
    SymbolicMatrixJacobian(M::AbstractArray{<:Num}, z::Num; simplify::Bool=true) = SymbolicMatrixJacobian(M, [z]; simplify=simplify)
    function SymbolicMatrixJacobian(M::AbstractArray{<:Num}, z::AbstractVector{<:Num}; simplify::Bool=true)
        reshape(Symbolics.jacobian(vec(M), z; simplify=simplify), size(M)..., length(z))
    end

    if deriv == :doublejacobian
        SymbolicDoubleJacobian(Fexpr, var; simplify=simplify)
    elseif deriv == :matrixjacobian
        SymbolicMatrixJacobian(Fexpr, var; simplify=simplify)
    else
        (@eval Symbolics.$deriv)(Fexpr, var; simplify=simplify)
    end
end
"""
    GetSymbolicDerivative(F::Function, inputdim::Int=GetArgLength(F), deriv::Symbol=:jacobian; timeout::Real=5, inplace::Bool=false, parallel::Bool=false)
Computes symbolic derivatives, including `:jacobian`, `:gradient`, `:hessian` and `:derivative` which are specified via `deriv`.
Special care has to be taken that the correct `inputdim` is specified! Silent errors may occur otherwise.
"""
function GetSymbolicDerivative(F::Function, inputdim::Int=GetArgLength(F), deriv::Symbol=:jacobian; timeout::Real=5, kwargs...)
    @assert deriv ∈ [:derivative, :gradient, :jacobian, :hessian, :doublejacobian, :matrixjacobian]
    @variables x X[1:inputdim]
    var = inputdim > 1 ? X : (try F(rand()); x catch; X end)
    Fexpr = KillAfter(F, var; timeout=timeout)
    # Warning already thrown in KillAfter
    isnothing(Fexpr) && return nothing
    GetSymbolicDerivative(Fexpr, var, deriv; kwargs...)
end
GetSymbolicDerivative(F::Function, deriv::Symbol; kwargs...) = GetSymbolicDerivative(F, GetArgLength(F), deriv; kwargs...)

function GetSymbolicDerivative(Fexpr::Union{<:AbstractArray{<:Num},<:Num}, var::Union{<:AbstractVector{<:Num},<:Num}, deriv::Symbol=:jacobian; simplify::Bool=true, inplace::Bool=false, parallel::Bool=false, kwargs...)
    derivative = SymbolicPassthrough(Fexpr, var, deriv; simplify=simplify)
    Builder(derivative, var; parallel=parallel, inplace=inplace, kwargs...)
end


text:
"""
        cocoso(decisionMat, weights, fns, lambda)

Apply CoCoSo (Combined Compromise Solution) method for a given matrix and weights.

# Arguments:
 - `decisionMat::DataFrame`: n × m matrix of objective values for n alternatives and m criteria
 - `weights::Array{Float64, 1}`: m-vector of weights that sum up to 1.0. If the sum of weights is not 1.0, it is automatically normalized.
 - `fns::Array{Function, 1}`: m-vector of functions to be applied on the columns.
 - `lambda::Float64`: joint criterion. 0<=lambda<=1, default=0.5.

# Description
cocoso() applies the CoCoSo method to rank n alterntives subject to m criteria which are supposed to be
either maximized or minimized.

# Output
- `::CoCoSoResult`: CoCoSoResult object that holds multiple outputs including scores, rankings, and best index.

# Examples
```julia-repl
julia> decmat = [3        12.5        2        120        14        3;
       5        15        3        110        38        4;
       3        13        2        120        19        3;
       4        14        2        100        31        4;
       3        15        1.5        125        40        4]
5×6 Array{Float64,2}:
 3.0  12.5  2.0  120.0  14.0  3.0
 5.0  15.0  3.0  110.0  38.0  4.0
 3.0  13.0  2.0  120.0  19.0  3.0
 4.0  14.0  2.0  100.0  31.0  4.0
 3.0  15.0  1.5  125.0  40.0  4.0

julia> df = makeDecisionMatrix(decmat)
5×6 DataFrame
 Row │ Crt1     Crt2     Crt3     Crt4     Crt5     Crt6
     │ Float64  Float64  Float64  Float64  Float64  Float64
─────┼──────────────────────────────────────────────────────
   1 │     3.0     12.5      2.0    120.0     14.0      3.0
   2 │     5.0     15.0      3.0    110.0     38.0      4.0
   3 │     3.0     13.0      2.0    120.0     19.0      3.0
   4 │     4.0     14.0      2.0    100.0     31.0      4.0
   5 │     3.0     15.0      1.5    125.0     40.0      4.0

julia> weights = [0.221, 0.159, 0.175, 0.127, 0.117, 0.201];

julia> fns = [maximum, minimum, minimum, maximum, minimum, maximum];

julia> lambda = 0.5;

julia> result = cocoso(df, weights, fns, lambda);

julia> result.scores
7-element Array{Float64,1}:
 2.0413128390265998
 2.787989783418825
 2.8823497955972495
 2.4160457689259287
 1.2986918936013303
 1.4431429073391682
 2.519094173200623

julia> result.bestIndex
3
```
# References

Yazdani, M., Zarate, P., Kazimieras Zavadskas, E. and Turskis, Z. (2019), "A combined compromise solution (CoCoSo) method for multi-criteria decision-making problems", Management Decision, Vol. 57 No. 9, pp. 2501-2519. https://doi.org/10.1108/MD-05-2017-0458

"""
function cocoso(decisionMat::DataFrame, weights::Array{Float64,1}, fns::Array{Function,1}, lambda::Float64=0.5):CoCoSoResult

    row, col = size(decisionMat)
    w = unitize(weights)
    colMax = colmaxs(decisionMat)
    colMin = colmins(decisionMat)

    A = similar(decisionMat)

    for i in 1:row
        for j in 1:col
            if fns[j] == maximum
                @inbounds A[i, j] = (decisionMat[i, j] - colMin[j]) / (colMax[j] - colMin[j])
            elseif fns[j] == minimum
                @inbounds A[i, j] = (colMax[j] - decisionMat[i, j]) / (colMax[j] - colMin[j])
            end
        end
    end

    scoreMat = similar(A)
    for i in 1:col
        scoreMat[:, i] = A[:, i].^w[i]
    end

    P = zeros(Float64, row)
    for i in 1:row
        P[i] = sum(scoreMat[i, :])
    end

    S = w * A |> rowsums

    scoreTable = [S P]

    kA = (S .+ P) ./ sum(scoreTable)

    kB = (S ./ minimum(S)) .+ (P ./ minimum(P))

    kC = ((lambda .* S) .+ ((1-lambda) .* P)) ./ ((lambda .* maximum(S)) .+ ((1-lambda) * maximum(P)))

    scores = (kA .+ kB .+ kC) ./ 3 .+ (kA .* kB .* kC) .^ (1/3)

    rankings = sortperm(scores)

    bestIndex = rankings |> last

    result = CoCoSoResult(
        decisionMat,
        w,
        scores,
        rankings,
        bestIndex
    )

    return result
end

max_stars_repo_name: bahadirfyildirim/JMcDM

max_issues_repo_name: bahadirfyildirim/JMcDM

max_forks_repo_name: bahadirfyildirim/JMcDM


text:
struct RelativisticPseudoPotential{T} <: AbstractPseudoPotential{T}
    name::String
    gst_config::Configuration{Orbital{Int}}
    Q::Int
    V₋::Vector{GaussianExpansion{T}} # For negative κ
    V₊::Vector{GaussianExpansion{T}} # For positive κ
    reference::String
end

AtomicPotentials.charge(pp::RelativisticPseudoPotential) = num_electrons(pp.gst_config)
AtomicPotentials.ground_state(pp::RelativisticPseudoPotential) = pp.gst_config

Base.show(io::IO, pp::RelativisticPseudoPotential{T}) where {T} =
    write(io, "Relativistic pseudo-potential for $(pp.name) ($(pp.gst_config)), Z = $(charge(pp))")

function κrange(pp::RelativisticPseudoPotential)
    κ₋max = length(pp.V₋)
    κ₊max = length(pp.V₊)

    sort(vcat(-1:-1:-κ₋max,1:κ₊max), lt=(a,b)->abs(a)<abs(b) || abs(a)==abs(b) && a<0)
end

function Base.show(io::IO, ::MIME"text/plain", pp::RelativisticPseudoPotential)
    show(io, pp)
    println(io)
    κ₋max = length(pp.V₋)
    κ₊max = length(pp.V₊)
    print(io, "Long-range Q = $(pp.Q), κ₋ ∈ -1:-1:$(-κ₋max)")
    κ₊max > 0 && print(io, ", κ₊ ∈ 1:$(κ₊max)")
    println(io, "\nData from \"$(pp.reference)\"")

    κs = κrange(pp)

    data = map(κs) do κ
        data = κ > 0 ? pp.V₊[κ] : pp.V₋[-κ]
        nk = length(data)
        ℓ = AtomicLevels.kappa_to_ℓ(κ)
        j = AtomicLevels.kappa_to_j(κ)
        hcat([spectroscopic_label(ℓ);repeat([""],nk-1)],
             [j;repeat([""],nk-1)],
             1:nk, data.n, data.B, data.β)
    end |> v -> vcat(v...)

    headers = ["ℓ", "j", "k", "n", "B", "β"]

    pretty_table(io, data, headers)
end

function (pp::RelativisticPseudoPotential{T})(orb::RelativisticOrbital, r::AbstractVector{T}) where T
    V = -pp.Q./r

    κ = orb.κ

    data = if κ < 0
        -κ > length(pp.V₋) && return V
        pp.V₋[-κ]
    else
        κ > length(pp.V₊) && return V
        pp.V₊[κ]
    end
    V += data(r)

    V
end

function spin_orbit_potential(pp::RelativisticPseudoPotential{T}, r::AbstractVector{T},
                              a::SpinOrbital{<:Orbital}, b::SpinOrbital{<:Orbital}) where T
    # This returns the (off-)diagonal, spin–orbit part of the
    # relativistic pseudopotential.
    s = half(1)
    ℓ,mℓa = jmⱼ(a)
    mas  = spin(a)
    ℓb,mℓb = jmⱼ(b)
    mbs = spin(b)

    V = zeros(T, length(r))

    ℓ == ℓb || return V

    κs = κrange(pp)

    for j = [ℓ+s, ℓ-s]
        κ = AtomicLevels.ℓj_to_kappa(ℓ,j)
        j ≥ abs(mℓa+mas) && κ ∈ κs || continue

        coeff = clebschgordan(
            ℓ,mℓa,s,mas,j
        )*clebschgordan(
            ℓ,mℓb,s,mbs,j
        )
        V += coeff*((κ > 0 ? pp.V₊[κ] : pp.V₋[-κ])(r))
    end

    V
end

function spin_orbit_potential(pp::RelativisticPseudoPotential{T}, r::AbstractVector{T},
                              a::SpinOrbital{<:RelativisticOrbital}, b::SpinOrbital{<:RelativisticOrbital}) where T
    # This returns the (off-)diagonal, spin–orbit part of the
    # relativistic pseudopotential.
    κ,κb = a.orb.κ,b.orb.κ

    V = zeros(T, length(r))

    a.orb.ℓ == b.orb.ℓ && a.orb.j == b.orb.j && a.m == b.m || return V

    κs = κrange(pp)

    if κ ∈ κs
        V += (κ > 0 ? pp.V₊[κ] : pp.V₋[-κ])(r)
    end

    V
end

function (pp::RelativisticPseudoPotential{T})(orb::SpinOrbital, r::AbstractVector{T}) where T
    V = -pp.Q./r

    V += spin_orbit_potential(pp, r, orb, orb)

    V
end

(pp::RelativisticPseudoPotential{T})(orb::AbstractOrbital, r::T) where T = pp(orb, [r])[1]


text:
s1 = "The quick brown fox jumps over the lazy dog α,β,γ"

# [search](http://docs.julialang.org/en/latest/stdlib/base/#Base.search) returns the first index of a char
i = search(s1, 'b')
println(i)
#> 11
# the second argument is equivalent to the second argument of split, see below

# or a range if called with another string
r = search(s1, "brown")
println(r)
#> 11:15

# string [replace](http://docs.julialang.org/en/latest/stdlib/base/#Base.replace) is done thus:
r = replace(s1, "brown", "red")
show(r); println()
#> "The quick red fox jumps over the lazy dog"

# search and replace can also take a regular expressions by preceding the string with 'r'
r = search(s1, r"b[\w]*n")
println(r)
#> 11:15

# again with a regular expression
r = replace(s1, r"b[\w]*n", "red")
show(r); println()
#> "The quick red fox jumps over the lazy dog"

# there are also functions for regular expressions that return RegexMatch types
# [match](http://julia.readthedocs.org/en/latest/stdlib/base/#Base.match) scans left to right for the first match (specified starting index optional)
r = match(r"b[\w]*n", s1)
println(r)
#> RegexMatch("brown")

# RegexMatch types have a property match that holds the matched string
show(r.match); println()
#> "brown"

# [matchall](http://julia.readthedocs.org/en/latest/stdlib/base/#Base.matchall) returns a vector with RegexMatches for each match
r = matchall(r"[\w]{4,}", s1)
println(r)
#> SubString{UTF8String}["quick","brown","jumps","over","lazy"]

# [eachmatch](http://julia.readthedocs.org/en/latest/stdlib/base/#Base.eachmatch) returns an iterator over all the matches
r = eachmatch(r"[\w]{4,}", s1)
for i in r print("\"$(i.match)\" ") end
println()
#> "quick" "brown" "jumps" "over" "lazy"

# a string can be repeated using the [repeat](http://julia.readthedocs.org/en/latest/manual/strings/#common-operations) function,
# or more succinctly with the [^ syntax](http://julia.readthedocs.org/en/latest/stdlib/base/#Base.^):
r = "hello "^3
show(r); println() #> "hello hello hello "

# the [strip](http://docs.julialang.org/en/latest/stdlib/base/#Base.strip) function works the same as python:
# e.g., with one argument it strips the outer whitespace
r = strip("hello ")
show(r); println() #> "hello"
# or with a second argument of an array of chars it strips any of them;
r = strip("hello ", ['h', ' '])
show(r); println() #> "ello"
# (note the array is of chars and not strings)

# similarly [split](http://docs.julialang.org/en/latest/stdlib/base/#Base.split) works in basically the same way as python:
r = split("hello, there,bob", ',')
show(r); println() #> ["hello"," there","bob"]
r = split("hello, there,bob", ", ")
show(r); println() #> ["hello","there,bob"]
r = split("hello, there,bob", [',', ' '], limit=0, keep=false)
show(r); println() #> ["hello","there","bob"]
# (the last two arguements are limit and include_empty, see docs)

# the opposite of split: [join](http://docs.julialang.org/en/latest/stdlib/base/#Base.join) is simply
r = join(collect(1:10), ", ")
println(r) #> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10

max_stars_repo_name: inkydragon/JuliaByExample

max_issues_repo_name: inkydragon/JuliaByExample

max_forks_repo_name: inkydragon/JuliaByExample