text:
@testset "Approximate" begin
    seq = dna"ACGTACG"

    @testset "forward" begin
        @test approxsearch(seq, dna"", 0) === 1:0
        @test approxsearch(seq, dna"AC", 0) === 1:2
        @test approxsearch(seq, dna"AC", 0, 2) === 5:6
        @test approxsearch(seq, dna"AC", 0, 2, 5) === 0:-1
        @test approxsearch(seq, dna"AT", 0) === 0:-1
        @test approxsearch(seq, dna"AT", 1) === 1:1
        @test approxsearch(seq, dna"AT", 1, 2) === 3:4
        @test approxsearch(seq, dna"AT", 1, 2, 3) === 0:-1
        @test approxsearch(seq, dna"NG", 0) === 2:3
        @test approxsearch(seq, dna"NG", 1) === 1:1
        @test approxsearch(seq, dna"GN", 0) === 3:4
        @test approxsearch(seq, dna"GN", 1) === 1:1
        @test approxsearch(seq, dna"ACG", 0) === 1:3
        @test approxsearch(seq, dna"ACG", 1) === 1:2
        @test approxsearch(seq, dna"ACG", 2) === 1:1
        @test approxsearch(seq, dna"ACG", 3) === 1:0
        @test approxsearch(seq, dna"ACG", 4) === 1:0

        @test approxsearchindex(seq, dna"", 0) === 1
        @test approxsearchindex(seq, dna"AC", 0) === 1
        @test approxsearchindex(seq, dna"AC", 0, 2) === 5
        @test approxsearchindex(seq, dna"AC", 0, 2, 5) === 0

        query = ApproximateSearchQuery(dna"ACG")
        @test approxsearch(seq, query, 1) === 1:2
        @test approxsearch(seq, query, 1, 2) === 2:3
        @test approxsearch(seq, query, 1, 2, 2) === 0:-1
        @test approxsearchindex(seq, query, 1) === 1
        @test approxsearchindex(seq, query, 1, 2) === 2
        @test approxsearchindex(seq, query, 1, 2, 2) === 0
    end

    @testset "backward" begin
        # TODO: maybe this should return 8:7 like rsearch
        @test approxrsearch(seq, dna"", 0) === 7:6
        @test approxrsearch(seq, dna"AC", 0) === 5:6
        @test approxrsearch(seq, dna"AC", 0, 5) === 1:2
        @test approxrsearch(seq, dna"AC", 0, 5, 2) === 0:-1
        @test approxrsearch(seq, dna"AT", 0) === 0:-1
        @test approxrsearch(seq, dna"AT", 1) === 5:6
        @test approxrsearch(seq, dna"AT", 1, 6) === 5:6
        @test approxrsearch(seq, dna"AT", 1, 5) === 5:5
        @test approxrsearch(seq, dna"AT", 1, 3, 2) === 0:-1
        @test approxrsearch(seq, dna"NG", 0) === 6:7
        @test approxrsearch(seq, dna"NG", 0, 6) === 2:3
        @test approxrsearch(seq, dna"GN", 0) === 3:4
        @test approxrsearch(seq, dna"GN", 1) === 7:7
        @test approxrsearch(seq, dna"ACG", 0) === 5:7
        @test approxrsearch(seq, dna"ACG", 1) === 6:7
        @test approxrsearch(seq, dna"ACG", 2) === 7:7
        @test approxrsearch(seq, dna"ACG", 3) === 7:6
        @test approxrsearch(seq, dna"ACG", 4) === 7:6

        # TODO: maybe this should return 8 like rsearchindex
        @test approxrsearchindex(seq, dna"", 0) === 7
        @test approxrsearchindex(seq, dna"AC", 0) === 5
        @test approxrsearchindex(seq, dna"AC", 0, 5) === 1
        @test approxrsearchindex(seq, dna"AC", 0, 5, 2) === 0

        query = ApproximateSearchQuery(dna"ACG")
        @test approxrsearch(seq, query, 1, 7) === 6:7
        @test approxrsearch(seq, query, 1, 6) === 5:6
        @test approxrsearch(seq, query, 1, 6, 6) === 0:-1
        @test approxrsearchindex(seq, query, 1, 7) === 6
        @test approxrsearchindex(seq, query, 1, 6) === 5
        @test approxrsearchindex(seq, query, 1, 6, 6) === 0
    end
end


text:
#!/usr/bin/env julia

"""
    get_transcripts_from_gff(fa::String, gff::String)

    Computes and returns array of transcripts from reference genome and gff3 file.

    # Arguments
    - `fa::String`: path to fasta file of reference genome
    - `gff::String`: path to gff3 transcripts file
"""
function get_transcripts_from_gff(fa::String, gff::String)::Array{FASTA.Record,1}

    # Reads Gff file generates dict chr => [transcriptsPositions]
    println("Reading GFF")
    reader = open(GFF3.Reader, gff)
    transdict = Dict{String,Array{String,1}}()
    ct = 0
    num = 0
    currtrans = ""
    coordex = ""
    strnd = ""
    chr = ""
    tic()

    for record in reader
        featureType = GFF3.featuretype(record)
        # Check if new feature found, then push exons
        if featureType == "transcript" || featureType == "mRNA"
            if currtrans != string(GFF3.attributes(record, "ID")[1]) && currtrans != ""
                if coordex == ""
                    println(STDERR,"ERROR! Exon is empty!")
                else
                    coordex = coordex * strnd
                end

                if !(in(coordex, transdict[chr]))
                    push!(transdict[chr],coordex)
                    num += 1
                else
                    ct += 1
                end
            end
            chr = string(GFF3.seqid(record))
            currtrans = string(GFF3.attributes(record, "ID")[1])
            strnd = string(GFF3.strand(record))
            # Check if chr exist in dict
            if !(haskey(transdict,chr))
                transdict[chr] = []
            end
            coordex = ""

        elseif featureType == "exon"
            strpos = string(GFF3.seqstart(record))
            endpos = string(GFF3.seqend(record))

            if currtrans == GFF3.attributes(record, "Parent")[1]
                coordex = coordex * strpos * "-" * endpos * ","
            else
                println(STDERR,"ERROR! Exon before transcript.")
                println(record)
            end
        end
    end

    # Adds If exon was last record in gff3 file.
    if coordex != ""
        coordex = coordex * strnd

        if !(in(coordex, transdict[chr]))
            push!(transdict[chr],coordex)
            num += 1
        else
            ct += 1
        end
    end

    toc()
    println(STDOUT,"Dublicated transcripts count in gff = ",ct)
    println(STDOUT,"Number of transcripts in gff = ",num)
    close(reader)

    println(STDOUT,"Collecting sequences from ref genome using data from GFF")
    tic()
    # Read ref genome FASTA File and optains transcripts sequences
    file_stream = open(fa,"r")
    result = @parallel (vcat) for record in collect(FASTA.Reader(file_stream))
        get_transcripts_from_dict(record, transdict)
    end
    toc()
    close(file_stream)

    return result
end

"""
    get_transcripts_from_dict(record::BioSequences.FASTA.Record, transdict::Dict)

    Function returns transcripts list as an array of FASTA records.

    # Arguments
    - `record::BioSequences.FASTA.Record`: fasta record of sequence
    - `transdict::Dict`: dictionary with transcripts coordinates
"""
function get_transcripts_from_dict(record::FASTA.Record,
    transdict::Dict)::Array{FASTA.Record,1}

    outrecords = Array{FASTA.Record,1}()
    chr = FASTA.identifier(record)

    if haskey(transdict,chr)
        # Reads transcripts in Chromosome from transcripts dict
        d = transdict[chr]

        for trans in d
            transseq = dna""
            transcoordsandstrand = split(trans,",")
            transcoords = transcoordsandstrand[1:end-1]
            strand = transcoordsandstrand[end]

            for coords in transcoords
                coords = split(coords, "-")
                coor1 = parse(Int, coords[1])
                coor2 = parse(Int, coords[2])
                try
                    transseq = transseq * FASTA.sequence(record,coor1:coor2)
                catch error
                    println(STDERR, "ERROR! $error")
                    continue
                end
            end
            # If strand - makes reverse complement
            if strand == "-"
                transseq = reverse_complement!(transseq)
            end
            if !(transseq == dna"")
                push!(outrecords, FASTA.Record("Seq", transseq))
            end
        end
    else
        println(STDERR,"ERROR! Chromosome '",chr, "' is not in GFF file.")
    end

    return outrecords
end

ext: jl

max_stars_repo_name: ThermofisherAndBabraham/polyAterminus

max_issues_repo_name: ThermofisherAndBabraham/polyAterminus

max_forks_repo_name: ThermofisherAndBabraham/polyAterminus


text:
using Images, FileIO
using ImageTransformations
using CairoMakie
using JLD2
using WindMixing: plot_profiles_fluxes_final

u_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\u.png"))))
v_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\v.png"))))
T_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\T.png"))))
uw_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\uw.png"))))
vw_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\vw.png"))))
wT_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\wT.png"))))
Ri_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\Ri.png"))))
z_img = load(assetpath(joinpath(pwd(), "images_plot\\z.png")))
T_3D_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\T_3D.png"))))

axis_images = (
    u = u_img,
    v = v_img,
    T = T_img,
    z = z_img,
    uw = uw_img,
    vw = vw_img,
    wT = wT_img,
    Ri = Ri_img,
    T_3D = T_3D_img,
)

RESULTS_DIR = "C:\\Users\\xinle\\Documents\\OceanParameterizations.jl\\final_results"
NDE_DIR = "18sim_old"

filename = "train_wind_-5e-4_heating_-3e-8_new"

file = jldopen(joinpath(RESULTS_DIR, NDE_DIR, filename, "profiles_fluxes_oceananigans.jld2"))

plot_data = file["NDE_profile"]

close(file)

frame = 1009

plot_profiles_fluxes_final(plot_data, frame, axis_images, "final_results\\profiles_fluxes_$(filename)_$(frame).png")


text:
const CALLBACK_DEBUG_SIGNAL = Ref(false)

function debug_callback(args...)
    CALLBACK_DEBUG_SIGNAL[] = true
    default_debug_callback(args...)
end

const debug_callback_c = @cfunction(debug_callback, UInt32, (DebugUtilsMessageSeverityFlagEXT, DebugUtilsMessageTypeFlagEXT, Ptr{VkCore.VkDebugUtilsMessengerCallbackDataEXT}, Ptr{Cvoid}))
const API_VERSION = v"1.2"
const VALIDATION_LAYER = "VK_LAYER_KHRONOS_validation"

const INSTANCE_LAYERS = String[
]
const INSTANCE_EXTENSIONS = String[
]
const DEVICE_EXTENSIONS = String[
]
const ENABLED_FEATURES = PhysicalDeviceFeatures(
)

let available_layers = unwrap(enumerate_instance_layer_properties())
    if VALIDATION_LAYER ∈ getproperty.(available_layers, :layer_name)
        push!(INSTANCE_LAYERS, VALIDATION_LAYER)
        true
    else
        @warn "Validation layer not found."
        false
    end
end

const WITH_DEBUG = let available_extensions = unwrap(enumerate_instance_extension_properties())
    if "VK_EXT_debug_utils" ∈ getproperty.(available_extensions, :extension_name)
        push!(INSTANCE_EXTENSIONS, "VK_EXT_debug_utils")
        true
    else
        @warn "VK_EXT_debug_utils not supported"
        false
    end
end

@testset "Vulkan tests" begin
    include("init.jl")

    @testset "Utilities" begin
        @testset "Function pointers" begin
            @test function_pointer(instance, "vkEnumeratePhysicalDevices") == get_instance_proc_addr("vkEnumeratePhysicalDevices"; instance)
            @test function_pointer(device, "vkQueueWaitIdle") == get_device_proc_addr(device, "vkQueueWaitIdle")
        end
    end

    @testset "Debugging" begin
        @testset "Validation" begin
            if WITH_DEBUG
                @test CALLBACK_DEBUG_SIGNAL[]
            end
        end
    end

    @testset "Bitmask flags" begin
        buffer = Buffer(device, 24, BUFFER_USAGE_INDIRECT_BUFFER_BIT, SHARING_MODE_EXCLUSIVE, [0]; flags=BUFFER_CREATE_PROTECTED_BIT & BUFFER_CREATE_SPARSE_ALIASED_BIT)
        buffer_2 = Buffer(device, 24, BUFFER_USAGE_SHADER_DEVICE_ADDRESS_BIT, SHARING_MODE_EXCLUSIVE, [0]; flags=0)
        @test buffer isa Buffer
        @test buffer_2 isa Buffer
    end

    @testset "Buffers" begin
        include("buffers.jl")
    end

    @testset "Introspection" begin
        @test Vk.hl_type(VkCore.VkPhysicalDeviceFeatures) == PhysicalDeviceFeatures
        @test Vk.intermediate_type(VkCore.VkPhysicalDeviceFeatures) == Vk.intermediate_type(PhysicalDeviceFeatures) == _PhysicalDeviceFeatures
        @test Vk.core_type(PhysicalDeviceFeatures) == VkCore.VkPhysicalDeviceFeatures
        @test Vk.structure_type(PhysicalDeviceFeatures2) ==
            Vk.structure_type(_PhysicalDeviceFeatures2) ==
            Vk.structure_type(VkCore.VkPhysicalDeviceFeatures2) ==
            VkCore.VK_STRUCTURE_TYPE_PHYSICAL_DEVICE_FEATURES_2
    end

    @testset "Next chains" begin
        # Initialization of chained empty structs.
        f1 = Vk.initialize(_PhysicalDeviceFeatures2, _PhysicalDeviceVulkan12Features, _PhysicalDeviceVulkanMemoryModelFeatures)
        @test f1.vks.pNext ≠ C_NULL
        @test Base.unsafe_load(Base.unsafe_convert(Ptr{VkCore.VkPhysicalDeviceVulkan12Features}, f1.vks.pNext)).pNext ≠ C_NULL

        f1 = Vk.initialize(PhysicalDeviceFeatures2, PhysicalDeviceVulkan12Features, PhysicalDeviceVulkanMemoryModelFeatures)
        @test f1.next isa PhysicalDeviceVulkan12Features
        @test f1.next.next isa PhysicalDeviceVulkanMemoryModelFeatures

        # Chaining/unchaining of structs.
        chain_args = [PhysicalDeviceFeatures2(PhysicalDeviceFeatures()), PhysicalDeviceVulkan12Features(), PhysicalDeviceVulkanMemoryModelFeatures(false, false, false)]
        chained = Vk.chain(chain_args...)
        @test Vk.unchain(chained) == chain_args
        @test f1 == chained

        # Conversion to core data structures and back.
        f2 = _PhysicalDeviceFeatures2(f1)
        f3 = Base.unsafe_convert(VkCore.VkPhysicalDeviceFeatures2, f2)
        f4 = PhysicalDeviceFeatures2(f3, PhysicalDeviceVulkan12Features, PhysicalDeviceVulkanMemoryModelFeatures)
        @test f1 == f4

        # Queries with optional chain members.
        feats = get_physical_device_features_2(device.physical_device, PhysicalDeviceVulkan12Features, PhysicalDeviceVulkanMemoryModelFeatures)
        @test feats isa PhysicalDeviceFeatures2
        @test feats.next isa PhysicalDeviceVulkan12Features
        @test feats.next.next isa PhysicalDeviceVulkanMemoryModelFeatures
        feats2 = get_physical_device_features_2(device.physical_device, PhysicalDeviceVulkan12Features)
        @test feats2 == @set(feats2.next.next = C_NULL)

        props = get_physical_device_properties_2(device.physical_device, PhysicalDeviceProtectedMemoryProperties, PhysicalDeviceProvokingVertexPropertiesEXT)
        @test props isa PhysicalDeviceProperties2
        @test props.next isa PhysicalDeviceProtectedMemoryProperties
        @test props.next.next isa PhysicalDeviceProvokingVertexPropertiesEXT
        props = get_physical_device_properties_2(device.physical_device, PhysicalDeviceVulkan11Properties, PhysicalDeviceVulkan12Properties)
        @test props.properties.api_version ≥ v"1.1"
    end

    @testset "Handle constructors" begin
        info = FenceCreateInfo()
        fence_1 = Fence(device, info)
        @test fence_1 isa Fence
        fence_2 = Fence(device)
        @test fence_2 isa Fence
        fence_3 = unwrap(create_fence(device, info))
        @test fence_3 isa Fence
        fence_4 = unwrap(create_fence(device))
        @test fence_4 isa Fence

        info = _FenceCreateInfo()
        fence_1 = Fence(device, info)
        @test fence_1 isa Fence
        fence_2 = Fence(device)
        @test fence_2 isa Fence
        fence_3 = unwrap(_create_fence(device, info))
        @test_throws MethodError unwrap(create_fence(device, info))
        @test fence_3 isa Fence
        fence_4 = unwrap(_create_fence(device))
        @test fence_4 isa Fence
    end
end

GC.gc()

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
#!/usr/bin/env julia

#
# Task ventilator
# Binds PUSH socket to tcp://localhost:5557
# Sends batch of tasks to workers via that socket
#

using ZMQ

context = Context()
sender = Socket(context, PUSH)
ZMQ.bind(sender, "tcp://*:5557")

# Socket with direct access to the sink: used to syncronize start of batch
sink = Socket(context, PUSH)
ZMQ.connect(sink, "tcp://localhost:5558")

println("Press Enter when the workers are ready: ")
_ = readline(STDIN)
println("Sending tasks to workers...")

# The first message is "0" and signals start of batch
ZMQ.send(sink, b"0")

# Initialize random number generator
srand(1)

# Send 100 tasks
total_msec = 0
for task_nbr in [1:100]
    # Random workload from 1 to 100 msecs
    workload = rand(1:100)
    total_msec += workload

    ZMQ.send(sender, "$workload")
end

println("Total expected cost: $total_msec msec")

# Give 0MQ time to deliver
sleep(1)

# Making a clean exit.
ZMQ.close(sender)
ZMQ.close(sink)
ZMQ.close(context)

max_stars_repo_name: JailbreakFox/LightWeightRepository

max_issues_repo_name: JailbreakFox/LightWeightRepository

max_forks_repo_name: JailbreakFox/LightWeightRepository


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

"Represents a Basic Block, in the DomTree"
struct DomTreeNode
    # How deep we are in the DomTree
    level::Int
    # The BB indices in the CFG for all Basic Blocks we immediately dominate
    children::Vector{Int}
end
DomTreeNode() = DomTreeNode(1, Vector{Int}())

"Data structure that encodes which basic block dominates which."
struct DomTree
    # Which basic block immediately dominates each basic block (ordered by BB indices)
    # Note: this is the inverse of the nodes, children field
    idoms::Vector{Int}

    # The nodes in the tree (ordered by BB indices)
    nodes::Vector{DomTreeNode}
end

"""
    Checks if bb1 dominates bb2.
    bb1 and bb2 are indexes into the CFG blocks.
    bb1 dominates bb2 if the only way to enter bb2 is via bb1.
    (Other blocks may be in between, e.g bb1->bbX->bb2).
"""
function dominates(domtree::DomTree, bb1::Int, bb2::Int)
    bb1 == bb2 && return true
    target_level = domtree.nodes[bb1].level
    source_level = domtree.nodes[bb2].level
    source_level < target_level && return false
    for _ in (source_level - 1):-1:target_level
        bb2 = domtree.idoms[bb2]
    end
    return bb1 == bb2
end

bb_unreachable(domtree::DomTree, bb::Int) = bb != 1 && domtree.nodes[bb].level == 1

function update_level!(domtree::Vector{DomTreeNode}, node::Int, level::Int)
    worklist = Tuple{Int, Int}[(node, level)]
    while !isempty(worklist)
        (node, level) = pop!(worklist)
        domtree[node] = DomTreeNode(level, domtree[node].children)
        foreach(domtree[node].children) do child
            push!(worklist, (child, level+1))
        end
    end
end

"Iterable data structure that walks though all dominated blocks"
struct DominatedBlocks
    domtree::DomTree
    worklist::Vector{Int}
end

"Returns an iterator that walks through all blocks dominated by the basic block at index `root`"
function dominated(domtree::DomTree, root::Int)
    doms = DominatedBlocks(domtree, Vector{Int}())
    push!(doms.worklist, root)
    doms
end

function iterate(doms::DominatedBlocks, state::Nothing=nothing)
    isempty(doms.worklist) && return nothing
    bb = pop!(doms.worklist)
    for dominated in doms.domtree.nodes[bb].children
        push!(doms.worklist, dominated)
    end
    return (bb, nothing)
end

function naive_idoms(cfg::CFG)
    nblocks = length(cfg.blocks)
    # The extra +1 helps us detect unreachable blocks below
    dom_all = BitSet(1:nblocks+1)
    dominators = BitSet[n == 1 ? BitSet(1) : copy(dom_all) for n = 1:nblocks]
    changed = true
    while changed
        changed = false
        for n = 2:nblocks
            if isempty(cfg.blocks[n].preds)
                continue
            end
            firstp, rest = Iterators.peel(Iterators.filter(p->p != 0, cfg.blocks[n].preds))
            new_doms = copy(dominators[firstp])
            for p in rest
                intersect!(new_doms, dominators[p])
            end
            push!(new_doms, n)
            changed = changed || (new_doms != dominators[n])
            dominators[n] = new_doms
        end
    end
    # Compute idoms
    idoms = fill(0, nblocks)
    for i = 2:nblocks
        if dominators[i] == dom_all
            idoms[i] = 0
            continue
        end
        doms = collect(dominators[i])
        for dom in doms
            i == dom && continue
            hasany = false
            for p in doms
                if p !== i && p !== dom && dom in dominators[p]
                    hasany = true; break
                end
            end
            hasany && continue
            idoms[i] = dom
        end
    end
    idoms
end

# Construct Dom Tree
function construct_domtree(cfg::CFG)
    idoms = SNCA(cfg)
    # Compute children
    nblocks = length(cfg.blocks)
    domtree = DomTreeNode[DomTreeNode() for _ = 1:nblocks]
    for (idx, idom) in Iterators.enumerate(idoms)
        (idx == 1 || idom == 0) && continue
        push!(domtree[idom].children, idx)
    end
    # Recursively set level
    update_level!(domtree, 1, 1)
    DomTree(idoms, domtree)
end

#================================ [SNCA] ======================================#
#
#   This section implements the Semi-NCA (SNCA) dominator tree construction from
#   described in Georgiadis' PhD thesis [LG05], which itself is a simplification
#   of the Simple Lenguare-Tarjan (SLT) algorithm [LG79]. This algorithm matches
#   the algorithm choice in LLVM and seems to be a sweet spot in implementation
#   simplicity and efficiency.
#
#   [LG05]  Linear-Time Algorithms for Dominators and Related Problems
#           Loukas Georgiadis, Princeton University, November 2005, pp. 21-23:
#           ftp://ftp.cs.princeton.edu/reports/2005/737.pdf
#
#   [LT79]  A fast algorithm for finding dominators in a flowgraph
#           Thomas Lengauer, Robert Endre Tarjan, July 1979, ACM TOPLAS 1-1
#           http://www.dtic.mil/dtic/tr/fulltext/u2/a054144.pdf
#
begin
    # We could make these real structs, but probably not worth the extra
    # overhead. Still, give them names for documentary purposes.
    const BBNumber = UInt
    const DFSNumber = UInt

    """
    Keeps the per-BB state of the Semi NCA algorithm. In the original
    formulation, there are three separate length `n` arrays, `label`, `semi` and
    `ancestor`. Instead, for efficiency, we use one array in a array-of-structs
    style setup.
    """
    struct Node
        semi::DFSNumber
        label::DFSNumber
    end

    struct DFSTree
        # Maps DFS number to BB number
        numbering::Vector{BBNumber}
        # Maps BB number to DFS number
        reverse::Vector{DFSNumber}
        # Records parent relationships in the DFS tree (DFS number -> DFS number)
        # Storing it this way saves a few lookups in the snca_compress! algorithm
        parents::Vector{DFSNumber}
    end
    length(D::DFSTree) = length(D.numbering)
    preorder(D::DFSTree) = OneTo(length(D))
    _drop(xs::AbstractUnitRange, n::Integer) = (first(xs)+n):last(xs)

    function DFSTree(nblocks::Int)
        DFSTree(
            Vector{BBNumber}(undef, nblocks),
            zeros(DFSNumber, nblocks),
            Vector{DFSNumber}(undef, nblocks))
    end

    function DFS(cfg::CFG, current_node::BBNumber)::DFSTree
        dfs = DFSTree(length(cfg.blocks))
        # TODO: We could reuse the storage in DFSTree for our worklist. We're
        # guaranteed for the worklist to be smaller than the remaining space in
        # DFSTree
        worklist = Tuple{DFSNumber, BBNumber}[(0, current_node)]
        dfs_num = 1
        parent = 0
        while !isempty(worklist)
            (parent, current_node) = pop!(worklist)
            dfs.reverse[current_node] != 0 && continue
            dfs.reverse[current_node] = dfs_num
            dfs.numbering[dfs_num] = current_node
            dfs.parents[dfs_num] = parent
            for succ in cfg.blocks[current_node].succs
                push!(worklist, (dfs_num, succ))
            end
            dfs_num += 1
        end
        # If all blocks are reachable, this is a no-op, otherwise,
        # we shrink these arrays.
        resize!(dfs.numbering, dfs_num - 1)
        resize!(dfs.parents, dfs_num - 1)
        dfs
    end

    """
    Matches the snca_compress algorithm in Figure 2.8 of [LG05], with the
    modification suggested in the paper to use `last_linked` to determine
    whether an ancestor has been processed rather than storing `0` in the
    ancestor array.
    """
    function snca_compress!(state::Vector{Node}, ancestors::Vector{DFSNumber},
                            v::DFSNumber, last_linked::DFSNumber)
        u = ancestors[v]
        @assert u < v
        if u >= last_linked
            snca_compress!(state, ancestors, u, last_linked)
            if state[u].label < state[v].label
                state[v] = Node(state[v].semi, state[u].label)
            end
            ancestors[v] = ancestors[u]
        end
        nothing
    end

    function snca_compress_worklist!(
            state::Vector{Node}, ancestors::Vector{DFSNumber},
            v::DFSNumber, last_linked::DFSNumber)
        # TODO: There is a smarter way to do this
        u = ancestors[v]
        worklist = Tuple{DFSNumber, DFSNumber}[(u,v)]
        @assert u < v
        while !isempty(worklist)
            u, v = last(worklist)
            if u >= last_linked
                if ancestors[u] >= last_linked
                    push!(worklist, (ancestors[u], u))
                    continue
                end
                if state[u].label < state[v].label
                    state[v] = Node(state[v].semi, state[u].label)
                end
                ancestors[v] = ancestors[u]
            end
            pop!(worklist)
        end
    end

    """
        SNCA(cfg::CFG)

    Determines a map from basic blocks to the block which immediately dominate them.
    Expressed as indexes into `cfg.blocks`.

    The main Semi-NCA algrithm. Matches Figure 2.8 in [LG05].
    Note that the pseudocode in [LG05] is not entirely accurate.
    The best way to understand what's happening is to read [LT79], then the
    description of SLT in in [LG05] (warning: inconsistent notation), then
    the description of Semi-NCA.
    """
    function SNCA(cfg::CFG)
        D = DFS(cfg, BBNumber(1))
        # `label` is initialized to the identity mapping (though
        # the paper doesn't make that clear). The rational for this is Lemma
        # 2.4 in [LG05] (i.e. Theorem 4 in ). Note however, that we don't
        # ever look at `semi` until it is fully initialized, so we could leave
        # it uninitialized here if we wanted to.
        state = Node[ Node(typemax(DFSNumber), w) for w in preorder(D) ]
        # Initialize idoms to parents. Note that while idoms are eventually
        # BB indexed, we keep it DFS indexed until a final post-processing
        # pass to avoid extra memory references during the O(N^2) phase below.
        idoms_dfs = copy(D.parents)
        # We abuse the parents array as the ancestors array.
        # Semi-NCA does not look at the parents array at all.
        # SLT would, but never simultaneously, so we could still
        # do this.
        ancestors = D.parents
        for w::DFSNumber ∈ reverse(_drop(preorder(D), 1))
            # LLVM initializes this to the parent, the paper initializes this to
            # `w`, but it doesn't really matter (the parent is a predecessor,
            # so at worst we'll discover it below). Save a memory reference here.
            semi_w = typemax(DFSNumber)
            for v ∈ cfg.blocks[D.numbering[w]].preds
                # For the purpose of the domtree, ignore virtual predecessors
                # into catch blocks.
                v == 0 && continue
                vdfs = D.reverse[v]
                # Ignore unreachable predecessors
                vdfs == 0 && continue
                last_linked = DFSNumber(w + 1)
                # N.B.: This conditional is missing from the psuedocode
                # in figure 2.8 of [LG05]. It corresponds to the
                # `ancestor[v] != 0` check in the `eval` implementation in
                # figure 2.6
                if vdfs >= last_linked
                    # For performance, if the number of ancestors is small
                    # avoid the extra allocation of the worklist.
                    if length(ancestors) <= 32
                        snca_compress!(state, ancestors, vdfs, last_linked)
                    else
                        snca_compress_worklist!(state, ancestors, vdfs, last_linked)
                    end
                end
                semi_w = min(semi_w, state[vdfs].label)
            end
            state[w] = Node(semi_w, semi_w)
        end
        for v ∈ _drop(preorder(D), 1)
            idom = idoms_dfs[v]
            vsemi = state[v].semi
            while idom > vsemi
                idom = idoms_dfs[idom]
            end
            idoms_dfs[v] = idom
        end
        # Reexpress the idom relationship in BB indexing
        idoms_bb = Int[ (i == 1 || D.reverse[i] == 0) ? 0 : D.numbering[idoms_dfs[D.reverse[i]]] for i = 1:length(cfg.blocks) ]
        idoms_bb
    end
end

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
module MultiScaleEntropy

using RecurrenceAnalysis

export multiscaleentropy,
       moving_average,
       arithmean, geommean, mixmean

# Local recurrences (excluding self-recurrence)
local_rr(x::AbstractMatrix{Bool}, exclude=1) = [countnz(x[t+exclude:end,t]) for t in (1:size(x)[2])]

# Note: in this function x is assumed to be previously normalized and embedded, radius is absolute
function sampled_local_rr(x, r, scale=1, delayed=false)
    n = div(size(x)[1], scale)
    !delayed && rem(size(x)[1], scale) > 0 && (n += 1)
    recurrences = zeros(n, delayed?scale:1)
    for s = 1:size(recurrences)[2]
        xs = x[s:scale:scale*n,:]
        rmat = recurrencematrix(xs, r, normalize=false)
        recurrences[:,s] = local_rr(rmat)
    end
    recurrences
end

function full_local_rr(x, radius, scale)
    rmat = recurrencematrix(x, radius, normalize=false)
    recurrences = hcat(local_rr(rmat, scale))
end

# Auxiliary functions
moving_average(x, ns) = filt(ones(ns)/ns, 1, x)[ns:end]
arithmean(x) = mean(x)
geommean(x) = prod(x)^(1/length(x))
mixmean(x) = geommean(mean(x, 1))

function multiscaleentropy(x, m, r, bounds; filtf=moving_average, normbyscale=false, rcount="simple", avg=arithmean)
    # Select methods
    countrecurrences = Dict(
        "simple"    => sampled_local_rr,
        "full"      => full_local_rr,
        "composite" => (x,r,s)->sampled_local_rr(x,r,s,true)
    )[rcount]
    ent = zeros(bounds[2]-bounds[1]+1)
    # Normalize signal
    x ./= std(x)
    for s in range(bounds...)
        # 1. Filter signal by scale
        xf = filtf(x, s)
        normbyscale && (xf ./=std(xf))
        # 2. Embed (and trim) for m and m+1
        xe = embed(xf, m, s)
        xe_plus = [xe[1:end-s,:] xe[1+s:end,end]]
        xe = xe[1:end-s,:]
        # 3. Count recurrences
        rc = countrecurrences(xe, r, s)
        rc_plus = countrecurrences(xe_plus, r, s)
        # 4. Calculate factors of entropy
        # (add self-recurrence if using geometric mean)
        (avg == geommean) && (rc += 1; rc_plus += 1)
        a = avg(rc_plus)
        b = avg(rc)
        ent[s] = -log(a/b)
    end
    ent
end

end


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

#goldschmidt-digits.jl

include("../unum.jl")
using Unums

#goldschmidt algorithm test on Float16, Float32, Float64.  The purpose of this is to
#determine if the understanding of where to put ULPs is correct.

#returns a Int16 "1" if it's negative, "0" if it's positive.
function signof(x::FloatingPoint)
  T = typeof(x)
  if (T == Float64)
    return ((reinterpret(UInt64, x) & 0x8000_0000_0000_0000) != 0) ? 1 : 0
  elseif (T == Float32)
    return ((reinterpret(UInt32, x) & 0x8000_0000) != 0) ? 1 : 0
  end
end

function exponentof(x::FloatingPoint)
  T = typeof(x)
  if (T == Float64)
    return Int64((reinterpret(Int64, x) & 0x7FF0_0000_0000_0000) >> 52 - 1023)
  elseif (T == Float32)
    return Int64((reinterpret(Int32, x) & 0x7F80_0000) >> 23 - 127)
  end
end

function castfrac(x::FloatingPoint)
  T = typeof(x)
  if (T == Float64)
    return (reinterpret(UInt64, x) & 0x000F_FFFF_FFFF_FFFF) << 12
  elseif (T == Float32)
    return UInt64(reinterpret(UInt32, x) & 0x007F_FFFF) << 41
  end
end

function maskfor(T::Type)
  if (T == Float64)
    return 0xFFFF_FFFF_FFFF_FF00  #56 digits suffices for Float64 (52 digits)
  elseif (T == Float32)
    return 0xFFFF_FFE0_0000_0000  #27 digits of precision seems to suffice for a Float32 (23 digits)
  end
end

o16 = one(UInt16)
z16 = zero(UInt16)
__clz_array=[0x0004,0x0003,0x0002,0x0002, o16, o16, o16, o16, z16,z16,z16,z16,z16,z16,z16,z16]
function leading_zeros(n)
  (n == 0) && return 64
  res::UInt16 = 0
  #use the binary search method
  (n & 0xFFFF_FFFF_0000_0000 == 0) && (n <<= 32; res += 0x0020)
  (n & 0xFFFF_0000_0000_0000 == 0) && (n <<= 16; res += 0x0010)
  (n & 0xFF00_0000_0000_0000 == 0) && (n <<= 8;  res += 0x0008)
  (n & 0xF000_0000_0000_0000 == 0) && (n <<= 4;  res += 0x0004)
  res + __clz_array[(n >> 60) + 1]
end

#simple fused - multiply - add.  Assumes num1 has hidden bits "carry" and num2
#has no hidden bits.
rm = 0x0000_0000_FFFF_FFFF
function sfma(carry, num1, num2)
  (fracprod, _) = Unums.__chunk_mult(num1, num2)
  (_carry, fracprod) = Unums.__carried_add(carry, num1, fracprod)
  ((carry & 0x1) != 0) && ((_carry, fracprod) = Unums.__carried_add(_carry, num2, fracprod))
  ((carry & 0x2) != 0) && ((_carry, fracprod) = Unums.__carried_add(_carry, lsh(num2, 1), fracprod))
  (_carry, fracprod)
end

#performs a simple multiply, Assumes that number 1 has a hidden bit of exactly one
#and number 2 has a hidden bit of exactly zero
#(1 + a)(0 + b) = b + ab
function smult(a::UInt64, b::UInt64)

  (fraction, _) = Unums.__chunk_mult(a, b)
  carry = one(UInt64)

  #only perform the respective adds if the *opposing* thing is not subnormal.
  ((carry, fraction) = Unums.__carried_add(carry, fraction, b))

  #carry may be as high as three!  So we must shift as necessary.
  (fraction, shift, is_ubit) = Unums.__shift_after_add(carry, fraction, _)
  fraction << 1
end

function reassemble(T::Type, sign, exp_f, number)
  if T == Float64
    number = (number >> 12) | (exp_f + 1023) << 52 | sign << 63
  elseif T == Float32
    number = (number >> 9) | (exp_f + 127) << 55 | sign << 63
  end
  last(reinterpret(T, [number]))
end

function calculate(carry, number)
  carry + big(number) / (big(1) << 64)
end

function maxexp(T)
  if T == Float64
    1023
  elseif T == Float32
    127
  end
end

function minexp(T, subnormalq)
  if T == Float64
    -1022 - (subnormalq ? 52 : 0)
  elseif T == Float32
    -126 - (subnormalq ? 23 : 0)
  end
end

function findshift(shift)
  1 << (64 - shift)
end

function fixsn(T, exp_f, frac)
  if T == Float64
    shift = -1022 - exp_f
    return (-1023, (frac >> shift | findshift(shift)))
  elseif T == Float32
    shift = -126 - exp_f
    return (-127, (frac >> shift | findshift(shift)))
  end
end

function exct(x::FloatingPoint, y::FloatingPoint)
  #exact floating point with the goldschmidt algorithm.
  T = typeof(x)

  if isnan(x) || isnan(y)
    return (nan(T), false)
  end

  #figure out the sign.
  sign::UInt64 = signof(x) $ signof(y)

  #calculate the exponent.
  exp_f::Int64 = exponentof(x) - exponentof(y) + (issubnormal(x) ? 1 : 0) - (issubnormal(y) ? 1 : 0)

  #figure the decimals.
  numerator::UInt64 = castfrac(x)

  if (issubnormal(x))
<<<<<<< HEAD:oddsandends/goldschmidt-digits.jl
    shift::UInt64 = clz(numerator) + 1
=======
    shift::UInt64 = leading_zeros(numerator) + 1
>>>>>>> directory-reorg:src/oddsandends/goldschmidt-digits.jl
    numerator = numerator << shift
    exp_f -= shift
  end
  #save the old numerator
  old_numerator = numerator
  #set the carry on the numerator
  carry::UInt64 = 1

  denominator::UInt64 = castfrac(y)
  #adjust the denominator in the case that it's a subnormal
  if issubnormal(y)
    shift = leading_zeros(denominator)
    denominator = denominator << shift
    exp_f += shift
  else
    denominator = (denominator >> 1) | 0x8000_0000_0000_0000
    exp_f -= 1
  end
  #and then save this old denominator.
  old_denominator = denominator

  (exp_f > maxexp(T)) && (return (sign == 1 ? convert(T,-Inf) : inf(T), false))
  (exp_f < minexp(T, true) - 2) && (return ((sign == 1 ? convert(T,-0.0) : 0.0), false))

  ourfrac_mask = maskfor(T)
  #do the goldschmidt algorithm
  for (idx = 1:32)
    factor::UInt64 = (-denominator)
    #simple-fused-multiply-add.
    (carry, numerator) = sfma(carry, numerator, factor)
    (_, denominator) = sfma(UInt64(0), denominator, factor)
    (~denominator & ourfrac_mask == 0) && break
    denominator &= ourfrac_mask
    numerator &= ourfrac_mask
  end
  if carry > 1
    numerator = numerator >> 1 | (carry & 0x1 << 63)
    exp_f += 1
  end

  (exp_f > maxexp(T)) && (return ((sign == 1 ? convert(T,-Inf) : inf(T)), false))
  (exp_f < minexp(T, true)) && (return ((sign == 1 ? convert(T,-0.0) : 0.0), false))

  #now, at this point, we have to go back and 'check our work.'
  #multiply the old_denominator times the current numerator, and should result
  #in the old_numerator

  numerator &= ourfrac_mask
  ans_subnormal = exp_f < minexp(T, false)
  is_ulp = true

  f_d = 0x0000_0000_0000_1000
  f_ma = 0xFFFF_FFFF_FFFF_F000
  #do a "remultiply" operation.  First attempt with the lower unit.

  reseq = smult(numerator & f_ma, old_denominator)
  resph = smult((numerator & f_ma + f_d), old_denominator)

  if (old_numerator < reseq)
    numerator = (numerator - f_d)
  elseif (old_numerator == reseq)
    #need to run an ulp check here.
    numerator = (numerator - f_d)
  elseif (old_numerator > (resph))
    numerator = (numerator + f_d)
  end

  (exp_f < minexp(T, false)) && ((exp_f, numerator) = fixsn(T, exp_f, numerator))

  #reassemble the value into the requisite floating point
  (reassemble(T, sign, exp_f, numerator), is_ulp)
end
#=
#one-time testing
#x = reinterpret(Float32, 0b00000110001001000111001101001111)
#y = reinterpret(Float32, 0b10000000010010100000111001111111)

x = reinterpret(Float64, rand(UInt64))
y = reinterpret(Float64, rand(UInt64))

#test exact divisions
#y = floor(rand() * 100000)
#q = floor(rand() * 100000)
#x = q * y
#println("theo. res: $(q)")

#x = 3.0
#y = 2.0

println("answer:    $(x/y)")
z = exct(x, y)
println("gscalc:    $z")
println("abits:     ", bits(x/y))
println("zbits:     ", bits(z))
#println("theo. bits:", bits(q))
=#

#continuous testing.
count = 0
errors = 0
while (true)
  #x = reinterpret(Float64, rand(UInt64))
  #y = reinterpret(Float64, rand(UInt64))
  x = 30.0
  y = 1.5
  (z, ulp) = exct(x, y)

  bigres = big(x) / big(y)
  bigguess = big(z)

  (isinf(z)) && continue

  if ulp
    nextfrac = big(reinterpret(Float64, (reinterpret(UInt64, z) + 1)))
    if (abs(bigguess) > abs(bigres))
      println("lower bound bad")
      println("gsans:", bits(z))
      println("fpans:", bits(x/y))
      println("count, $count")
      exit()
    end
    if (abs(nextfrac) < abs(bigres))
      println("upper bound bad")
      println(bits(z))
      println(bits(x/y))
      exit()
    end
  end
  count += 1
  #println(count)
end

#NB:  This technique still has a vanishingly small ~(0.005%) error rate in assigning
#the correct ULP to the quotient.  Will need to take a  more careful look at why this
#is occurring.

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using Comonicon, PyPlot
using DelimitedFiles

@cast function hamming(which::String, n::Int, alpha::Float64;
            degree::Int=3)
    for seed=1:100
        @show n, seed
        if which == "regular"
            folder = joinpath("data", "IndependentSet_Regular$(n)d$(degree)seed$(seed)")
        elseif which == "diag"
            folder = joinpath("data", "IndependentSet_Diag$(n)x$(n)f0.8seed$(seed)")
        else
            error("")
        end
        maxn = readdlm(joinpath(folder, "SizeMax1.dat"))[]
        K = ceil(Int, maxn * alpha)

        hamming = readdlm(joinpath(folder, "hamming-K$(K)-n10000.dat"))
        f = plot(collect(0:length(hamming)-1), hamming)

        #ylabel("hamming distance")
        PyPlot.axis("off")
        savefig(joinpath(folder, "hamming-K$(K)-n10000.png"))
        clf()
    end
    #ylabel("hamming distance")
    #savefig(joinpath("hamming-$(which)-size$(n)-n10000.png"))
    #clf()
end

@cast function grid(which::String, n::Int, alpha::Float64;
            degree::Int=3, graphsize::Int=3)
    f, axs = subplots(graphsize, graphsize)
    for seed=1:graphsize^2
        @show n, seed
        if which == "regular"
            folder = joinpath("data", "IndependentSet_Regular$(n)d$(degree)seed$(seed)")
        elseif which == "diag"
            folder = joinpath("data", "IndependentSet_Diag$(n)x$(n)f0.8seed$(seed)")
        else
            error("")
        end
        maxn = readdlm(joinpath(folder, "SizeMax1.dat"))[]
        K = ceil(Int, maxn * alpha)

        hamming = readdlm(joinpath(folder, "hamming-K$(K)-n10000.dat")) |> vec
        axs[seed].bar(collect(0:length(hamming)-1), hamming)
        axs[seed].axis("off")
    end
    fname = "grid-$(which)-size$(n)d$(degree)-k$(graphsize)-alpha$(alpha)-n10000.pdf"
    println("saving to: $fname")
    PyPlot.suptitle("\${\\rm ratio}=$(alpha)\$", fontsize=28, y=0.07)
    savefig(fname)
end

@cast function griddata(which::String, n::Int, alpha::Float64;
            degree::Int=3, graphsize::Int=3)
    hammings = zeros(Int, 1+(which == "regular" ? n : round(Int, n^2 * 0.8)), graphsize^2)
    for seed=1:graphsize^2
        @show n, seed
        if which == "regular"
            folder = joinpath("data", "IndependentSet_Regular$(n)d$(degree)seed$(seed)")
        elseif which == "diag"
            folder = joinpath("data", "IndependentSet_Diag$(n)x$(n)f0.8seed$(seed)")
        else
            error("")
        end
        maxn = readdlm(joinpath(folder, "SizeMax1.dat"))[]
        K = ceil(Int, maxn * alpha)
        @show K

        hamming = readdlm(joinpath(folder, "hamming-K$(K)-n10000.dat")) |> vec
        hammings[:,seed] .= hamming
    end
    fname = "grid-$(which)-size$(n)d$(degree)-k$(graphsize)-alpha$(alpha)-n10000.dat"
    println("saving to: $fname")
    writedlm(fname, hammings)
end

@main


text:
using Documenter, ReferenceTests

makedocs(
    modules = [ReferenceTests],
    clean = false,
    format = :html,
    sitename = "ReferenceTests.jl",
    authors = "Christof Stocker",
    linkcheck = !("skiplinks" in ARGS),
    pages = Any[
        "Home" => "index.md",
    ],
    html_prettyurls = !("local" in ARGS),
)

deploydocs(
    repo = "github.com/Evizero/ReferenceTests.jl.git",
    target = "build",
    julia = "1.0",
    deps = nothing,
    make = nothing,
)


text:
function HMMbackward(v,phghm,pvgh)
    #%HMMBACKWARD HMM Backward Pass
    #% beta=HMMbackward(v,phghm,pvgh)
    #%
    #% Inputs:
    #% v : visible (observation) sequence being a vector v=[2 1 3 3 1 ...]
    #% phghm : homogeneous transition distribution phghm(i,j)=p(h(t)=i|h(t-1)=j)
    #% pvgh : homogeneous emission distribution pvgh(i,j)=p(v(t)=i|h(t)=j)
    #%
    #% Outputs:
    #% beta:  beta messages: \propto p(v(t+1:T)|h(t))
    #           p(v(t:T)|h(t)) = sum_{h(t)} {p(v(t)| h(t)) * p(h(t)|h(t-1)) * p(v(t+1:T)|h(t))}
    #           beta(t-1) = p(v(t:T)|h(t)
    #           beta(t) = p(v(t+1:T)|h(t)
    #% See also HMMbackward.m, HMMviterbi.m, demoHMMinference.m

    T=length(v); H=size(phghm,1)

    pvghtrans=pvgh'
    phghmtrans=phghm'
    beta=ones(H,T)
    for t=T:-1:2
        tmp = phghmtrans * (beta[:,t] .* pvghtrans[:,v[t]])
        beta[:,t-1] = tmp ./ sum(tmp)
    end
    return beta
end


text:
module FastCombinations

export Combinations, MultiCombinations

import Base: start, next, done, length, eltype, collect

# Computes the next $K states of the iterator `itr`, starting from a given state.
# This function can handle cases where K is larger than the length of the given iterator.
@generated function initial_states{K}(itr, ::Type{Val{K}})
    K == 0 && return ()
    syms = [gensym(:state) for _ in 1:K]
    expr = Expr(:block, :($(syms[1]) = start(itr)))
    for i in 2:K
        sym = syms[i]
        push!(expr.args, quote
            (_, $sym) = next(itr, $(syms[i - 1]))
            done(itr, $sym) && return $(Expr(:tuple, repeated(sym, K)...))
        end)
    end
    push!(expr.args, Expr(:tuple, syms...))
    expr
end
# Computes the next $K states of the iterator `itr`, starting from a given state.
# This function assumes that the iterator has at least K elements left
# after the given state.
@generated function next_states{S, K}(itr, state::S, ::Type{Val{K}})
    K == 0 && return ()
    syms = [gensym(:state) for _ in 1:K]
    expr = Expr(:block, :($(syms[1]) = state))
    for i in 2:K
        push!(expr.args, quote
            (_, $(syms[i])) = next(itr, $(syms[i - 1]))
        end)
    end
    push!(expr.args, Expr(:tuple, syms...))
    expr
end

# A k-combination of a set is a way to select exactly k distinct elements from the set.
immutable Combinations{K, T}
    itr::T
    function Combinations(itr::T)
        0 < K || error("Received K <= 0")  # TODO error message
        new{K, T}(itr)
    end
end

Base.convert{K, T}(::Type{Combinations{K}}, itr::T) = Combinations{K, T}(itr)

eltype{K, T}(combs::Combinations{K, T}) = NTuple{K, eltype(combs.itr)}
collect(combs::Combinations) = collect(eltype(combs), combs)
@generated length{K, T}(combs::Combinations{K, T}) = :(binomial(length(combs.itr), $K))

start{K, T}(combs::Combinations{K, T}) = initial_states(combs.itr, Val{K})
done{K, T, State}(combs::Combinations{K, T}, state::NTuple{K, State}) = done(combs.itr, state[end])

@generated function next{K, T, State}(combs::Combinations{K, T}, state::NTuple{K, State})
    comb_elements = [gensym(:comb) for _ in 1:K]
    new_states = [gensym(:state) for _ in 1:K]
    labels = [gensym(Symbol(:label, l)) for l in 1:K]

    next_expr = Expr(:block)
    for i in 1:K
        push!(next_expr.args, :(($(comb_elements[i]), $(new_states[i])) = next(combs.itr, state[$i])))
    end

    goto_block = quote
        current_state = $(new_states[end])
        done(combs.itr, current_state) || @goto $(labels[1])
    end
    for i in 2:K
        push!(goto_block.args, quote
            next_state = $(new_states[end+1-i])
            next(combs.itr, next_state)[2] == current_state || @goto $(labels[i])
            current_state = next_state
        end)
    end
    push!(goto_block.args, :(@goto result))  # Final yield, no additional computation

    label_defs = Expr(:block)
    for i in 1:K
        lhs1 = Expr(:tuple, new_states[1:end-i]...)
        rhs1 = Expr(:tuple, [:(state[$i]) for i in 1:K-i]...)
        lhs2 = Expr(:tuple, new_states[end+2-i:end]...)
        rhs2 = :(next_states(combs.itr, next(combs.itr, $(new_states[end+1-i]))[2], Val{$(i-1)}))

        push!(label_defs.args, quote
            @label $(labels[i])
            $lhs1 = $rhs1
            $lhs2 = $rhs2
            @goto result
        end)
    end

    result_comb = Expr(:tuple, comb_elements...)
    result_state = Expr(:tuple, new_states...)
    result_expr = quote
        @label result
        $result_comb, $result_state
    end
    Expr(:block, next_expr, goto_block, label_defs, result_expr)
end

# a k-multicombination is a combination in which repetition of elements is allowed
# (though the elements must still be in sorted order)
immutable MultiCombinations{K, T}
    itr::T
    function MultiCombinations(itr::T)
        0 < K || error("Received K <= 0")  # TODO error message
        new{K, T}(itr)
    end
end

Base.convert{K, T}(::Type{MultiCombinations{K}}, itr::T) = MultiCombinations{K, T}(itr)

eltype{K, T}(combs::MultiCombinations{K, T}) = NTuple{K, eltype(combs.itr)}
collect(combs::MultiCombinations) = collect(eltype(combs), combs)
@generated length{K, T}(combs::MultiCombinations{K, T}) = :(binomial(length(combs.itr) + $(K - 1), $K))

@generated start{K, T}(combs::MultiCombinations{K, T}) = quote
    state = start(combs.itr)
    $(Expr(:tuple, repeated(:state, K)...))
end
done{K, T, State}(combs::MultiCombinations{K, T}, state::NTuple{K, State}) = done(combs.itr, state[1])

@generated function next{K, T, State}(combs::MultiCombinations{K, T}, state::NTuple{K, State})
    comb_elements = [gensym(:comb) for _ in 1:K]
    new_states = [gensym(:state) for _ in 1:K]
    labels = [gensym(Symbol(:label, l)) for l in 1:K]

    next_expr = Expr(:block)
    for i in 1:K
        push!(next_expr.args, :(($(comb_elements[i]), $(new_states[i])) = next(combs.itr, state[$i])))
    end

    goto_block = Expr(:block)
    for i in 1:K
        push!(goto_block.args, quote
            done(combs.itr, $(new_states[end+1-i])) || @goto $(labels[i])
        end)
    end
    push!(goto_block.args, :(@goto result))

    label_defs = Expr(:block)
    for i in 1:K
        lhs1 = Expr(:tuple, new_states[1:end-i]...)
        rhs1 = Expr(:tuple, [:(state[$i]) for i in 1:K-i]...)
        lhs2 = Expr(:tuple, new_states[end+2-i:end]...)
        rhs2 = Expr(:tuple, repeated(new_states[end+1-i], i-1)...)

        push!(label_defs.args, quote
            @label $(labels[i])
            $lhs1 = $rhs1
            $lhs2 = $rhs2
            @goto result
        end)
    end

    result_comb = Expr(:tuple, comb_elements...)
    result_state = Expr(:tuple, new_states...)
    result_expr = quote
        @label result
        $result_comb, $result_state
    end
    Expr(:block, next_expr, goto_block, label_defs, result_expr)
end

end # module


text:
module Searching

export binarysearch, BinarySearchTree, BinaryNode, delete!
export size, get, maximum, minimum

# Binary search
struct UnsortedException <: Exception end
Base.showerror(io::IO, e::UnsortedException) = print(io, "input vecotr v is not sorted")

midpoint(lo::Integer, hi::Integer) = lo + div(hi - lo, 2)

function binarysearch(key, v::AbstractVector; algo=:recursive)
    issorted(v) || throw(UnsortedException())
    lo, hi = 1, length(v)
    if algo == :recursive
        binarysearch_rec(key, v, lo, hi)
    elseif algo == :iterative
        binarysearch_iter(key, v, lo, hi)
    end
end

function binarysearch_rec(key, v::AbstractVector, lo::Integer, hi::Integer)
    i, j = lo, hi
    if i <= j
        mid = midpoint(i, j)
        if key < v[mid]
            binarysearch_rec(key, v, i, mid-1)
        elseif key > v[mid]
            binarysearch_rec(key, v, mid+1, j)
        else
            return mid
        end
    else
        return -1
    end
end

function binarysearch_iter(key, v::AbstractVector, lo::Integer, hi::Integer)
    i, j = lo, hi
    while i <= j
        mid = midpoint(i, j)
        if key < v[mid]
            j = mid - 1
        elseif key > v[mid]
            i = mid + 1
        else
            return mid
        end
    end
    return -1
end

# Binary search tree

mutable struct BinaryNode{T}
    val::T
    numnodes::Int
    left::Union{BinaryNode{T}, Nothing}
    right::Union{BinaryNode{T}, Nothing}

    # root constructor
    BinaryNode(val::T) where T = new{T}(val, 1, nothing, nothing)
end

abstract type SearchTree end

"""
Binary search tree.
"""
struct BinarySearchTree{T} <: SearchTree
    root::BinaryNode{T}
end

function BinarySearchTree(v::AbstractVector)
    !isempty(v) || error("empty vector")

    bst = BinarySearchTree(BinaryNode(v[1]))
    for val in v[2:end]
        put!(bst, val)
    end

    return bst
end

Searching.size(::Nothing) = 0
Searching.size(node::BinaryNode{T}) where T = node.numnodes
Searching.size(tree::BinarySearchTree{T}) where T = Searching.size(tree.root)

function updatesize!(node)
    node.numnodes = Searching.size(node.left) + Searching.size(node.right) + 1
    nothing
end

"""
    put!(tree, val)

Insert a value into the binary search tree.

    put!(node, val)

Insert a value following the node.
"""
function put!(node::BinaryNode{T}, val::T) where T
    if val < node.val
        node.left = put!(node.left, val)
    elseif val > node.val
        node.right = put!(node.right, val)
    else
        node.val = val
    end
    updatesize!(node)
    return node
end

put!(tree::BinarySearchTree{T}, val::T) where T = BinarySearchTree(put!(tree.root, val))
put!(::Nothing, val) = BinaryNode(val)

"""
    get(node, val)

Search for the node by value. Return `nothing` if not found.
"""
function get(node::BinaryNode, val)
    if val < node.val
        return get(node.left, val)
    elseif val > node.val
        return get(node.right, val)
    else
        return node
    end
end

get(tree::BinarySearchTree, val) = get(tree.root, val)
get(::Nothing, val)= nothing

"""
    maximum(BinarySearchTree)

Return node with maximum value.
"""
function maximum(node)
    isnothing(node) && return nothing

    while !isnothing(node.right)
        node = node.right
    end

    return node
end

maximum(tree::BinarySearchTree) = maximum(tree.root)

function minimum(node)
    isnothing(node) && return nothing

    while !isnothing(node.left)
        node = node.left
    end

    return node
end

minimum(tree::BinarySearchTree) = minimum(tree.root)

"""
    rank(node, val)

Find the rank for value `val`: number of nodes with values less
than `val` + 1.
"""
function rank(node::BinaryNode, val)
    if val < node.val
        valrank = rank(node.left, val)
    elseif val > node.val
        valrank = rank(node.right, val) + Searching.size(node.left) + 1
    else
        valrank = Searching.size(node.left)
    end

    return valrank
end

rank(::Nothing, val) = 0
rank(tree::BinarySearchTree, val) = rank(tree.root, val) + 1

"""
    select(node, rankval)

Select the node with ranking `rankval`.
We assume rankval is in range between 1 to the size of the tree.
"""
function select(node::BinaryNode, valrank)
    noderank = Searching.size(node.left) + 1
    if valrank < noderank
       node = select(node.left, valrank)
    elseif valrank > noderank
       node = select(node.right, valrank - noderank)
    else
        # found position with `valrank`
        # return node
    end

    return  node
end

select(::Nothing, valrank) = nothing
select(tree::BinarySearchTree, valrank) = select(tree.root, valrank)

"""
    deletemin!(node)
Find the minimum node and delete it from the tree.
We can find the minimum node by move left recursively, until
we find the node with no left child.

"""
function deletemin!(node)
    if isnothing(node.left)
        node = node.right
    else
        node.left = deletemin!(node.left)
    end
    updatesize!(node)
    return node
end

deletemin!(tree::BinarySearchTree) = BinarySearchTree(deletemin!(tree.root))

"""
    copyval!(tonode, fromnode)

Copy values from `fromnode` to `tonode`.
"""
function copyval!(tonode::BinaryNode, fromnode::BinaryNode)
    tonode.val = fromnode.val
    tonode.numnodes = fromnode.numnodes
    nothing
end

"""
    delete!(tree, val)

Delete the node from the binary search tree.
It returns the successor of the deleted node.

    delete!(node, val)

Delete the node with value `val` in the substree
starts with `node`.
"""
function delete!(node, val)
    if isnothing(node)
        return nothing
    elseif val < node.val
        node.left = delete!(node.left, val)
    elseif val > node.val
        node.right = delete!(node.right, val)
    # found the value and node has only one child node
    elseif isnothing(node.right)
        node = node.left
    elseif isnothing(node.left)
        node = node.right
    # node has both left and right child nodes
    else
        # replace current node with its successor
        # here I cheated using copy value (without links),
        # instead of replacing the object
        # This way we keep the original reference to make the
        # side effect consistent: the returned node and original "node"
        # refers to the same object (a substree starting from that node)
        copyval!(node, minimum(node.right))

        # left child remain the same
        # delete the the successor from the right subtree
        # link successor the right child to the pruned subtree
        node.right = deletemin!(node.right)
    end

    updatesize!(node)
    return node
end

delete!(tree::BinarySearchTree, val) = BinarySearchTree(delete!(tree.root, val))

end # module

max_stars_repo_name: Jetafull/AlgorithmsInJulia

max_issues_repo_name: Jetafull/algorithms-in-julia

max_forks_repo_name: Jetafull/algorithms-in-julia


text:
import Base.Dict

export Atoms

"""
`JData`: datatype for storing any data

some data which needs to be updated if the configuration (positions only!) has
changed too much.
"""
mutable struct JData{T}
   max_change::T     # how much X may change before recomputing
   accum_change::T   # how much has it changed already
   data::Any
end

struct LinearConstraint{T}
   C::Matrix{T}
   b::Vector{T}
   desc::String   # description of the constraint
end

mutable struct DofManager{T}
   variablecell::Bool
   xfree::Vector{Int}
   lincons::Vector{LinearConstraint{T}}
   # ----
   X0::Vector{JVec{T}}
   F0::JMat{T}
end

"""
`Atoms{T <: AbstractFloat} <: AbstractAtoms`

The main type storing information about atomic configurations. This is
not normally constructed directly by calling `Atoms`, but using one of the
utility functions such as `bulk`.

## Constructors

TODO

## Convenience functions

* `bulk`

# Getters and Setters

* `positions`, `set_positions!`
* `momenta`, `set_momenta!`
* `masses`, `set_masses!`
* `numbers`, `set_numbers!`
* `cell`, `set_cell!`
* `pbc`, `set_pbc!`
* `calculator`, `set_calculator!`
* `get_data`, `set_data!`
"""
mutable struct Atoms{T <: AbstractFloat} <: AbstractAtoms{T}
   X::Vector{JVec{T}}              # positions
   P::Vector{JVec{T}}              # momenta (or velocities?)
   M::Vector{T}                    # masses
   Z::Vector{Int16}                # atomic numbers
   cell::JMat{T}                   # cell
   pbc::JVec{Bool}                 # boundary condition
   calc::Union{Nothing, AbstractCalculator}
   dofmgr::DofManager
   data::Dict{Any,JData{T}}
end

Atoms(X::Vector{JVec{T}}, P::Vector{JVec{T}}, M::Vector{T},
      Z::Vector{Int16}, cell::JMat{T}, pbc::JVec{Bool},
      calc::Union{Nothing, AbstractCalculator} = nothing) where {T} =
   Atoms(X, P, M, Z, cell, pbc, calc,
         DofManager(length(X), T), Dict{Any,JData{T}}())

Base.eltype(::Atoms{T}) where {T} = T

# derived properties
length(at::Atoms) = length(at.X)

getindex(at::Atoms, i::Integer) = at.X[i]

function setindex!(at::Atoms{T}, val, i::Integer) where {T}
   u = norm(at.X[i] - val)
   at.X[i] = convert(JVec{T}, val)
   update_data!(at, u)
   return val
end

# ------------------------ access to struct fields ----------------------

# ----- getters

positions(at::Atoms) = copy(at.X)
momenta(at::Atoms) = copy(at.P)
masses(at::Atoms) = copy(at.M)
numbers(at::Atoms) = copy(at.Z)
atomic_numbers(at::Atoms) = copy(at.Z)
cell(at::Atoms) = at.cell
pbc(at::Atoms) = at.pbc
calculator(at::Atoms) = at.calc

chemical_symbols(at::Atoms) = Chemistry.chemical_symbol.(at.Z)

# ----- setters

function set_positions!(at::Atoms{T}, X::AbstractVector{JVec{T}})  where T
   update_data!(at, dist(at, X))
   at.X .= X
   return at
end
function set_momenta!(at::Atoms{T}, P::AbstractVector{JVec{T}})  where T
   at.P .= P
   return at
end
function set_masses!(at::Atoms{T}, M::AbstractVector{T})  where T
   update_data!(at, Inf)
   at.M .= M
   return at
end
function set_numbers!(at::Atoms, Z::AbstractVector{<:Integer})
   update_data!(at, Inf)
   at.Z .= Z
   return at
end
function set_cell!(at::Atoms{T}, C::AbstractMatrix) where T
   update_data!(at, Inf)
   at.cell = JMat{T}(C)
   return at
end
function set_pbc!(at::Atoms, p::Union{AbstractVector, Tuple})
   q = JVec{Bool}(p...)
   if at.pbc != q
      update_data!(at, Inf)
      at.pbc = q
   end
   return at
end
function set_calculator!(at::Atoms, calc::Union{AbstractCalculator, Nothing})
   at.calc = calc
   return at
end

# --------------- equality tests -------------------

import Base.==
==(at1::Atoms{T}, at2::Atoms{T}) where T = (
   isapprox(at1, at2, tol = 0.0)  &&
          (at1.data == at2.data)  &&
          (at1.calc == at2.calc)  &&
          (at1.dofmgr == at2.dofmgr) )

import Base.isapprox
function isapprox(at1::Atoms{T}, at2::Atoms{T}; tol = sqrt(eps(T))) where {T}
   if length(at1) != length(at2)
      return false
   end
   if tol > 0
      ndigits = floor(Int, abs(log10(tol)))
      X1 = [round.(x, ndigits) for x in at1.X]
      X2 = [round.(x, ndigits) for x in at2.X]
   else
      X1, X2 = at1.X, at2.X
   end
   p1 = sortperm(X1)
   p2 = sortperm(X2)
   return (maxdist(at1.X[p1], at2.X[p2]) <= tol) &&
          (maxdist(at1.P[p1], at2.P[p2]) <= tol)  &&
          (norm(at1.M[p1] - at2.M[p2], Inf) <= tol) &&
          (at1.Z[p1] == at2.Z[p2]) &&
          (norm(at1.cell - at2.cell, Inf) <= tol) &&
          (at1.pbc == at2.pbc)
end

# --------------- some aliases for easier / direct access --------------
#                 and some access not covered above

function neighbourlist(at::Atoms{T}, rcut::AbstractFloat;
            recompute=false, key="nlist:default",
            int_type = Int32, kwargs...) where {T}
   # we are allowed to use a stored list
   if !recompute
      if has_data(at, key)
         nlist = get_data(at, key)::PairList{T, int_type}
         if cutoff(nlist) >= rcut
            return nlist
         end
      end
   end
   # we are either forces to recompute, or no nlist exists or is not
   # suitable for the request.
   nlist = PairList(positions(at), rcut, cell(at), pbc(at);
                    int_type=int_type, kwargs...)::PairList{T, int_type}
   set_data!(at, key, nlist, 0.0)
   return nlist
end

# function neighbourlist(at::Atoms{T}, rcut::AbstractFloat;
#             recompute=false, key="nlist:default", kwargs...) where {T}
#    return PairList(positions(at), rcut, cell(at), pbc(at); kwargs...)
# end

neighbourlist(at::Atoms) = neighbourlist(at, cutoff(at))

# --------------- data Dict handling ------------------

"""
`has_data(at::Atoms, key)`: checks whether `at` is storing a datum with key `key`; see also `set_data!`
"""
has_data(at::Atoms, key) = haskey(at.data, key)

"""
`get_data(at::Atoms, key)`: retrieves a datum;
see also `set_data!`. Unlike `get_positions`, etc, this returns the original
datum, not a copy (unless it is immutable).
"""
get_data(at::Atoms, key) = (at.data[key]).data

"""
`set_data!(a::Atoms, key, value, max_change=Inf)`:

Stores `value` under an arbitrary key `key` inside `at`. The entry is paired
with an updatemeasure `chg`. Whenever atom positions or the cell change by a
distance `d`, the counter is incremented `chg += d`. When `chg > max_change`,
the entry is deleted from the dictionary.

The two most common uses are with `max_change = Inf` (default) - in this case
the entry will never be deleted; and `max_change = 0.0` - in this case the entry
will be deleted whenever atom positions or the cell change.

The change counter is used, e.g., for neighbourlists with buffer, and for
preconditioners so that the list need not be recomputed each time the
configuration changes.
"""
function set_data!(at::Atoms, key, value, max_change=Inf)
   at.data[key] = JData(max_change, 0.0, value)
   return at
end

"""
`function set_transient!(a::Atoms, key, value, max_change=0.0)`: this is an
alias for `set_data!` but with default `max_change = 0.0` instead of `max_change =
Inf`.
"""
set_transient!(at::Atoms, key, value, max_change=0.0) =
      set_data!(at, key, value, max_change)

"""
`update_data!(at::Atoms, r::Real)`:
increment the change counter in all stored data and delete the
expired ones.
"""
function update_data!(at::Atoms, r::Real)
   for (key, t) in at.data
      t.accum_change += r
      if t.accum_change >= t.max_change
         delete!(at.data, key)
      end
   end
   return at
end

# ------------------------ FIO  ----------------------

Dict(at::Atoms) =
   Dict( "__id__" => "JuLIP_Atoms",
         "X"      => at.X,
         "P"      => at.P,
         "M"      => at.M,
         "Z"      => at.Z,
         "cell"   => at.cell,
         "pbc"    => at.pbc,
         "calc"   => nothing,
         "data"   => nothing )

Atoms(D::Dict) = Atoms(D["X"], D["P"], D["M"], D["Z"], D["cell"], D["pbc"])
          # calc = D["calc"],
          # cons = D["cons"],
          # data = D["data"] )

Base.convert(::Val{:JuLIP_Atoms}, D) = Atoms(D)


text:
# using Pkg
# Pkg.add(PackageSpec(url="https://github.com/legend-exp/LegendHDF5IO.jl.git"))
using HDF5
using LegendHDF5IO: readdata, writedata
using SolidStateDetectors

hd5_sim_filename = "sim.h5"

sim = if isfile(hd5_sim_filename)
    @info "Reading from file $hd5_sim_filename"
    s = HDF5.h5open(hd5_sim_filename, "r") do h5f
        Simulation(readdata(h5f, "Simulation"))
    end
    set_charge_drift_model!(s, ADLChargeDriftModel())
    s
else
    s = Simulation(SSD_examples[:InvertedCoax])
    simulate!(s, max_refinements = 1)
    HDF5.h5open(hd5_sim_filename, "w") do h5f
        writedata(h5f, "Simulation", NamedTuple(s))
    end
end

println(sim)

# Storing and loading individual fields:

sim2 = Simulation(SSD_examples[:InvertedCoax])
calculate_weighting_potential!(sim2, 1) # Weighting potential of channel 1
calculate_weighting_potential!(sim2, 2) # Weighting potential of channel 1

hd5_wps_filename = "weighting_potentials.h5"

HDF5.h5open(hd5_wps_filename, "w") do h5f
    writedata(h5f, "WeightingPotential_001", NamedTuple(sim2.weighting_potentials[1]))
    writedata(h5f, "WeightingPotential_002", NamedTuple(sim2.weighting_potentials[2]))
end

# now read only the first one
sim2 = Simulation(SSD_examples[:InvertedCoax]) # empty Simulation object

sim2.weighting_potentials[1] = HDF5.h5open(hd5_wps_filename, "r") do h5f
    WeightingPotential(readdata(h5f, "WeightingPotential_001"))
end

max_stars_repo_name: rjtayl/NDP

max_issues_repo_name: rjtayl/NDP

max_forks_repo_name: rjtayl/NDP


text:
using svopt
using Test

ϵ = 1e-5

@testset "ratio helpers" begin
    a = -5
    b = 10
    @test svopt.fw1(a, b) - a == b - svopt.fw2(a, b)
end # testset

@testset "optfun" begin
    f1(x) = x^2
    f2(x) = (x-1)^2*10 + 10
    f3(x) = (x+10)^2 - 10
    f4(x) = x^2 - 1
    f5(x) = x^2/3 + 2x - sin(x)

    @test isapprox(optfun(f1, -3, 10), 0, atol=ϵ)
    @test isapprox(optfun(f2, -3, 10), 1, atol=ϵ)
    @test isapprox(optfun(f3, -30, 30), -10, atol=ϵ)
    @test isapprox(optfun(f4, -10, 30), 0, atol=ϵ)
    @test isapprox(optfun(f5, -10, 5), -3.99083, atol=ϵ)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiPolicyV1beta1PodSecurityPolicySpec <: SwaggerModel
    allowPrivilegeEscalation::Any # spec type: Union{ Nothing, Bool } # spec name: allowPrivilegeEscalation
    allowedCSIDrivers::Any # spec type: Union{ Nothing, Vector{IoK8sApiPolicyV1beta1AllowedCSIDriver} } # spec name: allowedCSIDrivers
    allowedCapabilities::Any # spec type: Union{ Nothing, Vector{String} } # spec name: allowedCapabilities
    allowedFlexVolumes::Any # spec type: Union{ Nothing, Vector{IoK8sApiPolicyV1beta1AllowedFlexVolume} } # spec name: allowedFlexVolumes
    allowedHostPaths::Any # spec type: Union{ Nothing, Vector{IoK8sApiPolicyV1beta1AllowedHostPath} } # spec name: allowedHostPaths
    allowedProcMountTypes::Any # spec type: Union{ Nothing, Vector{String} } # spec name: allowedProcMountTypes
    allowedUnsafeSysctls::Any # spec type: Union{ Nothing, Vector{String} } # spec name: allowedUnsafeSysctls
    defaultAddCapabilities::Any # spec type: Union{ Nothing, Vector{String} } # spec name: defaultAddCapabilities
    defaultAllowPrivilegeEscalation::Any # spec type: Union{ Nothing, Bool } # spec name: defaultAllowPrivilegeEscalation
    forbiddenSysctls::Any # spec type: Union{ Nothing, Vector{String} } # spec name: forbiddenSysctls
    fsGroup::Any # spec type: Union{ Nothing, IoK8sApiPolicyV1beta1FSGroupStrategyOptions } # spec name: fsGroup
    hostIPC::Any # spec type: Union{ Nothing, Bool } # spec name: hostIPC
    hostNetwork::Any # spec type: Union{ Nothing, Bool } # spec name: hostNetwork
    hostPID::Any # spec type: Union{ Nothing, Bool } # spec name: hostPID
    hostPorts::Any # spec type: Union{ Nothing, Vector{IoK8sApiPolicyV1beta1HostPortRange} } # spec name: hostPorts
    privileged::Any # spec type: Union{ Nothing, Bool } # spec name: privileged
    readOnlyRootFilesystem::Any # spec type: Union{ Nothing, Bool } # spec name: readOnlyRootFilesystem
    requiredDropCapabilities::Any # spec type: Union{ Nothing, Vector{String} } # spec name: requiredDropCapabilities
    runAsGroup::Any # spec type: Union{ Nothing, IoK8sApiPolicyV1beta1RunAsGroupStrategyOptions } # spec name: runAsGroup
    runAsUser::Any # spec type: Union{ Nothing, IoK8sApiPolicyV1beta1RunAsUserStrategyOptions } # spec name: runAsUser
    seLinux::Any # spec type: Union{ Nothing, IoK8sApiPolicyV1beta1SELinuxStrategyOptions } # spec name: seLinux
    supplementalGroups::Any # spec type: Union{ Nothing, IoK8sApiPolicyV1beta1SupplementalGroupsStrategyOptions } # spec name: supplementalGroups
    volumes::Any # spec type: Union{ Nothing, Vector{String} } # spec name: volumes

    function IoK8sApiPolicyV1beta1PodSecurityPolicySpec(;allowPrivilegeEscalation=nothing, allowedCSIDrivers=nothing, allowedCapabilities=nothing, allowedFlexVolumes=nothing, allowedHostPaths=nothing, allowedProcMountTypes=nothing, allowedUnsafeSysctls=nothing, defaultAddCapabilities=nothing, defaultAllowPrivilegeEscalation=nothing, forbiddenSysctls=nothing, fsGroup=nothing, hostIPC=nothing, hostNetwork=nothing, hostPID=nothing, hostPorts=nothing, privileged=nothing, readOnlyRootFilesystem=nothing, requiredDropCapabilities=nothing, runAsGroup=nothing, runAsUser=nothing, seLinux=nothing, supplementalGroups=nothing, volumes=nothing)
        o = new()
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("allowPrivilegeEscalation"), allowPrivilegeEscalation)
        setfield!(o, Symbol("allowPrivilegeEscalation"), allowPrivilegeEscalation)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("allowedCSIDrivers"), allowedCSIDrivers)
        setfield!(o, Symbol("allowedCSIDrivers"), allowedCSIDrivers)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("allowedCapabilities"), allowedCapabilities)
        setfield!(o, Symbol("allowedCapabilities"), allowedCapabilities)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("allowedFlexVolumes"), allowedFlexVolumes)
        setfield!(o, Symbol("allowedFlexVolumes"), allowedFlexVolumes)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("allowedHostPaths"), allowedHostPaths)
        setfield!(o, Symbol("allowedHostPaths"), allowedHostPaths)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("allowedProcMountTypes"), allowedProcMountTypes)
        setfield!(o, Symbol("allowedProcMountTypes"), allowedProcMountTypes)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("allowedUnsafeSysctls"), allowedUnsafeSysctls)
        setfield!(o, Symbol("allowedUnsafeSysctls"), allowedUnsafeSysctls)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("defaultAddCapabilities"), defaultAddCapabilities)
        setfield!(o, Symbol("defaultAddCapabilities"), defaultAddCapabilities)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("defaultAllowPrivilegeEscalation"), defaultAllowPrivilegeEscalation)
        setfield!(o, Symbol("defaultAllowPrivilegeEscalation"), defaultAllowPrivilegeEscalation)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("forbiddenSysctls"), forbiddenSysctls)
        setfield!(o, Symbol("forbiddenSysctls"), forbiddenSysctls)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("fsGroup"), fsGroup)
        setfield!(o, Symbol("fsGroup"), fsGroup)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("hostIPC"), hostIPC)
        setfield!(o, Symbol("hostIPC"), hostIPC)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("hostNetwork"), hostNetwork)
        setfield!(o, Symbol("hostNetwork"), hostNetwork)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("hostPID"), hostPID)
        setfield!(o, Symbol("hostPID"), hostPID)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("hostPorts"), hostPorts)
        setfield!(o, Symbol("hostPorts"), hostPorts)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("privileged"), privileged)
        setfield!(o, Symbol("privileged"), privileged)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("readOnlyRootFilesystem"), readOnlyRootFilesystem)
        setfield!(o, Symbol("readOnlyRootFilesystem"), readOnlyRootFilesystem)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("requiredDropCapabilities"), requiredDropCapabilities)
        setfield!(o, Symbol("requiredDropCapabilities"), requiredDropCapabilities)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("runAsGroup"), runAsGroup)
        setfield!(o, Symbol("runAsGroup"), runAsGroup)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("runAsUser"), runAsUser)
        setfield!(o, Symbol("runAsUser"), runAsUser)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("seLinux"), seLinux)
        setfield!(o, Symbol("seLinux"), seLinux)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("supplementalGroups"), supplementalGroups)
        setfield!(o, Symbol("supplementalGroups"), supplementalGroups)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicySpec, Symbol("volumes"), volumes)
        setfield!(o, Symbol("volumes"), volumes)
        o
    end
end # type IoK8sApiPolicyV1beta1PodSecurityPolicySpec

const _property_map_IoK8sApiPolicyV1beta1PodSecurityPolicySpec = Dict{Symbol,Symbol}(Symbol("allowPrivilegeEscalation")=>Symbol("allowPrivilegeEscalation"), Symbol("allowedCSIDrivers")=>Symbol("allowedCSIDrivers"), Symbol("allowedCapabilities")=>Symbol("allowedCapabilities"), Symbol("allowedFlexVolumes")=>Symbol("allowedFlexVolumes"), Symbol("allowedHostPaths")=>Symbol("allowedHostPaths"), Symbol("allowedProcMountTypes")=>Symbol("allowedProcMountTypes"), Symbol("allowedUnsafeSysctls")=>Symbol("allowedUnsafeSysctls"), Symbol("defaultAddCapabilities")=>Symbol("defaultAddCapabilities"), Symbol("defaultAllowPrivilegeEscalation")=>Symbol("defaultAllowPrivilegeEscalation"), Symbol("forbiddenSysctls")=>Symbol("forbiddenSysctls"), Symbol("fsGroup")=>Symbol("fsGroup"), Symbol("hostIPC")=>Symbol("hostIPC"), Symbol("hostNetwork")=>Symbol("hostNetwork"), Symbol("hostPID")=>Symbol("hostPID"), Symbol("hostPorts")=>Symbol("hostPorts"), Symbol("privileged")=>Symbol("privileged"), Symbol("readOnlyRootFilesystem")=>Symbol("readOnlyRootFilesystem"), Symbol("requiredDropCapabilities")=>Symbol("requiredDropCapabilities"), Symbol("runAsGroup")=>Symbol("runAsGroup"), Symbol("runAsUser")=>Symbol("runAsUser"), Symbol("seLinux")=>Symbol("seLinux"), Symbol("supplementalGroups")=>Symbol("supplementalGroups"), Symbol("volumes")=>Symbol("volumes"))
const _property_types_IoK8sApiPolicyV1beta1PodSecurityPolicySpec = Dict{Symbol,String}(Symbol("allowPrivilegeEscalation")=>"Bool", Symbol("allowedCSIDrivers")=>"Vector{IoK8sApiPolicyV1beta1AllowedCSIDriver}", Symbol("allowedCapabilities")=>"Vector{String}", Symbol("allowedFlexVolumes")=>"Vector{IoK8sApiPolicyV1beta1AllowedFlexVolume}", Symbol("allowedHostPaths")=>"Vector{IoK8sApiPolicyV1beta1AllowedHostPath}", Symbol("allowedProcMountTypes")=>"Vector{String}", Symbol("allowedUnsafeSysctls")=>"Vector{String}", Symbol("defaultAddCapabilities")=>"Vector{String}", Symbol("defaultAllowPrivilegeEscalation")=>"Bool", Symbol("forbiddenSysctls")=>"Vector{String}", Symbol("fsGroup")=>"IoK8sApiPolicyV1beta1FSGroupStrategyOptions", Symbol("hostIPC")=>"Bool", Symbol("hostNetwork")=>"Bool", Symbol("hostPID")=>"Bool", Symbol("hostPorts")=>"Vector{IoK8sApiPolicyV1beta1HostPortRange}", Symbol("privileged")=>"Bool", Symbol("readOnlyRootFilesystem")=>"Bool", Symbol("requiredDropCapabilities")=>"Vector{String}", Symbol("runAsGroup")=>"IoK8sApiPolicyV1beta1RunAsGroupStrategyOptions", Symbol("runAsUser")=>"IoK8sApiPolicyV1beta1RunAsUserStrategyOptions", Symbol("seLinux")=>"IoK8sApiPolicyV1beta1SELinuxStrategyOptions", Symbol("supplementalGroups")=>"IoK8sApiPolicyV1beta1SupplementalGroupsStrategyOptions", Symbol("volumes")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiPolicyV1beta1PodSecurityPolicySpec }) = collect(keys(_property_map_IoK8sApiPolicyV1beta1PodSecurityPolicySpec))
Swagger.property_type(::Type{ IoK8sApiPolicyV1beta1PodSecurityPolicySpec }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiPolicyV1beta1PodSecurityPolicySpec[name]))}
Swagger.field_name(::Type{ IoK8sApiPolicyV1beta1PodSecurityPolicySpec }, property_name::Symbol) =  _property_map_IoK8sApiPolicyV1beta1PodSecurityPolicySpec[property_name]

function check_required(o::IoK8sApiPolicyV1beta1PodSecurityPolicySpec)
    (getproperty(o, Symbol("fsGroup")) === nothing) && (return false)
    (getproperty(o, Symbol("runAsUser")) === nothing) && (return false)
    (getproperty(o, Symbol("seLinux")) === nothing) && (return false)
    (getproperty(o, Symbol("supplementalGroups")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiPolicyV1beta1PodSecurityPolicySpec }, name::Symbol, val)
end


text:
export _mult, _dot
export _cross, _cross_x, _cross_y, _cross_z
export _sum_cross, _sum_cross_x, _sum_cross_y, _sum_cross_z
export _outer
export slice_inv

"""
local multiplication, avoiding broadcast for AD performance
"""
# scalar-scalar muliplication
function _mult(s₁::AbstractArray{T1,2},s₂::AbstractArray{T2,2}) where {T1,T2}
	@tullio s₃[ix,iy] := s₁[ix,iy] * s₂[ix,iy]
end
function _mult(s₁::AbstractArray{T1,3},s₂::AbstractArray{T2,3}) where {T1,T2}
	@tullio s₃[ix,iy,iz] := s₁[ix,iy,iz] * s₂[ix,iy,iz]
end
# scalar-vector muliplication
function _mult(s::AbstractArray{T1,2},v₁::AbstractArray{T2,3}) where {T1,T2}
	@tullio v₂[a,ix,iy] := s[ix,iy] * v₁[a,ix,iy]
end
function _mult(v₁::AbstractArray{T2,3},s::AbstractArray{T1,2}) where {T1,T2}
	@tullio v₂[a,ix,iy] := s[ix,iy] * v₁[a,ix,iy]
end
function _mult(s::AbstractArray{T1,3},v₁::AbstractArray{T2,4}) where {T1,T2}
	@tullio v₂[a,ix,iy,iz] := s[ix,iy,iz] * v₁[a,ix,iy,iz]
end
function _mult(v₁::AbstractArray{T2,4},s::AbstractArray{T1,3}) where {T1,T2}
	@tullio v₂[a,ix,iy,iz] := s[ix,iy,iz] * v₁[a,ix,iy,iz]
end
# scalar-matrix muliplication
function _mult(s::AbstractArray{T1,2},v₁::AbstractArray{T2,4}) where {T1,T2}
	@tullio v₂[a,b,ix,iy] := s[ix,iy] * v₁[a,b,ix,iy]
end
function _mult(v₁::AbstractArray{T2,4},s::AbstractArray{T1,2}) where {T1,T2}
	@tullio v₂[a,b,ix,iy] := s[ix,iy] * v₁[a,b,ix,iy]
end
function _mult(s::AbstractArray{T1,3},v₁::AbstractArray{T2,5}) where {T1,T2}
	@tullio v₂[a,b,ix,iy,iz] := s[ix,iy,iz] * v₁[a,b,ix,iy,iz]
end
function _mult(v₁::AbstractArray{T2,5},s::AbstractArray{T1,3}) where {T1,T2}
	@tullio v₂[a,b,ix,iy,iz] := s[ix,iy,iz] * v₁[a,b,ix,iy,iz]
end

"""
dot products in the first few dimensions
"""
# first-order (linear) tensor-vector muliplication
function _dot(χ::AbstractArray{T,2},v₁::AbstractVector) where T<:Real
	@tullio v₂[i] := χ[i,j] * v₁[j]
end

function _dot(χ::AbstractArray{T,4},v₁::AbstractArray{Complex{T},3}) where T<:Real
	@tullio v₂[i,ix,iy] := χ[i,j,ix,iy] * v₁[j,ix,iy]
end

function _dot(χ::AbstractArray{T,5},v₁::AbstractArray{Complex{T},4}) where T<:Real
	@tullio v₂[i,ix,iy,iz] := χ[i,j,k,ix,iy,iz] * v₁[j,ix,iy,iz]
end

# first-order (linear) tensor-tensor muliplication

function _dot(A::AbstractArray{T,4},B::AbstractArray{T,4}) where T<:Real
	@tullio C[a,b,ix,iy] := A[a,j,ix,iy] * B[j,b,ix,iy]
end

function _dot(A::AbstractArray{T,4},B::AbstractArray{T,4},C::AbstractArray{T,4}) where T<:Real
	@tullio D[a,b,ix,iy] := A[a,j,ix,iy] * B[j,k,ix,iy] * C[k,b,ix,iy]
end

function _dot(A::AbstractArray{T,5},B::AbstractArray{T,5}) where T<:Real
	@tullio C[a,b,ix,iy,iz] := A[a,j,ix,iy,iz] * B[j,b,ix,iy,iz]
end

function _dot(A::AbstractArray{T,5},B::AbstractArray{T,5},C::AbstractArray{T,5}) where T<:Real
	@tullio D[a,b,ix,iy,iz] := A[a,j,ix,iy,iz] * B[j,k,ix,iy,iz] * C[k,b,ix,iy,iz]
end

# second-order nonlinear tensor muliplication
function _dot(χ::AbstractArray{T,3},v₁::AbstractVector,v₂::AbstractVector) where T<:Real
	@tullio v₃[i] := χ[i,j,k] * v₁[j] * v₂[k]
end

function _dot(χ::AbstractArray{T,5},v₁::AbstractArray{Complex{T},3},v₂::AbstractArray{Complex{T},3}) where T<:Real
	@tullio v₃[i,ix,iy] := χ[i,j,k,ix,iy] * v₁[j,ix,iy] * v₂[k,ix,iy]
end

function _dot(χ::AbstractArray{T,6},v₁::AbstractArray{Complex{T},4},v₂::AbstractArray{Complex{T},4}) where T<:Real
	@tullio v₃[i,ix,iy,iz] := χ[i,j,k,ix,iy,iz] * v₁[j,ix,iy,iz] * v₂[k,ix,iy,iz]
end

# third-order nonlinear tensor muliplication
function _dot(χ::AbstractArray{T,4},v₁::AbstractVector,v₂::AbstractVector,v₃::AbstractVector) where T<:Real
	@tullio v₄[i] := χ[i,j,k,l] * v₁[j] * v₂[k] * v₃[l]
end

"""
Tensor utilities
"""

function slice_inv(A::Array{T,4}) where {T}
    B = reinterpret(SArray{Tuple{3,3},T,2,9}, vec(A))
    C = map(inv, B)
    reshape(reinterpret(T, B), size(A))
end

function slice_inv(A::Array{T,5}) where {T}
    B = reinterpret(SArray{Tuple{3,3},T,2,9}, vec(A))
    C = map(inv, B)
    reshape(reinterpret(T, B), size(A))
end

"""
Spatially local linear algebra operations for flat (non-nested) arrays of 3-vectors and 3x3 tensors
"""
function _cross_x(v1::TA1,v2::TA2) where {TA1<:AbstractArray{<:Number,3},TA2<:AbstractArray{<:Number,3}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio v3_x[ix,iy] := v1[2,ix,iy] * v2[3,ix,iy] - v1[3,ix,iy] * v2[2,ix,iy]
end
function _cross_y(v1::TA1,v2::TA2) where {TA1<:AbstractArray{<:Number,3},TA2<:AbstractArray{<:Number,3}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio v3_y[ix,iy] := v1[1,ix,iy] * v2[3,ix,iy] - v1[3,ix,iy] * v2[1,ix,iy]
end
function _cross_z(v1::TA1,v2::TA2) where {TA1<:AbstractArray{<:Number,3},TA2<:AbstractArray{<:Number,3}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio v3_z[ix,iy] := v1[1,ix,iy] * v2[2,ix,iy] - v1[2,ix,iy] * v2[1,ix,iy]
end
function _cross(v1::TA1,v2::TA2) where {TA1<:AbstractArray{<:Number,3},TA2<:AbstractArray{<:Number,3}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio v3[i,ix,iy] := v1[mod(i-2),ix,iy] * v2[mod(i-1),ix,iy] - v1[mod(i-1),ix,iy] * v2[mod(i-2),ix,iy] (i in 1:3)
end
function _cross(v1::TA,v2::TV) where {TA<:AbstractArray{<:Number,3},TV<:AbstractVector{<:Number}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio v3[i,ix,iy] := v1[mod(i-2),ix,iy] * v2[mod(i-1)] - v1[mod(i-1),ix,iy] * v2[mod(i-2)] (i in 1:3)
end
function _cross(v1::TV,v2::TA) where {TA<:AbstractArray{<:Number,3},TV<:AbstractVector{<:Number}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio v3[i,ix,iy] := v1[mod(i-2)] * v2[mod(i-1),ix,iy] - v1[mod(i-1)] * v2[mod(i-2),ix,iy] (i in 1:3)
end
function _sum_cross_x(v1::TA1,v2::TA2) where {TA1<:AbstractArray{<:Number,3},TA2<:AbstractArray{<:Number,3}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio sum_v3_x := v1[2,ix,iy] * v2[3,ix,iy] - v1[3,ix,iy] * v2[2,ix,iy]
end
function _sum_cross_y(v1::TA1,v2::TA2) where {TA1<:AbstractArray{<:Number,3},TA2<:AbstractArray{<:Number,3}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio sum_v3_y := v1[1,ix,iy] * v2[3,ix,iy] - v1[3,ix,iy] * v2[1,ix,iy]
end
function _sum_cross_z(v1::TA1,v2::TA2) where {TA1<:AbstractArray{<:Number,3},TA2<:AbstractArray{<:Number,3}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio sum_v3_z := v1[1,ix,iy] * v2[2,ix,iy] - v1[2,ix,iy] * v2[1,ix,iy]
end
function _sum_cross(v1::TA1,v2::TA2) where {TA1<:AbstractArray{<:Number,3},TA2<:AbstractArray{<:Number,3}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio sum_v3[i] := v1[mod(i-2),ix,iy] * v2[mod(i-1),ix,iy] - v1[mod(i-1),ix,iy] * v2[mod(i-2),ix,iy] (i in 1:3)
end
function _cross_x(v1::TA1,v2::TA2) where {TA1<:AbstractArray{<:Number,4},TA2<:AbstractArray{<:Number,4}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio v3_x[ix,iy,iz] := v1[2,ix,iy,iz] * v2[3,ix,iy,iz] - v1[3,ix,iy,iz] * v2[2,ix,iy,iz]
end
function _cross_y(v1::TA1,v2::TA2) where {TA1<:AbstractArray{<:Number,4},TA2<:AbstractArray{<:Number,4}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio v3_y[ix,iy,iz] := v1[1,ix,iy,iz] * v2[3,ix,iy,iz] - v1[3,ix,iy,iz] * v2[1,ix,iy,iz]
end
function _cross_z(v1::TA1,v2::TA2) where {TA1<:AbstractArray{<:Number,4},TA2<:AbstractArray{<:Number,4}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio v3_z[ix,iy,iz] := v1[1,ix,iy,iz] * v2[2,ix,iy,iz] - v1[2,ix,iy,iz] * v2[1,ix,iy,iz]
end
function _cross(v1::TA1,v2::TA2) where {TA1<:AbstractArray{<:Number,4},TA2<:AbstractArray{<:Number,4}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio v3[i,ix,iy,iz] := v1[mod(i-2),ix,iy,iz] * v2[mod(i-1),ix,iy,iz] - v1[mod(i-1),ix,iy,iz] * v2[mod(i-2),ix,iy,iz] (i in 1:3)
end
function _cross(v1::TA,v2::TV) where {TA<:AbstractArray{<:Number,4},TV<:AbstractVector{<:Number}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio v3[i,ix,iy,iz] := v1[mod(i-2),ix,iy,iz] * v2[mod(i-1)] - v1[mod(i-1),ix,iy,iz] * v2[mod(i-2)] (i in 1:3)
end
function _cross(v1::TV,v2::TA) where {TA<:AbstractArray{<:Number,4},TV<:AbstractVector{<:Number}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio v3[i,ix,iy,iz] := v1[mod(i-2)] * v2[mod(i-1),ix,iy,iz] - v1[mod(i-1)] * v2[mod(i-2),ix,iy,iz] (i in 1:3)
end
function _sum_cross_x(v1::TA1,v2::TA2) where {TA1<:AbstractArray{<:Number,4},TA2<:AbstractArray{<:Number,4}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio sum_v3_x := v1[2,ix,iy,iz] * v2[3,ix,iy,iz] - v1[3,ix,iy,iz] * v2[2,ix,iy,iz]
end
function _sum_cross_y(v1::TA1,v2::TA2) where {TA1<:AbstractArray{<:Number,4},TA2<:AbstractArray{<:Number,4}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio sum_v3_y := v1[1,ix,iy,iz] * v2[3,ix,iy,iz] - v1[3,ix,iy,iz] * v2[1,ix,iy,iz]
end
function _sum_cross_z(v1::TA1,v2::TA2) where {TA1<:AbstractArray{<:Number,4},TA2<:AbstractArray{<:Number,4}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio sum_v3_z := v1[1,ix,iy,iz] * v2[2,ix,iy,iz] - v1[2,ix,iy,iz] * v2[1,ix,iy,iz]
end
function _sum_cross(v1::TA1,v2::TA2) where {TA1<:AbstractArray{<:Number,4},TA2<:AbstractArray{<:Number,4}}
        @assert size(v1,1)==3
        @assert size(v2,1)==3
        @tullio sum_v3[i] := v1[mod(i-2),ix,iy,iz] * v2[mod(i-1),ix,iy,iz] - v1[mod(i-1),ix,iy,iz] * v2[mod(i-2),ix,iy,iz] (i in 1:3)
end

function rrule(::typeof(_cross),v₁,v₂)
        v₃ = _cross(v₁,v₂)
        function _cross_pullback(v̄₃)
                return NoTangent(), @thunk(_cross(conj(v₂),v̄₃)), @thunk(_cross(v̄₃,conj(v₁)))
        end
        return v₃, _cross_pullback
end

function rrule(::typeof(_sum_cross),v₁::AbstractArray{<:Number,N},v₂::AbstractArray{<:Number,N}) where N
        sum_v₃ = _sum_cross(v₁,v₂)
        function _sum_cross_pullback(sum_v̄₃)
                return NoTangent(), @thunk(_cross(conj(v₂),repeat(sum_v̄₃,outer=(1,size(v₂)[2:N]...)))) , @thunk(_cross(repeat(sum_v̄₃,outer=(1,size(v₁)[2:N]...)),conj(v₁)))
        end
        return sum_v₃, _sum_cross_pullback
end

function rrule(::typeof(_sum_cross_x),v₁::AbstractArray{<:Number,N},v₂::AbstractArray{<:Number,N}) where N
        sum_v₃_x = _sum_cross_x(v₁,v₂)
        function _sum_cross_x_pullback(sum_v̄₃_x)
                return NoTangent(), @thunk(_cross(conj(v₂),repeat([sum_v̄₃_x,0.,0.],outer=(1,size(v₂)[2:N]...)))) , @thunk(_cross(repeat([sum_v̄₃_x,0.,0.],outer=(1,size(v₁)[2:N]...)),conj(v₁)))
        end
        return sum_v₃_x, _sum_cross_x_pullback
end

function rrule(::typeof(_sum_cross_y),v₁::AbstractArray{<:Number,N},v₂::AbstractArray{<:Number,N}) where N
        sum_v₃_y = _sum_cross_y(v₁,v₂)
        function _sum_cross_y_pullback(sum_v̄₃_y)
                return NoTangent(), @thunk(_cross(conj(v₂),repeat([0.,sum_v̄₃_y,0.],outer=(1,size(v₂)[2:N]...)))) , @thunk(_cross(repeat([0.,sum_v̄₃_y,0.],outer=(1,size(v₁)[2:N]...)),conj(v₁)))
        end
        return sum_v₃_y, _sum_cross_y_pullback
end

function rrule(::typeof(_sum_cross_z),v₁::AbstractArray{<:Number,N},v₂::AbstractArray{<:Number,N}) where N
        sum_v₃_z = _sum_cross_z(v₁,v₂)
        function _sum_cross_z_pullback(sum_v̄₃_z)
                return NoTangent(), @thunk(_cross(conj(v₂),repeat([0.,0.,sum_v̄₃_z],outer=(1,size(v₂)[2:N]...)))) , @thunk(_cross(repeat([0.,0.,sum_v̄₃_z],outer=(1,size(v₁)[2:N]...)),conj(v₁)))
        end
        return sum_v₃_z, _sum_cross_z_pullback
end

# Outer products
function _outer(v1::TA1,v2::TA2) where {TA1<:AbstractArray{<:Number,3},TA2<:AbstractArray{<:Number,3}}
		@tullio A[i,j,ix,iy] := v1[i,ix,iy] * conj(v2)[j,ix,iy]
		# v2_star = conj(v2)
		# @tullio A[i,j,ix,iy] := v1[i,ix,iy] * v2_star[j,ix,iy]
end

function _outer(v1::TA1,v2::TA2) where {TA1<:AbstractArray{<:Number,4},TA2<:AbstractArray{<:Number,4}}
        @tullio A[i,j,ix,iy,iz] := v1[i,ix,iy,iz] * conj(v2)[j,ix,iy,iz]
		# v2_star = conj(v2)
		# @tullio A[i,j,ix,iy,iz] := v1[i,ix,iy,iz] * v2_star[j,ix,iy,iz]
end


text:
using ADCME
using LinearAlgebra
using MAT
using AdFem
using PyPlot
using SparseArrays

SMALL_SIZE = 22
MEDIUM_SIZE = 24
BIGGER_SIZE = 24

plt.rc("font", size=SMALL_SIZE)          # controls default text sizes
plt.rc("axes", titlesize=SMALL_SIZE)     # fontsize of the axes title
plt.rc("axes", labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
plt.rc("xtick", labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc("ytick", labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc("legend", fontsize=SMALL_SIZE)    # legend fontsize
plt.rc("figure", titlesize=BIGGER_SIZE)  # fontsize of the figure title

# Example 1:

function u_exact(x,y)
    if y == 0.0
        return 1.0
    else
        return 0.0
    end
end

function v_exact(x,y)
    0.0
end

function p_exact(x,y)
    0.0
end

function ffunc(x, y)
    0.0
end
function gfunc(x, y)
    0.0
end

function hfunc(x,y)
    0.0
end

function nu_exact(x, y)
    1 + 6 * x^2 + x / (1 + 2 * y^2)
end

function nu_nn(x, y)
    1.5 + fc([x y], [20,20,20,1])|>squeeze
end

m = 20
n = 20
h = 1/n
# nu = 0.01

F1 = compute_fem_source_term1(eval_f_on_gauss_pts(ffunc, m, n, h), m, n, h)
F2 = compute_fem_source_term1(eval_f_on_gauss_pts(gfunc, m, n, h), m, n, h)
H = h^2*eval_f_on_fvm_pts(hfunc, m, n, h)
B = constant(compute_interaction_matrix(m, n, h))

# compute F
xy = fem_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
u0 = @. u_exact(x,y)
v0 = @. v_exact(x,y)

xy = fvm_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
p0 = @. p_exact(x,y)
# nu_gauss = Variable(ones(4*m*n))
# nu_fvm = Variable(ones(m*n))
nu_fvm = matread("sess_figures/ex1/ones50.mat")["Variablecolon0"]
nu_fvm = constant(nu_fvm)
nu_gauss = reshape(repeat(nu_fvm, 1, 4), (-1,))

Laplace = compute_fem_laplace_matrix1(nu_gauss, m, n, h)

# Laplace = constant(nu*compute_fem_laplace_matrix1(m, n, h))
function compute_residual(S)
    u, v, p = S[1:(m+1)*(n+1)], S[(m+1)*(n+1)+1:2(m+1)*(n+1)], S[2(m+1)*(n+1)+1:2(m+1)*(n+1)+m*n]
    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    interaction = compute_interaction_term(p, m, n, h) # julia kernel needed
    f1 = compute_fem_source_term1(ugauss.*ux, m, n, h)
    f2 = compute_fem_source_term1(vgauss.*uy, m, n, h)
    f3 = -interaction[1:(m+1)*(n+1)]
    f4 = Laplace*u
    f5 = -F1
    F = f1 + f2 + f3 + f4 + f5

    g1 = compute_fem_source_term1(ugauss.*vx, m, n, h)
    g2 = compute_fem_source_term1(vgauss.*vy, m, n, h)
    g3 = -interaction[(m+1)*(n+1)+1:end]
    g4 = Laplace*v
    g5 = -F2
    G = g1 + g2 + g3 + g4 + g5

    H0 = -B * [u;v] + H

    R = [F;G;H0]
    return R
end

function compute_jacobian(S)
    u, v, p = S[1:(m+1)*(n+1)], S[(m+1)*(n+1)+1:2(m+1)*(n+1)], S[2(m+1)*(n+1)+1:2(m+1)*(n+1)+m*n]
    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    M1 = constant(compute_fem_mass_matrix1(ux, m, n, h))
    M2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h)) # a julia kernel needed
    M3 = Laplace
    Fu = M1 + M2 + M3

    Fv = constant(compute_fem_mass_matrix1(uy, m, n, h))

    N1 = constant(compute_fem_mass_matrix1(vy, m, n, h))
    N2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h))
    N3 = Laplace
    Gv = N1 + N2 + N3

    Gu = constant(compute_fem_mass_matrix1(vx, m, n, h))

    J0 = [Fu Fv
          Gu Gv]
    J = [J0 -B'
        -B spdiag(zeros(size(B,1)))]
end

NT = 5
# S = zeros(m*n+2(m+1)*(n+1), NT+1)
# S[1:m+1, 1] = ones(m+1,1)

bd = bcnode("all", m, n, h)
# bd = [bd; bd .+ (m+1)*(n+1); ((1:m) .+ 2(m+1)*(n+1))]
bd = [bd; bd .+ (m+1)*(n+1); 2*(m+1)*(n+1)+(n-1)*m+1:2*(m+1)*(n+1)+(n-1)*m+2] # only apply Dirichlet to velocity

function solve_steady_cavityflow_one_step(S)
    residual = compute_residual(S)
    J = compute_jacobian(S)

    J, _ = fem_impose_Dirichlet_boundary_condition1(J, bd, m, n, h)
    # residual[bd] .= 0.0
    residual = scatter_update(residual, bd, zeros(length(bd)))

    d = J\residual
    S_new = S - d
    return S_new
end

function condition(i, S_arr)
    i <= NT + 1
end

function body(i, S_arr)
    S = read(S_arr, i-1)
    op = tf.print("i=",i)
    i = bind(i, op)
    S_new = solve_steady_cavityflow_one_step(S)
    S_arr = write(S_arr, i, S_new)
    return i+1, S_arr
end

S_arr = TensorArray(NT+1)
S_arr = write(S_arr, 1, [u0; v0; p0])

i = constant(2, dtype=Int32)

_, S = while_loop(condition, body, [i, S_arr])
S = set_shape(stack(S), (NT+1, 2*(m+1)*(n+1)+m*n))

S_true = matread("sess_figures/ex1/steady_cavity_data.mat")["V"]
loss = mean((S[end,1:2*(m+1)*(n+1)] - S_true[end,1:2*(m+1)*(n+1)])^2)
loss = loss * 1e10

# vmax = maximum(nu_exact.(x,y))
# vmin = minimum(nu_exact.(x,y))

sess = Session(); init(sess)
# ADCME.load(sess, "sess_figures/ex1/ones50.mat")
nu_fvm = matread("sess_figures/ex1/ones50.mat")["Variablecolon0"]

vmax = 6.8
vmin = -0.5

figure();visualize_scalar_on_fvm_points(nu_exact.(x,y), m, n, h, vmax=vmax, vmin=vmin); tight_layout()
savefig("sess_figures/ex1/ones/viscosity_exact.png")

# figure();visualize_scalar_on_gauss_points(run(sess, nu_gauss), m, n, h); #gca().invert_yaxis()
# savefig("sess_figures/ex1/ones/viscosity_prediction.png")

figure();visualize_scalar_on_gauss_points(run(sess, nu_gauss), m, n, h, vmax=vmax, vmin=vmin); tight_layout() ;#gca().invert_yaxis()
savefig("sess_figures/ex1/ones/viscosity_prediction.png")

# figure();visualize_scalar_on_fvm_points(nu_exact.(x,y).-run(sess, nu_gauss)[1:4:end], m, n, h);
# savefig("sess_figures/ex1/ones/viscosity_error.png")

figure();visualize_scalar_on_fvm_points(nu_exact.(x,y).-nu_fvm, m, n, h); tight_layout()
savefig("sess_figures/ex1/ones/viscosity_error.png")

vmax = 41
vmin = -39
figure();visualize_scalar_on_fvm_points(S_true[end, 2*(m+1)*(n+1)+1:end], m, n, h, vmax=vmax, vmin=vmin); tight_layout()
savefig("sess_figures/ex1/ones/pressure_data.png")

figure();visualize_scalar_on_fvm_points(run(sess, S[end,2*(m+1)*(n+1)+1:end]), m, n, h, vmax=vmax, vmin=vmin); tight_layout()
savefig("sess_figures/ex1/ones/pressure_prediction.png")

figure();visualize_scalar_on_fvm_points(S_true[end, 2*(m+1)*(n+1)+1:end] - run(sess, S[end,2*(m+1)*(n+1)+1:end]), m, n, h); tight_layout()
savefig("sess_figures/ex1/ones/pressure_error.png")

figure();visualize_scalar_on_fvm_points_nocontour(S_true[end, 2*(m+1)*(n+1)+1:end] - run(sess, S[end,2*(m+1)*(n+1)+1:end]), m, n, h); tight_layout()
savefig("sess_figures/ex1/ones/pressure_error_nocontour.png")

function visualize_scalar_on_fvm_points_nocontour(φ::Array{Float64, 1}, m::Int64, n::Int64, h::Float64)
    φ = reshape(φ, m, n)'|>Array
    m_ = mean(φ)
    s = std(φ)
    vmin, vmax = m_ - 2s, m_ + 2s
    x = (1:m)*h .- 0.5h
    y = (1:n)*h .- 0.5h
    ln = pcolormesh(x, y, φ, vmin= vmin, vmax=vmax)
    colorbar()
    # c = contour(φ[1,:,:], 10, cmap="jet", vmin=vmin,vmax=vmax)
    axis("scaled")
    xlabel("x")
    ylabel("y")
    ln = gca().pcolormesh(x, y, φ, vmin= vmin, vmax=vmax)
    # c = gca().contour(x, y, φ, 10, cmap="jet", vmin=vmin,vmax=vmax)
    gca().invert_yaxis()
end


text:
###############
# H-OBCA: Hierarchical Optimization-based Collision Avoidance - a path planner for autonomous parking
# Copyright (C) 2018
# Alexander LINIGER [liniger@control.ee.ethz.ch; Automatic Control Lab, ETH Zurich]
# Xiaojing ZHANG [xiaojing.zhang@berkeley.edu; MPC Lab, UC Berkeley]
# Atsushi SAKAI [atsushisakai@global.komatsu; Komatsu Ltd / MPC Lab]
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
###############
# The paper describing the theory can be found here:
# 	X. Zhang, A. Liniger and F. Borrelli; "Optimization-Based Collision Avoidance"; Technical Report, 2017, [https://arxiv.org/abs/1711.03449]
#   X. Zhang, A. Liniger, A. Sakai and F. Borrelli; "Autonomous  Parking  using  Optimization-Based  Collision  Avoidance"; Technical Report, 2018 [add URL]
###############

###############
# setup.jl: run this file before running main.jl
###############

##############################
# include JuMP -> Optimization modeling tool
# include IPOPT -> Interior point based nonlinear solver
# include PyPlot -> Ploting library (matplotlib python)
##############################
using PyPlot, NearestNeighbors, ControlSystems, PyCall

include("ParkingConstraints.jl")
include("ParkingSignedDist.jl")
include("obstHrep.jl")

include("plotTraj.jl")
include("hybrid_a_star.jl")

include("DualMultWS.jl")
include("veloSmooth.jl")

@pyimport scipy.interpolate as intp

# function that clears terminal output
clear() = run(@static is_unix() ? `clear` : `cmd /c cls`)

max_stars_repo_name: treyfortmuller/barc

max_issues_repo_name: treyfortmuller/barc