text:
# Printing, show and parse
# ------------------------

Base.summary(seq::BioSequence{<:DNAAlphabet}) = string(length(seq), "nt ", "DNA Sequence")
Base.summary(seq::BioSequence{<:RNAAlphabet}) = string(length(seq), "nt ", "RNA Sequence")
Base.summary(seq::BioSequence{<:AminoAcidAlphabet}) = string(length(seq), "aa ", "Amino Acid Sequence")
Base.summary(seq::BioSequence{<:CharAlphabet}) = string(length(seq), "char ", "Char Sequence")

# Buffer type. Not exposed to user, so code should be kept simple and performant.
# B is true if it is buffered and false if it is not
mutable struct SimpleBuffer{B, T} <: IO
    len::UInt
    arr::Vector{UInt8}
    io::T
end

SimpleBuffer(io::IO) = SimpleBuffer{true, typeof(io)}(0, Vector{UInt8}(undef, 1024), io)
SimpleBuffer(io::IO, len) = SimpleBuffer{false, typeof(io)}(0, Vector{UInt8}(undef, len), io)

function Base.write(sb::SimpleBuffer{true}, byte::UInt8)
    sb.len ≥ 1024 && flush(sb)
    sb.len += 1
    @inbounds sb.arr[sb.len] = byte
end

function Base.write(sb::SimpleBuffer{false}, byte::UInt8)
    len = sb.len + 1
    sb.len = len
    @inbounds sb.arr[len] = byte
end

# Flush entire buffer to its io
@noinline function Base.flush(sb::SimpleBuffer{true})
    arr = sb.arr
    GC.@preserve arr unsafe_write(sb.io, pointer(arr), UInt(1024))
    sb.len = 0
end

# Flush all unflushed bytes to its io. Does not close source or make buffer unwritable
function Base.close(sb::SimpleBuffer)
    iszero(sb.len) && return 0
    arr = sb.arr
    GC.@preserve arr unsafe_write(sb.io, pointer(arr), sb.len)
    sb.len = 0
end

Base.print(io::IO, seq::BioSequence; width::Integer = 0) = _print(SimpleBuffer(io), seq, width)

function padded_length(len::Integer, width::Integer)
    den = ifelse(width < 1, typemax(Int), width)
    return len + div(len-1, den)
end

# Generic method. The different name allows subtypes of BioSequence to
# selectively call the generic print despite being more specific type
function _print(buffer::SimpleBuffer, seq::BioSequence, width::Integer)
    col = 0
    for x in seq
        col += 1
        if (width > 0) & (col > width)
            write(buffer, '\n')
            col = 1
        end
        print(buffer, x)
    end
    close(buffer)
    return nothing
end

Base.show(io::IO, seq::BioSequence) = showcompact(io, seq)

function Base.show(io::IO, ::MIME"text/plain", seq::BioSequence)
    println(io, summary(seq), ':')
    showcompact(io, seq)
end

function showcompact(io::IO, seq::BioSequence)
    # don't show more than this many characters
    # to avoid filling the screen with junk
    if isempty(seq)
        print(io, "< EMPTY SEQUENCE >")
    else
        width = displaysize()[2]
        if length(seq) > width
            half = div(width, 2)
            for i in 1:half-1
                print(io, seq[i])
            end
            print(io, '…')
            for i in lastindex(seq)-half+2:lastindex(seq)
                print(io, seq[i])
            end
        else
            for x in seq
                print(io, convert(Char, x))
            end
        end
    end
end

function string_compact(seq::BioSequence)
    buf = IOBuffer()
    showcompact(buf, seq)
    return String(take!(buf))
end

Base.parse(::Type{S}, str::AbstractString) where {S<:BioSequence} = convert(S, str)


text:
"""
$(TYPEDEF)
"""
struct DAESolution{T,N,uType,duType,uType2,DType,tType,P,A,ID,DE} <: AbstractDAESolution{T,N,uType}
  u::uType
  du::duType
  u_analytic::uType2
  errors::DType
  t::tType
  prob::P
  alg::A
  interp::ID
  dense::Bool
  tslocation::Int
  destats::DE
  retcode::Symbol
end
(sol::DAESolution)(t,deriv::Type=Val{0};idxs=nothing,continuity=:left) = sol.interp(t,idxs,deriv,sol.prob.p,continuity)
(sol::DAESolution)(v,t,deriv::Type=Val{0};idxs=nothing,continuity=:left) = sol.interp(v,t,idxs,deriv,sol.prob.p,continuity)

function build_solution(prob::AbstractDAEProblem, alg, t, u, du = nothing;
                        timeseries_errors = length(u) > 2,
                        dense = false,
                        dense_errors = dense,
                        calculate_error = true,
                        k = nothing,
                        interp = du === nothing ? LinearInterpolation(t, u) : HermiteInterpolation(t, u, du),
                        retcode = :Default,
                        destats = nothing,
                        kwargs...)
  T = eltype(eltype(u))
  N = length((size(prob.u0)..., length(u)))

  if has_analytic(prob.f)
    u_analytic = Vector{typeof(prob.u0)}()
    errors = Dict{Symbol,real(eltype(prob.u0))}()

    sol = DAESolution{T,N,typeof(u),typeof(du),typeof(u_analytic),typeof(errors),typeof(t),
                      typeof(prob),typeof(alg),typeof(interp),typeof(destats)}(
                        u,du,u_analytic,errors,t,prob,alg,interp,dense,0,destats,retcode)

    if calculate_error
      calculate_solution_errors!(sol;timeseries_errors=timeseries_errors,dense_errors=dense_errors)
    end
    sol
  else
    DAESolution{T,N,typeof(u),typeof(du),Nothing,Nothing,typeof(t),
                typeof(prob),typeof(alg),typeof(interp),typeof(destats)}(
                  u,du,nothing,nothing,t,prob,alg,interp,dense,0,destats,retcode)
  end
end

function calculate_solution_errors!(sol::AbstractDAESolution;
                                    fill_uanalytic = true, timeseries_errors = true, dense_errors = true)
  prob = sol.prob
  f = prob.f

  if fill_uanalytic
    for i in 1:size(sol.u, 1)
      push!(sol.u_analytic, f.analytic(prob.du0, prob.u0, prob.p, sol.t[i]))
    end
  end

  save_everystep = length(sol.u) > 2
  if !isempty(sol.u_analytic)
    sol.errors[:final] = norm(recursive_mean(abs.(sol.u[end] - sol.u_analytic[end])))

    if save_everystep && timeseries_errors
      sol.errors[:l∞] = norm(maximum(vecvecapply(x -> abs.(x), sol.u - sol.u_analytic)))
      sol.errors[:l2] = norm(sqrt(recursive_mean(vecvecapply(x -> float.(x).^2, sol.u - sol.u_analytic))))
      if sol.dense && dense_errors
        densetimes = collect(range(sol.t[1]; stop = sol.t[end], length = 100))
        interp_u = sol(densetimes)
        interp_analytic = VectorOfArray([f.analytic(prob.du0, prob.u0, prob.p, t) for t in densetimes])
        sol.errors[:L∞] = norm(maximum(abs.(interp_u .- interp_analytic)))
        sol.errors[:L2] = norm(sqrt(mean((interp_u .- interp_analytic).^2)))
      end
    end
  end

  nothing
end

function build_solution(sol::AbstractDAESolution{T,N},u_analytic,errors) where {T,N}
  DAESolution{T,N,typeof(sol.u),typeof(sol.du),typeof(u_analytic),typeof(errors),typeof(sol.t),
                     typeof(sol.prob),typeof(sol.alg),typeof(sol.interp),typeof(sol.destats)}(
                     sol.u,sol.du,u_analytic,errors,sol.t,
              sol.prob,sol.alg,sol.interp,sol.dense,sol.tslocation,sol.destats,sol.retcode)
end

function solution_new_retcode(sol::AbstractDAESolution{T,N},retcode) where {T,N}
  DAESolution{T,N,typeof(sol.u),typeof(sol.du),typeof(sol.u_analytic),
              typeof(sol.errors),typeof(sol.t),
              typeof(sol.prob),typeof(sol.alg),typeof(sol.interp),typeof(sol.destats)}(
              sol.u,sol.du,sol.u_analytic,sol.errors,sol.t,
              sol.prob,sol.alg,sol.interp,sol.dense,sol.tslocation,sol.destats,retcode)
end

function solution_new_tslocation(sol::AbstractDAESolution{T,N},tslocation) where {T,N}
  DAESolution{T,N,typeof(sol.u),typeof(sol.du),typeof(sol.u_analytic),
              typeof(sol.errors),typeof(sol.t),
              typeof(sol.prob),typeof(sol.alg),typeof(sol.interp),typeof(sol.destats)}(
              sol.u,sol.du,sol.u_analytic,sol.errors,sol.t,
              sol.prob,sol.alg,sol.interp,sol.dense,tslocation,sol.destats,sol.retcode)
end

function solution_slice(sol::AbstractDAESolution{T,N},I) where {T,N}
  DAESolution{T,N,typeof(sol.u),typeof(sol.du),typeof(sol.u_analytic),
              typeof(sol.errors),typeof(sol.t),
              typeof(sol.prob),typeof(sol.alg),typeof(sol.interp),typeof(sol.destats)}(
              sol.u[I],sol.du[I],
              sol.u_analytic === nothing ? nothing : sol.u_analytic[I],
              sol.errors,sol.t[I],
              sol.prob,sol.alg,
              sol.interp,false,
              sol.tslocation,sol.destats,sol.retcode)
end

ext: jl


text:
using AES, Test

k = [0x2b, 0x7e, 0x15, 0x16, 0x28, 0xae, 0xd2, 0xa6, 0xab, 0xf7, 0x15, 0x88, 0x09, 0xcf, 0x4f, 0x3c]
key = AES128Key(k)

plaintext = "The quick brown fox jumps over the lazy dog."
plaintext8 = "The quick brown fox jumps over the 😴 dog."

cipher = AESCipher(;key_length=128, mode=AES.ECB, key=key)
ct = encrypt(plaintext, cipher)
pt = decrypt(ct, cipher)
ct8 = encrypt(plaintext8, cipher)
pt8 = decrypt(ct8, cipher)

ciphertext = [0x16, 0xfa, 0x65, 0x87, 0x31, 0x00, 0x2a, 0xd6,
			  0xe3, 0x4a, 0x2f, 0xa0, 0x0f, 0x29, 0x0d, 0x9f,
			  0x97, 0x4f, 0x7b, 0xac, 0x10, 0x45, 0x57, 0x4b,
			  0x74, 0xc2, 0x04, 0x9e, 0x65, 0xd2, 0xa8, 0x89,
			  0x41, 0x37, 0xa2, 0x65, 0xbb, 0xb0, 0x1b, 0x60,
			  0x88, 0x95, 0x89, 0x9d, 0x38, 0x1e, 0xe7, 0x36]

@testset "ASCII" begin
	@test pt == transcode(UInt8,plaintext)
	@test ct.data == ciphertext
end

@testset "UTF-8" begin
	@test pt8 == transcode(UInt8,plaintext8)
end


text:
using Test

@testset "SSPs Component" begin
    include("test_SSPs.jl")
end

@testset "RegionAggregatorSum Component" begin
end

@testset "Coupled" begin
end


text:
using NCDatasets
using Test
sz = (4,5)
filename = tempname()
#filename = "/tmp/mytest.nc"

NCDatasets.NCDataset(filename,"c") do ds

    ds.dim["lon"] = sz[1]
    ds.dim["lat"] = sz[2]

    v = NCDatasets.defVar(ds,"temperature",Float32,("lon","lat"),
                          attrib = ["long_name" => "Temperature",
                                    "test_vector_attrib" => [1,2,3]])

    # write attributes
    v.attrib["units"] = "degree Celsius"
    v.attrib["comment"] = "this is a string attribute with unicode Ω ∈ ∑ ∫ f(x) dx "

    # check presence of attribute
    @test NCDatasets.haskey(v.attrib,"comment")
    @test "comment" in v.attrib
    @test v.attrib["long_name"] == "Temperature"
    @test v.attrib[:long_name] == "Temperature"
    @test v.attrib["test_vector_attrib"] == [1,2,3]
    @test v.attrib["comment"] == "this is a string attribute with unicode Ω ∈ ∑ ∫ f(x) dx "

    @test NCDatasets.get(v.attrib,"does-not-exist","default") == "default"
    @test NCDatasets.get(v.attrib,"units","default") == "degree Celsius"

    # test deletion of attributes
    v.attrib["todelete"] = "foobar"
    @test haskey(v.attrib,"todelete")
    delete!(v.attrib,"todelete")
    @test !haskey(v.attrib,"todelete")

    for T in [UInt8,Int8,UInt16,Int16,UInt32,Int32,UInt64,Int64,Float32,Float64,
              String,Char]
        # scalar attribute
        name = "scalar-attrib-$T"

        refdata =
            if T == Char
                'a'
            elseif T == String
                "abc"
            else
                123
            end

        v.attrib[name] = T(refdata)
        attval = v.attrib[name]

        if T == Char
            # a single Char is returned as strings
            @test typeof(attval) == String
            @test attval[1] == refdata
        else
            if T == Int64
                # not supported in NetCDF, converted as Int32
                @test typeof(attval) in [Int32,Int64]
            else
                @test typeof(attval) == T
            end

            @test attval == refdata
        end

        # vector attribute
        name = "vector-attrib-$T"

        refvecdata =
            if T == Char
                ['a','b']
            elseif T == String
                ["abc","xyz"]
            else
                [1,2,3,4]
            end

        attval = T.(refvecdata)
        attval = attval

        if T == Char
            # vector of Char are returned as strings
            @test eltype(attval) == T
            @test Vector{Char}(attval) == refvecdata
        else
            if T == Int64
                # not supported in NetCDF, converted as Int32
                @test eltype(attval) in [Int32,Int64]
            else
                @test eltype(attval) == T
            end

            @test attval == refvecdata
        end

    end

    # arrays cannot be attributes
    @test_throws ErrorException v.attrib["error_attrib"] = zeros(2,2)

    # symbols in the attrib dict
    foo = NCDatasets.defVar(ds,"foovar",Int64,("lon","lat"),
                            attrib = [:long_name => "foo variable"])
    @test foo.attrib["long_name"] == "foo variable"

end

rm(filename)

filename = tempname()

NCDatasets.NCDataset(filename,"c",format = :netcdf3_classic) do ds
    # test deletion of attributes
    ds.attrib["todelete"] = "foobar"
end

NCDatasets.NCDataset(filename,"a") do ds
    @test haskey(ds.attrib,"todelete")
    delete!(ds.attrib,"todelete")
    @test !haskey(ds.attrib,"todelete")
end

rm(filename)

#filename = "/tmp/mytest.nc"


text:
"""
$(SIGNATURES)
Wait for the requests to complete, poll until done.

Args:
    requestIds (List[str]): The request IDs that should be polled.
    maxSeconds (int, optional): Maximum wait time. Defaults to 60.
    expectedStatus (str, optional): Expected status message per request.
        Defaults to "Complete".
"""
function waitForCompletion(
    navAbilityClient::NavAbilityClient,
    requestIds::AbstractVector{<:AbstractString};
    maxSeconds::Integer = 120,
    expectedStatuses::Union{Nothing,<:AbstractVector{<:AbstractString}} = nothing,
    exceptionMessage::AbstractString = "Requests did not complete in time")
#
    if expectedStatuses == Nothing
        expectedStatuses = ["Complete", "Failed"]
    end
    wait_time = maxSeconds
    tasksComplete = false
    while !tasksComplete
        statuses = values(fetch( getStatusesLatest(navAbilityClient, requestIds) ))
        tasksComplete = all(s["state"] in expectedStatuses for s in statuses)
        if tasksComplete
            return
        else
            sleep(2)
            wait_time -= 2
            wait_time <= 0 && throw(error(exceptionMessage))
        end
    end
end

# Dispatch for vector of Tasks
function waitForCompletion(
    navAbilityClient::NavAbilityClient,
    requestIds::AbstractVector{<:Task};
    kw...)
    #
    rids = fetch.(requestIds) .|> string
    waitForCompletion(navAbilityClient, rids; kw...)
end

# helper functions to construct for most likely user object
function GraphVizApp(ct::Client; variableStartsWith=nothing)
    suffix = ""
    if !(variableStartsWith isa Nothing)
        suffix *= "&variableStartsWith"
        if 0<length(variableStartsWith)
            suffix *= "=$variableStartsWith"
        end
    end
    GraphVizApp("https://app.navability.io/cloud/graph/?userId=$(ct.userId)&robotStartsWith=$(ct.robotId)&sessionStartsWith=$(ct.sessionId)"*suffix)
end
function MapVizApp(ct::Client; variableStartsWith=nothing)
    suffix = ""
    if !(variableStartsWith isa Nothing)
        suffix *= "&variableStartsWith"
        if 0<length(variableStartsWith)
            suffix *= "=$variableStartsWith"
        end
    end
    MapVizApp("https://app.navability.io/cloud/map/?userId=$(ct.userId)&robotStartsWith=$(ct.robotId)&sessionStartsWith=$(ct.sessionId)"*suffix)
end

# overload show dispatch for convenient static and interactive links to NavAbility App visualizations
Base.show(io::IO, ::MIME"text/plain", gv::Union{GraphVizApp,MapVizApp}) = println(gv.url)
Base.show(io::IO, ::MIME"text/markdown", gv::GraphVizApp) = display("text/markdown", "[![Navigate to Factor Graph]($assetGraphVizImg)]($(gv.url))")
Base.show(io::IO, ::MIME"text/markdown", gv::MapVizApp) = display("text/markdown", "[![Navigate to Factor Graph]($assetGeomVizImg)]($(gv.url))")


text:
include("query_macros.jl")

"""
A RQL datum, stored in `datum` below.

DATUM = 1;
"""
@reql_zero(1, datum,
  ReqlDatum)

"MAKE_ARRAY = 2; // DATUM... -> ARRAY"
@reql_onearr(2, make_array,
  ReqlDatum,
  ReqlArray)

# Compound types

"""
Takes an integer representing a variable and returns the value stored
in that variable.  It's the responsibility of the client to translate
from their local representation of a variable to a unique _non-negative_
integer for that variable.  (We do it this way instead of letting
clients provide variable names as strings to discourage
variable-capturing client libraries, and because it's more efficient

VAR = 10; // !NUMBER -> DATUM
"""
@reql_one(10, var,
  ReqlString,
  ReqlDatum)

"""
Takes some javascript code and executes it.

JAVASCRIPT = 11; // STRING {timeout: !NUMBER} -> DATUM |
                 // STRING {timeout: !NUMBER} -> Function(\*)
"""
@reql_one(11, javascript,
  ReqlString,
  ReqlDatum)

"alias for RethinkDB.javascript"
@reql_one(11, js,
  ReqlString,
  ReqlDatum)

"UUID = 169; // () -> DATUM"
@reql_zero(169, datum,
  ReqlDatum)

"""
Takes an HTTP URL and gets it.  If the get succeeds and
returns valid JSON, it is converted into a DATUM

HTTP = 153; // STRING {data: OBJECT | STRING,
            //         timeout: !NUMBER,
            //         method: STRING,
            //         params: OBJECT,
            //         header: OBJECT | ARRAY,
            //         attempts: NUMBER,
            //         redirects: NUMBER,
            //         verify: BOOL,
            //         page: FUNC | STRING,
            //         page_limit: NUMBER,
            //         auth: OBJECT,
            //         result_format: STRING,
            //         } -> STRING | STREAM
"""
@reql_one(153, http,
  ReqlString,
  ReqlString)

"""
Takes a string and throws an error with that message.
Inside of a `default` block, you can omit the first
argument to rethrow whatever error you catch (this is most
useful as an argument to the `default` filter optarg).

ERROR = 12; // STRING -> Error | -> Error
"""
@reql_one(12, error,
  ReqlString,
  ReqlError)

"""
Takes nothing and returns a reference to the implicit variable.

IMPLICIT_VAR = 13; // -> DATUM
"""
@reql_zero(13, implicit_var,
  ReqlDatum)

# Data Operators

"""
Returns a reference to a database.

DB = 14; // STRING -> Database
"""
@reql_one(14, db,
  ReqlString,
  ReqlDatabase)

"""
Returns a reference to a table.

TABLE = 15; // Database, STRING, {read_mode:STRING, identifier_format:STRING} -> Table
            // STRING, {read_mode:STRING, identifier_format:STRING} -> Table
"""
@reql_one_two(15, table,
  ReqlDatabase, ReqlString,
  ReqlTable)
@reql_one(15, table,
  ReqlString,
  ReqlTable)

"""
Gets a single element from a table by its primary or a secondary key.

GET = 16; // Table, STRING -> SingleSelection | Table, NUMBER -> SingleSelection |
          // Table, STRING -> NULL            | Table, NUMBER -> NULL |
"""
@reql_one_two(16, get,
  ReqlTable, ReqlString,
  ReqlSingleSelection)
@reql_one_two(16, get,
  ReqlTable, ReqlNumber,
  ReqlSingleSelection)

"GET_ALL = 78; // Table, DATUM..., {index:!STRING} => ARRAY"
@reql_one_twoarr(78, get_all,
  ReqlTable, ReqlDatum,
  ReqlArray)
@reql_one(78, get_all,
  ReqlTable,
  ReqlArray)

# Simple DATUM Ops

"EQ = 17; // DATUM... -> BOOL"
@reql_onearr(17, eq,
  ReqlDatum,
  ReqlBool)

"NE = 18; // DATUM... -> BOOL"
@reql_onearr(18, ne,
  ReqlDatum,
  ReqlBool)

"LT = 19; // DATUM... -> BOOL"
@reql_onearr(19, lt,
  ReqlDatum,
  ReqlBool)

"LE = 20; // DATUM... -> BOOL"
@reql_onearr(20, le,
  ReqlDatum,
  ReqlBool)

"GT = 21; // DATUM... -> BOOL"
@reql_onearr(21, gt,
  ReqlDatum,
  ReqlBool)

"GE = 22; // DATUM... -> BOOL"
@reql_onearr(22, ge,
  ReqlDatum,
  ReqlBool)

"NOT = 23; // BOOL -> BOOL"
@reql_one(23, not,
  ReqlBool,
  ReqlBool)

"""
ADD can either add two numbers or concatenate two arrays.

ADD = 24; // NUMBER... -> NUMBER | STRING... -> STRING
"""
@reql_onearr(24, add,
  ReqlNumber,
  ReqlNumber)
@reql_onearr(24, add,
  ReqlString,
  ReqlString)

"SUB = 25; // NUMBER... -> NUMBER"
@reql_onearr(25, sub,
  ReqlNumber,
  ReqlNumber)

"MUL = 26; // NUMBER... -> NUMBER"
@reql_onearr(26, mul,
  ReqlNumber,
  ReqlNumber)

"DIV = 27; // NUMBER... -> NUMBER"
@reql_onearr(27, div,
  ReqlNumber,
  ReqlNumber)

"MOD = 28; // NUMBER, NUMBER -> NUMBER"
@reql_one_two(28, mod,
  ReqlNumber, ReqlNumber,
  ReqlNumber)

"FLOOR = 183;    // NUMBER -> NUMBER"
@reql_one(183, floor,
  ReqlNumber,
  ReqlNumber)

"CEIL = 184;     // NUMBER -> NUMBER"
@reql_one(184, ceil,
  ReqlNumber,
  ReqlNumber)

"ROUND = 185;    // NUMBER -> NUMBER"
@reql_one(185, round,
  ReqlNumber,
  ReqlNumber)

# DATUM Array Ops

"""
Append a single element to the end of an array (like `snoc`).

APPEND = 29; // ARRAY, DATUM -> ARRAY
"""
@reql_one_two(29, append,
  ReqlArray, ReqlDatum,
  ReqlArray)

"""
Prepend a single element to the end of an array (like `cons`).

PREPEND = 80; // ARRAY, DATUM -> ARRAY
"""
@reql_one_two(80, prepend,
  ReqlArray, ReqlDatum,
  ReqlArray)

"""
Remove the elements of one array from another array.

DIFFERENCE = 95; // ARRAY, ARRAY -> ARRAY
"""
@reql_one_two(95, difference,
  ReqlArray, ReqlArray,
  ReqlArray)

# DATUM Set Ops

#=
Set ops work on arrays. They don't use actual sets and thus have
performance characteristics you would expect from arrays rather than
from sets. All set operations have the post condition that they
array they return contains no duplicate values.
=#

"SET_INSERT = 88; // ARRAY, DATUM -> ARRAY"
@reql_one_two(88, set_insert,
  ReqlArray, ReqlDatum,
  ReqlArray)

"SET_INTERSECTION = 89; // ARRAY, ARRAY -> ARRAY"
@reql_one_two(89, set_intersection,
  ReqlArray, ReqlArray,
  ReqlArray)

"SET_UNION = 90; // ARRAY, ARRAY -> ARRAY"
@reql_one_two(90, set_union,
  ReqlArray, ReqlArray,
  ReqlArray)

"SET_DIFFERENCE = 91; // ARRAY, ARRAY -> ARRAY"
@reql_one_two(91, set_difference,
  ReqlArray, ReqlArray,
  ReqlArray)

"SLICE = 30; // Sequence, NUMBER, NUMBER -> Sequence"
@reql_one_two_three(30, slice,
  ReqlSequence, ReqlNumber, ReqlNumber,
  ReqlSequence)

"SKIP = 70; // Sequence, NUMBER -> Sequence"
@reql_one_two(70, skip,
  ReqlSequence, ReqlNumber,
  ReqlSequence)

"LIMIT = 71; // Sequence, NUMBER -> Sequence"
@reql_one_two(71, limit,
  ReqlSequence, ReqlNumber,
  ReqlSequence)

"OFFSETS_OF = 87; // Sequence, DATUM -> Sequence | Sequence, Function(1) -> Sequence"
@reql_one_two(87, offsets_of,
  ReqlSequence, ReqlDatum,
  ReqlSequence)
@reql_one_two(87, offsets_of,
  ReqlSequence, ReqlFunction1,
  ReqlSequence)

# TODO: Function or Datum
"CONTAINS = 93; // Sequence, (DATUM | Function(1))... -> BOOL"
@reql_one_twoarr(93, contains,
  ReqlSequence, ReqlDatum,
  ReqlBool)

# Stream/Object Ops

"""
Get a particular field from an object, or map that over a
sequence.

GET_FIELD = 31; // OBJECT, STRING -> DATUM
                // | Sequence, STRING -> Sequence
"""
@reql_one_two(31, get_field,
  ReqlObject, ReqlString,
  ReqlDatum)
@reql_one_two(31, get_field,
  ReqlSequence, ReqlString,
  ReqlSequence)

"""
Return an array containing the keys of the object.

KEYS = 94; // OBJECT -> ARRAY
"""
@reql_one(94, keys,
  ReqlObject,
  ReqlArray)

"""
Return an array containing the values of the object.

VALUES = 186; // OBJECT -> ARRAY
"""
@reql_one(186, values,
  ReqlObject,
  ReqlArray)

#=
TODO
Creates an object

OBJECT = 143; // STRING, DATUM, ... -> OBJECT
=#

"""
Check whether an object contains all the specified fields,
or filters a sequence so that all objects inside of it
contain all the specified fields.

HAS_FIELDS = 32; // OBJECT, Pathspec... -> BOOL
"""
@reql_one_two(32, has_fields,
  ReqlObject, ReqlPathSpec,
  ReqlBool)

#=
TODO
x.with_fields(...) <=> x.has_fields(...).pluck(...)
WITH_FIELDS = 96; // Sequence, Pathspec... -> Sequence
=#

"""
Get a subset of an object by selecting some attributes to preserve,
or map that over a sequence.  (Both pick and pluck, polymorphic.)

PLUCK = 33; // Sequence, Pathspec... -> Sequence | OBJECT, Pathspec... -> OBJECT
"""
@reql_one_twoarr(33, pluck,
  ReqlSequence, ReqlPathSpec,
  ReqlSequence)
@reql_one_twoarr(33, pluck,
  ReqlObject, ReqlPathSpec,
  ReqlObject)

#=
TODO
Get a subset of an object by selecting some attributes to discard, or
map that over a sequence.  (Both unpick and without, polymorphic.)
WITHOUT = 34; // Sequence, Pathspec... -> Sequence | OBJECT, Pathspec... -> OBJECT
=#

#=
TODO
Merge objects (right-preferential)
MERGE = 35; // OBJECT... -> OBJECT | Sequence -> Sequence
=#

# Sequence Ops

#=
Get all elements of a sequence between two values.
Half-open by default, but the openness of either side can be
changed by passing 'closed' or 'open for `right_bound` or
`left_bound`.
=#

#=
BETWEEN_DEPRECATED = 36; // Deprecated version of between, which allows `null` to specify unboundedness
                         // With the newer version, clients should use `r.minval` and `r.maxval` for unboundedness
=#

#=
TODO
 BETWEEN   = 182; // StreamSelection, DATUM, DATUM, {index:!STRING, right_bound:STRING, left_bound:STRING} -> StreamSelection
=#

#=
TODO
REDUCE = 37; // Sequence, Function(2) -> DATUM
=#

#=
TODO
MAP = 38; // Sequence, Function(1) -> Sequence
          // The arity of the function should be
          // Sequence..., Function(sizeof...(Sequence)) -> Sequence
=#

"""
Filter a sequence with either a function or a shortcut
object (see API docs for details).  The body of FILTER is
wrapped in an implicit `.default(false)`, and you can
change the default value by specifying the `default`
optarg.  If you make the default `r.error`, all errors
caught by `default` will be rethrown as if the `default`
did not exist.

FILTER = 39; // Sequence, Function(1), {default:DATUM} -> Sequence |
             // Sequence, OBJECT, {default:DATUM} -> Sequence
"""
@reql_one_two(39, filter,
  ReqlSequence, ReqlFunction1,
  ReqlSequence)
@reql_one_two(39, filter,
  ReqlSequence, ReqlObject,
  ReqlSequence)

#=
TODO
Map a function over a sequence and then concatenate the results together.

CONCAT_MAP = 40; // Sequence, Function(1) -> Sequence
=#

#=
TODO
Order a sequence based on one or more attributes.

ORDER_BY = 41; // Sequence, (!STRING | Ordering)..., {index: (!STRING | Ordering)} -> Sequence
=#

"""
Get all distinct elements of a sequence (like `uniq`).

DISTINCT = 42; // Sequence -> Sequence
"""
@reql_one(42, distinct,
  ReqlSequence,
  ReqlSequence)

"""
Count the number of elements in a sequence, or only the elements that match
a given filter.

COUNT = 43; // Sequence -> NUMBER | Sequence, DATUM -> NUMBER | Sequence, Function(1) -> NUMBER
"""
@reql_one(43, count,
  ReqlSequence,
  ReqlNumber)
@reql_one_two(43, count,
  ReqlSequence, ReqlDatum,
  ReqlNumber)
@reql_one_two(43, count,
  ReqlSequence, ReqlFunction,
  ReqlNumber)

"IS_EMPTY = 86; // Sequence -> BOOL"
@reql_one(86, is_empty,
  ReqlSequence,
  ReqlBool)

#=
TODO
Take the union of multiple sequences (preserves duplicate elements! (use distinct)).
UNION = 44; // Sequence... -> Sequence
=#

"""
Get the Nth element of a sequence.

NTH = 45; // Sequence, NUMBER -> DATUM
"""
@reql_one_two(45, nth,
  ReqlSequence, ReqlNumber,
  ReqlDatum)

"""
do NTH or GET_FIELD depending on target object

BRACKET = 170; // Sequence | OBJECT, NUMBER | STRING -> DATUM
"""
@reql_one_two(170, bracket,
  ReqlSequence, ReqlNumber,
  ReqlDatum)
@reql_one_two(170, bracket,
  ReqlSequence, ReqlString,
  ReqlDatum)
@reql_one_two(170, bracket,
    ReqlDatum, ReqlString,
    ReqlDatum)
@reql_one_two(170, bracket,
  ReqlObject, ReqlString,
  ReqlDatum)

# OBSOLETE_GROUPED_MAPREDUCE = 46;

# OBSOLETE_GROUPBY = 47;

#=
TODO
INNER_JOIN = 48; // Sequence, Sequence, Function(2) -> Sequence
=#

#=
TODO
OUTER_JOIN = 49; // Sequence, Sequence, Function(2) -> Sequence
=#

#=
TODO
An inner-join that does an equality comparison on two attributes.

EQ_JOIN = 50; // Sequence, !STRING, Sequence, {index:!STRING} -> Sequence
=#

"ZIP = 72; // Sequence -> Sequence"
@reql_one(72, zip,
  ReqlSequence,
  ReqlSequence)

"""
RANGE = 173; // -> Sequence                        [0, +inf)
             // NUMBER -> Sequence                 [0, a)
             // NUMBER, NUMBER -> Sequence         [a, b)
"""
@reql_zero(173, range,
  ReqlSequence)
@reql_one(173, range,
  ReqlNumber,
  ReqlSequence)
@reql_one_two(173, range,
  ReqlNumber, ReqlNumber,
  ReqlSequence)

# Array Ops

"""
Insert an element in to an array at a given index.

INSERT_AT = 82; // ARRAY, NUMBER, DATUM -> ARRAY
"""
@reql_one_two_three(82, insert_at,
  ReqlArray, ReqlNumber, ReqlDatum,
  ReqlArray)

"""
Remove an element at a given index from an array.

DELETE_AT = 83; // ARRAY, NUMBER -> ARRAY |
                // ARRAY, NUMBER, NUMBER -> ARRAY
"""
@reql_one_two(83, delete_at,
  ReqlArray, ReqlNumber,
  ReqlArray)
@reql_one_two_three(83, delete_at,
  ReqlArray, ReqlNumber, ReqlNumber,
  ReqlArray)

"""
Change the element at a given index of an array.

CHANGE_AT = 84; // ARRAY, NUMBER, DATUM -> ARRAY
"""
@reql_one_two_three(84, change_at,
  ReqlArray, ReqlNumber, ReqlDatum,
  ReqlArray)

"""
Splice one array in to another array.

SPLICE_AT = 85; // ARRAY, NUMBER, ARRAY -> ARRAY
"""
@reql_one_two_three(85, splice_at,
  ReqlArray, ReqlNumber, ReqlArray,
  ReqlArray)

# Type Ops

"""
Coerces a datum to a named type (e.g. "bool").
If you previously used `stream_to_array`, you should use this instead
with the type "array".

COERCE_TO = 51; // Top, STRING -> Top
"""
@reql_one_two(51, coerce_to,
  ReqlTop, ReqlString,
  ReqlTop)

"""
Returns the named type of a datum (e.g. TYPE_OF(true) = "BOOL")

TYPE_OF = 52; // Top -> STRING
"""
@reql_one(52, type_of,
  ReqlTop,
  ReqlString)

# Write Ops (the OBJECTs contain data about number of errors etc.)

"""
TODO: func
Updates all the rows in a selection.  Calls its Function with the row
to be updated, and then merges the result of that call.

UPDATE = 53; // StreamSelection, Function(1), {non_atomic:BOOL, durability:STRING, return_changes:BOOL} -> OBJECT |
             // SingleSelection, Function(1), {non_atomic:BOOL, durability:STRING, return_changes:BOOL} -> OBJECT |
             // StreamSelection, OBJECT,      {non_atomic:BOOL, durability:STRING, return_changes:BOOL} -> OBJECT |
             // SingleSelection, OBJECT,      {non_atomic:BOOL, durability:STRING, return_changes:BOOL} -> OBJECT
"""
@reql_one_two(53, update,
  ReqlStreamSelection, ReqlFunction1,
  ReqlObject)
@reql_one_two(53, update,
  ReqlSingleSelection, ReqlFunction1,
  ReqlObject)
@reql_one_two(53, update,
  ReqlStreamSelection, ReqlObject,
  ReqlObject)
@reql_one_two(53, update,
  ReqlSingleSelection, ReqlObject,
  ReqlObject)

"""
Deletes all the rows in a selection.

DELETE = 54; // StreamSelection, {durability:STRING, return_changes:BOOL} -> OBJECT | SingleSelection -> OBJECT
"""
@reql_one(54, delete,
  ReqlStreamSelection,
  ReqlObject)
@reql_one(54, delete,
  ReqlSingleSelection,
  ReqlObject)

#=
TODO
Replaces all the rows in a selection.  Calls its Function with the row
to be replaced, and then discards it and stores the result of that
call.

REPLACE = 55; // StreamSelection, Function(1), {non_atomic:BOOL, durability:STRING, return_changes:BOOL} -> OBJECT | SingleSelection, Function(1), {non_atomic:BOOL, durability:STRING, return_changes:BOOL} -> OBJECT
=#

"""
Inserts into a table.  If `conflict` is replace, overwrites
entries with the same primary key.  If `conflict` is
update, does an update on the entry.  If `conflict` is
error, or is omitted, conflicts will trigger an error.

INSERT = 56; // Table, OBJECT, {conflict:STRING, durability:STRING, return_changes:BOOL} -> OBJECT |
             // Table, Sequence, {conflict:STRING, durability:STRING, return_changes:BOOL} -> OBJECT
"""
@reql_one_two(56, insert,
  ReqlTable, ReqlObject,
  ReqlObject)
@reql_one_two(56, insert,
  ReqlTable, ReqlSequence,
  ReqlObject)

# Administrative OPs

"""
Creates a database with a particular name.

DB_CREATE = 57; // STRING -> OBJECT
"""
@reql_one(57, db_create,
  ReqlString,
  ReqlObject)

"""
Drops a database with a particular name.

DB_DROP = 58; // STRING -> OBJECT
"""
@reql_one(58, db_drop,
  ReqlString,
  ReqlObject)

"""
Lists all the databases by name.  (Takes no arguments)

DB_LIST = 59; // -> ARRAY
"""
@reql_zero(59, db_list,
  ReqlArray)

"""
Creates a table with a particular name in a particular
database.  (You may omit the first argument to use the
default database.)

TABLE_CREATE = 60; // Database, STRING, {primary_key:STRING, shards:NUMBER, replicas:NUMBER, primary_replica_tag:STRING} -> OBJECT
                   // Database, STRING, {primary_key:STRING, shards:NUMBER, replicas:OBJECT, primary_replica_tag:STRING} -> OBJECT
                   // STRING, {primary_key:STRING, shards:NUMBER, replicas:NUMBER, primary_replica_tag:STRING} -> OBJECT
                   // STRING, {primary_key:STRING, shards:NUMBER, replicas:OBJECT, primary_replica_tag:STRING} -> OBJECT
"""
@reql_one_two(60, table_create,
  ReqlDatabase, ReqlString,
  ReqlObject)
@reql_one(60, table_create,
  ReqlString,
  ReqlObject)

"""
Drops a table with a particular name from a particular
database.  (You may omit the first argument to use the
default database.)

TABLE_DROP = 61; // Database, STRING -> OBJECT
                 // STRING -> OBJECT
"""
@reql_one_two(61, table_drop,
  ReqlDatabase, ReqlString,
  ReqlObject)
@reql_one(61, table_drop,
  ReqlString,
  ReqlObject)

"""
Lists all the tables in a particular database.  (You may
omit the first argument to use the default database.)

TABLE_LIST = 62; // Database -> ARRAY
                 //  -> ARRAY
"""
@reql_one(62, table_list,
  ReqlDatabase,
  ReqlArray)
@reql_zero(62, table_list,
  ReqlArray)

"""
Returns the row in the `rethinkdb.table_config` or `rethinkdb.db_config` table
that corresponds to the given database or table.
CONFIG  = 174; // Database -> SingleSelection
               // Table -> SingleSelection
"""
@reql_one(174, config,
  ReqlDatabase,
  ReqlTable)

"""
Returns the row in the `rethinkdb.table_status` table that corresponds to the
given table.

STATUS  = 175; // Table -> SingleSelection
"""
@reql_one(175, status,
  ReqlTable,
  ReqlSingleSelection)

"""
Called on a table, waits for that table to be ready for read/write operations.
Called on a database, waits for all of the tables in the database to be ready.
Returns the corresponding row or rows from the `rethinkdb.table_status` table.

WAIT    = 177; // Table -> OBJECT
               // Database -> OBJECT
"""
@reql_one(177, wait,
  ReqlTable,
  ReqlObject)
@reql_one(177, wait,
  ReqlDatabase,
  ReqlObject)

"""
Generates a new config for the given table, or all tables in the given database
The `shards` and `replicas` arguments are required. If `emergency_repair` is
specified, it will enter a completely different mode of repairing a table
which has lost half or more of its replicas.

RECONFIGURE   = 176; // Database|Table, {shards:NUMBER, replicas:NUMBER [,
                     //                  dry_run:BOOLEAN]
                     //                 } -> OBJECT
                     // Database|Table, {shards:NUMBER, replicas:OBJECT [,
                     //                  primary_replica_tag:STRING,
                     //                  nonvoting_replica_tags:ARRAY,
                     //                  dry_run:BOOLEAN]
                     //                 } -> OBJECT
                     // Table, {emergency_repair:STRING, dry_run:BOOLEAN} -> OBJECT
"""
@reql_one(176, reconfigure,
  ReqlDatabase,
  ReqlObject)
@reql_one(176, reconfigure,
  ReqlTable,
  ReqlObject)

"""
Balances the table's shards but leaves everything else the same. Can also be
applied to an entire database at once.

REBALANCE     = 179; // Table -> OBJECT
                     // Database -> OBJECT
"""
@reql_one(179, rebalance,
  ReqlTable,
  ReqlObject)
@reql_one(179, rebalance,
  ReqlTable,
  ReqlDatabase)

"""
Ensures that previously issued soft-durability writes are complete and
written to disk.

SYNC = 138; // Table -> OBJECT
"""
@reql_one(138, sync,
  ReqlTable,
  ReqlObject)

"""
Set global, database, or table-specific permissions

GRANT         = 188; //          -> OBJECT
                     // Database -> OBJECT
                     // Table    -> OBJECT
"""
@reql_zero(188, grant,
  ReqlObject)
@reql_one(188, grant,
  ReqlDatabase,
  ReqlObject)
@reql_one(188, grant,
  ReqlTable,
  ReqlObject)

# Secondary indexes OPs

#=
TODO
Creates a new secondary index with a particular name and definition.
INDEX_CREATE = 75; // Table, STRING, Function(1), {multi:BOOL} -> OBJECT
=#

"""
Drops a secondary index with a particular name from the specified table.

INDEX_DROP = 76; // Table, STRING -> OBJECT
"""
@reql_one_two(76, index_drop,
  ReqlTable, ReqlString,
  ReqlObject)

"""
Lists all secondary indexes on a particular table.

INDEX_LIST = 77; // Table -> ARRAY
"""
@reql_one(77, index_list,
  ReqlTable,
  ReqlArray)

"""
Gets information about whether or not a set of indexes are ready to
be accessed. Returns a list of objects that look like this:
{index:STRING, ready:BOOL[, progress:NUMBER]}

INDEX_STATUS = 139; // Table, STRING... -> ARRAY
"""
@reql_one_twoarr(139, index_status,
  ReqlTable, ReqlString,
  ReqlArray)

"""
Blocks until a set of indexes are ready to be accessed. Returns the
same values INDEX_STATUS.

INDEX_WAIT = 140; // Table, STRING... -> ARRAY
"""
@reql_one_twoarr(140, index_wait,
  ReqlTable, ReqlString,
  ReqlArray)

"""
Renames the given index to a new name

INDEX_RENAME = 156; // Table, STRING, STRING, {overwrite:BOOL} -> OBJECT
"""
@reql_one_two_three(156, index_rename,
  ReqlTable, ReqlString, ReqlString,
  ReqlObject)

# Control Operators

#=
TODO
Calls a function on data

FUNCALL = 64; // Function(*), DATUM... -> DATUM
=#

"""
Executes its first argument, and returns its second argument if it
got [true] or its third argument if it got [false] (like an `if`
statement).

BRANCH = 65; // BOOL, Top, Top -> Top
"""
@reql_one_two_three(65, branch,
  ReqlBool, ReqlTop, ReqlTop,
  ReqlTop)

"""
Returns true if any of its arguments returns true (short-circuits).

OR = 66; // BOOL... -> BOOL
"""
@reql_onearr(66, or,
  ReqlBool,
  ReqlBool)

"""
Returns true if all of its arguments return true (short-circuits).

AND = 67; // BOOL... -> BOOL
"""
@reql_onearr(67, and,
  ReqlBool,
  ReqlBool)

#=
TODO
Calls its Function with each entry in the sequence
and executes the array of terms that Function returns.

FOR_EACH = 68; // Sequence, Function(1) -> OBJECT
=#

# Special Terms

"""
An anonymous function.  Takes an array of numbers representing
variables (see [VAR] above), and a [Term] to execute with those in
scope.  Returns a function that may be passed an array of arguments,
then executes the Term with those bound to the variable names.  The
user will never construct this directly.  We use it internally for
things like `map` which take a function.  The "arity" of a [Function] is
the number of arguments it takes.
For example, here's what `_X_.map{|x| x+2}` turns into:

```
Term {
   type = MAP;
   args = [_X_,
           Term {
             type = Function;
             args = [Term {
                       type = DATUM;
                       datum = Datum {
                         type = R_ARRAY;
                         r_array = [Datum { type = R_NUM; r_num = 1; }];
                       };
                     },
                     Term {
                       type = ADD;
                       args = [Term {
                                 type = VAR;
                                 args = [Term {
                                           type = DATUM;
                                           datum = Datum { type = R_NUM;
                                                           r_num = 1};
                                         }];
                               },
                               Term {
                                 type = DATUM;
                                 datum = Datum { type = R_NUM; r_num = 2; };
                               }];
                     }];
           }];
```

FUNC = 69; // ARRAY, Top -> ARRAY -> Top
"""
@reql_one_two(69, func,
  ReqlArray, ReqlTop,
  ReqlTop)

"""
Indicates to ORDER_BY that this attribute is to be sorted in ascending order.

ASC = 73; // !STRING -> Ordering
"""
@reql_one(73, asc,
  ReqlNumber,
  ReqlOrdering)

"""
Indicates to ORDER_BY that this attribute is to be sorted in descending order.

DESC = 74; // !STRING -> Ordering
"""
@reql_one(74, desc,
  ReqlNumber,
  ReqlOrdering)

"""
Gets info about anything.  INFO is most commonly called on tables.

INFO = 79; // Top -> OBJECT
"""
@reql_one(79, info,
  ReqlTop,
  ReqlObject)

"""
`a.match(b)` returns a match object if the string `a`
matches the regular expression `b`.

MATCH = 97; // STRING, STRING -> DATUM
"""
@reql_one_two(97, match,
  ReqlString, ReqlString,
  ReqlDatum)

"""
Change the case of a string.

UPCASE   = 141; // STRING -> STRING
"""
@reql_one(141, upcase,
  ReqlString,
  ReqlString)

"""
Change the case of a string.

DOWNCASE = 142; // STRING -> STRING
"""
@reql_one(142, downcase,
  ReqlString,
  ReqlString)

"""
Select a number of elements from sequence with uniform distribution.

SAMPLE = 81; // Sequence, NUMBER -> Sequence
"""
@reql_one_two(81, sample,
  ReqlSequence, ReqlNumber,
  ReqlSequence)

"""
Evaluates its first argument.  If that argument returns
NULL or throws an error related to the absence of an
expected value (for instance, accessing a non-existent
field or adding NULL to an integer), DEFAULT will either
return its second argument or execute it if it's a
function.  If the second argument is a function, it will be
passed either the text of the error or NULL as its
argument.

DEFAULT = 92; // Top, Top -> Top
"""
@reql_one_two(92, default,
  ReqlTop, ReqlTop,
  ReqlTop)

"""
Parses its first argument as a json string and returns it as a
datum.

JSON = 98; // STRING -> DATUM
"""
@reql_one(98, json,
  ReqlString,
  ReqlDatum)

"""
Returns the datum as a JSON string.
N.B.: we would really prefer this be named TO_JSON and that exists as
an alias in Python and JavaScript drivers; however it conflicts with the
standard `to_json` method defined by Ruby's standard json library.
TO_JSON_STRING = 172; // DATUM -> STRING
"""
@reql_one(172, to_json_string,
  ReqlDatum,
  ReqlString)

"""
Parses its first arguments as an ISO 8601 time and returns it as a
datum.

ISO8601 = 99; // STRING -> PSEUDOTYPE(TIME)
"""
@reql_one(99, iso8601,
  ReqlString,
  ReqlTime)

"""
Prints a time as an ISO 8601 time.

TO_ISO8601 = 100; // PSEUDOTYPE(TIME) -> STRING
"""
@reql_one(100, to_iso8601,
  ReqlTime,
  ReqlString)

"""
Returns a time given seconds since epoch in UTC.

EPOCH_TIME = 101; // NUMBER -> PSEUDOTYPE(TIME)
"""
@reql_one(101, epoch_time,
  ReqlNumber,
  ReqlTime)

"""
Returns seconds since epoch in UTC given a time.

TO_EPOCH_TIME = 102; // PSEUDOTYPE(TIME) -> NUMBER
"""
@reql_one(102, to_epoch_time,
  ReqlTime,
  ReqlNumber)

"""
The time the query was received by the server.

NOW = 103; // -> PSEUDOTYPE(TIME)
"""
@reql_zero(103, now,
  ReqlTime)

"""
Puts a time into an ISO 8601 timezone.

IN_TIMEZONE = 104; // PSEUDOTYPE(TIME), STRING -> PSEUDOTYPE(TIME)
"""
@reql_one_two(104, in_timezone,
  ReqlTime, ReqlString,
  ReqlTime)

"""
a.during(b, c) returns whether a is in the range [b, c)

DURING = 105; // PSEUDOTYPE(TIME), PSEUDOTYPE(TIME), PSEUDOTYPE(TIME) -> BOOL
"""
@reql_one_two_three(105, during,
  ReqlTime, ReqlTime, ReqlTime,
  ReqlBool)

"""
Retrieves the date portion of a time.

DATE = 106; // PSEUDOTYPE(TIME) -> PSEUDOTYPE(TIME)
"""
@reql_one(106, date,
  ReqlTime,
  ReqlTime)

"""
x.time_of_day == x.date - x

TIME_OF_DAY = 126; // PSEUDOTYPE(TIME) -> NUMBER
"""
@reql_one(126, time_of_day,
  ReqlTime,
  ReqlNumber)

"""
Returns the timezone of a time.

TIMEZONE = 127; // PSEUDOTYPE(TIME) -> STRING
"""
@reql_one(127, timezone,
  ReqlTime,
  ReqlString)

# These access the various components of a time.

"YEAR = 128; // PSEUDOTYPE(TIME) -> NUMBER"
@reql_one(128, year,
  ReqlTime,
  ReqlNumber)

"MONTH = 129; // PSEUDOTYPE(TIME) -> NUMBER"
@reql_one(129, month,
  ReqlTime,
  ReqlNumber)

"DAY = 130; // PSEUDOTYPE(TIME) -> NUMBER"
@reql_one(130, day,
  ReqlTime,
  ReqlNumber)

"DAY_OF_WEEK = 131; // PSEUDOTYPE(TIME) -> NUMBER"
@reql_one(131, day_of_week,
  ReqlTime,
  ReqlNumber)

"DAY_OF_YEAR = 132; // PSEUDOTYPE(TIME) -> NUMBER"
@reql_one(132, day_of_year,
  ReqlTime,
  ReqlNumber)

"HOURS = 133; // PSEUDOTYPE(TIME) -> NUMBER"
@reql_one(133, hours,
  ReqlTime,
  ReqlNumber)

"MINUTES = 134; // PSEUDOTYPE(TIME) -> NUMBER"
@reql_one(134, minutes,
  ReqlTime,
  ReqlNumber)

"SECONDS = 135; // PSEUDOTYPE(TIME) -> NUMBER"
@reql_one(135, seconds,
  ReqlTime,
  ReqlNumber)

#=
TODO
Construct a time from a date and optional timezone or a
date+time and optional timezone.

TIME = 136; // NUMBER, NUMBER, NUMBER, STRING -> PSEUDOTYPE(TIME) |
            // NUMBER, NUMBER, NUMBER, NUMBER, NUMBER, NUMBER, STRING -> PSEUDOTYPE(TIME) |
=#

# Constants for ISO 8601 days of the week.

"MONDAY = 107;    // -> 1"
@reql_zero(107, monday,
  ReqlNumber)

"TUESDAY = 108;   // -> 2"
@reql_zero(108, tuesday,
  ReqlNumber)

"WEDNESDAY = 109; // -> 3"
@reql_zero(109, wednesday,
  ReqlNumber)

"THURSDAY = 110;  // -> 4"
@reql_zero(110, thursday,
  ReqlNumber)

"FRIDAY = 111;    // -> 5"
@reql_zero(111, friday,
  ReqlNumber)

"SATURDAY = 112;  // -> 6"
@reql_zero(112, saturday,
  ReqlNumber)

"SUNDAY = 113;    // -> 7"
@reql_zero(113, sunday,
  ReqlNumber)

# Constants for ISO 8601 months.

"JANUARY = 114;   // -> 1"
@reql_zero(114, january,
  ReqlNumber)

"FEBRUARY = 115;  // -> 2"
@reql_zero(115, february,
  ReqlNumber)

"MARCH = 116;     // -> 3"
@reql_zero(116, march,
  ReqlNumber)

"APRIL = 117;     // -> 4"
@reql_zero(117, april,
  ReqlNumber)

"MAY = 118;       // -> 5"
@reql_zero(118, may,
  ReqlNumber)

"JUNE = 119;      // -> 6"
@reql_zero(119, june,
  ReqlNumber)

"JULY = 120;      // -> 7"
@reql_zero(120, july,
  ReqlNumber)

"AUGUST = 121;    // -> 8"
@reql_zero(121, august,
  ReqlNumber)

"SEPTEMBER = 122; // -> 9"
@reql_zero(122, september,
  ReqlNumber)

"OCTOBER = 123;   // -> 10"
@reql_zero(123, october,
  ReqlNumber)

"NOVEMBER = 124;  // -> 11"
@reql_zero(124, november,
  ReqlNumber)

"DECEMBER = 125;  // -> 12"
@reql_zero(125, december,
  ReqlNumber)

#=
TODO
Indicates to MERGE to replace, or remove in case of an empty literal, the
other object rather than merge it.

LITERAL = 137; // -> Merging
               // JSON -> Merging
=#

#=
"GROUP = 144; // SEQUENCE, STRING -> GROUPED_SEQUENCE | SEQUENCE, FUNCTION -> GROUPED_SEQUENCE"
@reql_one_two(144, group,
  ReqlSequence, ReqlString,
  ReqlSequence)

"SUM = 145; // SEQUENCE, STRING -> GROUPED_SEQUENCE | SEQUENCE, FUNCTION -> GROUPED_SEQUENCE"
@reql_one_two(145, sum, ReqlArray, ReqlString)

"AVG = 146; // SEQUENCE, STRING -> GROUPED_SEQUENCE | SEQUENCE, FUNCTION -> GROUPED_SEQUENCE"
@reql_one_two(146, avg, ReqlArray, ReqlString)

"MIN = 147; // SEQUENCE, STRING -> GROUPED_SEQUENCE | SEQUENCE, FUNCTION -> GROUPED_SEQUENCE"
@reql_one_two(147, min, ReqlArray, ReqlString)

"MAX = 148; // SEQUENCE, STRING -> GROUPED_SEQUENCE | SEQUENCE, FUNCTION -> GROUPED_SEQUENCE"
@reql_one_two(148, max, ReqlArray, ReqlString)
=#

"""
`str.split()` splits on whitespace
`str.split(" ")` splits on spaces only
`str.split(" ", 5)` splits on spaces with at most 5 results
`str.split(nil, 5)` splits on whitespace with at most 5 results

SPLIT = 149; // STRING -> ARRAY | STRING, STRING -> ARRAY | STRING, STRING, NUMBER -> ARRAY | STRING, NULL, NUMBER -> ARRAY
"""
@reql_one(149, split,
  ReqlString,
  ReqlArray)
@reql_one_two(149, split,
  ReqlString, ReqlString,
  ReqlArray)
@reql_one_two_three(149, split,
  ReqlString, ReqlString, ReqlNumber,
  ReqlArray)
@reql_one_two_three(149, split,
  ReqlString, ReqlNull, ReqlNumber,
  ReqlArray)

#=
TODO
UNGROUP = 150; // GROUPED_DATA -> ARRAY
=#

"""
Takes a range of numbers and returns a random number within the range

RANDOM = 151; // NUMBER, NUMBER {float:BOOL} -> DATUM
"""
@reql_one_two(151, random,
  ReqlNumber, ReqlNumber,
  ReqlDatum)

#=
"CHANGES = 152; // TABLE -> STREAM"
@reql_one(152, changes, ReqlTerm)
=#

#=
"ARGS = 154; // ARRAY -> SPECIAL (used to splice arguments)"
@reql_one(154, args, ReqlArray)
=#

# BINARY is client-only at the moment, it is not supported on the server
# BINARY = 155; // STRING -> PSEUDOTYPE(BINARY)

#=
"GEOJSON = 157;           // OBJECT -> PSEUDOTYPE(GEOMETRY)"
@reql_one(157, geojson, ReqlObject)

"TO_GEOJSON = 158;        // PSEUDOTYPE(GEOMETRY) -> OBJECT"
@reql_one(158, to_geojson, ReqlTerm)

"POINT = 159;             // NUMBER, NUMBER -> PSEUDOTYPE(GEOMETRY)"
@reql_one_two(159, point, ReqlNumber, ReqlNumber)
=#

#=
TODO
LINE = 160;              // (ARRAY | PSEUDOTYPE(GEOMETRY))... -> PSEUDOTYPE(GEOMETRY)
=#

#=
TODO
POLYGON = 161;           // (ARRAY | PSEUDOTYPE(GEOMETRY))... -> PSEUDOTYPE(GEOMETRY)
=#

#=
"DISTANCE = 162;          // PSEUDOTYPE(GEOMETRY), PSEUDOTYPE(GEOMETRY) {geo_system:STRING, unit:STRING} -> NUMBER"
@reql_one_two(162, distance, ReqlTerm, ReqlTerm)

"INTERSECTS = 163;        // PSEUDOTYPE(GEOMETRY), PSEUDOTYPE(GEOMETRY) -> BOOL"
@reql_one_two(163, intersects, ReqlTerm, ReqlTerm)

"INCLUDES = 164;          // PSEUDOTYPE(GEOMETRY), PSEUDOTYPE(GEOMETRY) -> BOOL"
@reql_one_two(164, includes, ReqlTerm, ReqlTerm)

"CIRCLE = 165;            // PSEUDOTYPE(GEOMETRY), NUMBER {num_vertices:NUMBER, geo_system:STRING, unit:STRING, fill:BOOL} -> PSEUDOTYPE(GEOMETRY)"
@reql_one(165, distance, ReqlTerm)

"GET_INTERSECTING = 166;  // TABLE, PSEUDOTYPE(GEOMETRY) {index:!STRING} -> StreamSelection"
@reql_one_two(166, get_intersecting, ReqlTerm, ReqlTerm)

"FILL = 167;              // PSEUDOTYPE(GEOMETRY) -> PSEUDOTYPE(GEOMETRY)"
@reql_one(167, fill, ReqlTerm)

"GET_NEAREST = 168;       // TABLE, PSEUDOTYPE(GEOMETRY) {index:!STRING, max_results:NUM, max_dist:NUM, geo_system:STRING, unit:STRING} -> ARRAY"
@reql_one_two(168, get_nearest, ReqlTerm, ReqlTerm)

"POLYGON_SUB = 171;       // PSEUDOTYPE(GEOMETRY), PSEUDOTYPE(GEOMETRY) -> PSEUDOTYPE(GEOMETRY)"
@reql_one_two(171, polygon_sub, ReqlTerm, ReqlTerm)

# Constants for specifying key ranges

"MINVAL = 180;"
@reql_zero(180, minval)

"MAXVAL = 181;"
@reql_zero(181, maxval)
=#

function row(s::ByteString)
  implicit_var() |>
    d -> bracket(d, s)
end

@reql_one_two(69, func,
  ReqlArray, ReqlTop,
  ReqlTop)

function func0(f)
  o = func(make_array(), f)
  ReqlFunction1(-1, o)
end

max_stars_repo_name: mrcsparker/julia-rethinkdb

max_issues_repo_name: mrcsparker/julia-rethinkdb

max_forks_repo_name: mrcsparker/julia-rethinkdb


text:
# Set up function
function _set_up()
    path = tempdir()
    one_file  = joinpath(path, "one.tmp")
    two_file  = joinpath(path, "two.tmp")
    list_file = joinpath(path, "list.tmp")
    write(one_file,  "ONE")
    write(two_file,  "TWO")
    write(list_file, one_file, "\n", two_file, "\n")
    one_file, two_file, list_file
end

# Clean up functions
function _delete_file(path)
    if isfile(path)
        rm(path)
    end
end

function _clean_up_outputs()
    path = tempdir()
    _delete_file(joinpath(path, "one.mitos.tmp"))
    _delete_file(joinpath(path, "two.mitos.tmp"))
end

function _clean_up()
    path = tempdir()
    _delete_file(joinpath(path, "list.tmp"))
    _delete_file(joinpath(path, "one.tmp"))
    _delete_file(joinpath(path, "two.tmp"))
end

# Tests
@testset "Template" begin

    # julia bin
    julia = joinpath(Base.Sys.BINDIR, "julia")
    # ../../
    mitos_folder = splitdir(splitdir(dirname(@__FILE__))[1])[1]

    template = joinpath(mitos_folder, "scripts", "_template.jl")
    file_one, file_two, list_file = _set_up()
    one_out = joinpath(tempdir(), "one.mitos.tmp")
    two_out = joinpath(tempdir(), "two.mitos.tmp")

    @testset "STDIN -> STDOUT" begin
        out = read(pipeline(`$julia -e 'print("Hello hello")'`, `$julia $template`), String)

        @test occursin(r"RUN : 0", out)
        @test length(collect((m.match for m = eachmatch(r"[Hh]ello", out)))) == 2

        @static if Sys.isunix()

            @testset "--arg" begin
                out = read(pipeline(`cat $list_file`, `$julia $template --arg 42`), String)

                @test occursin(r"RUN : 42", out)
                @test length(collect((m.match for m = eachmatch(r"\.tmp", out)))) == 2
            end

            @testset "--list & -a" begin
                _clean_up_outputs()
                out = read(pipeline(`cat $list_file`, `$julia $template -a 42 --list`), String)

                @test occursin(r"ONE", out) # Printed into STDOUT
                @test occursin(r"TWO", out) # Printed into STDOUT
                @test !occursin(r"RUN : 42", out) # Printed into the file
                @test filesize(one_out) != 0 # Created file
                @test filesize(two_out) != 0 # Created file
                @test occursin(r"RUN : 42", read(one_out, String)) # Printed into the file
                @test occursin(r"RUN : 42", read(two_out, String)) # Printed into the file

                _clean_up_outputs()
            end
        end
    end

    @testset "File -> File" begin

        @testset "--list & -a" begin
            _clean_up_outputs()
            out = read(`$julia $template $list_file --a 42 --list`, String)

            @test occursin(r"ONE", out) # Printed into STDOUT
            @test occursin(r"TWO", out) # Printed into STDOUT
            @test !occursin(r"RUN : 42", out) # Printed into the file
            @test filesize(one_out) != 0 # Created file
            @test filesize(two_out) != 0 # Created file
            @test occursin(r"RUN : 42", read(one_out, String)) # Printed into the file
            @test occursin(r"RUN : 42", read(two_out, String)) # Printed into the file

            _clean_up_outputs()
        end
    end

    @testset "-p (parallel)" begin

        @testset "File -> File (--list)" begin
            _clean_up_outputs()
            out = read(`$julia $template $list_file -p 2 --a 42 --list`, String)

            @test occursin(r"ONE", out) # Printed into STDOUT
            @test occursin(r"TWO", out) # Printed into STDOUT
            @test !occursin(r"RUN : 42", out) # Printed into the file
            @test filesize(one_out) != 0 # Created file
            @test filesize(two_out) != 0 # Created file
            @test occursin(r"RUN : 42", read(one_out, String)) # Printed into the file
            @test occursin(r"RUN : 42", read(two_out, String)) # Printed into the file

            _clean_up_outputs()
        end

        @static if Sys.isunix()

            @testset "STDIN -> File (--list)" begin
                _clean_up_outputs()
                out = read(pipeline(`cat $list_file`,`$julia $template -p 2 --a 42 --list`), String)

                @test occursin(r"ONE", out) # Printed into STDOUT
                @test occursin(r"TWO", out) # Printed into STDOUT
                @test !occursin(r"RUN : 42", out) # Printed into the file
                @test filesize(one_out) != 0 # Created file
                @test filesize(two_out) != 0 # Created file
                @test occursin(r"RUN : 42", read(one_out, String)) # Printed into the file
                @test occursin(r"RUN : 42", read(two_out, String)) # Printed into the file
                _clean_up_outputs()
            end
        end
    end

    _clean_up()
end


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ d2c27263-8b8e-486e-9f85-a639ce3749c7
using MD5, PlutoTest

# ╔═╡ 0a80eec0-af92-11eb-2306-a1148fe38e04
md"""
# [Day 4: The Ideal Stocking Stuffer](https://adventofcode.com/2015/day/4)
Santa needs help [mining](https://en.wikipedia.org/wiki/Bitcoin#Mining) some AdventCoins (very similar to [bitcoins](https://en.wikipedia.org/wiki/Bitcoin)) to use as gifts for all the economically forward-thinking little girls and boys.

To do this, he needs to find [MD5](https://en.wikipedia.org/wiki/MD5) hashes which, in [hexadecimal](https://en.wikipedia.org/wiki/Hexadecimal), start with at least five zeroes. The input to the MD5 hash is some secret key (your puzzle input, given below) followed by a number in decimal. To mine AdventCoins, **you must find Santa the lowest positive number (no leading zeroes: `1`, `2`, `3`, ...) that produces such a hash**.

For example:

-   If your secret key is `abcdef`, the answer is `609043`, because the MD5 hash of `abcdef609043` starts with five zeroes (`000001dbbfa...`), and it is the lowest such number to do so.
-   If your secret key is `pqrstuv`, the lowest number it combines with to make an MD5 hash starting with five zeroes is `1048970`; that is, the MD5 hash of `pqrstuv1048970` looks like `000006136ef...`.
"""

# ╔═╡ 48e941ba-480f-4f31-b94c-b004130aa114
puzzle_input = "iwrupvqb"

# ╔═╡ 09b4af88-3d43-4136-8904-12b4f6aa3625
function answer(input, zero_count)
	n = 1
	search = "0" ^ zero_count
	while true
		hash = bytes2hex(md5(input * string(n)))
		if startswith(hash, search)
			return n
		end
		n += 1
	end
end

# ╔═╡ fb6daadc-d6b2-47dd-b7c5-c119d24c33cf
@assert answer("abcdef", 5) == 609043

# ╔═╡ b1b9e6bd-aa83-4c59-b524-189e5bf0d964
@assert answer("pqrstuv", 5) == 1048970

# ╔═╡ e9c559dc-544f-4c28-a9a8-d98477da244e
part1 = answer(puzzle_input, 5)

# ╔═╡ 1fe61914-4666-44dd-9fac-3d96c0484a50
@test part1 == 346386

# ╔═╡ 30832092-32ec-447b-8e66-6a114734629a
md"Your puzzle answer was `346386`."

# ╔═╡ f57a76e1-43db-47ad-971d-89d7a8300db0
md"""
# Part Two
Now find one that starts with six zeroes.
"""

# ╔═╡ 4747ba76-4eac-43ba-bd6a-911e54604c9f
part2 = answer(puzzle_input, 6)

# ╔═╡ 83b80af3-6153-4add-853f-c9d9eaae5268
@test part2 == 9958218

# ╔═╡ 0aa393b8-5b60-4aea-aa0a-025261708ee1
md"Your puzzle answer was `9958218`."

# ╔═╡ Cell order:
# ╠═d2c27263-8b8e-486e-9f85-a639ce3749c7
# ╟─0a80eec0-af92-11eb-2306-a1148fe38e04
# ╟─48e941ba-480f-4f31-b94c-b004130aa114
# ╠═09b4af88-3d43-4136-8904-12b4f6aa3625
# ╠═fb6daadc-d6b2-47dd-b7c5-c119d24c33cf
# ╠═b1b9e6bd-aa83-4c59-b524-189e5bf0d964
# ╠═e9c559dc-544f-4c28-a9a8-d98477da244e
# ╠═1fe61914-4666-44dd-9fac-3d96c0484a50
# ╟─30832092-32ec-447b-8e66-6a114734629a
# ╟─f57a76e1-43db-47ad-971d-89d7a8300db0
# ╠═4747ba76-4eac-43ba-bd6a-911e54604c9f
# ╠═83b80af3-6153-4add-853f-c9d9eaae5268
# ╟─0aa393b8-5b60-4aea-aa0a-025261708ee1


text:
using Dates

export IterationStatistics

"""
A simple structure containing performance stats from the iteration.
"""
mutable struct IterationStatistics
    startTimestamp::DateTime
    endTimestamp::DateTime
    numNodes::Int
    result::String
    additionalInfo::Dict{String, String}
    IterationStatistics(startTimestamp, endTimestamp, numNodes, result, additionalInfo::Dict{String, String} = Dict{String, String}()) = new(startTimestamp, endTimestamp, numNodes, result, additionalInfo)
    IterationStatistics() = new(Dates.now(), Dates.now(), 0, "N/A")
end


text:
# Generate new interactions to perturb ensembles

export
    repeatpocketpoints,
    perturbensemble,
    clusterligsite

"""
Find new interactions arising from fake modulator atoms and form a new `Constraints` object.
Arguments are atoms to calculate constraints with, `Constraints` object to add to, generated modulator
coordinates, distance cutoff and constraint tolerances for intra-modulator and modulator-protein
interactions, and minimum lower distance constraint.
Returns a new `Constraints` object including the modulator.
"""
function Constraints(atoms::Array{Atom,1},
                    old_constraints::Constraints,
                    mod_coords::Array{Float64};
                    intra_cutoff::Real=defaults["mod_intra_cutoff"], # If this is 0.0 there are no intra-modulator interactions
                    intra_tolerance::Real=defaults["mod_intra_tolerance"],
                    inter_cutoff::Real=defaults["mod_inter_cutoff"],
                    inter_tolerance::Real=defaults["mod_inter_tolerance"],
                    min_constraint_dist::Real=defaults["mod_min_constraint_dist"])
    n_atoms = length(atoms)
    @assert n_atoms > 0 "No atoms in atom list"
    constraint_atoms = old_constraints.atoms
    @assert n_atoms == length(constraint_atoms) "Different number of atoms in atom list and constraints"
    n_constraints = size(old_constraints.pres_inds, 1)
    @assert n_constraints == length(old_constraints.lower) && n_constraints == length(old_constraints.upper) "The number of lower and upper constraints do not correspond to the number of constraints present"
    n_mod_atoms = size(mod_coords, 2)
    @assert n_mod_atoms > 0 "No modulator coordinates in list"
    @assert intra_cutoff >= 0 "intra_cutoff cannot be negative"
    @assert intra_tolerance >= 0 "intra_tolerance cannot be negative"
    @assert inter_cutoff >= 0 "inter_cutoff cannot be negative"
    @assert inter_tolerance >= 0 "inter_tolerance cannot be negative"
    @assert min_constraint_dist >= 0 "min_constraint_dist cannot be negative"
    n_total = n_atoms + n_mod_atoms
    # Copy Constraints object and form larger constraints lists
    new_lower = copy(old_constraints.lower)
    new_upper = copy(old_constraints.upper)
    new_pres_inds = copy(old_constraints.pres_inds)
    new_atoms = deepcopy(constraint_atoms)
    # Add new atoms to atom array
    for i in 1:n_mod_atoms
        push!(new_atoms, Atom(mod_atom_info["atom_name"], mod_atom_info["res_name"], mod_atom_info["chain_id"], mod_atom_info["res_n"], mod_coords[:,i], mod_atom_info["element"]))
    end
    new_inds_i = Int[]
    new_inds_j = Int[]
    # Find interactions within modulator
    intra_counter = 0
    intra_cutoff_sq = intra_cutoff^2
    for i in 1:n_mod_atoms
        for j in 1:(i-1)
            sq_dist = 0.0
            for m = 1:3
                sq_dist += abs2(mod_coords[m,i] - mod_coords[m,j])
            end
            if sq_dist < intra_cutoff_sq
                dist = sqrt(sq_dist)
                push!(new_lower, max(dist-intra_tolerance, min_constraint_dist))
                push!(new_upper, dist+intra_tolerance)
                push!(new_inds_i, n_atoms+i)
                push!(new_inds_j, n_atoms+j)
                intra_counter += 1
            end
        end
    end
    # Find interactions between modulator and protein
    inter_counter = 0
    inter_cutoff_sq = inter_cutoff^2
    for i in 1:n_mod_atoms
        for j in 1:n_atoms
            sq_dist = 0.0
            for m = 1:3
                sq_dist += abs2(mod_coords[m,i] - atoms[j].coords[m])
            end
            if sq_dist < inter_cutoff_sq
                dist = sqrt(sq_dist)
                push!(new_lower, max(dist-inter_tolerance, min_constraint_dist))
                push!(new_upper, dist+inter_tolerance)
                push!(new_inds_i, n_atoms+i)
                push!(new_inds_j, j)
                inter_counter += 1
            end
        end
    end
    if length(new_pres_inds) > 0
        new_pres_inds = vcat(new_pres_inds, hcat(new_inds_i, new_inds_j))
    else
        new_pres_inds = hcat(new_inds_i, new_inds_j)
    end
    #println("Found ", intra_counter, " interactions within the modulator")
    println("Found ", inter_counter, " interactions between the modulator and the protein")
    return Constraints(new_atoms, new_lower, new_upper, new_pres_inds)
end

"""
Repeat pocket point to get a certain number of points.
Points are repeated as long as whole copies can be made, then the remainder are chosen randomly.
Returns coordinate array.
"""
function repeatpocketpoints(coords::Array{Float64}, n_out_points::Integer=defaults["mod_n_points"])
    n_pocket_points = Int(length(coords) / 3)
    inds_to_use = repeat(collect(1:n_pocket_points), outer=[Int(floor(n_out_points / n_pocket_points))])
    append!(inds_to_use, rand(1:n_pocket_points, n_out_points % n_pocket_points))
    return getindex(coords, collect(1:3), inds_to_use)
end

"""
Gets ensembles for multiple modulators.
Use pocket points to generate additional constraints.
"""
function perturbensemble(atoms::Array{Atom,1},
                    constraints::Constraints,
                    n_strucs::Integer,
                    mod_path::Union{AbstractString, Nothing},
                    n_mods::Integer)
    ensemble_mods = ModelledEnsemble[]
    if n_mods > 0 && mod_path != nothing
        pock_points = readpocketpoints(mod_path)
        n_pocks = length(pock_points)
        if n_mods > n_pocks
            println(n_mods, " modulators asked for but only ", n_pocks, " found")
        end
        n_mods_to_use = min(n_mods, n_pocks)
        for i in 1:n_mods_to_use
            println("Generating ensemble with modulator ", i, " of ", n_mods_to_use)
            mod_coords = repeatpocketpoints(pock_points[i])
            new_constraints = Constraints(atoms, constraints, mod_coords)
            push!(ensemble_mods, generateensemble(new_constraints, n_strucs))
        end
    end
    return ensemble_mods
end

"""
Assign pocket numbers to LIGSITEcs pocket_r.pdb file points.
Arguments are pocket_r.pdb file, pocket_all.pdb file and new pocket_r.pdb file.
New file has pocket number in place of residue number.
Apears to agree largely but not completely with the LIGSITEcs clustering.
Unassigned points given number 0 and should be fixed manually.
"""
function clusterligsite(point_filepath::AbstractString,
                    centre_filepath::AbstractString,
                    out_filepath::AbstractString;
                    cluster_dist::Real=defaults["cluster_dist"])
    point_lines = readpdblines(point_filepath)
    n_points = length(point_lines)
    point_coords = zeros(3, n_points)
    for (i, line) in enumerate(point_lines)
        point_coords[1,i] = parse(Float64, line[31:38])
        point_coords[2,i] = parse(Float64, line[39:46])
        point_coords[3,i] = parse(Float64, line[47:54])
    end
    println("Read ", n_points, " pocket points from LIGSITEcs pocket points PDB file")
    centre_coords, centre_vols = readligsite(centre_filepath)
    n_centres = length(centre_vols)
    assignments = zeros(Int, n_points)

    # Assign points within cluster_dist Angstroms of a centre to belong to the largest pocket centre
    for p in 1:n_points
        centre_closest = 0
        dist_closest = cluster_dist
        for c in n_centres:-1:1
            dist = norm(point_coords[:,p] - centre_coords[:,c])
            if dist <= cluster_dist
                centre_closest = c
                dist_closest = dist
            end
        end
        assignments[p] = centre_closest
    end

    # While some are still unassigned and have not stalled
    n_unassigned = 2 # Placeholder
    new_n_unassigned = 1 # Placeholder
    while new_n_unassigned != n_unassigned && new_n_unassigned > 0
        n_unassigned = sum(x -> x == 0, assignments)
        for p in 1:n_points
            # For each unassigned point
            if assignments[p] == 0
                point_closest = 0
                dist_closest = cluster_dist
                for a in 1:n_points
                    # Check each assigned point
                    if assignments[a] > 0
                        dist = norm(point_coords[:,p] - point_coords[:,a])
                        # Assign it if it within cluster_dist Angstroms and not already assigned to a larger pocket
                        if dist <= cluster_dist && (point_closest == 0 || assignments[a] < assignments[point_closest])
                            point_closest = a
                            dist_closest = dist
                        end
                    end
                end
                if point_closest > 0
                    assignments[p] = assignments[point_closest]
                end
            end
        end
        new_n_unassigned = sum(x -> x == 0, assignments)
    end
    if new_n_unassigned > 0
        println("Could not assign ", new_n_unassigned, " points")
    else
        println("Assigned all points")
    end

    writeclusterpoints(out_filepath, point_lines, assignments)
end


text:
using Test, Batman

## Ideas for tests
# 1. Simple coverage tests. Try to use all functions
# 2. Benchmarking tests; speed of objective functions
# 3. Bias / pull for known fits.

include("poisson.jl")


text:
# Execute a "benchmarkable" function to estimate the time required to perform
# a single evaluation of the benchmark's core expression. To do this reliably,
# we employ a series of estimates that allow us to decide on a sampling
# strategy and an estimation strategy for performing our benchmark.
#
# Arguments:
#
#     f!::Function: A "benchmarkable" function to be evaluated.
#
#     samples::Integer: The number of samples the user requests. Defaults to
#         100 samples.
#
#     budget::Integer: The time in seconds the user is willing to spend on
#         benchmarking. Defaults to 10 seconds.
#
#     τ::Float64: The minimum R² of the OLS model before the geometric search
#         procedure is considered to have converged. Defaults to 0.95.
#
#     α::Float64: The growth rate for the geometric search. Defaults to 1.1.
#
#     ols_samples::Integer: The number of samples per unique value of `n_evals`
#         when the geometric search procedure is employed. Defaults to 100.
#
#     verbose::Bool: Should the system print our verbose progress information?
#         Defaults to false.
#
# Returns:
#
#      r::Results: A Results object containing information about the
#         benchmark's full execution history.

function execute(
        f!::Function,
        samples::Integer = 100,
        budget::Integer = 10,
        τ::Float64 = 0.95,
        α::Float64 = 1.1,
        ols_samples::Integer = 100,
        verbose::Bool = false,
    )
    # We track the total time spent executing our benchmark.
    start_time = time()

    # Run the "benchmarkable" function once. This can require a compilation
    # step, which might bias the resulting estimates.
    s = Samples()
    f!(s, 1, 1)

    # Determine the elapsed time for the very first call to f!.
    biased_time_ns = s.elapsed_times[1]

    # We stop benchmarking f! if we've already exhausted our time budget.
    time_used = time() - start_time
    if time_used > budget
        return Results(false, false, false, s, time_used)
    end

    # We determine the maximum number of samples we could record given
    # our remaining time budget. We convert to nanoseconds before comparing
    # this number with our estimate of the per-evaluation time cost.
    remaining_time_ns = 10^9 * (budget - time_used)
    max_samples = floor(Integer, remaining_time_ns / biased_time_ns)

    # We stop benchmarking if running f! one more time would put us over
    # our time budget.
    if max_samples < 1
        return Results(false, false, false, s, time_used)
    end

    # Having reached this point, we can afford to record at least one more
    # sample without using up our time budget. The core question now is:
    #
    #     Is the expression being measured so fast that a single sample needs
    #     to evaluate the core expression multiple times before a single sample
    #     can provide an unbiased estimate of the expression's execution time?
    #
    # To determine this, we execute f! one more time. This provides
    # our first potentially unbiased estimate of the execution time, because
    # all compilations costs should now have been paid.

    # Before we execute f!, we empty our biased Samples object to discard
    # the values associated with our first execution of f!.
    empty!(s)

    # We evaluate f! to generate our first potentially unbiased sample.
    f!(s, 1, 1)

    # We can now improve our estimate of the expression's per-evaluation time.
    debiased_time_ns = s.elapsed_times[1]

    # If we've used up our time budget, we stop. We also stop if the user
    # only requested a single sample.
    time_used = time() - start_time
    if time_used > budget || samples == 1
        return Results(true, false, false, s, time_used)
    end

    # Now we determine if the function is so fast that we need to execute the
    # core expression multiple times per sample. We do this by determining if
    # the single-evaluation time is at least 1,000 times larger than the system
    # clock's resolution. If the function is at least that costly to execute,
    # then we determine how many single-evaluation samples we should employ.
    if debiased_time_ns > 1_000 * estimate_clock_resolution()
        remaining_time_ns = 10^9 * (budget - time_used)
        max_samples = floor(Integer, remaining_time_ns / debiased_time_ns)
        n_samples = min(max_samples, samples - 1)
        f!(s, n_samples, 1)
        return Results(true, true, false, s, time() - start_time)
    end

    # If we've reached this far, we are benchmarking a function that is so fast
    # that we need to be careful with our execution strategy. In particular,
    # we need to evaluate the core expression multiple times to generate a
    # single sample. To determine the correct number of times we should
    # evaluate the core expression per sample, we perform a geometric search
    # that starts at 2 evaluations per sample and increases by a factor of 1.1
    # evaluations on each iteration. Having generated data in this form, we
    # use an OLS regression to estimate the per-evaluation timing of our core
    # expression. We stop our geometric search when the OLS linear model is
    # almost perfect fit to our empirical data.

    # We start by executing two evaluations per sample.
    n_evals = 2.0

    # print header about the search progress
    verbose && @printf "%s\t%20s\t%8s\t%s\n" "time_used" "n_evals" "b" "r²"

    # Now we perform a geometric search.
    finished = false
    a, b = NaN, NaN
    while !finished
        # Gather many samples, each of which includes multiple evaluations.
        f!(s, ols_samples, ceil(Integer, n_evals))

        # Perform an OLS regression to estimate the per-evaluation time.
        a, b, r² = ols(s.evaluations, s.elapsed_times)

        # Stop our search when either:
        #  (1) The OLS fit is good enough; or
        #  (2) We've exhausted our time budget.
        time_used = time() - start_time
        min_time = 0.5
        if (r² > τ && time_used > min_time) || time_used > budget
            finished = true
        end

        # We optionally print out information about our search's progress.
        if verbose
            @printf(
                "%.1f\t%24.1f\t%12.2f\t%1.3f\n",
                time_used,
                n_evals,
                b,
                r²
            )
        end

        # We increase the number of evaluations per sample for the next round
        # of our search.
        n_evals *= α
    end

    return Results(true, true, true, s, time() - start_time)
end


text:
"""
    PaduaTransforms

an implementation of the Padua transform and its inverse via the fast Fourier transform.
"""
module PaduaTransforms

using StaticArrays: SVector
using FFTW
using LinearAlgebra: rmul!

export getpaduanum, getdegree, nextpaduanum
export getpaduapoints
export PaduaTransformPlan, paduatransform!
export InvPaduaTransformPlan, invpaduatransform!

## Number of Padua Points and Degree ##

"""
    getpaduanum(n)

calculates number of Padua points needed to approximate a function using Chebyshev polynomials
up to total degree `n`. This number is equal to the number of coefficients. The formula is

```math
N = (n + 1) ⋅ (n + 2) ÷ 2
```

# Examples
```jldoctest
julia> getpaduanum(13)
105
```
"""
getpaduanum(degree) = (degree + 1) * (degree + 2) ÷ 2

"""
    getdegree(N)

calculates total degree, given the number of coefficients or Padua points `N`.
Throws an error if `N` is not a possible number of Padua points. The formula is

```math
n = \\frac{\\sqrt{1 + 8N} - 3}{2}
```

# Examples
```jldoctest
julia> getdegree(105)
13

julia> getdegree(104)
ERROR: ArgumentError: number of Padua points or coeffs must be (n + 1) * (n + 2) ÷ 2
[...]
```
"""
function getdegree(N)
    d = (sqrt(1 + 8N) - 3) / 2
    isinteger(d) ? Int(d) : throw(ArgumentError(
        "number of Padua points or coeffs must be (n + 1) * (n + 2) ÷ 2"))
end

"""
    nextpaduanum(N)

get next valid number of Padua points ≥ `N`.

# Examples

```jldoctest
julia> nextpaduanum(104)
105
```
"""
function nextpaduanum(N)
    d = Int(cld(sqrt(1 + 8N) - 3, 2))
    getpaduanum(d)
end

## Padua Points ##

"""
    paduapoint(T::Type, j::Integer, i::Integer, n::Integer)

returns the Padua point ``z_{ij}``, where

```math
z_{ij} = (\\cos{\\frac{jπ}{n}}, \\cos{\\frac{iπ}{n+1}})
```

Note, that only points with ``i-j`` even are actually Padua points.
Check with [`ispadua`](@ref).

# Examples
```jldoctest
julia> [PaduaTransforms.paduapoint(Float32, x, y, 1) for y in 0:1+1, x in 0:1]
3×2 Matrix{Tuple{Float32, Float32}}:
 (1.0, 1.0)   (-1.0, 1.0)
 (1.0, 0.0)   (-1.0, 0.0)
 (1.0, -1.0)  (-1.0, -1.0)
```

"""
function paduapoint(::Type{T}, j, i, n) where T
    x = cospi(T(j) / T(n))
    y = cospi(T(i) / T(n + 1))
    return x, y
end

"""
    ispadua(i, j)

returns if [`paduapoint`](@ref) at position `(i, j)` is a Padua point.

# Examples
```jldoctest
julia> pointornothing(i, j, n) = PaduaTransforms.ispadua(i, j) ? PaduaTransforms.paduapoint(Float64, j, i, n) : nothing
pointornothing (generic function with 1 method)

julia> [pointornothing(y, x, 2) for y in 0:3, x in 0:2]
4×3 Matrix{Union{Nothing, Tuple{Float64, Float64}}}:
 (1.0, 1.0)   nothing      (-1.0, 1.0)
 nothing      (0.0, 0.5)   nothing
 (1.0, -0.5)  nothing      (-1.0, -0.5)
 nothing      (0.0, -1.0)  nothing
```

"""
ispadua(i, j) = iseven(i - j)

"""
    getpaduapoints([T=Float64,] n)

returns the Padua points

```math
\\textrm{Pad}_n = \\{(\\cos{\\frac{jπ}{n}}, \\cos{\\frac{iπ}{n + 1}}) \\; | \\;
    0 ≤ j ≤ n, \\; 0 ≤ i ≤ n + 1, \\; i - j \\; \\textrm{even} \\}
```

where each row is a point

# Examples

```jldoctest
julia> getpaduapoints(Float32, 1)
3×2 Matrix{Float32}:
  1.0   1.0
  1.0  -1.0
 -1.0   0.0
```
"""
function getpaduapoints(::Type{T}, n) where T
    out = Matrix{T}(undef, getpaduanum(n), 2)

    i = 1
    for x in 0:n
        for y in 0:n+1
            if ispadua(x, y)
                out[i, 1:2] .= paduapoint(T, x, y, n)
                i += 1
            end
        end
    end

    return out
end

getpaduapoints(n) = getpaduapoints(Float64, n)

"""
    getpaduapoints(f::Function, [T=Float64,] n)

evaluates the function `f` on the Padua points for degree `n`. If `f` returns a single value,
`getpaduapoints` returns a `Vector{T}`, else if `f` returns a tuple or other iterable
`getpaduapoints` returns a `Matrix{T}` where each row represents `f` applied to a Padua point.

# Examples
```jldoctest
julia> getpaduapoints(Float64, 2) do x, y; 3x - y, y^2; end
6×2 Matrix{Float64}:
  2.0  1.0
  3.5  0.25
 -0.5  0.25
  1.0  1.0
 -4.0  1.0
 -2.5  0.25
```
"""
function getpaduapoints(f::Function, ::Type{T}, n) where T
    D = length(f(zero(T), zero(T)))
    if D == 1
        out = Vector{T}(undef, getpaduanum(n))
    else
        out = Matrix{T}(undef, getpaduanum(n), D)
    end

    # Function barrier to alleviate issues from type instability
    _fillpoints!(out, f, T, n)

    return out
end

getpaduapoints(f::Function, n) = getpaduapoints(f, Float64, n)

function _fillpoints!(out::AbstractMatrix, f, ::Type{T}, n) where T
    i = 1
    for x in 0:n
        for y in 0:n+1
            if ispadua(x, y)
                v = paduapoint(T, x, y, n)
                out[i, :] .= f(v[1], v[2])
                i += 1
            end
        end
    end

    out
end

function _fillpoints!(out::AbstractVector, f, ::Type{T}, n) where T
    i = 1
    for x in 0:n
        for y in 0:n+1
            if ispadua(x, y)
                v = paduapoint(T, x, y, n)
                out[i] = f(v[1], v[2])
                i += 1
            end
        end
    end

    out
end

## Padua Transform ##

struct PaduaTransformPlan{T, P}
    degree::Int
    vals::Matrix{T}
    dctplan::P
end

"""
    PaduaTransformPlan{T}(n::Integer)

create plan to compute coefficients of Chebyshev polynomials in 2D up to total degree `n`
using the Padua transform.

"""
function PaduaTransformPlan{T}(degree::Integer) where T
    vals = Matrix{T}(undef, degree + 2, degree + 1)
    plan = FFTW.plan_r2r!(vals, FFTW.REDFT00)

    PaduaTransformPlan{T, typeof(plan)}(degree, vals, plan)
end

"""
    weight!(mat::AbstractMatrix, degree::Integer)

weight fourier coefficients to obtain Chebyshev coefficients as part of a [`paduatransform!`](@ref).
The weighting factor applied to the coefficients is

```math
w = \\frac{1}{n(n+1)} ⋅ \\begin{cases}
    \\frac{1}{2} & \\textrm{if on vertex}   \\\\
    1           & \\textrm{if on edge}    \\\\
    2           & \\textrm{if in interior} \\\\
\\end{cases}
```

# Examples
```julia-repl
julia> weight!(ones(4+2, 4+1), 4)
6×5 Matrix{Float64}:
 0.025  0.05  0.05  0.05  0.025
 0.05   0.1   0.1   0.1   0.05
 0.05   0.1   0.1   0.1   0.05
 0.05   0.1   0.1   0.1   0.05
 0.05   0.1   0.1   0.1   0.05
 0.025  0.05  0.05  0.05  0.025
```
"""
function weight!(mat::AbstractMatrix{T}, degree::Integer) where T
    rmul!(mat, T(2 / ( degree * (degree + 1) )))
    rmul!(@view(mat[1, :]), T(0.5))
    rmul!(@view(mat[end, :]), T(0.5))
    rmul!(@view(mat[:, 1]), T(0.5))
    rmul!(@view(mat[:, end]), T(0.5))

    mat
end

"""
    tovalsmat!(mat::Matrix, from::AbstractVector, degree::Integer)

write values of function evaluated at Padua points from `from` to matrix `mat`.

# Examples
```jldoctest
julia> PaduaTransforms.tovalsmat!(ones(3 + 2, 3 + 1), 1:getpaduanum(3), 3)
5×4 Matrix{Float64}:
 1.0  0.0  6.0   0.0
 0.0  4.0  0.0   9.0
 2.0  0.0  7.0   0.0
 0.0  5.0  0.0  10.0
 3.0  0.0  8.0   0.0

julia> PaduaTransforms.tovalsmat!(ones(2 + 2, 2 + 1), 1:getpaduanum(2), 2)
4×3 Matrix{Float64}:
 1.0  0.0  5.0
 0.0  3.0  0.0
 2.0  0.0  6.0
 0.0  4.0  0.0
```
"""
function tovalsmat!(mat::Matrix{T}, from::AbstractVector, degree::Integer) where T
    axes(from, 1) == 1:getpaduanum(degree) || error()
    size(mat) == (degree + 2, degree + 1) || error()

    if isodd(degree)
        # x 0
        # 0 x
        # x 0

        @inbounds for i in 1:length(from)
            mat[2i - 1] = from[i]
            mat[2i] = zero(T)
        end
    else
        @assert iseven(degree)
        # x 0 x
        # 0 x 0
        # x 0 x
        # 0 x 0

        valspercol = (degree + 2) ÷ 2

        # odd columns (j is column index)
        for j in 1:2:degree + 1, i in 1:valspercol
            k = (j - 1) * valspercol + i
            @inbounds mat[2i - 1, j] = from[k]
            @inbounds mat[2i, j] = zero(T)
        end

        # even columns
        for j in 2:2:degree + 1, i in 1:valspercol
            k = (j - 1) * valspercol + i
            @inbounds mat[2i - 1, j] = zero(T)
            @inbounds mat[2i, j] = from[k]
        end
    end

    mat
end

"""
    fromcoeffsmat!(to::AbstractVector, mat::Matrix, degree::Integer, ::Val{lex})

write Chebyshev coefficients from `mat` into vector `to`. `lex::Bool` determines whether
coefficients should be written in lexigraphical order or not. The lower right triangle does
not get written into `to`. These would represent greater polynomial degrees than `degree`.

If `lex` is `Val(true)` the coefficients correspond to the following basis polynomials

```math
T_0(x) T_0(y), T_1(x) T_0(y), T_0(x) T_1(y), T_2(x) T_0(y), T_1(x) T_1(y), T_0(x) T_2(y), ...
```

else if `lex` is `Val(false)` they correspond to

```math
T_0(x) T_0(y), T_0(x) T_1(y), T_1(x) T_0(y), T_0(x) T_2(y), T_1(x) T_1(y), T_2(x) T_0(y), ...
```

# Examples
```julia-repl
julia> mat = [(x, y) for y in 0:2+1, x in 0:2]
4×3 Matrix{Tuple{Int64, Int64}}:
 (0, 0)  (1, 0)  (2, 0)
 (0, 1)  (1, 1)  (2, 1)
 (0, 2)  (1, 2)  (2, 2)
 (0, 3)  (1, 3)  (2, 3)

julia> to1 = similar(mat, getpaduanum(2)); to2 = similar(mat, getpaduanum(2));

julia> fromcoeffsmat!(to1, mat, 2, Val(true))
6-element Vector{Tuple{Int64, Int64}}:
 (0, 0)
 (1, 0)
 (0, 1)
 (2, 0)
 (1, 1)
 (0, 2)

julia> fromcoeffsmat!(to2, mat, 2, Val(false))
6-element Vector{Tuple{Int64, Int64}}:
 (0, 0)
 (0, 1)
 (1, 0)
 (0, 2)
 (1, 1)
 (2, 0)
```
"""
function fromcoeffsmat!(to::AbstractVector, mat::Matrix, degree::Integer, ::Val{false})
    length(to) == getpaduanum(degree) || error()
    axes(mat) == (1:(degree + 2), 1:(degree + 1)) || error()

    n = firstindex(to)
    for d in 1:degree + 1
        for ix in 1:d
            iy = d - ix + 1
            @assert ix + iy == d + 1 "ix and iy must lie on d-th diagonal"

            to[n] = mat[iy, ix]
            n += 1
        end
    end

    to
end

function fromcoeffsmat!(to::AbstractVector, mat::Matrix, degree::Integer, ::Val{true})
    length(to) == getpaduanum(degree) || error()
    size(mat) == (degree + 2, degree + 1) || error()

    n = firstindex(to)
    for d in 1:degree + 1
        for iy in 1:d
            ix = d - iy + 1
            @assert ix + iy == d + 1 "ix and iy must lie on d-th diagonal"

            to[n] = mat[iy, ix]
            n += 1
        end
    end

    to
end

"""
    fromcoeffsmat!(to::AbstractMatrix, mat::Matrix, degree::Integer)

copy Chebyshev coefficients from `mat` to `to` without copying coefficients coresponding to
total degree greater than `degree`.

# Examples
```jldoctest
julia> PaduaTransforms.fromcoeffsmat!(zeros(4, 4), reshape(1:20, 5, 4), 3)
4×4 Matrix{Float64}:
 1.0  6.0  11.0  16.0
 2.0  7.0  12.0   0.0
 3.0  8.0   0.0   0.0
 4.0  0.0   0.0   0.0
```
"""
function fromcoeffsmat!(to::AbstractMatrix, mat::AbstractMatrix, degree::Integer)
    axes(to) == (1:degree + 1, 1:degree + 1) || error()
    axes(mat) == (1:degree + 2, 1:degree + 1) || error()

    for j in 1:degree + 1
        for i in 1:degree + 2 - j
            @inbounds to[i, j] = mat[i, j]
        end
    end

    to
end

"""
    paduatransform!(out, P::PaduaTransformPlan, vals[, lex])

obtain coefficients of Chebyshev polynomials on 2D via the Padua transform, given values
`vals` evaluated at the Padua points. Coefficients will be written into `out`, which should
either be a matrix or a vector.

if `out` is a matrix, make sure that all entries in the lower right diagonal are zero as
these will not get overwritten.

`lex` determines the order in which coefficeints are written into `out` if `out` is a vector.
See [`fromcoeffsmat!`](@ref) for details.

# Examples
```jldoctest
julia> plan = PaduaTransformPlan{Float64}(3);

julia> f(x, y) = 3 + 4x + 5 * y * (2x^2 - 1)
f (generic function with 1 method)

julia> vals = getpaduapoints(f, 3)
10-element Vector{Float64}:
 12.0
  7.0
  2.0
  3.232233047033631
  6.767766952966369
 -1.5000000000000004
  1.0000000000000004
  3.5000000000000013
  2.5355339059327378
 -4.535533905932738

julia> paduatransform!(zeros(4, 4), plan, vals)
4×4 Matrix{Float64}:
  3.0           4.0         0.0  1.4803e-16
 -5.92119e-16  -1.4803e-16  5.0  0.0
  0.0           0.0         0.0  0.0
 -2.96059e-16   0.0         0.0  0.0

julia> paduatransform!(zeros(getpaduanum(3)), plan, vals, Val(true))
10-element Vector{Float64}:
  3.0
  4.0
 -5.921189464667501e-16
  0.0
 -1.4802973661668753e-16
  0.0
  1.4802973661668753e-16
  5.0
  0.0
 -2.9605947323337506e-16

julia> paduatransform!(zeros(getpaduanum(3)), plan, vals, Val(false))
10-element Vector{Float64}:
  3.0
 -5.921189464667501e-16
  4.0
  0.0
 -1.4802973661668753e-16
  0.0
 -2.9605947323337506e-16
  0.0
  5.0
  1.4802973661668753e-16
```
"""
function paduatransform!(P::PaduaTransformPlan)
    coeffs = P.dctplan * P.vals
    weight!(coeffs, P.degree)

    coeffs
end

function paduatransform!(out, P::PaduaTransformPlan, vals, args...)
    tovalsmat!(P.vals, vals, P.degree)
    coeffs = paduatransform!(P)
    fromcoeffsmat!(out, coeffs, P.degree, args...)
end

"""
    paduatransform!(out::AbstractArray{<:Any, 3}, P::PaduaTransformPlan, vals::AbstractMatrix, args...)

transforms each column in `vals` and writes the resulting coefficients in a slice of `out`.
"""
function paduatransform!(out::AbstractArray{<:Any, 3}, P::PaduaTransformPlan, vals::AbstractMatrix, args...)
    axes(out, 3) == axes(vals, 2)|| error()

    @views for i in axes(out, 3)
        paduatransform!(out[:, :, i], P, vals[:, i], args...)
    end

    out
end

"""
    paduatransform!(out::Array{<:Any, 3}, P::PaduaTransformPlan, vals::AbstractVector{<:AbstractVector{T}}, args...)

transforms vector of vectors to Chebyshev coefficients and writes the resulting coefficients
in a slice of `out`. Each vector in the vector of vectors represents a point. Each slice of
`out` represents a transform of one dimension.
"""
function paduatransform!(out::Array{<:Any, 3}, P::PaduaTransformPlan, vals::AbstractVector{<:AbstractVector{T}}, args...) where T
    # Here, each column is a point and each row represents one dimension
    r = reinterpret(reshape, T, vals)
    axes(out, 3) == axes(r, 1) || error()

    @views for i in axes(out, 3)
        paduatransform!(out[:, :, i], P, r[i, :], args...)
    end

    out
end

## Inverse Padua Transform ##

struct InvPaduaTransformPlan{T, P}
    degree::Int
    coeffs::Matrix{T}
    dctplan::P
end

"""
    InvPaduaTransformPlan{T}(n::Integer)

create plan to compute values on Padua points, given coefficients of Chebyshev polynomials
up to total degree `n`.
"""
function InvPaduaTransformPlan{T}(degree::Integer) where T
    coeffs = Matrix{T}(undef, degree + 2, degree + 1)
    iplan = FFTW.plan_r2r!(coeffs, FFTW.REDFT00)

    InvPaduaTransformPlan{T, typeof(iplan)}(degree, coeffs, iplan)
end

"""
    tocoeffsmat!(mat::AbstractMatrix, coeffs::AbstractMatrix)

writes coefficients in `coeffs` into matrix `mat` for the [`invpaduatransform!`](@ref).

# Examples
```jldoctest
julia> PaduaTransforms.tocoeffsmat!(zeros(5, 4), reshape(1:16, 4, 4))
5×4 Matrix{Float64}:
 1.0  5.0   9.0  13.0
 2.0  6.0  10.0  14.0
 3.0  7.0  11.0  15.0
 4.0  8.0  12.0  16.0
 0.0  0.0   0.0   0.0
```
"""
function tocoeffsmat!(mat::AbstractMatrix{T}, coeffs::AbstractMatrix) where T
    mat[1:end-1, :]  = coeffs
    mat[    end, :] .= zero(T)

    mat
end

"""
    invweight!(coeffs::AbstractMatrix)

weight Chebyshev coefficients before the Fourier transform for the [`invpaduatransform!`](@ref).
using the weighting

```math
w = \\begin{cases}
    1            & \\textrm{if on vertex}   \\\\
    \\frac{1}{2} & \\textrm{if on edge}     \\\\
    \\frac{1}{4} & \\textrm{if in interior} \\\\
\\end{cases}
```

# Examples
```jldoctest
julia> PaduaTransforms.invweight!(ones(5, 5))
5×5 Matrix{Float64}:
 1.0  0.5   0.5   0.5   1.0
 0.5  0.25  0.25  0.25  0.5
 0.5  0.25  0.25  0.25  0.5
 0.5  0.25  0.25  0.25  0.5
 1.0  0.5   0.5   0.5   1.0
```
"""
function invweight!(coeffs::AbstractMatrix{T}) where T
    rmul!(@view(coeffs[:,2:end-1]), T(0.5))
    rmul!(@view(coeffs[2:end-1, :]), T(0.5))

    coeffs
end

"""
    fromvalsmat!(to::AbstractVector, mat::AbstractMatrix, n::Integer)

write values from `mat` into the vector `to` after an [`invpaduatransform!`](@ref) of total
degree `n`.
"""
function fromvalsmat!(to::AbstractVector, mat::AbstractMatrix, degree::Integer)
    axes(to, 1) == 1:getpaduanum(degree) || error()
    axes(mat) == (1:degree + 2, 1:degree + 1) || error()

    if isodd(degree)
        # x 0
        # 0 x
        # x 0

        @inbounds for i in 1:length(to)
            to[i] = mat[2i - 1]
        end
    else
        @assert iseven(degree)
        # x 0 x
        # 0 x 0
        # x 0 x
        # 0 x 0

        valspercol = (degree + 2) ÷ 2

        # odd columns (j is column index)
        for j in 1:2:degree + 1, i in 1:valspercol
            k = (j - 1) * valspercol + i
            @inbounds to[k] = mat[2i - 1, j]
        end

        # even columns
        for j in 2:2:degree + 1, i in 1:valspercol
            k = (j - 1) * valspercol + i
            @inbounds to[k] = mat[2i, j]
        end
    end

    to
end

function invpaduatransform!(IP::InvPaduaTransformPlan)
    invweight!(IP.coeffs)
    IP.dctplan * IP.coeffs

    IP.coeffs
end

"""
    invpaduatransform!(vals::AbstractVector, IP::InvPaduaTransformPlan, coeffs::AbstractMatrix)

evaluates the polynomial defined by the coefficients of Chebyshev polynomials `coeffs` on the
Padua points using the inverse transform plan `IP` and writes the resulting values into `vals`.
"""
function invpaduatransform!(vals::AbstractVector, IP::InvPaduaTransformPlan, coeffs::AbstractMatrix)
    tocoeffsmat!(IP.coeffs, coeffs)
    invpaduatransform!(IP)
    fromvalsmat!(vals, IP.coeffs, IP.degree)
end

function invpaduatransform!(vals::AbstractMatrix, IP::InvPaduaTransformPlan, coeffs::AbstractArray{<:Any, 3})
    axes(coeffs, 3) == axes(vals, 2) || error()

    @views for i in axes(coeffs, 3)
        invpaduatransform!(vals[:, i], IP, coeffs[:, :, i])
    end

    vals
end

end  # Padua


text:
function fib(n)
  x,y = (0,1)
  for i = 1:n x,y = (y, x+y) end
  x
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export NLfunForward_test1

"""
weighted forward function
"""
function NLfunForward_test1(A, x, g, params)

    e = params["numr"]*params["nr"]
    w = params["w"]

    e1 = params["numr"]*params["nr1"]
    W_data = params["weight"]
    U_t = copy(W_data[1:e1,1])
    V_t = copy(W_data[e1+1:end,1])
    U = reshape(copy(U_t),params["numr"], params["nr1"])
    V = reshape(copy(V_t),params["numc"], params["nr1"])

    PU = U*U'
    PV = V*V'
    m_U::Int,n_U::Int = size(PU)
    m_V::Int,n_V::Int = size(PV)

    PU_org = Matrix{Complex{Float64}}(I,m_U,n_U)-copy(PU)
    PV_org = Matrix{Complex{Float64}}(I,m_V,n_V)-copy(PV)

    QU = (1/w)*PU+PU_org
    QV = (1/w)*PV+PV_org

    L = x[1:e]
    R = x[(e+1):end]
    L = reshape(L, params["numr"], params["nr"])
    R = reshape(R, params["numc"], params["nr"])

    if isempty(g)
        f1 = params["afun"](QU*L*(QV*R)', params)
        f2 = 0.
    else
        fp = params["afunT"](g)
        f1 = [vec(QU*fp*(QV*R)); vec(QV*fp'*(QU*L))]
        f2 = vec(QU*fp*QV)
    end

    return f1,f2
end


text:
using CompScienceMeshes, BEAST, Base.Test

D, Δx = 1.0, 0.35
fn = joinpath(dirname(@__FILE__),"assets","sphere35.in")
#Γ = meshsphere(D, Δx)
Γ = readmesh(fn)
X = raviartthomas(Γ)
Y = subset(X,[1])

o, x, y, z = euclidianbasis(3)
direction, polarisation = z, x
Nt = 200

sol1 = 1.0
Δt1= 0.08/sol1
U1 = timebasisdelta(Δt1, Nt)
W1 = Y ⊗ U1
duration1, delay1, amplitude1 = 8.0/sol1, 12.0/sol1, 1.0
gaussian1 = creategaussian(duration1, delay1, duration1)
E1 = planewave(polarisation, direction, derive(gaussian1), sol1)
b1 = assemble(E1,W1)
taxis1 = collect((0:Nt-1)*Δt1)

sol2 = 36.0
Δt2 = 0.08/sol2
U2 = timebasisdelta(Δt2, Nt)
W2 = Y ⊗ U2
duration2, delay2, amplitude2 = 8.0/sol2, 12.0/sol2, 1.0
gaussian2 = creategaussian(duration2, delay2, duration2)
E2 = planewave(polarisation, direction, derive(gaussian2), sol2)
b2 = assemble(E2,W2)
taxis2 = collect((0:Nt-1)*Δt2)

gaussian1.(taxis1)[120]
gaussian2.(taxis2)[120]
@test all(gaussian1.(taxis1) .≈ gaussian2.(taxis2))

derive(gaussian1).(taxis1)[120]
derive(gaussian2).(taxis2)[120]/sol2
@test all(derive(gaussian1).(taxis1)/sol1 .≈ derive(gaussian2).(taxis2)/sol2)

timeels1, timead1 = BEAST.assemblydata(U1)
@test timeels1[100][1][1] ≈ Δt1*100
@test timeels1[100][2][1] ≈ Δt1*101

timeels2, timead2 = BEAST.assemblydata(U2)
@test timeels2[100][1][1] ≈ Δt2*100
@test timeels2[100][2][1] ≈ Δt2*101

b1[120]/sol1
b2[120]/sol2
@test all(b1/sol1 .≈ b2/sol2)


text:
abstract type Frame end

abstract type BodyFrame <: Frame end
struct Body <: BodyFrame end

abstract type WindFrame <: Frame end
struct Wind <: WindFrame end

abstract type GlobalFrame <: Frame end
struct Global <: GlobalFrame end

# function rotation_matrix(::Body, ::Global, model, X, freestream)

# end

@inline function R_from_wind(::Body, α, θ)
    R = ST.@SArray [-cos(α) sin(α); sin(α) cos(α)]
end

@inline function R_from_wind(::Wind, α, θ)
    R = ST.@SArray [1 0; 0 1]
end

@inline function R_from_wind(::Global, α, θ)
    δ = θ - α
    R = ST.@SArray [-cos(δ) -sin(δ); -sin(δ) cos(δ)]
end

@inline function R_from_body(::Body, α, θ)
    R = ST.@SArray [1 0; 0 1]
end

@inline function R_from_body(::Wind, α, θ)
    R = ST.@SArray [-cos(α) sin(α); sin(α) cos(α)]
end

@inline function R_from_body(::Global, α, θ)
    R = ST.@SArray [cos(θ) -sin(θ); sin(θ) cos(θ)]
end

@inline function R_from_global(::Body, α, θ)
    R = ST.@SArray [cos(θ) sin(θ); -sin(θ) cos(θ)]
end

@inline function R_from_global(::Wind, α, θ)
    δ = θ - α
    R = ST.@SArray [-cos(δ) -sin(δ); -sin(δ) cos(δ)]
end

@inline function R_from_global(::Global, α, θ)
    R = ST.@SArray [1 0; 0 1]
end

max_stars_repo_name: byuflowlab/MultifidelityTrajectories

max_issues_repo_name: byuflowlab/MultifidelityTrajectories

max_forks_repo_name: byuflowlab/MultifidelityTrajectories


text:
module TriDiagSolvers

include("TriDiagSolverFuncs.jl")

export solve_tdma!
export solve_tridiag_wrapper!

using LinearAlgebra
using ..FiniteDifferenceGrids
using ..StateVecs

struct DifussionImplicit end
function TridiagonalMatrix(q::StateVec, tmp::StateVec, Δt::T, ρ::S, a::S, K::S, grid::Grid{T}, ::DifussionImplicit, i::Int) where {S<:Symbol, T}
  a_ = [   -Δt*grid.Δzi2*     (tmp[ρ, Dual(k)] .* q[a, Dual(k), i] .* tmp[K, Dual(k), i])[2] for k in over_elems_real(grid)[1:end-1]]
  b_ = [1 + Δt*grid.Δzi2*(sum((tmp[ρ, Dual(k)] .* q[a, Dual(k), i] .* tmp[K, Dual(k), i])))  for k in over_elems_real(grid)]
  c_ = [   -Δt*grid.Δzi2*     (tmp[ρ, Dual(k)] .* q[a, Dual(k), i] .* tmp[K, Dual(k), i])[1] for k in over_elems_real(grid)[2:end]]
  return Tridiagonal(a_,b_,c_)
end

"""
    solve_tdma!(q::StateVec, tendencies::StateVec, tmp::StateVec, x::S, ρ::S, a_τp1::S, a_τ::S, K::S, grid::Grid{T}, Δt::T, i::Int=1)

Solve tri-diagonal system using Diffusion-implicit time-marching.
"""
function solve_tdma!(q::StateVec, tendencies::StateVec, tmp::StateVec, x::S, ρ::S, a_τp1::S, a_τ::S, K::S, grid::Grid{T}, Δt::T, i::Int=1) where {S<:Symbol, T}
  B = [tmp[a_τ, k, i]/q[a_τp1, k, i]*q[x, k, i] +
  Δt*tendencies[x, k, i]/(tmp[ρ, k]*q[a_τp1, k, i])
  for k in over_elems_real(grid)]
  A = TridiagonalMatrix(q, tmp, Δt, ρ, a_τp1, K, grid, DifussionImplicit(), i)
  X = inv(A)*B
  assign_real!(q, x, grid, X, i)
end

"""
    solve_tridiag_wrapper!(grid::Grid, sv::StateVec, ϕ::Symbol, i::Int, tri_diag::StateVec)

Solve the tri-diagonal system, given by state vector `tri_diag` for variable `ϕ` in sub-domain `i`.
"""
function solve_tridiag_wrapper!(grid::Grid, sv::StateVec, ϕ::Symbol, i::Int, tri_diag::StateVec)
  f = [tri_diag[:f, k] for k in over_elems_real(grid)]
  a = [tri_diag[:a, k] for k in over_elems_real(grid)]
  b = [tri_diag[:b, k] for k in over_elems_real(grid)]
  c = [tri_diag[:c, k] for k in over_elems_real(grid)]
  solve_tridiag_old(grid.n_elem_real, f, a, b, c)
  assign_real!(sv, ϕ, grid, f, i)
end

# TODO: Replace this with Tridiagonal/inv
function solve_tridiag_old(nz::I, x::Vector{T}, a::Vector{T}, b::Vector{T}, c::Vector{T}) where {T<:AbstractFloat,I<:Int}
  scratch = deepcopy(x)
  scratch[1] = c[1]/b[1]
  x[1] = x[1]/b[1]
  @inbounds for i in 2:nz
    m = 1/(b[i] - a[i] * scratch[i-1])
    scratch[i] = c[i] * m
    x[i] = (x[i] - a[i] * x[i-1])*m
  end
  @inbounds for i in nz-1:-1:1
    x[i] = x[i] - scratch[i] * x[i+1]
  end
end

end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA