text:
#=--- CpuId / CpuId.jl ----------------------------------------------------=#

"""
# Module CpuId

Query information about and directly from your CPU.
"""
module CpuId

export cpuvendor, cpubrand, cpumodel, cachesize, cachelinesize,
       simdbytes, simdbits, address_size, physical_address_size,
       cpu_base_frequency, cpu_max_frequency, cpu_bus_frequency,
       has_cpu_frequencies, hypervised, hvvendor, hvversion,
       hvinfo, cpucores, cputhreads, cpucycle, cpucycle_id,
       perf_revision, perf_fix_counters, perf_fix_bits, perf_gen_counters,
       perf_gen_bits, cpuinfo, cpufeature, cpufeatures, cpufeaturedesc,
       cpufeaturetable, cpuarchitecture, cacheinclusive

using Markdown: MD, Table, parse
const MarkdownString = MD     # Rename Markdown constructors
const MarkdownTable = Table   # to avoid deprecation warning in 0.7-beta
const parse_markdown = parse

using Base: @_noinline_meta, @_inline_meta

# Particular feature flag query is also externalized due to largeness of dicts.
include("cpufeature.jl")

# The low-level assembly functions are confined to a helper module.
include("CpuInstructions.jl")

using .CpuInstructions: cpuid, rdtsc, rdtscp

"""
Helper function, tagged noinline to not have detrimental effect on performance.
"""
@noinline _throw_unsupported_leaf(leaf) =
    error("This CPU does not provide information on cpuid leaf 0x$(string(leaf, base=16, pad=8)).")

"""
    hasleaf(leaf::UInt32) ::Bool

Helper function (not exported) to test whether the CPU claims to provide the
given leaf in a `cpuid` instruction call.

Note: It appears LLVM really know its gear: If this function is inlined, and
      just-in-time compiled, then this test is eliminated completly if the
      executing machine does support this feature. Yeah!
"""
function hasleaf(leaf::UInt32) ::Bool
    @_inline_meta
    eax, ebx, ecx, edx = cpuid(leaf & 0xffff_0000)
    eax >= leaf
end

"""
    cpucycle()

Read the CPU's [Time Stamp Counter, TSC](https://en.wikipedia.org/wiki/Time_Stamp_Counter),
directly with a `rdtsc` instruction.  This counter is increased for every CPU
cycle, until reset.  This function has, when inlined, practically no overhead
and is, thus, probably the fasted method to count how many cycles the CPU has
spent working since last read.

Note, the TSC runs at a constant rate if `hasfeature(:TSCINV)==true`;
otherwise, it is tied to the current CPU clock frequency.

Hint: This function is extremely efficient when inlined into your own code.
      Convince yourself by typing `@code_native CpuId.cpucycle()`.
      To use this for benchmarking, simply subtract the results of two calls.
      The result is the actual CPU clock cycles spent, independent of the
      current (and possible non-constant) CPU clock frequency.
"""
function cpucycle end
@eval cpucycle() = $(cpufeature(TSC)) ? rdtsc() : zero(UInt64)

"""
    cpucycle_id()

Read the CPU's [Time Stamp Counter, TSC](https://en.wikipedia.org/wiki/Time_Stamp_Counter),
and executing CPU id directly with a `rdtscp` instruction.  This function is
similar to the `cpucycle()`, but uses an instruction that also allows to
detect if the code has been moved to a different executing CPU.  See also the
comments for `cpucycle()` which equally apply.
"""
function cpucycle_id end
@eval cpucycle_id() = $(cpufeature(RDTSCP)) ? rdtscp() : (zero(UInt64),zero(UInt64))

"""
    cpumodel()

Obtain the CPU model information as a Dict of pairs of
`:Family`, `:Model`, `:Stepping`, and `:CpuType`.
"""
function cpumodel() ::Dict{Symbol, UInt8}
    #  Stepping:    0: 3
    #  Model:       4: 7
    #  Family:      8:11
    #  Processor:  12:13
    #  Ext.Model:  16:19
    #  Ext.Family: 20:27
    eax, ebx, ecx, edx = cpuid(0x01)
    Dict( :Family    => UInt8(((eax & 0x0000_0f00) >> 8)
                             +((eax & 0x0ff0_0000) >> (20-4)))
        , :Model     => UInt8(((eax & 0x0000_00f0) >> 4)
                             +((eax & 0x000f_0000) >> (16-4)))
        , :Stepping  =>  UInt8(eax & 0x0000_000f)
        , :CpuType   =>  UInt8((eax & 0x0000_3000) >> 12))
end

"""
    hypervised()

Check whether the CPU reports to run a hypervisor context, that is,
whether the current process runs in a virtual machine.

A positive answer may indicate that other information reported by the CPU
is fake, such as number of physical and logical cores.  This is because
the hypervisor is free to decide which information to pass.
"""
function hypervised() ::Bool
    # alternative: 0x8000_000a, eax bit 8 set if hv present.
    eax, ebx, ecx, edx = cpuid(0x0000_0001)
    ((ecx >> 31) & one(UInt32)) != zero(UInt32)
end

"""
    hvvendorstring()

Determine the hypervisor vendor string as provided by the cpu by executing a
`cpuid` instruction.  Note, this string has a fixed length of 12 characters.
Use `hvvendor()` if you prefer getting a parsed Julia symbol.  If the CPU is
not running a hypervisor, a string of undefined content will be returned.
"""
function hvvendorstring()
    eax, ebx, ecx, edx = cpuid(0x4000_0000)
    String( reinterpret(UInt8, [ebx, ecx, edx] ) )
end

"""
    hvversion()

Get a dictionary with additional information of the running hypervisor.
The dictionary is empty if no hypervisor is detected, and only tags that are
provided by the hypervisor are included.

Note, the data available is hypervisor vendor dependent.
"""
function hvversion()

    d = Dict{Symbol,Any}()
    !hypervised() && return d

    vendor = d[:vendor] = hvvendor()

    # Signature and version info appear to be vendor specific.

    if vendor == :Microsoft

        # Specs see e.g.
        # https://msdn.microsoft.com/en-us/library/windows/hardware/dn613994(v=vs.85).aspx

        leaf = 0x4000_0001
        if hasleaf(leaf)
            eax, ebx, ecx, edx = cpuid(leaf)
            eax != 0x00 && (d[:signature] = String( reinterpret(UInt8, ( [eax, ] ))))
        end

        leaf = 0x4000_0002
        if hasleaf(leaf)
            eax, ebx, ecx, edx = cpuid(leaf)
            if eax != 0x00 && ebx != 0x00
                d[:version] = VersionNumber(
                                    Int((ebx >> 16) & 0xffff),
                                    Int( ebx        & 0xffff),
                                    Int(eax))
            end
            if ecx != 0x00 && edx != 0x00
                d[:servicepack] = VersionNumber(
                                    Int( ecx ),                     # "servicepack"
                                    Int((edx >> 24) & 0x0000_00ff), # "branch"
                                    Int( edx        & 0x00ff_ffff)) # "number"
            end
        end

        return d

    elseif vendor == :Xen

        # Xen is specified e.g. here
        # https://xenbits.xen.org/docs/unstable/hypercall/x86_64/include,public,arch-x86,cpuid.h.html

        leaf = 0x4000_0001
        if hasleaf(leaf)
            eax, ebx, ecx, edx = cpuid(leaf)
            d[:version] = VersionNumber(Int((eax >> 16) & 0xffff), Int(eax & 0xffff))
        end

        return d

    elseif vendor == :VMware

        # Specs see
        # https://kb.vmware.com/selfservice/microsites/search.do?language=en_US&cmd=displayKC&externalId=1009458
        # VmWare seems to provide frequency information, but nothing else valuable

        leaf = 0x4000_0010
        if hasleaf(leaf)
            eax, ebx, ecx, edx = cpuid(leaf)
            eax != 0x00 && (d[:tscfreq] = Int(eax))
            ebx != 0x00 && (d[:busfreq] = Int(ebx))
        end

        return d

    elseif vendor == :KVM

        return d

    end

    return d
end

"""
    hvinfo() ::MarkdownString

Generate a markdown table of all the detected/available/supported tags of a
running hypervisor.  If there is no hosting hypervisor, an empty markdown
string is returned.
"""
function hvinfo()
    d = hvversion()
    isempty(d) && return MarkdownString()

    md = "| Hypervisor | Value |\n|:------|:------|\n"

    haskey(d, :vendor) &&
        (md *= string( "| Vendor      | '", d[:vendor], "' |\n"))

    haskey(d, :signature) &&
        (md *= string( "| Signature   | '", d[:signature], "' |\n"))

    haskey(d, :version) &&
        (md *= string( "| Version     | ", d[:version], " |\n"))

    haskey(d, :servicepack) &&
        (md *= string( "| Servicepack | ", d[:servicepack], " |\n"))

    haskey(d, :tscfreq) &&
        (md *= string( "| Frequencies | TSC = ", d[:tscfreq] ÷ 1000,
                       " MHz, bus = ", d[:busfreq] ÷ 1000, " MHz |\n"))

    parse_markdown(md)
end

"""
    cpuvendorstring()

Determine the cpu vendor string as provided by the cpu by executing a
`cpuid` instruction.  Note, this string has a fixed length of 12 characters.
Use `cpuvendor()` if you prefer getting a parsed Julia symbol.
"""
function cpuvendorstring()
    eax, ebx, ecx, edx = cpuid(0x00)
    String( reinterpret(UInt8, [ebx, edx, ecx] ) )
end

"Map vendor string of type 'char[12]' provided by `cpuid, eax=0x0` to a Julia symbol."
const _cpuid_vendor_id = Dict(
    "AMDisbetter!" => :AMD,
    "AuthenticAMD" => :AMD,
    "CentaurHauls" => :Centaur,
    "CyrixInstead" => :Cyrix,
    "GenuineIntel" => :Intel,
    "Geode by NSC" => :NSC,
    "NexGenDriven" => :NexGen,
    "RiseRiseRise" => :Rise,
    "SiS SiS SiS " => :SiS,
    "TransmetaCPU" => :Transmeta,
    "GenuineTMx86" => :Transmeta,
    "UMC UMC UMC " => :UMC,
    "VIA VIA VIA " => :VIA,
    "Vortex86 SoC" => :Vortex,
    # Hypervisors
    "KVMKVMKVM\0\0\0" => :KVM,   # KVM actually appends null characters...
    "Microsoft Hv" => :Microsoft,
    "prl hyperv  " => :Parallels,
    "VMwareVMware" => :VMware,
    "XenVMMXenVMM" => :Xen,
)

"""
    cpuvendor()

Determine the cpu vendor as a Julia symbol.  In case the CPU vendor
identification is unknown `:Unknown` is returned (then also consider raising
an issue on Github).
"""
cpuvendor() = get(_cpuid_vendor_id, cpuvendorstring(), :Unknown)

"""
    hvvendor()

Determine the hypervisor vendor as a Julia symbol or `:NoHypervisor` if not
running a hypervisor. In case the hypervisor vendor identification is unknown
`:Unknown` is returned (then also consider raising an issue on Github).
"""
hvvendor() = hypervised() ?
                get(_cpuid_vendor_id, hvvendorstring(), :Unknown) :
                :NoHypervisor

"""
    cpubrand()

Determine the cpu brand as a string as provided by the CPU through executing
the `cpuid` instruction.  This function throws if no CPU brand information is
available form the CPU, which should never be the case on recent hardware.
"""
function cpubrand() ::String

    leaf = 0x8000_0004
    hasleaf(leaf) || _throw_unsupported_leaf(leaf)

    # Extract the information from leaf 0x8000_0002..0x8000_0004
    rstrip( String( reinterpret(UInt8,
                    [cpuid(0x8000_0002)..., cpuid(0x8000_0003)..., cpuid(0x8000_0004)..., 0x0000_0000] )
                  )
          , '\0')
end

"""
    cpuarchitecture()

This function tries to infer the CPU microarchitecture with a call to the
`cpuid` instruction.  For now, only Intel CPUs are suppored according to the
following table.  Others are identified as `:Unknown`.

Table C-1 of Intel's Optimization Reference Manual:

| Family_Model                     | Microarchitecture   |
| :------------------------------- | :------------------ |
| 06_4EH, 06_5EH                   | Skylake             |
| 06_3DH, 06_47H, 06_56H           | Broadwell           |
| 06_3CH, 06_45H, 06_46H, 06_3FH   | Haswell             |
| 06_3AH, 06_3EH                   | Ivy Bridge          |
| 06_2AH, 06_2DH                   | Sandy Bridge        |
| 06_25H, 06_2CH, 06_2FH           | Westmere            |
| 06_1AH, 06_1EH, 06_1FH, 06_2EH   | Nehalem             |
| 06_17H, 06_1DH                   | Enhanced Intel Core |
| 06_0FH                           | Intel Core          |
"""
function cpuarchitecture() ::Symbol

    # See also the C++ library VC++ which has quite good detection.
    # https://github.com/VcDevel/Vc/blob/master/cmake/OptimizeForArchitecture.cmake
    #
    # See also Table 35-1 in Intel's Architectures Software Developer Manual.

    cpumod = cpumodel()
    family = cpumod[:Family]
    model  = cpumod[:Model]

    family == 0x06 &&
        return (model == 0x66) ? :Cannonlake :
               (model == 0x8e || model == 0x9e) ? :Kabylake :
               (model == 0x4e || model == 0x5e || model == 0x55) ? :Skylake :
               (model == 0x3d || model == 0x47 || model == 0x4f || model == 0x56) ? :Broadwell :
               (model == 0x3c || model == 0x3f || model == 0x45 || model == 0x46) ? :Haswell :
               (model == 0x3a || model == 0x3e) ? :IvyBridge :
               (model == 0x2a || model == 0x2d) ? :SandyBridge :
               (model == 0x25 || model == 0x2c || model == 0x2f) ? :Westmere :
               (model == 0x1a || model == 0x1e || model == 0x1f || model == 0x2e) ?  :Nehalem :
               (model == 0x17 || model == 0x1d) ?  :EnhancedIntelCore :
               (model == 0x0f || model == 0x1d) ?  :IntelCore :
               # Atom models
               (model == 0x5c) ? :Goldmont :
               (model == 0x5a) ? :Silvermont :
               (model == 0x1c) ? :Atom :
               # Xeon Phi models
               (model == 0x57) ? :KnightsLanding :
               # Well, this is awkward...
               :UnknownIntel

    # Xeon Phi family 0x07, model 0x01, or Itanium ?
    family == 0x07 && return :Itanium

    # AMD types
    family == 0x0f && return :K8
    family == 0x1f && return :K10
    family == 0x2f && return :Griffin       # not confirmed
    family == 0x3f && return :Llano         # not confirmed
    family == 0x5f && return :Bobcat
    family == 0x6f && return :Bulldozer
    family == 0x7f && return (model == 0x30) ? :Puma : :Jaguar
    family == 0x8f && return :Zen

    return :Unknown
end

"""
    cachelinesize()

Query the CPU about the L1 data cache line size in bytes.  This is typically
64 byte.  Returns zero if cache line size information is not available from
the CPU.
"""
function cachelinesize() ::Int

    # If the extended leaf is not available, same information
    # should be available through leaf 0x01, ebx[15:8]

    leaf = 0x8000_0006
    hasleaf(leaf) || return 0

    eax, ebx, ecx, edx = cpuid(leaf)
    (ecx & 0xff) % Int
end

"""
    simdbytes()

Query the CPU on the maximum supported SIMD vector size in bytes, or
`sizeof(Int)` if no SIMD capability is reported by the invoked `cpuid`
instruction.
"""
function simdbytes() ::Int
    @_inline_meta

    simd = 0

    if hasleaf(0x0000_0007)
        eax, ebx, ecx, edx = cpuid(0x0000_0007)
        simd = ebx & (1<<16) != 0 ? 512 ÷ 8 :  # AVX512F instruction set
               ebx & (1<< 5) != 0 ? 256 ÷ 8 : 0   # AVX2
    end

    if simd == 0
        eax, ebx, ecx, edx = cpuid(0x0000_0001)
        simd = ecx & (1<<28) != 0 ? 256 ÷ 8 :     # AVX
            edx & (1<<26) != 0 ? 128 ÷ 8 :     # SSE2
            edx & (1<<25) != 0 ? 128 ÷ 8 :     # SSE
            edx & (1<<23) != 0 ?  64 ÷ 8 :     # MMX
            sizeof(Int)
    end

    simd
end

"""
    simdbits()

Query the CPU on the maximum supported SIMD vector size in bits, or
`sizeof(Int)` in bits if no SIMD capability is reported by the invoked `cpuid`
instruction.
"""
simdbits() = simdbytes() * 8

"""
    cpucores()

Determine the number of physical cores on the current executing CPU by
invoking a `cpuid` instruction.  On systems with multiple CPUs, this only
gives information on the single CPU that is executing the code.
Returns zero if querying this feature is not supported, which may also be due
to a running hypervisor (as observed on hvvendor() == :Microsoft).

Also, this function does not take logical cores (aka hyperthreading) into
account, but determines the true number of physical cores, which typically
also share L3 caches and main memory bandwidth.

See also the Julia global variable `Base.Sys.CPU_THREADS`, which gives the
total count of all logical cores on the machine.
"""
function cpucores() ::Int

    leaf = 0x0000_000b
    hasleaf(leaf) || return zero(UInt32)

    # The number of cores reported by cpuid is actually already the total
    # number of cores at that level, including all of the lower levels.
    # Thus, we need to find the highest level...which is 0x02 == "Core"
    # on ecx[15:08] per the specs, and divide it by the number of
    # 0x01 == "SMT" logical cores.

    sl = zero(UInt32)
    nc = zero(UInt32) # "Core" count
    nl = zero(UInt32) # "SMT" count
    while (true)
        # ebx[15:0] must be non-zero according to manual
        eax, ebx, ecx, edx = cpuid(leaf, sl)
        ebx & 0xffff == 0x0000 && break
        sl += one(UInt32)
        lt = ((ecx >> 8) & 0xff) & 0x03
        # Let's assume for now there's only one valid entry for each level
        lt == 0x01 && ( nl = ebx & 0xffff; continue )
        lt == 0x02 && ( nc = ebx & 0xffff; continue )
        # others are invalid and shouldn't be considered..
    end

    return iszero(nc) ? # no cores detected? then maybe its AMD?
        # AMD
        ((cpuid(0x8000_0008)[3] & 0x00ff)+1) :
        # Intel, we need nonzero values of nc and nl
        (iszero(nl) ? nc : nc ÷ nl)
end

"""
    cpunodes() -> Int

Determine the number of core complexes, aka nodes, on this processor.
This notion is introduced by AMD, where L3 caches are shared among the
cores of a comples
"""
function cpunodes()

    # AMD gives the nodes per processor
    # in extended topology information.
    # CPUID[0x8000_001e][ECX][10:8]
    cpunodes_amd() = 1 + ((cpuid(0x8000_001e)[3] >> 8) & 0b0111)

    return 1

end

"""
    cputhreads_per_core() -> Int

Determine the of threads per hardware core on the currently executing CPU.
A value larger than one indicates simulatenous multithreading being enabled,
aka SMT, aka Hyperthreading.
"""
function cputhreads_per_core()

    cputhreads_per_core_amd() =
        # AMD gives the threads per physical core
        # in extended topology information.
        ((cpuid(0x8000_001e)[2] >> 8) & 0x00ff) + 1

    return 1

end

"""
    cputhreads()

Determine the number of logical cores on the current executing CPU by
invoking a `cpuid` instruction.  On systems with multiple CPUs, this only
gives information on the single CPU that is executing the code.
Returns zero if querying this feature is not supported, which may also be due
to a running hypervisor (as observed on hvvendor() == :Microsoft).

In contrast to `cpucores()`, this function also takes logical cores aka
hyperthreading into account.  For practical purposes, only I/O intensive code
should make use of these total number of cores; memory or computation bound
code will not benefit, but rather experience a detrimental effect.

See also Julia's global variable `Base.Sys.CPU_THREADS`, which gives the total
count of all logical cores on the machine.  Thus, `Base.Sys.CPU_THREADS ÷
CpuId.cputhreads()` should give you the number of CPUs (packages) in your
system.
"""
function cputhreads() ::Int

    function cputhreads_amd()
        # AMD stores the total number of threads
        # aka logical processors directly
        ((cpuid(0x0000_0001)[2] >> 16) & 0x00ff)
    end

    # 1) First try to detect whether we have legacy style core count encoding
    #    This is also correct for AMD, but not for modern Intel.
    #    nc = ((cpuid(0x0000_0001)[2] >> 16) % 8)
    #    if !iszero(nc) return nc
    # 2) Try the modern intel extended information

    leaf = 0x0000_000b
    hasleaf(leaf) || return zero(UInt32)

    # The number of cores reported by cpuid is actually already the total
    # number of cores at that level, including all of the lower levels.
    # Thus, we only need to find the highest level...which is 0x02 == "Core"
    # on ecx[15:08] per the specs.

    sl = zero(UInt32)
    nc = zero(UInt32)
    while (true)
        # ebx[15:0] must be non-zero according to manual
        eax, ebx, ecx, edx = cpuid(leaf, sl)
        ebx & 0xffff == 0x0000 && break
        sl += one(UInt32)
        lt = (ecx >> 8) & 0xff  # level type, 0x01 == "SMT", 0x02 == "Core", other are invalid.
        # Let's assume for now there's only one valid entry for each level
        lt != 0x02 && continue
        nc = ebx & 0xffff
    end

    return iszero(nc) ? # no cores detected? then maybe its AMD?
        # AMD
        ((cpuid(0x0000_0001)[2] >> 16) & 0x00ff) :
        # Intel
        (nc)

end

@deprecate cpucores_total() cputhreads()

"""
    address_size()

Determine the maximum virtual address size supported by this CPU as
reported by the `cpuid` instructions.

This information may be used to determine the number of high bits that can be
used in a pointer for tagging;  viz. `sizeof(Int) - address_size() ÷ 8`,
which gives on most 64 bit Intel machines 2 bytes = 16 bit for other purposes.
"""
function address_size() ::Int

    leaf = 0x8000_0008
    hasleaf(leaf) || _throw_unsupported_leaf(leaf)

    eax, ebx, ecx, edx = cpuid(leaf)
    (eax & 0xff00) >> 8
end

"""
    physical_address_size()

Determine the maximum phyiscal addresses size supported by this CPU as
reported by the `cpuid` instructions.  Prefer to make use of `address_size`
for practical purposes; use this only for diagnostic issues, such as
determining the theoretical maximum memory size.  Also note that this address
size is manipulated by a running hypervisor.
"""
function physical_address_size() ::Int

    leaf = 0x8000_0008
    hasleaf(leaf) || _throw_unsupported_leaf(leaf)

    eax, ebx, ecx, edx = cpuid(leaf)
    eax & 0xff
end

"""
    cachesize()
    cachesize(lvl::Integer)

Obtain information on the CPU's *data* cache sizes.

Determine the data cache size for each cache level as reported by the CPU
using a set of calls to the `cpuid` instruction.  Returns a tuple with the
tuple indices matching the cache levels; sizes are given in bytes.

If given an integer, then the data cache size of the respective cache level
will be returned.  This is significantly faster than the tuple version above.

Note that these are total cache sizes, where some cache levels are typically
shared by multiple cpu cores, the higher cache levels may include lower levels.
To print the cache levels in kbyte, use e.g. `CpuId.cachesize() .÷ 1024`.
"""
function cachesize end

"""
Helper function that performs the actual computation of the cache size with
register values retrieved from `cpuid` on leaf 0x04.

> Cache size information on leaf 0x04 is computed with
>     size in bytes = (ways+1) * (partitions+1) * (linesize+1) * (sets+1)
> where
>     ways = ebx[22:31], partitions = ebx[12:21], linesize = ebx[0:11]
>     sets = ecx[:]
"""
function __datacachesize(eax::UInt32, ebx::UInt32, ecx::UInt32) ::UInt32
    @_inline_meta
    (1 + (ebx>>22) & 0x03ff ) *    # ways
    (1 + (ebx>>12) & 0x03ff ) *    # partitions
    (1 +  ebx      & 0x0fff ) *    # linesize
    (1 +  ecx )                    # sets
end

@noinline function cachesize()

    # TODO: This function fails compilation if inlined.

    function cachesize_level(leaf, sl::UInt32)
        eax, ebx, ecx, edx = cpuid(leaf, sl)
        # if eax is zero in the lowest 5 bits, we've reached the sentinel.
        eax & 0x1f == 0 && return ()
        # could do a sanity check: cache level reported in eax bits 5:7
        # if lowest bit on eax is zero, then its not a data cache
        eax & 0x01 == 0 && return cachesize_level(leaf, sl + one(UInt32))
        # otherwise this should be a valid data or shared cache level
        (signed(__datacachesize(eax, ebx, ecx)), cachesize_level(leaf, sl + one(UInt32))...)
    end

    # TODO: This is awkwardly slow and requires some rework.
    #       Potential approach: Recurse to the last found cache level, there
    #       allocate a small array, then fill the array when leaving each
    #       recursion level.

    # AMD Exteneded Cache
    leaf = 0x8000_001d
    hasleaf(leaf) && cpufeature(TOPX) &&
        return (cachesize_level(leaf,zero(UInt32))...,)
    # Intel
    leaf = 0x0000_0004
    hasleaf(leaf) && return (cachesize_level(leaf,zero(UInt32))...,)
    # no cache data available
    ()
end

cachesize(lvl::Integer) = cachesize(UInt32(lvl))

function cachesize(lvl::UInt32) ::Int
    @_inline_meta

    lvl == 0 && return 0

    leaf = 0x0000_0004
    hasleaf(leaf) || _throw_unsupported_leaf(leaf)

    # Loop over all subleaves of leaf 0x04 until the cache level bits say
    # we've reached to target level, then check whether it's a data cache or
    # shared cache.

    # Assuming cache enumeration is linear, it is sufficient to start at
    # sub-leaf 'lvl'. This should limit looping to at most two iterations.
    # Though it might fail... Let's see if there are bug reports for weird
    # architectures.
    sl = lvl - one(UInt32)

    while (true)
        eax, ebx, ecx, edx = cpuid(leaf, sl)
        # still at a valid cache level ?
        eax & 0x1f == 0 && return 0
        # is this a data cache or shared cache level, eax[0:4]?
        # and is it the correct level, eax[5:7]?
        (eax & 0x01 == 0x01) && ((eax >> 5) & 0x07) == lvl &&
            return __datacachesize(eax, ebx, ecx)
        sl += one(UInt32)
    end
end

"""
    cacheinclusive()
    cacheinclusive(lvl::Integer)

Obtain information on the CPU's *data* cache inclusiveness. Returns `true`
for a cache that is inclusive of the lower cache levels, and `false` otherwise.

Determine the data cache size for each cache level as reported by the CPU
using a set of calls to the `cpuid` instruction.  Returns a tuple with the
tuple indices matching the cache levels.

If given an integer, then the data cache inclusiveness of the respective cache level
will be returned.  This is significantly faster than the tuple version above.
"""
function cacheinclusive end

function cacheinclusive()

    function cacheinclusive_level(leaf, sl::UInt32)
        eax, ebx, ecx, edx = cpuid(leaf, sl)
        # if eax is zero in the lowest 5 bits, we've reached the sentinel.
        eax & 0x1f == 0 && return ()
        # could do a sanity check: cache level reported in eax bits 5:7
        # if lowest bit on eax is zero, then its not a data cache
        eax & 0x01 == 0 && return cacheinclusive_level(leaf, sl + one(UInt32))
        # otherwise this should be a valid data or shared cache level
        ((edx & 0x02) != 0x00, cacheinclusive_level(leaf, sl + one(UInt32))...)
    end

    # TODO: This is awkwardly slow and requires some rework.
    #       Potential approach: Recurse to the last found cache level, there
    #       allocate a small array, then fill the array when leaving each
    #       recursion level.

    leaf = 0x0000_0004
    hasleaf(leaf) && return (cacheinclusive_level(leaf,zero(UInt32))...,)
    # no cache data available
    ()
end

cacheinclusive(lvl::Integer) = cacheinclusive(UInt32(lvl))

function cacheinclusive(lvl::UInt32) ::Int
    @_inline_meta

    lvl == 0 && return false

    leaf = 0x0000_0004
    hasleaf(leaf) || _throw_unsupported_leaf(leaf)

    # Loop over all subleaves of leaf 0x04 until the cache level bits say
    # we've reached to target level, then check whether it's a data cache or
    # shared cache.

    # Assuming cache enumeration is linear, it is sufficient to start at
    # sub-leaf 'lvl'. This should limit looping to at most two iterations.
    # Though it might fail... Let's see if there are bug reports for weird
    # architectures.
    sl = lvl - one(UInt32)

    while (true)
        eax, ebx, ecx, edx = cpuid(leaf, sl)
        # still at a valid cache level ?
        eax & 0x1f == 0 && return false
        # is this a data cache or shared cache level, eax[0:4]?
        # and is it the correct level, eax[5:7]?
        (eax & 0x01 == 0x01) && ((eax >> 5) & 0x07) == lvl &&
            return (edx & 0x02) != 0x00
        sl += one(UInt32)
    end
end

"""
    has_cpu_frequencies()

Determine whether the CPU provides clock frequency information.  If true, then
`cpu_base_frequency()`, `cpu_max_frequency()` and `cpu_bus_frequency()` should
be expected to return sensible information.
"""
function has_cpu_frequencies() ::Bool

    leaf = 0x0000_0016
    hasleaf(leaf) || return false

    # frequencies are provided if any of the bits in question are non-zero
    eax, ebx, ecx = cpuid(leaf)
    (eax & 0xffff) != zero(UInt32) ||
    (ebx & 0xffff) != zero(UInt32) ||
    (ecx & 0xffff) != zero(UInt32)
end

"""
    cpu_base_frequency()

Determine the CPU nominal base frequency in MHz as reported directly from the
CPU through a `cpuid` instruction call.  Returns zero if the CPU doesn't
provide base frequency information.

The actual cpu frequency might be lower due to throttling, or higher due to
frequency boosting (see `cpu_max_frequency`).
"""
function cpu_base_frequency() ::Int

    leaf = 0x0000_0016
    hasleaf(leaf) || return 0

    eax, ebx, ecx, edx = cpuid(leaf)
    eax & 0xffff
end

"""
    cpu_max_frequency()

Determine the maximum CPU frequency in MHz as reported directly from the CPU
through a `cpuid` instrauction call.  The maximum frequency typically refers
to the CPU's boost frequency.  Returns zero if the CPU doesn't provide maximum
frequency information.
"""
function cpu_max_frequency() ::Int

    leaf = 0x0000_0016
    hasleaf(leaf) || return 0

    eax, ebx, ecx, edx = cpuid(leaf)
    ebx & 0xffff
end

"""
    cpu_bus_frequency()

Determine the bus CPU frequency in MHz as reported directly from the CPU through
a `cpuid` instrauction call.  Returns zero if the CPU doesn't
provide bus frequency information.
"""
function cpu_bus_frequency() ::Int

    leaf = 0x0000_0016
    hasleaf(leaf) || return 0

    eax, ebx, ecx, edx = cpuid(leaf)
    ecx & 0xffff
end

"""
    perf_revision()

Determine the revision number of the performance monitoring unit.

This information is only available if `cpufeature(PDCM) == true`.
"""
function perf_revision() ::Int

    leaf = 0x0000_000a
    hasleaf(leaf) || return 0

    eax, ebx, ecx, edx = cpuid(leaf)
    eax & 0xff
end

"""
    perf_fix_counters()

Determine the number of fixed-function performance counters on the executing
machine.

This information is only available if `cpufeature(PDCM) == true`.
"""
function perf_fix_counters() ::Int

    leaf = 0x0000_000a
    hasleaf(leaf) || return 0

    eax, ebx, ecx, edx = cpuid(leaf)
    edx & 0x1f
end

"""
    perf_gen_counters()

Determine the number of general purpose counters performance counters on the
executing CPU.  Number of counters is given as per logical processor.

This information is only available if `cpufeature(PDCM) == true`.
"""
function perf_gen_counters() ::Int

    leaf = 0x0000_000a
    hasleaf(leaf) || return 0

    eax, ebx, ecx, edx = cpuid(leaf)
    (eax >> 8) & 0xff
end

"""
    perf_fix_bits()

Determine the number of bits fixed-function counters performance counters on
the executing CPU.

This information is only available if `cpufeature(PDCM) == true`.
"""
function perf_fix_bits() ::Int

    leaf = 0x0000_000a
    hasleaf(leaf) || return 0

    eax, ebx, ecx, edx = cpuid(leaf)
    (edx >> 5) & 0xff
end

"""
    perf_gen_bits()

Determine the number of bits general purpose counters performance counters on
the executing CPU.

This information is only available if `cpufeature(PDCM) == true`.
"""
function perf_gen_bits() ::Int

    leaf = 0x0000_000a
    hasleaf(leaf) || return 0

    eax, ebx, ecx, edx = cpuid(leaf)
    (eax >> 16) & 0xff
end

"""
    cpuinfo()

Generate a markdown table with the results of all of the CPU querying
functions provided by the module `CpuId`.  Intended to give a quick overview
for diagnostic purposes e.g. in log files.
"""
function cpuinfo()

    unsupported = "Not supported by CPU"
    cachesz = cachesize()
    modelfl = cpumodel()

    address = string(address_size(), " bits virtual, ", physical_address_size(), " bits physical")
    cache   = string("Level ", 1:length(cachesz), " : ", map(x->div(x,1024), cachesz), " kbytes")
    cachels = string(cachelinesize(), " byte cache line size")
    cores = string( CpuId.cpucores(), " physical cores, "
                  , CpuId.cputhreads(), " logical cores (on executing CPU)")
    frequencies = !has_cpu_frequencies() ? unsupported :
                        string(cpu_base_frequency(), " / ",
                               cpu_max_frequency(), " MHz (base/max), ",
                               cpu_bus_frequency(), " MHz bus")
    hyperthreading = (CpuId.cpucores() == CpuId.cputhreads() ?  "No " : "") * "Hyperthreading hardware capability detected"
    hypervisor = hypervised() ? "Yes, $(hvvendor())" : "No"
    model = string("Family: 0x",     string(modelfl[:Family],   base=16, pad=2),
                   ", Model: 0x",    string(modelfl[:Model],    base=16, pad=2),
                   ", Stepping: 0x", string(modelfl[:Stepping], base=16, pad=2),
                   ", Type: 0x",     string(modelfl[:CpuType],  base=16, pad=2))
    simd = string(simdbits(), " bit = ", simdbytes(), " byte max. SIMD vector size" )
    tsc = string("TSC is ", (cpufeature(TSC) ? "" : "not "), "accessible via `rdtsc`")
    tscinv = cpufeature(TSCINV) ? "TSC runs at constant rate (invariant from clock frequency)" :
                                   "TSC increased at every clock cycle (non-invariant TSC)"
    perfmon = !cpufeature(PDCM) ? "Performance Monitoring Counters (PMC) are not supported" :
                                  "Performance Monitoring Counters (PMC) revision $(perf_revision())"
    perfmon2 = !cpufeature(PDCM) ?  [] :
        [ ["", "Available hardware counters per logical core:"]
        , ["", "$(perf_fix_counters()) fixed-function counters of $(perf_fix_bits()) bit width"]
        , ["", "$(perf_gen_counters()) general-purpose counters of $(perf_gen_bits()) bit width"] ]
    ibs = !cpufeature(IBS) ? [] : [["", "CPU supports AMD's Instruction Based Sampling (IBS)"]]

    MarkdownString( MarkdownTable( [
        [ "Cpu Property",       "Value"              ],
        #----------------------------------------------
        [ "Brand",              strip(cpubrand())    ],
        [ "Vendor",             cpuvendor()          ],
        [ "Architecture",       cpuarchitecture()    ],
        [ "Model",              model                ],
        [ "Cores",              cores                ],
        [ "",                   hyperthreading       ],
        [ "Clock Frequencies",  frequencies          ],
        [ "Data Cache",         cache                ],
        [ "",                   cachels              ],
        [ "Address Size",       address              ],
        [ "SIMD",               simd                 ],
        [ "Time Stamp Counter", tsc                  ],
        [ "",                   tscinv               ],
        [ "Perf. Monitoring",   perfmon              ],
                                perfmon2...,
                                ibs...,
        [ "Hypervisor",         hypervisor           ],
       ], [:l, :l] ) )
end

end # module CpuId

#=--- end of file ---------------------------------------------------------=#


text:
function EquinoctialOrbit(arg0::Orbit)
    return EquinoctialOrbit((Orbit,), arg0)
end

function EquinoctialOrbit(arg0::PVCoordinates, arg1::Frame, arg2::AbsoluteDate, arg3::jdouble)
    return EquinoctialOrbit((PVCoordinates, Frame, AbsoluteDate, jdouble), arg0, arg1, arg2, arg3)
end

function EquinoctialOrbit(arg0::TimeStampedPVCoordinates, arg1::Frame, arg2::jdouble)
    return EquinoctialOrbit((TimeStampedPVCoordinates, Frame, jdouble), arg0, arg1, arg2)
end

function EquinoctialOrbit(arg0::jdouble, arg1::jdouble, arg10::jdouble, arg11::jdouble, arg12::PositionAngle, arg13::Frame, arg14::AbsoluteDate, arg15::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return EquinoctialOrbit((jdouble, jdouble, jdouble, jdouble, PositionAngle, Frame, AbsoluteDate, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg11, arg12, arg13, arg14, arg15, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function EquinoctialOrbit(arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::PositionAngle, arg7::Frame, arg8::AbsoluteDate, arg9::jdouble)
    return EquinoctialOrbit((jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, PositionAngle, Frame, AbsoluteDate, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function add_kepler_contribution(obj::EquinoctialOrbit, arg0::PositionAngle, arg1::jdouble, arg2::Vector{jdouble})
    return jcall(obj, "addKeplerContribution", void, (PositionAngle, jdouble, Vector{jdouble}), arg0, arg1, arg2)
end

function eccentric_to_mean(::Type{EquinoctialOrbit}, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(EquinoctialOrbit, "eccentricToMean", jdouble, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function eccentric_to_true(::Type{EquinoctialOrbit}, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(EquinoctialOrbit, "eccentricToTrue", jdouble, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function get_a(obj::EquinoctialOrbit)
    return jcall(obj, "getA", jdouble, ())
end

function get_a_dot(obj::EquinoctialOrbit)
    return jcall(obj, "getADot", jdouble, ())
end

function get_e(obj::EquinoctialOrbit)
    return jcall(obj, "getE", jdouble, ())
end

function get_e_dot(obj::EquinoctialOrbit)
    return jcall(obj, "getEDot", jdouble, ())
end

function get_equinoctial_ex(obj::EquinoctialOrbit)
    return jcall(obj, "getEquinoctialEx", jdouble, ())
end

function get_equinoctial_ex_dot(obj::EquinoctialOrbit)
    return jcall(obj, "getEquinoctialExDot", jdouble, ())
end

function get_equinoctial_ey(obj::EquinoctialOrbit)
    return jcall(obj, "getEquinoctialEy", jdouble, ())
end

function get_equinoctial_ey_dot(obj::EquinoctialOrbit)
    return jcall(obj, "getEquinoctialEyDot", jdouble, ())
end

function get_hx(obj::EquinoctialOrbit)
    return jcall(obj, "getHx", jdouble, ())
end

function get_hx_dot(obj::EquinoctialOrbit)
    return jcall(obj, "getHxDot", jdouble, ())
end

function get_hy(obj::EquinoctialOrbit)
    return jcall(obj, "getHy", jdouble, ())
end

function get_hy_dot(obj::EquinoctialOrbit)
    return jcall(obj, "getHyDot", jdouble, ())
end

function get_i(obj::EquinoctialOrbit)
    return jcall(obj, "getI", jdouble, ())
end

function get_i_dot(obj::EquinoctialOrbit)
    return jcall(obj, "getIDot", jdouble, ())
end

function get_l(obj::EquinoctialOrbit, arg0::PositionAngle)
    return jcall(obj, "getL", jdouble, (PositionAngle,), arg0)
end

function get_l_dot(obj::EquinoctialOrbit, arg0::PositionAngle)
    return jcall(obj, "getLDot", jdouble, (PositionAngle,), arg0)
end

function get_le(obj::EquinoctialOrbit)
    return jcall(obj, "getLE", jdouble, ())
end

function get_le_dot(obj::EquinoctialOrbit)
    return jcall(obj, "getLEDot", jdouble, ())
end

function get_lm(obj::EquinoctialOrbit)
    return jcall(obj, "getLM", jdouble, ())
end

function get_lm_dot(obj::EquinoctialOrbit)
    return jcall(obj, "getLMDot", jdouble, ())
end

function get_lv(obj::EquinoctialOrbit)
    return jcall(obj, "getLv", jdouble, ())
end

function get_lv_dot(obj::EquinoctialOrbit)
    return jcall(obj, "getLvDot", jdouble, ())
end

function get_type(obj::EquinoctialOrbit)
    return jcall(obj, "getType", OrbitType, ())
end

function interpolate(obj::EquinoctialOrbit, arg0::AbsoluteDate, arg1::Stream)
    return jcall(obj, "interpolate", EquinoctialOrbit, (AbsoluteDate, Stream), arg0, arg1)
end

function mean_to_eccentric(::Type{EquinoctialOrbit}, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(EquinoctialOrbit, "meanToEccentric", jdouble, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function shifted_by(obj::EquinoctialOrbit, arg0::jdouble)
    return jcall(obj, "shiftedBy", EquinoctialOrbit, (jdouble,), arg0)
end

function to_string(obj::EquinoctialOrbit)
    return jcall(obj, "toString", JString, ())
end

function true_to_eccentric(::Type{EquinoctialOrbit}, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(EquinoctialOrbit, "trueToEccentric", jdouble, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

ext: jl


text:
module TestExamples
using PatternDispatch

@pattern f(x) =  x
@pattern f(2) = 42

println(Any[f(x) for x=1:4])

println()
show_dispatch(f)
show_dispatch(f, (Int,))

@pattern f2((x,y::Int)) = x*y
@pattern f2([x,y::Int]) = x/y
@pattern f2(x)          = nothing

println()
@show f2((2,5))
@show f2((4,3))
@show f2([4,3])
@show f2((4,'a'))
@show f2(Any[4,'a'])
@show f2(1)
@show f2("foo")
@show f2((1,))
@show f2((1,2,3))

@pattern f3(v~[x::Int, y::Int]) = Any[v,x*y]
@pattern f3(v) = nothing

println()
@show f3([3, 2])
@show f3(Any[3, 2])
@show f3([3, 2.0])

@pattern f4($nothing) = 1
@pattern f4(x)        = 2

println()
@show f4(nothing)
@show f4(1)
@show f4(:x)
@show f4("foo")

println()
@pattern ambiguous((x,y),z) = 2
@pattern ambiguous(x,(1,z)) = 3

end # module


text:
module EverestCustomBirdTutorialTrigger

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Custom Bird Tutorial (Everest)" => Ahorn.EntityPlacement(
        Maple.EverestCustomBirdTutorialTrigger,
        "rectangle"
    )
)

end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
import MathProgBase
import MathProgBase.SolverInterface

export BaronSolver
struct BaronSolver <: MathProgBase.AbstractMathProgSolver
    options
end

BaronSolver(;kwargs...) = BaronSolver(kwargs)

mutable struct BaronMathProgModel <: MathProgBase.AbstractNonlinearModel
    inner::BaronModel
end

MathProgBase.NonlinearModel(s::BaronSolver) = BaronMathProgModel(;s.options...)
MathProgBase.LinearQuadraticModel(s::BaronSolver) = MathProgBase.NonlinearToLPQPBridge(MathProgBase.NonlinearModel(s))

function MathProgBase.loadproblem!(mpb_model::BaronMathProgModel,
    nvar, ncon, xˡ, xᵘ, gˡ, gᵘ, sense,
    d::MathProgBase.AbstractNLPEvaluator)

    model = mpb_model.inner
    MathProgBase.initialize(d, [:ExprGraph])

    @assert nvar == length(xˡ) == length(xᵘ)
    for i in 1:nvar
        push!(model.variable_info, VariableInfo(xˡ[i], xᵘ[i], :Cont, "x$i"))
    end

    @assert ncon == length(gˡ) == length(gᵘ)
    for i in 1:ncon
    	constr_expr = verify_support(MathProgBase.constr_expr(d, c))
        push!(model.constraint_info, ConstraintInfo(constr_expr, gˡ[i], gᵘ[i], "e$i"))
    end

    m.objective_info = ObjectiveInfo(verify_support(MathProgBase.obj_expr(d)), sense)

    m.probfile = joinpath(m.temp_dir, "baron_problem.bar")
    m.sumfile  = joinpath(m.temp_dir, "sum.lst")
    m.resfile  = joinpath(m.temp_dir, "res.lst")
    return m
end

function MathProgBase.setvartype!(m::BaronMathProgModel, cat::Vector{Symbol})
    @assert length(cat) == length(m.inner.variable_info)
    for i in 1:length(cat)
        var_cat = cat[i]
        @assert var_cat in (:Cont, :Bin, :Int)
        m.inner.variable_info[i].category = var_cat
    end
    return
end

function MathProgBase.setwarmstart!(m::BaronMathProgModel, v::Vector{Float64})
    @assert length(v) == length(m.inner.variable_info)
    for i in 1:length(v)
        m.inner.variable_info[i].start = v[i]
    end
    return
end

function MathProgBase.optimize!(mpb_model::BaronMathProgModel)
    write_bar_file(m.inner)
    run(`$baron_exec $(m.inner.prob_file)`)

    read_results(m.inner)
end

function MathProgBase.status(m::BaronMathProgModel)
    status = m.solution_info.status
    if status == NORMAL_COMPLETION
    	return :Optimal
    elseif status == INFEASIBLE
        return :Infeasible
    elseif status == UNBOUNDED
        return :Unbounded
    elseif status in (CPU_TIME_LIMIT, TIME_LIMIT, USER_INTERRUPTION)
        return :UserLimit
    else
        return :Error
    end
end

MathProgBase.numvar(m::BaronMathProgModel) = length(m.inner.variable_info)
MathProgBase.getsolution(m::BaronMathProgModel) = m.inner.solution_info.feasible_point
MathProgBase.getobjval(m::BaronMathProgModel) = m.inner.solution_info.objective_value
MathProgBase.getobjbound(m::BaronMathProgModel) = m.inner.solution_info.dual_bound
MathProgBase.getsolvetime(m::BaronMathProgModel) = m.inner.solution_info.wall_time


text:
import Base: tuple_type_cons, tuple_type_head, tuple_type_tail, tail

eltypes(::Type{T}) where {T} = map_types(eltype, T)

map_types(f, ::Type{Tuple{}}) = Tuple{}
function map_types(f, ::Type{T}) where {T<:Tuple}
    tuple_type_cons(f(tuple_type_head(T)), map_types(f, tuple_type_tail(T)))
end
map_types(f, ::Type{NamedTuple{names, types}}) where {names, types} =
    NamedTuple{names, map_types(f, types)}

all_types(f, ::Type{Tuple{}}, ::Type{T}) where {T<:Tuple} = true

function all_types(f, ::Type{S}, ::Type{T}) where {S<:Tuple, T<:Tuple}
    f(tuple_type_head(S), tuple_type_head(T)) && all_types(f, tuple_type_tail(S), tuple_type_tail(T))
end

all_types(f, ::Type{NamedTuple{n1, t1}}, ::Type{NamedTuple{n2, t2}}) where {n1, t1, n2, t2} =
    all_types(f, t1, t2)

map_params(f, ::Type{Tuple{}}) = ()
function map_params(f, ::Type{T}) where {T<:Tuple}
    (f(tuple_type_head(T)), map_params(f, tuple_type_tail(T))...)
end
map_params(f, ::Type{NamedTuple{names, types}}) where {names, types} =
    NamedTuple{names}(map_params(f, types))

buildfromschema(initializer, ::Type{T}) where {T} = buildfromschema(initializer, T, staticschema(T))

function buildfromschema(initializer, ::Type{T}, ::Type{NT}) where {T, NT<:NamedTuple}
    nt = map_params(initializer, NT)
    StructArray{T}(nt)
end

Base.@pure SkipConstructor(::Type) = false

@generated function foreachcolumn(f, x::StructArray{T, N, NamedTuple{names, types}}, xs...) where {T, N, names, types}
    exprs = Expr[]
    for (i, field) in enumerate(names)
        sym = QuoteNode(field)
        args = [Expr(:call, :getfieldindex, :(getfield(xs, $j)), sym, i) for j in 1:length(xs)]
        push!(exprs, Expr(:call, :f, Expr(:., :x, sym), args...))
    end
    push!(exprs, :(return nothing))
    Expr(:block, exprs...)
end

function createinstance(::Type{T}, args...) where {T}
    SkipConstructor(T) ? unsafe_createinstance(T, args...) : T(args...)
end

createinstance(::Type{T}, args...) where {T<:Union{Tuple, NamedTuple}} = T(args)

@generated function unsafe_createinstance(::Type{T}, args...) where {T}
    v = fieldnames(T)
    new_tup = Expr(:(=), Expr(:tuple, v...), :args)
    construct = Expr(:new, :T, (:(convert(fieldtype(T, $(Expr(:quote, sym))), $sym)) for sym in v)...)
    Expr(:block, new_tup, construct)
end


text:
module OnGridTests

using Interpolations, Base.Test

nx, ny, nz = 10, 8, 9
xg, yg, zg = 1:nx, 1:ny, 1:nz

f1(x) = sin(2pi/nx * (x-1))
f2(x,y) = f1(x) * cos(4pi/ny * (y-1))
f3(x,y,z) = f2(x,y) * log(z)

simpledegs = (Linear,Constant)
bcdegs = (Quadratic,)
gridbehvs = (OnCell,OnGrid)
bcs = (Flat,Line,Free,Periodic,Reflect)

function OneD()
    println("Testing evaluation on grid and boundary in 1D...")
    for deg in simpledegs, gb in gridbehvs
        fg = Float64[f1(x) for x in xg]

        itp = Interpolation(fg, deg(gb()), ExtrapError())

        for i in xg
            @test_approx_eq_eps itp[i] fg[i] sqrt(eps())
        end

        if gb == OnCell
            for xb in [.5, nx+.5]
                itp[xb] # if we're out of bounds, this throws
            end
        end
    end

    for deg in bcdegs, gb in gridbehvs, bc in bcs
        fg = Float64[f1(x) for x in xg]

        itp = Interpolation(fg, deg(bc(),gb()), ExtrapError())

        for i in xg
            @test_approx_eq_eps itp[i] fg[i] sqrt(eps())
        end

        if gb == OnCell
            for xb in [.5, nx+.5]
                itp[xb]
            end
        end
    end
end

function TwoD()
    println("Testing evaluation on grid and boundary in 2D...")
    for deg in simpledegs, gb in gridbehvs
        fg = Float64[f2(x,y) for x in xg, y in yg]
        itp = Interpolation(fg, deg(gb()), ExtrapError())

        for i in xg, j in yg
            @test_approx_eq_eps itp[i,j] fg[i,j] sqrt(eps())
        end

        if gb == OnCell
            for xb in [.5, nx+.5], yb in [.5,ny+.5]
                itp[xb,yb]
            end
        end
    end

    for deg in bcdegs, gb in gridbehvs, bc in bcs
        fg = Float64[f2(x,y) for x in xg, y in yg]
        itp = Interpolation(fg, deg(bc(),gb()), ExtrapError())

        for i in xg, j in yg
            @test_approx_eq_eps itp[i,j] fg[i,j] sqrt(eps())
        end

        if gb == OnCell
            for xb in [.5, nx+.5], yb in [.5,ny+.5]
                itp[xb,yb]
            end
        end
    end
end

function ThreeD()
    println("Testing evaluation on grid and boundary in 3D...")
    for deg in simpledegs, gb in gridbehvs
        fg = Float64[f3(x,y,z) for x in xg, y in yg, z in zg]
        itp = Interpolation(fg, deg(gb()), ExtrapError())

        for i in xg, j in yg, k in zg
            @test_approx_eq_eps itp[i,j,k] fg[i,j,k] sqrt(eps())
        end

        if gb == OnCell
            for xb in [.5, nx+.5], yb in [.5,ny+.5], zb in [.5, nz+.5]
                itp[xb,yb,zb]
            end
        end
    end

    for deg in bcdegs, gb in gridbehvs, bc in bcs
        fg = Float64[f3(x,y,z) for x in xg, y in yg, z in zg]
        itp = Interpolation(fg, deg(bc(),gb()), ExtrapError())

        for i in xg, j in yg, k in zg
            @test_approx_eq_eps itp[i,j,k] fg[i,j,k] sqrt(eps())
        end

        if gb == OnCell
            for xb in [.5, nx+.5], yb in [.5,ny+.5], zb in [.5, nz+.5]
                itp[xb,yb,zb]
            end
        end
    end
end

OneD()
TwoD()
ThreeD()

end


text:
# No primatives
function cgal_full(
        C_input,
        b_input,
        A!,
        A_adj!;
        n,
        d,
        scale_X = 1,
        scale_C = 1,
        max_iters = 1_000,
        tol = 1e-10,
        print_iter = 10,
)

    # Copy variables
    C = copy(C_input) .* scale_C
    b = copy(b_input) .* scale_X
    norm_b = norm(b_input)

    α = 1
    t = 1
    β0 = 1
    K = Inf
    dual_gap = Inf

    Xt = zeros(n, n)
    yt = zeros(d)

    # --- cache & other memory allocations ---
    qmax = min(n-1, Int(ceil(max_iters^0.25 * log(n))))
    cache = (
        A_X = zeros(d),
        dual_update = zeros(d),
        evec = init_cache_evec(n, qmax),
    )
    Dt = copy(C)
    v = zeros(n)

    # Keep track of things (credit: https://github.com/ZIB-IOL/FrankWolfe.jl)
    headers = ["Iteration", "Primal", "Dual Gap", "Infeas", "Time"]
    print_header(headers)
    time_start = time_ns()
    while t <= max_iters #&& dual_gap >= max(tol, eps())
        # --- Parameter updates (from Yurtsever et al Alg 6.1) ---
        βt = β0 * sqrt(t + 1)
        η = 2/(t + 1)

        # --- Gradient computation ---
        # * Non-allocating *
        A!(cache.A_X, Xt)
        cache.A_X .-= b
        cache.A_X .*= βt
        cache.A_X .+= yt
        for col = 1:size(Dt, 2)
            for k = SparseArrays.getcolptr(Dt)[col]:(SparseArrays.getcolptr(Dt)[col+1]-1)
                @views Dt[rowvals(Dt)[k], col] = C[rowvals(Dt)[k], col]
            end
        end
        A_adj!(Dt, cache.A_X)
        # Dt = C + A*(yt + βt*(diag(Xt) - b))

        # --- Eigenvector computation ---
        # * Custom Method *
        q = Int(ceil(t^0.25 * log(n)))
        ξ = approx_min_evec!(v, Dt, n=n, q=q, cache=cache.evec)
        # * ArnoldiMethod.jl (more stable than Lanzcos) *
        # F = partialschur(Dt, nev=1, which=SR(), tol=sqrt(n)*eps())
        # !F[2].converged && error("Eigvec computation did not converge")
        # ξ, vv = partialeigen(F[1])
        # v .= vv[:,1]
        # * Standard Julia function (dense matrix) *
        # d, V = eigen(Dt)
        # ξ = d[1]
        # v = V[:,1]

        # --- Primal update ---
        Xt .= Xt .- η.*Xt
        # TODO: Call directly to avoid allocation
        BLAS.ger!(η, v, v, Xt)

        # --- Dual update ---
        A!(cache.A_X, Xt)
        cache.dual_update .= cache.A_X .- b
        primal_infeas = norm(cache.dual_update) * 1/scale_X / (1 + norm_b)
        γ = min(β0, 4*α^2*βt*η^2 / primal_infeas^2)
        @. yt += γ*cache.dual_update
        obj_val = sum(C.*Xt)

        @. cache.dual_update = yt + βt*cache.dual_update
        dual_gap = obj_val + dot(cache.dual_update, cache.A_X) - ξ[1]

        if t == 1 || t % print_iter == 0
            print_iter_func((
                string(t),
                obj_val * 1/scale_C * 1/scale_X,
                dual_gap,
                primal_infeas,
                (time_ns() - time_start) / 1e9
            ))
        end
        t += 1
    end
    print_footer()

    return Xt, yt

end


text:
module GridTests

using Interpolations, Base.Test

# On-grid values
A = randn(4,10)
const EPS = sqrt(eps())
for it in (Constant(OnCell()), Linear(OnGrid()), Quadratic(Free(),OnCell()))
    for eb in (ExtrapNaN(), ExtrapError())
        itp = Interpolation(A, it, eb)
        for i = 1:size(A,1)
            for j = 1:size(A,2)
                @test_approx_eq_eps itp[i,j] A[i,j] EPS
            end
        end
        # v = itp[1:size(A,1), 1:size(A,2)]
        # @assert all(abs(v - A) .< EPS)
    end
end

A = randn(4,5,4)
for it in (Constant(OnCell()), Linear(OnGrid()), Quadratic(Free(),OnCell()))
    for eb in (ExtrapNaN(), ExtrapError())
        itp = Interpolation(A, it, eb)
        for k = 1:size(A,3), j = 1:size(A,2), i = 1:size(A,1)
            @test_approx_eq_eps itp[i,j,k] A[i,j,k] EPS
        end
        # v = itp[1:size(A,1), 1:size(A,2), 1:size(A,3)]
        # @assert all(abs(v - A) .< EPS)
    end
end
A = randn(4,5,4,3)
for it in (Constant(OnCell()), Linear(OnGrid()), Quadratic(Free(),OnCell()))
    for eb in (ExtrapNaN(), ExtrapError())
        itp = Interpolation(A, it, eb)
        for i = 1:size(A,1), j = 1:size(A,2), k = 1:size(A,3), l = 1:size(A,4)
            @test_approx_eq_eps itp[i,j,k,l] A[i,j,k,l] EPS
        end
        # v = itp[1:size(A,1), 1:size(A,2), 1:size(A,3), 1:size(A,4)]
        # @assert all(abs(v - A) .< EPS)
    end
end

A = float([1:4])
for it in (Constant(OnCell()), Linear(OnGrid()), Quadratic(Free(),OnCell()))
    itp = Interpolation(A, it, ExtrapError())
    @test_throws BoundsError itp[-0.8]
end
for it in (Constant(OnCell()), Linear(OnGrid()), Quadratic(Free(),OnCell()))
    itp = Interpolation(A, it, ExtrapNaN())
    @test isnan(itp[-0.8])
end

end


text:
module FooBroken

"""
This docstring contains a broken doctest:

```jldoctest
julia> 2 + 2
42
```
"""
function foo end

end


text:
using SimpleSMC
import SequentialMonteCarlo.SMCModel
using RNGPool
using SMCExamples.FiniteFeynmanKac
using Random

include("test.jl")

setRNGs(0)

d = 3
ffk = FiniteFeynmanKac.randomFiniteFK(d, 10)
model = FiniteFeynmanKac.makeSMCModel(ffk)

numParticles = 1024
numTrials = 2

## just run the algorithm a few times
testSMC(model, numParticles, numTrials)

testCSMC(model, numParticles, numTrials)

lM = FiniteFeynmanKac.makelM(ffk)
testCSMCBS(model, lM, numParticles, numTrials)


text:
export Simulator
export EndStateSimulator, RewardAccumulator, StateActionRecorder

"Abstract type for simulators."
abstract type Simulator end

(sim::Simulator)(sol::Solution, domain::Domain, state::State, spec::Specification) =
    run_simulator(sim, sol, domain, state, spec)
(sim::Simulator)(sol::Solution, domain::Domain, state::State, spec) =
	sim(sol, domain, state, Specification(spec))
(sim::Simulator)(sol::Solution, domain::Domain, state::State) =
    sim(sol, domain, state, NullSpecification())

"Simulator that returns end state of simulation."
@kwdef struct EndStateSimulator <: Simulator
    max_steps::Int = typemax(Int)
end

function run_simulator(sim::EndStateSimulator, sol::Solution,
					   domain::Domain, state::State, spec::Specification)
    for t in 1:sim.max_steps
		if is_goal(spec, domain, state) break end
        act = get_action(sol, t, state)
        state = transition(domain, state, act)
    end
    return state
end

"Simulator that returns accumulated reward."
@kwdef struct RewardAccumulator <: Simulator
    max_steps::Int = typemax(Int)
end

function run_simulator(sim::RewardAccumulator, sol::Solution,
					   domain::Domain, state::State, spec::Specification)
	reward = 0.0
	discount = get_discount(spec)
    for t in 1:sim.max_steps
		if is_goal(spec, domain, state) break end
        act = get_action(sol, t, state)
        next_state = transition(domain, state, act)
		reward += discount * get_reward(spec, domain, state, act, next_state)
		discount *= get_discount(spec)
		state = next_state
    end
    return reward
end

"Simulator that records the state-action trajectory, including the start state."
@kwdef struct StateActionRecorder <: Simulator
    max_steps::Int = typemax(Int)
end

function run_simulator(sim::StateActionRecorder, sol::Solution,
 					   domain::Domain, state::State, spec::Specification)
    actions = Term[]
    trajectory = State[state]
    for t in 1:sim.max_steps
		if is_goal(spec, domain, state) break end
        act = get_action(sol, t, state)
        state = transition(domain, state, act)
		push!(actions, act)
        push!(trajectory, state)
    end
    return (actions, trajectory)
end


text:
using SamplingHS3

## Create gaussian mixture
gm2_mix = []
for zi=1:10;
    found_es = false
    vd_pre   = NaN * ones(2,2)
    while(!found_es)
        vd = qr(rand(2,2))[1]
        vd_pre = vd*diagm(rand(2))*vd'
        found_es = issymmetric(vd_pre)
    end
    gpi = SamplingHS3.GaussParam( 8*rand(2) , vd_pre )
    gm2_mix = [ gm2_mix ;  gpi ]
end
gm2 = SamplingHS3.GaussMixture(gm2_mix,0.5+rand(10))
# bounding box
bb_G = [eye(2);-eye(2)]
bb_h = [10.;10;10;10]

# quick plot as sanity check
using Plots
xx  = SamplingHS3.draw_gaussmixture( gm2 , 100000 )
SamplingHS3.gr_hist_2d( xx , nbins=40 )

## Use SamplingHS3

# Input:
# Density function: (x) -> SamplingHS3.pdf_gaussmixture_nonnormalized(gm2,x)
# Bounding box, given by bb_G,bb_h

# Sampling Steps:
# 1. generate uniform samples over polytope
X_init = SamplingHS3.sample_unif_hr( bb_G , bb_h , zeros(2) , 1000 , 2000 )

# 2. create "cost" function, i.e. for a pdf given
#    by f(x), we use f_cost(x) = -log(f(x))
f_cost = (x) -> -log(SamplingHS3.pdf_gaussmixture_nonnormalized(gm2,x))

# Optional: evaluate cost of initial samples (can be supplied to the ns function)
# x_init_cost = [  f_cost(X_init[:,zi]) for zi in (1:2000) ]

# 3. Run NS:
# Required parameters are:
# (obviously): initial samples, cost function, polytope description
# (and):       the NSConfig struct
ns_result = SamplingHS3.nested_sampling( X_init[:,:] , f_cost , bb_G , bb_h , SamplingHS3.NSConfig(8,0.95,SamplingHS3.StopCriterion(max_time=10.0),SamplingHS3.MCMCConfig()) )

# 4. Plot evolution of live samples during ns run (Just because it looks cool..)
x_live_samples_1 = SamplingHS3.reconstruct_live_samples(ns_result.ns_iter_data)
x_live  = x_live_samples_1[1]
xf_live = x_live_samples_1[2]
for zi in 1:length(x_live)
    #zcolor_i = 1:size(x_live[zi],2)
    zcolor_i = xf_live[zi]
    display( scatter(x_live[zi][1,:] , x_live[zi][2,:]  , ma = 0.75 , zcolor=zcolor_i , m=2, msw=0 , c=:lightrainbow) )
end

# 5. Generate an actual sample of the probability distribution from the ns run:


text:
# api.jl

include("./api.jl")

using SQLite
using DataFrames
using Dates

# Connect
connect(path::String)::SQLite.DB = SQLite.DB(path)

connect()::SQLite.DB = SQLite.DB()

archive(db, table::String, domainItems::Array{T, 1} where {T <: Any}) = begin
   dbi = [ createDatabaseItem( i ; action="CREATE" ) for i in domainItems]
   DataFrame( dbi ) |> SQLite.load!(db, table)
end # create

store(db, table::String, items::Array{T, 1} where {T <: Any}) = begin
   DataFrame( items ) |> SQLite.load!(db, table)
end # tore

# retrieve all item from a table
retrieve(db, table::String)::DataFrame = begin
   r = SQLite.DBInterface.execute(db, "select * from $table") |> DataFrame
   return r
end # retrieve

# retrieve item form a table based on a sql condition
retrieve(db, table::String, condition::String )::DataFrame = begin
   r = SQLite.DBInterface.execute( db, "select * from $table where $condition")  |> DataFrame
   return r
end # retrieve

# run custom function
runfunct(funct, x, y, z) = funct(x, y, z)


text:
using Documenter, YaoLang

makedocs(;
    modules = [YaoLang],
    format = Documenter.HTML(prettyurls = !("local" in ARGS)),
    pages = ["Home" => "index.md", "Compilation" => "compilation.md"],
    repo = "https://github.com/QuantumBFS/YaoLang.jl",
    sitename = "YaoLang.jl",
)

deploydocs(; repo = "github.com/QuantumBFS/YaoLang.jl")


text:
export
ruleVBWishartOut,
ruleSPWishartOutNPP

ruleSPWishartOutNPP(msg_out::Nothing,
                    msg_v::Message{PointMass, MatrixVariate},
                    msg_nu::Message{PointMass, Univariate}) =
    Message(MatrixVariate, Wishart, v=deepcopy(msg_v.dist.params[:m]), nu=deepcopy(msg_nu.dist.params[:m]))

ruleVBWishartOut(   dist_out::Any,
                    dist_v::ProbabilityDistribution{MatrixVariate},
                    dist_nu::ProbabilityDistribution{Univariate}) =
    Message(MatrixVariate, Wishart, v=unsafeMean(dist_v), nu=unsafeMean(dist_nu))


text:
module JackPolynomials

import DynamicPolynomials
export JackPolynomial
export ZonalPolynomial
export ZonalQPolynomial
export SchurPolynomial
export Zonal
export ZonalQ
export Schur
export Jack

function isPartition(lambda::Vector{<:Integer})
  return all(diff(lambda) .<= 0) && all(lambda .>= 0)
end

function dualPartition(lambda::Vector{T}) where {T<:Integer}
  out = T[]
  if !isempty(lambda)
    for i = 1:lambda[1]
      push!(out, sum(lambda .>= i))
    end
  end
  return out
end

function betaratio(
  kappa::Vector{I},
  mu::Vector{I},
  k::I,
  alpha::T,
) where {T<:Real,I<:Integer}
  muk = mu[k]
  t = k - alpha * muk
  u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
  v = map(i -> t - i + alpha * mu[i], 1:(k-1))
  muPrime = dualPartition(mu)
  w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
  prod1 = prod(u ./ (u .+ alpha .- 1))
  prod2 = prod((v .+ alpha) ./ v)
  prod3 = prod((w .+ alpha) ./ w)
  return alpha * prod1 * prod2 * prod3
end

function _N(lambda::Vector{I}, mu::Vector{I}) where {I<:Integer}
  prods = map(i -> prod(Iterators.drop(lambda .+ 1, i)), 1:length(lambda))
  return sum(mu .* prods)
end

# ------------------------------------------------------------------------------
#~~ Jack polynomial ~~~~##
# ------------------------------------------------------------------------------
"""
    Jack(x, lambda, alpha)

Evaluates a Jack polynomial.

# Arguments
- `x`: vector of real or complex numbers
- `lambda`: partition of an integer
- `alpha`: alpha parameter
"""
function Jack(
  x::Vector{C},
  lambda::Vector{I},
  alpha::T,
) where {T<:Real,I<:Integer,C<:Number}
  if !isPartition(lambda)
    error("`lambda` must be a partition of an integer")
  end
  if alpha <= 0
    error("`alpha` must be positive")
  end
  function jac(m::I, k::I, mu::Vector{I}, nu::Vector{I}, beta::Real)
    if isempty(nu) || nu[1] == 0 || m == 0
      return C(1)
    end
    if length(nu) > m && nu[m+1] > 0
      return C(0)
    end
    if m == 1
      return x[1]^(nu[1]) * prod(1 .+ alpha .* collect(1:(nu[1]-1)))
    end
    v = S[_N(lambda, nu), m]
    if k == 0 && !ismissing(v)
      return v
    end
    i = max(1, k)
    s = jac(m - 1, 0, nu, nu, T(1)) * beta * x[m]^(sum(mu) - sum(nu))
    while length(nu) >= i && nu[i] > 0
      if length(nu) == i || nu[i] > nu[i+1]
        nuPrime = copy(nu)
        nuPrime[i, 1] -= 1
        gamma = beta * betaratio(mu, nu, i, alpha)
        if nu[i] > 1
          s = s + jac(m, i, mu, nuPrime, gamma)
        else
          s =
            s +
            jac(m - 1, 0, nuPrime, nuPrime, T(1)) * gamma * x[m]^(sum(mu) - sum(nuPrime))
        end
      end
      i += 1
    end
    if k == 0
      S[_N(lambda, nu), m] = s
    end
    return s
  end # end jac --------------------------------------------------------------
  lx = length(x)
  S = Array{Union{Missing,C}}(missing, _N(lambda, lambda), lx)
  jac(lx, 0, lambda, lambda, T(1))
end

# ------------------------------------------------------------------------------
#~~ Symbolic Jack polynomial ~~~~##
# ------------------------------------------------------------------------------
function JackPolynomial0(
  m::I, lambda::Vector{I}, alpha::T
) where {T<:Real,I<:Integer}
  function jac(m::I, k::I, mu::Vector{I}, nu::Vector{I}, beta::T)
    if isempty(nu) || nu[1] == 0 || m == 0
      return T(1)
    end
    lnu = length(nu)
    if lnu > m && nu[m+1] > 0
      return T(0)
    end
    if m == 1
      return x[1]^(nu[1]) * prod(1 .+ alpha .* collect(1:(nu[1]-1)))
    end
    v = S[_N(lambda, nu), m]
    if k == 0 && !ismissing(v)
      return v
    end
    i = max(1, k)
    s = jac(m - 1, 0, nu, nu, T(1)) * beta * x[m]^(sum(mu) - sum(nu))
    while lnu >= i && nu[i] > 0
      if lnu == i || nu[i] > nu[i+1]
        nuPrime = copy(nu)
        nuPrime[i, 1] -= 1
        gamma = beta * betaratio(mu, nu, i, alpha)
        if nu[i] > 1
          s += jac(m, i, mu, nuPrime, gamma)
        else
          s += jac(m - 1, 0, nuPrime, nuPrime, T(1)) * gamma * x[m]^(sum(mu) - sum(nuPrime))
        end
      end
      i += 1
    end
    if k == 0
      S[_N(lambda, nu), m] = s
    end
    return s
  end
  DynamicPolynomials.@polyvar x[1:m]
  S = Array{Union{Missing,DynamicPolynomials.Polynomial{true,T}}}(
    missing,
    _N(lambda, lambda),
    m,
  )
  jac(m, 0, lambda, lambda, T(1))
end

"""
    JackPolynomial(m, lambda, alpha)

Symbolic Jack polynomial. The coefficients of the polynomial will have the
same type as `alpha`.

# Arguments
- `m`: integer, the number of variables
- `lambda`: partition of an integer
- `alpha`: alpha parameter
"""
function JackPolynomial(
  m::I, lambda::Vector{I}, alpha::T
) where {T<:Real,I<:Integer}
  if !isPartition(lambda)
    error("`lambda` must be a partition of an integer")
  end
  if alpha <= 0
    error("`alpha` must be positive")
  end
  jack = JackPolynomial0(m, lambda, alpha)
  if(typeof(jack) == T)
    DynamicPolynomials.@polyvar x[1:m]
    jack = sum(T(0) * x) + jack
  end
  return jack
end

# ------------------------------------------------------------------------------
#~~ Zonal polynomial ~~~~##
# ------------------------------------------------------------------------------
function _i(lambda::Vector{T}) where {T<:Integer}
  out = T[]
  for i in 1:length(lambda)
    out = vcat(out, repeat([i], lambda[i]))
  end
  return out
end

function _j(lambda::Vector{T}) where {T<:Integer}
  out = T[]
  mu = filter(x -> x > 0, lambda)
  for m in mu
    out = vcat(out, 1:m)
  end
  return out
end

function hookLengths(lambda::Vector{<:Integer}, alpha::Real)
  i = _i(lambda)
  j = _j(lambda)
  lambdap = dualPartition(lambda)
  upper = lambdap[j] .- i .+ alpha .* (lambda[i] .- j .+ 1)
  lower = lambdap[j] .- i .+ 1 .+ alpha .* (lambda[i] .- j)
  return vcat(upper, lower)
end

"""
    Zonal(x, lambda)

Evaluates a zonal polynomial.

# Arguments
- `x`: vector of real or complex numbers
- `lambda`: partition of an integer
"""
function Zonal(
  x::Vector{<:Union{R,Complex{R}}},
  lambda::Vector{<:Integer},
) where {R<:Real}
  jack = Jack(x, lambda, R(2))
  jlambda = prod(hookLengths(lambda, R(2)))
  n = sum(lambda)
  return jack * 2^n * factorial(n) / jlambda
end

# ------------------------------------------------------------------------------
#~~ Symbolic zonal polynomial ~~~~##
# ------------------------------------------------------------------------------
"""
    ZonalPolynomial(m, lambda[, type])

Symbolic zonal polynomial.

# Arguments
- `m`: integer, the number of variables
- `lambda`: partition of an integer
- `type`: the type of the coefficients of the polynomial; default `Rational`
"""
function ZonalPolynomial(
  m::I,
  lambda::Vector{I},
  type::Type = Rational
) where {I<:Integer}
  jack = JackPolynomial(m, lambda, type(2))
  jlambda = prod(hookLengths(lambda, type(2)))
  n = sum(lambda)
  return jack * 2^n * factorial(n) / jlambda
end

# ------------------------------------------------------------------------------
#~~ Quaternionic zonal polynomial ~~~~##
# ------------------------------------------------------------------------------
"""
    ZonalQ(x, lambda)

Evaluates a quaternionic zonal polynomial.

# Arguments
- `x`: vector of real or complex numbers
- `lambda`: partition of an integer
"""
function ZonalQ(
  x::Vector{<:Union{R,Complex{R}}},
  lambda::Vector{<:Integer},
) where {R<:Real}
  jack = Jack(x, lambda, R(0.5))
  jlambda = prod(hookLengths(lambda, R(0.5)))
  n = sum(lambda)
  return jack * factorial(n) / 2^n / jlambda
end

# ------------------------------------------------------------------------------
#~~ Symbolic quaternionic zonal polynomial ~~~~##
# ------------------------------------------------------------------------------
"""
    ZonalQPolynomial(m, lambda[, type])

Symbolic quaternionic zonal polynomial.

# Arguments
- `m`: integer, the number of variables
- `lambda`: partition of an integer
- `type`: the type of the coefficients of the polynomial; default `Rational`
"""
function ZonalQPolynomial(
  m::I,
  lambda::Vector{I},
  type::Type = Rational
) where {I<:Integer}
  jack = JackPolynomial(m, lambda, type(0.5))
  jlambda = prod(hookLengths(lambda, type(0.5)))
  n = sum(lambda)
  return jack * factorial(n) / 2^n / jlambda
end

# ------------------------------------------------------------------------------
#~~ Schur polynomial ~~~~##
# ------------------------------------------------------------------------------
"""
    Schur(x, lambda)

Evaluates a Schur polynomial.

# Arguments
- `x`: vector of real or complex numbers
- `lambda`: partition of an integer
"""
function Schur(x::Vector{T}, lambda::Vector{I}) where {T<:Number,I<:Integer}
  if !isPartition(lambda)
    error("`lambda` must be a partition of an integer")
  end
  function sch(m::Int64, k::Int64, nu::Vector{I})
    if isempty(nu) || nu[1] == 0 || m == 0
      return T(1)
    end
    if length(nu) > m && nu[m+1] > 0
      return T(0)
    end
    if m == 1
      return x[1]^nu[1]
    end
    v = S[_N(lambda, nu), m]
    if !ismissing(v)
      return v
    end
    s = sch(m - 1, 1, nu)
    lnu = length(nu)
    i = k
    while lnu >= i && nu[i] > 0
      if lnu == i || nu[i] > nu[i+1]
        nup = copy(nu)
        nup[i] -= 1
        if nu[i] > 1
          s += x[m] * sch(m, i, nup)
        else
          s += x[m] * sch(m - 1, 1, nup)
        end
      end
      i += 1
    end
    if k == 1
      S[_N(lambda, lambda), m] = s
    end
    return s
  end # end sch --------------------------------------------------------------
  lx = length(x)
  S = Array{Union{Missing,T}}(missing, _N(lambda, lambda), lx)
  sch(lx, 1, lambda)
end

# ------------------------------------------------------------------------------
#~~ Symbolic Schur polynomial ~~~~##
# ------------------------------------------------------------------------------
function SchurPolynomial0(
  m::I,
  lambda::Vector{I},
  T::Type
) where {I<:Integer}
  function sch(m::Int64, k::Int64, nu::Vector{I})
    if isempty(nu) || nu[1] == 0 || m == 0
      return T(1)
    end
    if length(nu) > m && nu[m+1] > 0
      return T(0)
    end
    if m == 1
      return x[1]^nu[1]
    end
    v = S[_N(lambda, nu), m]
    if !ismissing(v)
      return v
    end
    s = sch(m - 1, 1, nu)
    lnu = length(nu)
    i = k
    while lnu >= i && nu[i] > 0
      if lnu == i || nu[i] > nu[i+1]
        nup = copy(nu)
        nup[i] -= 1
        if nu[i] > 1
          s += x[m] * sch(m, i, nup)
        else
          s += x[m] * sch(m - 1, 1, nup)
        end
      end
      i += 1
    end
    if k == 1
      S[_N(lambda, lambda), m] = s
    end
    return s
  end # end sch --------------------------------------------------------------
  DynamicPolynomials.@polyvar x[1:m]
  S = Array{Union{Missing,DynamicPolynomials.Polynomial{true,T}}}(
    missing,
    _N(lambda, lambda),
    m
  )
  sch(m, 1, lambda)
end

"""
    SchurPolynomial(m, lambda[, type])

Symbolic Schur polynomial.

# Arguments
- `m`: integer, the number of variables
- `lambda`: partition of an integer
- `type`: the type of the coefficients of the polynomial; default `Int64`
"""
function SchurPolynomial(
  m::I,
  lambda::Vector{I},
  type::Type = Int64
) where {I<:Integer}
  if !isPartition(lambda)
    error("`lambda` must be a partition of an integer")
  end
  schur = SchurPolynomial0(m, lambda, type)
  if(typeof(schur) == type)
    DynamicPolynomials.@polyvar x[1:m]
    schur = sum(type(0) * x) + schur
  end
  return schur
end

end # module


text:
import CUDA: CuVector, CuMatrix

export CuBitVector, AbstractBitVector,
    chunks, _msk_end

"Custom CUDA version of BitVector (lacking lots of functionality, just a container for now)."
struct CuBitVector <: AbstractVector{Bool}
    chunks::CuVector{UInt64}
    len::Int
end

"Retro-fitted super type of all bit vectors"
const AbstractBitVector = Union{BitVector, CuBitVector}

AbstractBitVector(chunks::CuVector{UInt64}, len) = CuBitVector(chunks, len)
AbstractBitVector(chunks::Vector{UInt64}, len) = begin
    v = BitVector()
    v.chunks = chunks
    v.len = len
    # v.dims unused by vectors
    return v
end

import Base: length, _msk_end, getindex, size, unsafe_bitgetindex #extend
@inline length(v::CuBitVector) = v.len
@inline size(v::CuBitVector) = (length(v),)
@inline _msk_end(v::CuBitVector) = _msk_end(length(v))
@inline function getindex(B::CuBitVector, i::Int)
    Bc = B.chunks
    i1, i2 = Base.get_chunks_id(i)
    u = UInt64(1) << i2
    r = (Bc[i1] & u) != 0
    return r
end


text: plot(randn(10, 3), colors=["violet", "navyblue", "orange"], lwidth=0.1)


text:
using MeshCat
using Base.Test
using GeometryTypes
using CoordinateTransformations
using Colors
using MeshIO, FileIO

@testset "MeshCat" begin
    include("video_rendering.jl")
    include("paths.jl")
    include("visualizer.jl")
    include("notebook.jl")

    @testset "deprecations" begin
        # Deprecated in v0.0.2
        @testset "cylinder" begin
            l = 1.0
            r = 2.0
            c = HyperCylinder(l, r)
            @test c.origin == Point(0., 0, 0)
            @test c.extremity == Point(0., 0, l)
            @test radius(c) == r
        end
    end
end

module ModuleTest
    # Test for any https://github.com/JuliaLang/julia/issues/21653
    # related issues when MeshCat is included in another module
    using MeshCat
end


text:
#############################################################################################

print("   running RosCtx test...")

before_time = time()

@context RosCtx

function rosenbrock(x::Vector{Float64})
    a = 1.0
    b = 100.0
    result = 0.0
    for i in 1:length(x)-1
        result += (a - x[i])^2 + b*(x[i+1] - x[i]^2)^2
    end
    return result
end

x = rand(2)
@inferred(overdub(RosCtx(), rosenbrock, x))

messages = String[]
Cassette.prehook(::RosCtx, f, args...) = push!(messages, string("calling ", f, args))
@test overdub(RosCtx(), rosenbrock, x) == rosenbrock(x)
@test length(messages) > 90

Cassette.prehook(ctx::RosCtx, f, args...) = push!(ctx.metadata, string("calling ", f, args))
messages2 = String[]
@test overdub(RosCtx(metadata=messages2), rosenbrock, x) == rosenbrock(x)
@test messages == messages2

Cassette.prehook(::RosCtx, f, args...) = nothing
Cassette.prehook(ctx::RosCtx, f, args::Number...) = push!(ctx.metadata, args)
argslog = Any[]
@test @overdub(RosCtx(metadata=argslog), rosenbrock(x)) == rosenbrock(x)
for args in argslog
    @test all(x -> isa(x, Number), args)
end

println("done (took ", time() - before_time, " seconds)")

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

print("   running HookCtx test...")

before_time = time()

@context HookCtx

Cassette.prehook(ctx::HookCtx, f, args...) = push!(ctx.metadata[1], (f, args))
Cassette.posthook(ctx::HookCtx, out, f, args...) = push!(ctx.metadata[2], (f, args))

ctx = HookCtx(metadata=(Any[], Any[]))
@overdub(ctx, 1 + 1 * 1)
# assumes the prehook trace is `[(*, ...), (mul_int, ...), (+, ...), (add_int, ...)]`
@test ctx.metadata[1][1] == ctx.metadata[2][2]
@test ctx.metadata[1][2] == ctx.metadata[2][1]
@test ctx.metadata[1][3] == ctx.metadata[2][4]
@test ctx.metadata[1][4] == ctx.metadata[2][3]

Cassette.execute(::HookCtx, ::typeof(+), args...) = +(args...)
empty!(ctx.metadata[1])
empty!(ctx.metadata[2])

@overdub(ctx, 1 + 1 * 1)
# assumes the prehook trace is `[(*, ...), (mul_int, ...), (+, ...)]`
@test ctx.metadata[1][1] == ctx.metadata[2][2]
@test ctx.metadata[1][2] == ctx.metadata[2][1]
@test ctx.metadata[1][3] == ctx.metadata[2][3]

println("done (took ", time() - before_time, " seconds)")

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

print("   running SinCtx test...")

before_time = time()

@context SinCtx
x = rand()
sin_plus_cos(x) = sin(x) + cos(x)
@test @overdub(SinCtx(), sin_plus_cos(x)) === sin_plus_cos(x)
Cassette.execute(::SinCtx, ::typeof(sin), x) = cos(x)
@test @overdub(SinCtx(), sin_plus_cos(x)) === (2 * cos(x))

println("done (took ", time() - before_time, " seconds)")

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

print("   running FoldCtx test...")

before_time = time()

@context FoldCtx
x = 2
foldmul(x, args...) = Core._apply(Base.afoldl, (*, x), args...)
@test @overdub(FoldCtx(), foldmul(x)) === foldmul(x)

println("done (took ", time() - before_time, " seconds)")

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

print("   running CountCtx test...")

before_time = time()

@context CountCtx
count1 = Ref(0)
Cassette.prehook(ctx::CountCtx, f, args::Number...) = (ctx.metadata[] += 1)
@overdub(CountCtx(metadata=count1), sin(1))
Cassette.prehook(ctx::CountCtx, f, args::Number...) = (ctx.metadata[] += 2)
count2 = Ref(0)
@overdub(CountCtx(metadata=count2), sin(1))
@test (2 * count1[]) === count2[]

println("done (took ", time() - before_time, " seconds)")

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

print("   running CountCtx2 test...")

before_time = time()

@context CountCtx2
mutable struct Count{T}
    count::Int
end
function Cassette.prehook(ctx::CountCtx2{Count{T}}, f, arg::T, args::T...) where {T}
    ctx.metadata.count += 1
end
mapstr(x) = map(string, x)
c = Count{Union{String,Int}}(0)
@test @overdub(CountCtx2(metadata=c), mapstr(1:10)) == mapstr(1:10)
@test c.count > 1000

println("done (took ", time() - before_time, " seconds)")

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

print("   running SqrCtx test...")

before_time = time()

@context SqrCtx
square_closure(x) = (y -> y * x)(x)
x = rand()
@test square_closure(x) == @overdub(SqrCtx(), square_closure(x))

println("done (took ", time() - before_time, " seconds)")

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

print("   running CompCtx test...")

before_time = time()

@context CompCtx
comprehension1(x) = [i for i in x]
comprehension2(f, x, y) = [f(x, i) for i in y]
f, x, y = hypot, rand(), rand(2)
@test comprehension1(x) == @overdub(CompCtx(), comprehension1(x))
@test comprehension1(y) == @overdub(CompCtx(), comprehension1(y))
@test comprehension2(f, x, y) == @overdub(CompCtx(), comprehension2(f, x, y))

println("done (took ", time() - before_time, " seconds)")

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

print("   running PassCtx test...")

before_time = time()

@context PassCtx
sig_collection = DataType[]
mypass = @pass (ctx, sig, cinfo) -> begin
    # TODO: check !in(sig, sig_collection) first to mimic a caching technique that will be used in practice
    push!(sig_collection, sig)
    return cinfo
end
@overdub(PassCtx(pass=mypass), sum(rand(3)))
@test !isempty(sig_collection) && all(T -> T <: Tuple, sig_collection)

println("done (took ", time() - before_time, " seconds)")

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

print("   running WorldCtx test...")

before_time = time()

@context WorldCtx
worldtest = 0
oldctx = WorldCtx()
overdub(oldctx, sin, 1)

Cassette.prehook(::WorldCtx, args...) = (global worldtest += 1)
overdub(WorldCtx(), sin, 1)
@test worldtest > 100

tmp = worldtest
overdub(oldctx, sin, 1)
@test tmp < worldtest

tmp = worldtest
Cassette.prehook(::WorldCtx, args...) = nothing
overdub(WorldCtx(), sin, 1)
@test tmp === worldtest

println("done (took ", time() - before_time, " seconds)")

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

print("   running TraceCtx test...")

before_time = time()

x, y, z = rand(3)
trtest(x, y, z) = x*y + y*z
trkwtest(x; _y = 1.0, _z = 2.0) = trtest(x, _y, _z)

@context TraceCtx

function Cassette.execute(ctx::TraceCtx, args...)
    subtrace = Any[]
    push!(ctx.metadata, args => subtrace)
    if canoverdub(ctx, args...)
        newctx = similarcontext(ctx, metadata = subtrace)
        return overdub(newctx, args...)
    else
        return fallback(ctx, args...)
    end
end
trace = Any[]
@test @overdub(TraceCtx(metadata = trace), trtest(x, y, z)) == trtest(x, y, z)
@test trace == Any[
    (trtest,x,y,z) => Any[
        (*,x,y) => Any[(Base.mul_float,x,y)=>Any[]]
        (*,y,z) => Any[(Base.mul_float,y,z)=>Any[]]
        (+,x*y,y*z) => Any[(Base.add_float,x*y,y*z)=>Any[]]
    ]
]

# jrevels/Cassette.jl#48
@context HookTraceCtx

mutable struct HookTrace
    current::Vector{Any}
    stack::Vector{Any}
    HookTrace() = new(Any[], Any[])
end

tracekw = Any[]
@overdub(TraceCtx(metadata = tracekw), trkwtest(x, _y = y, _z = z)) == trtest(x, y, z)
subtracekw = first(Iterators.filter(t -> t[1] === (Core.kwfunc(trkwtest), (_y = y, _z = z), trkwtest, x), tracekw))[2]
@test subtracekw == trace

function enter!(t::HookTrace, args...)
    pair = args => Any[]
    push!(t.current, pair)
    push!(t.stack, t.current)
    t.current = pair.second
    return nothing
end

function exit!(t::HookTrace)
    t.current = pop!(t.stack)
    return nothing
end

Cassette.prehook(ctx::HookTraceCtx, args...) = enter!(ctx.metadata, args...)
Cassette.posthook(ctx::HookTraceCtx, args...) = exit!(ctx.metadata)

htrace = HookTrace()
@overdub(HookTraceCtx(metadata = htrace), trtest(x, y, z))
@test htrace.current == trace

println("done (took ", time() - before_time, " seconds)")

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

print("   running NestedReflectCtx test...")

before_time = time()

@context NestedReflectCtx
r_pre = Cassette.reflect((typeof(sin), Int))
r_post = Cassette.reflect((typeof(overdub), typeof(NestedReflectCtx()), typeof(sin), Int))
@test isa(r_pre, Cassette.Reflection) && isa(r_post, Cassette.Reflection)
Cassette.overdub_pass!(r_pre, typeof(NestedReflectCtx()))
@test r_pre.code_info.code == r_post.code_info.code

println("done (took ", time() - before_time, " seconds)")

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

print("   running GemvCtx test...")

before_time = time()

# issue #51
@context GemvCtx
α, β = 1.0, 2.0
A, X, Y = rand(Float64, 10, 10), rand(Float64, 10), rand(Float64, 10)
Y_copy = copy(Y)
Y_out = overdub(GemvCtx(), LinearAlgebra.BLAS.gemv!, 'T', α, A, X, β, Y)
Y_copy_out = LinearAlgebra.BLAS.gemv!('T', α, A, X, β, Y_copy)
@test Y_out == Y_copy_out
@test Y_out === Y
@test Y_copy_out === Y_copy

println("done (took ", time() - before_time, " seconds)")

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

print("   running InferCtx test...")

before_time = time()

@context InferCtx

dispatchtupletest(::Type{T}) where {T} = Base.isdispatchtuple(Tuple{T}) ? T : Any
relu(x) = max(zero(x), x)
relulayer(W, x, b) = relu.(W*x .+ b)

@inferred(overdub(InferCtx(), typejoin, Float32, Float32, Float32))
@inferred(overdub(InferCtx(), dispatchtupletest, Float32))
@inferred(overdub(InferCtx(), (a, b) -> Core.apply_type(a, b), AbstractVector, Int))
@inferred(overdub(InferCtx(), eltype, rand(1)))
@inferred(overdub(InferCtx(), *, rand(1, 1), rand(1, 1)))
@inferred(overdub(InferCtx(), *, rand(Float32, 1, 1), rand(Float32, 1, 1)))
@inferred(overdub(InferCtx(), *, rand(Float32, 1, 1), rand(Float32, 1)))
@inferred(overdub(InferCtx(), rand, Float32, 1))
@inferred(overdub(InferCtx(), broadcast, +, rand(1), rand(1)))
@inferred(overdub(InferCtx(), relulayer, rand(Float64, 1, 1), rand(Float32, 1), rand(Float32, 1)))

println("done (took ", time() - before_time, " seconds)")

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

print("   running InModuleCtx test...")

before_time = time()

module DefineStuffInModule
    using Cassette
    Cassette.@context InModuleCtx
    x = 0
    incrpass = Cassette.@pass (ctx, sig, code) -> (global x += 1; code)
    f(x) = Cassette.overdub(InModuleCtx(pass = incrpass), sin, x)
end

@test DefineStuffInModule.f(1) === sin(1)
@test DefineStuffInModule.x > 0

println("done (took ", time() - before_time, " seconds)")

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

print("   running InvokeCtx test...")

before_time = time()

invoked(x::Int) = x + x
invoked(x::Number) = x^2
invoker(x) = invoke(invoked, Tuple{Number}, x)

@context InvokeCtx
Cassette.prehook(ctx::InvokeCtx, f, args...) = push!(ctx.metadata, f)
ctx = InvokeCtx(metadata=Any[])

@test overdub(ctx, invoker, 3) === 9
# This is kind of fragile and may break for unrelated reasons - the main thing
# we're testing here is that we properly trace through the `invoke` call.
@test ctx.metadata == Any[Core.apply_type, Core.invoke, Core.apply_type,
                          Val{2}, Core.apply_type, Base.literal_pow, *,
                          Base.mul_int]

println("done (took ", time() - before_time, " seconds)")


text:
immutable OffsetArray{T,N} <: AbstractArray{T,N}
    data::Array{T,N}
    startOffset::Dims{N}
end

Base.eltype(A::OffsetArray) = eltype(A.data)
Base.indices(A::OffsetArray) = map((x,y)->(x-1)+y, A.startOffset, indices(A.data))

@inline Base.getindex(A::OffsetArray, idxs...) = A.data[map((x,y) -> (x-y+1),idxs,A.startOffset)...]
@inline Base.setindex!(A::OffsetArray, X, idxs...) = A.data[map((x,y) -> (x-y+1),idxs,A.startOffset)...] = X

OffsetArray{T}(::Type{T},startOffset::Dims,dims::Dims) = OffsetArray(Array{T}(dims),startOffset)
OffsetArray{T}(::Type{T},startOffset::Dims,dims::Integer...) = OffsetArray(Array{T}(dims...),startOffset)
OffsetArray(startOffset::Dims,dims::Dims) = OffsetArray(Float64,startOffset,dims)
OffsetArray(startOffset::Dims,dims::Integer...) = OffsetArray(Float64,startOffset,dims...)
offsetarray_fill{N}(v,startOffset::Dims{N},dims::Dims{N}) = OffsetArray(fill(v,dims),startOffset)

offsetarray_fill{N}(v::AbstractArray,startOffset::Dims{N},dims::Dims{N}) = OffsetArray(map(x->copy(v), zeros(Bool,dims)),startOffset)
offsetarray_fill{N}(v::AbstractArray,startOffset::Dims{N},dims::Integer...) = OffsetArray(map(x->copy(v), zeros(Bool,dims...)),startOffset)

offsetarray_zeros{T,N}(t::Type{T},startOffset::Dims{N},dims::Dims{N}) = OffsetArray(zeros(t,dims),startOffset)
offsetarray_ones{T,N}(t::Type{T},startOffset::Dims{N},dims::Integer...) = OffsetArray(ones(t,dims...),startOffset)
offsetarray_zeros{N}(startOffset::Dims{N},dims::Dims{N}) = OffsetArray(zeros(dims),startOffset)
offsetarray_ones{N}(startOffset::Dims{N},dims::Integer...) = OffsetArray(ones(dims...),startOffset)

immutable AccelGrid{T,N}
    gridSize::Float64
    data::OffsetArray{T,N}
end

@inline Base.getindex(A::AccelGrid,idxs...) = A.data[map((x,y)->clamp(round(Int,x/A.gridSize),minimum(y),maximum(y)),idxs,indices(A.data))...]
@inline Base.setindex!(A::AccelGrid,X,idxs...) = A.data[map((x,y)->clamp(round(Int,x/A.gridSize),minimum(y),maximum(y)),idxs,indices(A.data))...] = X
Base.eltype(A::AccelGrid) = eltype(A.data)

function AccelGrid{T,N}(::Type{T}, gridSize::Real,lowerBound::NTuple{N,Real}, upperBound::NTuple{N,Real})
    startBound = map(x->floor(Int,x/gridSize),lowerBound)
    endBound = map(x->ceil(Int,x/gridSize),upperBound)

    if T <: AbstractArray
        eT = eltype(T)
        data = offsetarray_fill(eT[],startBound,map(-,endBound,startBound))
    elseif T <: Real
        data = offsetarray_zeros(T,startBound,map(-,endBound,startBound))
    else
        data = offsetarray_fill(Nullable{T}(),startBound,map(-,endBound,startBound))
    end
    AccelGrid(convert(Float64,gridSize),data)
end

AccelGrid{N}(gridSize::Real,lowerBound::NTuple{N,Real}, upperBound::NTuple{N,Real}) = AccelGrid(Float64,gridSize,lowerBound,upperBound)

# function Base.display(a::AccelGrid)
#   bold = "\x1b[1m"; default = "\x1b[0m";
#   println("$(bold)AccelGrid:")
#   println("\tgrid size: ", a.gridSize)
#   print("\tdata: $(default)")
#   display(a.data)
# end


text:
module TestStaticArrays

using Test
using BangBang: ismutable
using StaticArrays

@testset begin
    @test !ismutable(SVector(0))
    @test !ismutable(SMatrix{1, 1}(0))
    @test ismutable(MVector(0))
    @test ismutable(MMatrix{1, 1}(0))
end

end  # module


text:
module MoireHamiltonian

using Plots

include("lattice.jl")
include("folding.jl")
include("bilayer.jl")

end


text:
using ITensors, Test
using Combinatorics: permutations

@testset "Combiner" begin

i = Index(2,"i")
j = Index(3,"j")
k = Index(4,"k")
l = Index(5,"l")

A = randomITensor(i, j, k, l)

@testset "Basic combiner properties" begin
  C = combiner(i, j, k)
  @test eltype(storage(C)) === Number
  @test_throws ErrorException ITensors.data(C)
  @test NDTensors.uncombinedinds(NDTensors.tensor(C)) == (i, j, k)
end

@testset "Empty combiner" begin
  C = combiner()
  @test order(C) == 0
  @test isnothing(combinedind(C))
  AC = A * C
  @test A == AC
  AC = C * A
  @test A == AC

  R = ITensor(j, l, k, i)
  R .= A .* C
  @test R == A
end

@testset "Two index combiner" begin
  for inds_ij ∈ permutations([i,j])
    C = combiner(inds_ij...)
    c = combinedind(C)
    B = A*C
    @test hasinds(B, l, k, c)
    @test c == commonind(B, C)
    @test combinedind(C) == c
    @test isnothing(combinedind(A))
    D = B*C
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    @test isnothing(combinedind(D))
  end

  for inds_il ∈ permutations([i,l])
    C = combiner(inds_il...)
    c = combinedind(C)
    B = A*C
    @test hasinds(B, j, k)
    @test c == commonind(B, C)
    D = B*C
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
  end

  for inds_ik ∈ permutations([i,k])
    C = combiner(inds_ik...)
    c = combinedind(C)
    B = A*C
    @test hasinds(B, j, l)
    @test c == commonind(B, C)
    D = B*C
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
  end

  for inds_jk ∈ permutations([j,k])
    C = combiner(inds_jk...)
    c = combinedind(C)
    B = A*C
    @test hasinds(B, i, l)
    @test c == commonind(B, C)
    D = B*C
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    D = C*B
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    B = C*A
    @test hasinds(B, i, l)
    @test c == commonind(B, C)
    D = B*C
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    D = C*B
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
  end

  for inds_jl ∈ permutations([j,l])
    C = combiner(inds_jl...)
    c = combinedind(C)
    B = A*C
    @test hasinds(B, i, k)
    @test c == commonind(B, C)
    D = B*C
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    D = C*B
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    B = C*A
    @test hasinds(B, i, k)
    @test c == commonind(B, C)
    D = B*C
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    D = C*B
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
  end

  for inds_kl ∈ permutations([k,l])
    C = combiner(inds_kl...)
    c = combinedind(C)
    B = A*C
    @test hasinds(B, i, j)
    @test c == commonind(B, C)
    D = B*C
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    D = C*B
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    B = C*A
    @test hasinds(B, i, j)
    @test c == commonind(B, C)
    D = B*C
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    D = C*B
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
  end
end

@testset "Three index combiner" begin
  for inds_ijl ∈ permutations([i,j,l])
    C = combiner(inds_ijl...)
    c = combinedind(C)
    B = A*C
    @test hasind(B, k)
    @test c == commonind(B, C)
    D = B*C
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    D = C*B
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    B = C*A
    @test hasind(B, k)
    @test c == commonind(B, C)
    D = B*C
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    D = C*B
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
  end

  for inds_ijk ∈ permutations([i,j,k])
    C = combiner(inds_ijk...)
    c = combinedind(C)
    B = A*C
    @test hasind(B, l)
    @test c == commonind(B, C)
    D = B*C
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    D = C*B
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    B = C*A
    @test hasind(B, l)
    @test c == commonind(B, C)
    D = B*C
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    D = C*B
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
  end

  for inds_jkl ∈ permutations([j,k,l])
    C = combiner(inds_jkl...)
    c = combinedind(C)
    B = A*C
    @test hasind(B, i)
    @test c == commonind(B, C)
    D = B*C
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    D = C*B
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    B = C*A
    @test hasind(B, i)
    @test c == commonind(B, C)
    D = B*C
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
    D = C*B
    @test hasinds(D, i, j, k, l)
    @test D ≈ A
  end
end

@testset "SVD/Combiner should play nice" begin
  C = combiner(i, j, k)
  c = combinedind(C)
  Ac = A*C
  U,S,V,spec,u,v = svd(Ac, c)
  Uc = C*U
  Ua,Sa,Va,spec,ua,va = svd(A, i, j, k)
  replaceind!(Ua, ua, u)
  @test A ≈ C*Ac
  @test A ≈ Ac*C
  @test Ua*C ≈ U
  @test C*Ua ≈ U
  @test Ua ≈ Uc
  @test Uc*S*V ≈ A
  @test (C*Ua)*S*V ≈ Ac
  C = combiner(i, j)
  c = combinedind(C)
  Ac = A*C
  U,S,V,spec,u,v = svd(Ac, c)
  Uc = U*C
  Ua,Sa,Va,spec,ua,va = svd(A, i, j)
  replaceind!(Ua, ua, u)
  @test Ua ≈ Uc
  @test Ua*C ≈ U
  @test C*Ua ≈ U
  @test Uc*S*V ≈ A
  @test (C*Ua)*S*V ≈ Ac
end

@testset "mult/Combiner should play nice" begin
  C = combiner(i, j, k)
  Ac = A*C
  B = randomITensor(l)
  AB = Ac*B
  @test AB*C ≈ A*B
end

@testset "Replace index combiner" begin
  C = combiner(l, tags="nl")
  c = combinedind(C)
  B = A*C
  replaceind!(B, c, l)
  @test B == A
end

end

nothing


text:
# Extra file for fast-fixes for packages ReactiveMP.jl depends on

# Hotfix for https://github.com/JuliaStats/Distributions.jl/issues/1270
function score(::DifferentialEntropy, marginal::Marginal{ <: Dirichlet })
    d  = getdata(marginal)
    α  = d.alpha
    α0 = d.alpha0
    k  = length(α)
    en = d.lmnB + (α0 - k) * digamma(α0) - sum(αj -> (αj - 1) * digamma(αj), α)
end