text:
include("rrc_types.jl")

function getFileContent(fName::String)
  char_pointer = ccall(dlsym(rrlib, :getFileContent), cdecl, Ptr{UInt8}, (Ptr{UInt8},), fName)
  julia_str = unsafe_string(char_pointer)
  freeText(char_pointer)
  return julia_str
end

function createText(text::String)
  char_pointer = ccall(dlsym(rrlib, :createText), cdecl, Ptr{UInt8}, (Ptr{UInt8},), text)
  julia_str = unsafe_string(char_pointer)
  freeText(char_pointer)
  return julia_str
end

function createTextMemory(count::Int64)
  char_pointer = ccall(dlsym(rrlib, :createTextMemory), cdecl, Ptr{UInt8}, (Int64,), count)
  julia_str = unsafe_string(char_pointer)
  freeText(char_pointer)
  return julia_str
end

###############################################################################
#                            List Handling Routines                           #
###############################################################################

function createRRList()
  return ccall(dlsym(rrlib, :createRRList), cdecl, Ptr{RRList}, ())
end

function freeRRList(theList)
  ccall(dlsym(rrlib, :freeRRList), cdecl, Nothing, (Ptr{RRList},), theList)
end

function getListLength(myList)
  return ccall(dlsym(rrlib, :getListLength), cdecl, Int64, (Ptr{RRList},), myList)
end

function createIntegerItem(value::Int64)
  return ccall(dlsym(rrlib, :createIntegerItem), cdecl, Ptr{RRListItem}, (Int64,), value)
end

function createDoubleItem(value::Float64)
  return ccall(dlsym(rrlib, :createDoubleItem), cdecl, Ptr{RRListItem}, (Float64,), value)
end

function createStringItem(value::String)
  return ccall(dlsym(rrlib, :createStringItem), cdecl, Ptr{RRListItem}, (Ptr{UInt8},), value)
end

function createListItem(value::Ptr{RRList})
  return ccall(dlsym(rrlib, :createListItem), cdecl, Ptr{RRListItem}, (Ptr{RRList},), value)
end

function addItem(list::Ptr{RRList}, item::Ptr{RRListItem})
  return ccall(dlsym(rrlib, :addItem), cdecl, Int64, (Ptr{RRList}, Ptr{RRListItem}), list, item)
end

function getListItem(list::Ptr{RRList}, index::Int64)
  return ccall(dlsym(rrlib, :getListItem), cdecl, Ptr{RRListItem}, (Ptr{RRList}, Int64), list, index)
end

function isListItemInteger(item::Ptr{RRListItem})
  return ccall(dlsym(rrlib, :isListItemInteger ), cdecl, Bool, (Ptr{RRListItem},), item)
end

function isListItemDouble(item::Ptr{RRListItem})
  return ccall(dlsym(rrlib, :isListItemDouble ), cdecl, Bool, (Ptr{RRListItem},), item)
end

function isListItemString(item::Ptr{RRListItem})
  return ccall(dlsym(rrlib, :isListItemString), cdecl, Bool, (Ptr{RRListItem},), item)
end

function isListItemList(item::Ptr{RRListItem})
  return ccall(dlsym(rrlib, :isListItemList), cdecl, Bool, (Ptr{RRListItem},), item)
end

function isListItem(item::Ptr{RRListItem}, itemType)
  return ccall(dlsym(rrlib, :isListItem), cdecl, Bool, (Ptr{RRListItem}, Ptr{Nothing}), item, itemType)
end

function getIntegerListItem(item::Ptr{RRListItem})
  value = Array{Int64}(undef,1)
  status = ccall(dlsym(rrlib, :getIntegerListItem), cdecl, Bool, (Ptr{RRListItem}, Ptr{Int64}), item, value)
  if status == false
    error(getLastError())
  end
  return value[1]
end

function getDoubleListItem(item::Ptr{RRListItem})
  value = Array{Float64}(undef,1)
  status = ccall(dlsym(rrlib, :getDoubleListItem ), cdecl, Bool, (Ptr{RRListItem}, Ptr{Float64}), item, value)
  if status == false
    error(getLastError())
  end
  return value[1]
end

function getStringListItem(item::Ptr{RRListItem})
  return bytestring(ccall(dlsym(rrlib, :getStringListItem), cdecl, Ptr{UInt8}, (Ptr{RRListItem},), item))
end

function getList(item::Ptr{RRListItem})
  return ccall(dlsym(rrlib, :getList), cdecl, Ptr{RRList}, (Ptr{RRListItem},), item)
end

###############################################################################
#                              Helper Routines                                #
###############################################################################

function getVectorLength(vector::Ptr{RRVector})
  return ccall(dlsym(rrlib, :getVectorLength), cdecl, Int64, (Ptr{RRVector},), vector)
end

function createVector(size::Int64)
  return ccall(dlsym(rrlib, :createVector), cdecl, Ptr{RRVector}, (Int64,), size)
end

function getVectorElement(vector::Ptr{RRVector}, index::Int64)
  value = Array{Float64}(undef,1)
  status = ccall(dlsym(rrlib, :getVectorElement), cdecl, Bool, (Ptr{RRVector}, Int64, Ptr{Float64}), vector, index, value)
  if status == false
    error(getLastError())
  end
  return value[1]
end

function setVectorElement(vector::Ptr{RRVector}, index::Int64, value::Float64)
  status = ccall(dlsym(rrlib, :setVectorElement), cdecl, Bool, (Ptr{RRVector}, Int64, Float64), vector, index, value)
  if status == false
    error(getLastError())
  end
end

function createRRMatrix(r::Int64, c::Int64)
  return ccall(dlsym(rrlib, :createRRMatrix), cdecl, Ptr{RRDoubleMatrix}, (Int64, Int64), r, c)
end

function getMatrixNumRows(m::Ptr{RRDoubleMatrix})
  return ccall(dlsym(rrlib, :getMatrixNumRows), cdecl, Int64, (Ptr{RRDoubleMatrix},), m)
end

function getMatrixNumCols(m::Ptr{RRDoubleMatrix})
  return ccall(dlsym(rrlib, :getMatrixNumCols), cdecl, Int64, (Ptr{RRDoubleMatrix},), m)
end

function getMatrixElement(m::Ptr{RRDoubleMatrix}, r::Int64, c::Int64)
  value = Array{Float64}(undef,1)
  status = ccall(dlsym(rrlib, :getMatrixElement), cdecl, Bool, (Ptr{RRDoubleMatrix}, Int64, Int64, Ptr{Cdouble}), m, r, c, value)
  if status == false
    error(getLastError())
  end
  return value[1]
end

function setMatrixElement(m::Ptr{RRDoubleMatrix}, r::Int64, c::Int64, value::Float64)
  status = ccall(dlsym(rrlib, :setMatrixElement), cdecl, Bool, (Ptr{RRDoubleMatrix}, Int64, Int64, Float64), m, r, c, value)
  if status == false
    error(getLastError())
  end
end

function getComplexMatrixElement(m::Ptr{RRComplexMatrix}, r::Int64, c::Int64)
  value = Ptr{RRComplex}
  status = ccall(dlsym(rrlib, :getComplexMatrixElement), cdecl, Bool, (Ptr{RRComplexMatrix}, Int64, Int64, Ptr{RRComplex}), m, r, c, value)
  if status == false
    error(getLastError())
  end
  return value
end

function setComplexMatrixElement(m::Ptr{RRComplexMatrix}, r::Int64, c::Int64, value::Float64)
  status = ccall(dlsym(rrlib, :setComplexMatrixElement), cdecl, Bool, (Ptr{RRComplexMatrix}, Int64, Int64, Ptr{RRComplex}), m, r, c, value)
  if status == false
    error(getLastError())
  end
end

function getRRDataNumRows(rrData::Ptr{RRCData})
  return ccall(dlsym(rrlib, :getRRDataNumRows), cdecl, Int64, (Ptr{RRCData},), rrData)
end

function getRRDataNumCols(rrData::Ptr{RRCData})
  return ccall(dlsym(rrlib, :getRRDataNumCols), cdecl, Int64, (Ptr{RRCData},), rrData)
end

function getRRCDataElement(rrData::Ptr{RRCData}, r::Int64, c::Int64)
  value = Array{Float64}(undef,1)
  status = ccall(dlsym(rrlib, :getRRCDataElement), cdecl, Bool, (Ptr{RRCData}, Int64, Int64, Ptr{Float64}), rrData, r, c, value)
  if status == false
    error(getLastError())
  end
  return value[1]
end

function getRRDataColumnLabel(rrData::Ptr{RRCData}, column::Int64)
  char_pointer = ccall(dlsym(rrlib, :getRRDataColumnLabel), cdecl, Ptr{UInt8}, (Ptr{RRCData}, Int64), rrData, column)
  julia_str = unsafe_string(char_pointer)
  freeText(char_pointer)
  return julia_str
end

function writeRRData(rr::Ptr{Nothing}, fileNameAndPath::String)
  status = ccall(dlsym(rrlib, :writeRRData), cdecl, Bool, (Ptr{Nothing}, Ptr{UInt8}), rr, fileNameAndPath)
  if status == false
    error(getLastError())
  end
end

function compileSource(rr::Ptr{Nothing}, fName::String)
  status = ccall(dlsym(rrlib, :compileSource), cdecl, Bool, (Ptr{Nothing}, Ptr{UInt8}), rr, fName)
  if status == false
    error(getLastError())
  end
end

###############################################################################
#                            ToString Routines                                #
###############################################################################
function vectorToString(vecHandle::Ptr{RRVector})
   char_pointer = ccall(dlsym(rrlib, :vectorToString), cdecl, Ptr{UInt8}, (Ptr{RRVector},), vecHandle)
   julia_str = unsafe_string(char_pointer)
   freeText(char_pointer)
   return julia_str
end

function complexVectorToString(vecHandle::Ptr{RRComplexVector})
  char_pointer = ccall(dlsym(rrlib, :complexVectorToString), cdecl, Ptr{UInt8}, (Ptr{RRComplexVector},), vecHandle)
  julia_str = unsafe_string(char_pointer)
  freeText(char_pointer)
  return julia_str
end

function rrCDataToString(rrData::Ptr{RRCData})
  char_pointer = ccall(dlsym(rrlib, :rrCDataToString), cdecl, Ptr{UInt8}, (Ptr{RRCData},), rrData)
  julia_str = unsafe_string(char_pointer)
  freeText(char_pointer)
  return julia_str
end

function matrixToString(matrixHandle::Ptr{RRDoubleMatrix})
  char_pointer = ccall(dlsym(rrlib, :matrixToString), cdecl, Ptr{UInt8}, (Ptr{RRDoubleMatrix},), matrixHandle)
  julia_str = unsafe_string(char_pointer)
  freeText(char_pointer)
  return julia_str
end

function complexMatrixToString(matrixHandle::Ptr{RRComplexMatrix})
  char_pointer = ccall(dlsym(rrlib, :complexMatrixToString), cdecl, Ptr{UInt8}, (Ptr{RRComplexMatrix},), matrixHandle)
  julia_str = unsafe_string(char_pointer)
  freeText(char_pointer)
  return julia_str
end

function stringArrayToString(list::Ptr{RRStringArray})
  char_pointer = ccall(dlsym(rrlib, :stringArrayToString), cdecl, Ptr{UInt8}, (Ptr{RRStringArray},), list)
  julia_str = unsafe_string(char_pointer)
  freeText(char_pointer)
  return julia_str
end

function listToString(list::Ptr{RRList})
   char_pointer = ccall(dlsym(rrlib, :listToString), cdecl, Ptr{UInt8}, (Ptr{RRList},), list)
   julia_str = unsafe_string(char_pointer)
   freeText(char_pointer)
   return julia_str
end

###############################################################################
#                          StringArray Routines                               #
###############################################################################
function getNumberOfStringElements(list::Ptr{RRStringArray})
  return ccall(dlsym(rrlib, :getNumberOfStringElements), cdecl, Int64, (Ptr{RRStringArray},), list)
end

function getStringElement(list::Ptr{RRStringArray}, index::Int64)
  char_pointer = ccall(dlsym(rrlib, :getStringElement), cdecl, Ptr{UInt8}, (Ptr{RRStringArray}, Int64), list, index)
  julia_str = unsafe_string(char_pointer)
  freeText(char_pointer)
  return julia_str
end

###############################################################################
#                          Free Memory Routines                               #
###############################################################################
function freeRRCData(handle::Ptr{RRCData})
  status = ccall(dlsym(rrlib, :freeRRCData), cdecl, Bool, (Ptr{RRCData},), handle)
  if status == false
    (error(getLastError()))
  end
end

function freeText(text::Ptr{UInt8})
  status = ccall(dlsym(rrlib, :freeText), cdecl, Bool, (Ptr{UInt8},), text)
  if status == false
    (error(getLastError()))
  end
end

function freeStringArray(sl::Ptr{RRStringArray})
  status =  ccall(dlsym(rrlib, :freeStringArray), cdecl, Bool, (Ptr{RRStringArray},), sl)
  if status == false
    (error(getLastError()))
  end
end

function freeVector(vector::Ptr{RRVector})
  status = ccall(dlsym(rrlib, :freeVector), cdecl, Bool, (Ptr{RRVector},), vector)
  if status == false
    (error(getLastError()))
  end
end

function freeMatrix(matrix::Ptr{RRDoubleMatrix})
  status = ccall(dlsym(rrlib, :freeMatrix), cdecl, Bool, (Ptr{RRDoubleMatrix},), matrix)
  if status == false
    (error(getLastError()))
  end
end

function convertStringArrayToJuliaArray(list::Ptr{RRStringArray})
  julia_arr = String[]
  try
    num_elem = getNumberOfStringElements(list)
    for i = 1:num_elem
      push!(julia_arr, getStringElement(list, i -1))
    end
  catch e
    throw(e)
  finally
    freeStringArray(list)
  end
  return julia_arr
end

function convertRRVectorToJuliaArray(vector::Ptr{RRVector})
  julia_arr = Float64[]

  try
    num_elem = getVectorLength(vector)
    for i = 1:num_elem
      push!(julia_arr, getVectorElement(vector, i-1))
    end
  catch e
    throw(e)
  finally
    freeVector(vector)
  end
  return julia_arr
end


text:
module TestUtils

using Test
using YaoCompiler
using YaoCompiler.Intrinsics

@device function demo_circ_simp()
    1 => shift(7π / 4)
    1 => H
    1 => Rx(π / 4)
    4 => H
    @ctrl 1 4 => Z
    @ctrl 4 1 => X
    1 => H
    4 => H
    1 => T
    4 => shift(3π / 2)
    4 => X
    1 => H
    4 => S
    4 => X
    2 => S
    @ctrl 2 3 => X
    2 => H
    @ctrl 2 3 => X
    2 => T
    3 => S
    2 => H
    3 => H
    3 => S
    @ctrl 2 3 => X
end

@testset "gate_count" begin
    @test gate_count(demo_circ_simp()) == Dict(
        :ctrl => IdDict(Z => 1, X => 4),
        :gate => IdDict(
            T => 2,
            H => 8,
            Rx(π / 4) => 1,
            X => 2,
            S => 4,
            shift(7π / 4) => 1,
            shift(3π / 2) => 1,
        ),
    )
end

end # TestUtils

ext: jl


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>10.0, :alg=>TMJets(abstol=1e-15, orderT=10, orderQ=2, maxsteps=50_000)),
    "X0"   => Hyperrectangle(low=[-0.5, -0.5], high=[-0.4, -0.4]),
    "info" => Dict(
        "name" => "buckling_column",
        "dim" => 2,
        "linear" => false,
        "hybrid" => false
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using Printf, IterTools

import Base.IteratorSize, Base.iterate, Base.length

iscontseq(n::Integer) = count_zeros(n) == leading_zeros(n) + trailing_zeros(n)
iscontseq(n::BigInt)  = !ismatch(r"0", rstrip(bin(n), '0'))

function makeint2seq(n::Integer)
    idex = collect(1:n)
    function int2seq(m::Integer)
        d = digits(m, base=2, pad=n)
        idex[d .== 1]
    end
    return int2seq
end

struct NCSubSeq{T<:Integer}
    n::T
end

mutable struct NCSubState{T<:Integer}
    m::T
    m2s::Function
end

Base.IteratorSize(::NCSubSeq) = Base.HasLength()
Base.length(a::NCSubSeq) = 2 ^ a.n - a.n * (a.n + 1) ÷ 2 - 1
function Base.iterate(a::NCSubSeq, as::NCSubState=NCSubState(5, makeint2seq(a.n)))
    if 2 ^ a.n - 3 < as.m
        return nothing
    end
    s = as.m2s(as.m)
    as.m += 1
    while iscontseq(as.m)
        as.m += 1
    end
    return (s, as)
end

n = 4
println("Testing NCSubSeq for ", n, " items:\n ", join(NCSubSeq(n), " "))

s = "Rosetta"
cs = split(s, "")
m = 10
n = length(NCSubSeq(length(cs))) - m
println("\nThe first and last ", m, " NC sub-sequences of \"", s, "\":")
for (i, a) in enumerate(NCSubSeq(length(s)))
    i <= m || n < i || continue
    println(@sprintf "%6d %s" i join(cs[a], ""))
    i == m || continue
    println("    .. ......")
end

println("\nThe first and last ", m, " NC sub-sequences of \"", s, "\"")
for x in IterTools.Iterators.flatten([1:10; 20:10:40; big.(50:50:200)])
    @printf "%7d → %d\n" x length(NCSubSeq(x))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    IntegrationVolume(dataset::DataSet{T, I}, spvol::HyperRectVolume{T}, searchpts::Bool = true)::IntegrationVolume{T, I}

*AHMI-internal, not part of stable public API.*

creates an integration region by calculating the point cloud an the volume of the spatial volume.
"""
function IntegrationVolume(
    dataset::DataSet{T, I},
    spvol::HyperRectVolume{T},
    searchpts::Bool = true
)::IntegrationVolume{T, I} where {T<:AbstractFloat, I<:Integer}

    cloud = PointCloud(dataset, spvol, searchpts)
    vol = prod(spvol.hi - spvol.lo)

    IntegrationVolume(cloud, spvol, vol)
end

"""
    modify_integrationvolume!(intvol::IntegrationVolume{T, I}, dataset::DataSet{T, I}, spvol::HyperRectVolume{T}, searchpts::Bool = true)

*AHMI-internal, not part of stable public API.*

updates an integration volume with new boundaries. Recalculates the pointcloud and volume.
"""
function modify_integrationvolume!(
    intvol::IntegrationVolume{T, I},
    dataset::DataSet{T, I},
    spvol::HyperRectVolume{T},
    searchpts::Bool = true
) where {T<:AbstractFloat, I<:Integer}

    if ndims(intvol.spatialvolume) != ndims(spvol)
        intvol.spatialvolume = deepcopy(spvol)
    else
        copy!(intvol.spatialvolume, spvol)
    end

    PointCloud!(intvol.pointcloud, dataset, spvol, searchpts)

    intvol.volume = prod(spvol.hi - spvol.lo)

    nothing
end

function shrink_integrationvol!(
    volume::IntegrationVolume{T, I},
    dataset::DataSet{T, I},
    newrect::HyperRectVolume{T}) where {T<:AbstractFloat, I<:Integer}

    i = volume.pointcloud.points
    for j = 1:i
        inV = true
        for p = 1:dataset.P
            if dataset.data[p, volume.pointcloud.pointIDs[i]] < newrect.lo[p] || dataset.data[p, volume.pointcloud.pointIDs[i]] > newrect.hi[p]
                inV = false
                break
            end
        end
        if !inV
            deleteat!(volume.pointcloud.pointIDs, i)
        end
        i -= 1
    end
    copy!(volume.spatialvolume, newrect)
    volume.pointcloud.points = length(volume.pointcloud.pointIDs)
    volume.volume = prod(newrect.hi .- newrect.lo)

    nothing
end

function update!(
    volume::IntegrationVolume{T, I},
    dataset::DataSet{T, I}) where {T<:AbstractFloat, I<:Integer}

    res = search(dataset, volume.spatialvolume, true)

    _update!(volume, res, false, true)

    #volume.volume doesn't change
    volume.pointcloud.probfactor = exp(volume.pointcloud.maxLogProb - volume.pointcloud.minLogProb)
    volume.pointcloud.probweightfactor = exp(volume.pointcloud.maxWeightProb - volume.pointcloud.minWeightProb)

    @debug "Hyperrectangle updated. Points:\t$(volume.pointcloud.points)\tVolume:\t$(volume.volume)\tProb. Factor:\t$(volume.pointcloud.probfactor)"
    nothing
end

function _update!(
    volume::IntegrationVolume{T, I},
    searchres::SearchResult{T, I},
    addpts::Bool, #if true adds the points to the volume, if false assumes to replace the points in the volume
    searchpts::Bool) where {T<:AbstractFloat, I<:Integer}

    volume.pointcloud.points = searchres.points + (addpts ? volume.pointcloud.points : 0)

    volume.pointcloud.maxLogProb = volume.pointcloud.maxLogProb > searchres.maxLogProb && addpts ? volume.pointcloud.maxLogProb : searchres.maxLogProb
    volume.pointcloud.minLogProb = volume.pointcloud.minLogProb < searchres.minLogProb && addpts ? volume.pointcloud.minLogProb : searchres.minLogProb

    volume.pointcloud.maxWeightProb = volume.pointcloud.maxWeightProb > searchres.maxWeightProb && addpts ? volume.pointcloud.maxWeightProb : searchres.maxWeightProb
    volume.pointcloud.minWeightProb = volume.pointcloud.minWeightProb < searchres.minWeightProb && addpts ? volume.pointcloud.minWeightProb : searchres.minWeightProb

    if searchpts && searchres.points > 0
        start = addpts ? length(volume.pointcloud.pointIDs) + 1 : 1
        resize!(volume.pointcloud.pointIDs, volume.pointcloud.points)
        copyto!(volume.pointcloud.pointIDs, start, searchres.pointIDs, 1)
    end

    nothing
end

function resize_integrationvol(
    original::IntegrationVolume{T, I},
    dataset::DataSet{T, I},
    changed_dim::I,
    newrect::HyperRectVolume{T},
    searchpts::Bool = false)::IntegrationVolume{T, I} where {T<:AbstractFloat, I<:Integer}

    result = deepcopy(original)

    resize_integrationvol!(result, original, dataset, datatre, changed_dim, newrect, searchpts)
end

function resize_integrationvol!(
    result::IntegrationVolume{T, I},
    original::IntegrationVolume{T, I},
    dataset::DataSet{T, I},
    changed_dim::I,
    newrect::HyperRectVolume{T},
    searchpts::Bool,
    searchVol::HyperRectVolume{T}) where {T<:AbstractFloat, I<:Integer}

    copy!(searchVol, newrect)
    increase = true

    #increase
    if original.spatialvolume.lo[changed_dim] > newrect.lo[changed_dim]
        searchVol.hi[changed_dim] = original.spatialvolume.lo[changed_dim]
        searchVol.lo[changed_dim] = newrect.lo[changed_dim]
    elseif original.spatialvolume.hi[changed_dim] < newrect.hi[changed_dim]
        searchVol.lo[changed_dim] = original.spatialvolume.hi[changed_dim]
        searchVol.hi[changed_dim] = newrect.hi[changed_dim]
    else
        increase = false
        if original.spatialvolume.lo[changed_dim] < newrect.lo[changed_dim]
            searchVol.lo[changed_dim] = original.spatialvolume.lo[changed_dim]
            searchVol.hi[changed_dim] = newrect.lo[changed_dim]
        elseif original.spatialvolume.hi[changed_dim] > newrect.hi[changed_dim]
            searchVol.hi[changed_dim] = original.spatialvolume.hi[changed_dim]
            searchVol.lo[changed_dim] = newrect.hi[changed_dim]
        else
            #check if pts inside vol are corrected
            @error "resize_integrationvol(): Volume $original didn't change.", searchVol.lo[changed_dim], "\n", searchVol.hi[changed_dim], "\n", original.spatialvolume.lo[changed_dim], "\n", original.spatialvolume.hi[changed_dim]
        end
    end

    result.pointcloud.points = original.pointcloud.points
    if searchpts
        if result.pointcloud.pointIDs != original.pointcloud.pointIDs
            result.pointcloud.pointIDs = deepcopy(original.pointcloud.pointIDs)
        end
    end

    if increase
        res = search(dataset, searchVol, searchpts)

        _update!(result, res, true, searchpts)
    else
        res = search(dataset, searchVol, searchpts)
        result.pointcloud.points = original.pointcloud.points - res.points
        if searchpts
            newids = search(dataset, newrect, searchpts).pointIDs
            resize!(result.pointcloud.pointIDs, result.pointcloud.points)
            copyto!(result.pointcloud.pointIDs, newids)
        end
    end

    result.volume = prod(newrect.hi - newrect.lo)
    result.pointcloud.probfactor = exp(result.pointcloud.maxLogProb - result.pointcloud.minLogProb)
    result.pointcloud.probweightfactor = exp(result.pointcloud.maxWeightProb - result.pointcloud.minWeightProb)
    copy!(result.spatialvolume, newrect)

    nothing
end

function Base.copy!(
    target::IntegrationVolume{T, I},
    src::IntegrationVolume{T, I}) where {T<:AbstractFloat, I<:Integer}

    target.volume = src.volume

    copy!(target.spatialvolume, src.spatialvolume)
    target.pointcloud.points = src.pointcloud.points

    copy!(target.pointcloud, src.pointcloud)

    nothing
end


text:
"""
$(TYPEDEF)

Default prior on all dimensions of a variable node in the factor graph.  `Prior` is
not recommended when non-Euclidean dimensions are used in variables.
"""
struct Prior{T <: SamplableBelief} <: AbstractPrior
  Z::T
end
Prior(::UniformScaling) = Prior(Normal())

getSample(cf::CalcFactor{<:Prior}, N::Int=1) = (reshape(rand(cf.factor.Z,N),:,N), )

# basic default
(s::CalcFactor{<:Prior})(z, x1) = z .- x1

## packed types are still developed by hand.  Future versions would likely use a @packable macro to write Protobuf safe versions of factors

"""
$(TYPEDEF)

Serialization type for Prior.
"""
mutable struct PackedPrior <: PackedInferenceType
  Z::String
  PackedPrior() = new()
  PackedPrior(z::AS) where {AS <: AbstractString} = new(z)
end
function convert(::Type{PackedPrior}, d::Prior)
  PackedPrior(convert(PackedSamplableBelief, d.Z))
end
function convert(::Type{Prior}, d::PackedPrior)
  Prior(convert(SamplableBelief, d.Z))
end

"""
Converter: Prior -> PackedPrior::Dict{String, Any}

FIXME see DFG #590 for consolidation with Serialization and Marshaling
"""
function convert(::Type{Dict{String, Any}}, prior::IncrementalInference.Prior)
    z = convert(Type{Dict{String, Any}}, prior.Z)
    return Packed_Factor([z], "Prior")
end

"""
Converter: PackedPrior::Dict{String, Any} -> Prior

FIXME see DFG #590 for consolidation on Serialization and Marshaling
"""
function convert(::Type{<:Prior}, prior::Dict{String, Any})
    # Genericize to any packed type next.
    z = prior["measurement"][1]
    z = convert(_evalType(z["distType"]), z)
    return Prior(z)
end


text:
#=
if "../src" ∉  LOAD_PATH
    push!(LOAD_PATH,"../src")
end

import COSMIC

#const data_dir = "/home/guo/haplox/HapBrain/Tuhao/data"
const cosmic = "/home/guo/haplox/HapBrain/data/cosmic/GRch37/v74/CosmicCompleteExport.tsv"
const test_cosmic = "/home/guo/haplox/HapBrain/data/cosmic/10k.tsv"
#const result_dir = "/home/guo/haplox/HapBrain/Tuhao/result"
#const baseinfo = joinpath(result_dir,"baseinfo.txt")
const gene_census = "/home/guo/haplox/HapBrain/data/cosmic/cancer_gene_census.csv"
#const gene_tuhao = joinpath(data_dir, "tuhao_genes.csv")

#const debug = false
#const test = false

COSMIC.build_raw_samples(COSMIC.@test cosmic)
=#


text:
# LLVMBuilder -- reliable LLVM builds all the time.
using BinaryBuilder, Pkg, LibGit2

include("../../fancy_toys.jl")

# Everybody is just going to use the same set of platforms

const llvm_tags = Dict(
    v"6.0.1" => "d359f2096850c68b708bc25a7baca4282945949f",
    v"8.0.1" => "19a71f6bdf2dddb10764939e7f0ec2b98dba76c9",
    v"9.0.1" => "c1a0a213378a458fbea1a5c77b315c7dce08fd05",
    v"10.0.1" => "ef32c611aa214dea855364efd7ba451ec5ec3f74",
    v"11.0.0" => "176249bd6732a8044d457092ed932768724a6f06",
    v"11.0.1" => "43ff75f2c3feef64f9d73328230d34dac8832a91",
    v"12.0.0" => "d28af7c654d8db0b68c175db5ce212d74fb5e9bc",
    v"12.0.1" => "fed41342a82f5a3a9201819a82bf7a48313e296b",
)

const buildscript = raw"""
# We want to exit the program if errors occur.
set -o errexit

if [[ ${target} == *mingw32* ]]; then
    export CCACHE_DISABLE=true
fi

cd ${WORKSPACE}/srcdir/llvm-project/llvm
LLVM_SRCDIR=$(pwd)

# Apply all our patches
if [ -d $WORKSPACE/srcdir/llvm_patches ]; then
for f in $WORKSPACE/srcdir/llvm_patches/*.patch; do
    echo "Applying patch ${f}"
    atomic_patch -p1 ${f}
done
fi

if [ -d $WORKSPACE/srcdir/clang_patches ]; then
cd ${WORKSPACE}/srcdir/llvm-project/clang
for f in $WORKSPACE/srcdir/clang_patches/*.patch; do
    echo "Applying patch ${f}"
    atomic_patch -p1 ${f}
done
fi

if [ -d $WORKSPACE/srcdir/crt_patches ]; then
cd ${WORKSPACE}/srcdir/llvm-project/compiler-rt
for f in $WORKSPACE/srcdir/crt_patches/*.patch; do
    echo "Applying patch ${f}"
    atomic_patch -p1 ${f}
done
fi

if [ -d $WORKSPACE/srcdir/libcxx_patches ]; then
cd ${WORKSPACE}/srcdir/llvm-project/libcxx
for f in $WORKSPACE/srcdir/libcxx_patches/*.patch; do
    echo "Applying patch ${f}"
    atomic_patch -p1 ${f}
done
fi

# Patches from the monorepo
if [ -d $WORKSPACE/srcdir/patches ]; then
cd ${WORKSPACE}/srcdir/llvm-project
for f in $WORKSPACE/srcdir/patches/*.patch; do
    echo "Applying patch ${f}"
    atomic_patch -p1 ${f}
done
fi

# The very first thing we need to do is to build llvm-tblgen for x86_64-linux-muslc
# This is because LLVM's cross-compile setup is kind of borked, so we just
# build the tools natively ourselves, directly.  :/

# Build llvm-tblgen, clang-tblgen, and llvm-config
mkdir ${WORKSPACE}/bootstrap
pushd ${WORKSPACE}/bootstrap
CMAKE_FLAGS=()
CMAKE_FLAGS+=(-DLLVM_TARGETS_TO_BUILD:STRING=host)
CMAKE_FLAGS+=(-DLLVM_HOST_TRIPLE=${MACHTYPE})
CMAKE_FLAGS+=(-DCMAKE_BUILD_TYPE=Release)
if [[ "${LLVM_MAJ_VER}" -gt "11" ]]; then
    CMAKE_FLAGS+=(-DLLVM_ENABLE_PROJECTS='clang;compiler-rt;mlir')
else
    CMAKE_FLAGS+=(-DLLVM_ENABLE_PROJECTS='clang;compiler-rt')
fi
CMAKE_FLAGS+=(-DCMAKE_CROSSCOMPILING=False)
CMAKE_FLAGS+=(-DCMAKE_TOOLCHAIN_FILE=${CMAKE_HOST_TOOLCHAIN})

cmake -GNinja ${LLVM_SRCDIR} ${CMAKE_FLAGS[@]}
if [[ "${LLVM_MAJ_VER}" -gt "11" ]]; then
    ninja -j${nproc} llvm-tblgen clang-tblgen mlir-tblgen mlir-linalg-ods-gen llvm-config
else
    ninja -j${nproc} llvm-tblgen clang-tblgen llvm-config
fi
popd

# Let's do the actual build within the `build` subdirectory
mkdir ${WORKSPACE}/build && cd ${WORKSPACE}/build

# Accumulate these flags outside CMAKE_FLAGS,
# they will be added at the end.
CMAKE_CPP_FLAGS=""
CMAKE_CXX_FLAGS=""
CMAKE_C_FLAGS=""

CMAKE_FLAGS=()

# Release build for best perofrmance
CMAKE_FLAGS+=(-DCMAKE_BUILD_TYPE=Release)
if [[ "${ASSERTS}" == "1" ]]; then
    CMAKE_FLAGS+=(-DLLVM_ENABLE_ASSERTIONS:BOOL=ON)
fi

# build for our host arch and our GPU targets NVidia and AMD
TARGETS=(host NVPTX AMDGPU)
# Add WASM and BPF for LLVM >6
if [[ "${LLVM_MAJ_VER}" != "6" ]]; then
    TARGETS+=(WebAssembly BPF)
fi
LLVM_TARGETS=$(IFS=';' ; echo "${TARGETS[*]}")
CMAKE_FLAGS+=(-DLLVM_TARGETS_TO_BUILD:STRING=$LLVM_TARGETS)

# We mostly care about clang and LLVM
if [[ "${LLVM_MAJ_VER}" -gt "11" ]]; then
    CMAKE_FLAGS+=(-DLLVM_ENABLE_PROJECTS='clang;compiler-rt;lld;mlir')
else
    CMAKE_FLAGS+=(-DLLVM_ENABLE_PROJECTS='clang;compiler-rt;lld')
fi
CMAKE_FLAGS+=(-DLLVM_TOOL_CLANG_TOOLS_EXTRA_BUILD=OFF)

# We want a build with no bindings
CMAKE_FLAGS+=(-DLLVM_BINDINGS_LIST="" )

# Turn off ZLIB and XML2
CMAKE_FLAGS+=(-DLLVM_ENABLE_ZLIB=OFF)
CMAKE_FLAGS+=(-DLLVM_ENABLE_LIBXML2=OFF)

# Disable useless things like docs, terminfo, etc....
CMAKE_FLAGS+=(-DLLVM_INCLUDE_DOCS=Off)
CMAKE_FLAGS+=(-DLLVM_ENABLE_TERMINFO=Off)
CMAKE_FLAGS+=(-DHAVE_HISTEDIT_H=Off)
CMAKE_FLAGS+=(-DHAVE_LIBEDIT=Off)

# We want a shared library
CMAKE_FLAGS+=(-DLLVM_BUILD_LLVM_DYLIB:BOOL=ON)
CMAKE_FLAGS+=(-DLLVM_LINK_LLVM_DYLIB:BOOL=ON)
# set a SONAME suffix for FreeBSD https://github.com/JuliaLang/julia/issues/32462
CMAKE_FLAGS+=(-DLLVM_VERSION_SUFFIX:STRING="jl")

if [[ "${target}" == *linux* || "${target}" == *mingw* ]]; then
    # https://bugs.llvm.org/show_bug.cgi?id=48221
    CMAKE_CXX_FLAGS+="-fno-gnu-unique"
fi

# Install things into $prefix, and make sure it knows we're cross-compiling
CMAKE_FLAGS+=(-DCMAKE_INSTALL_PREFIX=${prefix})
CMAKE_FLAGS+=(-DCMAKE_CROSSCOMPILING=True)

# Julia expects the produced LLVM tools to be installed into tools and not bin
# We can't simply move bin to tools since on MingW64 it will also contain the shlib.
CMAKE_FLAGS+=(-DLLVM_TOOLS_INSTALL_DIR="tools")

# Also build and install utils, since we want FileCheck, and lit
CMAKE_FLAGS+=(-DLLVM_UTILS_INSTALL_DIR="tools")
CMAKE_FLAGS+=(-DLLVM_INCLUDE_UTILS=True -DLLVM_INSTALL_UTILS=True)

# Include perf/oprofile/vtune markers
if [[ ${target} == *linux* ]]; then
    CMAKE_FLAGS+=(-DLLVM_USE_PERF=1)
#     CMAKE_FLAGS+=(-DLLVM_USE_OPROFILE=1)
fi
# if [[ ${target} == *linux* ]] || [[ ${target} == *mingw32* ]]; then
if [[ ${target} == *linux* ]]; then # TODO only LLVM12
    CMAKE_FLAGS+=(-DLLVM_USE_INTEL_JITEVENTS=1)
fi

# Tell LLVM where our pre-built tblgen tools are
CMAKE_FLAGS+=(-DLLVM_TABLEGEN=${WORKSPACE}/bootstrap/bin/llvm-tblgen)
CMAKE_FLAGS+=(-DCLANG_TABLEGEN=${WORKSPACE}/bootstrap/bin/clang-tblgen)
CMAKE_FLAGS+=(-DLLVM_CONFIG_PATH=${WORKSPACE}/bootstrap/bin/llvm-config)
if [[ "${LLVM_MAJ_VER}" -gt "11" ]]; then
    CMAKE_FLAGS+=(-DMLIR_TABLEGEN=${WORKSPACE}/bootstrap/bin/mlir-tblgen)
    CMAKE_FLAGS+=(-DMLIR_LINALG_ODS_GEN=${WORKSPACE}/bootstrap/bin/mlir-linalg-ods-gen)
fi

# Explicitly use our cmake toolchain file
CMAKE_FLAGS+=(-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN})

# Manually set the host triplet, as otherwise on some platforms it tries to guess using
# `ld -v`, which is hilariously wrong.
CMAKE_FLAGS+=(-DLLVM_HOST_TRIPLE=${target})

# Most targets use the actual target string, but we disagree on `aarch64-darwin` and `arm64-darwin`
CMAKE_TARGET=${target}

if [[ "${target}" == *apple* ]]; then
    # On OSX, we need to override LLVM's looking around for our SDK
    CMAKE_FLAGS+=(-DDARWIN_macosx_CACHED_SYSROOT:STRING=/opt/${target}/${target}/sys-root)
    CMAKE_FLAGS+=(-DDARWIN_macosx_OVERRIDE_SDK_VERSION:STRING=10.8)

    # LLVM actually won't build against 10.8, so we bump ourselves up slightly to 10.9
    export MACOSX_DEPLOYMENT_TARGET=10.9
    export LDFLAGS=-mmacosx-version-min=10.9

    # We need to link against libc++ on OSX
    CMAKE_FLAGS+=(-DLLVM_ENABLE_LIBCXX=ON)

    # If we're building for Apple, CMake gets confused with `aarch64-apple-darwin` and instead prefers
    # `arm64-apple-darwin`.  If this issue persists, we may have to change our triplet printing.
    if [[ "${target}" == aarch64* ]]; then
        CMAKE_TARGET=arm64-${target#*-}
    fi
fi

if [[ "${target}" == *apple* ]] || [[ "${target}" == *freebsd* ]]; then
    # On clang-based platforms we need to override the check for ffs because it doesn't work with `clang`.
    export ac_cv_have_decl___builtin_ffs=yes

    # We don't use X-ray on BSD systems
    CMAKE_FLAGS+=(-DCOMPILER_RT_BUILD_XRAY=OFF)
fi

if [[ "${target}" == *mingw* ]]; then
    CMAKE_CPP_FLAGS="${CMAKE_CPP_FLAGS} -remap -D__USING_SJLJ_EXCEPTIONS__ -D__CRT__NO_INLINE"
    # Windows is case-insensitive and some dependencies take full advantage of that
    echo "BaseTsd.h basetsd.h" >> /opt/${target}/${target}/include/header.gcc
fi

if [[ "${target}" == *musl* ]]; then
    # Taken from https://git.alpinelinux.org/cgit/aports/tree/main/compiler-rt/APKBUILD
    CMAKE_FLAGS+=(-DCOMPILER_RT_INCLUDE_TESTS=ON)
    CMAKE_FLAGS+=(-DCOMPILER_RT_BUILD_SANITIZERS=OFF)
    CMAKE_FLAGS+=(-DCOMPILER_RT_BUILD_XRAY=OFF)
fi

if [[ "${target}" == *freebsd* ]]; then
    # On FreeBSD, we must force even statically-linked code to have -fPIC
    CMAKE_FLAGS+=(-DCMAKE_POSITION_INDEPENDENT_CODE=TRUE)
fi

# Tell LLVM which compiler target to use, because it loses track for some reason
CMAKE_FLAGS+=(-DCMAKE_C_COMPILER_TARGET=${CMAKE_TARGET})
CMAKE_FLAGS+=(-DCMAKE_CXX_COMPILER_TARGET=${CMAKE_TARGET})
CMAKE_FLAGS+=(-DCMAKE_ASM_COMPILER_TARGET=${CMAKE_TARGET})

cmake -GNinja ${LLVM_SRCDIR} ${CMAKE_FLAGS[@]} -DCMAKE_CXX_FLAGS="${CMAKE_CPP_FLAGS} ${CMAKE_CXX_FLAGS}" -DCMAKE_C_FLAGS="${CMAKE_CPP_FLAGS} ${CMAKE_CXX_FLAGS}"
ninja -j${nproc} -vv

# Install!
ninja install

# Life is harsh on Windows and dynamic libraries are
# expected to live alongside the binaries. So we have
# to copy the *.dll from bin/ to tools/ as well...
if [[ "${target}" == *mingw* ]]; then
    cp ${prefix}/bin/*.dll ${prefix}/tools/
fi

# Work around llvm-config bug by creating versioned symlink to libLLVM
# https://github.com/JuliaLang/julia/pull/30033
if [[ "${target}" == *darwin* ]]; then
    LLVM_VER=$(${WORKSPACE}/bootstrap/bin/llvm-config --version | cut -d. -f1-2)
    ln -s libLLVM.dylib ${prefix}/lib/libLLVM-${LLVM_VER}.dylib
fi

# Lit is a python dependency and there is no proper install target
cp -r ${LLVM_SRCDIR}/utils/lit ${prefix}/tools/

install_license ${WORKSPACE}/srcdir/llvm-project/llvm/LICENSE.TXT
"""

# Also define some scripts for extraction:
const libllvmscript = raw"""
# First, find (true) LLVM library directory in ~/.artifacts somewhere
LLVM_ARTIFACT_DIR=$(dirname $(dirname $(realpath ${prefix}/tools/opt${exeext})))

# Clear out our `${prefix}`
rm -rf ${prefix}/*

# Copy over `llvm-config`, `libLLVM` and `include`, specifically.
mkdir -p ${prefix}/include ${prefix}/tools ${libdir} ${prefix}/lib
mv -v ${LLVM_ARTIFACT_DIR}/include/llvm* ${prefix}/include/
mv -v ${LLVM_ARTIFACT_DIR}/tools/llvm-config* ${prefix}/tools/
mv -v ${LLVM_ARTIFACT_DIR}/$(basename ${libdir})/*LLVM*.${dlext}* ${libdir}/
mv -v ${LLVM_ARTIFACT_DIR}/lib/*LLVM*.a ${prefix}/lib
install_license ${LLVM_ARTIFACT_DIR}/share/licenses/LLVM_full*/*
"""

const clangscript = raw"""
# First, find (true) LLVM library directory in ~/.artifacts somewhere
LLVM_ARTIFACT_DIR=$(dirname $(dirname $(realpath ${prefix}/tools/opt${exeext})))

# Clear out our `${prefix}`
rm -rf ${prefix}/*

# Copy over `clang`, `libclang` and `include`, specifically.
mkdir -p ${prefix}/include ${prefix}/bin ${libdir} ${prefix}/lib
mv -v ${LLVM_ARTIFACT_DIR}/include/clang* ${prefix}/include/
mv -v ${LLVM_ARTIFACT_DIR}/bin/clang* ${prefix}/bin/
mv -v ${LLVM_ARTIFACT_DIR}/$(basename ${libdir})/libclang*.${dlext}* ${libdir}/
mv -v ${LLVM_ARTIFACT_DIR}/lib/libclang*.a ${prefix}/lib
mv -v ${LLVM_ARTIFACT_DIR}/lib/clang ${prefix}/lib/clang
install_license ${LLVM_ARTIFACT_DIR}/share/licenses/LLVM_full*/*
"""

const mlirscript = raw"""
# First, find (true) LLVM library directory in ~/.artifacts somewhere
LLVM_ARTIFACT_DIR=$(dirname $(dirname $(realpath ${prefix}/tools/opt${exeext})))

# Clear out our `${prefix}`
rm -rf ${prefix}/*

# Copy over `libMLIR` and `include`, specifically.
mkdir -p ${prefix}/include ${prefix}/tools ${libdir} ${prefix}/lib
mv -v ${LLVM_ARTIFACT_DIR}/include/mlir* ${prefix}/include/
mv -v ${LLVM_ARTIFACT_DIR}/tools/mlir* ${prefix}/tools/
mv -v ${LLVM_ARTIFACT_DIR}/$(basename ${libdir})/*MLIR*.${dlext}* ${libdir}/
mv -v ${LLVM_ARTIFACT_DIR}/$(basename ${libdir})/*mlir*.${dlext}* ${libdir}/
install_license ${LLVM_ARTIFACT_DIR}/share/licenses/LLVM_full*/*
"""

const llvmscript = raw"""
# First, find (true) LLVM library directory in ~/.artifacts somewhere
LLVM_ARTIFACT_DIR=$(dirname $(dirname $(realpath ${prefix}/tools/opt${exeext})))

# Clear out our `${prefix}`
rm -rf ${prefix}/*

# Copy over everything, but eliminate things already put inside `Clang_jll` or `libLLVM_jll`:
mv -v ${LLVM_ARTIFACT_DIR}/* ${prefix}/
rm -vrf ${prefix}/include/{clang*,llvm*,mlir*}
rm -vrf ${prefix}/bin/{clang*,llvm-config,mlir*}
rm -vrf ${libdir}/libclang*.${dlext}*
rm -vrf ${libdir}/*LLVM*.${dlext}*
rm -vrf ${libdir}/*MLIR*.${dlext}*
rm -vrf ${prefix}/lib/*LLVM*.a
rm -vrf ${prefix}/lib/libclang*.a
rm -vrf ${prefix}/lib/clang
rm -vrf ${prefix}/lib/mlir
"""

function configure_build(ARGS, version; experimental_platforms=false, assert=false)
    # Parse out some args
    if "--assert" in ARGS
        assert = true
        deleteat!(ARGS, findfirst(ARGS .== "--assert"))
    end
    sources = [
        GitSource("https://github.com/llvm/llvm-project.git", llvm_tags[version]),
        DirectorySource("./bundled"),
    ]

    platforms = expand_cxxstring_abis(supported_platforms(;experimental=experimental_platforms))
    products = [
        LibraryProduct("libclang", :libclang, dont_dlopen=true),
        LibraryProduct(["LLVM", "libLLVM"], :libllvm, dont_dlopen=true),
        LibraryProduct(["LTO", "libLTO"], :liblto, dont_dlopen=true),
        ExecutableProduct("llvm-config", :llvm_config, "tools"),
        ExecutableProduct("clang", :clang, "bin"),
        ExecutableProduct("opt", :opt, "tools"),
        ExecutableProduct("llc", :llc, "tools"),
    ]
    if version >= v"8"
        push!(products, ExecutableProduct("llvm-mca", :llvm_mca, "tools"))
    end
    if version >= v"12"
        push!(products, LibraryProduct(["MLIR", "libMLIR"], :mlir, dont_dlopen=true))
        push!(products, LibraryProduct(["MLIRPublicAPI", "libMLIRPublicAPI"], :mlir_public, dont_dlopen=true))
        push!(products, LibraryProduct("libclang-cpp", :libclang_cpp, dont_dlopen=true))
    end

    name = "LLVM_full"
    config = "LLVM_MAJ_VER=$(version.major)\n"
    if assert
        config *= "ASSERTS=1\n"
        name = "$(name)_assert"
    end
    # Dependencies that must be installed before this package can be built
    # TODO: Zlib, LibXML2
    dependencies = Dependency[]
    return name, version, sources, config * buildscript, platforms, products, dependencies
end

function configure_extraction(ARGS, LLVM_full_version, name, libLLVM_version=nothing; experimental_platforms=false, assert=false)
    if isempty(LLVM_full_version.build)
        error("You must lock an extracted LLVM build to a particular LLVM_full build number!")
    end
    if name != "libLLVM" && (libLLVM_version === nothing || isempty(libLLVM_version.build))
        error("You must lock an extracted LLVM build to a particular libLLVM build number!")
    end
    version = VersionNumber(LLVM_full_version.major, LLVM_full_version.minor, LLVM_full_version.patch)
    compat_version = "$(version.major).$(version.minor).$(version.patch)"
    if name == "libLLVM"
        script = libllvmscript
        products = [
            LibraryProduct(["LLVM", "libLLVM"], :libllvm, dont_dlopen=true),
            ExecutableProduct("llvm-config", :llvm_config, "tools"),
        ]
    elseif name == "Clang"
        script = clangscript
        products = [
            LibraryProduct("libclang", :libclang, dont_dlopen=true),
            LibraryProduct("libclang-cpp", :libclang_cpp, dont_dlopen=true),
            ExecutableProduct("clang", :clang, "bin"),
        ]
    elseif name == "MLIR"
        script = mlirscript
        products = [
            LibraryProduct("libMLIR", :libMLIR, dont_dlopen=true),
        ]
    elseif name == "LLVM"
        script = llvmscript
        products = [
            LibraryProduct(["LTO", "libLTO"], :liblto, dont_dlopen=true),
            ExecutableProduct("opt", :opt, "tools"),
            ExecutableProduct("llc", :llc, "tools"),
        ]
        if version >= v"8"
            push!(products, ExecutableProduct("llvm-mca", :llvm_mca, "tools"))
        end
    end
    platforms = expand_cxxstring_abis(supported_platforms(;experimental=experimental_platforms))

    dependencies = BinaryBuilder.AbstractDependency[]

    # Parse out some args
    if "--assert" in ARGS
        assert = true
        deleteat!(ARGS, findfirst(ARGS .== "--assert"))
    end

    if assert
        push!(dependencies, BuildDependency(get_addable_spec("LLVM_full_assert_jll", LLVM_full_version)))
        if name in ("Clang", "LLVM", "MLIR")
            push!(dependencies, Dependency("libLLVM_assert_jll", libLLVM_version, compat=compat_version))
        end
        name = "$(name)_assert"
    else
        push!(dependencies, BuildDependency(get_addable_spec("LLVM_full_jll", LLVM_full_version)))
        if name in ("Clang", "LLVM", "MLIR")
            push!(dependencies, Dependency("libLLVM_jll", libLLVM_version, compat=compat_version))
        end
    end

    return name, version, [], script, platforms, products, dependencies
end

max_stars_repo_name: xue-cheng/Yggdrasil

max_issues_repo_name: xue-cheng/Yggdrasil

max_forks_repo_name: xue-cheng/Yggdrasil


text:
export LazyInitialized

"""
    LazyInitialized{T}()

A thread-safe, lazily-initialized wrapper for a value of type `T`. Initialize and fetch the
value by calling `get!`. The constructor is ensured to only be called once.

This type is intended for lazy initialization of e.g. global structures, without using
`__init__`. It is similar to protecting accesses using a lock, but is much cheaper.

"""
struct LazyInitialized{T}
    # 0: uninitialized
    # 1: initializing
    # 2: initialized
    guard::Threads.Atomic{Int}
    value::Base.RefValue{T}

    LazyInitialized{T}() where {T} =
        new(Threads.Atomic{Int}(0), Ref{T}())
end

@inline function Base.get!(constructor, x::LazyInitialized; hook=nothing)
    while x.guard[] != 2
        initialize!(x, constructor, hook)
    end
    assume(isassigned(x.value)) # to get rid of the check
    x.value[]
end

@noinline function initialize!(x::LazyInitialized, constructor::F1, hook::F2) where {F1, F2}
    status = Threads.atomic_cas!(x.guard, 0, 1)
    if status == 0
        try
          x.value[] = constructor()
          x.guard[] = 2
        catch
          x.guard[] = 0
          rethrow()
        end

        if hook !== nothing
          hook()
        end
    else
        yield()
    end
    return
end


text:
function rising_cutoff(t)
    r = 0
    if t <= -1
        r = 0
    elseif t >= 1
        r = 1
    else
        r = sin(π / 4 * (1 + sin(π / 2 * t)))
    end
    return r
end

function standardize_fiedler(v)
    v1 = deepcopy(v)
    v1 ./= norm(v1, 2)
    if v1[1] < 0
        v1 = -v1
    end
    return round.(v1; digits = 14)
end

function Lrw_eigenvec(W; nev = 6)
    # N = size(W, 1)
    # D = Diagonal(sum(W, dims = 1)[:])
    # try
    #     val, vtmp = eigs(D - W, D,
    #                      nev = nev, sigma = eps(), v0 = ones(N) / sqrt(N))
    #     vtmp = vtmp[:, sortperm(val)[2:end]]
    #     vtmp ./= sqrt.(sum(vtmp.^2; dims = 1))
    #     vtmp *= Diagonal(1 .- (vtmp[1, :] .< 0) .* 2)
    # catch
    #
    # return round.(vtmp; digits = 14)
    deg = sum(W, dims = 1)[:]  # weighted degree vector
    if minimum(deg) <= 10^3 * eps()
        L = Diagonal(deg) - W
        𝛌, 𝚽 = eigen(Matrix(L))
        standardize_eigenvector_signs!(𝚽)
        return round.(𝚽[:, 2:nev]; digits = 14)
    end
    Lsym = Diagonal(deg.^(-1/2)) * (Diagonal(deg) - W) * Diagonal(deg.^(-1/2))
    𝛌sym, 𝚽sym = eigen(Matrix(Lsym))
    𝚽rw = Diagonal(deg.^(-1/2)) * 𝚽sym
    𝚽rw ./= sqrt.(sum(𝚽rw.^2; dims = 1))
    𝚽rw *= Diagonal(1 .- (𝚽rw[1, :] .< 0) .* 2)
    return round.(𝚽rw[:, 2:nev]; digits = 14)
end

function sortnodes_inmargin(active_region, Na, v, W, idx; sign = :positive)
    ind = sortperm(v[active_region]; rev = (sign == :positive))
    # if there is a tie, use the more dims eigenmaps and sort by lexical order
    if length(unique(v[ind])) < length(active_region)
        emb = Lrw_eigenvec(W[idx, idx]; nev = 3)'
        # for the second dimension coordinate we sort
        # the pos region and the neg region in the same order
        emb_tp = [Tuple(emb[:, i] .* [1, (sign == :positive) * 2 - 1]) for i in 1:length(idx)]
        ind = sortperm(emb_tp[active_region]; rev = (sign == :positive))
    end
    return ind[(end - Na + 1):end]
end

function find_pairinds(W; ϵ::Float64 = 0.2, idx = 1:size(W, 1), used_node = Set())
    v = partition_fiedler(W[idx, idx])[2]
    v = standardize_fiedler(v)
    vmax = norm(v, Inf)
    N = length(v)
    pos_active_region = (Int64)[]
    neg_active_region = (Int64)[]

    tol = 10^3 * eps()
    tol *= ((sum(v .>= tol) > sum(v .<= -tol)) * 2 - 1)

    for i in 1:N
        if idx[i] ∈ used_node
            continue
        end
        if tol <= v[i] < ϵ * vmax
            push!(pos_active_region, i)
        end
        if -ϵ * vmax < v[i] < tol
            push!(neg_active_region, i)
        end
    end
    Np = length(pos_active_region)
    Nn = length(neg_active_region)
    Na = min(Nn, Np)  # number of pair inds in action region

    # indp = sortperm(v[pos_active_region]; rev = true)[(end - Na + 1):end]
    # indn = sortperm(v[neg_active_region])[(end - Na + 1):end]
    indp = sortnodes_inmargin(pos_active_region, Na, v, W, idx; sign = :positive)
    indn = sortnodes_inmargin(neg_active_region, Na, v, W, idx; sign = :negative)
    pair_inds = vcat(idx[pos_active_region[indp]]',
                     idx[neg_active_region[indn]]')

    return pair_inds, v
end

function pair_inds_shadding(W, GP; ϵ = 0.3, J = 1)
    rs = GP.rs
    inds = GP.inds
    (N, jmax) = Base.size(inds)

    used_node = Set()
    shading = zeros(N)
    for j = 1:J
        regioncount = count(!iszero, rs[:, j]) - 1
        for r = 1:regioncount
            indr = rs[r, j]:(rs[r + 1, j] - 1)
            pair_inds = find_pairinds(W; ϵ = ϵ, idx = inds[indr, j], used_node = used_node)[1]
            for i in 1:size(pair_inds, 2)
                pv, nv = pair_inds[:, i]
                # if pv ∈ used_node || nv ∈ used_node
                #     continue
                # end
                if shading[pv] == 0
                    shading[pv] = J + 3 - j
                end
                if shading[nv] == 0
                    shading[nv] = -(J + 3 - j)
                end
            end
            union!(used_node, Set(pair_inds[:]))
        end
    end
    return shading
end

function keep_folding!(U, used_node, W, GP; ϵ = 0.3, j = 1)
    if j == 0
        return U
    end

    rs = GP.rs
    inds = GP.inds
    N = Base.size(inds, 1)

    regioncount = count(!iszero, rs[:, j]) - 1
    for r = 1:regioncount
        indr = rs[r, j]:(rs[r + 1, j] - 1)
        if length(indr) == 1
            continue
        end
        pair_inds, v = find_pairinds(W; ϵ = ϵ, idx = inds[indr, j], used_node = used_node)
        vmax = norm(v, Inf)
        for i in 1:size(pair_inds, 2)
            pv, nv = pair_inds[:, i]
            if pv ∈ used_node || nv ∈ used_node
                continue
            end
            # use the half distance between the 1D embeddings,
            # i.e., t = const⋅(v[pv] - v[nv] / 2), to compute the rising
            # cutoff function, which satisfy r(t)^2 + r(-t)^2 = 1
            t = (v[findfirst(inds[indr, j] .== pv)] - v[findfirst(inds[indr, j] .== nv)]) / (2 * ϵ * vmax)
            # t = v[findfirst(inds[indr, j] .== pv)] / (ϵ * vmax)
            U[pv, pv] = rising_cutoff(t)
            U[pv, nv] = rising_cutoff(-t)
            U[nv, pv] = -rising_cutoff(-t)
            U[nv, nv] = rising_cutoff(t)
        end
        union!(used_node, Set(pair_inds[:]))
    end
end

function unitary_folding_operator(W, GP; ϵ = 0.3, J = 1)
    rs = GP.rs
    inds = GP.inds
    (N, jmax) = Base.size(inds)

    U = Matrix{Float64}(I, N, N)
    used_node = Set()

    for j = 1:J
        keep_folding!(U, used_node, W, GP; ϵ = ϵ, j = j)
    end

    return U
end

"""
    lp_ngwp(𝚽::Matrix{Float64}, W_dual::SparseMatrixCSC{Float64, Int64},
            GP_dual::GraphPart; ϵ::Float64 = 0.3)

construct the lapped NGWP and GP.tag in place.

# Input Arguments
- `𝚽::Matrix{Float64}`: graph Laplacian eigenvectors
- `W_dual::SparseMatrixCSC{Float64, Int64}`: weight matrix of the dual graph
- `GP_dual::GraphPart`: GraphPart object of the dual graph

# Output Argument
- `wavelet_packet::Array{Float64,3}`: the lapped NGWP dictionary. The first
    index is for selecting wavelets at a fixed level; the second index is for
    selecting the level `j`; the third index is for selecting elements in the
    wavelet vector.

"""
function lp_ngwp(𝚽::Matrix{Float64}, W_dual::SparseMatrixCSC{Float64, Int64},
                 GP_dual::GraphPart; ϵ::Float64 = 0.3)
    rs = GP_dual.rs
    inds = GP_dual.inds
    (N, jmax) = Base.size(inds)

    GP_dual.tag = zeros(Int, N, jmax)
    GP_dual.tag[:, 1] = Vector{Int}(0:(N - 1))

    Uf = Matrix{Float64}(I, N, N)
    used_node = Set()

    wavelet_packet = zeros(N, jmax, N)
    wavelet_packet[:, 1, :] = Matrix{Float64}(I, N, N)
    for j = 2:jmax
        regioncount = count(!iszero, rs[:, j]) - 1
        # Uf = unitary_folding_operator(W_dual, GP_dual; ϵ = ϵ, J = j - 1)
        keep_folding!(Uf, used_node, W_dual, GP_dual; ϵ = ϵ, j = j - 1)
        for r = 1:regioncount
            indr = rs[r, j]:(rs[r + 1, j] - 1)
            GP_dual.tag[indr, j] = Vector{Int}(0:(length(indr) - 1))
            wavelet_packet[indr, j, :] = const_meyer_wavelets(𝚽, Uf; idx = inds[indr, j])'
        end
    end
    return wavelet_packet
end

function const_meyer_wavelets(𝚽, Uf; idx = 1:size(Uf, 1))
    N = size(𝚽, 1)
    # 1. nice interpretable way by smooth orthogonal projector
    # assemble smooth orthogonal projector custom to nodes `idx`
    P = Uf[idx, :]' * Uf[idx, idx]
    if diag(P) == χ(idx, N)
        B = 𝚽[:, idx]
    else
        # folding the eigenspace, i.e., 𝚽's column space
        Y = 𝚽 * P
        # find its column space's orthogonal basis
        B = svd(Y).U
    end
    # 2. alternative faster way by orthogonal unfolding operator
    # B = 𝚽 * Uf[:, idx]
    # perform varimax rotation to get the meyer_wavelets
    Wavelets = varimax(B)
    return Wavelets
end

"""
    lp_ngwp_analysis(G::GraphSig, 𝚽::Matrix{Float64}, W_dual::SparseMatrixCSC{Float64, Int64},
            GP_dual::GraphPart; ϵ::Float64 = 0.3)

perform the LP-NGWP transform of the graph signal(s) `G.f`.

# Input Arguments
- `G::GraphSig`: a GraphSig object
- `𝚽::Matrix{Float64}`: graph Laplacian eigenvectors
- `W_dual::SparseMatrixCSC{Float64, Int64}`: weight matrix of the dual graph
- `GP_dual::GraphPart`: GraphPart object of the dual graph

# Output Argument
- `wavelet_packet::Array{Float64,3}`: the lapped NGWP dictionary. The first
    index is for selecting wavelets at a fixed level; the second index is for
    selecting the level `j`; the third index is for selecting elements in the
    wavelet vector.

"""
function lp_ngwp_analysis(G::GraphSig, 𝚽::Matrix{Float64}, W_dual::SparseMatrixCSC{Float64, Int64},
                 GP_dual::GraphPart; ϵ::Float64 = 0.3)
    rs = GP_dual.rs
    inds = GP_dual.inds
    (N, jmax) = Base.size(inds)
    fcols = size(G.f, 2)

    GP_dual.tag = zeros(Int, N, jmax)
    GP_dual.tag[:, 1] = Vector{Int}(0:(N - 1))

    Uf = Matrix{Float64}(I, N, N)
    used_node = Set()

    dmatrix = zeros(N, jmax, fcols)
    dmatrix[:, 1, :] = G.f
    for j = 2:jmax
        regioncount = count(!iszero, rs[:, j]) - 1
        keep_folding!(Uf, used_node, W_dual, GP_dual; ϵ = ϵ, j = j - 1)
        for r = 1:regioncount
            indr = rs[r, j]:(rs[r + 1, j] - 1)
            GP_dual.tag[indr, j] = Vector{Int}(0:(length(indr) - 1))
            dmatrix[indr, j, :] = const_meyer_wavelets(𝚽, Uf; idx = inds[indr, j])' * G.f
        end
    end
    return dmatrix
end


text:
module Blocks

using Base.FS
using Compat

importall   Base
import      Base: peek, throwto, open

if isless(Base.VERSION, v"0.4.0-")
import      Base: AsyncStream, localpart
typealias LibuvStream AsyncStream
elseif isless(Base.VERSION, v"0.4.0-rc3+30")
import      Base: AsyncStream, |>
typealias LibuvStream AsyncStream
else
import      Base: LibuvStream, |>
end

export      Block, |>, prepare, @prepare, BlockableIO,
            blocks, affinities, localpart,
            as_it_is, as_io, as_recordio, as_wordio, as_lines, as_bufferedio, as_bytearray,

            map, mapreduce, pmap, pmapreduce,

            BlockIO, close, eof, read!, write, readbytes, peek,
            readall, flush, nb_available, position, filesize, seek, seekend, seekstart, skip

include("pmap.jl")
include("block_io.jl")
include("block_framework.jl")
include("mapreduce.jl")

# sub modules
#include("dataframe_blocks.jl")
#include("hdfs_blocks.jl")
include("matop_blocks.jl")

end


text:
# Date: July 1, 2018
# Author: Bernhard König
# Title: Short example on how categorical variables with high cardinalilty can be treated and preprocessed for these algorithms
# This snippet shall merely show the functionality of the two very simple Julia funcitons mapToOther! and getCounts
# We refer to the general literature on data preprocessing (specifically for high cardinality categorical data) and feature generation

using Distributed
# Distributed.@everywhere using Revise
Distributed.@everywhere import CSV
Distributed.@everywhere import DataFrames
Distributed.@everywhere import DataFrames: DataFrame

Distributed.@everywhere using DecisionTrees

elt = [Int,	Float64,	Float64,	Float64,	Float64,	Float64,	Int,	String,	String,	String,	Int,	String,	String,	String,	String,	String,	String,	String,	String,	String,	String,	Int,	String,	Int,	String,	String,	Int,	String,	String,	String,	String,	String,	Int,	String,	String,	String,	String,	Int,	String,	String,	Int,	Int,	String,	Int,	Int,	String,	Int,	Int,	String,	Int,	String,	String,	String,	Int,	String,	String,	String,	String,	Int,	Int,	Int,	String,	String,	String,	String,	String,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Int,	Float64,	Int,	Float64,	Float64,	Float64,	Float64,	Int,	Int,	Int,	Int,	Int,	Int]
datafile = joinpath("data", "GermanMotorPremiums", "data1small.csv")
@assert isfile(datafile);
@time df_tmp = CSV.read(datafile, DataFrame, types=elt, strict=true);

selected_explanatory_vars = ["ART_DES_WOHNEIGENTUM","GEBURTSDATUM","FAMILIENSTAND","NATIONALITAET","PLZ_WOHNORT","GESCHLECHT","FINANZIERUNGSART","STADT","KENNZEICHEN"]

for x in selected_explanatory_vars
    @show x, eltype(df_tmp[Symbol(x)])
end

# consider the unique values of PLZ_WOHNORT
uq = unique(df_tmp[:PLZ_WOHNORT])
length(uq)

# we can analyse the frequency with this helper function
counts, freqs, vals, keep30 = getCounts(df_tmp[:PLZ_WOHNORT], threshold=30)
# keep is the 30 largest values

# alternatively we can keep all values with are more frequent, than x%
# Let us keep all values more frequent than 1/1000
counts, freqs, vals, keep = getCounts(df_tmp[:PLZ_WOHNORT], threshold=-0.002)

# the we can easily modify the data and overwrite the infrequent values with a custom value

thiscol = df_tmp[:PLZ_WOHNORT]
DecisionTrees.mapToOther!(thiscol,keep30,"some rare PLZ")
unique(df_tmp[:PLZ_WOHNORT])


text:
module ControlSystems

export  LTISystem,
        AbstractStateSpace,
        StateSpace,
        HeteroStateSpace,
        TransferFunction,
        DelayLtiSystem,
        Continuous,
        Discrete,
        ss,
        tf,
        zpk,
        LQG,
        isproper,
        # Linear Algebra
        balance,
        care,
        dare,
        dlyap,
        lqr,
        dlqr,
        kalman,
        dkalman,
        lqg,
        lqgi,
        covar,
        norm,
        hinfnorm,
        linfnorm,
        gram,
        ctrb,
        obsv,
        place,
        luenberger,
        # Model Simplification
        reduce_sys,
        sminreal,
        minreal,
        balreal,
        baltrunc,
        similarity_transform,
        prescale,
        innovation_form,
        # Stability Analysis
        isstable,
        pole,
        tzero,
        dcgain,
        zpkdata,
        damp,
        dampreport,
        markovparam,
        margin,
        delaymargin,
        gangoffour,
        # Connections
        append,
        series,
        parallel,
        feedback,
        feedback2dof,
        starprod,
        lft,
        # Discrete
        c2d,
        c2d_x0map,
        d2c,
        # Time Response
        step,
        impulse,
        lsim,
        solve,
        Simulator,
        # Frequency Response
        freqresp, freqrespv,
        evalfr,
        bode, bodev,
        nyquist, nyquistv,
        sigma, sigmav,
        # delay systems
        delay,
        pade,
        # demo systems
        ssrand,
        DemoSystems, # A module containing some example systems
        # utilities
        num,    #Deprecated
        den,    #Deprecated
        numvec,
        denvec,
        numpoly,
        denpoly,
        iscontinuous,
        isdiscrete

# QUESTION: are these used? LaTeXStrings, Requires, IterTools
using Plots, LaTeXStrings, LinearAlgebra
import Polynomials
import Polynomials: Polynomial, coeffs
using OrdinaryDiffEq, DelayDiffEq
export Plots
import Base: +, -, *, /, (==), (!=), isapprox, convert, promote_op
import Base: getproperty, getindex
import Base: exp # for exp(-s)
import LinearAlgebra: BlasFloat
export lyap # Make sure LinearAlgebra.lyap is available
import Printf, Colors
import DSP: conv
using MacroTools

abstract type AbstractSystem end

include("types/TimeEvolution.jl")
## Added interface:
#   timeevol(Lti) -> TimeEvolution (not exported)

include("types/Lti.jl")

include("types/SisoTf.jl")

# Transfer functions and tranfer function elemements
include("types/TransferFunction.jl")
include("types/SisoTfTypes/SisoZpk.jl")
include("types/SisoTfTypes/SisoRational.jl")
include("types/SisoTfTypes/promotion.jl")
include("types/SisoTfTypes/conversion.jl")

include("types/StateSpace.jl")

# TODO Sample time
include("types/PartionedStateSpace.jl")
include("types/DelayLtiSystem.jl")

# Convenience constructors
include("types/tf.jl")
include("types/zpk.jl")

include("types/lqg.jl") # QUESTION: is it really motivated to have an LQG type?

include("utilities.jl")

include("types/promotion.jl")
include("types/conversion.jl")
include("connections.jl")

# Analysis
include("freqresp.jl")
include("timeresp.jl")

include("matrix_comps.jl")
include("simplification.jl")

include("discrete.jl")
include("analysis.jl")
include("synthesis.jl")

include("simulators.jl")
include("pid_design.jl")

include("demo_systems.jl")

include("delay_systems.jl")

include("plotting.jl")

@deprecate num numvec
@deprecate den denvec
@deprecate norminf hinfnorm
@deprecate diagonalize(s::AbstractStateSpace, digits) diagonalize(s::AbstractStateSpace)

function covar(D::Union{AbstractMatrix,UniformScaling}, R)
    @warn "This call is deprecated due to ambiguity, use covar(ss(D), R) or covar(ss(D, Ts), R) instead"
    D*R*D'
end

# The path has to be evaluated upon initial import
const __CONTROLSYSTEMS_SOURCE_DIR__ = dirname(Base.source_path())

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 RateCardQueryParameters <: SwaggerModel
    OfferDurableId::Any # spec type: Union{ Nothing, String } # spec name: OfferDurableId
    Currency::Any # spec type: Union{ Nothing, String } # spec name: Currency
    Locale::Any # spec type: Union{ Nothing, String } # spec name: Locale
    RegionInfo::Any # spec type: Union{ Nothing, String } # spec name: RegionInfo

    function RateCardQueryParameters(;OfferDurableId=nothing, Currency=nothing, Locale=nothing, RegionInfo=nothing)
        o = new()
        validate_property(RateCardQueryParameters, Symbol("OfferDurableId"), OfferDurableId)
        setfield!(o, Symbol("OfferDurableId"), OfferDurableId)
        validate_property(RateCardQueryParameters, Symbol("Currency"), Currency)
        setfield!(o, Symbol("Currency"), Currency)
        validate_property(RateCardQueryParameters, Symbol("Locale"), Locale)
        setfield!(o, Symbol("Locale"), Locale)
        validate_property(RateCardQueryParameters, Symbol("RegionInfo"), RegionInfo)
        setfield!(o, Symbol("RegionInfo"), RegionInfo)
        o
    end
end # type RateCardQueryParameters

const _property_map_RateCardQueryParameters = Dict{Symbol,Symbol}(Symbol("OfferDurableId")=>Symbol("OfferDurableId"), Symbol("Currency")=>Symbol("Currency"), Symbol("Locale")=>Symbol("Locale"), Symbol("RegionInfo")=>Symbol("RegionInfo"))
const _property_types_RateCardQueryParameters = Dict{Symbol,String}(Symbol("OfferDurableId")=>"String", Symbol("Currency")=>"String", Symbol("Locale")=>"String", Symbol("RegionInfo")=>"String")
Base.propertynames(::Type{ RateCardQueryParameters }) = collect(keys(_property_map_RateCardQueryParameters))
Swagger.property_type(::Type{ RateCardQueryParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_RateCardQueryParameters[name]))}
Swagger.field_name(::Type{ RateCardQueryParameters }, property_name::Symbol) =  _property_map_RateCardQueryParameters[property_name]

function check_required(o::RateCardQueryParameters)
    (getproperty(o, Symbol("OfferDurableId")) === nothing) && (return false)
    (getproperty(o, Symbol("Currency")) === nothing) && (return false)
    (getproperty(o, Symbol("Locale")) === nothing) && (return false)
    (getproperty(o, Symbol("RegionInfo")) === nothing) && (return false)
    true
end

function validate_property(::Type{ RateCardQueryParameters }, name::Symbol, val)
    if name === Symbol("OfferDurableId")
    end
end


text:
using QMTK
using Compat.Test
using Compat.Iterators

for shape in (8, 6, 10)
    chain = Chain(Fixed, shape)

    @test collect(1:shape) == collect(sites(chain))

    for K in (1, 3, 5)
        t = map(collect(1:shape-K)) do x
            x, x+K
        end
        itr = bonds(chain, K)
        @test chain == lattice(itr)
        @test t == collect(itr)
    end
end

for shape in (8, 6, 10)
    chain = Chain(Periodic, shape)

    @test collect(1:shape) == collect(sites(chain))

    for K in (1, 3, 5)
        t = map(collect(1:shape)) do x
            x, (x+K-1)%shape+1
        end
        itr = bonds(chain, K)
        @test chain == lattice(itr)
        @test t == collect(bonds(chain, K))
    end
end


text:
using Random

using StaticArrays:SizedVector

struct RandomAgent{O,A,E <: AbstractEnvironment{O,A}} <: AbstractAgent{O,A}
    env::E
end

# observe_first!(::RandomAgent) = nothing
function select_action!(agent::RandomAgent, args...; kargs...)
    SizedVector{action_length(agent.env)}(agent.env.env.action_space.sample())
end
observe!(::RandomAgent, observation) = nothing
observe!(::RandomAgent, action, reward, observation, terminal) = nothing
update!(::RandomAgent) = nothing


text:
# Fusion trees:
#==============================================================================#
"""
    struct FusionTree{I, N, M, L, T}

Represents a fusion tree of sectors of type `I<:Sector`, fusing (or splitting) `N` uncoupled
sectors to a coupled sector. (It actually represents a splitting tree, but fusion tree
is a more common term). The `isdual` field indicates whether an isomorphism is present
(if the corresponding value is true) or not. The field `uncoupled` contains the sectors
coming out of the splitting trees, before the possible 𝑍 isomorphism. This fusion tree
has `M=max(0, N-2)` inner lines. Furthermore, for `FusionStyle(I) isa GenericFusion`,
the `L=max(0, N-1)` corresponding vertices carry a label of type `T`. If `FusionStyle(I)
isa MultiplicityFreeFusion, `T = Nothing`.
"""
struct FusionTree{I<:Sector, N, M, L, T}
    uncoupled::NTuple{N, I}
    coupled::I
    isdual::NTuple{N, Bool}
    innerlines::NTuple{M, I} # M = N-2
    vertices::NTuple{L, T} # L = N-1
    function FusionTree{I, N, M, L, T}(uncoupled::NTuple{N, I},
                                            coupled::I,
                                            isdual::NTuple{N, Bool},
                                            innerlines::NTuple{M, I},
                                            vertices::NTuple{L, T}) where
                                            {I<:Sector, N, M, L, T}
        # if N == 0
        #     @assert coupled == one(coupled)
        # elseif N == 1
        #     @assert coupled == uncoupled[1]
        # elseif N == 2
        #     @assert coupled ∈ ⊗(uncoupled...)
        # else
        #     @assert innerlines[1] ∈ ⊗(uncoupled[1], uncoupled[2])
        #     for n = 2:N-2
        #         @assert innerlines[n] ∈ ⊗(innerlines[n-1], uncoupled[n+1])
        #     end
        #     @assert coupled ∈ ⊗(innerlines[N-2], uncoupled[N])
        # end
        new{I, N, M, L, T}(uncoupled, coupled, isdual, innerlines, vertices)
    end
end
function FusionTree{I}(uncoupled::NTuple{N, Any},
                        coupled,
                        isdual::NTuple{N, Bool},
                        innerlines,
                        vertices = ntuple(n->nothing, max(0, N-1))
                        ) where {I<:Sector, N}
    if FusionStyle(I) isa GenericFusion
        fusiontreetype(I, N)(map(s->convert(I, s), uncoupled),
            convert(I, coupled), isdual, map(s->convert(I, s), innerlines), vertices)
    else
        vertices′ = ntuple(n->nothing, max(0, N-1))
        if vertices == vertices′ || all(isone, vertices)
            fusiontreetype(I, N)(map(s->convert(I, s), uncoupled),
                convert(I, coupled), isdual, map(s->convert(I, s), innerlines), vertices)
        else
            throw(ArgumentError("Incorrect fusion vertices"))
        end
    end
end
function FusionTree(uncoupled::NTuple{N, I},
                        coupled::I,
                        isdual::NTuple{N, Bool},
                        innerlines,
                        vertices = ntuple(n->nothing, max(0, N-1))
                        ) where {I<:Sector, N}
    if FusionStyle(I) isa GenericFusion
        fusiontreetype(I, N)(uncoupled, coupled, isdual, innerlines, vertices)
    else
        vertices′ = ntuple(n->nothing, max(0, N-1))
        if vertices == vertices′ || all(isone, vertices)
            fusiontreetype(I, N)(uncoupled, coupled, isdual, innerlines, vertices′)
        else
            throw(ArgumentError("Incorrect fusion vertices"))
        end
    end
end

function FusionTree{I}(uncoupled::NTuple{N}, coupled = one(I),
                        isdual = ntuple(n->false, N)) where {I<:Sector, N}
    FusionStyle(I) isa UniqueFusion ||
        error("fusion tree requires inner lines if `FusionStyle(I) <: MultipleFusion`")
    FusionTree{I}(map(s->convert(I, s), uncoupled), convert(I, coupled), isdual,
                    _abelianinner(map(s->convert(I, s), (uncoupled..., dual(coupled)))))
end
function FusionTree(uncoupled::NTuple{N, I}, coupled::I = one(I),
                        isdual = ntuple(n->false, N)) where {I<:Sector, N}
    FusionStyle(I) isa UniqueFusion ||
        error("fusion tree requires inner lines if `FusionStyle(I) <: MultipleFusion`")
    FusionTree{I}(uncoupled, coupled, isdual, _abelianinner((uncoupled..., dual(coupled))))
end

# Properties
sectortype(::Type{<:FusionTree{I}}) where {I<:Sector} = I
FusionStyle(::Type{<:FusionTree{I}}) where {I<:Sector} = FusionStyle(I)
BraidingStyle(::Type{<:FusionTree{I}}) where {I<:Sector} = BraidingStyle(I)
Base.length(::Type{<:FusionTree{<:Sector, N}}) where {N} = N

sectortype(f::FusionTree) = sectortype(typeof(f))
FusionStyle(f::FusionTree) = FusionStyle(typeof(f))
BraidingStyle(f::FusionTree) = BraidingStyle(typeof(f))
Base.length(f::FusionTree) = length(typeof(f))

# Hashing, important for using fusion trees as key in a dictionary
function Base.hash(f::FusionTree{I}, h::UInt) where {I}
    h = hash(f.isdual, hash(f.coupled, hash(f.uncoupled, h)))
    if FusionStyle(I) isa MultipleFusion
        h = hash(f.innerlines, h)
    end
    if FusionStyle(I) isa GenericFusion
        h = hash(f.vertices, h)
    end
    return h
end
function Base.isequal(f1::FusionTree{I, N}, f2::FusionTree{I, N}) where {I<:Sector, N}
    f1.coupled == f2.coupled || return false
    @inbounds for i = 1:N
        f1.uncoupled[i] == f2.uncoupled[i] || return false
        f1.isdual[i] == f2.isdual[i] || return false
    end
    if FusionStyle(I) isa MultipleFusion
        @inbounds for i=1:N-2
            f1.innerlines[i] == f2.innerlines[i] || return false
        end
    end
    if FusionStyle(I) isa GenericFusion
        @inbounds for i=1:N-1
            f1.vertices[i] == f2.vertices[i] || return false
        end
    end
    return true
end
Base.isequal(f1::FusionTree, f2::FusionTree) = false

# Facilitate getting correct fusion tree types
function fusiontreetype(::Type{I}, N::Int) where {I<:Sector}
    if N === 0
        FusionTree{I, 0, 0, 0, vertex_labeltype(I)}
    elseif N === 1
        FusionTree{I, 1, 0, 0, vertex_labeltype(I)}
    else
        FusionTree{I, N, N-2, N-1, vertex_labeltype(I)}
    end
end

# converting to actual array
function Base.convert(A::Type{<:AbstractArray}, f::FusionTree{I, 0}) where {I}
    X = convert(A, fusiontensor(one(I), one(I), one(I)))[1, 1, :]
    return X
end
function Base.convert(A::Type{<:AbstractArray}, f::FusionTree{I, 1}) where {I}
    c = f.coupled
    if f.isdual[1]
        sqrtdc = sqrtdim(c)
        Zcbartranspose = sqrtdc * convert(A, fusiontensor(conj(c), c, one(c)))[:, :, 1, 1]
        X = conj!(Zcbartranspose) # we want Zcbar^†
    else
        X = convert(A, fusiontensor(c, one(c), c))[:, 1, :, 1, 1]
    end
    return X
end

function Base.convert(A::Type{<:AbstractArray}, f::FusionTree{I, 2}) where {I}
    a, b = f.uncoupled
    isduala, isdualb = f.isdual
    c = f.coupled
    da, db, dc = dim.((a, b, c))
    μ = (FusionStyle(I) isa GenericFusion) ? f.vertices[1] : 1
    C = convert(A, fusiontensor(a, b, c))[:, :, :, μ]
    X = C
    if isduala
        Xtemp = X
        X = similar(Xtemp)
        Za = convert(A, FusionTree((a,), a, (isduala,), ()))
        TO.contract!(1, Za, :N, Xtemp, :N, 0, X, (1,), (2,), (2, 3), (1,), (1,2,3))
    end
    if isdualb
        Xtemp = X
        X = similar(Xtemp)
        Zb = convert(A, FusionTree((b,), b, (isdualb,), ()))
        TO.contract!(1, Zb, :N, Xtemp, :N, 0, X, (1,), (2,), (1, 3), (2,), (2,1,3))
    end
    return X
end

function Base.convert(A::Type{<:AbstractArray}, f::FusionTree{I,N}) where {I,N}
    tailout = (f.innerlines[1], TupleTools.tail2(f.uncoupled)...)
    isdualout = (false, TupleTools.tail2(f.isdual)...)
    ftail = FusionTree(tailout, f.coupled, isdualout,
                        Base.tail(f.innerlines), Base.tail(f.vertices))
    Ctail = convert(A, ftail)
    f1 = FusionTree((f.uncoupled[1], f.uncoupled[2]), f.innerlines[1],
                    (f.isdual[1], f.isdual[2]), (), (f.vertices[1],))
    C1 = convert(A, f1)
    dtail = size(Ctail)
    d1 = size(C1)
    X = similar(C1, (d1[1], d1[2], Base.tail(dtail)...))
    trivialtuple = ntuple(identity, Val(N))
    TO.contract!(1, C1, :N, Ctail, :N, 0, X,
                    (1,2), (3,), Base.tail(trivialtuple), (1,), (trivialtuple..., N+1))
    return X
end

# Show methods
function Base.show(io::IO, t::FusionTree{I, N, M, K, Nothing}) where {I<:Sector, N, M, K}
    print(IOContext(io, :typeinfo => I), "FusionTree{", type_repr(I), "}(",
        t.uncoupled, ", ", t.coupled, ", ", t.isdual, ", ", t.innerlines, ")")
end
function Base.show(io::IO, t::FusionTree{I}) where {I<:Sector}
    print(IOContext(io, :typeinfo => I), "FusionTree{", type_repr(I), "}(",
        t.uncoupled, ", ", t.coupled, ", ", t.isdual, ",",
        t.innerlines, ", ", t.vertices, ")")
end

# Manipulate fusion trees
include("manipulations.jl")

# Fusion tree iterators
include("iterator.jl")

# auxiliary routines
# _abelianinner: generate the inner indices for given outer indices in the abelian case
_abelianinner(outer::Tuple{}) = ()
_abelianinner(outer::Tuple{I}) where {I<:Sector} =
    outer[1] == one(I) ? () : throw(SectorMismatch())
_abelianinner(outer::Tuple{I, I}) where {I<:Sector} =
    outer[1] == dual(outer[2]) ? () : throw(SectorMismatch())
_abelianinner(outer::Tuple{I, I, I}) where {I<:Sector} =
    first(⊗(outer...)) == one(I) ? () : throw(SectorMismatch())
function _abelianinner(outer::NTuple{N, I}) where {I<:Sector, N}
    c = first(outer[1] ⊗ outer[2])
    return (c, _abelianinner((c, TupleTools.tail2(outer)...))...)
end


text:
export read_settings

read_settings(x) = as_namedtuple(x)
read_settings(x::String) = as_namedtuple(TOML.parsefile(x))
as_namedtuple(xs) = xs
function as_namedtuple(xs::Dict{<:AbstractString,<:Any})
  kt = Tuple((Symbol(x) for x in keys(xs)))
  vt = Tuple(values(xs))
  NamedTuple{kt}(as_namedtuple.(vt))
end

read_params(x) = x

function read_params(x::DataFrameRow)
  Dict(k => x[k][1] for k in names(x))
end

function read_params(x::DataFrame)
  @assert size(x,1) == 1 "Only a single row of parameters can be provided."
  Dict(k => x[!,k][1] for k in names(x))
end

function checkunits(params::Dict,unit,key)
  try
    uconvert(unit,params[key])
  catch e
    error("parameter $key with value $(params[key]) incompatible with unit $unit.")
  end
end

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


text:
export segy_read

"""
block = segy_read(file::String)
"""
function segy_read(file::String; buffer::Bool = true, warn_user::Bool = true)

    if buffer
        s = IOBuffer(read(open(file)))
    else
        s = open(file)
    end

    read_file(s, warn_user)
end

"""
block = segy_read(file::String, keys::Array{String,1})
"""
function segy_read(file::String, keys::Array{String,1}; buffer::Bool = true, warn_user::Bool = true)

    if buffer
        s = IOBuffer(read(open(file)))
    else
        s = open(file)
    end

    read_file(s, keys, warn_user)
end


text:
include("testenv.jl")
addprocs_with_testenv(2)

using JuliaDB
using Base.Test
using TextParse
using IndexedTables
using NamedTuples
using PooledArrays

include("test_iteration.jl")
include("test_util.jl")
include("test_table.jl")
include("test_query.jl")
include("test_join.jl")
include("test_misc.jl")
include("test_rechunk.jl")
include("test_readwrite.jl")
include("test_onlinestats.jl")
include("test_ml.jl")
include("test_weakrefstrings.jl")


text:
# ---
# title: 1584. Min Cost to Connect All Points
# id: problem1584
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Union Find
# link: <https://leetcode.com/problems/min-cost-to-connect-all-points/description/>
# hidden: true
# ---
#
# You are given an array `points` representing integer coordinates of some
# points on a 2D-plane, where `points[i] = [xi, yi]`.
#
# The cost of connecting two points `[xi, yi]` and `[xj, yj]` is the **manhattan
# distance** between them: `|xi - xj| + |yi - yj|`, where `|val|` denotes the
# absolute value of `val`.
#
# Return  _the minimum cost to make all points connected._ All points are
# connected if there is **exactly one** simple path between any two points.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/26/d.png)
#
#
#
#     Input: points = [[0,0],[2,2],[3,10],[5,2],[7,0]]
#     Output: 20
#     Explanation:![](https://assets.leetcode.com/uploads/2020/08/26/c.png)
#     We can connect the points as shown above to get the minimum cost of 20.
#     Notice that there is a unique path between every pair of points.
#
#
# **Example 2:**
#
#
#
#     Input: points = [[3,12],[-2,5],[-4,1]]
#     Output: 18
#
#
# **Example 3:**
#
#
#
#     Input: points = [[0,0],[1,1],[1,0],[-1,1]]
#     Output: 4
#
#
# **Example 4:**
#
#
#
#     Input: points = [[-1000000,-1000000],[1000000,1000000]]
#     Output: 4000000
#
#
# **Example 5:**
#
#
#
#     Input: points = [[0,0]]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= points.length <= 1000`
#   * `-106 <= xi, yi <= 106`
#   * All pairs `(xi, yi)` are distinct.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Pkg; Pkg.activate("..")

using OceanTurb, Printf, StaticArrays, BenchmarkTools

Ns = (32,)

function atomic_forward_euler_benchmark(m, Δt)
    @printf "% 24s:" "prepare"
    @btime OceanTurb.prepare!($(m.bcs), $(m.timestepper.eqn.K), $(m.solution), $(m.timestepper.eqn.update!), $(m.grid.N), $m)

    @printf "% 24s:" "calc explicit rhs"
    @btime OceanTurb.calc_explicit_rhs!($(m.timestepper.rhs), $(m.timestepper.eqn), $(m.solution), $m)

    @printf "% 24s:" "update solution"
    @btime OceanTurb.forward_euler_update!($(m.timestepper.rhs), $(m.solution), $Δt)

    @printf "% 24s:" "tick"
    @btime OceanTurb.tick!($(m.clock), $Δt)
    return nothing
end

function atomic_backward_euler_benchmark(m, Δt)
    @printf "% 24s:" "prepare"
    @btime OceanTurb.prepare!($(m.bcs), $(m.timestepper.eqn.K), $(m.solution), $(m.timestepper.eqn.update!), $(m.grid.N), $m)

    @printf "% 24s:" "calc implicit rhs"
    @btime OceanTurb.calc_implicit_rhs!($(m.timestepper.rhs), $(m.timestepper.eqn), $(m.solution), $m)

    @printf "% 24s:" "calc diffusive lhs"
    @btime OceanTurb.calc_diffusive_lhs!($Δt, $(m.timestepper.lhs), $(m.timestepper.eqn.K), $(m.solution), $m)

    @printf "% 24s:" "update solution"
    @btime OceanTurb.backward_euler_update!($(m.timestepper.rhs), $(m.timestepper.lhs), $(m.solution), $Δt)

    @printf "% 24s:" "tick"
    @btime OceanTurb.tick!($(m.clock), $Δt)
    return nothing
end

@printf "\nAtomic timestepping benchmarking of default configuration...\n"
Δt = 1.1
for N in Ns

    @printf "N: % 4d, forward Euler\n" N
    m = Diffusion.Model(N=N, stepper=:ForwardEuler)
    atomic_forward_euler_benchmark(m, Δt)

    @printf "N: % 4d, backward Euler\n" N
    m = Diffusion.Model(N=N, stepper=:BackwardEuler)
    atomic_backward_euler_benchmark(m, Δt)
end

@printf "\nAtomic timestepping benchmarking of configuration with one constant boundary condition...\n"
Δt = 1.1
for N in Ns

    @printf "N: % 4d, forward Euler\n" N
    m = Diffusion.Model(N=N, stepper=:ForwardEuler)
    m.bcs.c.top = BoundaryCondition(Gradient, 0)
    atomic_forward_euler_benchmark(m, Δt)

    @printf "N: % 4d, backward Euler\n" N
    m = Diffusion.Model(N=N, stepper=:BackwardEuler)
    m.bcs.c.top = BoundaryCondition(Gradient, 0)
    atomic_backward_euler_benchmark(m, Δt)
end

@printf "\nAtomic timestepping benchmarking of configuration with one function boundary condition...\n"
Δt = 1.1
bc(model) = time(model)
for N in Ns

    @printf "N: % 4d, forward Euler\n" N
    m = Diffusion.Model(N=N, stepper=:ForwardEuler)
    m.bcs.c.top = BoundaryCondition(Gradient, bc)
    atomic_forward_euler_benchmark(m, Δt)

    @printf "N: % 4d, backward Euler\n" N
    m = Diffusion.Model(N=N, stepper=:BackwardEuler)
    m.bcs.c.top = BoundaryCondition(Gradient, bc)
    atomic_backward_euler_benchmark(m, Δt)
end


text:
using FEASTSolver
using LinearAlgebra
using Profile
using Random

function bf()
    A = diagm(-1 => fill(-1.0, 99), 0 => fill(2, 100), 1 => fill(-1.0, 99))
    # Random.seed!(1234)
    # A = randn(100,100)
    # A .= A + A'
    function linprob()
	f(z) = I*z - A
    end

    return linprob()
end

C, R, = complex(0.0, 0.0), 10.0/1033.660731699839
T = bf()
Y = ComplexF64.(rand(Float64,100,10))
e, v, res = nlfeast!(T, rand(ComplexF64,100,10), 2^3, 0, c=C, r=R)
e, v, res = @timev nlfeast!(T, Y, 2^2, 10, c=C, r=R, debug=true, ϵ=10e-8)
# e, v, res = @profile nlfeast!(T, rand(ComplexF64,64,30), 2^4, 20, c=complex(1.0,1.0), r=0.5, debug=true, ϵ=10e-16)
# Profile.print()
# e, v, res = nlfeast_it!(T, rand(ComplexF64,64,30), 2^3, 0, c=complex(1.0,1.0), r=0.5)
# e, v, res = @timev nlfeast_it!(T, rand(ComplexF64,64,30), 2^3, 20, c=complex(1.0,1.0), r=0.5, debug=true, ϵ=10e-14)
# display(e)
inside(x) = in_contour(x, C, R)
in_eig = e[inside.(e)]
in_res = res[inside.(e)]
in_res_conv = in_res[in_res .<= 1e-6]
in_eig_conv = in_eig[in_res .<= 1e-6]
print("\nmax res inside: ")
println(maximum(res[inside.(e)]))
print("max res inside non spurious: ")
println(maximum(in_res_conv))
print("number inside : ")
println(size(e[inside.(e)])[1])
print("number inside converged : ")
println(size(in_eig_conv)[1])
display(in_eig_conv)


text:
function patdatainterp(th::Float64,filename)
    # load patient data
    file = MAT.matread(filename);
    Qdata = file["Q"]; # Q should be in m3/s
    tdata = file["tgood"]; # times should be in s
    # Pdata = file["P"]; #  pressure should be in Pa
    Vdata = file["V"]; # volume should be in m3
    thdata = file["th"];
    # vectorize data
    Qdata = vec(Qdata);
    tdata = vec(tdata);
    # Pdata = vec(Pdata);
    Vdata = vec(Vdata);
    # create interpolation objects
    Qditp = Interpolations.interpolate(Qdata, Interpolations.BSpline(Interpolations.Linear()), Interpolations.OnGrid());
    # Pditp = Interpolations.interpolate(Pdata, Interpolations.BSpline(Interpolations.Linear()), Interpolations.OnGrid());
    Vditp = Interpolations.interpolate(Vdata, Interpolations.BSpline(Interpolations.Linear()), Interpolations.OnGrid());
    # perform interpolations
    t = linspace(0,th,length(Qdata));
    tq = tdata*th/thdata;
    sitp = Interpolations.scale(Qditp,t);
    qq = [sitp[j] for j in tq];
    # sitp = Interpolations.scale(Pditp,t);
    # pq = [sitp[j] for j in tq];
    sitp = Interpolations.scale(Vditp,t);
    vq = [sitp[j] for j in tq];
    return tq,qq,vq
end


text:
using Clang # this requires Clang 0.14+
using Clang.Generators
using libsixel_jll

cd(@__DIR__)

include_dir = joinpath(libsixel_jll.artifact_dir, "include") |> normpath

options = load_options(joinpath(@__DIR__, "generator.toml"))

# add compiler flags, e.g. "-DXXXXXXXXX"
args = get_default_args()
push!(args, "-I$include_dir")

headers = [joinpath(include_dir, header) for header in readdir(include_dir) if endswith(header, ".h")]
# there is also an experimental `detect_headers` function for auto-detecting top-level headers in the directory
# headers = detect_headers(include_dir, args)

# create context
ctx = create_context(headers, args, options)

# run generator
build!(ctx)


text:
@testset "Neighborsearch" begin
  @testset "NeighborhoodSearch" begin
    𝒟 = RegularGrid(10,10)

    S = NeighborhoodSearch(𝒟, BallNeighborhood(1.0))
    n = GeoStatsBase.search([0.,0.], S)
    @test Set(n) == Set([1,2,11])
    n = GeoStatsBase.search([9.,0.], S)
    @test Set(n) == Set([9,10,20])
    n = GeoStatsBase.search([0.,9.], S)
    @test Set(n) == Set([91,81,92])
    n = GeoStatsBase.search([9.,9.], S)
    @test Set(n) == Set([100,99,90])

    S = NeighborhoodSearch(𝒟, BallNeighborhood(√2))
    n = GeoStatsBase.search([0.,0.], S)
    @test Set(n) == Set([1,2,11,12])
    n = GeoStatsBase.search([9.,0.], S)
    @test Set(n) == Set([9,10,19,20])
    n = GeoStatsBase.search([0.,9.], S)
    @test Set(n) == Set([81,82,91,92])
    n = GeoStatsBase.search([9.,9.], S)
    @test Set(n) == Set([89,90,99,100])

# Non-MinkowskiMetric example
    𝒟 = RegularGrid((360, 180), (0., -90.), (1., 1.))
    S = NeighborhoodSearch(𝒟, BallNeighborhood(150., Haversine(6371.0)))
    n = GeoStatsBase.search([0.,0.], S)
    @test Set(n) == Set([32041, 32402, 32401, 32761, 32760])
  end

  @testset "KNearestSearch" begin
    𝒟 = RegularGrid(10,10)
    S = KNearestSearch(𝒟, 3)
    n = GeoStatsBase.search([0.,0.], S)
    @test Set(n) == Set([1,2,11])
    n = GeoStatsBase.search([9.,0.], S)
    @test Set(n) == Set([9,10,20])
    n = GeoStatsBase.search([0.,9.], S)
    @test Set(n) == Set([91,81,92])
    n = GeoStatsBase.search([9.,9.], S)
    @test Set(n) == Set([100,99,90])
  end

  @testset "KBallSearch" begin
    #create a regular grid to test
    𝒟 = RegularGrid(10,10)

    #test 1: all points should be found
    k = 10
    ball = BallNeighborhood(100.0)
    searcher = KBallSearch(𝒟, k, ball)

    xₒ = [5.0,5.0]
    neighbors = GeoStatsBase.search(xₒ, searcher)

    @test length(neighbors) == 10

    #test 2: radius of 1.0 m 5 points should be found
    k = 10
    ball = BallNeighborhood(1.0)
    searcher = KBallSearch(𝒟, k, ball)

    xₒ = [5.0,5.0]
    neighbors = GeoStatsBase.search(xₒ, searcher)
    @test length(neighbors) == 5
    @test neighbors[1] == 56

    #testing with masks
    mask = trues(nelms(𝒟))
    mask[56] = false #excluding location #56
    neighbors = GeoStatsBase.search(xₒ, searcher, mask=mask)
    @test length(neighbors) == 4

    #test 3: radius of 1.0 m but far just to find 1 point
    k = 10
    ball = BallNeighborhood(1.0)
    searcher = KBallSearch(𝒟, k, ball)

    xₒ = [-0.2,-0.2]
    neighbors = GeoStatsBase.search(xₒ, searcher)

    @test length(neighbors) == 1

    #test 4: radius of 1.0 m but far just to find no points
    k = 10
    ball = BallNeighborhood(1.0)
    searcher = KBallSearch(𝒟, k, ball)

    xₒ = [-10.0,-10.0]
    neighbors = GeoStatsBase.search(xₒ, searcher)

    @test length(neighbors) == 0
  end

  @testset "BoundedSearch" begin
    𝒟 = RegularGrid(10,10)
    S1 = NeighborhoodSearch(𝒟, BallNeighborhood(5.0))
    S2 = KNearestSearch(𝒟, 10)
    B1 = BoundedSearch(S1, 5)
    B2 = BoundedSearch(S2, 5)
    n = GeoStatsBase.search(coordinates(𝒟, rand(1:100)), B1)
    @test length(n) == 5
    n = GeoStatsBase.search(coordinates(𝒟, rand(1:100)), B2)
    @test length(n) == 5
  end
end