text:
module Pykspa
using JSON, Dates

function get_sld(length, seed)
    sld = ""
    modulo = 541 * length + 4
    a = length * length
    for i in 1:length
        index = (a + (seed*((seed % 5) + (seed % 123456) + i*((seed & 1) + (seed % 4567))) & 0xFFFFFFFF)) % 26
        a += length
        a &= 0xFFFFFFFFF
        sld *= 'a'+index
        seed += (((7837632 * seed * length) & 0xFFFFFFFF) + 82344) % modulo
    end
    return sld
end

function gen_one(nr=missing, date=missing, set_nr=2, add_tld=false)
    if date === missing
        date = Dates.Date(Dates.now())
    end
    if nr === missing
        nr = rand(Int8)
    end
    seeds = JSON.parsefile(string(@__DIR__) * "/../data/set$(set_nr)_seeds.json")
    dt = Dates.datetime2unix(Dates.DateTime(date))
    days = set_nr==1 ? 20 : 1
    index = trunc(Int64, dt/(days*3600*24))
    if !(string(index) in keys(seeds))
        return domain
    end
    seed = parse(Int64, seeds[string(index)], base=16)
    orig_seed = seed
    s = seed % (nr + 1)
    seed += (s+1)
    len = (seed + nr)%7 + 6
    second_level_domain = get_sld(len,seed)
    if add_tld
        tlds = ["com", "net", "org", "info", "cc"]
            top_level_domain = tlds[(seed & 3)+1]
            domain = second_level_domain * '.' * top_level_domain
    else
        domain = second_level_domain
    end
    return domain
end

function gen_many(cnt=100)
    res = []
    date = Dates.Date(Dates.now())
    seeds = JSON.parsefile(string(@__DIR__) * "/../data/set2_seeds.json")
    dt = Dates.datetime2unix(Dates.DateTime(date))
    days = 1
    index = trunc(Int64, dt/(days*3600*24))
    if !(string(index) in keys(seeds))
        return domain
    end
    seed = parse(Int64, seeds[string(index)], base=16)
    orig_seed = seed

    for dga_nr in 1:cnt
        s = seed % (dga_nr + 1)
        seed += (s+1)
        len = (seed + dga_nr)%7 + 6
        second_level_domain = get_sld(len,seed)
        domain = second_level_domain
        res = [res..., domain]
    end

    return res
end

end  # module Pykspa


text:
module chapters_getting_started

using Dash, DashHtmlComponents, DashCoreComponents
include("../../utils.jl")
include("../../styles.jl")

include("../../reusable_components/Header.jl")
include("../../reusable_components/Example.jl")

getting_started_layout_1 = LoadExampleCode("./dash_docs/chapters/getting_started/examples/getting_started_layout_1.jl")

getting_started_layout_2 = LoadExampleCode("./dash_docs/chapters/getting_started/examples/getting_started_layout_2.jl")

getting_started_layout_3 = LoadExampleCode("./dash_docs/chapters/getting_started/examples/getting_started_layout_3.jl")

getting_started_layout_4 = LoadExampleCode("./dash_docs/chapters/getting_started/examples/getting_started_layout_4.jl")

getting_started_layout_5 = LoadExampleCode("./dash_docs/chapters/getting_started/examples/getting_started_layout_5.jl")

getting_started_layout_6 = LoadExampleCode("./dash_docs/chapters/getting_started/examples/getting_started_layout_6.jl")

app = dash()

app.layout = html_div() do
    Header("Dash Layout"),
    html_blockquote(dcc_markdown("This is the 2nd chapter of the [Dash Tutorial](/).
    The previous chapter covered [installation](/installation) and the next chapter covers [Dash callbacks](/basic-callbacks).")),
    dcc_markdown("""
    This tutorial will walk you through a fundamental aspect of Dash apps, the app `layout`, through several self-contained apps.
    """),
    html_hr(),
    dcc_markdown("""
    Dash apps are composed of two parts. The first part is the `layout` of the app and it describes what the application looks like.
    The second part describes the interactivity of the application and will be covered in the [next chapter](/basic-callbacks).

    Dash provides Julia methods for all of the visual components of the application.
    We maintain a set of components in the `DashCoreComponents` and the `DashHtmlComponents` packages
    but you can also build your own with JavaScript and React.js

    Note: Throughout this documentation, Julia code examples are meant to be saved as files and executed using `julia app.jl`.
    These examples are not intended to run in Jupyter notebooks as-is,
    although most can be modified slightly to function in that environment.

    To get started, create a file called `app.jl` with the following code:
    """),
    getting_started_layout_1.source_code,
    Example(getting_started_layout_1.layout),

    dcc_markdown("""

    Run the app with
    ```
    \$ julia app.jl
    ```

    In the `julia` REPL you can run the following code, assuming that `app.jl` is in your current working directory.
    ```
    julia> include("app.jl")
    ```

    You can visit the app by visiting the URL http://127.0.0.1:8050.

    Note:

    1. The `layout` is composed of a tree of "components" like `html_div` and `dcc_graph`.
    2. The `DashHtmlComponents` package has a component for every HTML tag. The `html_h1("Hello Dash")` component
    generates a `<h1>Hello Dash</h1>` HTML element in your application.
    3. Not all components are pure HTML. The `DashCoreComponents` package describes higher level components
    that are interactive and genrated with JavaScript, HTML, and CSS through the React.js library.
    4. Each component is described entirely through keyword attributes. Dash is *declarative*: you will
    primarily describe your application through these attributes.
    5. The `children` property is special. By convention, it's always the first attribute which means that you can omit it;
    `html_div(children="Hello Dash")` is the same as `html_div("Hello Dash")`. Also, it can contain a string, a number, a single
    component, or a list of components.
    6. The fonts in your application will look a little bit different than what is displayed here. This application is using
    a custom CSS stylesheet to modify the default styles of the elements. If you would like to use this stylesheet,
    you can initialize your app with

    `app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])`
    """),
    Header("More About HTML"),
    dcc_markdown("""
    The `DashHtmlComponents` package contains a component class for every HTML tag as well as keyword arguments for all of the
    HTML arguments.

    Let's customize the text in our app by modifying the inline styles of the components:
    """),

    getting_started_layout_2.source_code,
    Example(getting_started_layout_2.layout),

    dcc_markdown("""

    In this example, we modified the inline styles of the `html_div` and
    `html_h1` components with the style property. `html_h1("Hello Dash",
    style = Dict("color" => "#7FDBFF", "textAlign" => "center"))` is rendered
    in the Dash application as `<h1 style="text-align: center; color: rgb(127, 219, 255);">Hello Dash</h1>`

    There are a few important differences between the `DashHtmlComponents` and HTML attributes:

    1. The `style` property in HTML is a semicolon separated string. In Dash, you can just supply a `Dict`.
    2. The keys in the `style` `Dict` are [camelCased]("https://en.wikipedia.org/wiki/Camel_case"). So instead of `text-align`, it's `textAlign`.
    3. The HTML `class` attribute is `className` in Dash.
    4. The children of the HTML tag is specified through the `children` keyword argument. By convention, this is
    always the first argument and so it is often omitted. Besides that, all of the available HTML attributes and
    tags are available to you within your Julia context.
    """),
    Header("Reusable Components"),
    dcc_markdown("""
    By writing our markup in Julia, we can create
    complex reusable components like tables without switching contexts or
    languages.

    Here's a quick example that generates a `Table` from a `DataFrame`.
    Create a file named `app.jl` with the following code:
    """),
    getting_started_layout_3.source_code,
    Example(getting_started_layout_3.layout),
    Header("More About Visualization"),
    dcc_markdown("""
    The `DashCoreComponents` package includes a component called `graph`.

    `graph` renders interactive data visualizations using the open source
    plotly.js JavaScript graphing library. Plotly.js supports over 35 chart types
    and renders charts in both vector-quality SVG and high-performance WebGL.

    The `figure` argument in the `dcc_graph` component is the same `figure` argument
    that is used by `plotly.py`, Plotly's open-source Python graphing library. Check out the
    plotly.py documentation and gallery to learn more.

    Here's an example that creates a scatter plot from a `DataFrame`. Create a file named `app.jl`
    with the following code:
    """),
    getting_started_layout_4.source_code,
    Example(getting_started_layout_4.layout),
    dcc_markdown("""

    These graphs are interactive and responsive. *Hover* over points to see their values, *click* on legend items to
    toggle traces, *click and drag* to zoom, *hold down shift and click and drag* to pan.

    """),
    Header("Markdown"),
    dcc_markdown("""
    While Dash exposes HTML through the `DashHtmlComponents` package, it can be tedious to write your
    copy in HTML. For writing blocks of text, you can use the `dcc_markdown` component in the `DashCoreComponents` package.
    """),
    getting_started_layout_5.source_code,
    Example(getting_started_layout_5.layout),
    Header("Core Components"),
    dcc_markdown("""
    The `DashCoreComponents` package includes a set of higher level components like
    dropdowns, graphs, markdown blocks, and more. Like all Dash components, they are described
    entirely declaratively. Every option that is configurable is available as a keyword
    argument to the component.
    """),
    getting_started_layout_6.source_code,
    Example(getting_started_layout_6.layout),
    dcc_markdown("""
    The next part of the Dash tutorial covers how to make these apps interactive.

    [Dash Tutorial Part 3: Basic Callbacks](/basic-callbacks)
    """),

    html_hr(),
    dcc_markdown("[Back to the table of contents](/)")

end

end

ext: jl

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
using helloIBox
using Sockets
using HTTP
using Base64
using JSON
using Test

localIP =  Sockets.localhost
port = 8080
pathToIBox = "E:\\box_fail\\IBox\\IBoxLauncher.exe"
pathToIBox = "C:/Temp/IBox/IBoxLauncher.exe"
# pathToIBox = "Y:\\Yuly\\IBox\\IBoxLauncher.exe"
@testset "Access to data through IBox" begin
server_path = "E:/box_fail" # "C:/Users/yzh/Downloads/TestServer"
server_path = "C:/Users/yzh/Downloads/TestServer"

start_server(server_path; localIP = localIP, port = port)
r = HTTP.request("GET", "http://$localIP:$port/api/runIBox?res=oxy115829.dat&IBox_port=8888&IBox_path=$pathToIBox&IBox_host=$localIP")
#ЗАПРОС ДАННЫХ
r = HTTP.request("GET", "http://$localIP:$port/api/getData?res=oxy115829.dat&dataName=Freq&index=0&from=0&count=1")
Freq = reinterpret(Int32, base64decode(r.body)) |> collect
@test Freq[1]==257

r = HTTP.request("GET", "http://$localIP:$port/api/getData?res=oxy115829.dat&dataName=Ecg_1&from=0&count=20")
ecg1 = reinterpret(Int32, base64decode(r.body)) |> collect
@test length(ecg1)==20 && ecg1[1]==3099 && ecg1[20]==3062

# r = HTTP.request("GET", "http://$localIP:$port/api/getStructData?res=oxy115829.dat&dataName=/Mark/QRS&from=5847867&to=5852730&fields=QPoint,WidthQRS")

r = HTTP.request("GET", "http://$localIP:$port/api/getData?dataName=Ecg_1&from=10&to=15")
ECG = reinterpret(Int32, base64decode(r.body)) |> collect
@test ECG == [3110, 3102, 3072, 3038, 3041]
r = HTTP.request("GET", "http://$localIP:$port/api/getStructData?dataName=QRS&fields=QPoint&from=100&to=400")
QPoint = reinterpret(Int32, base64decode(r.body)) |> collect
@test QPoint == [101, 221]

r = HTTP.request("GET", "http://$localIP:$port/api/getDataTree")
tree = JSON.parse(String(r.body))
@test tree[1]["nodes"][1]["nodes"][1]["name"] == "Ecg"

r = HTTP.request("GET", "http://$localIP:$port/api/getType?dataName=QPoint")
@test String(r.body)=="Int32"
r = HTTP.request("GET", "http://$localIP:$port/api/getData?dataName=Freq&index=0&all")
Amp = reinterpret(Int32, base64decode(r.body)) |> collect

r = HTTP.request("GET", "http://$localIP:$port/api/getStructData?dataName=/Mark/QRS&from=110684&to=110694&fields=QPoint,WidthQRS")
helloIBox.unpack_vec(base64decode(r.body), Int32, Int16) == [(Int32(110629), Int16(24))]

r = HTTP.request("GET", "http://$localIP:$port/api/getStructData?res=oxy115829.dat&dataName=/Mark/QRS&from=100&to=300&fields=QPoint,WidthQRS")
QPoint =helloIBox.unpack_vec(base64decode(r.body), Int32, Int16)

r = HTTP.request("GET", "http://$localIP:$port/api/getStructData?res=oxy115829.dat&dataName=/Mark/QRS&from=100&to=300&fields=QPoint")
QPoint =helloIBox.unpack_vec(base64decode(r.body), Int32)

r = HTTP.request("GET", "http://$localIP:$port/api/getTag?res=oxy115829.dat&dataName=Ecg&index=2")
@test String(r.body)=="C1"

r = HTTP.request("GET", "http://$localIP:$port/api/getAttributes?res=oxy115829.dat&dataName=/Mark/EcgChannals/Ecg&index=1")
attr = JSON.parse(String(r.body))
@test attr["tag"] == "F" && attr["dstype"] == "signal"

#другой вариант обращения к данным
r = HTTP.request("GET", "http://$localIP:$port/api/getAttributes?res=oxy115829.dat&dataName=Ecg_2")
attr = JSON.parse(String(r.body))
@test attr["tag"] == "C1" && attr["dstype"] == "signal"

r = HTTP.request("GET", "http://$localIP:$port/api/getStructData?dataName=Ecg&index=1&from=100&to=1500")
r = HTTP.request("GET", "http://$localIP:$port/api/getData?dataName=Ecg&index=1&from=100&to=1500")

r = HTTP.request("GET", "http://$localIP:$port/api/Close")

r = HTTP.request("GET", "http://$localIP:$port/api/closeServer")

end

@testset "Manual changes" begin
portIBox = 8888
start_server(""; localIP = localIP, port = port)
r = HTTP.request("GET", "http://$localIP:$port/api/runIBox?res=oxy115829.dat&IBox_port=$portIBox&IBox_path=$pathToIBox&IBox_host=$localIP")

pathTOjson = joinpath(Base.@__DIR__, "files","oxy115829.002","command_FT0.json")
r = HTTP.request("POST", "http://$localIP:$port/api/manualChange?res=oxy115829.dat", ["Content-Type" => "application/json"], read(pathTOjson))

r = HTTP.request("GET", "http://$localIP:$port/api/getData?dataName=QRS&fields=QPoint&index=1&from=1&to=3")
QPoint = reinterpret(Int32, base64decode(r.body)) |> collect
@test QPoint == [10, 30, 50]
r = HTTP.request("GET", "http://$localIP:$port/api/getStructData?dataName=QRS&fields=QPoint,WidthQRS&from=1&to=300")
QPoint =helloIBox.unpack_vec(base64decode(r.body), Int32, Int16)

r = HTTP.request("GET", "http://$localIP:$port/api/getStructData?dataName=QRS&fields=QPoint&from=1&to=100")
QPoint = reinterpret(Int32, base64decode(r.body)) |> collect
@test QPoint == [10, 30, 50]

r = HTTP.request("GET", "http://$localIP:$port/api/getData?dataName=Ecg_2&from=0&to=100")
QPoint = reinterpret(Int32, base64decode(r.body)) |> collect

r = HTTP.request("GET", "http://$localIP:$port/api/Close")

r = HTTP.request("GET", "http://$localIP:$port/api/closeServer")

end

@testset "Additional Info" begin
portIBox = 8888
start_server(""; localIP = localIP, port = port)
r = HTTP.request("GET", "http://$localIP:$port/api/runIBox?res=oxy115829.dat&IBox_port=$portIBox&IBox_path=$pathToIBox&IBox_host=$localIP")

r = HTTP.request("GET", "http://$localIP:$port/api/getData?dataName=RRm&all")
RRm = reinterpret(Int32, base64decode(r.body)) |> collect

r = HTTP.request("GET", "http://$localIP:$port/api/getData?dataName=ClassQRS&all")
classQRS = reinterpret(Int32, base64decode(r.body)) |> collect
unClasses = unique(classQRS[findall(classQRS.!=10)])

r = HTTP.request("GET", "http://$localIP:$port/api/getStructData?dataName=QRS&fields=QPoint,WidthQRS,ClassQRS,SubClassQRS&from=1&to=300")
QPoint =helloIBox.unpack_vec(base64decode(r.body), Int32, Int16, Int32, Int16)

r = HTTP.request("GET", "http://$localIP:$port/api/getType?dataName=SubClassQRS")

r = HTTP.request("GET", "http://$localIP:$port/api/getData?dataName=XKoef&all")
XKoefWald = reinterpret(Int16, base64decode(r.body)) |> collect

r = HTTP.request("GET", "http://$localIP:$port/api/getData?dataName=SubClassQRS&all")
subClass = reinterpret(Int8, base64decode(r.body)) |> collect
unique(subClass)

r = HTTP.request("GET", "http://$localIP:$port/api/Close")

r = HTTP.request("GET", "http://$localIP:$port/api/closeServer")

end

@testset "QRS tree" begin

#проверяем дерево классов QRS
server_path = "C:/Users/yzh/Downloads/TestServer"

start_server(server_path; localIP = localIP, port = port)

r = HTTP.request("GET", "http://$localIP:$port/api/runIBox?res=oxy115829.dat&IBox_port=8888&IBox_path=$pathToIBox&IBox_host=$localIP&config=IBOpen&resName=111&arg=-open")
r = HTTP.request("GET", "http://$localIP:$port/api/getQRStree")
res = String(r.body)
treeQRS = JSON.parse(res)
@test treeQRS["size"] == 91328 && treeQRS["nodes"][1]["name"]==10
r = HTTP.request("GET", "http://$localIP:$port/api/Close")

r = HTTP.request("GET", "http://$localIP:$port/api/closeServer")

end


text:
# host array

export MtlArray

mutable struct MtlArray{T,N} <: AbstractGPUArray{T,N}
  buffer::MtlBuffer{T}
  dims::Dims{N}

  dev::MtlDevice
end

device(A::MtlArray) = A.dev

## constructors

# type and dimensionality specified, accepting dims as tuples of Ints
function MtlArray{T,N}(::UndefInitializer, dims::Dims{N}; storage=Shared) where {T,N}
    dev = current_device()
    # Check that requested size is not larger than maximum buffer size allowed
    sizeof(T) * prod(dims) > dev.maxBufferLength && error("Too large of Metal buffer requested of size $(Base.format_bytes(sizeof(T) * prod(dims))) (Max: $(Base.format_bytes(dev.maxBufferLength)))")
    buf = alloc(T, dev, prod(dims); storage=storage)

    obj = MtlArray{T,N}(buf, dims, dev)
    finalizer(obj) do arr
        free(arr.buffer)
    end
    return obj
end

# type and dimensionality specified, accepting dims as series of Ints
MtlArray{T,N}(::UndefInitializer, dims::Integer...) where {T,N} = MtlArray{T,N}(undef, Dims(dims))

# type but not dimensionality specified
MtlArray{T}(::UndefInitializer, dims::Dims{N}) where {T,N} = MtlArray{T,N}(undef, dims)
MtlArray{T}(::UndefInitializer, dims::Integer...) where {T} =
    MtlArray{T}(undef, convert(Tuple{Vararg{Int}}, dims))

# empty vector constructor
MtlArray{T,1}() where {T} = MtlArray{T,1}(undef, 0)

Base.similar(a::MtlArray{T,N}) where {T,N} = MtlArray{T,N}(undef, size(a))
Base.similar(a::MtlArray{T}, dims::Base.Dims{N}) where {T,N} = MtlArray{T,N}(undef, dims)
Base.similar(a::MtlArray, ::Type{T}, dims::Base.Dims{N}) where {T,N} = MtlArray{T,N}(undef, dims)

function Base.copy(a::MtlArray{T,N}) where {T,N}
  b = similar(a)
  @inbounds copyto!(b, a)
end

## array interface

Base.elsize(::Type{<:MtlArray{T}}) where {T} = sizeof(T)

Base.size(x::MtlArray) = x.dims
Base.sizeof(x::MtlArray) = Base.elsize(x) * length(x)

# XXX: we can't actually get a true pointer, it's always a buffer
# with Metal APIs accepting an offset param. should we then disallow this conversion?
Base.pointer(x::MtlArray{T}) where {T} = Base.unsafe_convert(MTL.MTLBuffer, x)
Base.pointer(x::MtlArray, index) = error("Cannot compute a pointer with offset")

## interop with other arrays

@inline function MtlArray{T,N}(xs::AbstractArray{T,N}) where {T,N}
  A = MtlArray{T,N}(undef, size(xs))
  copyto!(A, xs)
  return A
end

MtlArray{T,N}(xs::AbstractArray{S,N}) where {T,N,S} = MtlArray{T,N}(map(T, xs))

# underspecified constructors
MtlArray{T}(xs::AbstractArray{S,N}) where {T,N,S} = MtlArray{T,N}(xs)
(::Type{MtlArray{T,N} where T})(x::AbstractArray{S,N}) where {S,N} = MtlArray{S,N}(x)
MtlArray(A::AbstractArray{T,N}) where {T,N} = MtlArray{T,N}(A)

# idempotency
MtlArray{T,N}(xs::MtlArray{T,N}) where {T,N} = xs

## conversions

Base.convert(::Type{T}, x::T) where T <: MtlArray = x

## interop with C libraries

Base.unsafe_convert(::Type{<:Ptr}, x::MtlArray) =
  throw(ArgumentError("cannot take the host address of a $(typeof(x))"))

Base.unsafe_convert(t::Type{MTL.MTLBuffer}, x::MtlArray) = x.buffer

## interop with GPU arrays

# TODO Figure out global

function Base.convert(::Type{MtlDeviceArray{T,N,AS.Device}}, a::MtlArray{T,N}) where {T,N}
    MtlDeviceArray{T,N,AS.Device}(a.dims, reinterpret(Core.LLVMPtr{T, 1}, pointer(a).handle))
end

Adapt.adapt_storage(::Adaptor, xs::MtlArray{T,N}) where {T,N} =
  convert(MtlDeviceArray{T,N,AS.Device}, xs)

# Adapt.adapt_storage(::Adaptor, xs::MtlArray{T,N}) where {T,N} =
#   convert(Core.LLVMPtr{T,AS.Device}, xs)

function Base.convert(::Type{Core.LLVMPtr{T,AS.Device}}, a::MtlArray{T}) where {T}
    reinterpret(Core.LLVMPtr{T, 1}, a.buffer.handle)
end
## interop with CPU arrays

Base.unsafe_wrap(t::Type{<:Array}, arr::MtlArray, dims; own=false) = unsafe_wrap(t, arr.buffer, dims; own=own)

# We don't convert isbits types in `adapt`, since they are already
# considered GPU-compatible.

Adapt.adapt_storage(::Type{MtlArray}, xs::AbstractArray) =
  isbits(xs) ? xs : convert(MtlArray, xs)

# if an element type is specified, convert to it
Adapt.adapt_storage(::Type{<:MtlArray{T}}, xs::AbstractArray) where {T} =
  isbits(xs) ? xs : convert(MtlArray{T}, xs)

Adapt.adapt_storage(::Type{Array}, xs::MtlArray) = convert(Array, xs)

Base.collect(x::MtlArray{T,N}) where {T,N} = copyto!(Array{T,N}(undef, size(x)), x)

function Base.copyto!(dest::MtlArray{T}, doffs::Integer, src::Array{T}, soffs::Integer,
                      n::Integer) where T
  n==0 && return dest
  @boundscheck checkbounds(dest, doffs)
  @boundscheck checkbounds(dest, doffs+n-1)
  @boundscheck checkbounds(src, soffs)
  @boundscheck checkbounds(src, soffs+n-1)
  unsafe_copyto!(device(dest), dest, doffs, src, soffs, n)
  return dest
end

Base.copyto!(dest::MtlArray{T}, src::Array{T}) where {T} =
    copyto!(dest, 1, src, 1, length(src))

function Base.copyto!(dest::Array{T}, doffs::Integer, src::MtlArray{T}, soffs::Integer,
                      n::Integer) where T
  n==0 && return dest
  @boundscheck checkbounds(dest, doffs)
  @boundscheck checkbounds(dest, doffs+n-1)
  @boundscheck checkbounds(src, soffs)
  @boundscheck checkbounds(src, soffs+n-1)
  unsafe_copyto!(device(src), dest, doffs, src, soffs, n)
  return dest
end

Base.copyto!(dest::Array{T}, src::MtlArray{T}) where {T} =
    copyto!(dest, 1, src, 1, length(src))

function Base.copyto!(dest::MtlArray{T}, doffs::Integer, src::MtlArray{T}, soffs::Integer,
                      n::Integer) where T
  n==0 && return dest
  @boundscheck checkbounds(dest, doffs)
  @boundscheck checkbounds(dest, doffs+n-1)
  @boundscheck checkbounds(src, soffs)
  @boundscheck checkbounds(src, soffs+n-1)
  # TODO: which device to use here?
  if device(dest) == device(src)
    unsafe_copyto!(device(dest), dest, doffs, src, soffs, n)
  else
    error("Copy between different devices not implemented")
  end
  return dest
end

Base.copyto!(dest::MtlArray{T}, src::MtlArray{T}) where {T} =
    copyto!(dest, 1, src, 1, length(src))

function Base.unsafe_copyto!(dev::MtlDevice, dest::MtlArray{T}, doffs, src::Array{T}, soffs, n) where T
  # these copies are implemented using pure memcpy's, not API calls, so aren't ordered.
  synchronize()

  GC.@preserve src dest begin
    unsafe_copyto!(dev, pointer(dest), doffs, pointer(src, soffs), n)
  end
  if Base.isbitsunion(T)
    # copy selector bytes
    error("Not implemented")
  end
  return dest
end

function Base.unsafe_copyto!(dev::MtlDevice, dest::Array{T}, doff, src::MtlArray{T}, soff, n) where T
  # these copies are implemented using pure memcpy's, not API calls, so aren't ordered.
  synchronize()

  GC.@preserve src dest begin
    unsafe_copyto!(dev, pointer(dest, doff), pointer(src), soff, n)
  end
  #GC.@preserve src dest unsafe_copyto!(dev, pointer(dest, doffs), pointer(src, soffs), n)
  if Base.isbitsunion(T)
    # copy selector bytes
    error("Not implemented")
  end
  return dest
end

function Base.unsafe_copyto!(dev::MtlDevice, dest::MtlArray{T}, doffs, src::MtlArray{T}, soffs, n) where T
  # these copies are implemented using pure memcpy's, not API calls, so aren't ordered.
  synchronize()

  GC.@preserve src dest unsafe_copyto!(dev, pointer(dest), doffs, pointer(src), soffs, n)
#  GC.@preserve src dest unsafe_copyto!(dev, pointer(dest, doffs), pointer(src, soffs), n)
  if Base.isbitsunion(T)
    # copy selector bytes
    error("Not implemented")
  end
  return dest
end

## utilities

zeros(T::Type, dims...) = fill!(MtlArray{T}(undef, dims...), 0)
ones(T::Type, dims...) = fill!(MtlArray{T}(undef, dims...), 1)
zeros(dims...) = zeros(Float32, dims...)
ones(dims...) = Mtls(Float32, dims...)
fill(v, dims...) = fill!(MtlArray{typeof(v)}(undef, dims...), v)
fill(v, dims::Dims) = fill!(MtlArray{typeof(v)}(undef, dims...), v)

# optimized implementation of `fill!` for types that are directly supported by fillbuffer
function Base.fill!(A::MtlArray{T}, val) where T <: Union{UInt8,Int8}
  B = convert(T, val)
  unsafe_fill!(device(A), pointer(A), B, length(A))
  A
end

## views

device(a::SubArray) = device(parent(a))

# we don't really want an array, so don't call `adapt(Array, ...)`,
# but just want MtlArray indices to get downloaded back to the CPU.
# this makes sure we preserve array-like containers, like Base.Slice.
struct BackToCPU end
Adapt.adapt_storage(::BackToCPU, xs::MtlArray) = convert(Array, xs)

@inline function Base.view(A::MtlArray, I::Vararg{Any,N}) where {N}
    J = to_indices(A, I)
    @boundscheck begin
        # Base's boundscheck accesses the indices, so make sure they reside on the CPU.
        # this is expensive, but it's a bounds check after all.
        J_cpu = map(j->adapt(BackToCPU(), j), J)
        checkbounds(A, J_cpu...)
    end
    J_gpu = map(j->adapt(MtlArray, j), J)
    Base.unsafe_view(Base._maybe_reshape_parent(A, Base.index_ndims(J_gpu...)), J_gpu...)
end

# pointer conversions
## contiguous
function Base.unsafe_convert(::Type{MTL.MTLBuffer}, V::SubArray{T,N,P,<:Tuple{Vararg{Base.RangeIndex}}}) where {T,N,P}
    return Base.unsafe_convert(MTL.MTLBuffer, parent(V)) +
           Base._memory_offset(V.parent, map(first, V.indices)...)
end

## reshaped
function Base.unsafe_convert(::Type{MTL.MTLBuffer}, V::SubArray{T,N,P,<:Tuple{Vararg{Union{Base.RangeIndex,Base.ReshapedUnitRange}}}}) where {T,N,P}
   return Base.unsafe_convert(MTL.MTLBuffer, parent(V)) +
          (Base.first_index(V)-1)*sizeof(T)
end

## PermutedDimsArray

device(a::Base.PermutedDimsArray) = device(parent(a))

Base.unsafe_convert(::Type{MTL.MTLBuffer}, A::PermutedDimsArray) where {T} =
    Base.unsafe_convert(MTL.MTLBuffer, parent(A))

## reshape

device(a::Base.ReshapedArray) = device(parent(a))

Base.unsafe_convert(::Type{MTL.MTLBuffer}, a::Base.ReshapedArray{T}) where {T} =
  Base.unsafe_convert(MTL.MTLBuffer, parent(a))

## reinterpret

device(a::Base.ReinterpretArray) = device(parent(a))

Base.unsafe_convert(::Type{MTL.MTLBuffer}, a::Base.ReinterpretArray{T,N,S} where N) where {T,S} =
  MTL.MTLBuffer(Base.unsafe_convert(ZePtr{S}, parent(a)))


text:
using Test
using ManifoldsBase
using ManifoldsBase: AbstractNumbers, ℝ, ℂ, ℍ, _unify_number_systems

struct NotImplementedNumbers <: ManifoldsBase.AbstractNumbers end

@testset "Number systems" begin
    @test_throws ErrorException real_dimension(NotImplementedNumbers())

    @test ℝ isa ManifoldsBase.RealNumbers
    @test ManifoldsBase.RealNumbers() === ℝ
    @test real_dimension(ℝ) == 1
    @test repr(ℝ) == "ℝ"

    @test ℂ isa ManifoldsBase.ComplexNumbers
    @test ManifoldsBase.ComplexNumbers() === ℂ
    @test real_dimension(ℂ) == 2
    @test repr(ℂ) == "ℂ"

    @test ℍ isa ManifoldsBase.QuaternionNumbers
    @test ManifoldsBase.QuaternionNumbers() === ℍ
    @test real_dimension(ℍ) == 4
    @test repr(ℍ) == "ℍ"

    @test (@inferred _unify_number_systems(ℝ)) === ℝ
    @test (@inferred _unify_number_systems(ℂ)) === ℂ
    @test (@inferred _unify_number_systems(ℍ)) === ℍ
    @test (@inferred _unify_number_systems(ℝ, ℝ)) === ℝ
    @test (@inferred _unify_number_systems(ℝ, ℂ)) === ℂ
    @test (@inferred _unify_number_systems(ℝ, ℍ)) === ℍ
    @test (@inferred _unify_number_systems(ℂ, ℝ)) === ℂ
    @test (@inferred _unify_number_systems(ℂ, ℂ)) === ℂ
    @test (@inferred _unify_number_systems(ℂ, ℍ)) === ℍ
    @test (@inferred _unify_number_systems(ℍ, ℝ)) === ℍ
    @test (@inferred _unify_number_systems(ℍ, ℂ)) === ℍ
    @test (@inferred _unify_number_systems(ℍ, ℍ)) === ℍ
    @test (@inferred _unify_number_systems(ℝ, ℂ, ℝ)) === ℂ
    @test (@inferred _unify_number_systems(ℝ, ℝ, ℝ)) === ℝ
end


text:
include("sum0.jl")

function  mean0(a)
    # create variable to store ouput
    out = 0.0

    # sum0
    s = sum0(a)

    # length
    l = length(a)

    # calculation
    out = s/l

    # return
    return out
end


text:
struct Cuboid
    x::UnitRange
    y::UnitRange
    z::UnitRange
end

struct RebootStep
    on::Bool
    cuboid::Cuboid
end

mutable struct Reactor
    on_cuboids::Vector{Cuboid}
    Reactor() = new(Cuboid[])
end

function Base.in(a::UnitRange, b::UnitRange)
    amin, amax = extrema(a)
    bmin, bmax = extrema(b)
    bmin <= amin <= bmax && bmin <= amax <= bmax
end

function Base.in(a::Cuboid, b::Cuboid)
    a.x in b.x && a.y in b.y && a.z in b.z
end

function Base.isdisjoint(a::Cuboid, b::Cuboid)
    isdisjoint(a.x, b.x) || isdisjoint(a.y, b.y) || isdisjoint(a.z, b.z)
end

function Base.split(a::UnitRange, b::UnitRange)
    (a in b || isdisjoint(a, b)) && return [a]
    amin, amax = extrema(a)
    bmin, bmax = extrema(b)
    amin < bmin && bmax < amax && return [amin:bmin-1, b, bmax+1:amax]
    amin < bmin && return [amin:bmin-1, bmin:amax]
    [amin:bmax, bmax+1:amax]
end

function Base.split(a::Cuboid, b::Cuboid)
    (a in b || isdisjoint(a, b)) && return [a]
    xs = split(a.x, b.x)
    ys = split(a.y, b.y)
    zs = split(a.z, b.z)
    [Cuboid(x, y, z) for (x, y, z) in Iterators.product(xs, ys, zs)]
end

function volume(a::Cuboid)
    length(a.x) * length(a.y) * length(a.z)
end

function reboot!(reactor::Reactor, step::RebootStep)
    step.on ? turnon!(reactor, step.cuboid) : turnoff!(reactor, step.cuboid)
end

function turnon!(reactor::Reactor, on_cuboid::Cuboid)
    reactor.on_cuboids = turnon(reactor.on_cuboids, on_cuboid)
end

function turnoff!(reactor::Reactor, off_cuboid::Cuboid)
    reactor.on_cuboids = turnoff(reactor.on_cuboids, off_cuboid)
end

function turnon(on_cuboids::Vector{Cuboid}, on_cuboid::Cuboid)
    push!(turnoff(on_cuboids, on_cuboid), on_cuboid)
end

function turnoff(on_cuboids::Vector{Cuboid}, off_cuboid::Cuboid)
    new_on_cuboids = empty(on_cuboids)
    for on_cuboid in on_cuboids
        if isdisjoint(on_cuboid, off_cuboid)
            push!(new_on_cuboids, on_cuboid)
            continue
        else
            split_cuboids = vec(split(on_cuboid, off_cuboid))
            for split_cuboid in split_cuboids
                if isdisjoint(split_cuboid, off_cuboid)
                    push!(new_on_cuboids, split_cuboid)
                end
            end
        end
    end
    new_on_cuboids
end

function Base.parse(::Type{RebootStep}, str::String)
    strs = split(replace(str, " x=" => ",", "y=" => "", "z=" => "", ".." => ","), ",")
    on = popfirst!(strs) == "on"
    xmin, xmax, ymin, ymax, zmin, zmax = parse.(Int, strs)
    RebootStep(on, Cuboid(xmin:xmax, ymin:ymax, zmin:zmax))
end

function isinit(step::RebootStep)
    step.cuboid in Cuboid(-50:50, -50:50, -50:50)
end

function part1(input)
    steps = parse.(RebootStep, eachline(input))
    init_steps = filter(isinit, steps)
    reactor = Reactor()
    for init_step in init_steps
        reboot!(reactor, init_step)
    end
    sum(volume.(reactor.on_cuboids))
end

function part2(input)
    steps = parse.(RebootStep, eachline(input))
    reactor = Reactor()
    for step in steps
        reboot!(reactor, step)
    end
    sum(volume.(reactor.on_cuboids))
end

@show part1("day22/steps.txt")
@show part2("day22/steps.txt")


text:
# Routines for overlaying a hierarchical rectangular grid onto a gmsh CAD model

function gmsh_overlay_rectangular_grid(bb::NTuple{4, Float64},
                                       material::String,
                                       x::Vector{Vector{Float64}},
                                       y::Vector{Vector{Float64}})

    @info "Overlaying rectangular grid"
    # Ensure that the material is of the form "MATERIAL_xxxx"
    if (length(material) < 9)
        @error "material argument must of of the form 'MATERIAL_xxxxxx'"
    end
    if material[1:9] != "MATERIAL_"
        @error "material argument must of of the form 'MATERIAL_xxxxxx'"
    end

    model_dim_tags = gmsh.model.get_entities(2)
    grid_tags = gmsh_generate_rectangular_grid(bb, x, y; material = material)
    grid_dim_tags = [(2, tag) for tag in grid_tags]
    union_of_dim_tags = vcat(model_dim_tags, grid_dim_tags)
    groups = gmsh.model.get_physical_groups()
    names = [gmsh.model.get_physical_name(grp[1], grp[2]) for grp in groups]
    material_indices = findall(x->occursin("MATERIAL", x), names)
    # material hierarchy with the grid material at the bottom.
    material_hierarchy = names[material_indices]
    push!(material_hierarchy, material)
    out_dim_tags = gmsh_group_preserving_fragment(
                        union_of_dim_tags,
                        union_of_dim_tags;
                        material_hierarchy = material_hierarchy
                   )

    return out_dim_tags
end

function gmsh_overlay_rectangular_grid(bb::NTuple{4, Float64},
                                       material::String,
                                       nx::Vector{Int64},
                                       ny::Vector{Int64})

    @info "Overlaying rectangular grid"
    # Ensure that the material is of the form "MATERIAL_xxxx"
    if (length(material) < 9)
        @error "material argument must of of the form 'MATERIAL_xxxxxx'"
    end
    if material[1:9] != "MATERIAL_"
        @error "material argument must of of the form 'MATERIAL_xxxxxx'"
    end

    model_dim_tags = gmsh.model.get_entities(2)
    grid_tags = gmsh_generate_rectangular_grid(bb, nx, ny; material = material)
    grid_dim_tags = [(2, tag) for tag in grid_tags]
    union_of_dim_tags = vcat(model_dim_tags, grid_dim_tags)
    groups = gmsh.model.get_physical_groups()
    names = [gmsh.model.get_physical_name(grp[1], grp[2]) for grp in groups]
    material_indices = findall(x->occursin("MATERIAL", x), names)
    # material hierarchy with the grid material at the bottom.
    material_hierarchy = names[material_indices]
    push!(material_hierarchy, material)
    out_dim_tags = gmsh_group_preserving_fragment(
                        union_of_dim_tags,
                        union_of_dim_tags;
                        material_hierarchy = material_hierarchy
                   )

    return out_dim_tags
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
"""
    WavePacketState

Abstract base class for all states.
"""
abstract type WavePacketState end

"""
    Coherent(α[, N_cutoff])

Describes a coherent state ``|α⟩``. The autmatic choice for `N_cutoff`
yields a normalisation of above 0.999 for α ∈ [0,10].
"""
struct Coherent <: WavePacketState
	α::Complex{Float64}
	N_cutoff::Int
end

coherent_cutoff(α)::Int = ceil(Int, abs(α) * (abs(α) + 3.5) + 0.45 * √(abs(α)))
Coherent(α)                    = Coherent(Complex{Float64}(α), coherent_cutoff(α))
Coherent(α::Float64, N_cutoff) = Coherent(Complex{Float64}(α), N_cutoff)

abstract type Displaced    <: WavePacketState end
abstract type NonDisplaced <: WavePacketState end

#-------------------------------------------------------
# NonDisplaced States

"""
    ArbitraryState(amplitudes[, N_cutoff])

Creates a state with specific amplitudes.
"""
struct ArbitraryState <: NonDisplaced
	amplitudes::Array{Complex{Float64}, 1}
	N_cutoff::Int
end
ArbitraryState(amps::Array{Complex{Float64}, 1}) = ArbitraryState(amps, length(amps) - 1)

"""
    Fock(n[, N_cutoff])

Decribes the Fock state ``|n⟩``.
"""
struct Fock <: NonDisplaced
	n::Int
	N_cutoff::Int
end
Fock(n) = Fock(n, n)

"""
    SqueezedVacuum(r, ϕ = 0.0[, N_cutoff])
		SqueezedVacuum(ξ::Complex{Float64}[, N_cutoff])

Decribes the a squeezed state with squeezing amplitude `r` and squeezing angle `ϕ`, i.e., the
state ``\\exp[(ξ^{*} a^2 - ξ a^{\\dagger 2})/2] |0⟩`` with ``ξ = r e^{i ϕ}``.
"""
struct SqueezedVacuum <: NonDisplaced
	r::Float64
	ϕ::Float64
	N_cutoff::Int
end
# magic number choosen so that for r ∈ [0,3] the normalisation of the state is above 0.99
SqueezedVacuum(r, ϕ = 0.0) = SqueezedVacuum(r, ϕ, 2 * ceil(Int, -1.7 / log(tanh(r))))
SqueezedVacuum(ξ::Complex{Float64}) = SqueezedVacuum(abs(ξ), angle(ξ))
SqueezedVacuum(ξ::Complex{Float64}, N_cutoff) = SqueezedVacuum(abs(ξ), angle(ξ), N_cutoff)

#-------------------------------------------------------
# Displaced States

"""
    DisplacedArbitraryState(α, amplitudes[, N_cutoff])

Describes a state with specific amplitudes, which then is displaced by ``D(α)``.
"""
struct DisplacedArbitraryState <: Displaced
	α::Complex{Float64}
	amplitudes::Array{Complex{Float64}, 1}
	N_cutoff::Int
end
DisplacedArbitraryState(α, amps::Array{Complex{Float64}, 1}) =
	DisplacedArbitraryState(Complex{Float64}(α), amps, length(amps) - 1)

"""
	  DisplacedFock(α, n[, N_cutoff])

Decribes the displaced Fock state ``D(α)|n⟩``.
"""
struct DisplacedFock <: Displaced
	α::Complex{Float64}
	n::Int
	N_cutoff::Int
end
DisplacedFock(α, n) = DisplacedFock(Complex{Float64}(α), n, n)

#-------------------------------------------------------
# Implementation of `createState`

createState(wpt::Union{ArbitraryState, DisplacedArbitraryState}) = wpt.amplitudes

function createState(wpt::Union{Fock, DisplacedFock})
	data = zeros(wpt.N_cutoff + 1)
	wpt.n ≤ wpt.N_cutoff && (data[wpt.n+1] = 1.0)
	return data
end

function createState(wpt::SqueezedVacuum)
	data = zeros(Complex{Float64}, wpt.N_cutoff + 1)

	r, ϕ = wpt.r, wpt.ϕ
	tanh_r, exp_i_ϕ = tanh(r), exp(1im * ϕ)
	for n ∈ 0:wpt.N_cutoff÷2
		data[2n+1] = ( -exp_i_ϕ * tanh_r )^n *
			( √(factorial(2.0*n)) / (2^n * factorial(1.0*n)) )
	end
	return data ./ √(sum(abs2, data))
end


text:
VERSION < v"0.7.0-beta2.199" && __precompile__()

module Combinatorics

using Polynomials

include("numbers.jl")
include("factorials.jl")
include("combinations.jl")
include("permutations.jl")
include("partitions.jl")
include("multinomials.jl")
include("youngdiagrams.jl")

end #module


text:
# Normalize all columns from data frame or array (not in-place)
function zscore(x)

    if typeof(x) == DataFrame
        for c = 1:length(x)
            x[:,c] = zscore(x[:,c])
        end
    else
        x = (x.-mean(skipmissing(x))) / std(skipmissing(x))
        x = convert(Array{Union{Float64, Missings.Missing},1},x)        # Because for some reason, the type of x becomes Any sometimes
    end
    x
end

max_stars_repo_name: caurnhammer/gated-rnns

max_issues_repo_name: caurnhammer/gated-rnns

max_forks_repo_name: caurnhammer/gated-rnns


text:
struct Bounds{T}
    min_val::T
    max_val::T
end

function population_bounds(flip::Flip{Int64}, plan::Plan{Int64}, bounds::Bounds{Int64})::Bool
    @inbounds return (plan.district_populations[flip.old_assignment] - flip.population >= bounds.min_val &&
                      plan.district_populations[flip.new_assignment] + flip.population <= bounds.max_val)
end

function contiguous(flip::Flip{Int64},
                    plan::Plan{Int64},
                    graph::GraphIndex{Int64},
                    neighbors::BitSet)::Bool
    source_node = iterate(neighbors)[1]
    pop!(neighbors, source_node)

    queue = Array{Int64}(undef, graph.n_nodes)
    @inbounds for target_node in neighbors
        visited = zeros(Bool, graph.n_nodes)
        queue[1] = target_node
        visited[target_node] = true
        queue_pos = 1
        queue_len = 2
        found = false
        while queue_len > queue_pos
            curr_node = queue[queue_pos]
            if curr_node == source_node
                found = true
                break
            end
            queue_pos += 1
            for index in 1:graph.neighbors_per_node[curr_node]
                neighbor = graph.node_neighbors[index, curr_node]
                if (!visited[neighbor] &&
                    plan.assignment[neighbor] == flip.old_assignment &&
                    neighbor != flip.node)
                    visited[neighbor] = true
                    queue[queue_len] = neighbor
                    queue_len += 1
                end
            end
        end
        if (queue_len == queue_pos && !found)
            return false
        end
    end
    return true
end

max_stars_repo_name: pjrule/gridwalk

max_issues_repo_name: pjrule/gridwalk

max_forks_repo_name: pjrule/gridwalk


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

include("bat_sample.jl")
include("mcmc/mcmc.jl")
include("sampled_density.jl")
include("importance/importance_sampler.jl")
include("nested_sampling/nested_sampling.jl")
include("partitioned_sampling/partitioned_sampling.jl")


text:
"""Composes and invokes one or more filters.
"""
type CompositeFilter <: Filter
    filters::Vector{Filter}
    state::LifeCycle.State
    CompositeFilter(filters::Vector{Filter}) = new(filters, LifeCycle.INITIALIZED)
end

"Make the filter available for use."
function start(flt::CompositeFilter)
    state!(flt, LifeCycle.STARTING)
    for f in flt.filters
        start(f)
    end
    state!(flt, LifeCycle.STARTED)
end

"Disable the filter."
function stop(flt::CompositeFilter)
    state!(flt, LifeCycle.STOPPING)
    for f in flt.filters
        stop(f)
    end
    state!(flt, LifeCycle.STOPPED)
end

function filter{E <: Event}(flt::CompositeFilter, evnt::E)
    result = FilterResult.NEUTRAL
    for f in flt.filters
        result = filter(f, evnt)
        (result == FilterResult.ACCEPT || result == FilterResult.DENY) && return result
    end
    return result
end

function filter(flt::CompositeFilter, level::Level.EventLevel, marker::MARKER, msg)
    result = FilterResult.NEUTRAL
    for f in flt.filters
        result = filter(f, level, marker, msg)
        (result == FilterResult.ACCEPT || result == FilterResult.DENY) && return result
    end
    return result
end


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

# hard code this for efficiency:
const one_over_half_pi_interval = inv(0.5 * pi_interval(Float64))

"""
Multiply an interval by a positive constant.
For efficiency, does not check that the constant is positive.
"""
multiply_by_positive_constant(α, x::Interval) = @round(α*x.lo, α*x.hi)

half_pi(::Type{Float64}) = multiply_by_positive_constant(0.5, pi_interval(Float64))
half_pi(::Type{T}) where {T} = 0.5 * pi_interval(T)
half_pi(x::T) where {T<:AbstractFloat} = half_pi(T)

two_pi(::Type{Float64}) = multiply_by_positive_constant(2.0, pi_interval(Float64))
two_pi(::Type{T}) where {T} = 2 * pi_interval(T)

range_atan(::Type{T}) where {T<:Real} = Interval(-(pi_interval(T).hi), pi_interval(T).hi)
half_range_atan(::Type{T}) where {T} = (temp = half_pi(T); Interval(-(temp.hi), temp.hi) )
pos_range_atan(::Type{T}) where {T<:Real} = Interval(zero(T), pi_interval(T).hi)

const halfpi = pi / 2.0

"""Finds the quadrant(s) corresponding to a given floating-point
number. The quadrants are labelled as 0 for x ∈ [0, π/2], etc.
For numbers very near a boundary of the quadrant, a tuple of two quadrants
is returned. The minimum or maximum must then be chosen appropriately.

This is a rather indirect way to determine if π/2 and 3π/2 are contained
in the interval; cf. the formula for sine of an interval in
Tucker, *Validated Numerics*.
"""
function find_quadrants(x::T) where {T}
    temp = atomic(Interval{T}, x) / half_pi(x)

    return SVector(floor(temp.lo), floor(temp.hi))
end

function sin(a::Interval{T}) where T
    isempty(a) && return a

    whole_range = Interval{T}(-1, 1)

    diam(a) > two_pi(T).lo && return whole_range

    # The following is equiavlent to doing temp = a / half_pi  and
    # taking floor(a.lo), floor(a.hi)
    lo_quadrant = minimum(find_quadrants(a.lo))
    hi_quadrant = maximum(find_quadrants(a.hi))

    if hi_quadrant - lo_quadrant > 4  # close to limits
        return whole_range
    end

    lo_quadrant = mod(lo_quadrant, 4)
    hi_quadrant = mod(hi_quadrant, 4)

    # Different cases depending on the two quadrants:
    if lo_quadrant == hi_quadrant
        a.hi - a.lo > pi_interval(T).lo && return whole_range  # in same quadrant but separated by almost 2pi
        lo = @round(sin(a.lo), sin(a.lo)) # Interval(sin(a.lo, RoundDown), sin(a.lo, RoundUp))
        hi = @round(sin(a.hi), sin(a.hi)) # Interval(sin(a.hi, RoundDown), sin(a.hi, RoundUp))
        return hull(lo, hi)

    elseif lo_quadrant==3 && hi_quadrant==0
        return @round(sin(a.lo), sin(a.hi)) # Interval(sin(a.lo, RoundDown), sin(a.hi, RoundUp))

    elseif lo_quadrant==1 && hi_quadrant==2
        return @round(sin(a.hi), sin(a.lo)) # Interval(sin(a.hi, RoundDown), sin(a.lo, RoundUp))

    elseif ( lo_quadrant == 0 || lo_quadrant==3 ) && ( hi_quadrant==1 || hi_quadrant==2 )
        return @round(min(sin(a.lo), sin(a.hi)), 1)
        # Interval(min(sin(a.lo, RoundDown), sin(a.hi, RoundDown)), one(T))

    elseif ( lo_quadrant == 1 || lo_quadrant==2 ) && ( hi_quadrant==3 || hi_quadrant==0 )
        return @round(-1, max(sin(a.lo), sin(a.hi)))
        # Interval(-one(T), max(sin(a.lo, RoundUp), sin(a.hi, RoundUp)))

    else #if( lo_quadrant == 0 && hi_quadrant==3 ) || ( lo_quadrant == 2 && hi_quadrant==1 )
        return whole_range
    end
end

function quadrant(x::Float64)
    x_mod2pi = rem2pi(x, RoundNearest)  # gives result in [-pi, pi]

    x_mod2pi < -halfpi && return (2, x_mod2pi)
    x_mod2pi < 0 && return (3, x_mod2pi)
    x_mod2pi <= halfpi && return (0, x_mod2pi)

    return 1, x_mod2pi
end

function sin(a::Interval{Float64})

    T = Float64

    isempty(a) && return a

    whole_range = Interval{T}(-1, 1)

    diam(a) > two_pi(T).lo && return whole_range

    lo_quadrant, lo = quadrant(a.lo)
    hi_quadrant, hi = quadrant(a.hi)

    lo, hi = a.lo, a.hi  # should be able to use the modulo version of a, but doesn't seem to work

    # Different cases depending on the two quadrants:
    if lo_quadrant == hi_quadrant
        a.hi - a.lo > pi_interval(T).lo && return whole_range  #

        if lo_quadrant == 1 || lo_quadrant == 2
            # negative slope
            return @round(sin(hi), sin(lo))
        else
            return @round(sin(lo), sin(hi))
        end

    elseif lo_quadrant==3 && hi_quadrant==0
        return @round(sin(lo), sin(hi))

    elseif lo_quadrant==1 && hi_quadrant==2
        return @round(sin(hi), sin(lo))

    elseif ( lo_quadrant == 0 || lo_quadrant==3 ) && ( hi_quadrant==1 || hi_quadrant==2 )
        return @round(min(sin(lo), sin(hi)), 1)

    elseif ( lo_quadrant == 1 || lo_quadrant==2 ) && ( hi_quadrant==3 || hi_quadrant==0 )
        return @round(-1, max(sin(lo), sin(hi)))

    else #if( lo_quadrant == 0 && hi_quadrant==3 ) || ( lo_quadrant == 2 && hi_quadrant==1 )
        return whole_range
    end
end

function sinpi(a::Interval{T}) where T
    isempty(a) && return a
    w = a * pi_interval(T)
    return(sin(w))
end

function cos(a::Interval{T}) where T
    isempty(a) && return a

    whole_range = Interval(-one(T), one(T))

    diam(a) > two_pi(T).lo && return whole_range

    lo_quadrant = minimum(find_quadrants(a.lo))
    hi_quadrant = maximum(find_quadrants(a.hi))

    if hi_quadrant - lo_quadrant > 4  # close to limits
        return Interval(-one(T), one(T))
    end

    lo_quadrant = mod(lo_quadrant, 4)
    hi_quadrant = mod(hi_quadrant, 4)

    # Different cases depending on the two quadrants:
    if lo_quadrant == hi_quadrant # Interval limits in the same quadrant
        a.hi - a.lo > pi_interval(T).lo && return whole_range
        lo = @round(cos(a.lo), cos(a.lo))
        hi = @round(cos(a.hi), cos(a.hi))
        return hull(lo, hi)

    elseif lo_quadrant == 2 && hi_quadrant==3
        return @round(cos(a.lo), cos(a.hi))

    elseif lo_quadrant == 0 && hi_quadrant==1
        return @round(cos(a.hi), cos(a.lo))

    elseif ( lo_quadrant == 2 || lo_quadrant==3 ) && ( hi_quadrant==0 || hi_quadrant==1 )
        return @round(min(cos(a.lo), cos(a.hi)), 1)

    elseif ( lo_quadrant == 0 || lo_quadrant==1 ) && ( hi_quadrant==2 || hi_quadrant==3 )
        return @round(-1, max(cos(a.lo), cos(a.hi)))

    else#if ( lo_quadrant == 3 && hi_quadrant==2 ) || ( lo_quadrant == 1 && hi_quadrant==0 )
        return whole_range
    end
end

function cos(a::Interval{Float64})

    T = Float64

    isempty(a) && return a

    whole_range = Interval{Float64}(-one(T), one(T))

    diam(a) > two_pi(T).lo && return whole_range

    lo_quadrant, lo = quadrant(a.lo)
    hi_quadrant, hi = quadrant(a.hi)

    lo, hi = a.lo, a.hi

    # Different cases depending on the two quadrants:
    if lo_quadrant == hi_quadrant # Interval limits in the same quadrant
        a.hi - a.lo > pi_interval(T).lo && return whole_range

        if lo_quadrant == 2 || lo_quadrant == 3
            # positive slope
            return @round(cos(lo), cos(hi))
        else
            return @round(cos(hi), cos(lo))
        end

    elseif lo_quadrant == 2 && hi_quadrant==3
        return @round(cos(a.lo), cos(a.hi))

    elseif lo_quadrant == 0 && hi_quadrant==1
        return @round(cos(a.hi), cos(a.lo))

    elseif ( lo_quadrant == 2 || lo_quadrant==3 ) && ( hi_quadrant==0 || hi_quadrant==1 )
        return @round(min(cos(a.lo), cos(a.hi)), 1)

    elseif ( lo_quadrant == 0 || lo_quadrant==1 ) && ( hi_quadrant==2 || hi_quadrant==3 )
        return @round(-1, max(cos(a.lo), cos(a.hi)))

    else #if ( lo_quadrant == 3 && hi_quadrant==2 ) || ( lo_quadrant == 1 && hi_quadrant==0 )
        return whole_range
    end
end

function cospi(a::Interval{T}) where T
    isempty(a) && return a
    w = a * pi_interval(T)
    return(cos(w))
end

function find_quadrants_tan(x::T) where {T}
    temp = atomic(Interval{T}, x) / half_pi(x)

    return SVector(floor(temp.lo), floor(temp.hi))
end

function tan(a::Interval{T}) where T
    isempty(a) && return a

    diam(a) > pi_interval(T).lo && return entireinterval(a)

    lo_quadrant = minimum(find_quadrants_tan(a.lo))
    hi_quadrant = maximum(find_quadrants_tan(a.hi))

    lo_quadrant_mod = mod(lo_quadrant, 2)
    hi_quadrant_mod = mod(hi_quadrant, 2)

    if lo_quadrant_mod == 0 && hi_quadrant_mod == 1
        # check if really contains singularity:
        if hi_quadrant * half_pi(T) ⊆ a
            return entireinterval(a)  # crosses singularity
        end

    elseif lo_quadrant_mod == hi_quadrant_mod && hi_quadrant > lo_quadrant
        # must cross singularity
        return entireinterval(a)

    end

    # @show a.lo, a.hi
    # @show tan(a.lo), tan(a.hi)

    return @round(tan(a.lo), tan(a.hi))
end

function tan(a::Interval{Float64})

    T = Float64

    isempty(a) && return a

    diam(a) > pi_interval(T).lo && return entireinterval(a)

    lo_quadrant, lo = quadrant(a.lo)
    hi_quadrant, hi = quadrant(a.hi)

    lo_quadrant_mod = mod(lo_quadrant, 2)
    hi_quadrant_mod = mod(hi_quadrant, 2)

    if lo_quadrant_mod == 0 && hi_quadrant_mod == 1
        return entireinterval(a)  # crosses singularity

    elseif lo_quadrant_mod == hi_quadrant_mod && hi_quadrant != lo_quadrant
        # must cross singularity
        return entireinterval(a)

    end

    return @round(tan(a.lo), tan(a.hi))
end

function asin(a::Interval{T}) where T

    domain = Interval(-one(T), one(T))
    a = a ∩ domain

    isempty(a) && return a

    return @round(asin(a.lo), asin(a.hi))
end

function acos(a::Interval{T}) where T

    domain = Interval(-one(T), one(T))
    a = a ∩ domain

    isempty(a) && return a

    return @round(acos(a.hi), acos(a.lo))
end

function atan(a::Interval{T}) where T
    isempty(a) && return a

    return @round(atan(a.lo), atan(a.hi))
end

#atan{T<:Real, S<:Real}(y::Interval{T}, x::Interval{S}) = atan(promote(y, x)...)

function atan(y::Interval{Float64}, x::Interval{Float64})
    (isempty(y) || isempty(x)) && return emptyinterval(Float64)

    atomic(Interval{Float64}, atan(big53(y), big53(x)))
end

function atan(y::Interval{BigFloat}, x::Interval{BigFloat})
    (isempty(y) || isempty(x)) && return emptyinterval(BigFloat)

    T = BigFloat

    # Prevent nonsense results when y has a signed zero:
    if y.lo == zero(T)
        y = Interval(zero(T), y.hi)
    end
    if y.hi == zero(T)
        y = Interval(y.lo, zero(T))
    end

    # Check cases based on x
    if x == zero(x)

        y == zero(y) && return emptyinterval(T)
        y.lo ≥ zero(T) && return half_pi(T)
        y.hi ≤ zero(T) && return -half_pi(T)
        return half_range_atan(T)

    elseif x.lo > zero(T)

        y == zero(y) && return y
        y.lo ≥ zero(T) &&
            return @round(atan(y.lo, x.hi), atan(y.hi, x.lo)) # refinement lo bound
        y.hi ≤ zero(T) &&
            return @round(atan(y.lo, x.lo), atan(y.hi, x.hi))
        return @round(atan(y.lo, x.lo), atan(y.hi, x.lo))

    elseif x.hi < zero(T)

        y == zero(y) && return pi_interval(T)
        y.lo ≥ zero(T) &&
            return @round(atan(y.hi, x.hi), atan(y.lo, x.lo))
        y.hi < zero(T) &&
            return @round(atan(y.hi, x.lo), atan(y.lo, x.hi))
        return range_atan(T)

    else # zero(T) ∈ x

        if x.lo == zero(T)
            y == zero(y) && return y

            y.lo ≥ zero(T) && return @round(atan(y.lo, x.hi), half_range_atan(BigFloat).hi)

            y.hi ≤ zero(T) && return @round(half_range_atan(BigFloat).lo, atan(y.hi, x.hi))
            return half_range_atan(T)

        elseif x.hi == zero(T)
            y == zero(y) && return pi_interval(T)
            y.lo ≥ zero(T) && return @round(half_pi(BigFloat).lo, atan(y.lo, x.lo))
            y.hi < zero(T) && return @round(atan(y.hi, x.lo), -(half_pi(BigFloat).lo))
            return range_atan(T)
        else
            y.lo ≥ zero(T) &&
                return @round(atan(y.lo, x.hi), atan(y.lo, x.lo))
            y.hi < zero(T) &&
                return @round(atan(y.hi, x.lo), atan(y.hi, x.hi))
            return range_atan(T)
        end

    end
end


text:
module Complementarity

using JuMP, ConditionalJuMP, Cbc
using Test
using Pkg

# A simple complementarity-based time-stepping rigid body simulation. All
# notation is taken from Stewart & Trinkle "An Implicit Time-Stepping Scheme for
# Rigid Body Dynamics with Coulomb Friction". This particular example solves
# for all N timesteps simultaneously. That's not actually necessary, but it makes
# the code a bit simpler to read.
#
# The model consists of a point mass (visualized as a brick) moving in two dimensions
# with gravity and a single planar surface at y = 0.

const h = 0.05
const μ = 0.5
const n = [0, 1]
const mass = 1.0
const g = [0, -9.81]
const D = [
    1 -1
    1  1
]
const k = size(D, 2)
const e = ones(k)

struct LCPUpdate{T}
    q::Vector{T}
    v::Vector{T}
    β::Vector{T}
    λ::T
    c_n::T
end

function JuMP.getvalue(up::LCPUpdate)
    LCPUpdate(getvalue.((up.q, up.v, up.β, up.λ, up.c_n))...)
end

function update(q, v, model)
    qnext = @variable(model, [1:length(q)], lowerbound=-10, basename="qnext", upperbound=10)
    vnext = @variable(model, [1:length(v)], lowerbound=-10, basename="vnext", upperbound=10)
    β = @variable(model,     [1:k],         lowerbound=0,   basename="β",     upperbound=100)
    λ = @variable(model,                    lowerbound=0,   basename="λ",     upperbound=100)
    c_n = @variable(model,                  lowerbound=0,   basename="c_n",   upperbound=100)

    @constraints model begin
        mass * (vnext - v) .== h * n * c_n .+ h * D * β .+ h * mass * g # (5)
        qnext - q .== h .* vnext # (6)
        n' * qnext >= 0 # (7)
        λ * e + D' * vnext .>= 0 # (8)
        μ * c_n - e' * β >= 0 # (9)
    end

    @disjunction(model, (n' * qnext == 0), (c_n == 0)) # (10)
    Dtv = D' * vnext
    for j in 1:k
        @disjunction(model, ((λ + Dtv[j]) == 0), (β[j] == 0)) # (11)
    end
    @disjunction(model, (μ * c_n - e' * β == 0), (λ == 0)) # (12)

    LCPUpdate(qnext, vnext, β, λ, c_n)
end

"""
Simulate the system using the LCP update. This creates and solves a small
model for every time step 1:N.
"""
function simulate(q0, v0, N)
    q, v = q0, v0
    results = LCPUpdate{Float64}[]
    for i in 1:N
        m = Model(solver=CbcSolver(logLevel=0))
        up = update(q, v, m)
        solve(m)
        push!(results, getvalue(up))
        q = results[end].q
        v = results[end].v
    end
    results
end

"""
Like simulate(), but solves for all N states in a single big
optimization, rather than solving each state separately. This
is more like what we would do if we had a control input we wanted
to optimize over
"""
function optimize(q0, v0, N)::Vector{LCPUpdate{Float64}}
    q, v = q0, v0
    m = Model(solver=CbcSolver(logLevel=0))
    results = []
    for i in 1:N
        up = update(q, v, m)
        push!(results, up)
        q = results[end].q
        v = results[end].v
    end
    solve(m)
    getvalue.(results)
end

q0 = [-1, 0.5]
v0 = [2, 0.5]
N = 30
results1 = simulate(q0, v0, N)
results2 = optimize(q0, v0, N)

@test all([r.q for r in results1] .≈ [r.q for r in results2])
@test all([r.v for r in results1] .≈ [r.v for r in results2])
@test all([r.λ for r in results1] .≈ [r.λ for r in results2])

q = [r.q for r in results1]
v = [r.v for r in results1]
@test q[end] ≈ [-0.16892500000000002, 0]
@test v[end] ≈ [0, 0]

if haskey(Pkg.installed(), "DrakeVisualizer")
    @eval using DrakeVisualizer;
    @eval using CoordinateTransformations
    @eval using GeometryTypes
    DrakeVisualizer.any_open_windows() || DrakeVisualizer.new_window()

    vis = Visualizer()[:block]
    setgeometry!(vis, HyperRectangle(Vec(-0.1, -0.1, 0), Vec(0.2, 0.2, 0.2)))

    for qi in q
        settransform!(vis, Translation(qi[1], 0, qi[2]))
        sleep(h)
    end
end

end


text:
using GLVisualize, GeometryTypes

if !isdefined(:runtests)
	window = glscreen()
    timesignal = loop(linspace(0f0,1f0,360))
end

# use the performance tips to speed this up
# (http://docs.julialang.org/en/release-0.4/manual/performance-tips/)
# the array is 512x512 after all
const N = 512
const range = linspace(-5f0, 5f0, N)

function contour_inner(i, x, y)
    Intensity{1,Float32}(sin(1.3*x*i)*cos(0.9*y)+cos(.8*x)*sin(1.9*y)+cos(y*.2*x))
end
const data = zeros(Intensity{1,Float32}, N, N)

function contourdata(t)
    for i=1:size(data, 1)
        @simd for j=1:size(data, 2)
            @inbounds data[i,j] = contour_inner(t, range[i], range[j])
        end
    end
    data
end

renderable = visualize(map(contourdata, timesignal))

_view(renderable, window, camera=:orthographic_pixel)

if !isdefined(:runtests)
	renderloop(window)
end


text:
module RTLSDRBindings

# --- Print radio config
include("../../Printing.jl");
using .Printing

# -- Loading driver bindings
include("./LibRtlsdr.jl")
using  .LibRtlsdr

using Printf

# Methods extension
import Base:close;

# Symbols exportation
export openRTLSDR;
export updateCarrierFreq!;
export updateSamplingRate!;
export updateGain!;
export recv;
export recv!;
export send;
export print;
export getError
export RTLSDRBinding;

# Utils to generate String containers
"""
Init an empty string of size n, filled with space. Usefull to have container to get string from UHD.
"""
initEmptyString(n) = String(ones(UInt8,n)*UInt8(32))
"""
Retrict the big string used as log container to its usefull part
"""
function truncate(s::String)
    indexes =  findfirst("\0",s)
    if isnothing(indexes)
        # No end of line, controls that it is not only whitespace
        if length(unique(s)) == 1
            # String is empty, abord
            sO = " "
        else
            # Return the string as it is
            s0 = s
        end
    else
        # Truncated output
        sO = s[1: indexes[1] - 1]
    end
    return s
end

mutable struct RTLSDRRxWrapper
    buffer::Vector{UInt8}
end
# --- Main Rx structure
mutable struct RTLSDRRx
    rtlsdr::RTLSDRRxWrapper
    carrierFreq::Float64
    samplingRate::Float64;
    gain::Union{Int,Float64};
    antenna::String;
    packetSize::Csize_t;
    released::Int;
end

# --- Main Tx structure
mutable struct RTLSDRTx
    carrierFreq::Float64
    samplingRate::Float64
    gain::Union{Int,Float64}
    antenna::String
    packetSize::Csize_t
    released::Int
end

# --- Complete structure
mutable struct RTLSDRBinding
    radio::Ref{Ptr{rtlsdr_dev_t}}
    rx::RTLSDRRx
    tx::RTLSDRTx
end

function openRTLSDR(carrierFreq,samplingRate,gain;agc_mode=0,tuner_gain_mode=0)
    # --- Instantiate a new RTLSDR device
    ptr_rtlsdr = Ref{Ptr{rtlsdr_dev_t}}()
    chan = 0
    rtlsdr_open(ptr_rtlsdr,chan)
    # --- Get the instantiate object
    rtlsdr  = ptr_rtlsdr[]
    rtlsdr_reset_buffer(rtlsdr)
    # --- Configure it based on what we want
    # Sampling rate (2MHz max ?)
    rtlsdr_set_sample_rate(rtlsdr,samplingRate)
    samplingRate = rtlsdr_get_sample_rate(rtlsdr)
    # Carrier freq
    rtlsdr_set_center_freq(rtlsdr,carrierFreq)
    carrierFreq = rtlsdr_get_center_freq(rtlsdr)
    # Gain mode
    rtlsdr_set_tuner_gain_mode(rtlsdr,tuner_gain_mode)
    rtlsdr_set_agc_mode(rtlsdr,agc_mode)
    # ----------------------------------------------------
    # --- Wrap all into a custom structure
    # ----------------------------------------------------
    # Instantiate a buffer to handle async receive. Size is arbritrary and will be modified afterwards
    buffer = zeros(UInt8,512)
    rtlsdrRx = RTLSDRRxWrapper(buffer)
    rx = RTLSDRRx(
                  rtlsdrRx,
                  carrierFreq,
                  samplingRate,
                  0,
                  "RX",
                  0,
                  0
                 );
    tx = RTLSDRTx(
                  carrierFreq,
                  samplingRate,
                  0,
                  "TX",
                  0,
                  0
                 );
    radio = RTLSDRBinding(
                          ptr_rtlsdr,
                          rx,
                          tx
    )
    return radio
end

function updateCarrierFreq!(radio::RTLSDRBinding,carrierFreq)
    # --- Set the carrier frequency
    rtlsdr_set_center_freq(radio.radio[],carrierFreq)
    # --- Get the carrier frequency
    carrierFreq = rtlsdr_get_center_freq(radio.radio[])
    # --- Update structure
    radio.rx.carrierFreq = carrierFreq;
    radio.tx.carrierFreq = carrierFreq;
    return carrierFreq;
end
function updateSamplingRate!(radio,samplingRate)
    # --- Set the sampling rate
    rtlsdr_set_sample_rate(radio.radio[],samplingRate)
    # --- Get the sampling rate
    samplingRate = rtlsdr_get_sample_rate(radio.radio[])
    # --- Update Fields
    radio.rx.samplingRate = samplingRate;
    radio.tx.samplingRate = samplingRate;
    return samplingRate
end
function updateGain!(radio,gain)
    # @warn "Analog gain update is not supported for RTLSDR";
end

function Base.print(rx::RTLSDRRx);
    strF  = @sprintf(" Carrier Frequency: %2.3f MHz\n Sampling Frequency: %2.3f MHz\n",rx.carrierFreq/1e6,rx.samplingRate/1e6)
    @inforx "Current Simulated Radio Configuration in Rx mode\n$strF";
end
function Base.print(tx::RTLSDRTx);
    strF  = @sprintf(" Carrier Frequency: %2.3f MHz\n Sampling Frequency: %2.3f MHz\n",tx.carrierFreq/1e6,tx.samplingRate/1e6)
    @infotx "Current Simulated Radio Configuration in Rx mode\n$strF";
end
function Base.print(radio::RTLSDRBinding)
    print(radio.rx);
    print(radio.tx);
end

function recv!(sig::Vector{ComplexF32},radio::RTLSDRBinding)
    nbSamples = length(sig)
    # --- Instantiation of UInt8 buffer
    if length(radio.rx.rtlsdr.buffer) !== 2*nbSamples
        # We need to redimension the size of the receive buffer
        buffer = zeros(UInt8,nbSamples * 2)
        radio.rx.rtlsdr.buffer = buffer
    end
    pointerSamples   = Ref{Cint}(0)
    # --- Call the receive method
    ptr = pointer(radio.rx.rtlsdr.buffer,1)
    rtlsdr_read_sync(radio.radio[],ptr,nbSamples*2,pointerSamples)
    # --- Get the number of received complex symbols
    nbElem = pointerSamples[] ÷ 2
    # --- From Bytes to complex
    byteToComplex!(sig,radio.rx.rtlsdr.buffer)
    return nbElem
end

function recv(radio::RTLSDRBinding,nbSamples)
    sig = zeros(ComplexF32,nbSamples)
    nbElem = recv!(sig,radio)
    (nbElem == 0) && @warnrx "No received samples !"
    return sig
end

""" Transform the input buffer of UInt8 into a complexF32 vector
"""
function byteToComplex!(sig::Vector{ComplexF32},buff::Vector{UInt8})
    nbS = length(sig)
    @assert (length(buff) == 2nbS) "Complex buffer size should be 2 times lower (here $(length(sig))  than input Buffersize (here $(length(buff))"
    @inbounds @simd  for n ∈ 1 : nbS
        sig[n] = _btoc(buff[2(n-1)+1]) + 1im*_btoc(buff[2(n-1)+2])
    end
end
function byteToComplex(buff::Vector{UInt8})
    N = length(buff)
    nbS = N ÷ 2
    sig = zeros(ComplexF32,nbS)
    byteToComplex!(sig,buff)
end

""" Byte to float conversion
"""
function _btoc(in::UInt8)
    return in / 128 -1
end

function close(radio::RTLSDRBinding)
    rtlsdr_close(radio.radio[])
    radio.rx.released = 1;
    radio.tx.released = 1;
    @info "RTL-SDR device is now close"
end

function send(sig::Vector{Complex{Cfloat}},radio::RTLSDRBinding)
    @warntx "Unsupported send method for RTLSDR";
end

function scan()
    # --- Using counting method to get the nbumber of devioes
    nE = LibRtlsdr.rtlsdr_get_device_count()
    if nE > 0
        # Found a RTL SDR dongle
        println("Found $nE RTL-SDR dongle")
        # manufact = initEmptyString(200)
        # product = initEmptyString(200)
        # serial = initEmptyString(200)
        # radio = openRTLSDR(800e6,1e6,20)
        # rtlsdr_get_usb_strings(radio.radio[],manufact,product,serial)
        # (!isempty(manufact)) && (println("Manufact = $(truncate(manufact))"))
        # (!isempty(product)) && (println("Product  = $(truncate(product))"))
        # (!isempty(serial)) && (println("Serial   = $(truncate(serial))"))
        # close(radio)
    end
    return nE
end

function getError(radio::RTLSDRBinding)
    return 0
end
function getMD(radio::RTLSDRBinding)
    return 0
end
end


text:
using AbstractPattern
using AbstractPattern.BasicPatterns

const backend = MK(RedyFlavoured)
"""
[a, b..., c] -> :vec3 => [a], b, [c]
[a, b, c]    -> :vec => [a, b, c]
"""
function ellipsis_split(args::AbstractArray{T, 1}) where T
    ellipsis_index = findfirst(args) do arg
        Meta.isexpr(arg, :...)
    end
    if isnothing(ellipsis_index)
        Val(:vec) => args
    else
        Val(:vec3) => (
            args[1:ellipsis_index-1],
            args[ellipsis_index].args[1],
            args[ellipsis_index+1:end]
        )
    end
end

function P_partial_struct_decons(t, partial_fields, ps, prepr::AbstractString="$t")
    function tcons(_...)
        t
    end
    comp = PComp(
        prepr, tcons;
    )
    function extract(sub, i::Int)
        :($sub.$(partial_fields[i]))
    end
    decons(comp, extract, ps)
end

basic_ex2tf(eval::Function, a) =
    isprimitivetype(typeof(a)) ? literal(a) : error("invalid literal $a")
basic_ex2tf(eval::Function, l::LineNumberNode) = wildcard
basic_ex2tf(eval::Function, q::QuoteNode) = literal(q)
basic_ex2tf(eval::Function, s::String) = literal(s)
basic_ex2tf(eval::Function, n::Symbol) =
    n === :_ ?  wildcard : P_capture(n)

Base.@pure function qt2ex(ex::Any)
    if ex isa Expr
        Meta.isexpr(ex, :$) && return ex.args[1]
        Expr(:call, Expr, QuoteNode(ex.head), Expr(:vect, (qt2ex(e) for e in ex.args)...))
    elseif ex isa Symbol
        QuoteNode(ex)
    else
        ex
    end
end

function basic_ex2tf(eval::Function, ex::Expr)
    !(x) = basic_ex2tf(eval, x)
    hd = ex.head; args = ex.args; n_args = length(args)
    if hd === :||
        @assert n_args === 2
        l, r = args
        or(!l, !r)
    elseif hd === :&&
        @assert n_args === 2
        l, r = args
        and(!l, !r)

    elseif hd === :if
        @assert n_args === 2
        cond = args[1]
        guard() do _, scope, _
            see_captured_vars(cond, scope)
        end
    elseif hd === :&
        @assert n_args === 1
        val = args[1]
        guard() do target, scope, _
            see_captured_vars(:($target == $val), scope)
        end
    elseif hd === :let
        bind = args[1]
        @assert bind isa Expr
        if bind.head === :(=)
            @assert bind.args[1] isa Symbol
            P_bind(bind.args[1], bind.args[2], see_capture=true)
        else
            @assert bind.head === :block
            binds = Function[P_bind(arg.args[1], arg.args[2], see_capture=true) for arg in bind.args]
            push!(binds, wildcard)
            and(binds)
        end
    elseif hd === :(::)
        if n_args === 2
            p, ty = args
            ty = eval(ty)::TypeObject
            and(P_type_of(ty), !p)
        else
            @assert n_args === 1
            p = args[1]
            ty = eval(ty)::TypeObject
            P_type_of(ty)
        end
    elseif hd === :vect
        tag, split = ellipsis_split(args)
        return tag isa Val{:vec} ?
            P_vector([!e for e in split]) :
            let (init, mid, tail) = split
                P_vector3(
                    [!e for e in init],
                    !mid,
                    [!e for e in tail]
                )
            end
    elseif hd === :tuple
        P_tuple([!e for e in args])
    elseif hd === :call
        let f = args[1],
            args′ = view(args, 2:length(args))
            n_args′ = n_args - 1
            t = eval(f)
            if t === Core.svec
                tag, split = ellipsis_split(args′ )
                return tag isa Val{:vec} ?
                    P_svec([!e for e in split]) :
                    let (init, mid, tail) = split
                        P_svec3(
                            [!e for e in init],
                            !mid,
                            [!e for e in tail]
                        )
                    end
            end
            all_field_ns = fieldnames(t)
            partial_ns = Symbol[]
            patterns = Function[]
            if n_args′ >= 1 && Meta.isexpr(args′[1], :parameters)
                kwargs = args′[1].args
                args′ = view(args′, 2:length(args′))
            else
                kwargs = []
            end
            if length(all_field_ns) === length(args′)
                append!(patterns, [!e for e in args′])
                append!(partial_ns, all_field_ns)
            elseif length(all_field_ns) !== 0
                error("count of positional fields should be 0 or the same as the fields($all_field_ns)")
            end
            for e in kwargs
                if e isa Symbol
                    e in all_field_ns || error("unknown field name $e for $t when field punnning.")
                    push!(partial_ns, e)
                    push!(patterns, P_capture(e))
                elseif Meta.isexpr(e, :kw)
                    key, value = e.args
                    key in all_field_ns || error("unknown field name $key for $t when field punnning.")
                    @assert key isa Symbol
                    push!(partial_ns, key)
                    push!(patterns, and(P_capture(key), !value))
                end
            end
            P_partial_struct_decons(t, partial_ns, patterns)
        end
    elseif hd === :quote
        !qt2ex(args[1])
    else
        error("not implemented expr=>pattern rule for '($hd)' Expr.")
    end
end

const case_sym = Symbol("@case")
"""a minimal implementation of switch
"""
macro switch(val, ex)
    @assert Meta.isexpr(ex, :block)
    clauses = Union{LineNumberNode, Pair{<:Function, Symbol}}[]
    body = Expr(:block)
    alphabeta = 'a':'z'
    base = gensym()
    k = 0
    for i in eachindex(ex.args)
        stmt = ex.args[i]
        if Meta.isexpr(stmt, :macrocall) &&
           stmt.args[1] === case_sym &&
           length(stmt.args) == 3

            pattern = basic_ex2tf(__module__.eval, stmt.args[3])
            br :: Symbol = Symbol(alphabeta[i % 26], i <= 26 ? "" : string(i), base)
            push!(clauses,  pattern => br)
            push!(body.args, :(@label $br))
        else
            if stmt isa LineNumberNode
                push!(clauses, stmt)
            end
            push!(body.args, stmt)
        end
    end

    match_logic = backend(val, clauses, __source__)
    esc(Expr(
        :block,
        match_logic,
        body
    ))
end

val = 1
@switch val begin
    @case 1
        println("case 1!")
        return
    @case _
        println("case 2!")
        return
end

val = [1, 2, 3, 4]
@switch val begin
    @case (1, 2, a)
        println("emmmm")
        return
    @case [1, (b::AbstractArray{Int, 1})..., 4]
        println("b:", b, " sum(b):", sum(b))
        return
end

@switch :(a = (1, 2, 3)) begin
    @case :(a = (1, $(b...)))
        println(b)
        return
end

@switch (1, 5) begin
    @case (1, 2)
        println(b)
        return
    @case (a, &(a + 4)) && let x = Int[a, 1], x = x .+ a end
        println("this run ", x)
        return
end

@switch Tuple{Int, Int}.parameters begin
    @case Core.svec(a, b)
        println(a, b)
        return
end

@switch :(&ooo) begin
    @case :(&$a)
        println(a)
        return
end

# import MLStyle
# import Match
# import Rematch

# function f_rematch(value)
#     Rematch.@match value begin
#         _::String => :string
#         (2, a, 3) => (:some_tuple, a)
#         [1, a..., 3, 4] => (:some_vector, a)
#     end
# end

# function f_match(value)
#     Match.@match value begin
#         _::String => :string
#         (2, a, 3) => (:some_tuple, a)
#         [1, a..., 3, 4] => (:some_vector, a)
#     end
# end

# function f_mlstyle(value)
#     Rematch.@match value begin
#         _::String => :string
#         (2, a, 3) => (:some_tuple, a)
#         [1, a..., 3, 4] => (:some_vector, a)
#     end
# end

# function f_this(value)
#     @switch value begin
#     @case _::String
#         return :string
#     @case (2, a, 3)
#         return (:some_tuple, a)
#     @case [1, a..., 3, 4]
#         return (:some_vector, a)
#     end
# end

# using BenchmarkTools
# data = [
#     "asda",
#     (2, 1 ,3),
#     [1, 2, 2, 3, 4]
# ]

# fs = [f_rematch, f_match, f_mlstyle, f_this]

# for datum in data
#     for f in fs
#         @info :testing f datum
#         println(@btime $f($datum))
#     end
#     println("===================")
# end

max_stars_repo_name: KristofferC/AbstractPattern

max_issues_repo_name: KristofferC/AbstractPattern

max_forks_repo_name: KristofferC/AbstractPattern


text:
function ViennaOneModel(arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jdouble)
    return ViennaOneModel((Vector{jdouble}, Vector{jdouble}, jdouble), arg0, arg1, arg2)
end

function ViennaOneModel(arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jdouble, arg3::TimeScale)
    return ViennaOneModel((Vector{jdouble}, Vector{jdouble}, jdouble, TimeScale), arg0, arg1, arg2, arg3)
end

function compute_zenith_delay(obj::ViennaOneModel, arg0::RealFieldElement, arg1::Vector{RealFieldElement}, arg2::FieldAbsoluteDate)
    return jcall(obj, "computeZenithDelay", Vector{RealFieldElement}, (RealFieldElement, Vector{RealFieldElement}, FieldAbsoluteDate), arg0, arg1, arg2)
end

function compute_zenith_delay(obj::ViennaOneModel, arg0::jdouble, arg1::Vector{jdouble}, arg2::AbsoluteDate)
    return jcall(obj, "computeZenithDelay", Vector{jdouble}, (jdouble, Vector{jdouble}, AbsoluteDate), arg0, arg1, arg2)
end

function get_parameters_drivers(obj::ViennaOneModel)
    return jcall(obj, "getParametersDrivers", List, ())
end

function mapping_factors(obj::ViennaOneModel, arg0::RealFieldElement, arg1::RealFieldElement, arg2::Vector{RealFieldElement}, arg3::FieldAbsoluteDate)
    return jcall(obj, "mappingFactors", Vector{RealFieldElement}, (RealFieldElement, RealFieldElement, Vector{RealFieldElement}, FieldAbsoluteDate), arg0, arg1, arg2, arg3)
end

function mapping_factors(obj::ViennaOneModel, arg0::jdouble, arg1::jdouble, arg2::Vector{jdouble}, arg3::AbsoluteDate)
    return jcall(obj, "mappingFactors", Vector{jdouble}, (jdouble, jdouble, Vector{jdouble}, AbsoluteDate), arg0, arg1, arg2, arg3)
end

function path_delay(obj::ViennaOneModel, arg0::RealFieldElement, arg1::RealFieldElement, arg2::Vector{RealFieldElement}, arg3::FieldAbsoluteDate)
    return jcall(obj, "pathDelay", RealFieldElement, (RealFieldElement, RealFieldElement, Vector{RealFieldElement}, FieldAbsoluteDate), arg0, arg1, arg2, arg3)
end

function path_delay(obj::ViennaOneModel, arg0::jdouble, arg1::jdouble, arg2::Vector{jdouble}, arg3::AbsoluteDate)
    return jcall(obj, "pathDelay", jdouble, (jdouble, jdouble, Vector{jdouble}, AbsoluteDate), arg0, arg1, arg2, arg3)
end


text:
for OP in (:(&), :(|), :(⊻))
    @eval begin
        @inline function $OP(x::T, y::T) where T<:SafeInteger
            ix = baseint(x)
            iy = baseint(y)
            result = $OP(ix, iy)
            return safeint(result)
        end
    end
end

for OP in (:(&), :(|), :(⊻))
    for (T1, T2) in MixedTypes
        @eval begin
            @inline function $OP(x::$T1, y::$T2)
                xx, yy = promote(x, y)
                return $OP(xx, yy)
            end
            @inline function $OP(x::$T2, y::$T1)
                xx, yy = promote(x, y)
                return $OP(xx, yy)
            end
        end
   end
end

for OP in (:(<), :(<=), :(>=), :(>), :(!=), :(==), :isless, :isequal)
    @eval begin
        @inline function $OP(x::T, y::T) where T<:SafeInteger
            ix = baseint(x)
            iy = baseint(y)
            result = $OP(ix, iy)
            return safeint(result)
        end
    end
end

for OP in (:(<), :(<=), :(>=), :(>), :(!=), :(==), :isless, :isequal)
    for (T1, T2) in MixedTypes
        @eval begin
            @inline function $OP(x::$T1, y::$T2)
                xx, yy = promote(x, y)
                return $OP(xx, yy)
            end
            @inline function $OP(x::$T2, y::$T1)
                xx, yy = promote(x, y)
                return $OP(xx, yy)
            end
        end
    end
end

for OP in (:(>>>), :(>>), :(<<))
    @eval begin
       @inline function $OP(x::T, y::T) where T<:SafeInteger
            r1 = baseint(x)
            r2 = baseint(y)
            bitsof(T) < abs(r2) && throw(OverflowError("cannot shift $T by $y"))
            result = $OP(r1, r2)
            return reinterpret(T, result)
        end

        @inline function $OP(x::T1, y::T2) where {T1<:SafeInteger, T2<:SafeInteger}
            xx = baseint(x)
            yy = baseint(y)
            bitsof(T1) < abs(yy) && throw(OverflowError("cannot shift $T1 by $yy"))
            return reinterpret(T1, $OP(xx, yy))
        end

        @inline function $OP(x::T1, y::Int) where {T1<:SafeInteger}
            xx = baseint(x)
            bitsof(T1) < abs(y) && throw(OverflowError("cannot shift $T1 by $y"))
            return reinterpret(T1, $OP(xx, y))
        end

        @inline function $OP(x::T1, y::T2) where {T1<:SafeInteger, T2<:Base.BitSigned}
            xx = baseint(x)
            bitsof(T1) < abs(y) && throw(OverflowError("cannot shift $T1 by $y"))
            return reinterpret(T1, $OP(xx, y))
        end
        @inline function $OP(x::T1, y::T2) where {T1<:Base.BitSigned, T2<:SafeInteger}
            yy = baseint(y)
            bitsof(T1) < abs(yy) && throw(OverflowError("cannot shift $T1 by $yy"))
            return reinterpret(T1, $OP(x, yy))
        end
        @inline function $OP(x::T1, y::T2) where {T1<:SafeInteger, T2<:Base.BitUnsigned}
            xx = baseint(x)
            bitsof(T1) < abs(y) && throw(OverflowError("cannot shift $T1 by $y"))
            return reinterpret(T1, $OP(xx, y))
        end
        @inline function $OP(x::T1, y::T2) where {T1<:Base.BitUnsigned, T2<:SafeInteger}
            yy = baseint(y)
            bitsof(T1) < abs(yy) && throw(OverflowError("cannot shift $T1 by $yy"))
            return reinterpret(T1, $OP(x, yy))
        end
   end
end


text:
"""Tunes `param` until `expr` evaluates to zero within δ.
`expr` must be monotonically increasing in terms of `param`."""
macro binary_opt(expr, param, min, max, δ)
    :(m = $min;M = $max;z=zero($δ);
    while(true)
        $(esc(param)) = (m+M)/2
        println("$((m+M)/2)")
        e = $(esc(expr))
        if norm(e) < $δ
            break
        elseif e < z
            m = $(esc(param))
        else
            M = $(esc(param))
        end
    end)
end

"""Returns the location and values of local maxima in `sig`."""
function local_maxima(sig)
    l = length(sig)
    res = Tuple{Int, T}[]
    if l < 3
        return res
    end
    @inbounds for i = 2:l-1
        if sig[i-1] < sig[i] && sig[i+1] < sig[i]
            append!(res, (i, sig[i]))
        end
    end
    res
end


text:
struct RejectionMethod{X,R,U,DG<:DependencyGraph} <: DecoupledSearchAlgorithm
  state::X
  model::R
  interval_lo::Vector{Int}
  interval_hi::Vector{Int}
  rates::Vector{U}
  total_rate::U
  spc_dg::DG
  rxn_dg::DG
end

# need to think about whether its possible to do this with cumulative rates
RatesCache(::RejectionMethod) = HasRates()

@inline function initialize!(algorithm::RejectionMethod, state, model)
  interval_lo = algorithm.interval_lo
  interval_hi = algorithm.interval_hi
  rates = algorithm.rates
  total_rate = algorithm.total_rate

  fill!(total_rate, zero(eltype(total_rate)))

  # build the upper and lower bounds on state
  for i in eachindex(state)
    interval_lo[i], interval_hi[i] = make_interval(state[i], 0.1)
  end

  # fill in the propensities

  for j in eachindex(rates)
    lo = rate(model, interval_lo, j)
    hi = rate(model, interval_hi, j)

    @inbounds rates[j][1] = lo
    @inbounds rates[j][2] = hi

    @inbounds total_rate[1] += lo
    @inbounds total_rate[2] += hi
  end
  return nothing
end

@inline function generate_jump(algorithm::RejectionMethod)
  j, logu = accept_reject(RatesCache(algorithm), algorithm.rates, algorithm.total_rate, algorithm.state, algorithm.model)
  # j, u = accept_reject(RatesCache(algorithm), algorithm.rates, algorithm.total_rate, algorithm.state, algorithm.model)
  s = logu / upper(algorithm.total_rate)
  # s = - log(u) / upper(algorithm.total_rate)

  return (j, s)
end

# rates should be an iterator over intervals_rates
@inline function accept_reject(ctype::S, rates, total_rate, state, model) where S
  logu = zero(eltype(total_rate))
  # r = one(eltype(total_rate))
  rejected = true
  j = 0

  while rejected
    # search jumps based on upper bounds
    j = search_jump_rates(ctype, iterate_upper(rates), upper(total_rate))

    lo = lower(rates, j)
    hi = upper(rates, j)
    u  = rand() * hi

    rejected = (lo > zero(lo)) && (lo < u)

    if rejected
      rate_j = rate(model, state, j)
      rejected = (rate_j > zero(rate_j)) && (rate_j < u)
    end

    logu += randexp()
  end

  return j, logu
end

@inline function update!(algorithm::RejectionMethod, state, model, k)
  spc_dg = algorithm.spc_dg
  rxn_dg = algorithm.rxn_dg
  interval_lo = algorithm.interval_lo
  interval_hi = algorithm.interval_hi
  rates = algorithm.rates
  total_rate = algorithm.total_rate

  # for i in dependents(rxn_dg, k)
  for i in eachindex(state)
    @inbounds x = state[i]
    @inbounds x_lo = interval_lo[i]
    @inbounds x_hi = interval_hi[i]

    # check for states that have left interval
    if (x == 0) || (x < x_lo) || (x > x_hi)
      @inbounds interval_lo[i], interval_hi[i] = make_interval(x, 0.1)

      # check which propensity bounds must be updated based on a species-reaction dependency graph
      for j in dependents(spc_dg, i)
        @inbounds total_rate[1] -= rates[j][1]
        @inbounds total_rate[2] -= rates[j][2]

        @inbounds rates[j][1] = rate(model, interval_lo, j)
        @inbounds rates[j][2] = rate(model, interval_hi, j)

        @inbounds total_rate[1] += rates[j][1]
        @inbounds total_rate[2] += rates[j][2]
      end
    end
  end
  return nothing
end

function make_interval(x, p)
  if x == zero(x)
    lo = zero(x)
    hi = zero(x)
  elseif x < 25
    lo = max(0, x-4)
    hi = x+4
  else
    lo = trunc(Int, (1-p)*x)
    hi = trunc(Int, (1+p)*x)
  end

  return (lo, hi)
end

lower(x) = first(x)
upper(x) = last(x)

@inbounds lower(x, j) = first(x[j])
@inbounds upper(x, j) = last(x[j])

@inbounds iterate_upper(x) = (upper(x, j) for j in eachindex(x))
@inbounds iterate_lower(x) = (lower(x, j) for j in eachindex(x))


text:
push!(LOAD_PATH,"simulator")
using EVQueues, PyPlot

lambda=100.0;
mu=1.0;
gamma=0.2;

Tfinal=100.0;

#@time sim = ev_parallel(lambda,mu,gamma,Tfinal,C)
@time sim = ev_exact(lambda,mu,gamma,Tfinal)
@time sim2 = ev_exact(lambda,mu,999999.0,Tfinal)

#@time compute_statistics!(sim)


text:
using SerializationCaches
using Documenter

makedocs(modules=[SerializationCaches],
         sitename="SerializationCaches",
         authors="Beacon Biosignals and other contributors",
         pages=["API Documentation" => "index.md"])

deploydocs(repo="github.com/beacon-biosignals/SerializationCaches.jl.git")


text:
using Test
using SafeTestsets

@safetestset "Quadrature rules" begin include("quadrule_test.jl") end

@safetestset "Singular quadrature rules" begin include("singularquadrule_test.jl") end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
module ContractionOrder
    using LightGraphs
    using LightGraphs: SimpleEdge

    function log2sumexp2(s)
        ms = maximum(s)
        return log2(sum(x->exp2(x - ms), s)) + ms
    end

    function indexof(x, v)
        @inbounds for i=1:length(v)
            if v[i] == x
                return i
            end
        end
        return 0
    end

    include("edgeinfo.jl")
    include("greedy.jl")
end

using .ContractionOrder: order_greedy, TensorNetworkLayout
import .ContractionOrder: abstract_contract, log2sumexp2
using LightGraphs: SimpleEdge, src, dst

export trees_greedy, abstract_contract, log2sumexp2

function build_trees(N::Int, order::AbstractVector)
    ids = collect(Any, 1:N)
    for i=1:length(order)
        ta, tb = src(order[i]), dst(order[i])
        ids[ta] = ContractionTree(ids[ta], ids[tb])
        ids[tb] = nothing
    end
    filter(x->x!==(nothing), ids)
end

function build_order(trees::AbstractVector)
    res = SimpleEdge[]
    for t in trees
        build_order!(t, res)
    end
    return res
end

function build_order!(tree, res)
    if tree isa Integer
        return tree
    else
        a = build_order!(tree.left, res)
        b = build_order!(tree.right, res)
        push!(res, SimpleEdge(a, b))
        return min(a, b)
    end
end

function layout(tn::TensorNetwork)
    graph = tn2graph(tn)
    log2shapes = [[log2.(size(tn.tensors[i]))...] for i=1:length(tn)]
    TensorNetworkLayout(graph, log2shapes)
end

"""
    trees_greedy(tn::TensorNetwork, strategy="min_dim")

Returns a tuple of `(time complexities, space complexities, trees)`,
where `trees` is a vector of `ContractionTree` objects.
For disconnected graphs, the number of trees can be greater than 1.

`strategy` can be "min_dim", "min_reduce", "max_reduce", "min_reduce_tri" or "max_reduce_tri".
"""
function trees_greedy(tn::TensorNetwork; kwargs...)
    tc, sc, order = order_greedy(layout(tn); kwargs...)
    tc, sc, build_trees(length(tn), order)
end

abstract_contract(tn::TensorNetwork, trees::ContractionTree) = abstract_contract(tn, [trees])

function abstract_contract(tn::TensorNetwork, trees::AbstractVector)
    abstract_contract(layout(tn), build_order(trees))
end


text:
using Revise, Test, ForwardDiff, Parameters, Setfield, Plots, LinearAlgebra
using BifurcationKit
const BK = BifurcationKit

D(f, x, p, dx) = ForwardDiff.derivative(t->f(x .+ t .* dx, p), 0.)
norminf = x -> norm(x, Inf)
####################################################################################################
function TMvf!(dz, z, p, t)
	@unpack J, α, E0, τ, τD, τF, U0 = p
	E, x, u = z
	SS0 = J * u * x * E + E0
	SS1 = α * log(1 + exp(SS0 / α))
	dz[1] = (-E + SS1) / τ
	dz[2] =	(1.0 - x) / τD - u * x * E
	dz[3] = (U0 - u) / τF +  U0 * (1.0 - u) * E
	dz
end

TMvf(z, p) = TMvf!(similar(z), z, p, 0)
dTMvf = (z,p) -> ForwardDiff.jacobian(x-> TMvf(x,p), z)

# we group the differentials together
jet  = BK.get3Jet(TMvf, dTMvf)

par_tm = (α = 1.5, τ = 0.013, J = 3.07, E0 = -2.0, τD = 0.200, U0 = 0.3, τF = 1.5, τS = 0.007) #2.87
z0 = [0.238616, 0.982747, 0.367876 ]

opts_br = ContinuationPar(pMin = -10.0, pMax = -0.9, ds = 0.04, dsmax = 0.125, nInversion = 8, detectBifurcation = 3, maxBisectionSteps = 25, nev = 3)
	opts_br = @set opts_br.newtonOptions.verbose = false
	br, = continuation(TMvf, dTMvf, z0, par_tm, (@lens _.E0), opts_br;
	printSolution = (x, p) -> (E = x[1], x = x[2], u = x[3]),
	tangentAlgo = BorderedPred(),
	plot = true, verbosity = 0, normC = norminf)

plot(br, plotfold=false, markersize=4, legend=:topleft)
####################################################################################################
hopfpt = computeNormalForm(jet..., br, 4)

# newton parameters
optn_po = NewtonPar(verbose = true, tol = 1e-8,  maxIter = 8)

# continuation parameters
opts_po_cont = ContinuationPar(dsmax = 0.1, ds= -0.0001, dsmin = 1e-4, pMax = 0., pMin=-5., maxSteps = 110, newtonOptions = (@set optn_po.tol = 1e-7), nev = 2, precisionStability = 1e-8, detectBifurcation = 3, plotEveryStep = 10, saveSolEveryStep=1)

Mt = 200 # number of sections
	br_potrap, utrap = continuation(
	jet..., br, 4, opts_po_cont,
	PeriodicOrbitTrapProblem(M = Mt);
	linearPO = :Dense,
	verbosity = 2,	plot = true,
	printSolution = (x, p) -> (xtt=reshape(x[1:end-1],3,Mt); return (max = maximum(xtt[1,:]), min = minimum(xtt[1,:]), period = x[end])),
	plotSolution = (x, p; k...) -> begin
		xtt = BK.getTrajectory(p.prob, x, p.p)
		plot!(xtt.t, xtt.u[1,:]; label = "E", k...)
		plot!(xtt.t, xtt.u[2,:]; label = "x", k...)
		plot!(xtt.t, xtt.u[3,:]; label = "u", k...)
		plot!(br,subplot=1, putbifptlegend = false)
		end,
	normC = norminf)

plot(br, br_potrap, markersize = 3)
	plot!(br_potrap.param, br_potrap.min, label = "")

####################################################################################################
# idem with Standard shooting
using DifferentialEquations

# this is the ODEProblem used with `DiffEqBase.solve`
probsh = ODEProblem(TMvf!, copy(z0), (0., 1000.), par_tm; atol = 1e-10, rtol = 1e-9)

opts_po_cont = ContinuationPar(dsmax = 0.05, ds= -0.0001, dsmin = 1e-4, pMax = 0., pMin=-5., maxSteps = 210, newtonOptions = (@set optn_po.tol = 1e-6), nev = 25, precisionStability = 1e-8, detectBifurcation = 0, plotEveryStep = 10, saveSolEveryStep=0)

br_posh, = @time continuation(
	jet..., br, 4,
	# arguments for continuation
	opts_po_cont,
	# this is where we tell that we want Standard Shooting
	ShootingProblem(15, par_tm, probsh, Rodas4(), parallel = true);
	ampfactor = 1.0, δp = 0.0005,
	usedeflation = true,
	linearPO = :autodiffDense,
	updateSectionEveryStep = 2,
	verbosity = 2,	plot = true,
	printSolution = (x, p) -> (return (max = getMaximum(p.prob, x, @set par_tm.E0 = p.p), period = getPeriod(p.prob, x, @set par_tm.E0 = p.p))),
	plotSolution = (x, p; k...) ->
		begin
			xtt = BK.getTrajectory(p.prob, x, @set par_tm.E0 = p.p)
			plot!(xtt; legend = false, k...);
			plot!(br, subplot=1, putspecialptlegend = false)
		end,
	normC = norminf)

plot(br_posh, br, markersize=3)
plot(br, br_potrap, br_posh, markersize=3)
####################################################################################################
# idem with Poincaré shooting

function TMvfExtended!(dz, z, p, t)
	# we write the first part
	TMvf!(dz, z, p, t)
	dz[4:end] .= @views ForwardDiff.derivative(t -> TMvf(z[1:3] .+ t .* z[4:end], p), 0)
end

probmono = ODEProblem(TMvfExtended!, vcat(u0, u0), (0., 1000.), par_tm; atol = 1e-10, rtol = 1e-9)

opts_po_cont = ContinuationPar(dsmax = 0.02, ds= -0.001, dsmin = 1e-5, pMax = 0., pMin=-5., maxSteps = 200, newtonOptions = NewtonPar(optn_po;tol = 1e-6, maxIter=15), nev = 25, precisionStability = 1e-8, detectBifurcation = 0, plotEveryStep = 5, saveSolEveryStep = 2)

br_popsh, = @time continuation(
	jet..., br, 4,
	# arguments for continuation
	opts_po_cont,
	# this is where we tell that we want Poincaré Shooting
	PoincareShootingProblem(9, par_tm, probsh, Rodas4P(), probmono, Rosenbrock23(), parallel = false);
	ampfactor = 1.0, δp = 0.0005,
	usedeflation = false,
	linearPO = :autodiffDense,
	updateSectionEveryStep = 2,
	verbosity = 2,	plot = true,
	printSolution = (x, p) -> (return (max = getMaximum(p.prob, x, @set par_tm.E0 = p.p), period = getPeriod(p.prob, x, @set par_tm.E0 = p.p))),
	plotSolution = (x, p; k...) ->
		begin
			xtt = BK.getTrajectory(p.prob, x, @set par_tm.E0 = p.p)
			plot!(xtt; legend = false, k...);
			plot!(br,subplot=1, putspecialptlegend = false)
			# plot!(br_potrap,subplot=1, putspecialptlegend = false)
		end,
	callbackN = (x, f, J, res, iteration, itlinear, options; kwargs...) -> (return res<1e16),
	finaliseSolution = (z, tau, step, contResult; prob=nothing, k...) ->
		begin
			isnothing(prob) && return true
			T = getPeriod(prob, z.u, @set par_tm.E0 = z.p)
			@show T
			T < 20.
		end,
	normC = norminf)

plot(br, br_popsh, markersize=3)

_po = br_popsh.sol[1].x
fl = br_popsh.functional.flow


text:
struct GPModel{TGP} <: AbstractMCMC.AbstractModel
    gp::TGP
end

struct GPSampler{TS} <: AbstractMCMC.AbstractSampler
    sampler::TS
end

function StatsBase.sample(model::MCGP, nSamples::Int; kwargs...)
    return sample(Random.GLOBAL_RNG, model, nSamples; kwargs...)
end
function StatsBase.sample(rng::Random.AbstractRNG, model::MCGP, nSamples::Int; kwargs...)
    sampler = inference(model)
    thinning = get(kwargs, :thinning, sampler.thinning)
    discard_initial = get(kwargs, :discard_initial, sampler.nBurnin)
    progressname = get(kwargs, :progressname, "Sampling with $sampler")
    return sample(
        rng,
        GPModel(model),
        GPSampler(sampler),
        nSamples;
        thinning,
        discard_initial,
        progressname,
        kwargs...,
    )
end
function AbstractMCMC.step(
    rng, model::GPModel, sampler::GPSampler{<:GibbsSampling}; kwargs...
)
    compute_kernel_matrices!(model.gp)
    sample_local!(likelihood(model.gp), yview(model.gp), means(model.gp))
    f = sample_global!.(∇E_μ(model.gp), ∇E_Σ(model.gp), Zviews(model.gp), getf(model.gp))
    sampler.sampler.nIter += 1
    return f, nothing
end

function AbstractMCMC.step(
    rng, model::GPModel, sampler::GPSampler{<:GibbsSampling}, state; kwargs...
)
    sample_local!(likelihood(model.gp), yview(model.gp), means(model.gp))
    f = sample_global!.(∇E_μ(model.gp), ∇E_Σ(model.gp), Zviews(model.gp), getf(model.gp))
    sampler.sampler.nIter += 1
    return f, nothing
end

function AbstractMCMC.bundle_samples(
    samples, model::GPModel, sampler::GPSampler, state, chain_type::Type; kwargs...
)
    resume = get(kwargs, :cat, true)
    if !resume || isempty(sampler.sampler.sample_store) # check if either one should restart sampling or if no samples was ever taken
        sampler.sampler.sample_store = samples
    else
        sampler.sampler.sample_store = vcat(sampler.sampler.sample_store, samples)
    end
    return samples
end


text:
# Undergraduate Student: Arturo Burgos
# Professor: João Rodrigo Andrade
# Federal University of Uberlândia - UFU, Fluid Mechanics Laboratory - MFLab, Block 5P, Uberlândia, MG, Brazil

using Printf # Just the same as import a module
using Statistics
using CPUTime # Look for Pkg.add("CPUTime")

println("Choose a N value: ")
N = chomp(readline()) # Here I get the char from the user
n = parse(Int32,N)  # Here I convert the char into an integer(32bits)

A = rand(1:20000,n,n,3)

println("\n\nThe A[1,1,1] is: ",A[1,1,1])
println("\nThe A[1,1,2] is: ",A[1,1,2])
println("\nThe A[1,1,3] is: ",A[1,1,3])
println("\n\nThe A[n,n,1] is: ",A[n,n,1])
println("The A[n,n,2] is: ",A[n,n,2])
println("The A[n,n,3] is: ",A[n,n,3])

function vectorized(mtx)

    mtx[:, :, 1] = mtx[:, :, 2]
    mtx[:, :, 3] = mtx[:, :, 1]
    sleep(1)
    return mtx # retornar algo à variavel alocada, no caso "A"
end

A = vectorized(A)

println("\n\nThe new A[1,1,1] is: ",A[1,1,1])
println("\nThe new A[1,1,2] is: ",A[1,1,2])
println("\nThe new A[1,1,3] is: ",A[1,1,3])
println("\n\nThe A[n,n,1] is: ",A[n,n,1])
println("The A[n,n,2] is: ",A[n,n,2])
println("The A[n,n,3] is: ",A[n,n,3])

@time @CPUtime vectorized(A)

# @time, @elapsed, time_ns() ,time()

max_stars_repo_name: bangyen/Comparing-Languages

max_issues_repo_name: bangyen/Comparing-Languages

max_forks_repo_name: bangyen/Comparing-Languages


text:
module TestMapreduce
include("preamble.jl")
using Transducers: air

@testset "mapreduce" begin
    xf = opcompose(Filter(isodd), Map(inc))

    ans = @test_deprecated mapreduce(xf, +, 0:4, init = 0)
    DEPWARN_ERROR || @test ans == 6

    ans = @test_deprecated mapreduce(
        xf,
        +,
        air.([1:5;] .- 1);
        init = 0,
        basesize = typemax(Int),  # Broadcasted is not parallelizable yet
    )
    DEPWARN_ERROR || @test ans == 6

    ans = @test_deprecated mapreduce(xf, +, skipmissing(Any[0:4;]), init = 0)
    DEPWARN_ERROR || @test ans == 6

    ans = @test_deprecated mapreduce(xf, +, 1, init = 0)
    DEPWARN_ERROR || @test ans == 2
end

end  # module


text:
@testset "DBSCAN" begin

    Random.seed!(42)

    df = CSV.read("data/blob_data.csv", DataFrame, drop=[1]);
    X = df[:,1:2];

    ϵ = 0.35;
    min_pts = 10;

    @testset "test errors" begin
        @test_throws ArgumentError dbscan(X[:, 2], ϵ, min_pts)
    end

    @testset "test Array{Float64, 2}" begin
        X_array = Matrix(X)
        res_ar = dbscan(X_array, ϵ, min_pts)
        n_clusters = unique(res_ar.labels) |> length
        @test length(res_ar.labels) == size(X_array, 1)
    end

    @testset "test dimensions results" begin
        res = dbscan(X, ϵ, min_pts)
        n_clusters = unique(res.labels) |> length
        @test length(res.labels) == size(X, 1)
        @test length(res.clusters) == n_clusters - 1
        @test length(res.clusters[1]) == 1010
        @test length(res.clusters[2]) == 999
        @test length(res.clusters[3]) == 994
        @test length(res.clusters[4]) == 1000
    end

    @testset "test Integer conversion" begin
        X_int = round.(Int, X.*100)
        res_int = dbscan(X_int, ϵ*100, min_pts)
        @test eltype(res_int.df) == Float64
        @test length(res_int.labels) == size(X_int, 1)
    end

end