text:
import Base: iterate
using Base: SizeUnknown, EltypeUnknown

abstract type AbstractOrder end
struct PreOrder <: AbstractOrder end
struct PostOrder <: AbstractOrder end
struct LevelOrder <: AbstractOrder end

# S indicates whether a tuple in ts is a single tree or a collection of trees
struct PredicateIterator{S, T, O <: AbstractOrder}
    f::Function
    ts::T
    complete::Bool
    order::O
    function PredicateIterator(f::Function, ts::T, single::Bool; complete::Bool=true,
                                  order::O=PreOrder()) where {T, O <: AbstractOrder}
        new{single, T, O}(f, ts, complete, order)
    end
end

PredicateIterator(f::Function, t; kwargs...) = PredicateIterator(f, (t,), true; kwargs...)
PredicateIterator(f::Function, ts...; kwargs...) = PredicateIterator(f, ts, false; kwargs...)

_leaf_predicate(ns::Tuple) = all(n -> isnothing(n) || isleaf(n), ns)
_leaf_predicate(n) = isleaf(n)
_node_predicate(ns::Tuple) = true
_node_predicate(n) = true
function _type_predicate_object(ts::Tuple, t::Tuple{Vararg{Type}})
    @assert length(ts) == length(t)
    (ns::Tuple) -> all([ns[i] isa t[i] for i in eachindex(ts)])
end
function _type_predicate_object(ts::Tuple, t::Type{T}) where T
    (ns::Tuple) -> all(n -> n isa T, ns)
end

LeafIterator(ts...; kwargs...) = PredicateIterator(_leaf_predicate, ts...; kwargs...)
NodeIterator(ts...; kwargs...) = PredicateIterator(_node_predicate, ts...; kwargs...)
function TypeIterator(t::Union{Type, Tuple{Vararg{Type}}}, ts...; kwargs...)
    PredicateIterator(_type_predicate_object(ts, t), ts...; kwargs...)
end
TypeIterator(t::Type{T}, ts; kwargs...) where T = PredicateIterator(n -> n isa T, ts; kwargs...)

Base.IteratorSize(::PredicateIterator) = SizeUnknown()
Base.IteratorEltype(::PredicateIterator) = EltypeUnknown()

# TODO implement with stack?
# TODO maybe implement faster version with types
traverse!(o::AbstractOrder, complete::Bool, ts::Nothing, res::Vector) = res

function traverse!(o::AbstractOrder, complete::Bool, ts::Tuple, res::Vector)
    if !all(_isnothing_iter(ts))
        o isa PreOrder && push!(res, ts)
        for chs in _children_pairs(ts, complete) |> _iter
            traverse!(o, complete, chs, res)
        end
        o isa PostOrder && push!(res, ts)
    end
    res
end

function traverse!(o::LevelOrder, complete::Bool, ts::Tuple, res::Vector)
    q = Queue{Any}()
    enqueue!(q, ts)
    while !isempty(q)
        ts = first(q)
        if !all(_isnothing_iter(ts))
            push!(res, ts)
            for chs in _children_pairs(ts, complete) |> _iter
                enqueue!(q, chs)
            end
        end
        dequeue!(q)
    end
    res
end

Base.iterate(it::PredicateIterator{S, Tuple{}}) where S = nothing

function Base.iterate(it::PredicateIterator)
    ns = []
    traverse!(it.order, it.complete, it.ts, ns)
    return iterate(it, (1, ns))
end

function Base.iterate(it::PredicateIterator{true}, (i, ns))
    while i <= length(ns) && !it.f(only(ns[i]))
        i += 1
    end
    i <= length(ns) || return nothing
    return only(ns[i]), (i+1, ns)
end

function Base.iterate(it::PredicateIterator{false}, (i, ns))
    while i <= length(ns) && !it.f(ns[i])
        i += 1
    end
    i <= length(ns) || return nothing
    return ns[i], (i+1, ns)
end

# TODO implement MultiIterator, so that nothing below doesn't clash with complete iterators
# struct MultiIterator
#     its::Tuple
#     MultiIterator(ts...) = new(ts)
#     MultiIterator(ts) = new(ts)
# end

# function iterate(it::MultiIterator)
#     r = collect(map(iterate, it.its))
#     any(isnothing.(r)) && return nothing
#     return collect(zip(r...))
# end

# function iterate(it::MultiIterator, ss)
#     r = [iterate(i, s) for (i, s) in zip(it.its, ss)]
#     any(isnothing.(r)) && return nothing
#     return collect(zip(r...))
# end

# defaultstack(it::IteratorTypes) = Any[it.t]
# defaultstack(it::ZipIterator) = [Any[t] for t in it.ts]

# function iterate(it::PredicateIterator, s=Stack{Tuple}(it.ts))
#     r = nextstate(it, s)
#     isnothing(r) && return nothing
#     return r, s
# end

# expand(n, s) = append!(s, reverse(collect(_children_sorted(n))))

# function nextstate(it, s)
#     isempty(s) && return nothing
#     n = pop!(s)
#     processnode(it, n, s)
# end

# function nextstate(it::ZipIterator, ss)
#     any(isempty.(ss)) && return nothing
#     ns = pop!.(ss)
#     if !any(isleaf.(ns))
#         for (n, s) in zip(ns, ss)
#             expand(n, s)
#         end
#     end
#     tuple(ns...)
# end

# function processnode(it::NodeIterator, n, s)
#     expand(n, s)
#     n
# end

# processnode(it::LeafIterator, n, s) = processnode(NodeType(n), it, n, s)
# processnode(::LeafNode, it, n, s) = n
# processnode(_, it, n, s) = nextstate(it, expand(n, s))

# function processnode(it::PredicateIterator, n, s)
#     expand(n, s)
#     it.f(n) ? n : nextstate(it, s)
# end

# function processnode(it::TypeIterator{T, U}, n::T, s) where {T, U}
#     expand(n, s)
#     n
# end
# processnode(it::TypeIterator, n, s) = nextstate(it, expand(n, s))


text:
using CUDA, GPUArrays
using NiLang, NiLang.AD

"""
A reversible swap kernel for GPU for SWAP gate in quantum computing.
See the irreversible version for comparison

http://tutorials.yaoquantum.org/dev/generated/developer-guide/2.cuda-acceleration/
"""
@i @inline function swap_kernel(state::AbstractVector{T}, mask1, mask2) where T
    @invcheckoff b ← (blockIdx().x-1) * blockDim().x + threadIdx().x
    @invcheckoff if (b < length(state), ~)
        if (b&mask1==0 && b&mask2==mask2, ~)
            NiLang.SWAP(state[b+1], state[b ⊻ (mask1|mask2) + 1])
        end
    end
end

# TODO: support ::Type like argument.
"""
SWAP gate in quantum computing.
"""
@i function instruct!(state::CuVector, gate::Val{:SWAP}, locs::Tuple{Int,Int})
    mask1 ← 1 << (locs[1]-1)
    mask2 ← 1 << (locs[2]-1)
    @cuda threads=256 blocks=ceil(Int,length(state)/256) swap_kernel(state, mask1, mask2)
end

using Test
@testset "swap gate" begin
    v = cu(randn(128))
    v1 = instruct!(copy(v), Val(:SWAP), (3,4))[1]
    v2 = instruct!(copy(v1), Val(:SWAP), (3,4))[1]
    v3 = (~instruct!)(copy(v1), Val(:SWAP), (3,4))[1]
    @test !(v ≈ v1)
    @test v ≈ v2
    @test v ≈ v3
end

@i function loss(out!, state::CuVector)
    instruct!(state, Val(:SWAP), (3,4))
    out! += state[4]
end

loss(0.0, CuArray(randn(128)))
Grad(loss)(Val(1), 0.0, CuArray(randn(128)))

####################### A different loss ###############
@i function loss(out!, state::CuVector, target::CuVector)
    instruct!(state, Val(:SWAP), (3,4))
    out! += state' * target
end

# requires defining a new primitive, we don't how to parallelize a CUDA program automatically yet.
using LinearAlgebra: Adjoint
function (_::MinusEq{typeof(*)})(out!::GVar, x::Adjoint{<:Any, <:CuVector{<:GVar}}, y::CuVector{<:GVar})
    chfield(out!, value, value(out!)-(value.(x) * value.(y))[]),
    chfield.(parent(x), grad, grad.(parent(x)) .+ grad(out!)' .* conj.(value.(y)))',
    chfield.(y, grad, grad.(y) .+ grad(out!) .* conj.(value.(x')))
end

function (_::PlusEq{typeof(*)})(out!::GVar, x::Adjoint{<:Any, <:CuVector{<:GVar}}, y::CuVector{<:GVar})
    chfield(out!, value, value(out!)+(value.(x) * value.(y))[]),
    chfield.(parent(x), grad, grad.(parent(x)) .- grad(out!)' .* conj.(value.(y)))',
    chfield.(y, grad, grad.(y) .- grad(out!) .* conj.(value.(x')))
end

function (_::PlusEq{typeof(*)})(out!, x, y)
    out! += x * y
    out!, x, y
end

function (_::MinusEq{typeof(*)})(out!, x, y)
    out! -= x * y
    out!, x, y
end

loss(0.0, CuArray(randn(128)), CuArray(randn(128)))
Grad(loss)(Val(1), 0.0, CuArray(randn(128)), CuArray(randn(128)))

ext: jl


text:
using Roassal
using Random

function example01()
    canvas = RCanvas()
    box = RBox()
    add!(canvas, box)
    extent!(box, 25, 15)
    translate_to!(box, -20, -20)
    rshow(canvas)
end

function example02()
    canvas = RCanvas()
    Random.seed!(42)
    for x in 1:3
        for y in 1:3
            box = RBox(color=RColor(rand(), rand(), rand()))
            translate_to!(box, x * 20, y * 20)
            add!(canvas, box)
        end
    end
    rshow(canvas)
end

function example03()
    canvas = RCanvas()
    Random.seed!(42)
    for x in 1:3
        for y in 1:3
            box = RBox(color=RColor(rand(), rand(), rand()))
            translate_to!(box, x * 20, y * 20)
            add!(canvas, box)
        end
    end
    highlightable(get_shapes(canvas))
    #add_callback!(get_shapes(canvas), Callback(:mouseClick, (event, shape) -> remove_shape!(canvas, shape)))
    #highlightable(get_shapes(canvas)[1])
    #add_callback!(firstBox, Callback(:mouseMove, (event, shape) -> println("CALLBACK mouse move $(shape)")))
    rshow(canvas)
end

function example04()
    c = RCanvas()
    box1 = RBox(;width=30, height=30)
    box2 = RBox(;width=30, height=30)
    add!(c, box1)
    add!(c, box2)
    translate_to!(box2, (40, 0))

    set_color!(box1, RColor_RED)
    set_color!(box2, RColor_GREEN)

    highlightable(box1)
    highlightable(box2)
    popup(box2)
    rshow(c)
end

function example05()
    c = RCanvas()
    box1 = RBox()
    box2 = RBox()

    add!(c, box1)
    add!(c, box2)

    translate_to!(box2, 50, 30)

    line = RLine(box1, box2)
    add!(c, line)

    rshow(c)
end

example05()

print("Press Enter to exit")
readline()


text:
import Base
import PyPlot
export show

function show(stat::Array{CPUStat,1})

    for i = 1:length(stat);
        ax = PyPlot.gca()
        ax[:set_ylim]([0,100])
        PyPlot.plot(stat[i].cpu)
        PyPlot.title(stat[i].pid)
    end

end


text:
using BenchmarkTools

fullstrip(str::AbstractString, chars::AbstractString) = replace(str, Regex("[$chars]") => "")

punctuation_characters = ASCIIString("~\`!@#$%^&*()_-+=[{]}\\|;:',<.>/?1234567890")
#punctuation_characters = '!"#$%&\'()*+,-./:;<=>?@[\\]^_`{|}~ \t\n\r\x0b\x0c'

"""
    countwords(filename::String)

Opens the given file and makes a collection of its unique words.
"""
function countwords(filename::String)
    openedfile = open(filename)
    wordlist = String[]
    for line in eachline(openedfile)
        println(line)
        words = split( fullstrip(line, punctuation_characters))
        map(word -> push!(wordlist, lowercase(word)), words)
    end
    close(openedfile)
    filter!(!isempty, wordlist)
    println("Number of unique words: ", length(wordlist))
end

# Get the text file name from the command line.
n, = size(ARGS)
if n < 1
    println("Usage: ")
    println("    julia " + ARG[0] + " filename")
    println("       --->      Please specify the filename.")
    exit()
end

filename = ARGS[1]
@btime count_word(filename)

println(" ")


text:
using TemporalGPs:
    x0,
    fill_in_missings,
    replace_observation_noise_cov,
    transform_model_and_obs

println("missings:")
@testset "missings" begin

    rng = MersenneTwister(123456)

    storages = (
        dense=(name="dense storage Float64", val=ArrayStorage(Float64)),
        static=(name="static storage Float64", val=SArrayStorage(Float64)),
    )
    emission_types = (
        small_output=(name="small output", val=SmallOutputLGC),
        large_output=(name="large output", val=LargeOutputLGC),
        scalar_output=(name="scalar output", val=ScalarOutputLGC),
    )
    settings = [
        (tv=:time_varying, N=5, Dlat=3, Dobs=2, storage=storages.dense),
        (tv=:time_varying, N=5, Dlat=3, Dobs=2, storage=storages.static),
        (tv=:time_invariant, N=5, Dlat=3, Dobs=2, storage=storages.dense),
        (tv=:time_invariant, N=5, Dlat=3, Dobs=2, storage=storages.static),
    ]
    orderings = [
        Forward(),
    ]

    # We can test this stuff by analytically dropping particular elements of the LGSSM
    # and computing the marginalised dynamics.
    @testset "($tv, $N, $Dlat, $Dobs, $(storage.name), $(emission.name), $order)" for
        (tv, N, Dlat, Dobs, storage) in settings,
        emission in emission_types,
        order in orderings

        # Print current iteration to prevent CI timing out.
        println(
            "(time_varying=$tv, N=$N, Dlat=$Dlat, Dobs=$Dobs, " *
            "storage=$(storage.name), emissions=$(emission.val), ordering=$order)",
        )

        # Build LGSSM.
        model = if emission.val ∈ (SmallOutputLGC, LargeOutputLGC)
            random_lgssm(
                rng, order, Val(tv), emission.val, Dlat, Dobs, N, Val(:dense), storage.val,
            )
        elseif emission.val == ScalarOutputLGC
            random_lgssm(rng, order, Val(tv), emission.val, Dlat, Dobs, N, storage.val)
        else
            throw(error("Unrecognised storage $(emission.val)"))
        end

        # Verify the correct output types has been obtained.
        @test eltype(model.emissions) <: emission.val

        # Verify that model properties are as requested.
        @test storage_type(model) == storage.val

        @test length(model) == N
        @test getindex(model, N) isa TemporalGPs.ElementOfLGSSM

        y = rand(model)
        x = TemporalGPs.x0(model)

        # Make some of the data missing.
        missings_idx = [3, 5, 9]
        missings_idx = [2, 4]
        presents_idx = setdiff(eachindex(y), missings_idx)
        y_missing = Vector{Union{eltype(y), Missing}}(undef, length(y))
        y_missing[missings_idx] .= missing
        y_missing[presents_idx] .= y[presents_idx]

        # Construct a modified LGSSM that is shorter than the original with analytically-
        # marginalised dynamics where there are missing data.
        As = model.transitions.As
        as = model.transitions.as
        Qs = model.transitions.Qs
        transitions_missing = map(eachindex(y_missing)) do n
            if n - 1 ∈ missings_idx
                A = As[n] * As[n-1]
                a = As[n] * as[n-1] + as[n]
                Q = As[n] * Qs[n-1] * As[n]' + Qs[n]
                return A, a, Q
            else
                return As[n], as[n], Qs[n]
            end
        end

        new_As = getindex.(transitions_missing, 1)[presents_idx]
        new_as = getindex.(transitions_missing, 2)[presents_idx]
        new_Qs = getindex.(transitions_missing, 3)[presents_idx]
        new_model = LGSSM(
            GaussMarkovModel(order, new_As, new_as, new_Qs, x0(model)),
            model.emissions[presents_idx],
        )
        new_y = y[presents_idx]

        @testset "logpdf" begin
            @test logpdf(new_model, new_y) ≈ logpdf(model, y_missing)
        end

        @testset "filter" begin
            marginal_filters = _filter(new_model, new_y)
            filters = _filter(model, y_missing)
            @test mean.(marginal_filters) ≈ mean.(filters)[presents_idx]
            @test cov.(marginal_filters) ≈ cov.(filters)[presents_idx]
        end

        @testset "posterior" begin
            new_posterior = posterior(new_model, new_y)
            post = posterior(model, y_missing)

            new_post_marginals = marginals(new_posterior)
            post_marginals = marginals(post)
            @test mean.(new_post_marginals) ≈ mean.(post_marginals)[presents_idx]
            @test cov.(new_post_marginals) ≈ cov.(post_marginals)[presents_idx]

            @test logpdf(new_posterior, new_y) ≈ logpdf(post, y_missing) rtol=1e-4
        end

        # Only test the bits of AD that we haven't tested before.
        @testset "AD: transform_model_and_obs" begin
            fdm = central_fdm(2, 1)
            adjoint_test(fill_in_missings, (model.emissions.Q, y_missing); fdm=fdm)
            adjoint_test(replace_observation_noise_cov, (model, model.emissions.Q))
            adjoint_test(transform_model_and_obs, (model, y_missing); fdm=fdm)
        end
    end

    storages = (
        dense=(name="dense storage Float64", val=ArrayStorage(Float64)),
    )
    emission_types = (
        small_output=(name="small output", val=SmallOutputLGC),
        large_output=(name="large output", val=LargeOutputLGC),
    )
    settings = [
        (tv=:time_varying, N=5, Dlat=3, Dobs=2, storage=storages.dense),
        (tv=:time_invariant, N=5, Dlat=3, Dobs=2, storage=storages.dense),
    ]
    orderings = [
        Forward(),
    ]

    @testset "missing-in-obs (tv=$(tv), N=$N, Dlat=$Dlat, emission=$emission" for
        (tv, N, Dlat, Dobs, storage) in settings,
        emission in emission_types,
        order in orderings

        println("missing-in-obs (tv=$(tv), N=$N, Dlat=$Dlat, emission=$emission")

        # Build LGSSM.
        model = random_lgssm(
            rng, order, Val(tv), emission.val, Dlat, Dobs, N, Val(:diag), storage.val,
        )

        # Verify the correct output types has been obtained.
        @test eltype(model.emissions) <: emission.val

        # Verify that model properties are as requested.
        @test storage_type(model) == storage.val

        @test length(model) == N
        @test getindex(model, N) isa TemporalGPs.ElementOfLGSSM

        # Generate some missing data.
        y = rand(model)
        y_missing = map(y) do yn
            yn_missing = Vector{Union{Missing, eltype(yn)}}(undef, length(yn))
            yn_missing .= yn
            yn_missing[randperm(length(yn))[1]] = missing
            return yn_missing
        end

        # Check logpdf and inference run, infer, and play nicely with AD.
        @inferred logpdf(model, y_missing)
        adjoint_test(y_missing -> logpdf(model, y_missing), (y_missing, ))
        @inferred posterior(model, y_missing)
    end
end


text:
Base.promote_rule(::Type{Decimal}, ::Type{<:Real}) = Decimal

# override definitions in Base
Base.promote_rule(::Type{BigFloat}, ::Type{Decimal}) = Decimal
Base.promote_rule(::Type{BigInt}, ::Type{Decimal}) = Decimal

# Addition
# To add, convert both decimals to the same exponent.
# (If the exponents are different, use the smaller exponent
# to make sure we're adding integers.)
function +(x::Decimal, y::Decimal)
    cx = (-1)^x.s * x.c * 10^max(x.q - y.q, 0)
    cy = (-1)^y.s * y.c * 10^max(y.q - x.q, 0)
    s = (abs(cx) > abs(cy)) ? x.s : y.s
    c = BigInt(cx) + BigInt(cy)
    normalize(Decimal(s, abs(c), min(x.q, y.q)))
end

# Negation
-(x::Decimal) = Decimal((x.s == 1) ? 0 : 1, x.c, x.q)

# Subtraction
-(x::Decimal, y::Decimal) = +(x, -y)

# Multiplication
function *(x::Decimal, y::Decimal)
    s = (x.s == y.s) ? 0 : 1
    normalize(Decimal(s, BigInt(x.c) * BigInt(y.c), x.q + y.q))
end

# Inversion
function Base.inv(x::Decimal)
    str = string(x)
    if str[1] == '-'
        str = str[2:end]
    end
    b = ('.' in str) ? length(split(str, '.')[1]) : 0
    c = round(BigInt(10)^(-x.q + DIGITS) / x.c)
    q = (x.q < 0) ? 1 - b - DIGITS : -b - DIGITS
    normalize(Decimal(x.s, c, q))
end

# Division
/(x::Decimal, y::Decimal) = x * inv(y)

# Exponentiation, please test

function ^(x::Decimal, y::Decimal)
    #(-1.23)^(something) code to catch errors prior to the routine
    x.s == 1 && begin
        (y.q < 0  && throw(DomainError(y),"Exponentiation yielding a complex result requires a complex argument."))
    end
    y.s == 1 && (return ^(inv(x),-y))
    y.q < 0 && (return Decimal(BigFloat(x)^BigFloat(y))) #add correction to precision, its too much

    if iseven(y.c) #squared number
        zs = 0
    else
        if x.s == 0
            zs =0
        else
            zs = 1
        end
    end

        zc = x.c^((y.c)*(10^y.q))
        zq = (x.q)*((y.c)*(10^y.q))
        return  Decimal(zs,zc,zq)

end


text:
export Attribute, Attributes, Fill, Fills, Stroke, Strokes, Linestyle, Linewidth, Linewidths, Markersize, Markersizes, Marker, Font, Visible
export Clip, Textfill, Opacity, Operator, Invisible

Base.Broadcast.broadcastable(a::Attribute) = Ref(a)

Fill(s::String) = Fill(parse(Colors.Colorant, s))
Fill(s::String, r::Real) = Fill(Colors.RGBA(parse(Colors.Colorant, s), r))

Textfill(s::String) = Textfill(parse(Colors.Colorant, s))

Stroke(s::String) = Stroke(parse(Colors.Colorant, s))
Stroke(s::String, r::Real) = Stroke(Colors.RGBA(parse(Colors.Colorant, s), r))

Attributes() = Attributes(Dict{Type{<:Attribute}, Union{Function, Attribute}}())

Operator() = Operator(:over)

const Invisible = Visible(false)

parameterlesstypeof(a::Attribute) = typeof(a).name.wrapper

function Attributes(att::Attribute, varargs::Vararg{Attribute, N}) where N
    attributes = Attributes()
    for a in (att, varargs...)
        insert!(attributes, a)
    end
    attributes
end

function Base.insert!(attrs::Attributes, a::Attribute)
    pltype = parameterlesstypeof(a)
    attrs[pltype] = a
end

Base.haskey(a::Attributes, key) = haskey(a.attrs, key)
Base.setindex!(a::Attributes, value, key) = setindex!(a.attrs, value, key)
Base.getindex(a::Attributes, key) = getindex(a.attrs, key)

ismultiattr(a::Attribute) = false
ismultiattr(f::Fill{T}) where T <: Union{Array, Pair{Symbol, <:Function}} = true
ismultiattr(s::Stroke{T}) where T <: Union{Array, Pair{Symbol, <:Function}} = true


text:
using BinaryBuilder

name = "ZeroMQ"
version = v"4.3.2"

# Collection of sources required to build ZMQ
sources = [
    GitSource("https://github.com/zeromq/libzmq.git", "a84ffa12b2eb3569ced199660bac5ad128bff1f0"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libzmq
if [[ "${target}" == *-mingw* ]]; then
    # Apply patch from
    # https://github.com/msys2/MINGW-packages/blob/350ace4617661a4df7b9474c573b08325fa716c3/mingw-w64-zeromq/001-mingw-__except-fixes.patch
    atomic_patch -p1 ../patches/001-mingw-__except-fixes.patch
elif [[ "${target}" == *86*-linux-musl* ]]; then
    pushd /opt/${target}/lib/gcc/${target}/*/include
    # Fix bug in Musl C library, see
    # https://github.com/JuliaPackaging/BinaryBuilder.jl/issues/387
    atomic_patch -p0 $WORKSPACE/srcdir/patches/mm_malloc.patch
    popd
fi
sh autogen.sh
./configure --prefix=$prefix \
    --host=${target} \
    --without-docs \
    --enable-drafts \
    --with-libsodium \
    --disable-libunwind \
    --disable-perf \
    --disable-Werror \
    --disable-eventfd \
    --without-gcov \
    --disable-static \
    CXXFLAGS="-O2 -fms-extensions"
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libzmq", :libzmq),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
    Dependency("libsodium_jll")
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
# Slices
if write_slices
fields = Dict(
    "u" => model.velocities.u,
    "v" => model.velocities.v,
    "w" => model.velocities.w,
    "b" => model.tracers.b
)

surface_output_writer =
    NetCDFOutputWriter(model, fields, filename= filename_1 * "_surface.nc",
			           time_interval= slice_output_interval, zC=Nz, zF=Nz)

middepth_output_writer =
    NetCDFOutputWriter(model, fields, filename= filename_1 * "_middepth.nc",
                       time_interval = slice_output_interval, zC=Int(Nz/2), zF=Int(Nz/2))

zonal_output_writer =
    NetCDFOutputWriter(model, fields, filename= filename_1 * "_zonal.nc",
                       time_interval=slice_output_interval, yC=Int(Ny/2), yF=Int(Ny/2))

meridional_output_writer =
    NetCDFOutputWriter(model, fields, filename= filename_1 * "_meridional.nc",
                       time_interval=slice_output_interval, xC=Int(Nx/2), xF=Int(Nx/2))
end
## Zonal and Time Averages
if write_zonal
u, v, w = model.velocities
b = getproperty(model.tracers, :b)

u_scratch = XFaceField(model.architecture, model.grid)
v_scratch = YFaceField(model.architecture, model.grid)
w_scratch = ZFaceField(model.architecture, model.grid)
b_scratch =  CellField(model.architecture, model.grid)

zonal_fields = Dict(
    :u => AveragedField( u, dims=(1,)),
    :v => AveragedField( v, dims=(1,)),
    :w => AveragedField( w, dims=(1,)),
    :b => AveragedField( b, dims=(1,)),

    :uu => AveragedField(u * u, dims= 1),
    :vv => AveragedField(v * v, dims= 1),
    :ww => AveragedField(w * w, dims= 1),
    :uv => AveragedField(u * v, dims= 1),
    :vw => AveragedField(w * v, dims= 1),
    :uw => AveragedField(w * u, dims= 1),
    :ub => AveragedField(u * b, dims= 1),
    :vb => AveragedField(v * b, dims= 1),
    :wb => AveragedField(w * b, dims= 1),
)

zonal_statistics = JLD2OutputWriter(model, zonal_fields,
                     time_averaging_window = time_avg_window,
                             time_interval = zonal_output_interval,
                                    prefix = filename_1 * "_zonal_averages",
                                     force = true)
end

## Checkpointer
checkpointer = Checkpointer(model, prefix = filename_1 * "_checkpoint", time_interval = checkpoint_interval, force = true)

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
_class_params = """    outlier_fraction::Float64
The fraction of outliers (number between 0 and 1) in the data is used to determine the score threshold to classify the
samples into inliers and outliers.

    classify::Union{Function, Nothing}
A function to transform a vector of scores to a vector of bits, where 1 represents an outlier and 0 represents a normal
instance. *Hint:* Sometimes you don't want to transform your scores to classes, e.g. in ROC AUC evaluation, where you
can use `no_classify` to pass along the reduced (raw) scores. See [`classify`](@ref) for a specific implementation."""

_score_params = """    combine::Function
A function to reduce a matrix, where each row represents an instance and each column represents a score of specific
detector, to a vector of scores for each instance. See `combine` for a specific implementation. *Note:* This function
is not called if the input to the evaluator consists of a single train/test scores tuple.

    normalize::Union{Function, Nothing}
A function to reduce a matrix, where each row represents an instance and each column a score of specific detector, to a
vector of scores for each instance. See [`normalize`](@ref) for a specific implementation."""

"""
    Class(outlier_fraction = 0.1,
          classify = classify,
          combine = combine,
          normalize = normalize)

A flexible, quantile-thresholding classifier that maps the outlier scores of a single or multiple outlier detection
models to binary classes, where `1` represents inliers and `-1` represents outliers.

Parameters
----------
$_class_params

$_score_params

Examples
----------
$_evaluator
"""
MMI.@mlj_model mutable struct Class <: Evaluator
    outlier_fraction::Float64 = 0.1::(0 < _ < 1)
    classify::Function = classify
    combine::Function = combine
    normalize::Union{Nothing, Function} = normalize
end

"""
    Score(combine = combine,
          normalize = normalize)

Transform the results of a single or multiple outlier detection models to combined and normalized scores.

Parameters
----------
$_score_params

Examples
----------
```julia
using OutlierDetection: Score, KNN, fit, score
detector = KNN()
X = rand(10, 100)
model = fit(detector, X)
train_scores, test_scores = score(detector, model, X)
ŷ = detect(Score(), train_scores, test_scores)
```
"""
MMI.@mlj_model mutable struct Score <: Evaluator
    combine::Function = combine
    normalize::Union{Nothing, Function} = normalize
end

function detect(ev::Class, scores::Result...)::Labels
    _detect_helper(ev.outlier_fraction, ev.classify, ev.normalize, ev.combine, scores...)
end

function detect(ev::Score, scores::Result...)::Scores
    _detect_helper(nothing, nothing, ev.normalize, ev.combine, scores...)
end

function _detect_helper(outlier_fraction::Union{Real,Nothing},
                        classify::Union{Function,Nothing},
                        normalize::Union{Function,Nothing},
                        combine::Function,
                        scores::Result...)::Scores
    # transforms a variable number of equal-length scores into classes, where each input tuple represents the
    # train scores and test scores of a detector.

    # make sure that we have at least one tuple of train and test scores
    n_scores = length(scores)
    @assert n_scores > 0

    # conditionally normalize all scores if not nothing
    scores = isnothing(normalize) ? scores : map(tup -> normalize(tup...), scores)

    # [(train1, test1), (train2, test2)] -> matrix of scores for train [train1 train2]' and test [test1 test2]'
    # Note: The matrices contain one observation per row!
    reduce_cat = idx -> reduce(hcat, getfield.(scores, idx))
    scores_train, scores_test = reduce_cat(1), reduce_cat(2)

    # return an identity function that returns the existing scores if classify is nothing
    classify = isnothing(classify) ? (_, _, scores) -> scores : classify

    # scores_train and scores_test can be either a vector (if n_scores == 1) or a matrix (if n_scores > 1)
    n_scores == 1 ? classify(outlier_fraction, scores_train, scores_test) :
                    classify(outlier_fraction, combine(scores_train), combine(scores_test))
end


text:
struct UselessTable
        mat::Matrix{Any}
        heads::Vector{Any}
        foots::Vector{Any}
        useful::Any
end

import Base.show
function show(io::IO, tb::UselessTable)
        compact = get(io, :compact, false)
        dsize = get(io, :displaysize, displaysize())
        limit = get(io, :limit, true)

        offset = "  "
        nrow, ncol = size(tb.mat, 1), size(tb.mat, 2)
        mat = fill("", (nrow, ncol))
        for i in 1:nrow, j in 1:ncol
                if ismissing(tb.mat[i,j]) mat[i,j] = "N/A"
                elseif isnothing(tb.mat[i,j]) mat[i,j] = ""
                else mat[i,j] = string(tb.mat[i,j])
                end
                #todo compact
        end
        colwidth = [maximum(map(textwidth, mat[:,i])) for i in 1:ncol] .+ 1
        for i in 1:nrow, j in 1:ncol
                pad = ifelse(j == 1, rpad, lpad)
                val = mat[i,j]
                val = pad(val, colwidth[j] - textwidth(val) + length(val))
                if j == 1 val = offset * val end
                mat[i,j] = val
        end

        function middle(s::String, width::Int)
                r = offset * s
                n = textwidth(r)
                if n < width r = " " ^ div((width - n), 2) * r end
                return r
        end
        function right(s::String, width::Int)
                r = offset * s
                n = textwidth(r)
                if n < width r = " " ^ (width - n) * r end
                return r
        end

        totw = sum(colwidth)
        for i in 1:nrow
                #print headers
                if i == 1 && length(tb.heads) > 0
                        for head in tb.heads
                                println(io, middle(string(head), totw))
                        end
                        println(io, middle("-" ^ totw, totw))
                end
                for j in 1:ncol
                        print(io, mat[i,j])
                end
                println(io)
                #print foots
                if i == nrow && length(tb.foots) > 0
                        println(io, middle("=" ^ totw, totw))
                        for foot in tb.foots
                                println(io, right(string(foot), totw))
                        end
                end
        end
end

function UselessTable(mat::Matrix{Tv}; rnames=nothing, cnames=nothing, topleft=nothing, heads=[], foots=[], useful=nothing) where {Tv}
        ncol = size(mat, 2); nrow = size(mat, 1)
        if isnothing(cnames) cnames = 1:ncol end
        if isnothing(rnames) rnames = 1:nrow end
        ncol += 1; nrow += 1
        mat2 = Matrix{Any}(missing, (nrow, ncol))
        mat2[1,1] = topleft
        mat2[1,2:end] .= cnames
        mat2[2:end,1] .= rnames
        for i in 2:nrow, j in 2:ncol
                mat2[i,j] = mat[i-1,j-1]
        end
        return UselessTable(mat2, heads, foots, useful)
end

==(a::UselessTable, b::UselessTable) = a.mat == b.mat

#=
function UselessTable(dict::Dict{Tcname,Dict{Trname, Tv}}; topleft="") where {Tcname, Trname, Tv}
        cnames = Dict{Tcname, Int}(map(reverse, enumerate(keys(dict))))
        rnames = Dict{Trname, Int}()
        for col in values(dict), (rname, cell) in col
                if !haskey(rnames, rname) rnames[rname] = length(rnames) + 1 end
        end
        rnames2 = sort(collect(keys(rnames)), by=x->rnames[x])
        cnames2 = sort(collect(keys(cnames)), by=x->cnames[x])
        ncol = length(cnames) + 1
        nrow = length(rnames) + 1
        mat = Matrix{Any}(missing, (nrow, ncol))
        mat[1,1] = topleft
        mat[1,2:end] .= cnames2
        mat[2:end,1] .= rnames2
        for (cname, col) in dict, (rname, cell) in col
                mat[rnames[rname]+1,cnames[cname]+1] = cell
        end
        return UselessTable(mat, dict)
end
=#


text:
@testset "Contrast Tests" begin
    @testset "Single-Factor Contrast" begin
        @testset "Single Contrast" begin
            data = [AnovaValue(     "Total", 4617.60, 39),
                    AnovaResult("Treatment", 3497.60,  4, 874.40, 27.33, 0, 0), # p value and effect size unnecessary for this test
                    AnovaFactor(    "Error", 1120.00, 35,  32.00)]

            result = AnovaData(data, [AnovaFactor(r.name, r.ss, r.df, r.ms) for r ∈ data[3:3]], [5], 8, [4.00; 10.00; 11.00; 24.00; 29.00])

            expected = AnovaContrastResult(2640.2222, 35, 99.008334, 9.6573624e-12, 0.85954768)

            @test all(contrast(result, [1, 1, 1, 2, 2]) .≈ expected)
        end

        @testset "Difference Contrasts" begin
            data = [AnovaValue(     "Total", 54340.4, 39),
                    AnovaResult("Treatment", 26135.9,  3, 8711.97, 11.1199, 0, 0), # p value and effect size unnecessary for this test
                    AnovaFactor(    "Error", 28204.4, 36,  783.457)]
            result = AnovaData(data, [AnovaFactor(r.name, r.ss, r.df, r.ms) for r ∈ data[3:3]], [4], 10, [280.90, 240.75, 208.89, 239.96])

            expected = AnovaContrastResults([AnovaContrastResult(26044.011,   36, 24.931819,   1.53404354e-5, 0.63966842),
                                             AnovaContrastResult( 2665.05625, 36,  2.26777496, 0.140814326,   0.2434352),
                                             AnovaContrastResult( 9653.449,   36,  6.1608033,  0.017862323,   0.38226473)])

            @test all(differencecontrasts(result) ≈ expected)
        end

        @testset "Difference Contrasts (Reverse)" begin
            data = [AnovaValue(     "Total", 54340.4, 39),
                    AnovaResult("Treatment", 26135.9,  3, 8711.97, 11.1199, 0, 0), # p value and effect size unnecessary for this test
                    AnovaFactor(    "Error", 28204.4, 36,  783.457)]

            result = AnovaData(data, [AnovaFactor(r.name, r.ss, r.df, r.ms) for r ∈ data[3:3]], [4], 10, [280.90, 240.75, 208.89, 239.96])

            expected = AnovaContrastResults([AnovaContrastResult(  126.261778, 36,  0.120869854, 0.73011755,   0.057846875),
                                             AnovaContrastResult(26972.44225,  36, 22.9516466,   2.8424295e-5, 0.62396317),
                                             AnovaContrastResult(16120.225,    36, 10.28788115,  0.0028101042, 0.47144314)])

            @test all(differencecontrasts(result, true) ≈ expected)
        end

        @testset "Simple Contrasts (Default)" begin
            data = [AnovaValue(     "Total", 54340.4, 39),
                    AnovaResult("Treatment", 26135.9,  3, 8711.97, 11.1199, 0, 0), # p value and effect size unnecessary for this test
                    AnovaFactor(    "Error", 28204.4, 36,  783.457)]

            result = AnovaData(data, [AnovaFactor(r.name, r.ss, r.df, r.ms) for r ∈ data[3:3]], [4], 10, [280.90, 240.75, 208.89, 239.96])

            expected = AnovaContrastResults([AnovaContrastResult(16120.225, 36, 10.28788115, 0.0028101042, 0.47144314),
                                             AnovaContrastResult(51854.401, 36, 33.093329,   1.4877478e-6, 0.69207342),
                                             AnovaContrastResult(16760.836, 36, 10.6967172,  0.002368843,  0.47861035)])

            @test all(simplecontrasts(result) ≈ expected)
        end

        @testset "Simple Contrasts (Factor Level 3)" begin
            data = [AnovaValue(     "Total", 54340.4, 39),
                    AnovaResult("Treatment", 26135.9,  3, 8711.97, 11.1199, 0, 0), # p value and effect size unnecessary for this test
                    AnovaFactor(    "Error", 28204.4, 36,  783.457)]

            result = AnovaData(data, [AnovaFactor(r.name, r.ss, r.df, r.ms) for r ∈ data[3:3]], [4], 10, [280.90, 240.75, 208.89, 239.96])

            expected = AnovaContrastResults([AnovaContrastResult(51854.401, 36, 33.093329,  1.4877478e-6, 0.69207342),
                                             AnovaContrastResult(10150.596, 36,  6.4780811, 0.0153496883, 0.39051772),
                                             AnovaContrastResult( 9653.449, 36,  6.1608033, 0.017862323,  0.38226473)])

            @test all(simplecontrasts(result, 3) ≈ expected)
        end

        @testset "Repeated Contrasts" begin
            data = [AnovaValue(     "Total", 54340.4, 39),
                    AnovaResult("Treatment", 26135.9,  3, 8711.97, 11.1199, 0, 0), # p value and effect size unnecessary for this test
                    AnovaFactor(    "Error", 28204.4, 36,  783.457)]

            result = AnovaData(data, [AnovaFactor(r.name, r.ss, r.df, r.ms) for r ∈ data[3:3]], [4], 10, [280.90, 240.75, 208.89, 239.96])

            expected = AnovaContrastResults([AnovaContrastResult(16120.225, 36, 10.28788115, 0.0028101042, 0.47144314),
                                             AnovaContrastResult(10150.596, 36,  6.4780811,  0.0153496883, 0.39051772),
                                             AnovaContrastResult( 9653.449, 36,  6.1608033,  0.017862323,   0.38226473)])

            @test all(repeatedcontrasts(result) ≈ expected)
        end
    end
end


text:
using Distributions, Random, Plots; pyplot()
Random.seed!(4)

function simulateMM1DoobGillespie(lambda,mu,Q0,T)
    t, Q = 0.0 , Q0
    tValues, qValues = [0.0], [Q0]
    while t<T
        if Q == 0
            t += rand(Exponential(1/lambda))
            Q = 1
        else
            t += rand(Exponential(1/(lambda+mu)))
            Q += 2(rand() < lambda/(lambda+mu)) -1
        end
        push!(tValues,t)
        push!(qValues,Q)
    end
    return [tValues, qValues]
end

function stichSteps(epochs,q)
    n = length(epochs)
    newEpochs  = [ epochs[1] ]
    newQ = [ q[1] ]
    for i in 2:n
        push!(newEpochs,epochs[i])
        push!(newQ,q[i-1])
        push!(newEpochs,epochs[i])
        push!(newQ,q[i])
    end
    return [newEpochs, newQ]
end

lambda, mu = 0.7, 1.0
Tplot, Testimation = 200, 10^7
Q0 = 20

eL,qL = simulateMM1DoobGillespie(lambda, mu ,Q0, Testimation)
meanQueueLength = (eL[2:end]-eL[1:end-1])'*qL[1:end-1]/last(eL)
rho = lambda/mu
println("Estimated mean queue length: ", meanQueueLength )
println("Theoretical mean queue length: ", rho/(1-rho) )

epochs, qValues = simulateMM1DoobGillespie(lambda, mu, Q0,Tplot)
epochsForPlot, qForPlot = stichSteps(epochs,qValues)
plot(epochsForPlot,qForPlot,
	c=:blue, xlims=(0,Tplot), ylims=(0,25), xlabel="Time",
	ylabel="Customers in queue", legend=:none)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
"""
    `opensock(ipaddr, port, timeout)`
    `opensock(fun, ipaddr, port, timeout)`
    `opensock(fun, dev, timeout)`

Opens a TCP socket to a DTC Initium at ip address `ipaddr` and port `port`.
The user can also specify the timeout.

There are versions that accept a function as the first argument and opens calls the
function and closes the port. It is not very useful since when the DTC Initium closes
the port it loses current configurations.

"""
function opensock(ipaddr::IPv4, port=8400, timeout=5)

    sock = TCPSocket()
    t = Timer(_ -> close(sock), timeout)
    try
        connect(sock, ipaddr, port)
    catch e
        error("Could not connect to $ipaddr ! Turn on the device or set the right IP address!")
    finally
        close(t)
    end

    return sock
end

function opensock(fun::Function, ipaddr::IPv4, port=8400, timeout=5)

    io = opensock(ipaddr, port, timeout)
    try
        fun(io)
    finally
        close(io)
    end

end

opensock(fun::Function, dev::Initium, timeout=5) =
    opensock(fun, ipaddr(dev), portnum(dev), timeout)

"""
`SD1(dev, scnlst)`

Implements the `SD1` command which configures the scanners connected to the DTC Initium.
This function actually calls [`SD1cmd`](@ref) to build the command message.

The `scnlst` is vector of `Tuple{Int, Int, Int}` where the first element of the tuple
corresponds to the number of the connection where the scanner is located, the second
element is the number of pressure ports available and the third member of the tuple
is the logical range number (see the manual for further information).

```
"""
function SD1(dev::Initium, scnlst)

    io = socket(dev)
    isopen(io) || throw(ArgumentError("Socket not open!"))

    scn = Tuple{Int,Int,Int}[]
    #scn = scanners(dev)

    #if length(scn) > 0
    #    throw(ArgumentError("Scanners have already been added!"))
    #end

    for s in scnlst
        push!(scn, s)
    end

    cmd = SD1cmd(scn, crs=getcrs(dev))
    println(io, cmd)
    resp = read(io, 8)

    ispackerr(resp)  && throw(DTCInitiumError(resperr(resp)))

    return respconf(resp)
end

"""
`SD2(dev;  stbl=1, nfr=64, nms=1, msd=100, trm=0, scm=1, ocf=2)`

This function will configures the data acquisition parameters. The parameters are described below. For further information, see the user manual.

Parameters
 * `stbl` Chooses a DA Setup table. Should be an integer 1-5
 * `nfr` Integer that specifies the number of samples that should be averaged to create a measurement set.
 * `nms` Default number of pressure measurement sets that should be acquired. 0 if continuous data acquisition is desired.
 * `msd`  Delay between each pressure set. If 0, the scanners will acquire pressure as fast as possible
 * `trm` Trigger mode.
     - 0 no trigger
     - 1 Initial measurement is triggered
     - 2 Each pressure measurement is triggered
 * `scm` Defines scam modes. Just use 1. See the manual for further information
 * `ocf` Output format for the measurement. Just use 2 in general
     - 1 Raw measurements
     - 2 Receive pressure measurements in EU units
     - 3 EU press-sets and ez-set and temp-sets

"""
function SD2(dev::Initium; stbl=1, nfr=64, nms=1, msd=100, trm=0, scm=1, ocf=2)

    io = socket(dev)
    isopen(io) || throw(ArgumentError("Socket not open!"))

    params = setdaqparams(stbl=stbl, nfr=nfr, nms=nms, msd=msd, trm=trm, scm=scm, ocf=ocf)
    cmd = SD2cmd(params, crs=getcrs(dev))

    println(io, cmd)
    resp = read(io, 8)
    ispackerr(resp) && throw(DTCInitiumError(resperr(resp)))

    dev.params = params

    dev.isconfigured = true

    return respconf(resp)

end

"""
    `SD3(dev, stbl, ports)`

Configures the ports that should be used during data acquisition. Different port
configurations can be used, one for each `stbl` configured. This is useful where different measurements sets are necessary.

The ports are specified as a vector of [`PortRange`](@ref) objects. They can also be
specified as text ranges ("101-164") or `UnitRange{Int}` (101:164).

There are 8 scanner connections on the DTC Initium. The ports are specified by numbers
where de hundreds digit corresponds to the connection where scanner is located on the
DTC Initium. To refer to the first 32 ports on the second scanner, use "201-232".
If multiple scanners are specified, and all ports should be used, something
like "101-864" should work where all ports in 8 scanners are specified.

The function tries to check  if the specified ports are compatible with the scanners
specified in the [`SD1`](@ref) command ([`checkportlist`](@ref) ).

To get a vector with all port numbers configured, use function [`defscanlist`](@ref).

For further information, see the User Manual.

"""
function SD3(dev::Initium, stbl, ports::AbstractVector{PortRange})

    if !checkportlist(scanners(dev), ports)
        throw(ArgumentError("Invalid pressure ports"))
    end

    io = socket(dev)
    isopen(io) || throw(ArgumentError("Socket not open!"))

    pcmd = strportlist(ports)

    cmd = SD3cmd(stbl, pcmd, crs=getcrs(dev))

    println(io, cmd)
    resp = read(io, 8)
    ispackerr(resp) && throw(DTCInitiumError(resperr(resp)))

    return respconf(resp)

end

SD3(dev::Initium, stbl, ports...) = SD3(dev, stbl, portlist(ports...))

"""
`SD5(dev, actx)`

Implements the `SD5` command. This command does a lot of things but for now, only the control option is implemented (it calls [`SD5cmd`](@ref) with `stbl`=-1).

The default mode of operation uses `actx`=1. In this case, on each measurement, both the pressure and the temperature are measured simultaneously. If high speed is necessary, the temperature can be measured only once (at the beginning) and this is done setting `actx=0`. Positive values greater than 1 of `actx` specify the intermittent rate at which
temp-sets should be acquired.

"""
function SD5(dev::Initium, actx)

    stbl = -1
    io = socket(dev)

    isopen(io) || throw(ArgumentError("Socket not open!"))

    cmd = SD5cmd(stbl, actx, crs=getcrs(dev))

    println(io, cmd)
    resp = read(io, 8)
    ispackerr(resp) && throw(DTCInitiumError(resperr(resp)))

    dev.actx = actx

    return respconf(resp)

end

"""
`PC4(unx, fct=0; lrn=1)`

Change pressure units. There are predefined units for `unx` in 1-12. If `unx` is 0 or 13, a unit conversion factor from psi is employed specified by parameter `fct`.

| unx | Unit |   factor  |
| --- | ---- | --------- |
| 0   | user |  ?        |
| 1   | psi  | 1.0       |
| 2   | inH2O| 27.673    |
| 3   | Pa   | 6894.757  |
| 4   | kG/m2| 703.0696  |
| 5   | G/cm2| 70.30696  |
| 6   | ATM  | 0.068046  |
| 7   | mmHg | 51.71493  |
| 8   | mmH2O| 703.08    |
| 9   | bar  | 0.0689475 |
| 10  | kPa  | 6.894757  |
| 11  | mBar | 68.94757  |
| 12  | PSF  | 144.0     |
| 13  | user |   ?       |
"""
function PC4(dev, unx, fct=0; lrn=1)

    cmd = PC4cmd(unx, fct, lrn=lrn)

    io = socket(dev)
    isopen(io) || throw(ArgumentError("Socket not open!"))

    println(io, cmd)
    resp = read(io, 8)
    ispackerr(resp) && throw(DTCInitiumError(resperr(resp)))

    return respconf(resp)

end

"""
    `CA2(dev)`

Perform a zero calibration, basically measure the pressure and use it as a reference.
If the DTC Initium is provided with the appropriate pressure connections, the scanners
can be zeroed while measuring.

"""
function CA2(dev; lrn=1)
    cmd = CA2cmd(lrn)

    io = socket(dev)
    isopen(io) || throw(ArgumentError("Socket not open!"))

    println(io, cmd)
    resp = read(io, 8)
    ispackerr(resp) && throw(DTCInitiumError(resperr(resp)))

    return respconf(resp)

end

"""
    `AD0(dev)`

Stop data acquisition.
"""
function AD0(dev)
    cmd = AD0cmd()

    io = socket(dev)
    isopen(io) || throw(ArgumentError("Socket not open!"))

    # This is kind of tricky.
    # There might be more frames already coming in.
    # so will will try to read every packet until we get a confirmatio/error packet

    println(io, cmd)
    resp = UInt8[]

    while true
        resp = readresponse(io)
        rtype = resptype(resp)

        if rtype==4 || rtype==128
            break
        end

    end

    ispackerr(resp) && throw(DTCInitiumError(resperr(resp)))

    return respconf(resp)
end

"""
    `AD2(dev, stbl, nms)`

Acquire data using `stbl`. If `nms` is not provided it will use
the default value specified in the [`SD2`](@ref) command. This function acquires
all responses sequentially and therefore it can use any value of `ocf` (see [`SD2`](@ref))

"""
function AD2(dev, stbl, nms=-1)

    if stbl < 1 || stbl > 5
        throw(DomainError(stbl, "STBL should be between 1 and 5!"))
    end

    params = daqparams(dev)
    if !haskey(params, stbl)
        throw(DomainError(stbl, "STBL $stbl is not yet configured!"))
    end

    if nms < 0
        cmd = AD2cmd(stbl)
    else
        cmd = AD2cmd(stbl, nms)
    end
    println(cmd)
    io = socket(dev)
    isopen(io) || throw(ArgumentError("Socket not open!"))
    println(io, cmd)

    resp = Vector{UInt8}[]
    t1 = time_ns()
    t2 = UInt64(0)
    while true
        r = readresponse(io)
        rtype = resptype(r)
        t2 = time_ns()
        if rtype==4 || rtype==128
            break
        end

        push!(resp, r)
    end
    return resp, Float64(t2-t1)/1e6
end

"""
    `readresponse(io)`

Reads a response from the DTC Initium. It is capable of reading any type of response.
It will just return the bytes acquired, no processing or parsing is done.
"""
function readresponse(io)
    # Read header
    hdr = read(io, 8)
    msglen = resplen(hdr)
    if msglen > 8
        rest = read(io, msglen-8)
        return [hdr;rest]
    else
        return  hdr
    end

end
readresponse(dev::Initium) = readresponse(socket(dev))


text:
immutable IIRFilter
	filterobjPtr::Ptr{Void}
end

function Base.display(filt::IIRFilter)
	ccall((:iirfilt_rrrf_print,"libliquid"), Void, (Ptr{Void}, ), filt.filterobjPtr);
end

function createiirfilter(b::Array{Float32}, a::Array{Float32})
	ptr = ccall((:iirfilt_rrrf_create, "libliquid"), Ptr{Void}, (Ptr{Float32}, UInt, Ptr{Float32}, UInt),
		b,length(b), a, length(a));
	return IIRFilter(ptr);
end

@enum IIRFType BUTTER=0 CHEBY1 CHEBY2 ELLIP BESSEL
@enum IIRBType LPF=0 HPF BPF BSF

function getftype(input::AbstractString)
	if(input in ["lo", "low", "lowpass", "lpf"])
		return LPF
	elseif(input in ["hi","high", "hpf","highpass"])
		return HPF
	elseif(input in ["band", "bandpass", "bpf"])
		return BPF
	elseif(input in ["bandstop", "bsf"])
		return BSF
	else
		error("invalid band-type: " , input)
	end
end

function iir_filt(fc::Float64, f0::Float64, dtype::IIRFType,	band::AbstractString = "lpf",
	order::Int = 10,   Ap::Float64=1.0, As::Float64=60.0)
    btype = getftype(band);
	N = order;
	if(btype == BPF || btype == BSF)
		N = 2*order;
	end
	# r = N % 2;
	# L = (N-r) / 2;
	B = Array(Cfloat, Int(N+1));
	A = Array(Cfloat, Int(N+1));
	if(fc<=0.0 || fc>=0.5)
		error("cut-off frequency must be normalized and in range (0,0.5)")
	end
	if(f0<=0.0 || f0>=0.5)
		error("center frequency must be normalized and in range (0,0.5)")
	end
	ccall((:liquid_iirdes, "libliquid"), Void, (UInt, UInt, UInt, UInt, Cfloat,Cfloat,Cfloat,Cfloat,Ptr{Cfloat},Ptr{Cfloat}),
		dtype, btype, 1, N, Cfloat(fc), Cfloat(f0), Cfloat(Ap), Cfloat(As), B,A);
	return B,A
end

# liquid_iirdes(_ftype, _btype, _format, _n, _fc, _f0, _Ap, _As, *_B, *_A);
function butterworth(fc::Float64, f0::Float64, 	band::AbstractString = "lpf",
	order::Int = 10,   Ap::Float64=1.0, As::Float64=60.0)
	B,A = iir_filt(fc, f0, BUTTER, band, order, Ap, As )
end

function cheby1(fc::Float64, f0::Float64, 	band::AbstractString = "lpf",
	order::Int = 10,   Ap::Float64=1.0, As::Float64=60.0)
	B,A = iir_filt(fc, f0, CHEBY1, band, order, Ap, As )
end

function cheby2(fc::Float64, f0::Float64, 	band::AbstractString = "lpf",
	order::Int = 10,   Ap::Float64=1.0, As::Float64=60.0)
	B,A = iir_filt(fc, f0, CHEBY2, band, order, Ap, As )
end

function ellip(fc::Float64, f0::Float64, 	band::AbstractString = "lpf",
	order::Int = 10,   Ap::Float64=1.0, As::Float64=60.0)
	B,A = iir_filt(fc, f0, ELLIP, band, order, Ap, As )
end

function bessel(fc::Float64, f0::Float64, 	band::AbstractString = "lpf",
	order::Int = 10,   Ap::Float64=1.0, As::Float64=60.0)
	B,A = iir_filt(fc, f0, BESSEL, band, order, Ap, As )
end


text:
module TestRules

using  Test
import PlayingCards: Card
import Cribbage.Rules: cardvalue, validhands, validplays

@testset "Test cardvalue" begin

    @test cardvalue(Card("5♥")) == 5
    @test cardvalue(Card("A♠")) == 1
    @test cardvalue(Card("10♣")) == 10
    @test cardvalue(Card("J♦")) == 10
    @test cardvalue(Card("Q♥")) == 10
    @test cardvalue(Card("K♦")) == 10

end

@testset "Test validhands" begin

    originalhand = tuple(
        Card("5♥"), Card("A♠"), Card("10♣"), Card("J♦"), Card("Q♥")
    )
    hands = validhands(originalhand)
    @test length(hands) == 5
    for hand in hands
        @test length(hand) == 4
        for card in hand
            @test card ∈ originalhand
        end
    end

    originalhand = tuple(
        Card("5♥"), Card("A♠"), Card("10♣"), Card("J♦"), Card("Q♥"), Card("K♦")
    )
    hands = validhands(originalhand)
    @test length(hands) == 15
    for hand in hands
        @test length(hand) == 4
        for card in hand
            @test card ∈ originalhand
        end
    end
end

@testset "Test validplays" begin

    prior = tuple()
    hand = tuple(
        Card("5♥"), Card("10♣"), Card("J♦"), Card("Q♥")
    )
    possibleplays = validplays(prior, hand)
    @test length(possibleplays) == 4
    for card in possibleplays
        @test card ∈ hand
    end

    prior = tuple(Card("5♥"), Card("K♦"))
    hand = tuple(
        Card("10♣"), Card("J♦"), Card("Q♥")
    )
    possibleplays = validplays(prior, hand)
    @test length(possibleplays) == 3
    for card in possibleplays
        @test card ∈ hand
    end

    prior = tuple(Card("5♥"), Card("K♦"), Card("10♣"), Card("4♣"))
    hand = tuple(
        Card("J♦"), Card("Q♥")
    )
    possibleplays = validplays(prior, hand)
    @test length(possibleplays) == 0

    prior = tuple(Card("5♥"), Card("K♦"), Card("10♣"), Card("4♣"))
    hand = tuple()
    possibleplays = validplays(prior, hand)
    @test length(possibleplays) == 0
end

end  # module TestRules


text:
using DoubleFloats
using LinearAlgebra

using SemiLagrangian:
    nosplit,
    Advection,
    sizeall,
    AdvectionData,
    getdata,
    advection!,
    UniformMesh,
    AbstractInterpolation,
    Lagrange,
    B_SplineLU,
    B_SplineFFT,
    gettabmod,
    CachePrecal,
    interpolate!,
    compute_charge!,
    compute_elfield!,
    compute_ee,
    compute_ke,
    getenergy,
    getenergyall,
    dotprod,
    PoissonVar,
    getpoissonvar,
    TypePoisson,
    StdPoisson,
    StdPoisson2d,
    stdtomesh,
    TimeAlgorithm,
    NoTimeAlg,
    ABTimeAlg_ip,
    ABTimeAlg_new,
    ABTimeAlg_init,
    hamsplit_3_11,
    TimeOptimization,
    NoTimeOpt,
    SimpleThreadsOpt,
    SplitThreadsOpt,
    standardsplit,
    strangsplit,
    triplejumpsplit,
    order6split,
    hamsplit_3_11

function test_poisson2d(
    sz::NTuple{2,Int},
    interp::Vector{I},
    t_max::T,
    nbdt::Int,
) where {T,I<:AbstractInterpolation{T}}
    spmin, spmax, nsp = T(0), 4T(pi), sz[1]
    vmin, vmax, nv = T(-6), T(6), sz[2]

    mesh_sp = UniformMesh(spmin, spmax, nsp)
    mesh_v = UniformMesh(vmin, vmax, nv)

    dt = t_max / nbdt

    # dec1 = -2tan(dt/2) / step(mesh_sp) * mesh_v.points
    # dec2 = 2tan(dt/2) / step(mesh_v) * mesh_sp.points
    dec1 = zeros(T, sz)
    dec2 = zeros(T, sz)
    coef = 1 / sqrt(2T(pi))
    tabref = zeros(T, sz)
    for i = 1:sz[1], j = 1:sz[2]
        x = mesh_sp.points[i]
        y = mesh_v.points[j]
        dec1[i, j] = -dt / step(mesh_sp) * y
        tabref[i, j] = coef * exp(-0.5 * y^2) * (1 + 0.001 * cos(x / 2))
    end

    #    @show dec1, dec2

    cache = CachePrecal(interp, t_max)

    rho = zeros(T, sz[1])
    elf = zeros(T, sz[1])

    data = zeros(T, sz)

    copyto!(data, tabref)

    compute_charge!(rho, (mesh_sp,), data)
    compute_elfield!(elf, mesh_sp, rho)
    ee = compute_ee((mesh_sp,), (elf,))
    ke = compute_ke((mesh_sp,), (mesh_v,), data)

    println("0\t$ee\t$ke")

    buf = zeros(T, sz)
    diffmax = 0
    tabmod = gettabmod.(sz)
    for ind = 1:2nbdt
        dec2 = [dt / step(mesh_v) * elf[i] for i = 1:sz[1], j = 1:sz[2]]
        interpolate!(
            buf,
            data,
            ind -> (dec1[ind], dec2[ind]),
            interp,
            tabmod = tabmod,
            cache = cache,
        )
        copyto!(data, buf)

        compute_charge!(rho, (mesh_sp,), data)
        compute_elfield!(elf, mesh_sp, rho)

        ee = compute_ee((mesh_sp,), (elf,))
        ke = compute_ke((mesh_sp,), (mesh_v,), data)

        println("$(Float32(ind*dt))\t$ee\t$ke")
    end
    return diffmax
end
verif(pv::PoissonVar{T}, advd::AdvectionData) where {T} = missing
# function verif(pv::PoissonVar{T,N,Nsp,Nv,StdPoisson2dTry}, advd::AdvectionData) where{T,N,Nsp,Nv}
#     adv = advd.adv
#     sz = sizeall(adv)
#     d = zeros(T, sz)
#     coef = 1 / sqrt(2T(pi))
#     for ind in CartesianIndices(sz)
#         sp = stdtomesh(adv.t_mesh[1],pv.t_sp[ind])
#         v = stdtomesh(adv.t_mesh[2],pv.t_v[ind])
#         d[ind] = coef * exp(T(-0.5) * v^2) * (1 + T(big"0.001") * cos( sp / 2))
#     end
#     # println("t_sp=$(pv.t_sp[1:10,1:10])")
#     # println("t_v=$(pv.t_v[1:10,1:10])")
#     if !ismissing(pv.bufcur_sp)
#         # println("bufc_sp=$(pv.bufcur_sp[1:10])")
#         # println("bufc_v=$(pv.bufcur_v[1:10])")
#         println("extrema sp : $(extrema(pv.bufcur_sp))")
#         println("extrema v : $(extrema(pv.bufcur_v))")
#     end
#     res = norm(d-advd.data)
#     res2 = norm(d-advd.data, Inf)
#     println("verif : res=$res res2=$res2")
#     @show T
# end
function initdata(adv::Advection{T}) where {T}
    sz = sizeall(adv)
    coef = 1 / sqrt(2T(pi))
    tabref = zeros(T, sz)
    for i = 1:sz[1], j = 1:sz[2]
        x = adv.t_mesh[1].points[i]
        y = adv.t_mesh[2].points[j]
        tabref[i, j] = coef * exp(T(-0.5) * y^2) * (1 + T(big"0.5") * cos(x / 2))
    end
    return tabref
end

function get_init(advdref::AdvectionData{T}, nb::Int) where {T}

    dt = advdref.adv.dt_base

    tabst = [([2, 1], 1, 1, true), ([1, 2], 1, 2, true)]

    lag = Lagrange(19, T)

    adv = Advection(advdref.adv.t_mesh, [lag, lag], dt, tabst; tab_coef = hamsplit_3_11(dt))

    pvar = getpoissonvar(adv)

    advd = AdvectionData(adv, advdref.data, pvar)

    t_res = []
    for i = 1:nb
        while advection!(advd)
        end
        push!(t_res, copy(advd.data))
    end

    return t_res
end

function test_poisson2dadv(
    sz::NTuple{2,Int},
    interp::Vector{I},
    t_max::T,
    nbdt::Int,
    type::TypePoisson = StdPoisson2d,
    typeadd = 0,
    timealg::TimeAlgorithm = NoTimeAlg,
    ordalg = 0;
    timeopt::TimeOptimization=NoTimeOpt,
    flbiginit::Bool=false,
) where {T,I<:AbstractInterpolation{T}}
    spmin, spmax, nsp = T(0), 4T(pi), sz[1]
    vmin, vmax, nv = T(-9), T(9), sz[2]

    mesh_sp = UniformMesh(spmin, spmax, nsp)
    mesh_v = UniformMesh(vmin, vmax, nv)

    dt = t_max / nbdt

    tabst = [([1, 2], 2, 1, false)]

    adv = Advection(
        (mesh_sp, mesh_v),
        interp,
        dt,
        tabst,
        tab_coef = nosplit(dt),
        timealg = timealg,
        ordalg = ordalg,
        timeopt=timeopt,
    )

    tabref = initdata(adv)

    data = zeros(T, sz)
    copyto!(data, tabref)

    pvar = getpoissonvar(adv, type = type, typeadd = typeadd)

    # initdatas = timealg == ABTimeAlg_init ? map(x -> tabref, 1:(ordalg-1)) : missing

    # advd = AdvectionData(adv, data, pvar; initdatas = initdatas)

    advd = AdvectionData(adv, data, pvar)

    if timealg == ABTimeAlg_init
        nb = flbiginit ? 3ordalg-1 : ordalg-1
        advd.time_cur -= dt*nb
        advd.initdatas = get_init(advd, nb)
    end

    elenergy, kinenergy, energyall = getenergy(advd)

    verif(pvar, advd)

    enmax = enmin = energyall
    @show enmax, enmin

    diff = 0
    diffprec = 0

    borne_t = t_max - dt/2

    while advd.time_cur < borne_t
        while advection!(advd)
        end
        elenergy, kinenergy, energyall = getenergy(advd)
        enmax = max(energyall, enmax)
        enmin = min(energyall, enmin)

        t = advd.time_cur
        diff = enmax - enmin
        delta = diff - diffprec
        diffprec = diff
        println("t=$(Float32(t)) diff=$diff delta=$delta")
        println(
            "$(Float32(t))\t$(Float64(elenergy))\t$(Float64(kinenergy))\t$(Float64(energyall))",
        )
        verif(pvar, advd)
    end
    @show enmax, enmin
    return (enmax - enmin), advd.data
end

function test_poisson2d2d_adv(
    sz::NTuple{4,Int},
    interp::Vector{I},
    t_max::T,
    nbdt::Int,
    split=standardgsplit
) where {T,I<:AbstractInterpolation{T}}

    mesh_sp1 = UniformMesh(T(0), 4T(pi), sz[1])
    mesh_sp2 = UniformMesh(T(0), 4T(pi), sz[2])
    mesh_v1 = UniformMesh(T(-10), T(10), sz[3])
    mesh_v2 = UniformMesh(T(-10), T(10), sz[4])

    dt = t_max / nbdt

    tabst = [([1,2,3,4], 2, 1, true,),
    ([3,4,1,2,], 2, 2, true,),]
#    tabst = [([3, 4, 1, 2], 2, 1, true,),
#    ([1, 2, 3, 4], 2, 2, true,),]
#    tabst = [([3, 4, 1, 2], 1, 1, true,),
#    ([4, 3, 1, 2], 1, 1, true,),
#    ([1, 2, 3, 4], 1, 2, true,),
#    ([2,1, 3, 4], 1, 2, true,),]

    adv = Advection((mesh_sp1, mesh_sp2, mesh_v1, mesh_v2), interp, dt, tabst, tab_coef= split(dt))
    println("trace0")

    epsilon = T(0.5)
    fct_sp(x) = epsilon * cos(x / 2) + 1
    fct_v(v) = exp(-v^2 / 2) / sqrt(2T(pi))
    lgn1_sp = fct_sp.(mesh_sp1.points)
    lgn1_v = fct_v.(mesh_v1.points)
    lgn2_sp = fct_sp.(mesh_sp2.points)
    lgn2_v = fct_v.(mesh_v2.points)
    println("trace1")
    data = dotprod((lgn1_sp, lgn2_sp, lgn1_v, lgn2_v))
    println("trace2")

    pvar = getpoissonvar(adv)

    advd = AdvectionData(adv, data, pvar)

    elenergy, kinenergy, energyall = getenergy(advd)

    enmax = enmin = energyall

    diff = 0
    diffprec = 0
    for i = 1:nbdt
        while advection!(advd)
        end
        elenergy, kinenergy, energyall = getenergy(advd)
        enmax = max(energyall, enmax)
        enmin = min(energyall, enmin)
        t = i * dt
        diff = enmax - enmin
        delta = diff - diffprec
        diffprec = diff
        println("t=$(Float32(t)) diff=$diff delta=$delta")
        println(
            "$(Float32(t))\t$(Float64(elenergy))\t$(Float64(kinenergy))\t$(Float64(energyall))",
        )
    end
    @show enmax, enmin
    return enmax - enmin
end
function test_split4d(T, nbdt, split, order)
    sz = (32,32,32,32)
    lag = Lagrange(9,T)
    t_max = T(0.1)
    interp = [lag,lag,lag,lag]

    ret1 = test_poisson2d2d_adv(sz, interp, t_max, nbdt, split)
    ret2 = test_poisson2d2d_adv(sz, interp, t_max, 2nbdt, split)
    @show split, ret1, ret2, ret1/ret2

    @test 1.2ret1/ret2 > 2^order

end
function test_timealg(interp::Vector{I},nbdt, timealg,ordalg; flbiginit=false) where {T,I<:AbstractInterpolation{T}}
    sz = (128,100)
   t_max = T(big"0.1")
    ret1, _ = test_poisson2dadv(
        sz,
        interp,
        t_max,
        nbdt,
        StdPoisson2d,
        0,
        timealg,ordalg, flbiginit=flbiginit)
    nbdt *= 2
    ret2, _ = test_poisson2dadv(
        sz,
        interp,
        t_max,
        nbdt,
        StdPoisson2d,
        0,
        timealg,ordalg, flbiginit=flbiginit)

    @show ret1,ret2, ret1/ret2

    @test 1.25*ret1/ret2 > 2 ^ ordalg
end

@testset "test poisson2d ABTimeAlg_new" begin

    T = Double64
    interp = [B_SplineLU(11, 128, T), B_SplineLU(11, 100, T)]
    @time test_timealg(interp, 5, ABTimeAlg_new, 2 )

end

@testset "test poisson2d ABTimeAlg_ip" begin
    T = Double64
    interp = [Lagrange(7,T),Lagrange(7,T)]
    @time test_timealg(interp,5,ABTimeAlg_ip,2)
    @time test_timealg(interp,5,ABTimeAlg_ip,3)
    @time test_timealg(interp,5,ABTimeAlg_ip,4)
end

@testset "test poisson2d ABTimeAlg_init" begin
    T = Double64
    interp = [Lagrange(7,T),Lagrange(7,T)]
    @time test_timealg(interp,5,ABTimeAlg_init,2)
    @time test_timealg(interp,5,ABTimeAlg_init,3)
    @time test_timealg(interp,5,ABTimeAlg_init,3)
    @time test_timealg(interp,5,ABTimeAlg_init,4)
end
@testset "test poisson2d ABTimeAlg_init 2" begin
    T = Double64
    interp = [Lagrange(9,T),Lagrange(9,T)]
    @time test_timealg(interp,5,ABTimeAlg_init,2, flbiginit=true)
    @time test_timealg(interp,5,ABTimeAlg_init,3, flbiginit=true)
    @time test_timealg(interp,5,ABTimeAlg_init,4, flbiginit=true)
    @time test_timealg(interp,10,ABTimeAlg_init,5, flbiginit=true)
end

@testset "test poisson2d thread" begin
    T = Float64
    interp = [Lagrange(5,T),Lagrange(5,T)]
    sz = (64,64)
    flag=false
    for t_alg in (NoTimeAlg, ABTimeAlg_ip, ABTimeAlg_init, ABTimeAlg_init)
        res = []
        for t_opt in (NoTimeOpt, SimpleThreadsOpt, SplitThreadsOpt)
            @show t_alg, t_opt
            ordalg = t_alg == NoTimeAlg ? 0 : 2
            _, d1 = test_poisson2dadv(
                sz,
                interp,
                T(0.01),
                5,
                StdPoisson2d,
                0,
                t_alg, ordalg,
                timeopt=t_opt,
                flbiginit=flag)
            push!(res, d1)
            if length(res) >= 2
                @test res[end-1] == res[end]
            end
        end
        flag = t_alg == ABTimeAlg_init
    end
end

@testset "test split 4d" begin
    T = Float64
    @time test_split4d(T, 5, standardsplit, 1)
    # @time test_split4d(T, 5, strangsplit, 2)
    # @time test_split4d(T, 5, triplejumpsplit, 4)
    # @time test_split4d(T, 5, order6split, 6)
    # @time test_split4d(T, 5, hamsplit_3_11, 6)
end


text:
using ACEds.equivRotations3D: rmatrices

print("==========================\n")
for m in -1:1:1
    for mu in -1:1:1
        print("indices: ", m, " ", mu, "\n")
        print("\n")
        display(rmatrices[(m,mu)])
        print("-------------------------------\n")
    end
end


text:
# # Sampling from a GP
#
# ## Preliminary steps

# ### Loading necessary packages

using Plots
using AugmentedGaussianProcesses
using Distributions
using LinearAlgebra

# ### Generating some random data
kernel = SqExponentialKernel()
x = range(0, 10; length=50)
K = kernelmatrix(kernel, x)
f = rand(MvNormal(K + 1e-8I)) # Sample a random GP
y = rand.(Bernoulli.(AGP.logistic.(f)))
y_sign = Int.(sign.(y .- 0.5))

# ### We create a function to visualize the data

function plot_data(x, y; size=(300, 500))
    return Plots.scatter(x, y; alpha=0.2, markerstrokewidth=0.0, lab="", size=size)
end
plot_data(x, y; size=(500, 500))

# ### Run the variational gaussian process approximation
@info "Running full model"
mfull = VGP(x, y_sign, kernel, LogisticLikelihood(), AnalyticVI(); optimiser=false)
@time train!(mfull, 5)

# ### We can also create a sampling based model
@info "Sampling from model"
mmcmc = MCGP(x, y, kernel, LogisticLikelihood(), GibbsSampling(); optimiser=false)
m = mmcmc
@time samples = sample(mmcmc, 1000)

# ### We can now visualize the results of both models

# ### We first plot the latent function f (truth, the VI estimate, the samples)
p1 = plot(x, f; label="true f")
plot!(x, samples; label="", color=:black, alpha=0.02, lab="")
plot!(x, mean(mfull[1]); ribbon=sqrt.(var(mfull[1])), label="VI")
# ### And we can also plot the predictions vs the data
p2 = plot_data(x, y; size=(600, 400))
μ_vi, σ_vi = proba_y(mfull, x)
plot!(x, μ_vi; ribbon=σ_vi, label="VI")
μ_mcmc, σ_mcmc = proba_y(mmcmc, x)
plot!(x, μ_mcmc; ribbon=σ_mcmc, label="MCMC")


text:
import Colors: Colorant, Color, RGB, red, green, blue, parse, convert
import PyCall: @pyimport, PyObject
@pyimport cycler
import PyPlot: plt

# I need a way to convert the types from Colors.jl (RGB{8}(r, g, b) -> python tuple)
# this code is directly from PyCall.jl -> conversions.jl for tuple conversion
function PyObject(t::Color)
    trgb = convert(RGB, t)
    ctup = map(float, (red(trgb), green(trgb), blue(trgb)))
    o = PyObject(ctup)
    return o
end

#TODO: add RGBA versions

# string color names
# the order will matter ... or maybe I should have it return multiple
color_defs = [svg_rgb, crayons, xkcd_rgb]

function colorname2rgb(name::String)
    for def in color_defs
        if haskey(def, name)
            return parse(Colorant, def[name])
        end
    end
    throw(ArgumentError("Not a Color name: $name"))
end

#NC for Named Color
macro NC_str(name::String)
    c = colorname2rgb(name)
    return :($c)
end

#HC for Hex Color
macro HC_str(hexstr::String)
    #TODO: test to ensure actually a hex string
    c = parse(Colorant, hexstr)
    return :($c)
end

# Named color palettes
SEABORN_PALETTES = Dict(
:deep => [HC"#4C72B0", HC"#55A868", HC"#C44E52", HC"#8172B2", HC"#CCB974", HC"#64B5CD"],
:muted => [HC"#4878CF", HC"#6ACC65", HC"#D65F5F", HC"#B47CC7", HC"#C4AD66", HC"#77BEDB"],
:pastel => [HC"#92C6FF", HC"#97F0AA", HC"#FF9F9A", HC"#D0BBFF", HC"#FFFEA3", HC"#B0E0E6"],
:bright => [HC"#003FFF", HC"#03ED3A", "#E8000B", HC"#8A2BE2", HC"#FFC400", HC"#00D7FF"],
:dark => [HC"#001C7F", HC"#017517", "#8C0900", HC"#7600A1", HC"#B8860B", HC"#006374"],
:colorblind => [HC"#0072B2", HC"#009E73", HC"#D55E00", HC"#CC79A7", HC"#F0E442", HC"#56B4E9"])

color_cycle(carr::AbstractArray) = cycler.cycler("color", carr)

#TODO: this needs a more a more julian name something like style!, palette!
function set_palette(name::Symbol)
    colprop = color_cycle(SEABORN_PALETTES[name])
    plt[:rc]("axes", prop_cycle=colprop)
end

max_stars_repo_name: gabrielgellner/OdinSon

max_issues_repo_name: gabrielgellner/OdinSon

max_forks_repo_name: gabrielgellner/OdinSon


text:
export sample!
using StatsBase: sample, Weights

"""
    sample!(model::ABM, n [, weight]; kwargs...)

Replace the agents of the `model` with a random sample of the current agents with
size `n`.

Optionally, choose an agent property `weight` (Symbol) to weight the sampling.
This means that the higher the `weight` of the agent, the higher the probability that
this agent will be chosen in the new sampling.

# Keywords
* `replace = true` : whether sampling is performed with replacement, i.e. all agents can
  be chosen more than once.
* `rng = GLOBAL_RNG` : a random number generator to perform the sampling with.

See the Wright-Fisher example in the documentation for an application of `sample!`.
"""
function sample!(model::ABM, n::Int, weight=nothing; replace=true,
    rng::AbstractRNG=Random.GLOBAL_RNG)

    if weight != nothing
        weights = Weights([getproperty(a, weight) for a in values(model.agents)])
        newids = sample(rng, collect(keys(model.agents)), weights, n, replace=replace)
    else
        newids = sample(rng, collect(keys(model.agents)), n, replace=replace)
    end

    for (index, id) in enumerate(newids) # add new agents while adjusting id
        model.agents[index] = deepcopy(model.agents[id])
        model.agents[index].id = index
    end
    # kill extra agents
    if n < nagents(model)
        genocide!(model, n)
    end
    # Clean space
    clean_space!(model)

    return model
end

"""
Remove all IDs from space and add agent ids again.
"""
function clean_space!(model::ABM)
    if model.space != nothing
        for node in 1:nv(model)
            model.space.agent_positions[node] = Int[]
        end
        for (k, v) in model.agents
            push!(get_node_contents(v,model), v.id)
        end
    end
end


text:
module Instruments

export Instrument, GenericInstrument, connect!, disconnect!, write, read, query
export ResourceManager
export find_resources
export @scpifloat
export @scpibool

import Base: write, read, readavailable

# load the binary dependency path
if isfile(joinpath(dirname(dirname(@__FILE__)),"deps","deps.jl"))
    include("../deps/deps.jl")
else
    #error("Instruments.jl not properly installed. Please run Pkg.build(\"Instruments\")")
    @warn "No VISA libraries found, please check VISA installation for visa64 library!"
end

include("visa/VISA.jl")

include("instrument.jl")

include("scpi.jl")

ResourceManager() = viOpenDefaultRM()

# Helper functions to find instruments
find_resources(rm, expr::AbstractString="?*::INSTR") = Instruments.viFindRsrc(rm, expr)

end # module


text:
module LabelMe

import Base.==

export Annotation, Source, Object, Polygon, XML

type Polygon
    username::AbstractString
    points::Vector{Tuple{Float32,Float32}}

    Polygon() = new("", Vector())
    Polygon(username, points) = new(username, points)
end

==(a::Polygon, b::Polygon) = a.username == b.username && a.points == b.points

type Object
    id::Nullable{Integer}
    name::AbstractString
    deleted::Nullable{Bool}
    verified::Nullable{Bool}
    occluded::Nullable{Bool}
    attributes::AbstractString
    parts_parent::Nullable{Integer}
    parts_children::Vector{Integer}
    date::DateTime
    polygon::Polygon

    Object() = new(nothing, "", nothing, nothing, nothing, "", nothing,
        Vector(), DateTime(), Polygon())
    Object(id, name, deleted, verified, occluded, attributes, parts_parent, parts_children, date, polygon) =
        new(id, name, deleted, verified, occluded, attributes, parts_parent, parts_children, date, polygon)
    Object(polygon::Polygon) = new(nothing, "", nothing, nothing, nothing, "", nothing,
        Vector(), DateTime(), polygon)
end

function ==(a::Object, b::Object)
    return isequal(a.id, b.id) &&
        a.name == b.name &&
        isequal(a.deleted, b.deleted) &&
        isequal(a.verified, b.verified) &&
        isequal(a.occluded, b.occluded) &&
        a.attributes == b.attributes &&
        isequal(a.parts_parent, b.parts_parent) &&
        a.parts_children == b.parts_children &&
        a.date == b.date &&
        a.polygon == b.polygon
end

type Source
    sourceImage::AbstractString
    sourceAnnotation::AbstractString

    Source() = new("", "")
    Source(sourceImage, sourceAnnotation) = new(sourceImage, sourceAnnotation)
end

==(a::Source, b::Source) = a.sourceImage == b.sourceImage && a.sourceAnnotation == b.sourceAnnotation

type Annotation
    filename::AbstractString
    folder::AbstractString
    source::Source
    image_size::Nullable{Tuple{Integer,Integer}}
    objects::Vector{Object}

    Annotation() = new("", "", Source(), nothing, Vector())
    Annotation(filename, folder, source, image_size, objects) = new(filename, folder, source, image_size, objects)
end

function ==(a::Annotation, b::Annotation)
    return a.filename == b.filename &&
        a.folder == b.folder &&
        isequal(a.source, b.source) &&
        isequal(a.image_size, b.image_size) &&
        a.objects == b.objects
end

include("XML.jl")

end # module


text:
"""
t1.jl - do some calculations for fluxonium t1
"""

using Dierckx
using HDF5
using Interpolations
using LaTeXStrings
using LinearAlgebra
using TrajectoryOptimization
import Plots
using Printf
using StaticArrays

# Construct paths.
EXPERIMENT_META = "spin"
EXPERIMENT_NAME = "spin15"
WDIR = ENV["ROBUST_QOC_PATH"]
SAVE_PATH = joinpath(WDIR, "out", EXPERIMENT_META, EXPERIMENT_NAME)
T1_SAVE_FILE_PATH = joinpath(SAVE_PATH, "t1.h5")
T1CMP_PLOT_FILE_PATH = joinpath(SAVE_PATH, "t1splinecmp.png")

# Plotting configuration.
ENV["GKSwstype"] = "nul"
Plots.gr()
DPI = 300

# misc. constants
SAMPLE_SIZE = Int(1e3)

# Define experimental constants.
# E / h
EC = 0.479e9
EL = 0.132e9
EJ = 3.395e9
# Q_CAP = 1 / 8e-6
Q_CAP = 1.25e5
T_CAP = 0.042
H = 6.62607015e-34
HBAR = 1.05457148e-34
KB = 1.3806503e-23
# HBAR_BY_KB = HBAR / KB
HBAR_BY_KB = 7.638225390210369e-12
E_CHARGE = 1.60217646e-19
# FLUX_QUANTUM = H / (2 * E_CHARGE)
# INVERSE_FLUX_QUANTUM = 2 * E_CHARGE / H
INVERSE_FLUX_QUANTUM = 4.835977958971655e14
# FLUX_FRUSTRATION = FLUX_QUANTUM / 2
FLUX_FRUSTRATION = 1.033917036516689e-15
# raw T1s reported in units of microseconds
T1_ARRAY = [
    1597.923, 1627.93, 301.86, 269.03, 476.33, 1783.19, 2131.76, 2634.50,
    4364.68, 2587.82, 1661.915, 1794.468, 2173.88, 1188.83,
    1576.493, 965.183, 560.251, 310.88
] * 1e3
FBFQ_ARRAY = [
    0.26, 0.28, 0.32, 0.34, 0.36, 0.38, 0.4,
    0.42, 0.44, 0.46, 0.465, 0.47, 0.475,
    0.48, 0.484, 0.488, 0.492, 0.5
]

# Define the system.
# FLUXONIUM_STATE_COUNT is the state count used in the T1 calculations
FLUXONIUM_STATE_COUNT = 110
FLUXONIUM_LEVELS = Array(range(0., stop=FLUXONIUM_STATE_COUNT - 1,
                               length=FLUXONIUM_STATE_COUNT))
SQRT_FLUXONIUM_LEVELS_TRUNC = map(sqrt, FLUXONIUM_LEVELS[2:FLUXONIUM_STATE_COUNT])
ANNIHILATE = diagm(1 => SQRT_FLUXONIUM_LEVELS_TRUNC)
CREATE = diagm(-1 => SQRT_FLUXONIUM_LEVELS_TRUNC)
GAMMA_CAP_PREFACTOR = (16 * pi * EC * Q_CAP)^(-1)
E_PLASMA = (8 * EL * EC)^(0.5)
PHI_OSC = (8 * EC / EL)^(0.25)
PHI_OP = PHI_OSC * 2^(-0.5) * (CREATE + ANNIHILATE)
H_EXP_RAW = exp(1im * PHI_OP)
H_LC = diagm(E_PLASMA * FLUXONIUM_LEVELS)

### ANALYTICAL CALCULATIONS ###
function get_hamiltonian_oscillator_basis(flux :: Float64)
    reduced_flux = 2 * pi * flux * INVERSE_FLUX_QUANTUM
    h_exp = H_EXP_RAW * exp(1im * reduced_flux)
    h_cos = 0.5 * (h_exp + h_exp')
    h_fluxonium = real(H_LC - EJ * h_cos)
    return h_fluxonium
end

function get_t1(flux :: Float64)
    # cap
    h_fluxonium = get_hamiltonian_oscillator_basis(flux)
    evals, evecs = eigen(Hermitian(h_fluxonium))
    g0_phi_g1 = evecs[:, 1]' * PHI_OP * evecs[:, 2]
    omega = (evals[2] - evals[1]) * (2 * pi)
    gamma_cap = (GAMMA_CAP_PREFACTOR * omega^2 * coth(HBAR_BY_KB * omega / (2 * T_CAP))
                 * g0_phi_g1 * conj(g0_phi_g1))

    gamma = gamma_cap
    t1 = gamma^(-1)

    return t1
end

### EXPERIMENTAL DATA FIT ###
function compare_splines()
    # fit
    t1_spline_dierckx = Spline1D(FBFQ_ARRAY, T1_ARRAY)
    t1_spline_itp = extrapolate(interpolate((FBFQ_ARRAY,), T1_ARRAY, Gridded(Linear())), Flat())
    # fbfq_axis = range(minimum(FBFQ_ARRAY), stop=maximum(FBFQ_ARRAY), length=SAMPLE_SIZE)
    fbfq_axis = range(0, stop=2 * maximum(FBFQ_ARRAY), length=SAMPLE_SIZE)

    # plot
    fig = Plots.plot(dpi=DPI)
    Plots.plot!(fbfq_axis, map(t1_spline_dierckx, fbfq_axis), label="dierckx")
    Plots.plot!(fbfq_axis, map(t1_spline_itp, fbfq_axis), label="itp")
    Plots.scatter!(FBFQ_ARRAY, T1_ARRAY, label="data")
    Plots.xlabel!(L"\Phi / \Phi_{0}")
    Plots.ylabel!(L"T_{1}")
    Plots.savefig(fig, T1CMP_PLOT_FILE_PATH)
end

max_stars_repo_name: SchusterLab/rbqoc

max_issues_repo_name: SchusterLab/rbqoc

max_forks_repo_name: SchusterLab/rbqoc


text:
"""
    AmoebaFiber3D(f::AbstractPolynomial, w=(0., 0., 0.))

Construct the fiber `AmoebaFiber3D` of the bivariate polynomial `f` at `w`.
"""
function AmoebaFiber3D(p::MP.AbstractPolynomialLike{S}, w=(0.0, 0.0, 0.0)) where S
    vars = MP.variables(p)
    @assert (length(vars) == 3) "Expected a trivariate polynomial, but got $(p)."

    T = realified_coefficient_type(S)
    f_re, f_im = SP.Polynomial.(1.0 .* realify(p))

    v = exp.(w)
    U = zeros(T, 2, 6)
    x = SVector{6, Float64}(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

    exponents = SP.exponents.((f_re, f_im))
    original_coefficients = copy.(SP.coefficients.((f_re, f_im)))

    AmoebaFiber3D(f_re, f_im, v, U, x, exponents, original_coefficients)
end

function update_fiber!(F::AmoebaFiber3D, w)
    w1, w2, w3 = w
    F.v = (exp(w1), exp(w2), exp(w3))
    condition!(F)
    F
end

function condition!(fiber::AmoebaFiber3D)
    fiber.f_re.coefficients .= fiber.original_coefficients[1]
    fiber.f_im.coefficients .= fiber.original_coefficients[2]

    n1 = normalizer_amoeba_3d(fiber.exponents[1], fiber.v)
    n2 = normalizer_amoeba_3d(fiber.exponents[2], fiber.v)
    normalizer = inv(max(n1, n2))
    SP.scale_coefficients!(fiber.f_re, normalizer)
    SP.scale_coefficients!(fiber.f_im, normalizer)
    nothing
end

function normalizer_amoeba_3d(E, v)
    v1, v2, v3 = v
    normalizer = sum(1:size(E, 2)) do j
        m1 = E[1,j]+E[4,j]
        m2 = E[2,j]+E[5,j]
        m3 = E[3,j]+E[6,j]
        Base.FastMath.pow_fast(v1, m1) *
        Base.FastMath.pow_fast(v2, m2) *
        Base.FastMath.pow_fast(v3, m3)
    end
end

function condition_amoeba_fiber_3d!(f, E, v)
    v1, v2, v3 = v
    normalizer = sum(1:size(E, 2)) do j
        m1 = E[1,j]+E[4,j]
        m2 = E[2,j]+E[5,j]
        m3 = E[3,j]+E[6,j]
        Base.FastMath.pow_fast(v1, m1) *
        Base.FastMath.pow_fast(v2, m2) *
        Base.FastMath.pow_fast(v3, m3)
    end
    SP.scale_coefficients!(f, 1 / normalizer)
end

function phi!(fiber::AmoebaFiber3D, θ)
    v1, v2, v3 = fiber.v
    x = SVector(
        v1 * cos(θ[1]), v2 * cos(θ[2]), v3 * cos(θ[3]),
        v1 * sin(θ[1]), v2 * sin(θ[2]), v3 * sin(θ[3])
        )
    fiber.x = x
    x
end

function evaluate(fiber::AmoebaFiber3D, θ)
    x = phi!(fiber, θ)
    SVector(SP.evaluate(fiber.f_re, x), SP.evaluate(fiber.f_im, x))
end

function jacobian(fiber::AmoebaFiber3D, θ, precomputed=false)
    x = precomputed ? fiber.x : phi!(fiber, θ)
    U = fiber.U
    @inbounds begin
        ∇f_re = SP.gradient(fiber.f_re, x)
        for i=1:6
            U[1, i] = ∇f_re[i]
        end
        U[2, 4] = U[1, 1]
        U[2, 5] = U[1, 2]
        U[2, 6] = U[1, 3]
        U[2, 1] = -U[1, 4]
        U[2, 2] = -U[1, 5]
        U[2, 3] = -U[1, 6]
    end

    @inbounds out = begin
        x4, x5, x6 = -x[4], -x[5], -x[6]
        ∂re∂θ1 = U[1, 1] * x4
        ∂re∂θ2 = U[1, 2] * x5
        ∂re∂θ3 = U[1, 3] * x6
        ∂im∂θ1 = U[2, 1] * x4
        ∂im∂θ2 = U[2, 2] * x5
        ∂im∂θ3 = U[2, 3] * x6

        ∂re∂θ1 = muladd(U[1, 4], x[1], ∂re∂θ1)
        ∂re∂θ2 = muladd(U[1, 5], x[2], ∂re∂θ2)
        ∂re∂θ3 = muladd(U[1, 6], x[3], ∂re∂θ3)
        ∂im∂θ1 = muladd(U[2, 4], x[1], ∂im∂θ1)
        ∂im∂θ2 = muladd(U[2, 5], x[2], ∂im∂θ2)
        ∂im∂θ3 = muladd(U[2, 6], x[3], ∂im∂θ3)

        # SMatrix has column-major storage
        SMatrix{2,3}(∂re∂θ1, ∂im∂θ1, ∂re∂θ2, ∂im∂θ2, ∂re∂θ3, ∂im∂θ3)
    end

    out
end


text:
###
#   Utilities and show
###

left(x :: Interval) = x.lo
right(x :: Interval) = x.hi

left(x :: IntervalU) = left.(x.v)
right(x :: IntervalU) = right.(x.v)

isequal(x :: IntervalU, y ::IntervalU) = all(isequal.(sort(x.v),sort(y.v)))

function in( x:: IntervalU, y :: Vector{IntervalU})
    for i =1:length(y)
        if x == y[i]
            return true
        end
    end
    return false
end

function Base.show(io::IO, this::IntervalU)
    v = sort(this.v)
    print(io, join(v, " ∪ "));
end


text:
include("opcodes.jl")
include("macro.jl")
include("show.jl")

#=
The general idea here is that a macro will be created that can take
pseudo-assembly much like the C++ version of the simulator and turn it into
a vector of type AsapInstruction. Core models will then have this vector of
instructions as their program and keep a local PC to point to which instruction
they are executing on a given clock cycle.

Lables will be encoded using the standard Julia @label macro. During source code
tansformation, labels will be converted to indices of instructions in the
final vector of instructions.

- Alternatives:

This could be encoded using something like a ResumableFunction from the
ResumableFunctions package, but I feel like it would be difficult to bind the
functions to their respecitve cores in a clear and mostly type-stable manner.

Plus, using the vector of instructions format would allow an assembler to take
assembly files generated by the compiler and execute them without needing
to recompile any code.
=#

################################################################################
# ASSEMBLER
################################################################################
struct AsapProgram
    instructions :: Vector{AsapInstruction}
    labels       :: Dict{Symbol, Int}
end

AsapProgram() = AsapProgram(AsapInstruction[], Dict{Symbol,Int}())

Base.getindex(p::AsapProgram, i) = p.instructions[i]
Base.length(p::AsapProgram) = length(p.instructions)

function Base.show(io::IO, p::AsapProgram, pc = 0)
    # Reverse the label dict to go from istruction numbers to symbols.
    label_rev = Dict(v => k for (k,v) in p.labels)
    for i in 1:length(p.instructions)
        # Check if this instruction is the start of a label.
        if haskey(label_rev, i)
            println(io, "\n", label_rev[i], ":")
        end

        # Check if PC is here. If so, mark it with an arrow.
        leader = i == pc ? "--> " : "    "
        # Print out instruction.
        println(io, leader, "$i : $(p.instructions[i])")
    end
end

function assemble(temp_program :: Vector{InstructionLabelTarget})
    # First - need to take care of END_RPT() instructions.
    handle_rpt!(temp_program)

    # With repeats handled, we aren't removing any more instructions, so we can
    # split apart the labels and the instructions.
    program = [i.instruction for i in temp_program]
    labels = [i.label for i in temp_program]
    branch_targets = [i.branch_target for i in temp_program]

    # Make a dictionary of the first index where each label is seen.
    labeldict = makelabeldict(labels)

    # Set branch targets
    for (index, inst) in enumerate(program)
        if inst.op == BR || inst.op == BRL
            # Get the branch target symbol for this instruction
            branch_target = branch_targets[index]

            # If this is a return instruction, don't do anything.
            branch_target == :back && continue

            # Otherwise, get the index of this label
            target_index = labeldict[branch_target]

            program[index] = set(inst, set_branch_target(NULL, target_index))
        end
    end

    return AsapProgram(program, labeldict)
end

function handle_rpt!(program :: Vector{InstructionLabelTarget})
    # Iterate through each element of the intermediate instructions.
    index_of_last_rpt = 0
    index = 1
    while index <= length(program)
        # Get the instruction here
        instruction = program[index].instruction
        if instruction.op == END_RPT
            # Save the index of the instruction before this as the end of the
            # repeat block and delete this op from the instruction vector.
            last_rpt = program[index_of_last_rpt].instruction
            program[index_of_last_rpt].instruction = set(last_rpt, set_repeat_end(index - 1))

            # In this case, don't increment "index" because it will
            # automatically point to the next instruction since we deleted
            # :END_RPT
            deleteat!(program, index)
        elseif instruction.op == RPT
            # Set the start point for this repeat. Since we've removed all
            # END_RPT before this, we don't have to worry about this index
            # getting messed up.
            #
            # Add 1 so the start points to the next instruction.
            program[index].instruction = set(instruction, set_repeat_start(index + 1))

            # Mark the index of the last repeat function for fast setting.
            index_of_last_rpt = index
            index += 1
        else
            index += 1
        end
    end
end

function makelabeldict(labels)
    d = Dict{Symbol,Int}()
    for (index, label) in enumerate(labels)
        # Skip non-set labels
        label isa Nothing && continue

        # Check if a symbol is already in the dictionary. If so, do nothing.
        # Otherwise, mark this index.
        if !in(label, keys(d))
            d[label] = index
        end
    end
    return d
end


text:
abstract type AbstractPressureForce <: AbstractPhysicsComponent end

Base.@kwdef struct PressureDivergence{T} <: AbstractPressureForce
    eos::T
end

@inline calc_flux!(flux, ::Nothing, state, _...) = nothing
@inline calc_flux!(flux, ::AbstractPressureForce, state, _...) = nothing

@inline function calc_flux!(flux, pressure::PressureDivergence, state, aux, _...)
    eos = pressure.eos
    p = calc_pressure(eos, state, aux)

    flux.ρu += p * I

    return nothing
end

max_stars_repo_name: CliMA/ClimaCoupler


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

function compute_vel(a,v0,psi,sigma,tau,eta)
    compute_vel_ = load_op_and_grad("./build/libComputeVel","compute_vel")
a,v0,psi,sigma,tau,eta = convert_to_tensor([a,v0,psi,sigma,tau,eta], [Float64,Float64,Float64,Float64,Float64,Float64])
    compute_vel_(a,v0,psi,sigma,tau,eta)
end

# verify
n = 10
a = 2.0*ones(n)|>constant
x = 4.0*ones(n)|>constant
v0 = 0.001|>constant
Ψ = 2000.0*ones(n)|>constant
σ = 3.0*ones(n)|>constant
η = 2.0|>constant
s = x/2v0
inv_o = exp(-Ψ/a)
τ = a * (Ψ/a + log(s + sqrt(s*s+inv_o*inv_o)))*σ + η*x
# τ = a * asinh(*exp(Ψ/a))*σ + η*(x-u)/Δt

# TODO: specify your input parameters
u = compute_vel(a,v0,Ψ,σ,τ,η)
sess = Session(); init(sess)
@show run(sess, u)
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(m)
    # return sum(compute_vel(a,v0,Ψ,σ,τ,η)^2)
    return sum(compute_vel(a,v0,Ψ,σ,τ,m)^2)
end

# TODO: change `m_` and `v_` to appropriate values
# m_ = constant(rand(n))
# v_ = rand(n)

# m_ = τ
# v_ = rand(n)

m_ = constant(rand())
v_ = rand()

y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
module ODE

export evolve_ode, Midpoint

"""    stepsize_adjust(y, yerr, h, order; rtol=1e-8, atol=1e-8)

Returns `(hnew, retry_step)` giving a stepsize that will try to achieve the
given absolute and relative tolerances on the next step (or a retry of this
one).

`retry_step` is `true` if the error exceeds the bound by more than a factor of
two, and `false` otherwise.
"""
function stepsize_adjust(y, yerr, h, order; rtol=1e-8, atol=1e-8)
    scale_max = zero(y[1])
    for i in 1:length(y)
        scale = yerr[i]/(atol + rtol*abs(y[i]))
        scale_max = max(scale_max, scale)
    end

    retry_step = (scale_max > 2 ? true : false)

    hnew = h/(scale_max)^(1/(order+1))

    if hnew > 2*h
        hnew = 2*h
    elseif hnew < 1/2*h
        hnew = 1/2*h
    end

    (hnew, retry_step)
end

struct Stepper
    order::Int
    step
end

"""    evolve_ode(stepper, f_rhs, x0, y0, x1; h_init, atol, rtol)

Returns the state of the ODE represented by `f_rhs` evolved from its initial
condition `(x0, y0)` to `x1` using the single-step algorithm in `stepper`.

`h_init` is used to guess the first step.

The evolution is adjusted to control the the absolute and relative error in the
solution according to the given absolute and relative tolerances; steps may be
re-tried until an acceptable error is achieved.

The final step will land *exactly* on `x1`.
"""
function evolve_ode(stepper, f_rhs, x0, y0, x1; h_init=(x1-x0)/100, atol=1e-8, rtol=1e-8)
    x = x0
    y = y0
    h = h_init
    not_done = true

    while not_done
        if x + h > x1
            h = x1-x
            not_done = false
        end

        y_new, yerr = stepper.step(f_rhs, x, y, h)
        h_new, retry = stepsize_adjust(y_new, yerr, h, stepper.order; rtol=rtol, atol=atol)

        if retry
            h = h_new
            not_done = true
        else
            y = y_new
            x = x + h
            h = h_new
        end
    end

    y
end

function step_midpoint(f_rhs, x, y, h)
    dydx0 = f_rhs(x, y)
    h_half = h/2
    y_half = y .+ dydx0 .* h_half

    dydx_half = f_rhs(x+h_half, y_half)

    y_new = y .+ dydx_half .* h
    yerr = h .* abs.(dydx_half .- dydx0)

    (y_new, yerr)
end

Midpoint = Stepper(2, step_midpoint)

end # module


text:
type Sweet{T<:Real}
    sweetness :: T
end

issweet(A) = false
issweet(A::Sweet) = A.sweetness > 0


text:
# Adapted from ClimateMachine.jl
# Original code by jkozdon

module ClimaCoreVTK

export writevtk

using WriteVTK
import ClimaCore: Fields, Geometry, Spaces, Topologies, Operators

include("space.jl")
include("addfield.jl")

"""
    vtk_grid(basename, gridspace::ClimaCore.Spaces.AbstractSpace;
        basis=:cell, vtkargs...)

Construct a VTK grid from a `ClimaCore.Spaces.AbstractSpace`. If
`basis=:lagrange`, it will construct a mesh made of Lagrange cells (valid only
for spectral element spaces), otherwise it will it subdivide the space into
quads, with vertices at nodal points.
"""
function WriteVTK.vtk_grid(
    basename::String,
    gridspace::Spaces.AbstractSpace;
    basis,
    vtkargs...,
)

    if basis == :lagrange
        cells = vtk_cells_lagrange(gridspace)
    else
        cells = vtk_cells_linear(gridspace)
    end
    cart_coords =
        Geometry.Cartesian123Point.(
            Fields.coordinate_field(gridspace),
            Ref(gridspace.global_geometry),
        )
    return vtk_grid(
        basename,
        vec(parent(cart_coords.x1)),
        vec(parent(cart_coords.x2)),
        vec(parent(cart_coords.x3)),
        cells;
        vtkargs...,
    )
end

function vtk_file(
    basename::String,
    fields,
    ;
    gridspace::Spaces.AbstractSpace = vtk_grid_space(fields),
    basis = :cell,
    dataspace::Spaces.AbstractSpace = basis == :cell ?
                                      vtk_cell_space(gridspace) : gridspace,
    vtkargs...,
)
    vtkfile = vtk_grid(basename, gridspace; basis, vtkargs...)
    addfield!(vtkfile, nothing, fields, dataspace)
    return vtkfile
end

"""
    writevtk(
        basename::String,
        fields;
        basis=:cell,
        vtkargs...
    )

Write `fields` to as an unstructured mesh VTK file named `$(basename).vtu`.

`fields` can be either:
- a ClimaCore `Field` object,
- a `FieldVector` object,
- a `NamedTuple` of `Field`s.

The `basis` keyword option determines the type of cells used to write.:
- `:cell` (default): output values at cell centers (interpolating where
  necessary).
- `:point` (default): output values at cell vertices.
- `:lagrange`: output values at Lagrange nodes (valid only for spectral element
  spaces), using Use VTK Lagrange cells to accurately represent high-order
  elements.

Any additional keyword arguments are passed to
[`WriteVTK.vtk_grid`](https://jipolanco.github.io/WriteVTK.jl/stable/grids/syntax/#Supported-options).
"""
function writevtk(basename::String, fields; vtkargs...)
    vtkfile = vtk_file(basename, fields; vtkargs...)
    vtk_save(vtkfile)
end

"""
    writevtk(
        basename::String,
        times,
        fields;
        vtkargs...
    )

Write a sequence of fields at `times` as a Paraview collection.

`fields` can be either be an iterable collection of fields, or a `NamedTuple` of collections.
"""
function writevtk(basename::String, times, fields; vtkargs...)
    npad = ndigits(length(times); pad = 3)

    paraview_collection(basename) do pvd
        for (n, time) in enumerate(times)
            pvd[time] = vtk_file(
                basename * "_" * string(n, pad = npad),
                fields[n];
                vtkargs...,
            )
        end
    end
end
function writevtk(basename::String, times, fields::NamedTuple, vtkargs...)
    npad = ndigits(length(times); pad = 3)

    paraview_collection(basename) do pvd
        for (n, time) in enumerate(times)
            field = NamedTuple(key => val[n] for (key, val) in pairs(fields))
            pvd[time] = vtk_file(
                basename * "_" * string(n, pad = npad),
                field;
                vtkargs...,
            )
        end
    end
end

end # module


text:
include("../utility.jl")
using Distributions, Turing, Test

N = 5
setchunksize(4*N)
alg = HMC(1000, 0.2, 4)

@model vdemo() = begin
  v = Vector{Matrix{Real}}(undef, N)
  v ~ [Wishart(7, [1 0.5; 0.5 1])]
end

t_vec = @elapsed res_vec = sample(vdemo(), alg)

@model vdemo() = begin
  v = Vector{Matrix{Real}}(undef, N)
  for i = 1:N
    v[i] ~ Wishart(7, [1 0.5; 0.5 1])
  end
end

t_loop = @elapsed res = sample(vdemo(), alg)

println("Time for")
println("  Loop : $t_loop")
println("  Vec  : $t_vec")


text:
using IcanProj
using NetCDF
using DataFrames
using ProgressMeter
using PyCall
using CSV

function project_results(df::DataFrame, on::Symbol)

    map_values = fill(NaN, 1550, 1195)

    map_values[df[:ind_julia]] = df[on]

    return map_values

end

function add_variability(df_all)

    variables = ["ilwr", "iswr", "pres", "rainf", "rhum", "snowf", "tair", "wind", "prec"]

    for variable in variables

        namemean = Symbol(variable * "_mean")

        namestd = Symbol(variable * "_std")

        df_tmp = by(df_all, :ind_50km, d -> DataFrame(data_std = std(d[namemean])))

        names!(df_tmp, [:ind_50km, namestd])

        df_all = join(df_all, df_tmp, on = :ind_50km)

    end

    return df_all

end

function plot_map(df, variable, cb_label, cb_unit, title_label, figpath)

    resmap = project_results(df, variable)

    file = joinpath(figpath, String(variable) * ".png")

    cb_text =  cb_label * " " * cb_unit

    py"""

    import matplotlib.pyplot as plt

    fig, ax = plt.subplots()

    #fig.figsize = (8, 12)

    ax.set_title($(title_label))

    img_plot = ax.imshow($(resmap), cmap = 'jet')

    cbar = fig.colorbar(img_plot)

    cbar.set_label($(cb_text))

    ax.grid(linestyle='dotted')

    plt.tick_params(axis='both', left='off', top='off', right='off', bottom='off',
                    labelleft='off', labeltop='off', labelright='off', labelbottom='off')

    fig.savefig($(file), dpi = 200)

    plt.close('all')

    """

end

df_all = CSV.File(joinpath(dirname(pathof(IcanProj)), "..", "data", "forcings_summary.txt")) |> DataFrame

df_all[:prec_mean] = min.(365*8(df_all[:rainf_mean] + df_all[:snowf_mean]), 4500)

df_all = add_variability(df_all)

figpath = joinpath(dirname(pathof(IcanProj)), "..", "plots", "forcings")

variable = ["ilwr",
            "iswr",
            "pres",
            "rainf",
            "rhum",
            "snowf",
            "tair",
            "wind",
            "prec"]

cb_unit = ["(W/m2)",
           "(W/m2)",
           "(kPa)",
           "(mm/timestep)",
           "(%)",
           "(mm/timestep)",
           "(C)",
           "(m/s)",
           "(mm/year)"]

title_label = ["Incoming longwave radiation",
               "Incoming shortwave radiation",
               "Air pressure",
               "Rainfall",
               "Relative humidity",
               "Snowfall",
               "Air temperature",
               "Wind speed",
               "Precipitation"]

# Plot standard deviations

cb_label = "Standard deviation"

variable_std = convert.(Symbol, variable .* "_std")

for (v, cu, t) in zip(variable_std, cb_unit, title_label)

    plot_map(df_all, v, cb_label, cu, t, figpath)

end

# Plot averages

cb_label = "Average"

variable_mean = convert.(Symbol, variable .* "_mean")

for (v, cu, t) in zip(variable_mean, cb_unit, title_label)

    plot_map(df_all, v, cb_label, cu, t, figpath)

end