text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ADIOS2_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ADIOS2")
JLLWrappers.@generate_main_file("ADIOS2", UUID("44b1415a-b224-5c99-9168-ff4febb5f37f"))
end  # module ADIOS2_jll


text:
using MVApp
using Test
using GeometryBasics
using Statistics: mean, var, std

include("test_utils_movement.jl")
include("test_smoothing.jl")
include("test_shape.jl")

data_folder = "../data"
@testset "MovingFrame.jl" begin
    test_standardize()
    test_gradient()
    test_smoothing()
    test_shape()
end

ext: jl


text:
function Fract_Maxwell_G_reduce()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fm_model = Fract_Maxwell.Ga(t, [1.0, 1.0, 1.0, 0.0])
    maxwell = Maxwell.Ga(t, [1.0, 1.0])

    all(i -> isapprox(fm_model[i], maxwell[i]), eachindex(t))
end
@test Fract_Maxwell_G_reduce()

function FractS_Maxwell_G_reduce()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fms_model = FractS_Maxwell.Ga(t, [1.0, 1.0, 1.0])
    maxwell = Maxwell.Ga(t, [1.0, 1.0])

    all(i -> isapprox(fms_model[i], maxwell[i]), eachindex(t))
end
@test FractS_Maxwell_G_reduce()

function FractD_Maxwell_G_reduce()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fmd_model = FractD_Maxwell.Ga(t, [1.0, 1.0, 0.0])
    maxwell = Maxwell.Ga(t, [1.0, 1.0])

    all(i -> isapprox(fmd_model[i], maxwell[i]), eachindex(t))
end
@test FractD_Maxwell_G_reduce()

function Fract_Maxwell_J_reduce()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fm_model = Fract_Maxwell.Ja(t, [1.0, 1.0, 1.0, 0.0])
    maxwell = Maxwell.Ja(t, [1.0, 1.0])

    all(i -> isapprox(fm_model[i], maxwell[i]), eachindex(t))
end
@test Fract_Maxwell_J_reduce()

function FractS_Maxwell_J_reduce()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fms_model = FractS_Maxwell.Ja(t, [1.0, 1.0, 1.0])
    maxwell = Maxwell.Ja(t, [1.0, 1.0])

    all(i -> isapprox(fms_model[i], maxwell[i]), eachindex(t))
end
@test FractS_Maxwell_J_reduce()

function FractD_Maxwell_J_reduce()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fmd_model = FractD_Maxwell.Ja(t, [1.0, 1.0, 0.0])
    maxwell = Maxwell.Ja(t, [1.0, 1.0])

    all(i -> isapprox(fmd_model[i], maxwell[i]), eachindex(t))
end
@test FractD_Maxwell_J_reduce()

function Fract_Maxwell_Gp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fm_model = Fract_Maxwell.Gpa(chirp, [1.0, 1.0, 1.0, 0.0])
    maxwell = Maxwell.Gpa(chirp, [1.0, 1.0])

    all(i -> isapprox(fm_model[i], maxwell[i]), eachindex(chirp))
end
@test Fract_Maxwell_Gp_reduce()

function FractS_Maxwell_Gp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fms_model = FractS_Maxwell.Gpa(chirp, [1.0, 1.0, 1.0])
    maxwell = Maxwell.Gpa(chirp, [1.0, 1.0])

    all(i -> isapprox(fms_model[i], maxwell[i]), eachindex(chirp))
end
@test FractS_Maxwell_Gp_reduce()

function FractD_Maxwell_Gp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fmd_model = FractD_Maxwell.Gpa(chirp, [1.0, 1.0, 0.0])
    maxwell = Maxwell.Gpa(chirp, [1.0, 1.0])

    all(i -> isapprox(fmd_model[i], maxwell[i]), eachindex(chirp))
end
@test FractD_Maxwell_Gp_reduce()

function Fract_Maxwell_Gp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fm_model = Fract_Maxwell.Gppa(chirp, [1.0, 1.0, 1.0, 0.0])
    maxwell = Maxwell.Gppa(chirp, [1.0, 1.0])

    all(i -> isapprox(fm_model[i], maxwell[i]), eachindex(chirp))
end
@test Fract_Maxwell_Gp_reduce()

function FractS_Maxwell_Gp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fms_model = FractS_Maxwell.Gppa(chirp, [1.0, 1.0, 1.0])
    maxwell = Maxwell.Gppa(chirp, [1.0, 1.0])

    all(i -> isapprox(fms_model[i], maxwell[i]), eachindex(chirp))
end
@test FractS_Maxwell_Gp_reduce()

function FractD_Maxwell_Gp_reduce()
    ω_step = 0.01;
    chirp = collect(0.0:ω_step:10000.0);

    fmd_model = FractD_Maxwell.Gppa(chirp, [1.0, 1.0, 0.0])
    maxwell = Maxwell.Gppa(chirp, [1.0, 1.0])

    all(i -> isapprox(fmd_model[i], maxwell[i]), eachindex(chirp))
end
@test FractD_Maxwell_Gp_reduce()


text:
#= iterative improvement to a solution of a linear system A x = b

cf. Numerical Recipes (2007), sec 2.5

this is the same procedure as gsl_linalg_LU_refine:
https://github.com/ampl/gsl/blob/master/linalg/lu.c

the algorithm is:

1. compute the residual = A * x - b

2. solve for the correction delta, A * delta = residual

3. subtract this from the old solution for the improved solution
   xnew = x - delta

Note: in Numerical Recipes it is suggested to use higher precision to compute
the residual, since this can be very small. from my (very limited) testing, this
didn't seem crucial, and i'm also not sure what would be the best way to do that
here. the issue is that we just want to compute the residual in higher
precision, not performing the left division. so, if we were to expect the caller
to convert to higher precision upon calling this function, we'd still have to
convert back to lower precision before calling ldiv! (for efficiency reasons and
also since most ldiv methods don't work for higher precision types anyway). so
for now it stays like this.
=#
function refine_solution!(sol::Vector, A_fact, A_mat::AbstractMatrix, b_vec::Vector)
    mul!(b_vec, A_mat, sol, 1, -1)
    ldiv!(A_fact, b_vec)
    sol .-= b_vec
end


text:
using CSV, DataFrames
using Random, Statistics
using ChemistryFeaturization
using AtomicGraphNets
using Formatting
using CalibrationErrors

round4(x) = round(x, digits = 4)

function get_data(num_pts = 100, data_dir = joinpath(@__DIR__, "data"), verbose = true)
    println("Setting things up...")

    # data-related options
    train_frac = 0.8 # what fraction for training?
    num_train = Int32(round(train_frac * num_pts))
    num_test = num_pts - num_train
    prop = "formation_energy_per_atom"
    id = "task_id" # field by which to label each input material

    # set up the featurization
    featurization = GraphNodeFeaturization([
        "Group",
        "Row",
        "Block",
        "Atomic mass",
        "Atomic radius",
        "X",
    ])
    num_features =
        sum(ChemistryFeaturization.FeatureDescriptor.output_shape.(featurization.features)) # TODO: update this with cleaner syntax once new version of CF is tagged that has it
    # dataset...first, read in outputs
    info = CSV.read(string(data_dir, "/", prop, ".csv"), DataFrame)
    y = Array(Float32.(info[!, Symbol(prop)]))

    # shuffle data and pick out subset
    indices = shuffle(1:size(info, 1))[1:num_pts]
    info = info[indices, :]
    output = y[indices]

    # next, make and featurize graphs
    if verbose
        println("Building graphs and feature vectors from structures...")
    end
    inputs = FeaturizedAtoms[]

    # for r in eachrow(info)
    cifpaths = [
        joinpath(data_dir, format("{}_cifs", prop), string(r[Symbol(id)], ".cif")) for
        r in eachrow(info)
    ]
    graphs = skipmissing(AtomGraph.(cifpaths))
    inputs = [featurize(gr, featurization) for gr in graphs]

    # pick out train/test sets
    if verbose
        println("Dividing into train/test sets...")
    end
    train_output = output[1:num_train]
    test_output = output[num_train+1:end]
    train_input = inputs[1:num_train]
    test_input = inputs[num_train+1:end]
    train_data = zip(train_input, train_output)
    test_data = zip(test_input, test_output)

    # build the model
    if verbose
        println("Building the network...")
    end

    return num_features, train_input, train_output, test_input, test_output
end


text:
module Cthulhu

using TerminalMenus
using InteractiveUtils

using Core: MethodInstance
const Compiler = Core.Compiler

include("callsite.jl")

export descend, @descend, descend_code_typed, descend_code_warntype, @descend_code_typed, @descend_code_warntype

"""
  @descend_code_typed

  Evaluates the arguments to the function or macro call, determines their
types, and calls `code_typed` on the resulting expression.
"""
macro descend_code_typed(ex0...)
    InteractiveUtils.gen_call_with_extracted_types_and_kwargs(__module__, :descend_code_typed, ex0)
end

"""
  @descend_code_warntype

  Evaluates the arguments to the function or macro call, determines their
types, and calls `code_warntype` on the resulting expression.
"""
macro descend_code_warntype(ex0...)
    InteractiveUtils.gen_call_with_extracted_types_and_kwargs(__module__, :descend_code_warntype, ex0)
end

"""
  @descend

  Shortcut for [`@descend_code_typed`](@ref).
"""
macro descend(ex0...)
    esc(:(@descend_code_typed($(ex0...))))
end

"""
    descend_code_typed(f, tt; kwargs...)

Given a function and a tuple-type, interactively explore the output of
`code_typed` by descending into `invoke` statements. Type enter to select an
`invoke` to descend into, select ↩  to ascend, and press q or control-c to
quit.

# Usage:
```julia
function foo()
    T = rand() > 0.5 ? Int64 : Float64
    sum(rand(T, 100))
end

descend_code_typed(foo, Tuple{})
```
"""
descend_code_typed(f, @nospecialize(tt); kwargs...) =
    _descend_with_error_handling(f, tt; iswarn=false, kwargs...)

"""
    descend_code_warntype(f, tt)

Given a function and a tuple-type, interactively explore the output of
`code_warntype` by descending into `invoke` statements. Type enter to select an
`invoke` to descend into, select ↩  to ascend, and press q or control-c to
quit.

# Usage:
```julia
function foo()
    T = rand() > 0.5 ? Int64 : Float64
    sum(rand(T, 100))
end

descend_code_warntype(foo, Tuple{})
```
"""
descend_code_warntype(f, @nospecialize(tt)) =
    _descend_with_error_handling(f, tt; iswarn=true)

function _descend_with_error_handling(f, @nospecialize(tt); kwargs...)
    try
        _descend(f, tt; kwargs...)
    catch x
        if x isa InterruptException
            return nothing
        else
            rethrow(x)
        end
    end
    return nothing
end

unwrap_type(T) = T
unwrap_type(T::Core.Compiler.Const) = typeof(T.val)

function find_callsites(CI, mi, slottypes; params=current_params(), kwargs...)
    spvals = Core.Compiler.spvals_from_meth_instance(mi)
    callsites = Callsite[]
    for (id, c) in enumerate(CI.code)
        if c isa Expr
            callsite = nothing
            if c.head === :invoke
                rt = CI.ssavaluetypes[id]
                callsite = Callsite(id, c.args[1], rt)
            elseif c.head === :call
                rt = CI.ssavaluetypes[id]
                types = map(arg -> unwrap_type(Compiler.argextype(arg, CI, spvals, slottypes)), c.args)

                # Filter out builtin functions and intrinsic function
                if types[1] <: Core.Builtin || types[1] <: Core.IntrinsicFunction
                    continue
                end

                # Filter out abstract signatures
                # otherwise generated functions get crumpy
                if any(isabstracttype, types) || any(T->(T isa Union || T isa UnionAll), types)
                    continue
                end

                mi = first_method_instance(Tuple{types...})
                mi == nothing && continue
                callsite = Callsite(id, mi, rt)
            end

            if callsite !== nothing
                push!(callsites, callsite)
            end
        end
    end
    return callsites
end

if VERSION >= v"1.1.0-DEV.215"
function dce!(ci, mi)
    argtypes = Core.Compiler.matching_cache_argtypes(mi, nothing)[1]
    ir = Compiler.inflate_ir(ci, Core.Compiler.spvals_from_meth_instance(mi),
                             argtypes)
    compact = Core.Compiler.IncrementalCompact(ir, true)
    # Just run through the iterator without any processing
    Core.Compiler.foreach(x -> nothing, compact)
    ir = Core.Compiler.finish(compact)
    Core.Compiler.replace_code_newstyle!(ci, ir, length(argtypes)-1)
end
else
function dce!(ci, mi)
end
end

"""
  descend

  Shortcut for [`descend_code_typed`](@ref).
"""
const descend = descend_code_typed

function show_as_line(el)
    reduced_displaysize = displaysize(stdout) .- (0, 3)
    buf = IOBuffer()
    show(IOContext(buf, :limit=>true, :displasize=>reduced_displaysize), el)
    String(take!(buf))
end

function do_typeinf_slottypes(mi::Core.Compiler.MethodInstance, run_optimizer::Bool, params::Core.Compiler.Params)
    ccall(:jl_typeinf_begin, Cvoid, ())
    result = Core.Compiler.InferenceResult(mi)
    frame = Core.Compiler.InferenceState(result, false, params)
    frame === nothing && return (nothing, Any)
    if Compiler.typeinf(frame) && run_optimizer
        opt = Compiler.OptimizationState(frame)
        Compiler.optimize(opt, result.result)
        opt.src.inferred = true
    end
    ccall(:jl_typeinf_end, Cvoid, ())
    frame.inferred || return (nothing, Any)
    return (frame.src, result.result, frame.slottypes)
end

function preprocess_ci!(ci, mi, optimize)
    if optimize
        # if the optimizer hasn't run, the IR hasn't been converted
        # to SSA form yet and dce is not legal
        dce!(ci, mi)
        dce!(ci, mi)
    end
end

current_params() = Core.Compiler.Params(ccall(:jl_get_world_counter, UInt, ()))
function _descend(mi::MethodInstance; iswarn::Bool, params=current_params(), optimize::Bool=true, kwargs...)
    (CI, rt, slottypes) = do_typeinf_slottypes(mi, optimize, params)

    callsites = find_callsites(CI, mi, slottypes; params=params, kwargs...)
    while true
        println()
        println("│ ─ $(string(Callsite(-1, mi, rt)))")
        iswarn ? code_warntype(F, TT) : display(CI=>rt)
        println()
        TerminalMenus.config(cursor = '•', scroll = :wrap)
        menu = RadioMenu(vcat(map(show_as_line, callsites), ["↩ "]))
        println("In `$(mi.def.name)` select a call to descend into or ↩ to ascend. [q] to quit.")
        cid = request(menu)
        if cid == length(callsites) + 1
            break
        end
        if cid == -1
            throw(InterruptException())
        end
        callsite = callsites[cid]
        _descend(callsite.mi; iswarn=iswarn, kwargs...)
    end
end

function first_method_instance(F, TT; params=current_params())
    sig = Tuple{typeof(F), TT.parameters...}
    first_method_instance(sig; params=params)
end

function first_method_instance(sig; params=current_params())
    methds = Base._methods_by_ftype(sig, 1, params.world)
    methds === false && return nothing
    x = methds[1]
    meth = x[3]
    if isdefined(meth, :generator) && !isdispatchtuple(Tuple{sig.parameters[2:end]...})
        return nothing
    end
    mi = Compiler.code_for_method(meth, sig, x[2], params.world)
end

function _descend(@nospecialize(F), @nospecialize(TT); params=current_params(), kwargs...)
    mi = first_method_instance(F, TT; params=params)
    _descend(mi; params=params, kwargs...)
end

end


text:
module Allocs

using Base.StackTraces: StackTrace, StackFrame, lookup
using Base: InterpreterIP

# --- Raw results structs, originally defined in C ---

# The C jl_bt_element_t object contains either an IP pointer (size_t) or a void*.
const BTElement = Csize_t;

# matches jl_raw_backtrace_t on the C side
struct RawBacktrace
    data::Ptr{BTElement} # in C: *jl_bt_element_t
    size::Csize_t
end

# matches jl_raw_alloc_t on the C side
struct RawAlloc
    type::Ptr{Type}
    backtrace::RawBacktrace
    size::Csize_t
end

# matches jl_profile_allocs_raw_results_t on the C side
struct RawResults
    allocs::Ptr{RawAlloc}
    num_allocs::Csize_t
end

"""
    Profile.Allocs.@profile [sample_rate=0.0001] expr

Profile allocations that happen during `expr`, returning
both the result and and AllocResults struct.

A sample rate of 1.0 will record everything; 0.0 will record nothing.

```julia
julia> Profile.Allocs.@profile sample_rate=0.01 peakflops()
1.03733270279065e11

julia> results = Profile.Allocs.fetch()

julia> last(sort(results.allocs, by=x->x.size))
Profile.Allocs.Alloc(Vector{Any}, Base.StackTraces.StackFrame[_new_array_ at array.c:127, ...], 5576)
```
"""
macro profile(opts, ex)
    _prof_expr(ex, opts)
end
macro profile(ex)
    _prof_expr(ex, :(sample_rate=0.0001))
end

# globals used for tracking how many allocs we're missing
# vs the alloc counters used by @time
const _g_gc_num_before = Ref{Base.GC_Num}()
const _g_sample_rate = Ref{Real}()
const _g_expected_sampled_allocs = Ref{Float64}(0)

function _prof_expr(expr, opts)
    quote
        $start(; $(esc(opts)))
        try
            $(esc(expr))
        finally
            $stop()
        end
    end
end

"""
    Profile.Allocs.start(sample_rate::Real)

Begin recording allocations with the given sample rate
A sample rate of 1.0 will record everything; 0.0 will record nothing.
"""
function start(; sample_rate::Real)
    ccall(:jl_start_alloc_profile, Cvoid, (Cdouble,), Float64(sample_rate))

    _g_sample_rate[] = sample_rate
    _g_gc_num_before[] = Base.gc_num()
end

"""
    Profile.Allocs.stop()

Stop recording allocations.
"""
function stop()
    ccall(:jl_stop_alloc_profile, Cvoid, ())

    # increment a counter of how many allocs we would expect
    # the memory profiler to see, based on how many allocs
    # actually happened.
    gc_num_after = Base.gc_num()
    gc_diff = Base.GC_Diff(gc_num_after, _g_gc_num_before[])
    alloc_count = Base.gc_alloc_count(gc_diff)
    expected_samples = alloc_count * _g_sample_rate[]
    _g_expected_sampled_allocs[] += expected_samples
end

"""
    Profile.Allocs.clear()

Clear all previously profiled allocation information from memory.
"""
function clear()
    ccall(:jl_free_alloc_profile, Cvoid, ())

    _g_expected_sampled_allocs[] = 0
    return nothing
end

"""
    Profile.Allocs.fetch()

Retrieve the recorded allocations, and decode them into Julia
objects which can be analyzed.
"""
function fetch()
    raw_results = ccall(:jl_fetch_alloc_profile, RawResults, ())
    decoded_results = decode(raw_results)

    # avoid divide-by-0 errors
    if _g_expected_sampled_allocs[] > 0
        missed_allocs = max(0, _g_expected_sampled_allocs[] - length(decoded_results.allocs))
        missed_percentage = max(0, round(Int, missed_allocs / _g_expected_sampled_allocs[] * 100))
        if missed_percentage > 0
            @warn("The allocation profiler is not fully implemented, and missed approximately" *
            " $(missed_percentage)% (estimated $(round(Int, missed_allocs)) / $(round(Int,
            _g_expected_sampled_allocs[]))) " *
                    "of sampled allocs in the last run. " *
                    "For more info see https://github.com/JuliaLang/julia/issues/43688")
        else
            @warn("The allocation profiler is not fully implemented, and may have missed" *
            " some of the allocs. " *
                    "For more info see https://github.com/JuliaLang/julia/issues/43688")
        end
    end
    return decoded_results
end

# decoded results

struct Alloc
    type::Any
    stacktrace::StackTrace
    size::Int
end

struct AllocResults
    allocs::Vector{Alloc}
end

# Without this, the Alloc's stacktrace prints for lines and lines and lines...
function Base.show(io::IO, a::Alloc)
    stacktrace_sample = length(a.stacktrace) >= 1 ? "$(a.stacktrace[1]), ..." : ""
    print(io, "$Alloc($(a.type), $StackFrame[$stacktrace_sample], $(a.size))")
end

const BacktraceCache = Dict{BTElement,Vector{StackFrame}}

# copied from julia_internal.h
const JL_BUFF_TAG = UInt(0x4eadc000)

struct CorruptType end
struct BufferType end

function load_type(ptr::Ptr{Type})
    if UInt(ptr) < UInt(4096)
        return CorruptType
    elseif UInt(ptr) == JL_BUFF_TAG
        return BufferType
    end
    return unsafe_pointer_to_objref(ptr)
end

function decode_alloc(cache::BacktraceCache, raw_alloc::RawAlloc)::Alloc
    Alloc(
        load_type(raw_alloc.type),
        stacktrace_memoized(cache, load_backtrace(raw_alloc.backtrace)),
        UInt(raw_alloc.size)
    )
end

function decode(raw_results::RawResults)::AllocResults
    cache = BacktraceCache()
    allocs = [
        decode_alloc(cache, unsafe_load(raw_results.allocs, i))
        for i in 1:raw_results.num_allocs
    ]
    return AllocResults(allocs)
end

function load_backtrace(trace::RawBacktrace)::Vector{BTElement}
    out = Vector{BTElement}()
    for i in 1:trace.size
        push!(out, unsafe_load(trace.data, i))
    end

    return out
end

function stacktrace_memoized(
    cache::BacktraceCache,
    trace::Vector{BTElement},
    c_funcs::Bool=true
)::StackTrace
    stack = StackTrace()
    for ip in trace
        frames = get(cache, ip) do
            res = lookup(ip)
            cache[ip] = res
            return res
        end
        for frame in frames
            # Skip frames that come from C calls.
            if c_funcs || !frame.from_c
                push!(stack, frame)
            end
        end
    end
    return stack
end

# Precompile once for the package cache.
@assert precompile(start, ())
@assert precompile(stop, ())
@assert precompile(fetch, ())

end

max_stars_repo_name: jipolanco/julia

max_issues_repo_name: jipolanco/julia

max_forks_repo_name: jipolanco/julia


text:
# T is the bitpattern type for logging, i is the logbook number
struct Sherlog16{T<:AbstractFloat,i} <: AbstractSherlog
    val::T
end

# conversions back from Sherlog16
Base.Float64(x::Sherlog16) = Float64(x.val)
Base.Float32(x::Sherlog16) = Float32(x.val)
Base.Float16(x::Sherlog16) = Float16(x.val)
Base.Int64(x::Sherlog16) = Int64(x.val)
Base.Int32(x::Sherlog16) = Int32(x.val)
Base.Int16(x::Sherlog16) = Int16(x.val)

# otherwise conversion to itself is ambigious
Sherlog16{T,i}(x::Sherlog16) where {T<:AbstractFloat,i} = Sherlog16{T,i}(x.val)
Sherlog16(x::Sherlog16{T,i}) where {T<:AbstractFloat,i} = Sherlog16{T,i}(x.val)

# generator functions from Float, Int and Bool
Sherlog16(x::AbstractFloat) = Sherlog16{Float16,1}(x)   # use Float16 and logbook #1 as standard
Sherlog16{T}(x::AbstractFloat) where T = Sherlog16{T,1}(x)      # use logbook #1 if no logbook provided
Sherlog16(x::Integer) = Sherlog16{Float16,1}(x)           # same for integers
Sherlog16{T}(x::Integer) where T = Sherlog16{T,1}(x)
Sherlog16{T}(x::Bool) where T = if x Sherlog16{T,1}(1) else Sherlog16{T,1}(0) end

Base.promote_rule(::Type{Int64},::Type{Sherlog16{T,i}}) where {T,i} = Sherlog16{T,i}
Base.promote_rule(::Type{Int32},::Type{Sherlog16{T,i}}) where {T,i} = Sherlog16{T,i}
Base.promote_rule(::Type{Int16},::Type{Sherlog16{T,i}}) where {T,i} = Sherlog16{T,i}

Base.promote_rule(::Type{Float64},::Type{Sherlog16{T,i}}) where {T,i} = Sherlog16{T,i}
Base.promote_rule(::Type{Float32},::Type{Sherlog16{T,i}}) where {T,i} = Sherlog16{T,i}
Base.promote_rule(::Type{Float16},::Type{Sherlog16{T,i}}) where {T,i} = Sherlog16{T,i}

Base.bitstring(x::Sherlog16) = bitstring(x.val)
Base.show(io::IO,x::Sherlog16) = print(io,"Sherlog16(",string(x.val),")")

Base.eps(::Type{Sherlog16{T,i}}) where {T,i} = eps(Float64)
Base.eps(x::Sherlog16) = eps(x.val)

Base.typemin(::Sherlog16{T,i}) where {T,i} = Sherlog16{T,i}(typemin(T))
Base.typemax(::Sherlog16{T,i}) where {T,i} = Sherlog16{T,i}(typemax(T))
Base.floatmin(::Sherlog16{T,i}) where {T,i} = Sherlog16{T,i}(floatmin(T))
Base.floatmax(::Sherlog16{T,i}) where {T,i} = Sherlog16{T,i}(floatmax(T))
Base.precision(::Sherlog16{T,i}) where {T,i} = Sherlog16{T,i}(precision(T))

function +(x::Sherlog16{T,i},y::Sherlog16{T,i}) where {T,i}
    r = x.val + y.val
    log_it(T,r,i)
    return Sherlog16{T,i}(r)
end

function -(x::Sherlog16{T,i},y::Sherlog16{T,i}) where {T,i}
    r = x.val - y.val
    log_it(T,r,i)
    return Sherlog16{T,i}(r)
end

function *(x::Sherlog16{T,i},y::Sherlog16{T,i}) where {T,i}
    r = x.val * y.val
    log_it(T,r,i)
    return Sherlog16{T,i}(r)
end

function /(x::Sherlog16{T,i},y::Sherlog16{T,i}) where {T,i}
    r = x.val / y.val
    log_it(T,r,i)
    return Sherlog16{T,i}(r)
end

function ^(x::Sherlog16{T,i},y::Sherlog16{T,i}) where {T,i}
    r = x.val ^ y.val
    log_it(T,r,i)
    return Sherlog16{T,i}(r)
end

for O in ( :(-), :(+),
           :sign,
           :prevfloat, :nextfloat,
           :round, :trunc, :ceil, :floor,
           :inv, :abs, :sqrt, :cbrt,
           :exp, :expm1, :exp2, :exp10,
           :log, :log1p, :log2, :log10,
           :rad2deg, :deg2rad, :mod2pi, :rem2pi,
           :sin, :cos, :tan, :csc, :sec, :cot,
           :asin, :acos, :atan, :acsc, :asec, :acot,
           :sinh, :cosh, :tanh, :csch, :sech, :coth,
           :asinh, :acosh, :atanh, :acsch, :asech, :acoth,
           :sinc, :sinpi, :cospi,
           :sind, :cosd, :tand, :cscd, :secd, :cotd,
           :asind, :acosd, :atand, :acscd, :asecd, :acotd
          )
    @eval begin
        function Base.$O(x::Sherlog16{T,i}) where {T,i}
            r = $O(x.val)
            log_it(T,r,i)
            return Sherlog16{T,i}(r)
        end
    end
end

for O in ( :(<), :(<=))
    @eval begin
        Base.$O(x::Sherlog16, y::Sherlog16) = $O(x.val, y.val)
    end
end


text:
using StatGeochemBase
using NaNStatistics
using Test

@testset "Math" begin include("testMath.jl") end
@testset "Images" begin include("testImages.jl") end
@testset "ArrayStats" begin include("testArrayStats.jl") end
@testset "Interpolations" begin include("testInterpolations.jl") end
@testset "Geochronology" begin include("testGeochronology.jl") end

max_stars_repo_name: brenhinkeller/StatGeochemBase

max_issues_repo_name: brenhinkeller/StatGeochemBase

max_forks_repo_name: brenhinkeller/StatGeochemBase


text:
#=
installIJulia:
- Julia version:
- Author: logan
- Date: 2019-09-13
=#

import Pkg
Pkg.add("IJulia")
Pkg.build("IJulia")

max_stars_repo_name: btrspg/jupyterlab-docker

max_issues_repo_name: btrspg/jupyterlab-docker

max_forks_repo_name: btrspg/jupyterlab-docker


text:
module RomuRandoms

export Romu, romu, state, setstate,
    RomuQuadState, RomuTrioState, RomuDuoState, RomuDuoJrState, RomuQuad32State, RomuTrio32State, RomuMono32State,
    romu_quad, romu_trio, romu_duo, romu_duojr, romu_quad32, romu_trio32, romu_mono32

struct RomuException <: Exception
    msg::AbstractString
end

include("primitives.jl")  # rotate bits (ror, rol)
include("romu_states.jl") # setstate, state structs
include("romu64s.jl")     # romu_quad, romu_trio, romu_duo, romu_duojr
include("romu32s.jl")     # romu_quad32, romu_trio32, romu_mono32

const Romu = RomuTrioState
const romu = romu_trio

end #  RomuRandoms


text:
using OmegaCore, Distributions, Test, OmegaTest, OmegaDistributions

a ⟹ b = !a || b

function test_hi()
  μ = 1 ~ Normal(0, 1)
  σ = 2 ~ Uniform(1, 3)
  y = 3 ~ Normal.(μ, σ)
  changeμ = 4 ~ Bernoulli(0.5)
  choice(ω) = ifelse((changeμ)(ω), μ, σ)
  int_dist(ω) = ValueIntervention(choice(ω), 5.0)

  # with a 40% chance, either μ or σ will be intervened to 0.5
  joint = @joint(y, μ, σ, changeμ)
  joint_ = joint |ᵈ hi(int_dist)
  y_, μ_, σ_, changeμ_ = randsample(joint_)
  if changeμ_
    @test μ_ == 5.0
  else
    @test σ_ == 5.0
  end
end

function test_hi_self()
  # Intervenes one of the variables to the negation of itself
  μ = 1 ~ Poisson(10)
  σ = 2 ~ Uniform(1, 3)
  y = Variable(ω -> (μ(ω), σ(ω)))
  choice(ω) = ifelse((4 ~ Bernoulli(0.5))(ω), μ, σ)
  int_dist(ω) = let x_ = choice(ω)
    Intervention(x_, 0.0 .- x_)
  end

  joint = @joint(y, μ, σ, choice)
  joint_ = joint |ᵈ hi(int_dist)
  y_, μ_, σ_, choice_ = randsample(joint_)
  @test (choice_ == μ) ⟹ μ_ <= 0.0
  @test (choice_ == σ) ⟹ σ_ <= 0.0
  # @test isinferred(randsample, joint_)
end

function test_cd()
  # A model where x --causes-> y
  function m1()
    x = 1 ~ Normal(0, 1)
    y(ω) = (2 ~ Normal(x(ω), 1))(ω)
    (x = x, y = y)
  end

  # A model where y --causes-> x
  function m2()
    x = 3 ~ Normal(0, 1)
    y(ω) = (4 ~ Normal(x(ω), 1))(ω)
    (x = y, y = x)
  end

  is_model_1 = 3 ~ Bernoulli(0.5)
  model = pw(ifelse, is_model_1, m1(), m2())
  mx(ω) = model(ω).x
  mxω(ω) = mx(ω)(ω)
  my(ω) = model(ω).y

  x_ = 1.3
  y_ = 2.0

  # This won't do what we want!
  condition1(ω) = (mxω |ᵈ (my => 100.0))(ω) >= 101.0

  condition(ω) = (mxω |ᵈ (my =>ˡ y_))(ω) >= x_

  # condition = (mxω |ᵈ (my =>ˡ y_)) >=ₚ x_

  prob = mean(randsample(is_model_1 |ᶜ condition, 1000))
  println("""The model1 is correct (x causes y), given after forcing y to 2 x,
  we observed that x was 1.3 is""", prob)
end

@testset "Higher Intervention" begin
  test_hi()
  test_hi_self()
  test_cd()
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# Important note: even if one or more tests fail, we will still exit with status code 0.
#
# The reason for this is that we always want to upload code coverage, even if some of the
# tests fail. Therefore, even if the `coverage_linux64` builder passes, you should not
# assume that all of the tests passed. If you want to know if all of the tests are passing,
# please look at the status of the `tester_*` builders (e.g. `tester_linux64`).

const ncores = Sys.CPU_THREADS
@info "" Sys.CPU_THREADS
@info "" ncores

script_native_yes = """
    Base.runtests(["cmdlineargs"]; ncores = $(ncores))
"""
script_native_no = """
    Base.runtests(["all", "--skip", "cmdlineargs"]; ncores = $(ncores))
"""

base_cmd       = `$(Base.julia_cmd()) --code-coverage=all`
cmd_native_yes = `$(base_cmd) --sysimage-native-code=yes -e $(script_native_yes)`
cmd_native_no  = `$(base_cmd) --sysimage-native-code=no  -e $(script_native_no)`

@info "Running command" cmd_native_yes
p1 = run(pipeline(cmd_native_yes; stdin, stdout, stderr); wait = false)
wait(p1)

@info "Running command" cmd_native_no
p2 = run(pipeline(cmd_native_no; stdin, stdout, stderr); wait = false)
wait(p2)

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
module LevelSets

## Included files
include("levels.jl")
include("grids.jl")
include("utils.jl")

export regulargrid2d
export regulargrid3d
export get_coords
export levelcurve
export compute_subfaces
export levelsurface
export load_mesh

end # module


text:
# test null hypothesis cases

using Test
using IncrementalInference
using TensorCast
# going to introduce two new constraint types

##

@testset "Post 237 without nullhypo" begin

##

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)

addFactor!(fg, [:x0;], Prior(Normal()))
addFactor!(fg, [:x0;:x1], LinearRelative(Normal(10,1)))

solveTree!(fg)

pts_ = getBelief(fg, :x1) |> getPoints
@cast pts[i,j] := pts_[j][i]

@test 80 < sum(2 .< pts .< 18)

##

end

@testset "Post 237 nullhypo on prior" begin

##

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(10,1)), nullhypo=0.5)

solveTree!(fg)

pts_ = getBelief(fg, :x0) |> getPoints
@cast pts[i,j] := pts_[j][i]

@test 10 < sum(-15 .< pts .< 4) < 60
@test 30 < sum(4 .< pts .< 16) < 85
# @test 10 < sum(16 .< pts .< 60) < 60

##

end

@testset "Post 237 with nullhypo test" begin

##

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)

addFactor!(fg, [:x0;], Prior(Normal()))
addFactor!(fg, [:x0;:x1], LinearRelative(Normal(10,1)), nullhypo=0.5)

pts_ = approxConv(fg, :x0x1f1, :x1)
@cast pts[i,j] := pts_[j][i]

@test 20 < sum(pts .< 5)
@test 20 < sum(5 .< pts .< 15)

solveTree!(fg)

pts2_ = getBelief(fg, :x1) |> getPoints
@cast pts2[i,j] := pts2_[j][i]

@test 30 < sum(8 .< pts2 .< 12) <= 80
@test 80 < sum(-10 .< pts2 .< 30)

##

end

#


text:
module ExactDiagonalization

using Arpack: eigs
using LinearAlgebra: Eigen
using Printf: @printf, @sprintf
using Base.Iterators: product
using SparseArrays: SparseMatrixCSC
using QuantumLattices: EnumerativeVectorSpace, VectorSpace, DulPermutations, Combinations
using QuantumLattices: OperatorUnit, OperatorProd, OperatorSum, Operator, Operators, MatrixRepresentation, Transformation
using QuantumLattices: AbstractLattice, Bonds, Hilbert, Fock, Spin, FockTerm, SpinTerm, Term, Metric, OIDToTuple, Table, Boundary, Generator, Image, Engine
using QuantumLattices: reparameter, id, idtype, expand, plain, creation, rank

import LinearAlgebra: eigen
import QuantumLattices: add!, ⊕, ⊗, parameternames, contentnames, getcontent, dtype, kind, update!, statistics, matrix, dimension, Parameters

export Sector, BinaryBasis, BinaryBasisRange, BinaryBases, TargetSpace
export EDMatrix, EDMatrixRepresentation, SectorFilter, EDKind, ED
export productable, sumable

"""
    abstract type Sector <: OperatorUnit

A sector of the Hilbert space which form the bases of an irreducible representation of the Hamiltonian of a quantum lattice system.
"""
abstract type Sector <: OperatorUnit end

"""
    TargetSpace{S<:Sector} <: EnumerativeVectorSpace{S}

The target Hilbert space in which the exact diagonalization method is performed, which could be the direct sum of several sectors.
"""
struct TargetSpace{S<:Sector} <: EnumerativeVectorSpace{S}
    sectors::Vector{S}
end
@inline contentnames(::Type{<:TargetSpace}) = (:table,)
@inline getcontent(target::TargetSpace, ::Val{:table}) = target.sectors
function add!(target::TargetSpace, sector::Sector)
    @assert all(map(sec->sumable(sec, sector), target.sectors)) "add! error: could not be direct summed."
    push!(target.sectors, sector)
    return target
end
function add!(target::TargetSpace, another::TargetSpace)
    for sector in another
        add!(target, sector)
    end
    return target
end

"""
    ⊕(sector::Sector, sectors::Union{Sector, TargetSpace}...) -> TargetSpace
    ⊕(target::TargetSpace, sectors::Union{Sector, TargetSpace}...) -> TargetSpace

Get the direct sum of sectors and target spaces.
"""
@inline function ⊕(sector::Sector, sectors::Union{Sector, TargetSpace}...)
    result = TargetSpace(sector)
    map(op->add!(result, op), sectors)
    return result
end
@inline function ⊕(target::TargetSpace, sectors::Union{Sector, TargetSpace}...)
    result = TargetSpace(copy(target.sectors))
    map(op->add!(result, op), sectors)
    return result
end

"""
    TargetSpace(sector::Sector, sectors::Sector...)

Construct a target space from sectors.
"""
@inline function TargetSpace(sector::Sector, sectors::Sector...)
    result = TargetSpace([sector])
    for sec in sectors
        add!(result, sec)
    end
    return result
end

# Binary bases commonly used in canonical fermionic and hardcore bosonic quantum lattice systems
"""
    BinaryBasis{I<:Unsigned}

Binary basis represented by an unsigned integer.
"""
struct BinaryBasis{I<:Unsigned}
    rep::I
    BinaryBasis(i::Integer) = (rep = Unsigned(i); new{typeof(rep)}(rep))
end
@inline Base.:(==)(basis₁::BinaryBasis, basis₂::BinaryBasis) = basis₁.rep == basis₂.rep
@inline Base.isequal(basis₁::BinaryBasis, basis₂::BinaryBasis) = isequal(basis₁.rep, basis₂.rep)
@inline Base.:<(basis₁::BinaryBasis, basis₂::BinaryBasis) = basis₁.rep < basis₂.rep
@inline Base.isless(basis₁::BinaryBasis, basis₂::BinaryBasis) = isless(basis₁.rep, basis₂.rep)
@inline Base.one(basis::BinaryBasis) = one(typeof(basis))
@inline Base.one(::Type{BinaryBasis{I}}) where {I<:Unsigned} = BinaryBasis(one(I))
@inline Base.zero(basis::BinaryBasis) = zero(typeof(basis))
@inline Base.zero(::Type{BinaryBasis{I}}) where {I<:Unsigned} = BinaryBasis(zero(I))
@inline Base.show(io::IO, basis::BinaryBasis) = @printf io "%s" string(basis.rep, base=2)
@inline Base.eltype(basis::BinaryBasis) = eltype(typeof(basis))
@inline Base.eltype(::Type{<:BinaryBasis}) = Int
@inline Base.IteratorSize(::Type{<:BinaryBasis}) = Base.SizeUnknown()

"""
    iterate(basis::BinaryBasis, state=nothing)

Iterate over the numbers of the occupied single-particle orbitals.
"""
function Base.iterate(basis::BinaryBasis, state=nothing)
    (pos, rep) = isnothing(state) ? (0, basis.rep) : (state[1], state[2])
    while rep>0
        pos += 1
        isodd(rep) && return (pos, (pos, rep÷2))
        rep ÷= 2
    end
    return nothing
end

"""
    one(basis::BinaryBasis, state::Integer) -> BinaryBasis

Get a new basis with the specified single-particle state occupied.
"""
@inline Base.one(basis::BinaryBasis, state::Integer) = BinaryBasis(basis.rep | one(basis.rep)<<(state-1))

"""
    isone(basis::BinaryBasis, state::Integer) -> Bool

Judge whether the specified single-particle state is occupied for a basis.
"""
@inline Base.isone(basis::BinaryBasis, state::Integer) = (basis.rep & one(basis.rep)<<(state-1))>0

"""
    zero(basis::BinaryBasis, state::Integer) -> BinaryBasis

Get a new basis with the specified single-particle state unoccupied.
"""
@inline Base.zero(basis::BinaryBasis, state::Integer) = BinaryBasis(basis.rep & ~(one(basis.rep)<<(state-1)))

"""
    iszero(basis::BinaryBasis, state::Integer) -> Bool

Judge whether the specified single-particle state is unoccupied for a basis.
"""
@inline Base.iszero(basis::BinaryBasis, state::Integer) = !isone(basis, state)

"""
    count(basis::BinaryBasis) -> Int
    count(basis::BinaryBasis, start::Integer, stop::Integer) -> Int

Count the number of occupied single-particle states.
"""
@inline Base.count(basis::BinaryBasis) = count(basis, 1, ndigits(basis.rep, base=2))
@inline function Base.count(basis::BinaryBasis, start::Integer, stop::Integer)
    result = 0
    for i = start:stop
        isone(basis, i) && (result += 1)
    end
    return result
end

"""
    ⊗(basis₁::BinaryBasis, basis₂::BinaryBasis) -> BinaryBasis

Get the direct product of two binary bases.
"""
@inline ⊗(basis₁::BinaryBasis, basis₂::BinaryBasis) = BinaryBasis(basis₁.rep|basis₂.rep)

"""
    BinaryBasis(states; filter=index->true)
    BinaryBasis{I}(states; filter=index->true) where {I<:Unsigned}

Construct a binary basis with the given occupied orbitals.
"""
@inline BinaryBasis(states; filter=index->true) = BinaryBasis{typeof(Unsigned(first(states)))}(states; filter=filter)
function BinaryBasis{I}(states; filter=index->true) where {I<:Unsigned}
    rep, eye = zero(I), one(I)
    for (index, state) in enumerate(states)
        filter(index) && (rep += eye<<(state-1))
    end
    return BinaryBasis(rep)
end

"""
    BinaryBasisRange{I<:Unsigned} <: VectorSpace{BinaryBasis{I}}

A continuous range of binary basis.
"""
struct BinaryBasisRange{I<:Unsigned} <: VectorSpace{BinaryBasis{I}}
    slice::UnitRange{I}
end
@inline Base.issorted(bbr::BinaryBasisRange) = true
@inline dimension(bbr::BinaryBasisRange) = length(bbr.slice)
@inline Base.getindex(bbr::BinaryBasisRange, i::Integer) = BinaryBasis(bbr.slice[i])

"""
    BinaryBases{B<:BinaryBasis, T<:AbstractVector{B}} <: Sector

A set of binary bases.
"""
struct BinaryBases{B<:BinaryBasis, T<:AbstractVector{B}} <: Sector
    id::Vector{Tuple{B, Float64}}
    table::T
end
@inline Base.issorted(bs::BinaryBases) = true
@inline Base.length(bs::BinaryBases) = dimension(bs)
@inline dimension(bs::BinaryBases) = length(bs.table)
@inline Base.:(==)(bs₁::BinaryBases, bs₂::BinaryBases) = isequal(bs₁.id, bs₂.id)
@inline Base.isequal(bs₁::BinaryBases, bs₂::BinaryBases) = isequal(bs₁.id, bs₂.id)
@inline Base.getindex(bs::BinaryBases, i::Integer) = bs.table[i]
@inline Base.eltype(bs::BinaryBases) = eltype(typeof(bs))
@inline Base.eltype(::Type{<:BinaryBases{B}}) where {B<:BinaryBasis} = B
@inline Base.iterate(bs::BinaryBases, state=1) = state>dimension(bs) ? nothing : (bs.table[state], state+1)
function Base.repr(bs::BinaryBases)
    result = String[]
    for (rep, nparticle) in bs.id
        if isnan(nparticle)
            push!(result, @sprintf "2^%s" count(rep))
        else
            push!(result, @sprintf "C(%s, %s)" count(rep) Int(nparticle))
        end
    end
    return join(result, " ⊗ ")
end
function Base.show(io::IO, bs::BinaryBases)
    @printf io "%s:\n" repr(bs)
    for i = 1:dimension(bs)
        @printf io "  %s\n" bs[i]
    end
end
@inline Base.findfirst(b::BinaryBasis, bs::BinaryBases) = searchsortedfirst(bs.table, b)
@inline Base.findfirst(b::BinaryBasis, bs::BinaryBases{<:BinaryBasis, <:BinaryBasisRange}) = Int(b.rep+1)

"""
    ⊗(bs₁::BinaryBases, bs₂::BinaryBases) -> BinaryBases

Get the direct product of two sets of binary bases.
"""
function ⊗(bs₁::BinaryBases, bs₂::BinaryBases)
    @assert productable(bs₁, bs₂) "⊗ error: the input two sets of bases cannot be direct producted."
    table = Vector{promote_type(eltype(bs₁), eltype(bs₂))}(undef, dimension(bs₁)*dimension(bs₂))
    count = 1
    for (b₁, b₂) in product(bs₁, bs₂)
        table[count] = b₁⊗b₂
        count += 1
    end
    return BinaryBases(sort!([bs₁.id; bs₂.id]; by=first), sort!(table))
end

"""
    productable(bs₁::BinaryBases, bs₂::BinaryBases) -> Bool

Judge whether two sets of binary bases could be direct producted.
"""
function productable(bs₁::BinaryBases, bs₂::BinaryBases)
    for (irr₁, irr₂) in product(bs₁.id, bs₂.id)
        isequal(irr₁[1].rep & irr₂[1].rep, 0) || return false
    end
    return true
end

"""
    sumable(bs₁::BinaryBases, bs₂::BinaryBases) -> Bool

Judge whether two sets of binary bases could be direct summed.

Strictly speaking, two sets of binary bases could be direct summed if and only if they have no intersection. The time complexity to check the intersection is O(n log n), which costs a lot when the dimension of the binary bases is huge. It is also possible to judge whether they could be direct summed by close investigations on their ids, i.e. the single-particle states and occupation number. It turns out that this is a multi-variable pure integer linear programming problem. In the future, this function would be implemented based on this observation. At present, the direct summability should be handled by the users in priori.
"""
@inline sumable(bs₁::BinaryBases, bs₂::BinaryBases) = true

"""
    BinaryBases(nstate::Integer)
    BinaryBases(states)

Construct a set of binary bases that does not preserve the particle number conservation.
"""
@inline BinaryBases(nstate::Integer) = BinaryBases([(BinaryBasis(1:nstate), NaN)], BinaryBasisRange(UInt(0):UInt(2^nstate-1)))
@inline BinaryBases(states) = BinaryBases(Tuple(states))
function BinaryBases(states::NTuple{N, Integer}) where N
    states = NTuple{N, eltype(states)}(sort!(collect(states); rev=true))
    com = DulPermutations{N}((false, true))
    table = Vector{BinaryBasis{typeof(Unsigned(first(states)))}}(undef, length(com))
    for (i, poses) in enumerate(com)
        table[i] = BinaryBasis(states; filter=index->poses[index])
    end
    return BinaryBases([(BinaryBasis(states), NaN)], table)
end

"""
    BinaryBases(nstate::Integer, nparticle::Integer)
    BinaryBases(states, nparticle::Integer)

Construct a set of binary bases that preserves the particle number conservation.
"""
@inline BinaryBases(nstate::Integer, nparticle::Integer) = BinaryBases(1:nstate, Val(nparticle))
@inline BinaryBases(states, nparticle::Integer) = BinaryBases(states, Val(nparticle))
function BinaryBases(states, ::Val{N}) where N
    com = Combinations{N}(sort!(collect(states); rev=true))
    table = Vector{BinaryBasis{typeof(Unsigned(first(states)))}}(undef, length(com))
    for (i, poses) in enumerate(com)
        table[end+1-i] = BinaryBasis(poses)
    end
    return BinaryBases([(BinaryBasis(states), Float64(N))], table)
end

# CSC-formed sparse matrix representation of an operator
"""
    matrix(op::Operator, braket::NTuple{2, BinaryBases}, table; dtype=valtype(op)) -> SparseMatrixCSC{dtype, Int}

Get the CSC-formed sparse matrix representation of an operator.

Here, `table` specifies the order of the operator ids.
"""
function matrix(op::Operator, braket::NTuple{2, BinaryBases}, table; dtype=valtype(op))
    bra, ket = braket[1], braket[2]
    ndata, intermediate = 1, zeros(ket|>eltype, rank(op)+1)
    data, indices, indptr = zeros(dtype, dimension(ket)), zeros(Int, dimension(ket)), zeros(Int, dimension(ket)+1)
    sequence = NTuple{rank(op), Int}(table[op[i]] for i in reverse(1:rank(op)))
    iscreation = NTuple{rank(op), Bool}(oid.index.iid.nambu==creation for oid in reverse(id(op)))
    for i = 1:dimension(ket)
        flag = true
        indptr[i] = ndata
        intermediate[1] = ket[i]
        for j = 1:rank(op)
            isone(intermediate[j], sequence[j])==iscreation[j] && (flag = false; break)
            intermediate[j+1] = iscreation[j] ? one(intermediate[j], sequence[j]) : zero(intermediate[j], sequence[j])
        end
        if flag
            nsign = 0
            statistics(eltype(op))==:f && for j = 1:rank(op)
                nsign += count(intermediate[j], 1, sequence[j]-1)
            end
            indices[ndata] = findfirst(intermediate[end], bra)
            data[ndata] = op.value*(-1)^nsign
            ndata += 1
        end
    end
    indptr[end] = ndata
    return SparseMatrixCSC(dimension(bra), dimension(ket), indptr, indices[1:ndata-1], data[1:ndata-1])
end

# Generic exact diagonalization method
"""
    EDMatrix{S<:Sector, M<:SparseMatrixCSC} <: OperatorProd{M, Tuple{S, S}}

Matrix representation of quantum operators between a ket and a bra Hilbert space.
"""
struct EDMatrix{S<:Sector, M<:SparseMatrixCSC} <: OperatorProd{M, Tuple{S, S}}
    bra::S
    ket::S
    matrix::M
end
@inline parameternames(::Type{<:EDMatrix}) = (:sector, :value)
@inline getcontent(m::EDMatrix, ::Val{:id}) = (m.bra, m.ket)
@inline getcontent(m::EDMatrix, ::Val{:value}) = m.matrix
@inline dtype(::Type{<:EDMatrix{<:Sector, M}}) where {M<:SparseMatrixCSC} = eltype(M)
@inline EDMatrix(m::SparseMatrixCSC, braket::NTuple{2, Sector}) = EDMatrix(braket[1], braket[2], m)
@inline Base.promote_rule(M::Type{<:EDMatrix}, N::Type{<:Number}) = reparameter(M, :value, reparameter(valtype(M), 1, promote_type(dtype(M), N)))

"""
    EDMatrix(sector::Sector, m::SparseMatrixCSC)

Construct a matrix representation when the ket and bra spaces share the same bases.
"""
@inline EDMatrix(sector::Sector, m::SparseMatrixCSC) = EDMatrix(sector, sector, m)

"""
    EDMatrixRepresentation{S<:Sector, T} <: MatrixRepresentation

Exact matrix representation of a quantum lattice system on a target Hilbert space.
"""
struct EDMatrixRepresentation{S<:Sector, T} <: MatrixRepresentation
    brakets::Vector{NTuple{2, S}}
    table::T
end
@inline function Base.valtype(::Type{<:EDMatrixRepresentation{S}}, M::Type{<:Operator}) where {S<:Sector}
    M = EDMatrix{S, SparseMatrixCSC{valtype(M), Int}}
    return OperatorSum{M, idtype(M)}
end
@inline Base.valtype(R::Type{<:EDMatrixRepresentation}, M::Type{<:Operators}) = valtype(R, eltype(M))
function (representation::EDMatrixRepresentation)(m::Operator; kwargs...)
    result = zero(valtype(representation, m))
    for braket in representation.brakets
        add!(result, EDMatrix(matrix(m, braket, representation.table; kwargs...), braket))
    end
    return result
end

"""
    EDMatrixRepresentation(target::TargetSpace, table)

Construct a exact matrix representation.
"""
@inline EDMatrixRepresentation(target::TargetSpace, table) = EDMatrixRepresentation([(sector, sector) for sector in target], table)

"""
    SectorFilter{S} <: Transformation

Filter the target bra and ket Hilbert spaces.
"""
struct SectorFilter{S} <: Transformation
    brakets::Set{NTuple{2, S}}
end
@inline Base.valtype(::Type{<:SectorFilter}, M::Type{<:OperatorSum{<:EDMatrix}}) = M
@inline (sectorfileter::SectorFilter)(m::EDMatrix) = id(m)∈sectorfileter.brakets ? m : 0
@inline SectorFilter(sector::Sector, sectors::Sector...) = SectorFilter((sector, sector), map(op->(op, op), sectors)...)
@inline SectorFilter(braket::NTuple{2, Sector}, brakets::NTuple{2, Sector}...) = SectorFilter(push!(Set{typeof(braket)}(), braket, brakets...))

"""
    EDKind{K}

The kind of the exact diagonalization method applied to a quantum lattice system.
"""
struct EDKind{K} end
@inline EDKind(K::Symbol) = EDKind{K}()
@inline EDKind(::Type{T}) where {T<:Term} = error("EDKind error: not defined for $(kind(T)).")
@inline EDKind(::Type{<:FockTerm}) = EDKind(:FED)
@inline EDKind(::Type{<:SpinTerm}) = EDKind(:SED)
@inline @generated function EDKind(::Type{TS}) where {TS<:Tuple{Vararg{Term}}}
    exprs = []
    for i = 1:fieldcount(TS)
        push!(exprs, :(typeof(EDKind(fieldtype(TS, $i)))))
    end
    return Expr(:call, Expr(:call, :reduce, :promote_type, Expr(:tuple, exprs...)))
end

"""
    Metric(::EDKind{:FED}, hilbert::Hilbert{<:Fock}) -> OIDToTuple
    Metric(::EDKind{:SED}, hilbert::Hilbert{<:Spin}) -> OIDToTuple

Get the oid-to-tuple metric for a free fermionic/bosonic system or a free phononic system.
"""
@inline @generated Metric(::EDKind{:FED}, hilbert::Hilbert{<:Fock}) = OIDToTuple(:spin, fieldnames(keytype(hilbert))..., :orbital)
@inline @generated Metric(::EDKind{:SED}, hilbert::Hilbert{<:Spin}) = OIDToTuple(fieldnames(keytype(hilbert))..., :orbital)

"""
    ED{L<:AbstractLattice, G<:Generator, M<:Image} <: Engine

Exact diagonalization method of a quantum lattice system.
"""
struct ED{K, L<:AbstractLattice, G<:Generator, M<:Image} <: Engine
    lattice::L
    H::G
    Hₘ::M
    function ED{K}(lattice::AbstractLattice, H::Generator, mr::EDMatrixRepresentation) where K
        @assert isa(K, EDKind) "ED error: wrong kind."
        Hₘ = mr(H)
        new{K, typeof(lattice), typeof(H), typeof(Hₘ)}(lattice, H, Hₘ)
    end
end
@inline kind(ed::ED) = kind(typeof(ed))
@inline kind(::Type{<:ED{K}}) where K = K
@inline Base.valtype(::Type{<:ED{K, <:AbstractLattice, G} where K}) where {G<:Generator} = valtype(eltype(G))
@inline statistics(ed::ED) = statistics(typeof(ed))
@inline statistics(::Type{<:ED{K, <:AbstractLattice, G} where K}) where {G<:Generator} = statistics(eltype(eltype(G)))
@inline function update!(ed::ED; kwargs...)
    if length(kwargs)>0
        update!(ed.H; kwargs...)
        update!(ed.Hₘ, ed.H)
    end
    return ed
end
@inline Parameters(ed::ED) = Parameters(ed.H)

"""
    ED(lattice::AbstractLattice, hilbert::Hilbert, terms::Tuple{Vararg{Term}}, targetspace::TargetSpace; boundary::Boundary=plain)

Construct the exact diagonalization method for a quantum lattice system.
"""
function ED(lattice::AbstractLattice, hilbert::Hilbert, terms::Tuple{Vararg{Term}}, targetspace::TargetSpace; boundary::Boundary=plain)
    K = EDKind(typeof(terms))
    H = Generator(terms, Bonds(lattice), hilbert; half=false, boundary=boundary)
    mr = EDMatrixRepresentation(targetspace, Table(hilbert, Metric(K, hilbert)))
    return ED{K}(lattice, H, mr)
end

"""
    matrix(ed::ED, sector::Sector; kwargs...) -> EDMatrix
    matrix(ed::ED, sector=first(ed.Hₘ.transformation.brakets); kwargs...) -> EDMatrix

Get the sparse matrix representation of a quantum lattice system in a sector of the target space.
"""
@inline matrix(ed::ED, sector::Sector; kwargs...) = matrix(ed, (sector, sector); kwargs...)
@inline matrix(ed::ED, braket::NTuple{2, Sector}=first(ed.Hₘ.transformation.brakets); kwargs...) = expand(SectorFilter(braket)(ed.Hₘ))[braket]

"""
    eigen(m::EDMatrix; nev=6, which=:SR, tol=0.0, maxiter=300, sigma=nothing, v₀=dtype(m)[]) -> Eigen

Solve the eigen problem by the restarted Lanczos method provided by the Arpack package.
"""
@inline function eigen(m::EDMatrix; nev=6, which=:SR, tol=0.0, maxiter=300, sigma=nothing, v₀=dtype(m)[])
    eigensystem = eigs(m.matrix; nev=nev, which=which, tol=tol, maxiter=maxiter, sigma=sigma, ritzvec=true, v0=v₀)
    return Eigen(eigensystem[1], eigensystem[2])
end

end # module


text:
"""
    X = sylvckr(A,B,C)

Solve the continuous Sylvester matrix equation

                AX + XB = C

using the Kronecker product expansion of equations. `A` and `B` are
square matrices, and `A` and `-B` must not have common eigenvalues.
This function is not recommended for large order matrices.
"""
function sylvckr(A,B,C)

   m, n = size(C);
   [m; n] == LinearAlgebra.checksquare(A,B) ||
             throw(DimensionMismatch("A, B and Q have incompatible dimensions"))
   reshape((kron(Array{eltype(A),2}(I,n,n),A)+
            kron(transpose(B),Array{eltype(B),2}(I,m,m)))\(C[:]),m,n)
end
"""
    X = sylvdkr(A,B,C)

Solve the discrete Sylvester matrix equation

                AXB + X = C

using the Kronecker product expansion of equations. `A` and `B` are
square matrices, and `A` and `-B` must not have common reciprocal eigenvalues.
This function is not recommended for large order matrices.
"""
function sylvdkr(A,B,C)

    m, n = size(C);
    [m; n] == LinearAlgebra.checksquare(A,B) ||
              throw(DimensionMismatch("A, B and C have incompatible dimensions"))
     reshape((kron(transpose(B),A)+I)\(C[:]),m,n)
end
"""
    X = gsylvkr(A,B,C,D,E)

Solve the generalized Sylvester matrix equation

                AXB + CXD = E

using the Kronecker product expansion of equations. `A`, `B`, `C` and `D` are
square matrices. The pencils `A-λC` and `D+λB` must be regular and
must not have common eigenvalues.
This function is not recommended for large order matrices.
"""
function gsylvkr(A,B,C,D,E)

    m, n = size(E);
    [m; n; m; n] == LinearAlgebra.checksquare(A,B,C,D) ||
                    throw(DimensionMismatch("A, B, C, D and E have incompatible dimensions"))
    reshape((kron(transpose(B),A)+kron(transpose(D),C))\(E[:]),m,n)
end
"""
    sylvsyskr(A,B,C,D,E,F) -> (X,Y)

Solve the Sylvester system of matrix equations

                AX + YB = C
                DX + YE = F

using the Kronecker product expansion of equations. `(A,D)`, `(B,E)` are
pairs of square matrices of the same size.
The pencils `A-λD` and `-B+λE` must be regular and must not have common eigenvalues.
This function is not recommended for large order matrices.
"""
function sylvsyskr(A,B,C,D,E,F)

    m, n = size(C);
    (m == size(F,1) && n == size(F,2)) ||
          throw(DimensionMismatch("C and F must have the same dimensions"))
    [m; n; m; n] == LinearAlgebra.checksquare(A,B,D,E) ||
                    throw(DimensionMismatch("A, B, C, D, E and F have incompatible dimensions"))
     z = [ kron(Array{eltype(A),2}(I,n,n),A) kron(transpose(B),Array{eltype(B),2}(I,m,m)) ;
          kron(Array{eltype(D),2}(I,n,n),D) kron(transpose(E),Array{eltype(E),2}(I,m,m))] \ [C[:];F[:]]
    (reshape(z[1:m*n],m,n),reshape(z[m*n+1:end],m,n))
end
"""
    dsylvsyskr(A,B,C,D,E,F) -> (X,Y)

Solve the dual Sylvester system of matrix equations

       AX + DY = C
       XB + YE = F

using the Kronecker product expansion of equations. `(A,D)`, `(B,E)` are
pairs of square matrices of the same size.
The pencils `A-λD` and `-B+λE` must be regular and must not have common eigenvalues.
This function is not recommended for large order matrices.
"""
function dsylvsyskr(A,B,C,D,E,F)

    m, n = size(C);
    (m == size(F,1) && n == size(F,2)) ||
          throw(DimensionMismatch("C and F must have the same dimensions"))
    [m; n; m; n] == LinearAlgebra.checksquare(A,B,D,E) ||
                    throw(DimensionMismatch("A, B, C, D, E and F have incompatible dimensions"))
    z = [ kron(Array{eltype(A),2}(I,n,n),A) kron(Array{eltype(D),2}(I,n,n),D);
          kron(transpose(B),Array{eltype(B),2}(I,m,m)) kron(transpose(E),Array{eltype(E),2}(I,m,m))] \ [C[:];F[:]]
    (reshape(z[1:m*n],m,n),reshape(z[m*n+1:end],m,n))
end


text:
using WAVI
function MISMIP_PLUS_Ice1r(initial_thickness)
    #Grid and boundary conditions
    nx = 80
    ny = 10
    nσ = 4
    x0 = 0.0
    y0 = -40000.0
    dx = 8000.0
    dy = 8000.0
    h_mask=trues(nx,ny)
    u_iszero = falses(nx+1,ny); u_iszero[1,:].=true
    v_iszero=falses(nx,ny+1); v_iszero[:,1].=true; v_iszero[:,end].=true
    grid = Grid(nx = nx,
                ny = ny,
                nσ = nσ,
                x0 = x0,
                y0 = y0,
                dx = dx,
                dy = dy,
                h_mask = h_mask,
                u_iszero = u_iszero,
                v_iszero = v_iszero)

    #Bed
    bed = WAVI.mismip_plus_bed #function definition

    #solver parameters
    maxiter_picard = 30
    solver_params = SolverParams(maxiter_picard = maxiter_picard)

    #Physical parameters
    accumulation_rate = 0.3
    params = Params(default_thickness = default_thickness,
                    accumulation_rate = accumulation_rate)

    #initial conditions
    initial_conditions = InitialConditions(initial_thickness = initial_thickness)

    #make the model
    model = Model(grid = grid,
                     bed_elevation = bed,
                     params = params,
                     solver_params = solver_params,
                     initial_conditions = initial_conditions)

    #embed the model with melt rate model
    function m1(h, b)
        draft = -(918.0 / 1028.0) * h
        cavity_thickness = draft .- b
        cavity_thickness = max.(cavity_thickness, 0)
        m =  0.2*tanh.(cavity_thickness./75).*max.((-100 .- draft), 0)
        return m
    end
    arguments = (h = model.fields.gh.h, b = model.fields.gh.b)
    melt_rate_model = AnalyticMeltRate(melt_partial_cell = true,
                                    melt_rate_function = m1,
                                    function_arguments = arguments)
    add_melt_rate_model!(model,melt_rate_model)

    #timestepping parameters
    dt = 0.1
    end_time = 100.
    timestepping_params = TimesteppingParams(dt = dt, end_time = end_time)

    simulation = Simulation(model = model,
                        timestepping_params = timestepping_params)

    #perform the simulation
    run_simulation!(simulation)
    return simulation
end

#@time simulation = MISMIP_PLUS_test();


text:
using Test
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

@testset "GLShader.jl" begin

   # function GLShader(vertex, fragment)
   @testset "GLShader" begin
      @test
      @test
      @test
      @test
   end

   # function releaseGpuResources(shader::GLShader)
   @testset "releaseGpuResources" begin
      @test
      @test
      @test
      @test
   end

   # function createShader(type,source)
   @testset "createShader" begin
      @test
      @test
      @test
      @test
   end

   # function enableProgram(shader)
   @testset "enableProgram" begin
      @test
      @test
      @test
      @test
   end

   # function disableProgram(shader)
   @testset "disableProgram" begin
      @test
      @test
      @test
      @test
   end

end


text:
## Code dealing with Face
mutable struct Face
    nvertex::Int64
    X::Array{Float64,1}
    Y::Array{Float64,1}
    Z::Array{Float64,1}
    I::Array{Int64,1}
    J::Array{Int64,1}
    K::Array{Int64,1}
end

"""Default Constructor for Face
"""
function Face()
    return Face(0,Float64[], Float64[],Float64[],Int64[],Int64[],Int64[])
end

""" Add a vertex to a face
    Args
        f: Face
        x: x-vertex
        y: y-vertex
        z: z-vertex
        i: index of x
        j: index of y
        k: index of z
"""
function add_face_vertex(f::Face, x::Float64, y::Float64, z::Float64,i::Int64,j::Int64,k::Int64)
    push!(f.X,x)
    push!(f.Y,y)
    push!(f.Z,z)
    push!(f.I,i)
    push!(f.J,j)
    push!(f.K,k)
    f.nvertex+=1
end

""" Check to see if two faces are the same by matching vertices
    Args
        f: another face
    Returns
        match_indices e.g. [[1,2],[1,4]] Face 1 matches Face 2, Face 1 matches Face 4
"""
function match_face_indicies(f1::Face,f2::Face)

    tol = 1E-6
    matchedIndices = Any[]
    for i in 1:f1.nvertex
       for j in 1:f1.nvertex
           dx = abs(f1.X[i] - f2.X[j])
           dy = abs(f1.Y[i] - f2.Y[j])
           dz = abs(f1.Z[i] - f2.Z[j])
           if ((dx<tol) && (dy<tol) && (dz<tol))
               push!(matchedIndices,[i j])
           end
        end
    end
    return matchedIndices
end

Base.:(==)(lhs::Face, rhs::Face) = (length(match_face_indicies(lhs,rhs)) == rhs.nvertex)
Base.:(!=)(lhs::Face, rhs::Face) = (length(match_face_indicies(lhs,rhs)) != lhs.nvertex)
## End Face


text:
module Day10

using AdventOfCode2019
using LinearAlgebra

function day10(input::String = readInput(joinpath(@__DIR__, "..", "data", "day10.txt")))
    asteroidMap = _generate_map(input)
    station, maxAsteroids = _evaluate_best_asteroid(asteroidMap)
    vapo = _vaporize(asteroidMap, station, 200)
    vapo_puzzle_coords = _to_puzzle_coordinates(vapo)
    solution_part2 = 100vapo_puzzle_coords[1] + vapo_puzzle_coords[2]
    return [maxAsteroids, solution_part2]
end

function _to_puzzle_coordinates(c::CartesianIndex)
    return CartesianIndex(c[2]-1, c[1]-1)
end

function _evaluate_best_asteroid(asteroidMap::Array{Bool,2})
    maxCoord = CartesianIndex(-1, -1)
    maxCount = 0
    asteroidCoords = findall(asteroidMap)
    for a in asteroidCoords
        count = 0
        cp = copy(asteroidCoords)
        filter!(x -> x != a, cp)
        while length(cp) != 0
            b = pop!(cp)
            list = Array{CartesianIndex,1}()
            for ast in cp
                if _on_line(ast, a, b)
                    push!(list, ast)
                end
            end
            filter!(x -> !(x in list), cp)
            if _in_between(a, b, list)
                count += 2
            else
                count += 1
            end
        end
        if count > maxCount
            maxCount = count
            maxCoord = a
        end
    end
    return maxCoord, maxCount
end

function _generate_map(input::String)
    sinp = split(input)
    nrows = length(sinp)
    ncols = length(sinp[1])
    asteroidMap = zeros(Bool, (nrows, ncols))
    for (i, row) in enumerate(sinp)
        for (j, col) in enumerate(row)
            if col == '#'
                asteroidMap[i, j] = 1
            end
        end
    end
    return asteroidMap
end

function _on_line(p::CartesianIndex, a::CartesianIndex, b::CartesianIndex)
    return dot(b.I .- a.I, p.I .- a.I)^2 == sum((b.I .- a.I).^2) * sum((p.I .- a.I).^2)
end

function _in_between(p::CartesianIndex, b::CartesianIndex, list::Array{CartesianIndex,1})
    m = min(p, b, list...)
    M = max(p, b, list...)
    if p[1] == m[1]
        (p[2] == m[2] || p[2] == M[2]) && return false
    elseif p[1] == M[1]
        (p[2] == m[2] || p[2] == M[2]) && return false
    end
    return true
end

function _vaporize(asteroidMap::Array{Bool,2}, station::CartesianIndex, n::Int)
    asteroidCoords = findall(asteroidMap)
    filter!(x->x!=station, asteroidCoords)

    asteroidCoordsRight = copy(asteroidCoords)
    filter!(x -> (x[2] > station[2] || (x[2] == station[2] && x[1] < station[1])), asteroidCoordsRight)
    asteroidCoordsLeft = copy(asteroidCoords)
    filter!(x->!(x in asteroidCoordsRight), asteroidCoordsLeft)

    dirsRight = Array{CartesianIndex,1}()
    while length(asteroidCoordsRight) != 0
        ast = pop!(asteroidCoordsRight)
        filter!(x->!_on_line(x, station, ast), asteroidCoordsRight)
        push!(dirsRight, ast)
    end

    up = CartesianIndex(-1, 0)
    sort!(dirsRight, by=x->dot(up.I, (x-station).I)/(norm(up.I)*norm((x-station).I)), rev=true)

    dirsLeft = Array{CartesianIndex,1}()
    while length(asteroidCoordsLeft) != 0
        ast = pop!(asteroidCoordsLeft)
        filter!(x->!_on_line(x, station, ast), asteroidCoordsLeft)
        push!(dirsLeft, ast)
    end

    down = CartesianIndex(1, 0)
    sort!(dirsLeft, by=x->dot(down.I, (x-station).I)/(norm(down.I)*norm((x-station).I)), rev=true)

    dirs = union(dirsRight, dirsLeft)

    i = 1
    c = 0
    filter!(x->x!=station, asteroidCoords)
    length(asteroidCoords) < n && throw(AssertionError("not enough asteroids"))
    while true
        dir = dirs[i]
        candidates = []
        for ast in asteroidCoords
            if _on_line(ast, station, dir) && dot((ast-station).I, (dir-station).I) > 0
                push!(candidates, ast)
            end
        end
        if length(candidates) > 0
            ind = argmin([norm((x-station).I) for x in candidates])
            vapo = candidates[ind]
            filter!(x->x!=vapo, asteroidCoords)
            c += 1
            if c == n
                return vapo
            end
        end
        i = mod1(i + 1, length(dirs))
    end
end

end # module


text:
import Base: *, ∈, isempty

export CartesianProduct,
       CartesianProduct!,
       swap

"""
    CartesianProduct{N, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N}

Type that represents the Cartesian product of two sets, that is the set

```math
Z = \\{ z ∈ \\mathbb{R}^{n + m} : z = (x, y),\\qquad x ∈ X, y ∈ Y \\}.
```
If ``X ⊆ \\mathbb{R}^n`` and ``Y ⊆ \\mathbb{R}^m``, then ``Z`` is ``n+m``-dimensional.

### Fields

- `X` -- first set
- `Y` -- second set

### Notes

The Cartesian product of three elements is obtained recursively.
See also [`CartesianProductArray`](@ref) for an implementation of a Cartesian product of
many sets without recursion, instead using an array, making the operations more efficient.

The `EmptySet` is the absorbing element for `CartesianProduct`.

The Cartesian product preserves convexity: if the set arguments are convex, then
their Cartesian product is convex as well.

In some docstrings the word "block" is used to denote each wrapped set, with the natural
order, i.e. we say that the first block of a Cartesian product `cp` is `cp.X` and
the second block is `cp.Y`.

### Examples

The Cartesian product between two sets `X` and `Y` can be constructed either
using `CartesianProduct(X, Y)` or the short-cut notation `X × Y` (to enter the times
symbol, write `\times[TAB]`).

```jldoctest cartesianproduct_constructor
julia> I1 = Interval(0, 1);

julia> I2 = Interval(2, 4);

julia> I12 = I1 × I2;

julia> typeof(I12)
CartesianProduct{Float64, Interval{Float64, IntervalArithmetic.Interval{Float64}}, Interval{Float64, IntervalArithmetic.Interval{Float64}}}
```
A hyperrectangle is the cartesian product of intervals, so we can convert `I12`
exactly to a `Hyperrectangle` type:

```jldoctest cartesianproduct_constructor
julia> convert(Hyperrectangle, I12)
Hyperrectangle{Float64, Vector{Float64}, Vector{Float64}}([0.5, 3.0], [0.5, 1.0])
```
"""
struct CartesianProduct{N, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N}
    X::S1
    Y::S2
    function CartesianProduct(X::LazySet{N}, Y::LazySet{N}) where {N}
        return new{N, typeof(X), typeof(Y)}(X, Y)
    end
end

isoperationtype(::Type{<:CartesianProduct}) = true
isconvextype(::Type{CartesianProduct{N, S1, S2}}) where {N, S1, S2} = isconvextype(S1) && isconvextype(S2)

# EmptySet is the absorbing element for CartesianProduct
@absorbing(CartesianProduct, EmptySet)

"""
```
    *(X::LazySet, Y::LazySet)
```

Alias for the binary Cartesian product.
"""
*(X::LazySet, Y::LazySet) = CartesianProduct(X, Y)

"""
    ×

Unicode alias constructor × (`times`) for the binary Cartesian product operator.

### Notes

Write `\\times[TAB]` to enter this symbol.
"""
×(X::LazySet, Y::LazySet) = CartesianProduct(X, Y)

"""
    swap(cp::CartesianProduct)

Return a new `CartesianProduct` object with the arguments swapped.

### Input

- `cp` -- Cartesian product of two sets

### Output

A new `CartesianProduct` object with the arguments swapped.
"""
function swap(cp::CartesianProduct)
    return CartesianProduct(cp.Y, cp.X)
end

"""
    dim(cp::CartesianProduct)

Return the dimension of a Cartesian product.

### Input

- `cp` -- Cartesian product

### Output

The ambient dimension of the Cartesian product.
"""
function dim(cp::CartesianProduct)
    return dim(cp.X) + dim(cp.Y)
end

"""
    σ(d::AbstractVector, cp::CartesianProduct)

Return the support vector of a Cartesian product.

### Input

- `d`  -- direction
- `cp` -- Cartesian product

### Output

The support vector in the given direction.
If the direction has norm zero, the result depends on the wrapped sets.
"""
function σ(d::AbstractVector, cp::CartesianProduct)
    n1 = dim(cp.X)
    return [σ(d[1:n1], cp.X); σ(d[n1+1:length(d)], cp.Y)]
end

"""
    ρ(d::AbstractVector, cp::CartesianProduct)

Return the support function of a Cartesian product.

### Input

- `d`  -- direction
- `cp` -- Cartesian product

### Output

The support function in the given direction.
If the direction has norm zero, the result depends on the wrapped sets.
"""
function ρ(d::AbstractVector, cp::CartesianProduct)
    n1 = dim(cp.X)
    return ρ(d[1:n1], cp.X) + ρ(d[n1+1:length(d)], cp.Y)
end

"""
    isbounded(cp::CartesianProduct)

Determine whether a Cartesian product is bounded.

### Input

- `cp` -- Cartesian product

### Output

`true` iff both wrapped sets are bounded.
"""
function isbounded(cp::CartesianProduct)
    return isbounded(cp.X) && isbounded(cp.Y)
end

"""
    ∈(x::AbstractVector, cp::CartesianProduct)

Check whether a given point is contained in a Cartesian product.

### Input

- `x`  -- point/vector
- `cp` -- Cartesian product

### Output

`true` iff ``x ∈ cp``.
"""
function ∈(x::AbstractVector, cp::CartesianProduct)
    @assert length(x) == dim(cp)

    n1 = dim(cp.X)
    return view(x, 1:n1) ∈ cp.X &&
           view(x, n1+1:length(x)) ∈ cp.Y
end

"""
    isempty(cp::CartesianProduct)

Return if a Cartesian product is empty or not.

### Input

- `cp` -- Cartesian product

### Output

`true` iff any of the sub-blocks is empty.
"""
function isempty(cp::CartesianProduct)
    return isempty(cp.X) || isempty(cp.Y)
end

"""
    constraints_list(cp::CartesianProduct)

Return the list of constraints of a (polyhedral) Cartesian product.

### Input

- `cp` -- Cartesian product

### Output

A list of constraints.
"""
function constraints_list(cp::CartesianProduct)
    return constraints_list(CartesianProductArray([cp.X, cp.Y]))
end

"""
    vertices_list(cp::CartesianProduct{N}) where {N}

Return the list of vertices of a (polytopic) Cartesian product.

### Input

- `cp` -- Cartesian product

### Output

A list of vertices.

### Algorithm

We assume that the underlying sets are polytopic.
Then the high-dimensional set of vertices is just the Cartesian product of the
low-dimensional sets of vertices.
"""
function vertices_list(cp::CartesianProduct{N}) where {N}
    # collect low-dimensional vertices lists
    vlist_low = (vertices_list(cp.X), vertices_list(cp.Y))

    # create high-dimensional vertices list
    vlist = Vector{Vector{N}}()
    m = length(vlist_low[1]) * length(vlist_low[2])
    sizehint!(vlist, m)
    for v1 in vlist_low[1]
        for v2 in vlist_low[2]
            push!(vlist, vcat(v1, v2))
        end
    end

    return vlist
end

"""
    linear_map(M::AbstractMatrix, cp::CartesianProduct)

Concrete linear map of a (polyhedral) Cartesian product.

### Input

- `M`  -- matrix
- `cp` -- Cartesian product of two sets

### Output

A polytope if `cp` is bounded and a polyhedron otherwise.

### Algorithm

We convert the Cartesian product to constraint representation and then call
`linear_map` on the corresponding polyhedron.

This is a fallback implementation and it will fail if the wrapped sets are not polyhedral.
"""
function linear_map(M::AbstractMatrix, cp::CartesianProduct)
    return linear_map_cartesian_product(M, cp)
end

function linear_map_cartesian_product(M, cp)
    @assert dim(cp) == size(M, 2) "a linear map of size $(size(M)) cannot " *
                                  "be applied to a set of dimension $(dim(cp))"

    # use constraint representation
    T = isbounded(cp) ? HPolytope : HPolyhedron
    P = T(constraints_list(cp))
    return linear_map(M, P)
end

function concretize(cp::CartesianProduct)
    return cartesian_product(concretize(cp.X), concretize(cp.Y))
end

function project(cp::CartesianProduct, block::AbstractVector{Int}; kwargs...)
    n1 = dim(cp.X)
    if block[end] <= n1
        # projection completely in the first block
        return project(cp.X, block; kwargs...)
    elseif block[1] > n1
        # projection completely in the second block
        return project(cp.Y, block .- n1; kwargs...)
    end
    # projection is a new Cartesian product of the block-wise projections
    for (i, bi) in enumerate(block)
        if bi > n1
            X = project(cp.X, block[1:i-1]; kwargs...)
            Y = project(cp.Y, block[i:end] .- n1; kwargs...)
            return CartesianProduct(X, Y)
        end
    end
end

"""
    project(cp::CartesianProduct{N, IT, HT}, block::AbstractVector{Int};
            [kwargs...]) where {N, IT<:Interval, HT<:AbstractHyperrectangle{N}}

Concrete projection of a Cartesian product between an interval and a hyperrectangle.

### Input

- `cp`       -- Cartesian product between an interval and a hyperrectangle
- `block`    -- block structure, a vector with the dimensions of interest

### Output

A hyperrectangle representing the projection of the cartesian product `cp` on the
dimensions specified by `block`.
"""
function project(cp::CartesianProduct{N, IT, HT}, block::AbstractVector{Int};
                 kwargs...) where {N, IT<:Interval, HT<:AbstractHyperrectangle{N}}
    I = cp.X
    H = cp.Y
    block_vec = collect(block)
    if 1 ∉ block_vec
        block_vec .-= 1
        cH = center(H)
        rH = radius_hyperrectangle(H)
    else
        cH = vcat(center(I), center(H))
        rH = vcat(radius_hyperrectangle(I), radius_hyperrectangle(H))
    end
    return Hyperrectangle(cH[block_vec], rH[block_vec], check_bounds=false)
end

"""
    project(cp::CartesianProduct{N, IT, ZT}, block::AbstractVector{Int};
            [kwargs...]) where {N, IT<:Interval, ZT<:AbstractZonotope{N}}

Concrete projection of the Cartesian product between an interval and a zonotopic set.

### Input

- `cp`       -- Cartesian product between an interval and a zonotopic set
- `block`    -- block structure, a vector with the dimensions of interest

### Output

A zonotope representing the projection of the cartesian product `cp` on the
dimensions specified by `block`.
"""
function project(cp::CartesianProduct{N, IT, ZT}, block::AbstractVector{Int};
                 kwargs...) where {N, IT<:Interval, ZT<:AbstractZonotope{N}}
    block_vec = collect(block)
    Z = cp.Y
    if 1 ∉ block_vec
        block_vec .-= 1
    else
        Z = convert(Zonotope, cp)
    end
    M = projection_matrix(block_vec, dim(Z), N)
    return linear_map(M, Z)
end

"""
    project(cp::CartesianProduct{N, IT, Union{VP1, VP2}},
            block::AbstractVector{Int};
            [kwargs...]) where {N, IT<:Interval, VP1<:VPolygon{N}, VP2<:VPolytope{N}}

Concrete projection of the Cartesian product between an interval and a set in vertex representation.

### Input

- `cp`       -- Cartesian product between an interval and a `VPolygon` or a `VPolytope`
- `block`    -- block structure, a vector with the dimensions of interest

### Output

A `VPolytope` representing the projection of the cartesian product `cp` on the
dimensions specified by `block`.
"""
function project(cp::CartesianProduct{N, IT, Union{VP1, VP2}},
                 block::AbstractVector{Int};
                 kwargs...) where {N, IT<:Interval, VP1<:VPolygon{N}, VP2<:VPolytope{N}}
    I = cp.X
    P = cp.Y
    block_vec = collect(block)
    if 1 ∉ block_vec
        Pout = project(P, block_vec .- 1; kwargs...)
    else
        out = cartesian_product(I, P)
        Pout = project(out, block_vec; kwargs...)
    end
    return Pout
end

function volume(cp::CartesianProduct)
    return volume(cp.X) * volume(cp.Y)
end


text:
using SimJulia
using Dates

@resumable function datetimetest(sim::Simulation)
  println(nowDatetime(sim))
  @yield timeout(sim, Day(2))
  println(nowDatetime(sim))
end

datetime = now()
sim = Simulation(datetime)
@process datetimetest(sim)
run(sim, datetime+Month(3))
println(nowDatetime(sim))


text:
struct InputListener
    pipeline::Vector{Channel}
end

function InputListener(; size::Real=Inf)
    sequence_queue = Channel{String}(size, spawn=true) do ch
        while true
            sequence = T.read_stream()
            put!(ch, sequence)
        end
    end
    event_queue = Channel{T.Event}(size, spawn=true) do ch
        while true
            sequence = take!(sequence_queue)
            put!(ch, T.parse_sequence(sequence))
        end
    end

    return InputListener([sequence_queue, event_queue])
end

struct QuitEvent <: T.Event end
Base.show(io::IO, ::QuitEvent) = Base.print(io, "QuitEvent")

struct App
    size::Tuple{Int, Int}
    listener::InputListener
    model::Model
    view::View
end

event_queue(app::App) = app.listener.pipeline[end]

function init_term(app::App)
    T.cshow(false)
    Sys.iswindows() ? T.clear() : T.alt_screen(true)
    h, w = app.size; T.displaysize(h+3, w)
    T.raw!(true)
end

function reset_term(::App)
    T.raw!(false)
    sleep(1)
    !(Sys.iswindows()) && T.alt_screen(false)
    T.cshow()

    return
end

emit_quit_event(app::App) = put!(event_queue(app), QuitEvent())

function handle_lose(app::App)
    T.cmove_line_last()
    paint(app.view, state=:lose)
    emit_quit_event(app)
end

function handle_quit(app::App)
    keep_running = false
    foreach(close, app.listener.pipeline)

    T.cmove_line_last()
    T.cmove_up(2)
    T.println("Shutted down...")

    return keep_running
end

function handle_event(app::App)
    snake = app.model
    auto_move(snake)

    is_running = true
    while is_running
        e = take!(event_queue(app))
        if e === QuitEvent()
            is_running = handle_quit(app)
        elseif e === LossEvent()
            handle_lose(app)
        elseif e === UpdateEvent()
            paint(app.view)
        elseif e === T.KeyPressedEvent(T.UP)
            move(snake, :up)
        elseif e === T.KeyPressedEvent(T.DOWN)
            move(snake, :down)
        elseif e === T.KeyPressedEvent(T.RIGHT)
            move(snake, :right)
        elseif e === T.KeyPressedEvent(T.LEFT)
            move(snake, :left)
        elseif e === T.KeyPressedEvent(T.ESC)
            emit_quit_event(app)
        end
    end
end

function Base.run(app::App)
    init_term(app)
    paint(app.view)
    handle_event(app)
    reset_term(app)
end


text:
@testset "macro" begin
    @testset "private name" begin
        @system SPrivateName(Controller) begin
            _a: _aa => 1 ~ preserve
            __b: __bb => 2 ~ preserve
        end
        s = instance(SPrivateName)
        @test_throws ErrorException s._a
        @test s.__SPrivateName__a' == 1
        @test s.__SPrivateName__aa === s.__SPrivateName__a
        @test s.__b' == 2
        @test s.__bb === s.__b
    end

    @testset "private name args" begin
        @system SPrivateNameArgs(Controller) begin
            _x: _xx => 1 ~ preserve
            a(_x) ~ preserve
            b(_x) => x ~ preserve
            c(y=_x) => y ~ preserve
            d(_y=_x) => _y ~ preserve
            aa(_xx) ~ preserve
            bb(_xx) => xx ~ preserve
            cc(y=_xx) => y ~ preserve
            dd(_y=_xx) => _y ~ preserve
        end
        s = instance(SPrivateNameArgs)
        @test s.a' == 1
        @test s.b' == 1
        @test s.c' == 1
        @test s.d' == 1
        @test s.aa' == 1
        @test s.bb' == 1
        @test s.cc' == 1
        @test s.dd' == 1
    end

    @testset "private name tags" begin
        @system SPrivateNameTags(Controller) begin
            _t: _true => true ~ flag
            _f: _false => false ~ flag
            a => 1 ~ track(when=_t)
            b => 1 ~ track(when=_t|_f)
            c => 1 ~ track(when=_t&_f)
            d => 1 ~ track(when=_false)
        end
        s = instance(SPrivateNameTags)
        @test s.a' == 1
        @test s.b' == 1
        @test s.c' == 0
        @test s.d' == 0
    end

    @testset "private name mixin" begin
        @system SPrivateNameMixin1 begin
            _a: aa1 => 1 ~ preserve
            b(_a) ~ track
        end
        @system SPrivateNameMixin2 begin
            _a: aa2 => 2 ~ preserve
            c(_a) ~ track
        end
        @eval @system SPrivateNameMixed1(SPrivateNameMixin1, SPrivateNameMixin2, Controller)
        s1 = instance(SPrivateNameMixed1)
        @test_throws ErrorException s1._a
        @test s1.__SPrivateNameMixin1__a' == 1
        @test s1.__SPrivateNameMixin2__a' == 2
        @test s1.b' == 1
        @test s1.c' == 2
        @eval @system SPrivateNameMixed11(SPrivateNameMixed1, Controller)
        s11 = instance(SPrivateNameMixed11)
        @test_throws ErrorException s11._a
        @test s11.__SPrivateNameMixin1__a' == 1
        @test s11.__SPrivateNameMixin2__a' == 2
        @test s11.b' == 1
        @test s11.c' == 2
        @test_throws ErrorException @eval @system SPrivateNameMixed2(SPrivateNameMixin1, Controller) begin
            _a: aa1 => 3 ~ preserve
        end
        @eval @system SPrivateNameMixed3(SPrivateNameMixin1, Controller) begin
            d(_a) ~ track
        end
        @test_throws UndefVarError instance(SPrivateNameMixed3)
        @eval @system SPrivateNameMixed4(SPrivateNameMixin1, SPrivateNameMixin2, Controller) begin
            d(a=__SPrivateNameMixin1__a) ~ track
            e(a=__SPrivateNameMixin2__a) ~ track
        end
        s4 = instance(SPrivateNameMixed4)
        @test s4.d' == 1
        @test s4.e' == 2
    end

    @testset "alias" begin
        @system SAlias(Controller) begin
            a: aa => 1 ~ track
            b(a, aa) => a + aa ~ track
        end
        s = instance(SAlias)
        @test s.a' == s.aa' == 1
        @test s.b' == 2
    end

    @testset "single arg without key" begin
        @system SSingleArgWithoutKey(Controller) begin
            a => 1 ~ track
            b(a) ~ track
            c(x=a) ~ track
        end
        s = instance(SSingleArgWithoutKey)
        @test s.a' == 1
        @test s.b' == 1
        @test s.c' == 1
    end

    @testset "bool arg" begin
        @system SBoolArg(Controller) begin
            t => true ~ preserve::Bool
            f => false ~ preserve::Bool
            a(t & f) ~ track::Bool
            b(t | f) ~ track::Bool
            c(t & !f) ~ track::Bool
            d(x=t&f, y=t|f) => x | y ~ track::Bool
        end
        s = instance(SBoolArg)
        @test s.a' == false
        @test s.b' == true
        @test s.c' == true
        @test s.d' == true
    end

    @testset "type alias" begin
        @system STypeAlias(Controller) begin
            a => -1 ~ preserve::int
            b => 1 ~ preserve::uint
            c => 1 ~ preserve::float
            d => true ~ preserve::bool
            e => :a ~ preserve::sym
            f => "A" ~ preserve::str
            g => nothing ~ preserve::∅
            h => missing ~ preserve::_
            i => Cropbox.Dates.Date(2021) ~ preserve::date
            j => Cropbox.TimeZones.ZonedDateTime(2021, Cropbox.TimeZones.tz"UTC") ~ preserve::datetime
        end
        s = instance(STypeAlias)
        @test s.a' isa Int64 && s.a' === -1
        @test s.b' isa UInt64 && s.b' == 1
        @test s.c' isa Float64 && s.c' === 1.0
        @test s.d' isa Bool && s.d' === true
        @test s.e' isa Symbol && s.e' === :a
        @test s.f' isa String && s.f' === "A"
        @test s.g' isa Nothing && s.g' === nothing
        @test s.h' isa Missing && s.h' === missing
        @test s.i' isa Cropbox.Dates.Date && s.i' === Cropbox.Dates.Date(2021)
        @test s.j' isa Cropbox.TimeZones.ZonedDateTime && s.j' === Cropbox.TimeZones.ZonedDateTime(2021, Cropbox.tz"UTC")
    end

    @testset "type union nothing" begin
        @system STypeUnionNothing(Controller) begin
            a ~ preserve::{sym|∅}(parameter)
        end
        a1 = :hello
        s1 = instance(STypeUnionNothing; config=:0 => :a => a1)
        @test s1.a' isa Union{Symbol,Nothing}
        @test s1.a' === a1
        a2 = nothing
        s2 = instance(STypeUnionNothing; config=:0 => :a => a2)
        @test s2.a' isa Union{Symbol,Nothing}
        @test s2.a' === a2
    end

    @testset "type union missing" begin
        @system STypeUnionMissing(Controller) begin
            a ~ preserve::{int|_}(parameter)
        end
        a1 = 0
        s1 = instance(STypeUnionMissing; config=:0 => :a => a1)
        @test s1.a' isa Union{Int64,Missing}
        @test s1.a' === a1
        a2 = missing
        # missing filtered out by genparameter, leaving no return value (nothing)
        @test_throws MethodError instance(STypeUnionMissing; config=:0 => :a => a2)
    end

    @testset "type vector single" begin
        @system STypeVectorSingle(Controller) begin
            a ~ preserve::int[](parameter)
        end
        a = [1, 2, 3]
        s = instance(STypeVectorSingle; config=:0 => :a => a)
        @test s.a' isa Vector{Int64}
        @test s.a' === a
    end

    @testset "type vector union" begin
        @system STypeVectorUnion(Controller) begin
            a ~ preserve::{sym|∅}[](parameter)
        end
        a = [:a, nothing, :c]
        s = instance(STypeVectorUnion; config=:0 => :a => a)
        @test s.a' isa Vector{Union{Symbol,Nothing}}
        @test s.a' === a
    end

    @testset "patch type" begin
        @system SPatchType{A => Int, Int => Float32, int => Float64}(Controller) begin
            a => 1 ~ preserve::A
            b => 1 ~ preserve::Int
            c => 1 ~ preserve::int
        end
        s = instance(SPatchType)
        @test s.a' isa Int
        @test s.b' isa Float32
        @test s.c' isa Float64
    end

    @testset "patch type cascade" begin
        @system SPatchTypeCascade0{A => Int} begin
            a => 1 ~ preserve::A
        end
        @eval @system SPatchTypeCascade(SPatchTypeCascade0, Controller) begin
            b => 1 ~ preserve::A
        end
        s = instance(SPatchTypeCascade)
        @test s.a' isa Int
        @test s.b' isa Int
    end

    @testset "patch const" begin
        @system SPatchConst{x = 1, y = u"m"}(Controller) begin
            a => 1 ~ accumulate(init=x, unit=y)
        end
        s = instance(SPatchConst)
        @test Cropbox.unittype(s.a) === u"m"
        @test s.a' == 1u"m"
        update!(s)
        @test s.a' == 2u"m"
    end

    @testset "dynamic type" begin
        @system SDynamicTypeBase(Controller)
        @system SDynamicTypeChild(Controller) <: SDynamicTypeBase
        @test SDynamicTypeChild <: SDynamicTypeBase
        @system SDynamicTypeBaseDynamic(Controller) begin
            x ~ <:SDynamicTypeBase(override)
        end
        @system SDynamicTypeBaseStatic(Controller) begin
            x ~ ::SDynamicTypeBase(override)
        end
        b = instance(SDynamicTypeBase)
        c = instance(SDynamicTypeChild)
        @test instance(SDynamicTypeBaseDynamic; options=(; x = b)).x === b
        @test instance(SDynamicTypeBaseDynamic; options=(; x = c)).x === c
        @test instance(SDynamicTypeBaseStatic; options=(; x = b)).x === b
        @test_throws MethodError instance(SDynamicTypeBaseStatic; options=(; x = c))
    end

    @testset "duplicate variable" begin
        @test_logs (:warn, "duplicate variable") @eval @system SDuplicateVariable begin
            x => 1 ~ preserve
            x => 2 ~ preserve
        end
    end

    @testset "body replacement" begin
        @system SBodyReplacement1(Controller) begin
            a => 1 ~ preserve
        end
        @eval @system SBodyReplacement2(SBodyReplacement1) begin
            a => 2
        end
        s1 = instance(SBodyReplacement1)
        s2 = instance(SBodyReplacement2)
        @test s1.a isa Cropbox.Preserve
        @test s1.a' == 1
        @test s2.a isa Cropbox.Preserve
        @test s2.a' == 2
    end

    @testset "replacement with different alias" begin
        @system SReplacementDifferentAlias1 begin
            x: aaa => 1 ~ preserve
        end
        @test_logs (:warn, "variable replaced with inconsistent alias") @eval @system SReplacementDifferentAlias2(SReplacementDifferentAlias1) begin
            x: bbb => 2 ~ preserve
        end
    end

    @testset "custom system" begin
        abstract type SAbstractCustomSystem <: System end
        @system SCustomSystem <: SAbstractCustomSystem
        @test SCustomSystem <: System
        @test SCustomSystem <: SAbstractCustomSystem
    end

    @testset "return" begin
        #TODO: reimplement more robust return checking
        @test_skip @test_throws LoadError @eval @system SReturn begin
            x => begin
                return 1
            end ~ track
        end
    end
end


text:
@testset "OpenCL.Event" begin
    @testset "OpenCL.Event status" begin
        for platform in cl.platforms()
            if contains(platform[:name], "Portable")
                msg = "$(platform[:name]) does not implement User Events"
                warn(msg)
                continue
            end

            for device in cl.devices(platform)
                ctx = cl.Context(device)
                evt = cl.UserEvent(ctx)
                evt[:status]
                @test evt[:status] == :submitted
                cl.complete(evt)
                @test evt[:status] == :complete
                finalize(evt)
            end
        end
    end

    @testset "OpenCL.Event wait" begin
        for platform in cl.platforms()
            if contains(platform[:name], "Portable") ||
               contains(platform[:name], "Intel Gen OCL")
                msg = "$(platform[:name]) does not implement User Events or shows other problems"
                warn(msg)
                continue
            end

            for device in cl.devices(platform)
                ctx = cl.Context(device)
                # create user event
                usr_evt = cl.UserEvent(ctx)
                q = cl.CmdQueue(ctx)
                cl.enqueue_wait_for_events(q, usr_evt)

                # create marker event
                mkr_evt = cl.enqueue_marker(q)

                @test usr_evt[:status] == :submitted
                @test mkr_evt[:status] in (:queued, :submitted)

                cl.complete(usr_evt)
                @test usr_evt[:status] == :complete

                cl.wait(mkr_evt)
                @test mkr_evt[:status] == :complete

                @test cl.cl_event_status(:running) == cl.CL_RUNNING
                @test cl.cl_event_status(:submitted) == cl.CL_SUBMITTED
                @test cl.cl_event_status(:queued) == cl.CL_QUEUED
                @test cl.cl_event_status(:complete) == cl.CL_COMPLETE
            end
        end
    end

    @testset "OpenCL.Event callback" begin
        for platform in cl.platforms()
            v = cl.opencl_version(platform)
            if v.major == 1 && v.minor < 1
                info("Skipping OpenCL.Event callback for $(platform[:name]) version < 1.1")
                continue
            end

            if contains(platform[:name], "Portable") ||
               contains(platform[:name], "Intel Gen OCL")
                msg = "$(platform[:name]) does not implement User Events or shows other problems."
                warn(msg)
                continue
            end

            for device in cl.devices(platform)
                callback_called = false

                function test_callback(evt, status)
                    callback_called = true
                    println("Test Callback")
                end

                ctx = cl.Context(device)
                usr_evt = cl.UserEvent(ctx)
                queue = cl.CmdQueue(ctx)

                cl.enqueue_wait_for_events(queue, usr_evt)

                mkr_evt = cl.enqueue_marker(queue)
                cl.add_callback(mkr_evt, test_callback)

                @test usr_evt[:status] == :submitted
                @test mkr_evt[:status] in (:queued, :submitted)
                @test callback_called == false

                cl.complete(usr_evt)
                @test usr_evt[:status] == :complete

                cl.wait(mkr_evt)

                # Give callback some time to finish
                yield()
                sleep(0.5)

                @test mkr_evt[:status] == :complete
                @test callback_called
            end
        end
    end
end


text:
import Pkg
Pkg.add("Pluto")
Pkg.add("PlutoUI")
Pkg.add("Plots")
Pkg.add("PlotlyJS")
Pkg.precompile()
Pkg.gc()


text:
abstract type RegressionLikelihood{T<:Real} <: Likelihood{T} end

include("gaussian.jl")
include("studentt.jl")
include("laplace.jl")
include("heteroscedastic.jl")

""" Return the labels in a vector of vectors for multiple outputs"""
function treat_labels!(y::AbstractArray{T,N},likelihood::L) where {T,N,L<:RegressionLikelihood}
    @assert T<:Real "For regression target(s) should be real valued"
    @assert N <= 2 "Target should be a matrix or a vector"
    if N == 1
        return [y],1,likelihood
    else
        return [y[:,i] for i in 1:size(y,2)],size(y,2),likelihood
    end
end


text:
export find_point_group
import Molecules: center_of_mass, translate

function find_point_group(mol::Molecule)
    # Center Molecule on COM
    mol = translate(mol, center_of_mass(mol))
    # Calculate moment of inertia tensor (moit) eigenvalues
    Ia_mol, Ib_mol, Ic_mol = eigenmoit(calcmoit(mol))[1]

    if Ia_mol == 0.0
        # Linear Molecule
        if Molecules.isequivalent(Molecules.transform(mol, Molecules.inversion_matrix()), mol)
            return("Dinfh")
        else
            return("Cinfv")
        end
    elseif isapprox(Ia_mol, Ib_mol, atol=tol) && isapprox(Ia_mol, Ic_mol, atol = tol)
        # Molecule with high symmetry
        D = buildD(mol)
        SEAs = findSEA(D, 5)
        n = num_C2(mol, SEAs)
        invertable = Molecules.isequivalent(Molecules.transform(mol, Molecules.inversion_matrix()), mol)
        if n == 15
            return("Ih")
        elseif n == 9
            if invertable
                return("Oh")
            else
                return("O")
            end
        elseif n == 3
            if invertable
                return("Th")
            else
                return("Td")
            end
        end
    else
        # Build distance matrix and find symmetry equivalent atoms
        D = buildD(mol)
        SEAs = findSEA(D, 5)
        # Find sets of rotation elements
        rot_set = find_rotation_sets(mol, SEAs)
        if size(rot_set)[1] < 1
            # Path followed if no rotation elements found by 'find_rotation_sets'
            # Therefore, maximum rotation order of 2
            c2 = find_c2(mol, SEAs)
            if c2 !== nothing
                c2_ortho = is_there_ortho_c2(mol, c2, SEAs)
                σh = is_there_sigmah(mol, c2)
                if c2_ortho
                    if σh
                        return("D2h")
                    else
                        σv = is_there_sigmav(mol, SEAs, c2)
                        if σv
                            return("D2d")
                        else
                            return("D2")
                        end
                    end
                else
                    if σh
                        return("C2h")
                    else
                        σv = is_there_sigmav(mol, SEAs, c2)
                        if σv
                            return("C2v")
                        else
                            S4 = Molecules.Sn(c2, 4)
                            molB = Molecules.transform(mol, S4)
                            check = Molecules.isequivalent(mol, molB)
                            if check
                                return("S4")
                            else
                                return("C2")
                            end
                        end
                    end
                end
            else
                if Molecules.isequivalent(Molecules.transform(mol, Molecules.inversion_matrix()), mol)
                    return("Ci")
                else
                    σv = is_there_sigmav(mol, SEAs, [0.0,0.0,0.0])
                    if σv
                        return("Cs")
                    else
                        return("C1")
                    end
                end
            end
        else
            # Path followed if rotation elements found by 'find_rotation_sets'
            rots = find_rotations(mol, rot_set)
            Cn = highest_ordered_axis(rots)
            paxis = rots[1].axis
            q1 = is_there_ortho_c2(mol, paxis, SEAs)
            q2 = is_there_sigmah(mol, paxis)
            q3 = is_there_sigmav(mol, SEAs, paxis)
            if q1
                if q2
                    return("D"*string(Cn)*"h")
                else
                    if q3
                        return("D"*string(Cn)*"d")
                    else
                        return("D"*string(Cn))
                    end
                end
            else
                if q2
                    return("C"*string(Cn)*"h")
                else
                    if q3
                        return("C"*string(Cn)*"v")
                        return("C",Cn,"v")
                    else
                        S2n = Molecules.Sn(paxis, Cn*2)
                        molB = Molecules.transform(mol, S2n)
                        check = Molecules.isequivalent(mol, molB)
                        if check
                        return("S"*string(Cn*2))
                            return("S",Cn*2)
                        else
                        return("C"*string(Cn))
                            return("C",Cn)
                        end
                    end
                end
            end
        end
    end
end


text:
using FlightSims
import FlightSims: State, Dynamics!
using DataFrames
using ComponentArrays
using UnPack
using Transducers
using Plots
using SciMLBase
using Test

struct LinearSystem_ZOH_Input <: AbstractEnv
    linear_env::LinearSystem
end

function State(env::LinearSystem_ZOH_Input)
    State(env.linear_env)
end

function Dynamics!(env::LinearSystem_ZOH_Input)
    @Loggable function dynamics!(dx, x, input, t)
        @nested_log Dynamics!(env.linear_env)(dx, x, nothing, t; u=input)
    end
end

function main()
    A = [0 1;
         0 0]
    B = [0 1]'
    linear_env = LinearSystem(A, B)
    env = LinearSystem_ZOH_Input(linear_env)
    x0_state = [1, 2]
    input = zeros(1)
    x0 = State(env)(x0_state)
    t0 = 0.0
    tf = 20.0
    simulator = Simulator(x0, Dynamics!(env), input; tf=tf)
    # interactive sim
    Δt = 0.1  # save period
    df = DataFrame()
    @time for (i, t) in enumerate(t0:Δt:tf)
        # To perform interactive simulation,
        # you should be aware of the integrator interface
        # provided by DifferentialEquations.jl;
        # see https://diffeq.sciml.ai/stable/basics/integrator/#integrator
        # e.g., DO NOT directly change the integrator state
        state = simulator.integrator.u
        input = simulator.integrator.p
        if (i-1) % 10 == 0  # update input period
            input .= -sum(state)
        end
        step_until!(simulator, t, df)
    end
    # plot
    ts = df.time
    states = df.sol |> Map(datum -> datum.state) |> collect
    inputs = df.sol |> Map(datum -> datum.input) |> collect
    p_x = plot(ts, hcat(states...)';
               title="state variable", label=["x1" "x2"],
              )
    p_u = plot(ts, hcat(inputs...)';
               title="input variable", label="u",
               linetype=:steppost,  # to plot zero-order-hold input appropriately
              )
    fig = plot(p_x, p_u; layout=(2, 1))
    savefig("figures/interactive_sim.png")
    display(fig)
end

@testset "linear_system_zoh_input" begin
    main()
end