text:
order(crs::AbstractAlgebraicCoarsening) = crs.order

function _add_to_neighborhood!(N, k, v, maxsize)
    push!(N, (k, v)) # TODO Smarter implementation
    if length(N) > maxsize
        sort!(N; by=x->x[1], rev=true)
        pop!(N)
    end
    return nothing
end

function coarse_neighborhoods(W, C, invseed, order)
    nc = length(C)
    nf = length(invseed) - length(C)
    # TODO Other strategies for sparsifying neighborhood?
    N = [Vector{Tuple{Float64, Int}}() for i in 1:nf]
    rows = rowvals(W)
    vals = nonzeros(W)

    for (j, c) in enumerate(C)
        for idx in nzrange(W, c)
            row = rows[idx]
            val = vals[idx]
            if val != 0 && invseed[row] > nc
                idxf = invseed[row] - nc
                _add_to_neighborhood!(N[idxf], val, j, order)
            end
        end
    end
    return N
end

function _fillcoarseop_coarse!(Pi, Pj, Pv, C, idx)
    for (i, c) in enumerate(C)
        idx += 1
        Pi[idx] = c
        Pj[idx] = i
        Pv[idx] = 1.0
    end
    return idx
end

function _fillcoarseop_fine!(Pi, Pj, Pv, F, N, idx)
    for (i, f) in enumerate(F)
        total = sum(x -> x[1], N[i])
        for (w, c) in N[i]
            idx += 1
            Pi[idx] = f
            Pj[idx] = c
            Pv[idx] = w / total
        end
    end
    return idx
end

function _fillcoarseop!(Pi, Pj, Pv, C, F, N)
    idx = 0
    idx = _fillcoarseop_coarse!(Pi, Pj, Pv, C, idx)
    idx = _fillcoarseop_fine!(Pi, Pj, Pv, F, N, idx)
    return nothing
end

function formcoarseop(crs::AbstractAlgebraicCoarsening, C, F, N)
    nc = length(C)
    nf = length(F)
    k = sum(length, N) + nc
    Pi = Vector{Int}(undef, k)
    Pj = Vector{Int}(undef, k)
    Pv = Vector{Float64}(undef, k)
    _fillcoarseop!(Pi, Pj, Pv, C, F, N)
    return sparse(Pi, Pj, Pv, nc + nf, nc)
end

function fix_adjacency!(A)
    for i in 1:size(A, 1)
        A[i, i] = 0
    end
    return A
end

function fix_laplacian!(L)
    rows = rowvals(L)
    vals = nonzeros(L)
    for c in axes(L, 2)
        acc = 0
        for idx in nzrange(L, c)
            if rows[idx] != c
                acc -= vals[idx]
            end
        end
        L[c, c] = acc
    end
    return L
end


text:
# import Base.Dates: UTInstant, DateTime, TimeZone, Millisecond
using Base.Dates
import Base.Dates: value
import Base: ==, isequal, isless

abstract TimeError <: Exception

type AmbiguousTimeError <: TimeError
    local_dt::DateTime
    timezone::TimeZone
end

type NonExistentTimeError <: TimeError
    local_dt::DateTime
    timezone::TimeZone
end

function Base.showerror(io::IO, e::AmbiguousTimeError)
    print(io, "AmbiguousTimeError: ")
    print(io, "Local DateTime $(e.local_dt) is ambiguous within $(string(e.timezone))")
end

function Base.showerror(io::IO, e::NonExistentTimeError)
    print(io, "NonExistentTimeError: ")
    print(io, "Local DateTime $(e.local_dt) does not exist within $(string(e.timezone))")
end

# Using type Symbol instead of AbstractString for name since it
# gets us ==, and hash for free.

"""
    FixedTimeZone

A `TimeZone` with a constant offset for all of time.
"""
immutable FixedTimeZone <: TimeZone
    name::Symbol
    offset::UTCOffset
end

"""
    FixedTimeZone(name, utc_offset, dst_offset=0) -> FixedTimeZone

Constructs a `FixedTimeZone` with the given `name`, UTC offset (in seconds), and DST offset
(in seconds).
"""
function FixedTimeZone(name::AbstractString, utc_offset::Integer, dst_offset::Integer=0)
    FixedTimeZone(Symbol(name), UTCOffset(utc_offset, dst_offset))
end

function FixedTimeZone(name::AbstractString, utc_offset::Second, dst_offset::Second=Second(0))
    FixedTimeZone(Symbol(name), UTCOffset(utc_offset, dst_offset))
end

"""
    FixedTimeZone(::AbstractString) -> FixedTimeZone

Constructs a `FixedTimeZone` with a UTC offset from a string. Resulting `FixedTimeZone` will
be named like \"UTC±HH:MM[:SS]\".

# Examples
```julia
julia> FixedTimeZone(\"UTC+6\")
UTC+06:00

julia> FixedTimeZone(\"-1330\")
UTC-13:30

julia> FixedTimeZone(\"15:45:21\")
UTC+15:45:21
```
"""
function FixedTimeZone(s::AbstractString)
    const regex = r"""
    ^(?|
        UTC([+-]\d{1,2})?
        |
        (?:UTC(?=[+-]))?
        ([+-]?\d{2})
        (?|
            (\d{2})
            |
            \:(\d{2})
            (?:\:(\d{2}))?
        )
    )$
    """x

    m = match(regex, s)
    m == nothing && throw(ArgumentError("Unrecognized time zone: $s"))

    values = map(n -> n == nothing ? 0 : Base.parse(Int, n), m.captures)

    if values == [0, 0, 0]
        name = "UTC"
    elseif values[3] == 0
        name = @sprintf("UTC%+03d:%02d", values[1:2]...)
    else
        name = @sprintf("UTC%+03d:%02d:%02d", values...)
    end

    if values[1] < 0
        for i in 2:length(values)
            values[i] = -values[i]
        end
    end

    offset = values[1] * 3600 + values[2] * 60 + values[3]
    return FixedTimeZone(name, offset)
end

immutable Transition
    utc_datetime::DateTime  # Instant where new zone applies
    zone::FixedTimeZone
end

Base.isless(x::Transition,y::Transition) = isless(x.utc_datetime,y.utc_datetime)

"""
    VariableTimeZone

A `TimeZone` with an offset that changes over time.
"""
immutable VariableTimeZone <: TimeZone
    name::Symbol
    transitions::Vector{Transition}
    cutoff::Nullable{DateTime}
end

function VariableTimeZone(name::AbstractString, transitions::Vector{Transition}, cutoff::Nullable{DateTime})
    return VariableTimeZone(Symbol(name), transitions, cutoff)
end

function VariableTimeZone(name::AbstractString, transitions::Vector{Transition}, cutoff::DateTime)
    return VariableTimeZone(Symbol(name), transitions, Nullable(cutoff))
end

function VariableTimeZone(name::AbstractString, transitions::Vector{Transition})
    return VariableTimeZone(Symbol(name), transitions, Nullable{DateTime}())
end

type UnhandledTimeError <: TimeError
    tz::VariableTimeZone
end
Base.showerror(io::IO, e::UnhandledTimeError) = print(io, "TimeZone $(string(e.tz)) does not handle dates on or after $(get(e.tz.cutoff)) UTC")

# """
#     ZonedDateTime

# A `DateTime` that includes `TimeZone` information.
# """

immutable ZonedDateTime <: TimeType
    utc_datetime::DateTime
    timezone::TimeZone
    zone::FixedTimeZone  # The current zone for the utc_datetime.

    function ZonedDateTime(utc_datetime::DateTime, timezone::TimeZone, zone::FixedTimeZone)
        return new(utc_datetime, timezone, zone)
    end

    function ZonedDateTime(utc_datetime::DateTime, timezone::VariableTimeZone, zone::FixedTimeZone)
        if utc_datetime >= get(timezone.cutoff, typemax(DateTime))
            throw(UnhandledTimeError(timezone))
        end

        return new(utc_datetime, timezone, zone)
    end
end

"""
    ZonedDateTime(dt::DateTime, tz::TimeZone; from_utc=false) -> ZonedDateTime

Construct a `ZonedDateTime` by applying a `TimeZone` to a `DateTime`. When the `from_utc`
keyword is true the given `DateTime` is assumed to be in UTC instead of in local time and is
converted to the specified `TimeZone`.  Note that when `from_utc` is true the given
`DateTime` will always exists and is never ambiguous.
"""
function ZonedDateTime(dt::DateTime, tz::VariableTimeZone; from_utc::Bool=false)
    possible = interpret(dt, tz, from_utc ? UTC : Local)

    num = length(possible)
    if num == 1
        return first(possible)
    elseif num == 0
        throw(NonExistentTimeError(dt, tz))
    else
        throw(AmbiguousTimeError(dt, tz))
    end
end

function ZonedDateTime(dt::DateTime, tz::FixedTimeZone; from_utc::Bool=false)
    utc_dt = from_utc ? dt : dt - tz.offset
    return ZonedDateTime(utc_dt, tz, tz)
end

"""
    ZonedDateTime(dt::DateTime, tz::VariableTimeZone, occurrence::Integer) -> ZonedDateTime

Construct a `ZonedDateTime` by applying a `TimeZone` to a `DateTime`. If the `DateTime` is
ambiguous within the given time zone you can set `occurrence` to a positive integer to
resolve the ambiguity.
"""
function ZonedDateTime(dt::DateTime, tz::VariableTimeZone, occurrence::Integer)
    possible = interpret(dt, tz, Local)

    num = length(possible)
    if num == 1
        return first(possible)
    elseif num == 0
        throw(NonExistentTimeError(dt, tz))
    elseif occurrence > 0
        return possible[occurrence]
    else
        throw(AmbiguousTimeError(dt, tz))
    end
end

"""
    ZonedDateTime(dt::DateTime, tz::VariableTimeZone, is_dst::Bool) -> ZonedDateTime

Construct a `ZonedDateTime` by applying a `TimeZone` to a `DateTime`. If the `DateTime` is
ambiguous within the given time zone you can set `is_dst` to resolve the ambiguity.
"""
function ZonedDateTime(dt::DateTime, tz::VariableTimeZone, is_dst::Bool)
    possible = interpret(dt, tz, Local)

    num = length(possible)
    if num == 1
        return first(possible)
    elseif num == 0
        throw(NonExistentTimeError(dt, tz))
    elseif num == 2
        mask = [isdst(zdt.zone.offset) for zdt in possible]

        # Mask is expected to be unambiguous.
        !($)(mask...) && throw(AmbiguousTimeError(dt, tz))

        occurrence = is_dst ? findfirst(mask) : findfirst(!mask)
        return possible[occurrence]
    else
        throw(AmbiguousTimeError(dt, tz))
    end
end

# Convenience constructors
@doc """
    ZonedDateTime(y, [m, d, h, mi, s, ms], tz, [amb]) -> DateTime

Construct a `ZonedDateTime` type by parts. Arguments `y, m, ..., ms` must be convertible to
`Int64` and `tz` must be a `TimeZone`. If the given `DateTime` is ambiguous in the given
`TimeZone` then `amb` can be supplied to resolve ambiguity.
""" ZonedDateTime

@optional function ZonedDateTime(y::Integer, m::Integer=1, d::Integer=1, h::Integer=0, mi::Integer=0, s::Integer=0, ms::Integer=0, tz::VariableTimeZone, amb::Union{Integer,Bool})
    ZonedDateTime(DateTime(y,m,d,h,mi,s,ms), tz, amb)
end

@optional function ZonedDateTime(y::Integer, m::Integer=1, d::Integer=1, h::Integer=0, mi::Integer=0, s::Integer=0, ms::Integer=0, tz::TimeZone)
    ZonedDateTime(DateTime(y,m,d,h,mi,s,ms), tz)
end

function ZonedDateTime(parts::Union{Period,TimeZone}...)
    periods = Period[]
    timezone = Nullable{TimeZone}()
    for part in parts
        if isa(part, Period)
            push!(periods, part)
        elseif isnull(timezone)
            timezone = Nullable{TimeZone}(part)
        else
            throw(ArgumentError("Multiple time zones found"))
        end
    end

    isnull(timezone) && throw(ArgumentError("Missing time zone"))
    return ZonedDateTime(DateTime(periods...), get(timezone))
end

# Equality
==(a::ZonedDateTime, b::ZonedDateTime) = a.utc_datetime == b.utc_datetime
isless(a::ZonedDateTime, b::ZonedDateTime) = isless(a.utc_datetime, b.utc_datetime)

# Note: a more efficient version of: `hash(a) == hash(b)`. Assumes that the "zone" of the
# ZonedDateTimes is not being set incorrectly.
function isequal(a::ZonedDateTime, b::ZonedDateTime)
    isequal(a.utc_datetime, b.utc_datetime) && isequal(a.timezone, b.timezone)
end

function ==(a::VariableTimeZone, b::VariableTimeZone)
    a.name == b.name && a.transitions == b.transitions
end

ext: jl


text:
# # Logistic regression
using DataFrames
using Plots
using RDatasets
using Convex
using SCS

#-

## we'll use iris data
## predict whether the iris species is versicolor using the sepal length and width and petal length and width
iris = dataset("datasets", "iris")
## outcome variable: +1 for versicolor, -1 otherwise
Y = [species == "versicolor" ? 1.0 : -1.0 for species in iris.Species]
## create data matrix with one column for each feature (first column corresponds to offset)
X = hcat(ones(size(iris, 1)), iris.SepalLength, iris.SepalWidth, iris.PetalLength, iris.PetalWidth);

#-

## solve the logistic regression problem
n, p = size(X)
beta = Variable(p)
problem = minimize(logisticloss(-Y.*(X*beta)))

solve!(problem, SCSSolver(verbose=false))

#-

## let's see how well the model fits
using Plots
logistic(x::Real) = inv(exp(-x) + one(x))
perm = sortperm(vec(X*beta.value))
plot(1:n, (Y[perm] .+ 1)/2, st=:scatter)
plot!(1:n, logistic.(X*beta.value)[perm])


text:
__precompile__()

module JUnitTestSets

export JUnitTestSet, write_report, check_test_results, @jtestset, @jtest

using Test, Dates
include("JTestResult.jl")

struct JUnitTestSet
    id::String
    name::String
    tests::Int
    failures::Int
    time::Dates.Nanosecond
    results::Vector{JTestResult}
end

function evaluate_tests(arr)::JUnitTestSet
    evals::Vector{JTestResult} = []
    for a in arr
        if typeof(a) == JUnitTestSets.JTestResult
            push!(evals, eval(a))
        end
    end
    tests = length(evals)
    failures = length(filter(x -> !x.result, evals))
    etime = sum(map(x -> x.time, evals))
    return JUnitTestSet("tests", "Julia Tests", tests, failures, etime, evals)
end

function check_test_results(ts::JUnitTestSet)
    # Return 1 if any tests have failed, otherwise return 0. Handy for CIgit status

    for t in ts.results
        if !t.result
            return 1
        end
    end
    return 0
end

macro jtestset(args...)
    isempty(args) && error("No arguments to @testset")
    tests = Array(args[end].args)
    quote
        evaluate_tests([$(esc.(tests)...)])
    end
end

include("xml.jl")

end # module


text:
export HSegment, RandomHSegment, endpoints, collinear, bisector

"""
`HSegment(A,B)` creates a new line segment with endpoints `A` and `B`.
"""
struct HSegment <: HObject
    A::HPoint
    B::HPoint
    attr::Dict{Symbol,Any}
    function HSegment(a::HPoint,b::HPoint)
        @assert a!=b "The end points of a segment must be distinct"
        S = new(a,b,Dict{Symbol,Any}())
        set_color(S)
        set_thickness(S)
        set_line_style(S)
        return S
    end
end
HSegment(S::HSegment) = HSegment(S.A, S.B)  # copy constructor
HSegment(a::Number, B::HPoint) = HSegment(HPoint(a), B)
HSegment(A::HPoint, b::Number) = HSegment(A, HPoint(b))
HSegment(a::Number, b::Number) = HSegment(HPoint(a), HPoint(b))
HSegment(A::HPoint) = HSegment(A, HPoint())
HSegment(a::Number) = HSegment(HPoint(a), HPoint())

(+)(A::HPoint, B::HPoint) = HSegment(A,B)

function show(io::IO, L::HSegment)
    p,q = endpoints(L)
    a = getz(p)
    b = getz(q)
    print(io,"HSegment($a,$b)")
end

endpoints(L::HSegment) = [L.A, L.B]

RandomHSegment() = HSegment(RandomHPoint(), RandomHPoint())

function (==)(L::HSegment, LL::HSegment)
    return (endpoints(L)==endpoints(LL)) || (endpoints(L)==(LL.B,LL.A))
end

midpoint(L::HSegment) = midpoint(endpoints(L)...)
length(L::HSegment) = dist(endpoints(L)...)

"""
`collinear` checks if the arguments are collinear. Arguments are:
+ `a,b,c`: three points
+ `a,L`: point and segment (in either order)
+ `L,LL`: two segments
"""
function collinear(a::HPoint, b::HPoint, c::HPoint)::Bool
    if a==b || b==c || a==c
        return true
    end
    aa = reflect_across(a,HSegment(b,c))
    return aa == a
end

collinear(a::HPoint, L::HSegment) = collinear(a,endpoints(L)...)
collinear(L::HSegment, a::HPoint) = collinear(a,L)

function collinear(L::HSegment, LL::HSegment)
    a,b = endpoints(LL)
    return collinear(a,L) && collinear(b,L)
end

function in(a::HPoint, L::HSegment)
    x,y = endpoints(L)
    return between(x,a,y)
end

"""
`bisector(S::HSegment)` yields an `HLine` that's the perpendicular
bisector of the segment `S`. May also be invoked `bisector(A,B)`
where `A` and `B` are points.
"""
function bisector(S::HSegment)
    A,B = endpoints(S)
    @assert A != B "The end points of the segment must be distinct"
    P = midpoint(S)
    f = move2xplus(P,B)
    g = inv(f)
    zs = g(im)
    zt = g(-im)

    return HLine(angle(zs), angle(zt))
end

bisector(A::HPoint, B::HPoint) = bisector(A+B)


text:
module Ariel

export read_QuestaalXYZ

println("Ariel.jl")
println("Of sulphurous roaring the most mighty Neptune")
println("Seem to besiege and make his bold waves tremble,")

"""
 read_QuestaalXYZ(f::IOStream)

 Read in Questaal .XYZ file format to a trajectory object.

XYZ format looks like (128 H2O):
384 PLAT:   29.73674  0  0  0  29.73674  0  0  0  29.73674
t=103460fs T=285K V=-113.146Ry K.E.+V=-112.1104Ry alat=1; velocities times 1000
O            2.69595    -2.92790    13.17485    0.54262  -0.46064  -0.19684    0.00826   0.03806  -0.06875    1.01293   0.03968   0.25365  -0.09923   0.00000   0.00000   0.00000   0.00000   0.00000
"""
function read_QuestaalXYZ(f::IOStream)
    h=split(readline(f))
    natoms=parse(Int,h[1])
    println("natoms: $natoms")

    PLAT=reshape(map(x->parse(Float64,x), h[3:11] ), 3,3) # lattice parameters
    show(PLAT)

    readline(f) # don't do anything with the comment line, currently

#    traj=[]
#    while !eof(f)
#        append(traj,read_next_frame(f),natoms)
#    end
    A=reshape(readdlm(f, skipstart=2), natoms, 19, 2)
    show(A)
end

"Gobble natoms worth of frame"
function read_next_frame(f::IOStream,natoms)
    println("Oh lordy.")
end

end # module


text:
module CreateTableBitempWebApp
import SearchLight: AbstractModel, DbId, query, save!
import SearchLight.Migrations:
    create_table, column, columns, primary_key, add_index, drop_table, add_indices
import SearchLightPostgreSQL
import Base: @kwdef
import Intervals, Dates, TimeZones
using Intervals,
    Dates, SearchLight, SearchLight.Transactions, SearchLightPostgreSQL, TimeZones
export up, down

function createRevisionsTriggerAndConstraint(
    trigger::Symbol,
    constraint::Symbol,
    table::Symbol,
)
    SearchLight.query("""
                      CREATE TRIGGER $trigger
                      BEFORE INSERT OR UPDATE ON $table
                      FOR EACH ROW EXECUTE PROCEDURE f_versionrange();
                      """)
    SearchLight.query(
        """
        ALTER TABLE $table
        ADD CONSTRAINT $constraint EXCLUDE USING GIST (ref_component WITH =, ref_valid WITH &&)
        """,
    )
end

function up()

    create_table(:contracts) do
        [
            column(:id, :bigserial, "PRIMARY KEY")
            column(:ref_history, :bigint, "REFERENCES histories(id) ON DELETE CASCADE")
            column(:ref_version, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
        ]
    end

    create_table(:contractRevisions) do
        [
            column(:id, :bigserial, "PRIMARY KEY")
            column(:ref_component, :bigint, "REFERENCES contracts(id) ON DELETE CASCADE")
            column(:ref_validfrom, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
            column(:ref_invalidfrom, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
            column(:ref_valid, :int8range)
            column(:description, :string)
        ]
    end

    createRevisionsTriggerAndConstraint(
        :cr_versions_trig,
        :contractsversionrange,
        :contractRevisions,
    )

    create_table(:productItems) do
        [
            column(:id, :bigserial, "PRIMARY KEY")
            column(:ref_history, :bigint, "REFERENCES histories(id) ON DELETE CASCADE")
            column(:ref_version, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
            column(:ref_super, :bigint, "REFERENCES contracts(id) ON DELETE CASCADE")
        ]
    end

    create_table(:productItemRevisions) do
        [
            column(:id, :bigserial, "PRIMARY KEY")
            column(
                :ref_component,
                :bigint,
                "REFERENCES productitems(id) ON DELETE CASCADE",
            )
            column(:position, :bigint)
            column(:ref_validfrom, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
            column(:ref_invalidfrom, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
            column(:ref_valid, :int8range)
            column(:description, :string)
        ]
    end

    createRevisionsTriggerAndConstraint(
        :pi_versions_trig,
        :pi_versionrange,
        :productitemRevisions,
    )

    create_table(:partners) do
        [
            column(:id, :bigserial, "PRIMARY KEY")
            column(:ref_history, :bigint, "REFERENCES histories(id) ON DELETE CASCADE")
            column(:ref_version, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
        ]
    end

    create_table(:partnerRevisions) do
        [
            column(:id, :bigserial, "PRIMARY KEY")
            column(:ref_component, :bigint, "REFERENCES partners(id) ON DELETE CASCADE")
            column(:ref_validfrom, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
            column(:ref_invalidfrom, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
            column(:ref_valid, :int8range)
            column(:description, :string)
        ]
    end

    createRevisionsTriggerAndConstraint(
        :pr_versions_trig,
        :pr_versionrange,
        :partnerRevisions,
    )

    create_table(:contractPartnerRefs) do
        [
            column(:id, :bigserial, "PRIMARY KEY")
            column(:ref_history, :bigint, "REFERENCES histories(id) ON DELETE CASCADE")
            column(:ref_version, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
            column(:ref_super, :bigint, "REFERENCES contracts(id) ON DELETE CASCADE")
        ]
    end

    create_table(:contractPartnerRefRevisions) do
        [
            column(:id, :bigserial, "PRIMARY KEY")
            column(
                :ref_component,
                :bigint,
                "REFERENCES contractPartnerRefs(id) ON DELETE CASCADE",
            )
            column(:ref_validfrom, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
            column(:ref_invalidfrom, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
            column(:ref_valid, :int8range)
            column(:description, :string)
            column(:ref_partner, :bigint, "REFERENCES partners(id) ON DELETE CASCADE")
        ]
    end

    createRevisionsTriggerAndConstraint(
        :cprr_versions_trig,
        :cprr_versionrange,
        :contractPartnerRefRevisions,
    )

    create_table(:tariffs) do
        [
            column(:id, :bigserial, "PRIMARY KEY")
            column(:ref_history, :bigint, "REFERENCES histories(id) ON DELETE CASCADE")
            column(:ref_version, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
        ]
    end

    create_table(:tariffRevisions) do
        [
            column(:id, :bigserial, "PRIMARY KEY")
            column(:ref_component, :bigint, "REFERENCES tariffs(id) ON DELETE CASCADE")
            column(:ref_validfrom, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
            column(:ref_invalidfrom, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
            column(:ref_valid, :int8range)
            column(:description, :string)
        ]
    end

    createRevisionsTriggerAndConstraint(
        :tr_versions_trig,
        :tr_versionrange,
        :tariffRevisions,
    )

    create_table(:productItemTariffRefs) do
        [
            column(:id, :bigserial, "PRIMARY KEY")
            column(:ref_history, :bigint, "REFERENCES histories(id) ON DELETE CASCADE")
            column(:ref_version, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
            column(:ref_super, :bigint, "REFERENCES productitems(id) ON DELETE CASCADE")
        ]
    end

    create_table(:productItemTariffRefRevisions) do
        [
            column(:id, :bigserial, "PRIMARY KEY")
            column(
                :ref_component,
                :bigint,
                "REFERENCES productItemTariffRefs(id) ON DELETE CASCADE",
            )
            column(:ref_validfrom, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
            column(:ref_invalidfrom, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
            column(:ref_valid, :int8range)
            column(:description, :string)
            column(:ref_tariff, :bigint, "REFERENCES tariffs(id) ON DELETE CASCADE")
        ]
    end

    createRevisionsTriggerAndConstraint(
        :pitrr_versions_trig,
        :pitrr_versionrange,
        :productItemTariffRefRevisions,
    )

    create_table(:productItemPartnerRefs) do
        [
            column(:id, :bigserial, "PRIMARY KEY")
            column(:ref_history, :bigint, "REFERENCES histories(id) ON DELETE CASCADE")
            column(:ref_version, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
            column(:ref_super, :bigint, "REFERENCES productitems(id) ON DELETE CASCADE")
        ]
    end

    create_table(:productItemPartnerRefRevisions) do
        [
            column(:id, :bigserial, "PRIMARY KEY")
            column(
                :ref_component,
                :bigint,
                "REFERENCES productItemPartnerRefs(id) ON DELETE CASCADE",
            )
            column(:ref_validfrom, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
            column(:ref_invalidfrom, :bigint, "REFERENCES versions(id) ON DELETE CASCADE")
            column(:ref_valid, :int8range)
            column(:description, :string)
            column(:ref_partner, :bigint, "REFERENCES partners(id) ON DELETE CASCADE")
        ]
    end

    createRevisionsTriggerAndConstraint(
        :piprr_versions_trig,
        :piprr_versionrange,
        :productItemPartnerRefRevisions,
    )

end

function down()
    BitemporalPostgres.DDL.down()
    drop_table(:contractRevisions)
    drop_table(:contracts)
    drop_table(:contractProductItems)
    drop_table(:contractProductItemRevisions)
    drop_table(:contractPartnerRefRevisions)
    drop_table(:contractPartnerRefs)
    drop_table(:productItemTariffRefRevisions)
    drop_table(:productItemTariffRefs)
    drop_table(:productItemPartnerRefRevisions)
    drop_table(:productItemPartnerRefs)
    drop_table(:partnerRevisions)
    drop_table(:partners)
    drop_table(:tariffRevisions)
    drop_table(:tariffs)

end
end

max_stars_repo_name: michaelfliegner/BitempWebApp

max_issues_repo_name: michaelfliegner/BitempWebApp


text:
using ArgParse
using BioSequences
using CodecZlib
using FASTX

function skew1(chrom, seq, w)
	for i = 1:length(seq) -w
		g = 0
		c = 0
		for j = i:i+w
			if     seq[j] == DNA_C c += 1
			elseif seq[j] == DNA_G g += 1
			end
		end
		if (g+c) == 0 println("$chrom\t$i\t", 0)
		else          println("$chrom\t$i\t", (g-c)/(g+c))
		end
	end
end

function skew2(chrom, seq, w)
	for i = 1:length(seq) -w
		comp = composition(seq[i:i+w])
		c = comp[DNA_C]
		g = comp[DNA_G]
		if (g+c) == 0 println("$chrom\t$i\t", 0)
		else          println("$chrom\t$i\t", (g-c)/(g+c))
		end
	end
end

function skew3(chrom, seq, w)
	# first window
end

cli = ArgParseSettings()
@add_arg_table cli begin
	"--window"
		help = "window size"
		arg_type = Int
		default = 100
	"--algorithm"
		help = "1-3"
		arg_type = Int
		default = 1
	"fasta"
		help = "fasta file, required"
		required = true
end
arg = parse_args(ARGS, cli)

reader = nothing
if splitext(arg["fasta"])[2] == ".gz"
	reader = FASTA.Reader(GzipDecompressorStream(open(arg["fasta"])))
else
	reader = FASTA.Reader(open(arg["fasta"]))
end

## Main loop ##
w = arg["window"]
a = arg["algorithm"]
for record in reader
	c = FASTA.identifier(record)
	s = FASTA.sequence(record)
	if     a == 1 skew1(c, s, w)
	elseif a == 2 skew2(c, s, w)
	elseif a == 3 skew3(c, s, w)
	end
end

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
# Bernoulli distribution

export Bernoulli
import Base

@parameterized Bernoulli(p) ≃ CountingMeasure(ℤ[0:1])

function logdensity(d::Bernoulli{(:p,)}, y)
    p = d.p
    return y * log(p) + (1 - y) * log(1 - p)
end

function density(d::Bernoulli{(:p,)}, y)
    p = d.p
    return 2*p*y - p - y + 1
end

function logdensity(d::Bernoulli{(:logitp,)}, y)
    x = d.logitp
    return y * x - log1pexp(x)
end

function density(d::Bernoulli{(:logitp,)}, y)
    exp_x = exp(d.logitp)
    return exp_x ^ y / (1 + exp_x)
end

sampletype(::Bernoulli) = Bool

Base.rand(rng::AbstractRNG, T::Type, d::Bernoulli{(:p,)}) = rand(rng, T) < d.p

Base.rand(rng::AbstractRNG, T::Type, d::Bernoulli{(:logitp,)}) = rand(rng, T) < logistic(d.logitp)

≪(::Bernoulli, ::IntegerRange{lo,hi}) where {lo, hi} = lo ≤ 0 && 1 ≤ hi

distproxy(d::Bernoulli{(:p,)}) = Dists.Bernoulli(d.p)
distproxy(d::Bernoulli{(:logitp,)}) = Dists.Bernoulli(logistic(d.logitp))


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct Explainable <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Explainable(; kwargs...)
        obj = new(meta(Explainable), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Explainable
const __meta_Explainable = Ref{ProtoMeta}()
function meta(::Type{Explainable})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Explainable)
            __meta_Explainable[] = target = ProtoMeta(Explainable)
            allflds = Pair{Symbol,Union{Type,String}}[:conjunction => AbstractString, :id => Int64]
            meta(target, Explainable, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Explainable[]
    end
end
function Base.getproperty(obj::Explainable, name::Symbol)
    if name === :conjunction
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :id
        return (obj.__protobuf_jl_internal_values[name])::Int64
    else
        getfield(obj, name)
    end
end

mutable struct Explainables_Owned_OwnedEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Explainables_Owned_OwnedEntry(; kwargs...)
        obj = new(meta(Explainables_Owned_OwnedEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Explainables_Owned_OwnedEntry (mapentry)
const __meta_Explainables_Owned_OwnedEntry = Ref{ProtoMeta}()
function meta(::Type{Explainables_Owned_OwnedEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Explainables_Owned_OwnedEntry)
            __meta_Explainables_Owned_OwnedEntry[] = target = ProtoMeta(Explainables_Owned_OwnedEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => Explainable]
            meta(target, Explainables_Owned_OwnedEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Explainables_Owned_OwnedEntry[]
    end
end
function Base.getproperty(obj::Explainables_Owned_OwnedEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::Explainable
    else
        getfield(obj, name)
    end
end

mutable struct Explainables_RelationsEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Explainables_RelationsEntry(; kwargs...)
        obj = new(meta(Explainables_RelationsEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Explainables_RelationsEntry (mapentry)
const __meta_Explainables_RelationsEntry = Ref{ProtoMeta}()
function meta(::Type{Explainables_RelationsEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Explainables_RelationsEntry)
            __meta_Explainables_RelationsEntry[] = target = ProtoMeta(Explainables_RelationsEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => Explainable]
            meta(target, Explainables_RelationsEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Explainables_RelationsEntry[]
    end
end
function Base.getproperty(obj::Explainables_RelationsEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::Explainable
    else
        getfield(obj, name)
    end
end

mutable struct Explainables_AttributesEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Explainables_AttributesEntry(; kwargs...)
        obj = new(meta(Explainables_AttributesEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Explainables_AttributesEntry (mapentry)
const __meta_Explainables_AttributesEntry = Ref{ProtoMeta}()
function meta(::Type{Explainables_AttributesEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Explainables_AttributesEntry)
            __meta_Explainables_AttributesEntry[] = target = ProtoMeta(Explainables_AttributesEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => Explainable]
            meta(target, Explainables_AttributesEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Explainables_AttributesEntry[]
    end
end
function Base.getproperty(obj::Explainables_AttributesEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::Explainable
    else
        getfield(obj, name)
    end
end

mutable struct Explainables_Owned <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Explainables_Owned(; kwargs...)
        obj = new(meta(Explainables_Owned), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Explainables_Owned
const __meta_Explainables_Owned = Ref{ProtoMeta}()
function meta(::Type{Explainables_Owned})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Explainables_Owned)
            __meta_Explainables_Owned[] = target = ProtoMeta(Explainables_Owned)
            allflds = Pair{Symbol,Union{Type,String}}[:owned => Base.Dict{AbstractString,Explainable}]
            meta(target, Explainables_Owned, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Explainables_Owned[]
    end
end
function Base.getproperty(obj::Explainables_Owned, name::Symbol)
    if name === :owned
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{AbstractString,Explainable}
    else
        getfield(obj, name)
    end
end

mutable struct Explainables_OwnershipsEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Explainables_OwnershipsEntry(; kwargs...)
        obj = new(meta(Explainables_OwnershipsEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Explainables_OwnershipsEntry (mapentry)
const __meta_Explainables_OwnershipsEntry = Ref{ProtoMeta}()
function meta(::Type{Explainables_OwnershipsEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Explainables_OwnershipsEntry)
            __meta_Explainables_OwnershipsEntry[] = target = ProtoMeta(Explainables_OwnershipsEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => Explainables_Owned]
            meta(target, Explainables_OwnershipsEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Explainables_OwnershipsEntry[]
    end
end
function Base.getproperty(obj::Explainables_OwnershipsEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::Explainables_Owned
    else
        getfield(obj, name)
    end
end

mutable struct Explainables <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Explainables(; kwargs...)
        obj = new(meta(Explainables), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Explainables
const __meta_Explainables = Ref{ProtoMeta}()
function meta(::Type{Explainables})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Explainables)
            __meta_Explainables[] = target = ProtoMeta(Explainables)
            allflds = Pair{Symbol,Union{Type,String}}[:relations => Base.Dict{AbstractString,Explainable}, :attributes => Base.Dict{AbstractString,Explainable}, :ownerships => Base.Dict{AbstractString,Explainables_Owned}]
            meta(target, Explainables, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Explainables[]
    end
end
function Base.getproperty(obj::Explainables, name::Symbol)
    if name === :relations
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{AbstractString,Explainable}
    elseif name === :attributes
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{AbstractString,Explainable}
    elseif name === :ownerships
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{AbstractString,Explainables_Owned}
    else
        getfield(obj, name)
    end
end

mutable struct ConceptMap_MapEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ConceptMap_MapEntry(; kwargs...)
        obj = new(meta(ConceptMap_MapEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct ConceptMap_MapEntry (mapentry)
const __meta_ConceptMap_MapEntry = Ref{ProtoMeta}()
function meta(::Type{ConceptMap_MapEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ConceptMap_MapEntry)
            __meta_ConceptMap_MapEntry[] = target = ProtoMeta(ConceptMap_MapEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => Concept]
            meta(target, ConceptMap_MapEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ConceptMap_MapEntry[]
    end
end
function Base.getproperty(obj::ConceptMap_MapEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::Concept
    else
        getfield(obj, name)
    end
end

mutable struct ConceptMap <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ConceptMap(; kwargs...)
        obj = new(meta(ConceptMap), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct ConceptMap
const __meta_ConceptMap = Ref{ProtoMeta}()
function meta(::Type{ConceptMap})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ConceptMap)
            __meta_ConceptMap[] = target = ProtoMeta(ConceptMap)
            allflds = Pair{Symbol,Union{Type,String}}[:map => Base.Dict{AbstractString,Concept}, :explainables => Explainables]
            meta(target, ConceptMap, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ConceptMap[]
    end
end
function Base.getproperty(obj::ConceptMap, name::Symbol)
    if name === :map
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{AbstractString,Concept}
    elseif name === :explainables
        return (obj.__protobuf_jl_internal_values[name])::Explainables
    else
        getfield(obj, name)
    end
end

mutable struct ConceptMapGroup <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ConceptMapGroup(; kwargs...)
        obj = new(meta(ConceptMapGroup), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct ConceptMapGroup
const __meta_ConceptMapGroup = Ref{ProtoMeta}()
function meta(::Type{ConceptMapGroup})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ConceptMapGroup)
            __meta_ConceptMapGroup[] = target = ProtoMeta(ConceptMapGroup)
            allflds = Pair{Symbol,Union{Type,String}}[:owner => Concept, :concept_maps => Base.Vector{ConceptMap}]
            meta(target, ConceptMapGroup, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ConceptMapGroup[]
    end
end
function Base.getproperty(obj::ConceptMapGroup, name::Symbol)
    if name === :owner
        return (obj.__protobuf_jl_internal_values[name])::Concept
    elseif name === :concept_maps
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{ConceptMap}
    else
        getfield(obj, name)
    end
end

mutable struct Numeric <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Numeric(; kwargs...)
        obj = new(meta(Numeric), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Numeric
const __meta_Numeric = Ref{ProtoMeta}()
function meta(::Type{Numeric})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Numeric)
            __meta_Numeric[] = target = ProtoMeta(Numeric)
            allflds = Pair{Symbol,Union{Type,String}}[:long_value => Int64, :double_value => Float64, :nan => Bool]
            oneofs = Int[1,1,1]
            oneof_names = Symbol[Symbol("value")]
            meta(target, Numeric, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_Numeric[]
    end
end
function Base.getproperty(obj::Numeric, name::Symbol)
    if name === :long_value
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :double_value
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :nan
        return (obj.__protobuf_jl_internal_values[name])::Bool
    else
        getfield(obj, name)
    end
end

mutable struct NumericGroup <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function NumericGroup(; kwargs...)
        obj = new(meta(NumericGroup), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct NumericGroup
const __meta_NumericGroup = Ref{ProtoMeta}()
function meta(::Type{NumericGroup})
    ProtoBuf.metalock() do
        if !isassigned(__meta_NumericGroup)
            __meta_NumericGroup[] = target = ProtoMeta(NumericGroup)
            allflds = Pair{Symbol,Union{Type,String}}[:owner => Concept, :number => Numeric]
            meta(target, NumericGroup, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_NumericGroup[]
    end
end
function Base.getproperty(obj::NumericGroup, name::Symbol)
    if name === :owner
        return (obj.__protobuf_jl_internal_values[name])::Concept
    elseif name === :number
        return (obj.__protobuf_jl_internal_values[name])::Numeric
    else
        getfield(obj, name)
    end
end

export ConceptMap_MapEntry, ConceptMap, Explainables_RelationsEntry, Explainables_AttributesEntry, Explainables_OwnershipsEntry, Explainables_Owned_OwnedEntry, Explainables_Owned, Explainables, Explainable, ConceptMapGroup, Numeric, NumericGroup
# mapentries: "Explainables_Owned_OwnedEntry" => ("AbstractString", "Explainable"), "Explainables_RelationsEntry" => ("AbstractString", "Explainable"), "Explainables_OwnershipsEntry" => ("AbstractString", "Explainables_Owned"), "Explainables_AttributesEntry" => ("AbstractString", "Explainable"), "ConceptMap_MapEntry" => ("AbstractString", "Concept")


text:
using Autumn
include("scene.jl")
include("construct_observation_data.jl")

"""
USAGE: particles

a = parseautumn(" particles program from interface ")
mod = eval(compiletojulia(a))
observations, _ = generate_observations_particles(mod)
images = convert_render_output_to_image_seq(observations)

For other models, replace particles program string and replace generate_observations_particles
with one of the other generate_observations_[model name] in "construct_observation_data."

"""

function convert_render_output_to_image_seq(render_output, grid_size = 16)
  map(r -> convert_render_output_to_image(r, grid_size), render_output)
end

function convert_render_output_to_image(render_output, grid_size = 16)
  image = [RGBA(1.0, 0.0, 0.0, 1.0) for x in 1:grid_size, y in 1:grid_size]

  for cell in render_output
    color = rgb(cell.color)
    image[cell.position.y, cell.position.x] = RGBA(color.r, color.g, color.b, 1.0) # 0.6
  end
  image
end

# function particles_images(mod::Module)
#   observations, _ = generate_observations_particles(mod)
#   images = convert_render_output_to_image_seq(observations)
# end

# function ice_images(mod::Module)
#   observations, _ = generate_observations_ice(mod)
#   images = convert_render_output_to_image_seq(observations)
# end

# function ants_images(mod::Module)
#   a = au"""(program
#   (= GRID_SIZE 16)

#   (object Ant (Cell 0 0 "gray"))
#   (object Food (Cell 0 0 "red"))

#   (: ants (List Ant))
#   (= ants (initnext (map Ant (randomPositions GRID_SIZE 6))
#                     (updateObj (prev ants) (--> obj (move obj (unitVector obj (closest obj Food)))))))

#   (: foods (List Food))
#   (= foods (initnext (list)
#                      (updateObj (prev foods) (--> obj (if (intersects obj ants)
#                                                        then (removeObj obj)
#                                                        else obj)))))

#   (on clicked (= foods (addObj (prev foods) (map Food (randomPositions GRID_SIZE 4)))))
# )"""

#   global mod = eval(compiletojulia(a))
#   observations, _ = generate_observations_ants(mod)
#   images = convert_render_output_to_image_seq(observations)
# end

# function space_invaders_images(mod::Module)
#   a = au"""(program
#   (= GRID_SIZE 16)

#   (object Enemy (Cell 0 0 "blue"))
#   (object Hero (Cell 0 0 "black"))
#   (object Bullet (Cell 0 0 "red"))
#   (object EnemyBullet (Cell 0 0 "orange"))

#   (: enemies1 (List Enemy))
#   (= enemies1 (initnext (mapPositions
#                           Enemy
#                           GRID_SIZE
#                           (--> pos (& (== (.. pos y) 1) (== (% (.. pos x) 2) 0
#                                                             ))))
#                         (prev enemies1)))

#   (: enemies2 (List Enemy))
#   (= enemies2 (initnext (mapPositions
#                           Enemy
#                           GRID_SIZE
#                           (--> pos (& (== (.. pos y) 3) (== (% (.. pos x) 2) 1
#                                                             ))))
#                         (prev enemies2)))

#   (: hero Hero)
#   (= hero (initnext (Hero (Position 8 15)) (prev hero)))

#   (: enemyBullets (List EnemyBullet))
#   (= enemyBullets (initnext (list) (updateObj (prev enemyBullets) (--> obj (move obj 0 1)))))

#   (: bullets (List Bullet))
#   (= bullets (initnext (list) (updateObj (prev bullets) (--> obj (move obj 0 -1)))))

#   (: time Int)
#   (= time (initnext 0 (+ (prev time) 1)))

#   (on left (= hero (moveLeftNoCollision (prev hero))))
#   (on right (= hero (moveRightNoCollision (prev hero))))
#   (on (& up (.. (prev hero) alive)) (= bullets (addObj (prev bullets) (Bullet (.. (prev hero) origin)))))

#   (on (== (% time 10) 5) (= enemies1 (updateObj (prev enemies1) (--> obj (moveLeft obj)))))
#   (on (== (% time 10) 0) (= enemies1 (updateObj (prev enemies1) (--> obj (moveRight obj)))))

#   (on (== (% time 10) 5) (= enemies2 (updateObj (prev enemies2) (--> obj (moveRight obj)))))
#   (on (== (% time 10) 0) (= enemies2 (updateObj (prev enemies2) (--> obj (moveLeft obj)))))

#   (on (intersects (prev bullets) (prev enemies1))
#   	(let ((= bullets (removeObj (prev bullets) (--> obj (intersects obj (prev enemies1)))))
#           (= enemies1 (removeObj (prev enemies1) (--> obj (intersects obj (prev bullets)))))))
#   )

#   (on (intersects (prev bullets) (prev enemies2))
#   	(let ((= bullets (removeObj (prev bullets) (--> obj (intersects obj (prev enemies2)))))
#           (= enemies2 (removeObj (prev enemies2) (--> obj (intersects obj (prev bullets)))))))
#   )

#   (on (== (% time 5) 2) (= enemyBullets (addObj (prev enemyBullets) (EnemyBullet (uniformChoice (map (--> obj (.. obj origin)) (prev enemies2)))))))
#   (on (intersects (prev hero) (prev enemyBullets)) (= hero (removeObj (prev hero))))

#   (on (intersects (prev bullets) (prev enemyBullets))
#     (let
#        ((= bullets (removeObj (prev bullets) (--> obj (intersects obj (prev enemyBullets)))))
#         (= enemyBullets (removeObj (prev enemyBullets) (--> obj (intersects obj (prev bullets))))))))
# )"""

#   global mod = eval(compiletojulia(a))
#   observations, _ = generate_observations_space_invaders(mod)
#   images = convert_render_output_to_image_seq(observations)
# end

# function disease_images(mod::Module)
#   a = """(program
#   (= GRID_SIZE 16)

#   (object Particle (: health Bool) (Cell 0 0 (if health then "gray" else "darkgreen")))

#   (: inactiveParticles (List Particle))
#   (= inactiveParticles (initnext (list (Particle true (Position 5 3)) (Particle true (Position 11 14)) (Particle true (Position 4 8)) (Particle true (Position 9 9)) (Particle true (Position 1 13)) (Particle true (Position 14 5)))
#                        (updateObj (prev inactiveParticles) (--> obj (if (! (.. activeParticle health))
#                                                                      then (updateObj obj "health" false)
#                                                                      else obj))
# 														   (--> obj (adjacent (.. obj origin) (.. (prev activeParticle) origin))))))

#   (: activeParticle Particle)
#   (= activeParticle (initnext (Particle false (Position 0 0)) (prev activeParticle)))

#   (on (!= (length (filter (--> obj (! (.. obj health))) (adjacentObjs activeParticle))) 0) (= activeParticle (updateObj (prev activeParticle) "health" false)))
#   (on (clicked (prev inactiveParticles))
#       (let ((= inactiveParticles (addObj (prev inactiveParticles) (prev activeParticle)))
#             (= activeParticle (objClicked click (prev inactiveParticles)))
#             (= inactiveParticles (removeObj inactiveParticles (objClicked click (prev inactiveParticles))))
#            )))
#   (on left (= activeParticle (moveNoCollision (prev activeParticle) -1 0)))
#   (on right (= activeParticle (moveNoCollision (prev activeParticle) 1 0)))
#   (on up (= activeParticle (moveNoCollision (prev activeParticle) 0 -1)))
#   (on down (= activeParticle (moveNoCollision (prev activeParticle) 0 1)))
# )"""
#   global mod = eval(compiletojulia(a))
#   observations, _ = generate_observations_disease(mod)
#   images = convert_render_output_to_image_seq(observations)
# end

# function water_plug_images(mod::Module)
#   a = """(program
#   (= GRID_SIZE 16)

#   (object Button (: color String) (Cell 0 0 color))
#   (object Vessel (Cell 0 0 "purple"))
#   (object Plug (Cell 0 0 "orange"))
#   (object Water (Cell 0 0 "blue"))

#   (: vesselButton Button)
#   (= vesselButton (Button "purple" (Position 2 0)))
#   (: plugButton Button)
#   (= plugButton (Button "orange" (Position 5 0)))
#   (: waterButton Button)
#   (= waterButton (Button "blue" (Position 8 0)))
#   (: removeButton Button)
#   (= removeButton (Button "black" (Position 11 0)))
#   (: clearButton Button)
#   (= clearButton (Button "red" (Position 14 0)))

#   (: vessels (List Vessel))
#   (= vessels (initnext (list (Vessel (Position 6 15)) (Vessel (Position 6 14)) (Vessel (Position 6 13)) (Vessel (Position 5 12)) (Vessel (Position 4 11)) (Vessel (Position 3 10)) (Vessel (Position 9 15)) (Vessel (Position 9 14)) (Vessel (Position 9 13)) (Vessel (Position 10 12)) (Vessel (Position 11 11)) (Vessel (Position 12 10))) (prev vessels)))
#   (: plugs (List Plug))
#   (= plugs (initnext (list (Plug (Position 7 15)) (Plug (Position 8 15)) (Plug (Position 7 14)) (Plug (Position 8 14)) (Plug (Position 7 13)) (Plug (Position 8 13))) (prev plugs)))
#   (: water (List Water))
#   (= water (initnext (list) (updateObj (prev water) nextLiquid)))

#   (= currentParticle (initnext "vessel" (prev currentParticle)))

#   (on (& clicked (& (isFree click) (== currentParticle "vessel"))) (= vessels (addObj (prev vessels) (Vessel (Position (.. click x) (.. click y))))))
#   (on (& clicked (& (isFree click) (== currentParticle "plug"))) (= plugs (addObj (prev plugs) (Plug (Position (.. click x) (.. click y))))))
#   (on (& clicked (& (isFree click) (== currentParticle "water"))) (= water (addObj (prev water) (Water (Position (.. click x) (.. click y))))))
#   (on (clicked vesselButton) (= currentParticle "vessel"))
#   (on (clicked plugButton) (= currentParticle "plug"))
#   (on (clicked waterButton) (= currentParticle "water"))
#   (on (clicked removeButton) (= plugs (removeObj plugs (--> obj true))))
#   (on (clicked clearButton) (let ((= vessels (removeObj vessels (--> obj true))) (= plugs (removeObj plugs (--> obj true))) (= water (removeObj water (--> obj true))))))
# )"""
#   global mod = eval(compiletojulia(a))
#   observations, _ = generate_observations_water_plug(mod)
#   images = convert_render_output_to_image_seq(observations)
# end

# function paint_images(mod::Module)
#   a = """(program
#   (= GRID_SIZE 16)

#   (object Particle (: color String) (Cell 0 0 color))

#   (: particles (List Particle))
#   (= particles (initnext (list) (prev particles)))

#   (: currColor String)
#   (= currColor (initnext "red" (prev currColor)))

#   (on (& clicked (isFree click)) (= particles (addObj (prev particles) (Particle currColor (Position (.. click x) (.. click y))))))
#   (on (& up (== (prev currColor) "red")) (= currColor "gold"))
#   (on (& up (== (prev currColor) "gold")) (= currColor "green"))
#   (on (& up (== (prev currColor) "green")) (= currColor "blue"))
#   (on (& up (== (prev currColor) "blue")) (= currColor "purple"))
#   (on (& up (== (prev currColor) "purple")) (= currColor "red"))
# )"""
#   global mod = eval(compiletojulia(a))
#   observations, _ = generate_observations_paint(mod)
#   images = convert_render_output_to_image_seq(observations)
# end

# function gravity_images(mod::Module)
#   a = """(program
#   (= GRID_SIZE 16)
#   (= background "black")

#   (object Button (: color String) (Cell 0 0 color))
#   (object Blob (list (Cell 0 0 "blue") (Cell 0 1 "blue") (Cell 1 0 "blue") (Cell 1 1 "blue")))

#   (: leftButton Button)
#   (= leftButton (initnext (Button "red" (Position 0 7)) (prev leftButton)))

#   (: rightButton Button)
#   (= rightButton (initnext (Button "darkorange" (Position 15 7)) (prev rightButton)))

#   (: upButton Button)
#   (= upButton (initnext (Button "gold" (Position 7 0)) (prev upButton)))

#   (: downButton Button)
#   (= downButton (initnext (Button "green" (Position 7 15)) (prev downButton)))

#   (: blobs (List Blob))
#   (= blobs (initnext (list) (prev blobs)))

#   (: gravity String)
#   (= gravity (initnext "down" (prev gravity)))

#   (on (== gravity "left") (= blobs (updateObj (prev blobs) (--> obj (moveLeftNoCollision obj)))))
#   (on (== gravity "right") (= blobs (updateObj (prev blobs) (--> obj (moveRightNoCollision obj)))))
#   (on (== gravity "up") (= blobs (updateObj (prev blobs) (--> obj (moveUpNoCollision obj)))))
#   (on (== gravity "down") (= blobs (updateObj (prev blobs) (--> obj (moveDownNoCollision obj)))))

#   (on (& clicked (isFree click)) (= blobs (addObj (prev blobs) (Blob (Position (.. click x) (.. click y))))) )

#   (on (clicked leftButton) (= gravity "left"))

#   (on (clicked rightButton) (= gravity "right"))

#   (on (clicked upButton) (= gravity "up"))

#   (on (clicked downButton) (= gravity "down"))
# )"""
#   global mod = eval(compiletojulia(a))
#   observations, _ = generate_observations_gravity(mod)
#   images = convert_render_output_to_image_seq(observations)
# end

# function sand_images(mod::Module)
#   a = """(program
#   (= GRID_SIZE 16)

#   (object Button (: color String) (Cell 0 0 color))
#   (object Sand (: liquid Bool) (Cell 0 0 (if liquid then "sandybrown" else "tan")))
#   (object Water (Cell 0 0 "skyblue"))

#   (: sandButton Button)
#   (= sandButton (initnext (Button "tan" (Position 2 0)) (prev sandButton)))

#   (: waterButton Button)
#   (= waterButton (initnext (Button "skyblue" (Position 7 0)) (prev waterButton)))

#   (: sand (List Sand))
#   (= sand (initnext (list)
#             (updateObj (prev sand) (--> obj (if (.. obj liquid)
# 											 then (nextLiquid obj)
# 											 else (nextSolid obj))))))

#   (: water (List Water))
#   (= water (initnext (list) (updateObj (prev water) (--> obj (nextLiquid obj)))))

#   (: clickType String)
#   (= clickType (initnext "sand" (prev clickType)))

#   (on true (= sand (updateObj (prev sand) (--> obj (updateObj obj "liquid" true)) (--> obj (& (! (.. obj liquid)) (intersects (adjacentObjs obj) (prev water)))))))

#   (on (clicked sandButton) (= clickType "sand"))
#   (on (clicked waterButton) (= clickType "water"))
#   (on (& (& clicked (isFree click)) (== clickType "sand"))  (= sand (addObj sand (Sand false (Position (.. click x) (.. click y))))))
#   (on (& (& clicked (isFree click)) (== clickType "water")) (= water (addObj water (Water (Position (.. click x) (.. click y))))))

# )"""
#   global mod = eval(compiletojulia(a))
#   observations, _ = generate_observations_sand_simple(mod)
#   images = convert_render_output_to_image_seq(observations, 16)
# end


text:
module GPs

# https://gregorygundersen.com/blog/2019/09/12/practical-gp-regression/
# NOTE: If A * x = b. Then A \ b solves for x.
# NOTE: cholesky(A) \ b = L' \ (L \ b) = x, where L = cholesky(A).L
# NOTE: A \ b and cholesky(A) \ b are respectively implemented via variants of
# `solve` and `cholesky solve` in other languages.

using Distributions
using Distances

include(joinpath(@__DIR__, "Kernels.jl"))
using .Kernels

import LinearAlgebra
import LinearAlgebra: cholesky, Symmetric

export GP, posterior, Kernels

struct GP
  mu
  kernel
  y
  X
  sigma
  C
  d
  a
end

function GP(y, X, kernel; mu=0, sigma=0)
  Ddata = pairwise(metric(kernel), X, dims=1)
  Kdata = kernel(Ddata)
  C = cholesky(Symmetric(Kdata))
  d = y .- mu
  a = C \ d  # = L' \ (L \ d), where L = cholesky(Kdata).L
  return GP(mu, kernel, y, X, sigma, C, d, a)
end

function posterior(gp, Xnew)
  Dnew = pairwise(metric(gp.kernel), Xnew, dims=1)
  D_new_data = pairwise(metric(gp.kernel), Xnew, gp.X, dims=1)

  Knew = gp.kernel(Dnew)
  K_new_data = gp.kernel(D_new_data, jitter=false)

  mu = gp.mu .+ K_new_data * gp.a
  S = Knew - K_new_data * (gp.C \ K_new_data')

  return MvNormal(mu, Symmetric(S) + LinearAlgebra.I * gp.sigma)
end

quantiles(X, q; dims) = mapslices(v -> quantile(v, q), X, dims=dims)

function ci(preds, alpha=0.05, dims=2)
    q_upper = 1 - alpha / 2
    q_lower = alpha / 2
    lower = vec(quantiles(preds, q_lower, dims=dims))
    upper = vec(quantiles(preds, q_upper, dims=dims))
    return lower, upper
end

end

max_stars_repo_name: luiarthur/TuringBnpBenchmarks

max_issues_repo_name: luiarthur/TuringBnpBenchmarks

max_forks_repo_name: luiarthur/TuringBnpBenchmarks


text:
"""
    mixerblock(planes, npatches; mlp_ratio = (0.5, 4.0), mlp_layer = mlp_block,
               dropout = 0., drop_path_rate = 0., activation = gelu)

Creates a feedforward block for the MLPMixer architecture.
([reference](https://arxiv.org/pdf/2105.01601))

# Arguments:

  - `planes`: the number of planes in the block
  - `npatches`: the number of patches of the input
  - `mlp_ratio`: number(s) that determine(s) the number of hidden channels in the token mixing MLP
    and/or the channel mixing MLP as a ratio to the number of planes in the block.
  - `mlp_layer`: the MLP layer to use in the block
  - `dropout`: the dropout rate to use in the MLP blocks
  - `drop_path_rate`: Stochastic depth rate
  - `activation`: the activation function to use in the MLP blocks
"""
function mixerblock(planes, npatches; mlp_ratio = (0.5, 4.0), mlp_layer = mlp_block,
                    dropout = 0.0, drop_path_rate = 0.0, activation = gelu)
    tokenplanes, channelplanes = [Int(r * planes) for r in mlp_ratio]
    return Chain(SkipConnection(Chain(LayerNorm(planes),
                                      swapdims((2, 1, 3)),
                                      mlp_layer(npatches, tokenplanes; activation, dropout),
                                      swapdims((2, 1, 3)),
                                      DropPath(drop_path_rate)), +),
                 SkipConnection(Chain(LayerNorm(planes),
                                      mlp_layer(planes, channelplanes; activation, dropout),
                                      DropPath(drop_path_rate)), +))
end

"""
    mlpmixer(block, imsize::Dims{2} = (224, 224); inchannels = 3, norm_layer = LayerNorm,
             patch_size::Dims{2} = (16, 16), embedplanes = 512, drop_path_rate = 0.,
             depth = 12, nclasses = 1000, kwargs...)

Creates a model with the MLPMixer architecture.
([reference](https://arxiv.org/pdf/2105.01601)).

# Arguments

  - `block`: the type of mixer block to use in the model - architecture dependent
    (a constructor of the form `block(embedplanes, npatches; drop_path_rate, kwargs...)`)
  - `imsize`: the size of the input image
  - `inchannels`: the number of input channels
  - `norm_layer`: the normalization layer to use in the model
  - `patch_size`: the size of the patches
  - `embedplanes`: the number of channels after the patch embedding (denotes the hidden dimension)
  - `drop_path_rate`: Stochastic depth rate
  - `depth`: the number of blocks in the model
  - `nclasses`: number of output classes
  - `kwargs`: additional arguments (if any) to pass to the mixer block. Will use the defaults if
    not specified.
"""
function mlpmixer(block, imsize::Dims{2} = (224, 224); inchannels = 3,
                  norm_layer = LayerNorm,
                  patch_size::Dims{2} = (16, 16), embedplanes = 512, drop_path_rate = 0.0,
                  depth = 12, nclasses = 1000, kwargs...)
    npatches = prod(imsize .÷ patch_size)
    dp_rates = LinRange{Float32}(0.0, drop_path_rate, depth)
    layers = Chain(PatchEmbedding(imsize; inchannels, patch_size, embedplanes),
                   Chain([block(embedplanes, npatches; drop_path_rate = dp_rates[i],
                                kwargs...)
                          for i in 1:depth]))

    classification_head = Chain(norm_layer(embedplanes), seconddimmean,
                                Dense(embedplanes, nclasses))
    return Chain(layers, classification_head)
end

# Configurations for MLPMixer models
mixer_configs = Dict(:small => Dict(:depth => 8, :planes => 512),
                     :base => Dict(:depth => 12, :planes => 768),
                     :large => Dict(:depth => 24, :planes => 1024),
                     :huge => Dict(:depth => 32, :planes => 1280))

struct MLPMixer
    layers::Any
end

"""
    MLPMixer(size::Symbol = :base; patch_size::Dims{2} = (16, 16),
             imsize::Dims{2} = (224, 224), drop_path_rate = 0., nclasses = 1000)

Creates a model with the MLPMixer architecture.
([reference](https://arxiv.org/pdf/2105.01601)).

# Arguments

  - `size`: the size of the model - one of `small`, `base`, `large` or `huge`
  - `patch_size`: the size of the patches
  - `imsize`: the size of the input image
  - `drop_path_rate`: Stochastic depth rate
  - `nclasses`: number of output classes

See also [`Metalhead.mlpmixer`](#).
"""
function MLPMixer(size::Symbol = :base; patch_size::Dims{2} = (16, 16),
                  imsize::Dims{2} = (224, 224), drop_path_rate = 0.0, nclasses = 1000)
    @assert size in keys(mixer_configs) "`size` must be one of $(keys(mixer_configs))"
    depth = mixer_configs[size][:depth]
    embedplanes = mixer_configs[size][:planes]
    layers = mlpmixer(mixerblock, imsize; patch_size, embedplanes, depth, drop_path_rate,
                      nclasses)
    return MLPMixer(layers)
end

@functor MLPMixer

(m::MLPMixer)(x) = m.layers(x)

backbone(m::MLPMixer) = m.layers[1]
classifier(m::MLPMixer) = m.layers[2]

"""
    resmixerblock(planes, npatches; dropout = 0., drop_path_rate = 0., mlp_ratio = 4.0,
                  activation = gelu, λ = 1e-4)

Creates a block for the ResMixer architecture.
([reference](https://arxiv.org/abs/2105.03404)).

# Arguments

  - `planes`: the number of planes in the block
  - `npatches`: the number of patches of the input
  - `mlp_ratio`: ratio of the number of hidden channels in the channel mixing MLP to the number
    of planes in the block
  - `mlp_layer`: the MLP block to use
  - `dropout`: the dropout rate to use in the MLP blocks
  - `drop_path_rate`: Stochastic depth rate
  - `activation`: the activation function to use in the MLP blocks
  - `λ`: initialisation constant for the LayerScale
"""
function resmixerblock(planes, npatches; mlp_ratio = 4.0, mlp_layer = mlp_block,
                       dropout = 0.0, drop_path_rate = 0.0, activation = gelu, λ = 1e-4)
    return Chain(SkipConnection(Chain(Flux.Scale(planes),
                                      swapdims((2, 1, 3)),
                                      Dense(npatches, npatches),
                                      swapdims((2, 1, 3)),
                                      LayerScale(planes, λ),
                                      DropPath(drop_path_rate)), +),
                 SkipConnection(Chain(Flux.Scale(planes),
                                      mlp_layer(planes, Int(mlp_ratio * planes); dropout,
                                                activation),
                                      LayerScale(planes, λ),
                                      DropPath(drop_path_rate)), +))
end

struct ResMLP
    layers::Any
end

"""
    ResMLP(size::Symbol = :base; patch_size::Dims{2} = (16, 16), imsize::Dims{2} = (224, 224),
           drop_path_rate = 0., nclasses = 1000)

Creates a model with the ResMLP architecture.
([reference](https://arxiv.org/abs/2105.03404)).

# Arguments

  - `size`: the size of the model - one of `small`, `base`, `large` or `huge`
  - `patch_size`: the size of the patches
  - `imsize`: the size of the input image
  - `drop_path_rate`: Stochastic depth rate
  - `nclasses`: number of output classes

See also [`Metalhead.mlpmixer`](#).
"""
function ResMLP(size::Symbol = :base; patch_size::Dims{2} = (16, 16),
                imsize::Dims{2} = (224, 224), drop_path_rate = 0.0, nclasses = 1000)
    @assert size in keys(mixer_configs) "`size` must be one of $(keys(mixer_configs))"
    depth = mixer_configs[size][:depth]
    embedplanes = mixer_configs[size][:planes]
    layers = mlpmixer(resmixerblock, imsize; mlp_ratio = 4.0, patch_size, embedplanes,
                      drop_path_rate, depth, nclasses)
    return ResMLP(layers)
end

@functor ResMLP

(m::ResMLP)(x) = m.layers(x)

backbone(m::ResMLP) = m.layers[1]
classifier(m::ResMLP) = m.layers[2]

"""
    SpatialGatingUnit(norm, proj)

Creates a spatial gating unit as described in the gMLP paper.
([reference](https://arxiv.org/abs/2105.08050))

# Arguments

  - `norm`: the normalisation layer to use
  - `proj`: the projection layer to use
"""
struct SpatialGatingUnit{T, F}
    norm::T
    proj::F
end

"""
    SpatialGatingUnit(planes::Integer, npatches::Integer; norm_layer = LayerNorm)

Creates a spatial gating unit as described in the gMLP paper.
([reference](https://arxiv.org/abs/2105.08050))

# Arguments

  - `planes`: the number of planes in the block
  - `npatches`: the number of patches of the input
  - `norm_layer`: the normalisation layer to use
"""
function SpatialGatingUnit(planes::Integer, npatches::Integer; norm_layer = LayerNorm)
    gateplanes = planes ÷ 2
    norm = norm_layer(gateplanes)
    proj = Dense(2 * eps(Float32) .* rand(Float32, npatches, npatches), ones(npatches))
    return SpatialGatingUnit(norm, proj)
end

@functor SpatialGatingUnit

function (m::SpatialGatingUnit)(x)
    u, v = chunk(x, 2; dims = 1)
    v = m.norm(v)
    v = m.proj(permutedims(v, (2, 1, 3)))
    return u .* permutedims(v, (2, 1, 3))
end

"""
    spatial_gating_block(planes, npatches; mlp_ratio = 4.0, mlp_layer = gated_mlp_block,
                         norm_layer = LayerNorm, dropout = 0.0, drop_path_rate = 0.,
                         activation = gelu)

Creates a feedforward block based on the gMLP model architecture described in the paper.
([reference](https://arxiv.org/abs/2105.08050))

# Arguments

  - `planes`: the number of planes in the block
  - `npatches`: the number of patches of the input
  - `mlp_ratio`: ratio of the number of hidden channels in the channel mixing MLP to the number
    of planes in the block
  - `norm_layer`: the normalisation layer to use
  - `dropout`: the dropout rate to use in the MLP blocks
  - `drop_path_rate`: Stochastic depth rate
  - `activation`: the activation function to use in the MLP blocks
"""
function spatial_gating_block(planes, npatches; mlp_ratio = 4.0, norm_layer = LayerNorm,
                              mlp_layer = gated_mlp_block, dropout = 0.0,
                              drop_path_rate = 0.0,
                              activation = gelu)
    channelplanes = Int(mlp_ratio * planes)
    sgu = inplanes -> SpatialGatingUnit(inplanes, npatches; norm_layer)
    return SkipConnection(Chain(norm_layer(planes),
                                mlp_layer(sgu, planes, channelplanes; activation, dropout),
                                DropPath(drop_path_rate)), +)
end

struct gMLP
    layers::Any
end

"""
    gMLP(size::Symbol = :base; patch_size::Dims{2} = (16, 16),
         imsize::Dims{2} = (224, 224), drop_path_rate = 0., nclasses = 1000)

Creates a model with the gMLP architecture.
([reference](https://arxiv.org/abs/2105.08050)).

# Arguments

  - `size`: the size of the model - one of `small`, `base`, `large` or `huge`
  - `patch_size`: the size of the patches
  - `imsize`: the size of the input image
  - `drop_path_rate`: Stochastic depth rate
  - `nclasses`: number of output classes

See also [`Metalhead.mlpmixer`](#).
"""
function gMLP(size::Symbol = :base; patch_size::Dims{2} = (16, 16),
              imsize::Dims{2} = (224, 224), drop_path_rate = 0.0, nclasses = 1000)
    @assert size in keys(mixer_configs) "`size` must be one of $(keys(mixer_configs))"
    depth = mixer_configs[size][:depth]
    embedplanes = mixer_configs[size][:planes]
    layers = mlpmixer(spatial_gating_block, imsize; mlp_layer = gated_mlp_block,
                      patch_size, embedplanes, drop_path_rate, depth, nclasses)
    return gMLP(layers)
end

@functor gMLP

(m::gMLP)(x) = m.layers(x)

backbone(m::gMLP) = m.layers[1]
classifier(m::gMLP) = m.layers[2]


text:
function _mergesort!(a::Vector, tmp::Vector, leftindex, rightindex, isless)
    if (rightindex + 1 - leftindex)  <= 9
        insertionsort!(view(a, leftindex:rightindex), isless)
        return
    end
    midindex = (leftindex + rightindex) ÷ 2
    _mergesort!(a, tmp, leftindex, midindex, isless)
    _mergesort!(a, tmp, midindex+1, rightindex, isless)

    for i = midindex:-1:leftindex
        tmp[i] = a[i]
    end
    for j = 1:rightindex-midindex
        tmp[rightindex-j+1] = a[midindex+j]
    end
    let i = k = leftindex, j = rightindex
        while k <= rightindex
            if isless(tmp[i], tmp[j])
                a[k] = tmp[i]
                i += 1
            else
                a[k] = tmp[j]
                j -= 1
            end
            k += 1
        end
    end
end

function mergesort!(a::Vector, isless=isless)
    _mergesort!(a, similar(a), 1, length(a), isless)
    a
end


text:
import Base: copy, eltype, iterate, IteratorSize
import Base: getindex, setindex!, delete!

export BBSearch, BreadthFirstBBSearch, DepthFirstBBSearch
export data
export copy, eltype, iterate, IteratorSize

abstract type AbstractBBNode end

"""
    BBNode <: AbstractBBNode

Intermediate node of a `BBTree`. Does not contain any data by itself,
only redirect toward its children.
"""
struct BBNode <: AbstractBBNode
    parent::Int
    children::Vector{Int}
end

"""
    BBLeaf{DATA} <: AbstractBBLeaf

Leaf node of a `BBTree` that contains some data. Its status is either
    - `:working`: the leaf will be further processed.
    - `:final`: the leaf won't be touched anymore.
"""
struct BBLeaf{DATA} <: AbstractBBNode
    data::DATA
    parent::Int
    status::Symbol
end

function BBLeaf(data::DATA, parent::Int) where {DATA}
    BBLeaf{DATA}(data, parent, :working)
end

function BBNode(leaf::BBLeaf, child1::Int, child2::Int)
    BBNode(leaf.parent, Int[child1, child2])
end

"""
    BBTree{DATA}

Tree storing the data used and produced by a branch and bound search in a
structured way.

Nodes and leaves can be accessed using their index using the bracket syntax
`wt[node_id]`. However this is slow, as nodes and leaves are stored separately.

Support the iterator interface. The element yielded by the iteration are
tuples `(node_id, lvl)` where `lvl` is the depth of the node in the tree.
"""
struct BBTree{DATA}
    nodes::Dict{Int, BBNode}
    leaves::Dict{Int, BBLeaf{DATA}}
    working_leaves::Vector{Int}
end

function BBTree(rootdata::DATA) where {DATA}
    rootleaf = BBLeaf(rootdata, 0)
    BBTree{DATA}(Dict{Int, BBNode}(), Dict(1 => rootleaf), Int[1])
end

show(io::IO, wn::BBNode) = print(io, "Node with children $(wn.children)")

function show(io::IO, wl::BBLeaf)
    print(io, "Leaf (:$(wl.status)) with data $(wl.data)")
end

function show(io::IO, wt::BBTree{DATA}) where {DATA}
    println(io, "Working tree with $(nnodes(wt)) elements of type $DATA")

    if nnodes(wt) > 0
        println(io, "Indices: ", vcat(collect(keys(wt.nodes)), collect(keys(wt.leaves))) |> sort)
        println(io, "Structure:")
        for (id, lvl) in wt
            println(io, "  "^lvl * "[$id] $(wt[id])")
        end
    end
end

# Root node has id 1 and parent id 0
root(wt::BBTree) = wt[1]
is_root(wt::BBTree, id::Int) = (id == 1)

"""
    nnodes(wt::BBTree)

Number of nodes (including leaves) in a `BBTree`.
"""
nnodes(wt::BBTree) = length(wt.nodes) + length(wt.leaves)

"""
    data(leaf::BBLeaf)

Return the data stored in the leaf.
"""
data(leaf::BBLeaf) = leaf.data

"""
    data(wt::BBTree)

Return all the data stored in a `BBTree` as a list. The ordering of the elements
is arbitrary.
"""
data(wt::BBTree) = data.(values(wt.leaves))

function newid(wt::BBTree)
    k1 = keys(wt.nodes)
    k2 = keys(wt.leaves)

    if length(k1) > 0
        m1 = maximum(k1)
    else
        m1 = 0
    end

    if length(k2) > 0
        m2 = maximum(k2)
    else
        m2 = 0
    end

    return max(m1, m2) + 1
end

# Index operations (slower than manipulating the node directly in the correct
# dictionary)
function getindex(wt::BBTree, id)
    haskey(wt.nodes, id) && return wt.nodes[id]
    haskey(wt.leaves, id) && return wt.leaves[id]
    error("getindex failed: no index $id")  # TODO: make better error
end

setindex!(wt::BBTree, val::BBNode, id) = setindex!(wt.nodes, val, id)
setindex!(wt::BBTree, val::BBLeaf, id) = setindex!(wt.leaves, val, id)

function delete!(wt::BBTree, id)
    if haskey(wt.nodes, id)
        delete!(wt.nodes, id)
    elseif haskey(wt.leaves, id)
        delete!(wt.leaves, id)
    else
        error("delete! failed: no index $id")  # TODO: make better error
    end
end

"""
    discard_leaf!(wt::BBTree, id::Int)

Delete the `BBLeaf` with index `id` and all its ancestors to which it is
the last descendant.
"""
function discard_leaf!(wt::BBTree, id::Int)
    leaf = wt.leaves[id]
    delete!(wt.leaves, id)
    recursively_delete_parent!(wt, leaf.parent, id)
end

function recursively_delete_parent!(wt, id_parent, id_child)
    if !is_root(wt, id_child)
        parent = wt.nodes[id_parent]
        siblings = parent.children
        if length(parent.children) == 1  # The child has no siblings, so delete the parent
            delete!(wt.nodes, id_parent)
            recursively_delete_parent!(wt, parent.parent, id_parent)
        else  # The child has siblings so remove it from the children list
            deleteat!(parent.children, searchsortedfirst(parent.children, id_child))
        end
    end
end

function iterate(wt::BBTree, (id, lvl)=(0, 0))
    id, lvl = next_id(wt, id, lvl)
    lvl == 0 && return nothing
    return (id, lvl), (id, lvl)
end

function next_id(wt::BBTree, id, lvl)
    lvl == 0 && return (1, 1)
    node = wt[id]
    isa(node, BBNode) && return (node.children[1], lvl + 1)
    return next_sibling(wt, id, lvl)
end

function next_sibling(wt::BBTree, sibling, lvl)
    parent = wt[sibling].parent
    parent == 0 && return (0, 0)
    children = wt[parent].children
    maximum(children) == sibling && return next_sibling(wt, parent, lvl - 1)
    id = minimum(filter(x -> x > sibling, children))
    return (id, lvl)
end

"""
    BBSearch{DATA}

Branch and bound search interface in element of type DATA.

This interface provide an iterable that perform the search.

There is currently three types of search supported `BreadFirstBBSearch`,
`DepthFirstBBSearch` and `KeyBBSearch`, each one processing the element of the
tree in a different order. When subtyping one of these, the following methods
must be implemented:
    - `root_element(::BBSearch)`: return the element with which the search is started
    - `process(::BBSearch, elem::DATA)`: return a symbol representing the action
        to perform with the element `elem` and an object of type `DATA` reprensenting
        the state of the element after processing (may return `elem` unchanged).
    - `bisect(::BBSearch, elem::DATA)`: return two elements of type `DATA` build
        by bisecting `elem`

Subtyping `BBSearch` directly allows to have control over the order in which
the elements are process. To do this the following methods must be implemented:
    - `root_element(::BBSearch)`: return the first element to be processed. Use
        to build the initial tree.
    - `get_leaf_id!(::BBSearch, wt::BBTree)`: return the id of the next leaf that
        will be processed and remove it from the list of working leaves of `wt`.
    - `insert_leaf!(::BBSearch, wt::BBTree, leaf::BBLeaf)`: insert a leaf in the
        list of working leaves.

# Valid symbols returned by the process function
    - `:store`: the element is considered as final and is stored, it will not be
        further processed
    - `:bisect`: the element is bisected and each of the two resulting part will
        be processed
    - `:discard`: the element is discarded from the tree, allowing to free memory
"""
abstract type BBSearch{DATA} end

abstract type BreadthFirstBBSearch{DATA} <: BBSearch{DATA} end
abstract type DepthFirstBBSearch{DATA} <: BBSearch{DATA} end

"""
    KeyBBSearch{DATA} <: BBSearch{DATA}

Interface to a branch and bound search that use a key function to decide which
element to process first. The search process first the element with the largest
key as computed by `keyfunc(ks::KeyBBSearch, elem)`.

!!WARNING: Untested.
"""
abstract type KeyBBSearch{DATA} <: BBSearch{DATA} end

"""
    root_element(search::BBSearch)

Return the initial element of the search. The `BBTree` will be build around it.

Can be define for custom searches that are direct subtype of `BBSearch`, default
behavior is to fetch the field `initial` of the search.
"""
root_element(search::BBSearch) = search.initial

"""
    get_leaf_id!(::BBSearch, wt::BBTree)

Return the id of the next leaf that will be processed and remove it from the
list of working leaves.

Must be define for custom searches that are direct subtype of `BBSearch`.
"""
get_leaf_id!(::BreadthFirstBBSearch, wt::BBTree) = popfirst!(wt.working_leaves)
get_leaf_id!(::DepthFirstBBSearch, wt::BBTree) = pop!(wt.working_leaves)
get_leaf_id!(::KeyBBSearch, wt::BBTree) = popfirst!(wt.working_leaves)

"""
    insert_leaf!(::BBSearch, wt::BBTree, leaf::BBLeaf)

Insert the id of a new leaf that has been produced by bisecting an older leaf
into the list of working leaves.

Must be define for custom searches that are direct subtype of `BBSearch`.
"""
function insert_leaf!(::Union{BreadthFirstBBSearch{DATA}, DepthFirstBBSearch{DATA}},
                      wt::BBTree{DATA}, leaf::BBLeaf{DATA}) where {DATA}
    id = newid(wt)
    wt.leaves[id] = leaf
    push!(wt.working_leaves, id)
    return id
end

function insert_leaf!(::KS, wt::BBTree{DATA}, leaf::BBLeaf{DATA}) where {DATA, KS <: KeyBBSearch{DATA}}
    id = newid(wt)
    wt.leaves[id] = leaf
    keys = keyfunc.(KS, wt.working_leaves)
    current = keyfunc(KS, leaf)

    # Keep the working_leaves sorted
    insert!(wt.working_leaves, searchsortedfirst(keys, current), id)
    return id
end

eltype(::Type{BBS}) where {DATA, BBS <: BBSearch{DATA}} = BBTree{DATA}
IteratorSize(::Type{BBS}) where {BBS <: BBSearch} = Base.SizeUnknown()

function iterate(search::BBSearch{DATA},
                 wt::BBTree=BBTree(root_element(search))) where {DATA}

    isempty(wt.working_leaves) && return nothing

    id = get_leaf_id!(search, wt)
    X = wt.leaves[id]
    action, newdata = process(search, data(X))
    if action == :store
        wt.leaves[id] = BBLeaf(newdata, X.parent, :final)
    elseif action == :bisect
        child1, child2 = bisect(search, newdata)
        leaf1 = BBLeaf(child1, id, :working)
        leaf2 = BBLeaf(child2, id, :working)
        id1 = insert_leaf!(search, wt, leaf1)
        id2 = insert_leaf!(search, wt, leaf2)
        wt.nodes[id] = BBNode(X, id1, id2)
        delete!(wt.leaves, id)
    elseif action == :discard
        discard_leaf!(wt, id)
    else
        error("Branch and bound: process function of the search object return " *
              "unknown action: $action for element $X. Valid actions are " *
              ":store, :bisect and :discard.")
    end
    return wt, wt
end


text:
export RatApp, f0, f10, f20

import Base: +, -, *, /, \, inv, conj, sqrt
import Base: zero, one, oneunit, promote_rule

struct RatApp{T<:Number} <: Number
    a::T
    b::T
    c::T
    d::T
    RatApp(a::T, b::T, c::T, d::T) where T<:Number = new{T}(a, b, c, d)
    RatApp{T}(a::T, b::T, c::T, d::T) where T<:Number = new{T}(a, b, c, d)
end

function RatApp{T}(a::T...) where T
    RatApp(a...)
end
function RatApp{T}(a...) where T
    RatApp{T}(convert.(T, a)...)
end
function RatApp(a::Number, b::Number, c::Number, d::Number)
    RatApp(promote(a, b, c, d)...)
end
function RatApp(a::Number, b::Number)
    RatApp(promote(a, b, oftype(a, 1), oftype(a, 0))...)
end
function RatApp(a::Number)
    RatApp(promote(a, oftype(a, 0), oftype(a, 1), oftype(a, 0))...)
end
function convert(R::Type{RatApp{T}}, x::RatApp{S}) where {T,S}
    R(convert(T, x.a), convert(T, x.b), convert(T, x.c), convert(T, x.d))
end

function promote_rule(A::Type{RatApp{T}}, B::Type{RatApp{S}}) where {S,T}
    RatApp{promote_type(T, S)}
end
function promote_rule(A::Type{RatApp{T}}, B::Type{S}) where {S,T}
    RatApp{promote_type(T, S)}
end

zero(::Type{R}) where {T,R<:RatApp{T}} = R(zero(T), zero(T), oneunit(T), zero(T))
oneunit(::Type{R}) where {T,R<:RatApp{T}} = R(oneunit(T), zero(T), oneunit(T), zero(T))
one(::Type{R}) where R<:RatApp = oneunit(R)

-(x::T) where T<:RatApp = T(-x.a, -x.b, x.c, x.d)
-(x::T, y::T) where T<:RatApp = normalized(minus(x, y))
+(x::T, y::T) where T<:RatApp = normalized(minus(x, y))
inv(x::T) where T<:RatApp = T(x.c, x.d, x.a, x.b)
*(x::T, y::T) where T<:RatApp = normalized(times(x, y))
/(x::T, y::T) where T<:RatApp = normalized(divideby(x, y))
conj(x::RatApp{<:Real}) = x
conj(x::R) where R<:RatApp = R(conj(x.a), conj(x.b), conj(x.c), conj(x.d))

function f0(x::T) where T<:RatApp
    x.a / x.c
end
function f10(x::T) where T<:RatApp
    (x.b * x.c - x.a * x.d) / x.a * x.c
end
function f20(x::T) where T<:RatApp
    -2 * f10(x) * x.d / x.c
end

function plus(x::T, y::T) where T<:RatApp
    a, b, c, d = x.a, x.b, x.c, x.d
    A, B, C, D = y.a, y.b, y.c, y.d
    aa = a * C + c * A
    bb = a * D + b * C + c * B + d * A
    cc = c * C
    dd = c * D + d * C
    T(aa, bb, cc, dd)
end
function minus(x::T, y::T) where T<:RatApp
    a, b, c, d = x.a, x.b, x.c, x.d
    A, B, C, D = y.a, y.b, y.c, y.d
    aa = a * C - c * A
    bb = a * D + b * C - c * B - d * A
    cc = c * C
    dd = c * D + d * C
    T(aa, bb, cc, dd)
end

function times(x::T, y::T) where T<:RatApp
    a, b, c, d = x.a, x.b, x.c, x.d
    A, B, C, D = y.a, y.b, y.c, y.d
    aa = a * A
    bb = a * B + b * A
    cc = c * C
    dd = c * D + d * C
    T(aa, bb, cc, dd)
end
function divideby(x::T, y::T) where T<:RatApp
    a, b, c, d = x.a, x.b, x.c, x.d
    A, B, C, D = y.c, y.d, y.a, y.b
    aa = a * A
    bb = a * B + b * A
    cc = c * C
    dd = c * D + d * C
    T(aa, bb, cc, dd)
end

function sqrt(x::R) where R<:RatApp
    aa = sqrt(x.a)
    bb = x.b / 2aa
    cc = sqrt(x.c)
    dd = x.d / 2cc
    R(aa, bb, cc, dd)
end

@inline function normalized(x::T) where T<:RatApp
    a, b, c, d = x.a, x.b, x.c, x.d
    n1 = max(abs(a), abs(b))
    n2 = max(abs(c), abs(d))
    n = max(n1, n2)
    if iszero(n)
        x
    elseif iszero(n1)
        T(a, b, one(a), a)
    elseif iszero(n2)
        T(one(c), c, c, c)
    elseif a * d == b * c
        if iszero(a)
            if abs(d) >= abs(b)
                T(b / d, zero(a), one(a), zero(a))
            else
                T(one(a), zero(a), d / b, zero(a))
            end
        else
            if abs(c) >= abs(a)
                T(a / c, zero(a), one(a), zero(a))
            else
                T(one(a), zero(a), c / a, zero(a))
            end
        end
    elseif min(n1, n2) > 2^20
        n = oftype(n, 2^20)
        if real(c) < 0 || c == 0 && d < 0
            n = -n
        end
        T(a / n, b / n, c / n, d / n)
    else
        x
    end
end


text:
# A standard implementation of the Conjugate Gradient method.
# The only non-standard point about it is that it does not check
# that the operator is definite.
# It is possible to check that the system is inconsistent by
# monitoring ‖p‖, which would cost an extra norm computation per
# iteration.
#
# Dominique Orban, <dominique.orban@gerad.ca>
# Salt Lake City, UT, March 2015.

export cg

"""The conjugate gradient method to solve the symmetric linear system Ax=b.

The method does _not_ abort if A is not definite.

A preconditioner M may be provided in the form of a linear operator and is
assumed to be symmetric and positive definite.
"""
function cg(A :: AbstractLinearOperator, b :: AbstractVector{T};
            M :: AbstractLinearOperator=opEye(size(A,1)), atol :: Float64=1.0e-8,
            rtol :: Float64=1.0e-6, itmax :: Int=0, radius :: Float64=0.0,
            verbose :: Bool=false) where T <: Number

  n = size(b, 1);
  (size(A, 1) == n & size(A, 2) == n) || error("Inconsistent problem size");
  verbose && @printf("CG: system of %d equations in %d variables\n", n, n);

  # Initial state.
  x = zeros(T, n);
  r = copy(b)
  z = M * r
  p = copy(z)
  γ = @kdot(n, r, z)
  γ == 0 && return x, SimpleStats(true, false, [0.0], [], "x = 0 is a zero-residual solution")

  iter = 0;
  itmax == 0 && (itmax = 2 * n);

  rNorm = sqrt(γ);
  rNorms = [rNorm;];
  ε = atol + rtol * rNorm;
  verbose && @printf("%5d  %8.1e  ", iter, rNorm);

  solved = rNorm <= ε;
  tired = iter >= itmax;
  on_boundary = false;
  status = "unknown";

  while ! (solved || tired)
    Ap = A * p;
    pAp = @kdot(n, p, Ap)

    α = γ / pAp;

    # Compute step size to boundary if applicable.
    σ = radius > 0.0 ? maximum(to_boundary(x, p, radius)) : α

    verbose && @printf("%8.1e  %7.1e  %7.1e\n", pAp, α, σ);

    # Move along p from x to the boundary if either
    # the next step leads outside the trust region or
    # we have nonpositive curvature.
    if (radius > 0.0) & ((pAp <= 0.0) | (α > σ))
      α = σ
      on_boundary = true
    end

    @kaxpy!(n,  α,  p, x)
    @kaxpy!(n, -α, Ap, r)
    z = M * r
    γ_next = @kdot(n, r, z)
    rNorm = sqrt(γ_next);
    push!(rNorms, rNorm);

    solved = (rNorm <= ε) | on_boundary;
    if !solved
      β = γ_next / γ;
      γ = γ_next;

      @kscal!(n, β, p)
      @kaxpy!(n, 1.0, z, p)
    end

    iter = iter + 1;
    tired = iter >= itmax;
    verbose && @printf("%5d  %8.1e  ", iter, rNorm);
  end
  verbose && @printf("\n");

  status = on_boundary ? "on trust-region boundary" : (tired ? "maximum number of iterations exceeded" : "solution good enough given atol and rtol")
  stats = SimpleStats(solved, false, rNorms, T[], status);
  return (x, stats);
end


text:
# Testing SDR backends
using Test
using AbstractSDRs

# ----------------------------------------------------
# --- Testing radioSim back end
# ----------------------------------------------------
# This backend requires no SDR connected so we can use it by default
# Most of the test are also sone in the comon API call but set a specific test call here
# to test everything is OK on the API side and that specific RadioSim properties are OK (especially control that emulated received data is the one desired)
include("test_RadioSim.jl")

# ----------------------------------------------------
# --- Test functions
# ----------------------------------------------------
# Define here key behaviour common for all radio
# You can have a look here to find pratical examples of radio use
include("test_functions.jl")

# ----------------------------------------------------
# --- Common API call for all backend
# ----------------------------------------------------
# Issue is that we need a connected device to perform all test so we will do as follows
# => Try to scan for device. If detected then proceed to the test, otherwise drop the test
# For UHD, we need to hardcode a potential IP address in case of non broadcast ehternet address. If we detect the UHD device with scan(:uhd) the we use the associated IP address. Otherwise we use UHD_ADDRESS. It means that in case of non broadcast ethernet link, you should modify the following line according to your UHD IP address
# --- USRP address
global UHD_ADDRESS = "192.168.10.16"
#  This address may vary from time to time, not really sure how to handle dynamic testing

# --- Define test backend
# backends = [:radiosim;:uhd;:pluto;:rtlsdr]
backends = [:radiosim;:rtlsdr]
# backends = [:radiosim;:pluto]
for sdr ∈ backends
    # --- Flaging test
    println("######################################")
    println("# --- Testing $sdr backend ---")
    println("######################################")
    # --- Test scan
    str = scan(sdr)
    if isempty(str)
        println("No device found, rescan with IP address $UHD_ADDRESS")
        # If we have nothing, try with the USRP address
        str = scan(sdr;args="addr=$UHD_ADDRESS")
    end
    if isempty(str)
        # ----------------------------------------------------
        # --- We have find nothing, drop the rest of the tests
        # ----------------------------------------------------
        @warn "Unable to detect any SDR devices based on backend $sdr\n Abandon rest of tests"
    else
        # For UHD, we update str to be sure we can use the SDR latter on
        global UHD_ADDRESS = str[1]
        # ----------------------------------------------------
        # --- Testing stuff
        # ----------------------------------------------------
        @testset "Opening" begin
            # check_scan(sdr;UHD_ADDRESS)
            check_open(sdr)
        end

        @testset "Radio configuration" begin
            check_carrierFreq(sdr)
            check_samplingRate(sdr)
            check_gain(sdr)
        end

        @testset "Checking data retrieval"  begin
            check_recv(sdr)
            check_recv_preAlloc(sdr)
            check_recv_iterative(sdr)
        end

        @testset "Checking data transmission"  begin
            check_send(sdr)
        end
    end
end