text:
"Abstract type for GraphQL objects"
abstract type GraphQLType end

"Object type"
struct Object <: GraphQLType
    name::String
end
"InputObject type"
struct InputObject <: GraphQLType
    name::String
end

const GQL_DEFAULT_SCALAR_TO_JULIA_TYPE = Dict(
    "String" => String,
    "Float" => Float64,
    "Boolean" => Bool,
    "Int" => Int64,
    "ID" => String,
)

"""
    getjuliatype(field; level=:top, scalar_types::Dict{String, DataType}=GQL_DEFAULT_SCALAR_TO_JULIA_TYPE)

Returns type of a GraphQL field by determining the kind of the field and
acting accordingly.

# Arguments
- `field`: introspected field dictionary.
- `scalar_types::Dict{String, DataType}=GQL_DEFAULT_SCALAR_TO_JULIA_TYPE`: optional dictionary
    of scalar type name to Julia type which defaults to internal dictionary containing default
    scalars only.
"""
function getjuliatype(field; scalar_types::Dict{String, DataType}=GQL_DEFAULT_SCALAR_TO_JULIA_TYPE)
    if isobject(field)
        return Object
    elseif isinputobject(field)
        return InputObject
    elseif isscalar(field)
        name = getsubfield(field)["name"]
        if !haskey(scalar_types, name)
            if scalar_types == GQL_DEFAULT_SCALAR_TO_JULIA_TYPE
                throw(ArgumentError("Scalar type $(name) is not a default GraphQL SCALAR type. Use the scalar_types kwarg."))
            else
                throw(ArgumentError("Scalar type $(name) is not defined in scalar_types."))
            end
        end
        return scalar_types[name]
    elseif isenum(field)
        return String
    elseif islist(field)
        return Vector{getjuliatype(getsubfield(field); scalar_types=scalar_types)}
    elseif isnonnull(field)
        return getjuliatype(getsubfield(field); scalar_types=scalar_types)
    else
        throw(GraphQLClientException("getjuliatype function does not handle types of \"$(gettype(field))\""))
    end
end

"""
    get_field_type_string(field)

Returns string of the type of a GQL field by determining the kind of the field and
acting accordingly.
"""
function get_field_type_string(field)
    if isscalar(field) || isinputobject(field) || isenum(field) || isobject(field)
        return getsubfield(field)["name"]
    elseif islist(field)
        # Add square brackets and get type from next level down
        return "[" * get_field_type_string(getsubfield(field)) * "]"
    elseif isnonnull(field)
        # Add ! to show that it is required and get type from next level down
        return get_field_type_string(getsubfield(field)) * "!"
    else
        throw(GraphQLClientException("get_field_type_string function does not handle types of \"$(gettype(field))\""))
    end
end

"""
    getsubfield(field)

Return `field["type"]` or `field["ofType"]`, throwing an error if neither exist.
"""
function getsubfield(field)
    haskey(field, "type") && return field["type"]
    haskey(field, "ofType") && return field["ofType"]
    throw(GraphQLClientException("$field is not a field, cannot return subfield."))
end

gettype(field) = getsubfield(field)["kind"]

"""
    istype(field, comparison)

Checks whether kind of `field` (which has either `type` or `ofType` as a key)
is `comparison`. Error if both `type` and `ofType` are not keys.
"""
istype(field, comparison) = gettype(field) == comparison

# Generate the simple cases
for type in ("OBJECT", "LIST", "SCALAR", "ENUM")
    f = Symbol("is", lowercase(type))
    @eval begin
        $f(field) = istype(field, $type)
    end
end

# Handle these separately due to underscores and more complicated logic
isnonnull(field) = istype(field, "NON_NULL")
isinputobject(field) = istype(field, "INPUT_OBJECT")
is_nonnull_input_object(field) = isnonnull(field) && istype(getsubfield(field), "INPUT_OBJECT")

"""
    getroottypefield(field)

Get the root type field of a NON_NULL or LIST field.
"""
function getroottypefield(field)
    if islist(field) || isnonnull(field)
        getroottypefield(getsubfield(field))
    else
        getsubfield(field)
    end
end

"""
    getroottype(field)

Get the name of the root type of a field.
"""
getroottype(field) = getroottypefield(field)["name"]

"""
    isroottypeenum(field)

Return `true` if field is an arbitrarily nested enum (nested in both NON_NULL and/or LISTs).
"""
isroottypeenum(field) = isenum(field) || ((islist(field) || isnonnull(field)) && isroottypeenum(getsubfield(field)))

"""
    isroottypeobject(field)

Return `true` if field is an arbitrarily nested enum (nested in both NON_NULL and/or LISTs).
"""
isroottypeobject(field) = isobject(field) || ((islist(field) || isnonnull(field)) && isroottypeobject(getsubfield(field)))

"""
    isroottypescalar(field)

Return `true` if field is an arbitrarily nested scalr (nested in both NON_NULL and/or LISTs).
"""
isroottypescalar(field) = isscalar(field) || ((islist(field) || isnonnull(field)) && isroottypescalar(getsubfield(field)))


text:
@testset "Karp Minimum Cycle Mean" begin
     w = [0. 0.  0. 0. 0. 0. 1. 0.  0.  0.0  0.  0.
          0. 0.  0. 0. 0. 0. 1. 0.  0.  0.0  0.  0.
          0. 0.  0. 0. 0. 0. 0. 0.  1.1 0.0  0.  0.
          0. 1.2 0. 0. 0. 0. 1. 0.  0.  0.0  0.  0.
          0. 1.  0. 0. 0. 0. 0. 0.  1.3 0.0  0.  0.
          0. 1.2 0. 0. 0. 0. 1. 0.  0.  0.0  0.  0.
          0. 1.  0. 0. 0. 0. 0. 0.  1.3 0.0  0.  0.
          0. 0.  0. 0. 0. 0. 0. 0.  0.  0.0  0.  0.
          0. 0.  0. 0. 0. 0. 0. 0.  0.  0.0  0.4 0.
          0. 0.  0. 0. 0. 0. 1. 0.  0.  0.0  0.  0.
          0. 0.  0. 0. 0. 0. 1. 0.  0.  0.0  0.  0.
          0. 0.  0. 0. 0. 0. 0. 0.  1.2 0.0  0.  0.]

    w2 = [Inf  Inf    Inf  Inf  Inf  Inf   -0.09 Inf  Inf    Inf  Inf    Inf
          Inf  Inf    Inf  Inf  Inf  Inf   -0.06 Inf  Inf    Inf  Inf    Inf
          Inf  Inf    Inf  Inf  Inf  Inf  Inf    Inf   -0.01 Inf  Inf    Inf
          Inf   -0.19 Inf  Inf  Inf  Inf   -0.09 Inf  Inf    Inf  Inf    Inf
          Inf   -0.02 Inf  Inf  Inf  Inf  Inf    Inf   -0.03 Inf  Inf    Inf
          Inf   -0.19 Inf  Inf  Inf  Inf   -0.09 Inf  Inf    Inf  Inf    Inf
          Inf   -0.02 Inf  Inf  Inf  Inf  Inf    Inf   -0.03 Inf  Inf    Inf
          Inf  Inf    Inf  Inf  Inf  Inf  Inf    Inf  Inf    Inf  Inf    Inf
          Inf  Inf    Inf  Inf  Inf  Inf  Inf    Inf  Inf    Inf    0.39 Inf
          Inf  Inf    Inf  Inf  Inf  Inf   -0.09 Inf  Inf    Inf  Inf    Inf
          Inf  Inf    Inf  Inf  Inf  Inf   -0.06 Inf  Inf    Inf  Inf    Inf
          Inf  Inf    Inf  Inf  Inf  Inf  Inf    Inf   -0.01 Inf  Inf    Inf]

     g1 = SimpleDiGraph(12)
     add_edge!(g1, 1, 7)
     add_edge!(g1, 2, 7)
     add_edge!(g1, 3, 9)
     add_edge!(g1, 4, 2)
     add_edge!(g1, 4, 7)
     add_edge!(g1, 5, 2)
     add_edge!(g1, 5, 9)
     add_edge!(g1, 6, 2)
     add_edge!(g1, 6, 7)
     add_edge!(g1, 7, 2)
     add_edge!(g1, 7, 9)
     add_edge!(g1, 9, 11)
     add_edge!(g1, 10, 7)
     add_edge!(g1, 11, 7)
     add_edge!(g1, 12, 9)

     @testset "simple digraphs" for g in testgraphs(g1)
        c, λ = karp_minimum_cycle_mean(g, w)
        @test c == [9, 11, 7]
        @test λ == 0.9

        c2, λ2 = karp_minimum_cycle_mean(g, w2)
        @test λ2 == -0.04
        @test c2 == [2, 7]
     end

    # Tricky test case
    # Backward walk from 3 is good:
    # 1 - 2 - 3 - 4 - 3
    # Only the first of the two backward walk from 4 is good:
    # 1 - 3 - 4 - 3 - 4
    # 1 - 2 - 2 - 3 - 4
    # We will do the backward walk from 3 in the test since 0 = F_4(3) < F_4(4) = 1
    tricky = SimpleDiGraph(4)
    add_edge!(tricky, 1, 3)
    add_edge!(tricky, 1, 2)
    add_edge!(tricky, 3, 4)
    add_edge!(tricky, 4, 1) # Make it strongly connected
    add_edge!(tricky, 4, 3)
    add_edge!(tricky, 2, 3)
    add_edge!(tricky, 2, 2)
    distmx = [Inf  0.  1. Inf
              Inf  1.  0. Inf
              Inf Inf Inf  0.
               1. Inf  0. Inf]

    @testset "tricky case" for g in testgraphs(tricky)
        c, λ = karp_minimum_cycle_mean(g, distmx)
        @test λ == 0.
        @test sort(c) == [3, 4]

        c, λ = karp_minimum_cycle_mean(g, distmx .- 2)
        @test λ == -2.
        @test sort(c) == [3, 4]
    end

    # Multiple strongly connected components and one component with one node and no edge
    multi = SimpleDiGraph(3)
    add_edge!(multi, 1, 1)
    add_edge!(multi, 2, 2)
    distmx = [0. Inf
              Inf -1]

    @testset "multiple SCCs" for g in testgraphs(multi)
        c, λ = karp_minimum_cycle_mean(g, distmx)
        @test λ == -1.
        @test c == [2]
    end
end

ext: jl


text:
using DataFrames
using CSV
using HorseML.Clustering
using HorseML.Clustering: fit!

@testset "Clustering" begin
    sy, st = rand(), rand()
    y, t = rand(20), rand(20)
    my1, my2, my3 = fill!(Array{Float64}(undef, 1, 1), sy), y[:, :], rand(20, 2)

    #Load the data
    df = CSV.read("clustering.csv", DataFrame)
    K = 3
    x = Array(df)

    #Test for Kmeans
    model = Kmeans(K)
    @test_nowarn fit!(model, x)
    @test_nowarn model(x)

    μ = model.μ
    @test_throws DomainError dm(sy, st, μ)
    @test typeof(dm(y, t, μ)) == Float64
    @test typeof(dm(y, t, μ, reduction = "sum")) == Float64
    @test size(dm(y, t, μ, reduction = "none")) == (20, )
    @test_throws ArgumentError dm(y, t, μ, reduction="NaN")
    @test_throws DomainError dm(my1, st, μ)
    @test dm(my2, st, μ) == dm(vec(my2), st, μ)
    @test_throws DimensionMismatch dm(my3, st, μ)
    @test dm(my2, t, μ) == dm(y, t, μ)
    @test_throws DimensionMismatch dm(my3, t, μ)

    #Test for GMM(Gaussian Mixture Model)
    model = GMM(K)
    @test_nowarn fit!(model, x)
    @test_nowarn model(x)

    π, μ, Σ = model.π, model.μ, model.Σ
    @test nlh(x, π, μ, Σ) ≈ 215.78213657665862
end


text:
# const ABI_UNIT_NAMES = [lowercase(s) for s in [
#         "au",
#         "Angstr", "Angstrom", "Angstroms", "Bohr", "Bohrs",
#         "eV", "Ha", "Hartree", "Hartrees", "K", "Ry", "Rydberg", "Rydbergs",
#         "T", "Tesla"]]
#
# const abi_conversions = Dict{Symbol,Any}(:ev => 1 / 27.2113845,
#                                          :ha => 1.0,
#                                          :ry => 0.5,
#                                          :ang => 1.889716164632,
#                                          :bohr => 1.0,
#                                          :au => 1.0,
#                                          :angstr => 1.889716164632,
#                                          :angstrom => 1.889716164632,
#                                          :bohrs => 1.0,
#                                          :hartree => 1.0,
#                                          :hartrees => 1.0,
#                                          :k => 1/3.1577513e5,
#                                          :rydberg => 0.5,
#                                          :rydbergs => 0.5,
#                                          :t =>1/2.35e5,
#                                          :tesla => 1/2.35e5)
#
# convert_2abi(value, s::String) = value * abi_conversions[Symbol(s)]
# convert_2abi(value, s::Symbol) = value * abi_conversions[s]
#
# function construct_abi_flags()
#     out      = Dict{Symbol, Type}()
#     open(joinpath(assets_dir,"inputs/abinit/input_variables.txt"), "r") do f
#         while !eof(f)
#             spl = split(readline(f))
#             out[Meta.parse(spl[1])] = eval(Meta.parse(spl[2]))
#         end
#     end
#     return out
# end
#

include(joinpath(depsdir, "abinitflags.jl"))
const AbinitFlags  = _ABINITFLAGS()
#
# flagtype(input::DFInput{Abinit}, flag) = haskey(AbinitFlags, flag) ? AbinitFlags[flag] : Nothing


text:
using SafeTestsets

@safetestset "Wald" begin
    using Test, SequentialSamplingModels, KernelDensity, Random
    Random.seed!(22158)
    d = Wald(2, 1, .1)
    @test mean(d) ≈ (1/2) + .1 atol = 1e-5

    function simulate(υ, α, θ)
        noise = 1.0
        #Time Step
        Δt = .0005
        #Evidence step
        Δe = noise*sqrt(Δt)
        e = 0.0
        t = θ
        p = .5*(1 + υ*sqrt(Δt)/noise)
        while (e < α)
            t += Δt
            e += rand() ≤ p ? Δe : -Δe
        end
        return t
    end
    rts = map(_->simulate(3, 1, .2), 1:10^5)
    approx_pdf = kde(rts)
    x = .2:.01:1.5
    y′ = pdf(approx_pdf, x)
    y = pdf.(Wald(3,1,.2), x)
    @test y′ ≈ y rtol = .03
    @test mean(rts) ≈ mean(Wald(3,1,.2)) atol = 5e-3
    @test std(rts) ≈ std(Wald(3,1,.2)) atol = 1e-3
end

@safetestset "Wald Mixture" begin
    using Test, SequentialSamplingModels, KernelDensity, Random
    Random.seed!(22198)
    d = WaldMixture(2, .2, 1, .1)
    @test mean(d) ≈ (1/2) + .1 atol = 1e-5
    rts = rand(d, 100000)
    approx_pdf = kde(rts)
    x = .2:.01:1.5
    y′ = pdf(approx_pdf, x)
    y = pdf.(d, x)
    @test y′ ≈ y rtol = .03
    @test mean(rts) ≈ mean(d) atol = 5e-3

    y′ = @. logpdf(d, x) |> exp
    @test y′ ≈ y
end

@safetestset "LogNormal Race Tests" begin
    using SequentialSamplingModels, Test, Random, KernelDensity, Distributions
    Random.seed!(54054)
    d1 = LNR(;μ=[1.0], σ=1.0, ϕ=.1)
    v1 = .3
    p1 = pdf(d1, 1, v1)
    p2 = pdf(LogNormal(1, 1), v1-.1)
    @test p1 ≈ p2
    d2 = LNR(;μ=[1.0,0.0],σ=1.0,ϕ=.1)
    d3 = LNR(;μ=[1.0,0.0,0.0], σ=1.0, ϕ=.1)
    p1 = pdf(d2, 1, v1)
    p2 = pdf(d3, 1 ,v1)
    @test p1 > p2
    @test p1 ≈ pdf(LogNormal(1, 1), v1-.1)*(1-cdf(LogNormal(0, 1), v1-.1))

    m1,m2=-2,-1
    σ = .9
    ϕ = 0.0
    d = LNR(μ=[m1,m2], σ=σ, ϕ=ϕ)
    data = rand(d, 10^4)
    x = range(m1*.8, m1*1.2, length=100)
    y = map(x->sum(logpdf.(LNR(μ=[x,m2], σ=σ, ϕ=ϕ), data)), x)
    mv,mi = findmax(y)
    @test m1 ≈ x[mi] atol = .05

    x = range(m2*.8, m2*1.2, length=100)
    y = map(x->sum(logpdf.(LNR(μ=[m1,x], σ=σ, ϕ=ϕ), data)), x)
    mv,mi = findmax(y)
    @test m2 ≈ x[mi] atol = .05

    x = range(σ*.8, σ*1.2, length=100)
    y = map(x->sum(logpdf.(LNR(μ=[m1,m2], σ=x, ϕ=ϕ), data)), x)
    mv,mi = findmax(y)
    @test σ ≈ x[mi] atol = .05

    x = range(ϕ*.8, ϕ*1.2, length=100)
    y = map(x->sum(logpdf.(LNR(μ=[m1,m2], σ=σ, ϕ=x), data)), x)
    mv,mi = findmax(y)
    @test ϕ ≈ x[mi] atol = .0005

    d = LNR(μ=[m1,m2], σ=σ, ϕ=ϕ)
    data = rand(d, 10^4)
    y1 = map(x->logpdf(d, x), data)
    y2 = map(x->log(pdf(d, x)), data)
    @test y1 ≈ y2 atol = .00001

    d = LNR(μ=[1.,.5], σ=.4, ϕ=.5)
    data = rand(d, 10^4)
    y1 = map(x->logpdf(d, x), data)
    y2 = map(x->log(pdf(d, x)), data)
    @test y1 ≈ y2 atol = .00001

    dist = LNR(μ=[-1.5,-.9], σ=.5, ϕ=.3)
    data = rand(dist, 10^5)
    data1 = filter(x->x[1] == 1, data)
    p1 = mean(x->x[1] == 1, data)
    p2 = 1 - p1
    rt1 = map(x->x[2], data1)
    approx_pdf = kde(rt1)
    x = .2:.01:1.5
    y′ = pdf(approx_pdf, x)*p1
    y = pdf.(dist, (1,), x)
    @test y′ ≈ y rtol = .03

    data2 = filter(x->x[1] == 2, data)
    rt2 = map(x->x[2], data2)
    approx_pdf = kde(rt2)
    x = .2:.01:1.5
    y′ = pdf(approx_pdf, x)*p2
    y = pdf.(dist, (2,), x)
    @test y′ ≈ y rtol = .03
end

@safetestset "LBA Tests" begin
    using SequentialSamplingModels, Test, KernelDensity, Random
    Random.seed!(10542)
    dist = LBA(ν=[3.0,2.0], A = .8, k = .2, τ = .3)
    choice,rt = rand(dist, 10^5)
    rt1 = rt[choice .== 1]
    p1 = mean(x->x==1, choice)
    p2 = 1 - p1
    approx_pdf = kde(rt1)
    x = .2:.01:1.5
    y′ = pdf(approx_pdf, x)*p1
    y = pdf.(dist, (1,), x)
    @test y′ ≈ y rtol = .03

    rt2 = rt[choice .== 2]
    approx_pdf = kde(rt2)
    x = .2:.01:1.5
    y′ = pdf(approx_pdf, x)*p2
    y = pdf.(dist, (2,), x)
    @test y′ ≈ y rtol = .03
end

@safetestset "Racing Diffusion Model" begin
    # cd(@__DIR__)
    # using Pkg
    # Pkg.activate("..")

    using SequentialSamplingModels, Test, KernelDensity, QuadGK, Random
    using Interpolations, Distributions
    #import KernelDensity: kernel_dist
    import SequentialSamplingModels: WaldA
    kernel_dist(::Type{Epanechnikov}, w::Float64) = Epanechnikov(0.0, w)
    kernel(data) = kde(data; kernel=Epanechnikov)
    Random.seed!(741)

    dist = WaldA(ν=.5, k=.3, A=.7, θ=.2)
    rts = map(_->rand(dist), 1:10^6)
    approx_pdf = kernel(rts)
    x = .201:.01:2.5
    y′ = pdf(approx_pdf, x)
    y = pdf.(dist, x)
    @test mean(abs.(y .- y′)) < .02
    #@test std(abs.(y .- y′)) < .04

    p′ = quadgk(x->pdf(dist, x), .2, Inf)[1]
    @test p′ ≈ 1 rtol = .001

    p = cdf(dist, .25)
    p′ = quadgk(x->pdf(dist, x), .2, .25)[1]
    @test p′ ≈ p rtol = .001
    @test p ≈ mean(rts .< .25) rtol = .01

    p = cdf(dist, .5)
    p′ = quadgk(x->pdf(dist, x), .2, .5)[1]
    @test p′ ≈ p rtol = .001
    @test p ≈ mean(rts .< .5) rtol = .01

    p = cdf(dist, .6) - cdf(dist, .3)
    p′ = quadgk(x->pdf(dist, x), .3, .6)[1]
    @test p′ ≈ p rtol = .001
    @test p ≈ mean((rts .< .6 ).& (rts .> .3)) rtol = .01

    dist = WaldA(ν=1.0, k=.3, A=1.0, θ=.2)
    rts = map(_->rand(dist), 1:10^6)
    approx_pdf = kernel(rts)
    x = .201:.01:2.5
    y′ = pdf(approx_pdf, x)
    y = pdf.(dist, x)
    @test mean(abs.(y .- y′)) < .02
   # @test std(abs.(y .- y′)) < .04

    p′ = quadgk(x->pdf(dist, x), .2, Inf)[1]
    @test p′ ≈ 1 rtol = .001

    p = cdf(dist, .25)
    p′ = quadgk(x->pdf(dist, x), .2, .25)[1]
    @test p′ ≈ p rtol = .001
    @test p ≈ mean(rts .< .25) rtol = .01

    p = cdf(dist, .5)
    p′ = quadgk(x->pdf(dist, x), .2, .5)[1]
    @test p′ ≈ p rtol = .001
    @test p ≈ mean(rts .< .5) rtol = .01

    p = cdf(dist, 1.0) - cdf(dist, .6)
    p′ = quadgk(x->pdf(dist, x), .6, 1.0)[1]
    @test p′ ≈ p rtol = .001
    @test p ≈ mean((rts .< 1).& (rts .> .6)) rtol = .01

    p = cdf(dist, 1.5) - cdf(dist, 1.4)
    p′ = quadgk(x->pdf(dist, x), 1.4, 1.5)[1]
    @test p′ ≈ p rtol = .001
    @test p ≈ mean((rts .< 1.5).& (rts .> 1.4)) rtol = .02

    dist = DiffusionRace(;ν=[1.0,.5], k=0.5, A=1.0, θ=.2)
    data = rand(dist, 10^6)
    data1 = filter(x->x[1] == 1, data)
    p1 = mean(x->x[1] == 1, data)
    p2 = 1 - p1
    rt1 = map(x->x[2], data1)
    # approx_pdf = kernel(rt1)
    # x = .201:.01:2.5
    # y′ = pdf(approx_pdf, x) * p1
    # y = pdf.(dist, (1,), x)
    # @test y′ ≈ y rtol = .03

    data2 = filter(x->x[1] == 2, data)
    rt2 = map(x->x[2], data2)
    approx_pdf = kernel(rt2)
    # x = .201:.01:1.5
    # y′ = pdf(approx_pdf, x) * p2
    # y = pdf.(dist, (2,), x)
    # @test y′ ≈ y rtol = .03

    p′ = quadgk(x->pdf(dist, 1, x), .2, Inf)[1]
    @test p′ ≈ p1 rtol = .001

    p = quadgk(x->pdf(dist, 1, x), .2, .3)[1]
    @test p ≈ mean(rt1 .< .3) * p1 atol = .01

    p = quadgk(x->pdf(dist, 1, x), .2, .5)[1]
    @test p ≈ mean(rt1 .< .5) * p1 atol = .01

    p = quadgk(x->pdf(dist, 1, x), .2, 1)[1] - quadgk(x->pdf(dist, 1, x), .2, .6)[1]
    @test p ≈ mean((rt1 .< 1) .& (rt1 .> .6)) * p1 atol = .01

    p = quadgk(x->pdf(dist, 1, x), .2, 1.5)[1] - quadgk(x->pdf(dist, 1, x), .2, 1.4)[1]
    @test p ≈ mean((rt1 .< 1.5).& (rt1 .> 1.4)) atol = .01

    p′ = quadgk(x->pdf(dist, 2, x), .2, Inf)[1]
    @test p′ ≈ p2 rtol = .001

    p = quadgk(x->pdf(dist, 2, x), .2, .3)[1]
    @test p ≈ mean(rt2 .< .3) * p2 atol = .01

    p = quadgk(x->pdf(dist, 2, x), .2, .5)[1]
    @test p ≈ mean(rt2 .< .5) * p2 atol = .02

    p = quadgk(x->pdf(dist, 2, x), .2, 1)[1] - quadgk(x->pdf(dist, 2, x), .2, .6)[1]
    @test p ≈ mean((rt2 .< 1) .& (rt2 .> .6)) * p2 atol = .01

    p = quadgk(x->pdf(dist, 2, x), .2, 1.5)[1] - quadgk(x->pdf(dist, 2, x), .2, 1.4)[1]
    @test p ≈ mean((rt2 .< 1.5).& (rt2 .> 1.4)) * p2 atol = .01
end

@safetestset "Attentional Diffusion" begin
    using Test, SequentialSamplingModels, StatsBase, Random
    Random.seed!(5511)

    mutable struct Transition
        state::Int
        n::Int
        mat::Array{Float64,2}
    end

    function Transition(mat)
        n = size(mat,1)
        state = rand(1:n)
        return Transition(state, n, mat)
    end

    function attend(t)
        w = t.mat[t.state,:]
        next_state = sample(1:t.n, Weights(w))
        t.state = next_state
        return next_state
    end

    model = AttentionalDiffusion(;ν1=5.0, ν2=4.0)

    tmat = Transition(
        [.98 .015 .005;
        .015 .98 .005;
        .45 .45 .1]
        )

    rts1 = rand(model, 1000, x->attend(x), tmat)
    n1,n2 = length.(rts1)
    @test n1 > n2

    model = AttentionalDiffusion(;ν1=4.0, ν2=5.0)
    rts2 = rand(model, 1000, x->attend(x), tmat)
    n1,n2 = length.(rts2)
    @test n1 < n2

    μ_rts1 = mean(vcat(rts1...))
    model = AttentionalDiffusion(;ν1=5.0, ν2=5.0)
    rts3 = rand(model, 1000, x->attend(x), tmat)
    μ_rts3 = mean(vcat(rts3...))
    @test μ_rts1 < μ_rts3
end


text:
"""
    abstract type AbstractCXXRecordDecl <: AbstractRecordDecl
Supertype for `CXXRecordDecl`s.
"""
abstract type AbstractCXXRecordDecl <: AbstractRecordDecl end

"""
    struct CXXRecordDecl <: AbstractCXXRecordDecl
Hold a pointer to a `clang::CXXRecordDecl` object.
"""
struct CXXRecordDecl <: AbstractCXXRecordDecl
    ptr::CXCXXRecordDecl
end

Base.unsafe_convert(::Type{CXCXXRecordDecl}, x::CXXRecordDecl) = x.ptr
Base.cconvert(::Type{CXCXXRecordDecl}, x::CXXRecordDecl) = x


text:
export currentFramePatchPeriod, getRepetitionsOfSameState, manualCorrection, FT1A05,
       determineWindow, getMeasurementsMotionCompFD, getPeriod

struct LocationInMeasurement
  frame::Int
  patch::Int
  period::Float64
end

"""
    currentFramePatchPeriod(currentTimeInms, numPeriodsPerFrame, numPeriodsPerPatch,
                            numPatches, tsc[; DFCyclesSwitch==7])
Determine frame, patch, and period from time within the measurement sequence.
"""
function currentFramePatchPeriod(currentTimeInms, numPeriodsPerFrame, numPeriodsPerPatch,
                                 numPatches, tsc; DFCyclesSwitch=7)

  frame = ceil(Int, currentTimeInms/((numPeriodsPerFrame+DFCyclesSwitch*numPatches)*tsc))
  rest = currentTimeInms-(frame-1)*(numPeriodsPerFrame+DFCyclesSwitch*numPatches)*tsc
  patch = ceil(Int, rest/((numPeriodsPerPatch+DFCyclesSwitch)*tsc))
  rest2 = rest-(patch-1)*(numPeriodsPerPatch+DFCyclesSwitch)*tsc
  period = rest2/tsc
  if currentTimeInms == 0
    frame = 1
    patch = 1
    period = 1.0
  end
  return LocationInMeasurement(frame, patch, period)
end

"""
    calculateRepetitionsInFramePatchPeriod(totalduration, motFreq, currentTimeInms,
                      endtime, numPeriodsPerPatch, numPeriodsPerFrame, numPatches,tsc)

Based on the starting point (currentTimeInms) in time, the repetitions are
calculated from the motion frequency at that moment (entries in motFreq) until
the endtime is reached. The times are converted to frame, patch, and period.

"""
function calculateRepetitionsInFramePatchPeriod(totalduration, motFreq, currentTimeInms, endtime,
                                        numPeriodsPerPatch, numPeriodsPerFrame, numPatches, tsc)

  #N = ceil(Int, totalduration*maximum(motFreq)) + 1
  tmot = Vector{LocationInMeasurement}()

  while (currentTimeInms < endtime)
    loc = currentFramePatchPeriod(currentTimeInms, numPeriodsPerFrame,
                                            numPeriodsPerPatch, numPatches,tsc)

    currentTimeInms += 1 / motFreq[loc.frame,loc.patch]

     if (loc.period+1/(motFreq[loc.frame,loc.patch])/tsc < numPeriodsPerPatch-1 &&
         loc.period >= 2)
       push!(tmot, loc)
     end
  end
  return tmot
end

function calculateBeginningInms(firstFrame, tsc, numPeriodsPerPatch, numPatches, DFCyclesSwitch)
  return (firstFrame-1)*tsc*(numPeriodsPerPatch*numPatches+DFCyclesSwitch*numPatches)
end

"""
    getRepetitionsOfSameState(bMeas::MPIFile,motFreq,frames[;DFCyclesSwitch==7])

Based on the motion frequency the repetitions of the same state are calculated. I.e. the frame, patch, and period of the repetition

#Arguments:
- `motFreq::Float`: Assumed motion frequency in Hz
"""
function getRepetitionsOfSameState(bMeas::MPIFile, motFreq, frames::UnitRange;
                                   DFCyclesSwitch=7)

  numPeriodsPerFrame = acqNumPeriodsPerFrame(bMeas)
  numPatches = acqNumPatches(bMeas)
  numPeriodsPerPatch = acqNumPeriodsPerPatch(bMeas)
  endtime = last(frames)*dfCycle(bMeas)*(numPeriodsPerPatch*numPatches+DFCyclesSwitch*numPatches)
  currentTimeInms = calculateBeginningInms(first(frames), dfCycle(bMeas), numPeriodsPerPatch,
                                           numPatches, DFCyclesSwitch)
  totalduration = endtime-currentTimeInms

  tmot = calculateRepetitionsInFramePatchPeriod(totalduration, motFreq, currentTimeInms,
                          endtime, numPeriodsPerPatch, numPeriodsPerFrame, numPatches, dfCycle(bMeas))

  return tmot
end

function FT1A05(windowLength)
  window = zeros(windowLength)
  n = 0:windowLength-1
  window = 0.2769-0.5261.*cos.(2*pi.*n./windowLength)+0.1971.*cos.(4*pi.*n./windowLength)
  return window
end

"""
    determineWindow(totalLength, filledLength, window)

Calculate window function with different window types of length filledLength.
If totalLength > filledLength, window function is filled with zeros.
"""
function determineWindow(totalLength, filledLength, windowType)
  window = zeros(totalLength)
  empty = 0#floor(Int,(totalLength-filledLength)/2)
  if windowType == 1
    window[empty+1:empty+filledLength] = hanning(filledLength)
  elseif windowType == 2
    window[empty+1:empty+filledLength] = FT1A05(filledLength)
  elseif windowType == 3
    window[empty+1:empty+filledLength] = ones(filledLength)
  end
  return window
end

"""
    numDFPeriodsInMotionCycle(motFreq, frames, tsc)

Based on the motion frequency in every frame and patch, the maximum number of
complete DF cycles is caluclated for the completion of one motion cycle

"""
function numDFPeriodsInMotionCycle(motFreq::Matrix, frames, tsc)
  return floor(Int, 1.0 / (minimum(motFreq[frames,:])*tsc) )
end

function numDFPeriodsInMotionCycle(motFreq::Number, tsc)
  return floor(Int, 1.0 / (motFreq*tsc) )
end

"""
   calculateSubDFPeriodShift(currentPeriod, incr, samplesInOneDFCycle)

currentPeriod is given as floating number (not an Integer although this seems to be right).
The decimal digits give the shift < DF duration.
"""
function calculateSubDFPeriodShift(currentPeriod, incr, samplesInOneDFCycle, samplingPrecision)
  delta = floor(Int,(currentPeriod+incr-floor(Int,currentPeriod+incr))*samplesInOneDFCycle)
  if samplingPrecision == false
    delta = 0
  end
  return delta
end

function fillDataIntoVirtualFrame!(ui, ufinal, normalization, numMotPeriods,
                   incrementPerPeriod, samplesInOneDFCycle, currentPeriod,
                   currentPatch, window, samplingPrecision)

  numDFCyclesPerWindowLength = div(length(window), size(ui,1))
  N = size(ui,1)

  for period = 1:numMotPeriods
    incr = incrementPerPeriod*(period-1)
    delta = calculateSubDFPeriodShift(currentPeriod,incr,samplesInOneDFCycle,samplingPrecision)

    centralPeriod = floor(Int, currentPeriod+incr)
    shiftedWindow = circshift(window, delta)

    for recChan = 1:size(ui,2)
      p = (centralPeriod-1):(centralPeriod-1+numDFCyclesPerWindowLength)
      buf = vec(ui[:,recChan,p,currentPatch])
      buf2 = circshift(buf[delta+1:delta+N*numDFCyclesPerWindowLength,:], delta)
      buf3 = buf2 .* shiftedWindow

      for l=1:numDFCyclesPerWindowLength
        ufinal[:,recChan,period,currentPatch] += buf3[((l-1)*N+1):(l*N)]
        normalization[:,recChan,period,currentPatch] +=
                             shiftedWindow[((l-1)*N+1):(l*N)]
      end
    end
  end
  return nothing
end

"""
    normalizeSignalLevelandFormatData(ufinal,numPatches,freq,count)

Normalize data based on number of repetitions to ensure same signal level for all patches and reformatting
"""
function normalizeSignalLevelandFormatData(ufinal, normalization, freq)
  if any(normalization .== 0)
    error("Could not construct virtual frame!")
  end
  ufinal ./= normalization

  ufinal = rfft(ufinal,1)
  ufinal = permutedims(reshape(ufinal,size(ufinal)[1]*size(ufinal)[2],
                               size(ufinal)[3],size(ufinal)[4])[freq,:,:],[1,3,2])
  return ufinal
end

"""
	getMeasurementsMotionCompFD(bMeas::MPIFile, motFreq, tmot, freq, frames, sigma,
                      samplingPrecision, windowType)
	Averages over raw data corresponding to the same motion state and creates the virtual frames

- bMeas:		Meas data
- motFreq:		Array containing dominant motion frequencies for each patch and frame
- tmot:			Array containing repetitions of the same state (frame,patch,period)
- freq:			Selected frequencies for reconstruction
- frames:       Selected frames
- alpha:        Window width relative to the DF cycle.
- samplingPrecision:	true: rounding motion period to sampling precision, false: rounding to DF period precision
- window:		1: Hann window, 2:FT1A05, 3: Rectangle

"""
function getMeasurementsMotionCompFD(bMeas::MPIFile, motFreq, tmot, freq, frames::UnitRange,
                                     alpha, samplingPrecision, windowType)
  oldFrame = first(frames)

  numMotPeriods = numDFPeriodsInMotionCycle(motFreq, frames, dfCycle(bMeas))

  ui_ = getMeasurements(bMeas, frames=first(frames), spectralLeakageCorrection=false)
  ui = reshape(ui_, size(ui_,1), 3, div(size(ui_,3), acqNumPatches(bMeas)), acqNumPatches(bMeas))
  ufinal = zeros(Float32, size(ui,1), size(ui,2), numMotPeriods, size(ui,4))
  normalization = zeros(Float32, size(ui,1), size(ui,2), numMotPeriods, size(ui,4))
  count = zeros(acqNumPatches(bMeas))
  numDFCyclesPerWindowLength = ceil(Int, alpha)
  window = determineWindow(size(ui,1)*numDFCyclesPerWindowLength,
                           floor(Int, alpha*size(ui,1)), windowType)

  # collect data from all repetitions from first till last frame and fill into virtual frame ufinal
  for i = 1:length(tmot)
    currentFrame = tmot[i].frame
    if currentFrame > first(frames)-1 && currentFrame < last(frames)+1
      # loading new data if frame has changed
      if currentFrame != oldFrame
        ui[:] = vec(getMeasurements(bMeas, frames=currentFrame, spectralLeakageCorrection=false))
        oldFrame = currentFrame
      end

      currentPatch = tmot[i].patch

      currentMotFreq = motFreq[currentFrame,currentPatch]
      # The number of full DF cycles within one motion cycle varies for different patches and frames
      # For consistent images, the same number is required for all patches.
      # Thus for higher motion in a patch, the increment per Period is lowered
      incrementPerPeriod = numDFPeriodsInMotionCycle(currentMotFreq, dfCycle(bMeas)) / numMotPeriods

      currentPeriod = tmot[i].period
      fillDataIntoVirtualFrame!(ui, ufinal, normalization, numMotPeriods, incrementPerPeriod,
                               size(ui,1), currentPeriod, currentPatch, window,
                               samplingPrecision)
    end
  end

  return normalizeSignalLevelandFormatData(ufinal, normalization, freq)
end


text:
#=
Setup
-----

Packages only need to be installed once per project. Installing packages can
take a few moments.
=#

Pkg.add("RDatasets")

#=
Basic data manipulation
-----------------------

The RDatasets package provides some sample datasets and makes it easy to
import them into Julia as [DataFrames](http://dataframesjl.readthedocs.org/en/latest/).
=#
using RDatasets
mac = dataset("Zelig", "macro")
typeof(mac)

#=
We can extract and compute statistical summaries from individual columns, or
get a statistical summary of the entire dataset.
=#
mean(mac[:Trade])
median(mac[:Trade])
describe(mac)

#=
Subsetting allows us to keep only rows that satisfy a criterion, and datasets
provide many other sophisticated ways to split, modify and combine datasets.
=#
mac[mac[:Year] .> 1985, :]

#=
Plotting
--------

We can use the plotting package Gadfly to visualize the dataset.
=#
using Gadfly
plot(mac, x=:Year, y=:Unem, color=:Country, Geom.line)

## Create a histogram with bin count
plot(dataset("ggplot2", "diamonds"), x=:Price, color=:Cut,
     Geom.histogram(bincount=40))


text:
### Operations and Strategies
include("options/payoff_base.jl")
include("options/engines.jl")
include("options/payoffs.jl")


text:
# !! auto-generated code, do not edit !!
module actionlib_msgs
using RobotOSData.Messages
using RobotOSData.StdMsgs
using RobotOSData.CommonMsgs
struct GoalID <: Readable
    stamp::ROSTime
    id::String
end
struct GoalStatus <: Readable
    goal_id::GoalID
    status::UInt8
    text::String
end
struct GoalStatusArray <: Readable
    header::Header
    status_list::Vector{GoalStatus}
end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module AdjustersTest

using Test
using Dates

@testset "trunc" begin
    dt = DatesPlus.Date(2012, 12, 21)
    @test trunc(dt, DatesPlus.Year) == DatesPlus.Date(2012)
    @test trunc(dt, DatesPlus.Quarter) == DatesPlus.Date(2012, 10)
    @test trunc(dt, DatesPlus.Month) == DatesPlus.Date(2012, 12)
    @test trunc(dt, DatesPlus.Day) == DatesPlus.Date(2012, 12, 21)
    dt = DatesPlus.DateTime(2012, 12, 21, 16, 30, 20, 200)
    @test trunc(dt, DatesPlus.Year) == DatesPlus.DateTime(2012)
    @test trunc(dt, DatesPlus.Quarter) == DatesPlus.DateTime(2012, 10)
    @test trunc(dt, DatesPlus.Month) == DatesPlus.DateTime(2012, 12)
    @test trunc(dt, DatesPlus.Day) == DatesPlus.DateTime(2012, 12, 21)
    @test trunc(dt, DatesPlus.Hour) == DatesPlus.DateTime(2012, 12, 21, 16)
    @test trunc(dt, DatesPlus.Minute) == DatesPlus.DateTime(2012, 12, 21, 16, 30)
    @test trunc(dt, DatesPlus.Second) == DatesPlus.DateTime(2012, 12, 21, 16, 30, 20)
    @test trunc(dt, DatesPlus.Millisecond) == DatesPlus.DateTime(2012, 12, 21, 16, 30, 20, 200)
    t = DatesPlus.Time(1, 2, 3, 4, 5, 6)
    @test trunc(t, DatesPlus.Hour) == DatesPlus.Time(1)
    @test trunc(t, DatesPlus.Minute) == DatesPlus.Time(1, 2)
    @test trunc(t, DatesPlus.Second) == DatesPlus.Time(1, 2, 3)
    @test trunc(t, DatesPlus.Millisecond) == DatesPlus.Time(1, 2, 3, 4)
    @test trunc(t, DatesPlus.Microsecond) == DatesPlus.Time(1, 2, 3, 4, 5)
    @test trunc(t, DatesPlus.Nanosecond) == DatesPlus.Time(1, 2, 3, 4, 5, 6)
end

Jan = DatesPlus.DateTime(2013, 1, 1) #Tuesday
Feb = DatesPlus.DateTime(2013, 2, 2) #Saturday
Mar = DatesPlus.DateTime(2013, 3, 3) #Sunday
Apr = DatesPlus.DateTime(2013, 4, 4) #Thursday
May = DatesPlus.DateTime(2013, 5, 5) #Sunday
Jun = DatesPlus.DateTime(2013, 6, 7) #Friday
Jul = DatesPlus.DateTime(2013, 7, 7) #Sunday
Aug = DatesPlus.DateTime(2013, 8, 8) #Thursday
Sep = DatesPlus.DateTime(2013, 9, 9) #Monday
Oct = DatesPlus.DateTime(2013, 10, 10) #Thursday
Nov = DatesPlus.DateTime(2013, 11, 11) #Monday
Dec = DatesPlus.DateTime(2013, 12, 11) #Wednesday

@testset "Date functions" begin
    @testset "lastdayofmonth" begin
        @test DatesPlus.lastdayofmonth(Jan) == DatesPlus.DateTime(2013, 1, 31)
        @test DatesPlus.lastdayofmonth(Feb) == DatesPlus.DateTime(2013, 2, 28)
        @test DatesPlus.lastdayofmonth(Mar) == DatesPlus.DateTime(2013, 3, 31)
        @test DatesPlus.lastdayofmonth(Apr) == DatesPlus.DateTime(2013, 4, 30)
        @test DatesPlus.lastdayofmonth(May) == DatesPlus.DateTime(2013, 5, 31)
        @test DatesPlus.lastdayofmonth(Jun) == DatesPlus.DateTime(2013, 6, 30)
        @test DatesPlus.lastdayofmonth(Jul) == DatesPlus.DateTime(2013, 7, 31)
        @test DatesPlus.lastdayofmonth(Aug) == DatesPlus.DateTime(2013, 8, 31)
        @test DatesPlus.lastdayofmonth(Sep) == DatesPlus.DateTime(2013, 9, 30)
        @test DatesPlus.lastdayofmonth(Oct) == DatesPlus.DateTime(2013, 10, 31)
        @test DatesPlus.lastdayofmonth(Nov) == DatesPlus.DateTime(2013, 11, 30)
        @test DatesPlus.lastdayofmonth(Dec) == DatesPlus.DateTime(2013, 12, 31)

        @test DatesPlus.lastdayofmonth(Date(Jan)) == DatesPlus.Date(2013, 1, 31)
        @test DatesPlus.lastdayofmonth(Date(Feb)) == DatesPlus.Date(2013, 2, 28)
        @test DatesPlus.lastdayofmonth(Date(Mar)) == DatesPlus.Date(2013, 3, 31)
        @test DatesPlus.lastdayofmonth(Date(Apr)) == DatesPlus.Date(2013, 4, 30)
        @test DatesPlus.lastdayofmonth(Date(May)) == DatesPlus.Date(2013, 5, 31)
        @test DatesPlus.lastdayofmonth(Date(Jun)) == DatesPlus.Date(2013, 6, 30)
        @test DatesPlus.lastdayofmonth(Date(Jul)) == DatesPlus.Date(2013, 7, 31)
        @test DatesPlus.lastdayofmonth(Date(Aug)) == DatesPlus.Date(2013, 8, 31)
        @test DatesPlus.lastdayofmonth(Date(Sep)) == DatesPlus.Date(2013, 9, 30)
        @test DatesPlus.lastdayofmonth(Date(Oct)) == DatesPlus.Date(2013, 10, 31)
        @test DatesPlus.lastdayofmonth(Date(Nov)) == DatesPlus.Date(2013, 11, 30)
        @test DatesPlus.lastdayofmonth(Date(Dec)) == DatesPlus.Date(2013, 12, 31)
    end
    @testset "firstdayofmonth" begin
        @test DatesPlus.firstdayofmonth(Jan) == DatesPlus.DateTime(2013, 1, 1)
        @test DatesPlus.firstdayofmonth(Feb) == DatesPlus.DateTime(2013, 2, 1)
        @test DatesPlus.firstdayofmonth(Mar) == DatesPlus.DateTime(2013, 3, 1)
        @test DatesPlus.firstdayofmonth(Apr) == DatesPlus.DateTime(2013, 4, 1)
        @test DatesPlus.firstdayofmonth(May) == DatesPlus.DateTime(2013, 5, 1)
        @test DatesPlus.firstdayofmonth(Jun) == DatesPlus.DateTime(2013, 6, 1)
        @test DatesPlus.firstdayofmonth(Jul) == DatesPlus.DateTime(2013, 7, 1)
        @test DatesPlus.firstdayofmonth(Aug) == DatesPlus.DateTime(2013, 8, 1)
        @test DatesPlus.firstdayofmonth(Sep) == DatesPlus.DateTime(2013, 9, 1)
        @test DatesPlus.firstdayofmonth(Oct) == DatesPlus.DateTime(2013, 10, 1)
        @test DatesPlus.firstdayofmonth(Nov) == DatesPlus.DateTime(2013, 11, 1)
        @test DatesPlus.firstdayofmonth(Dec) == DatesPlus.DateTime(2013, 12, 1)

        @test DatesPlus.firstdayofmonth(Date(Jan)) == DatesPlus.Date(2013, 1, 1)
        @test DatesPlus.firstdayofmonth(Date(Feb)) == DatesPlus.Date(2013, 2, 1)
        @test DatesPlus.firstdayofmonth(Date(Mar)) == DatesPlus.Date(2013, 3, 1)
        @test DatesPlus.firstdayofmonth(Date(Apr)) == DatesPlus.Date(2013, 4, 1)
        @test DatesPlus.firstdayofmonth(Date(May)) == DatesPlus.Date(2013, 5, 1)
        @test DatesPlus.firstdayofmonth(Date(Jun)) == DatesPlus.Date(2013, 6, 1)
        @test DatesPlus.firstdayofmonth(Date(Jul)) == DatesPlus.Date(2013, 7, 1)
        @test DatesPlus.firstdayofmonth(Date(Aug)) == DatesPlus.Date(2013, 8, 1)
        @test DatesPlus.firstdayofmonth(Date(Sep)) == DatesPlus.Date(2013, 9, 1)
        @test DatesPlus.firstdayofmonth(Date(Oct)) == DatesPlus.Date(2013, 10, 1)
        @test DatesPlus.firstdayofmonth(Date(Nov)) == DatesPlus.Date(2013, 11, 1)
        @test DatesPlus.firstdayofmonth(Date(Dec)) == DatesPlus.Date(2013, 12, 1)
    end
    @testset "firstdayofweek" begin
        # 2014-01-06 is a Monday = 1st day of week
        a = DatesPlus.Date(2014, 1, 6)
        b = DatesPlus.Date(2014, 1, 7)
        c = DatesPlus.Date(2014, 1, 8)
        d = DatesPlus.Date(2014, 1, 9)
        e = DatesPlus.Date(2014, 1, 10)
        f = DatesPlus.Date(2014, 1, 11)
        g = DatesPlus.Date(2014, 1, 12)
        @test DatesPlus.firstdayofweek(a) == a
        @test DatesPlus.firstdayofweek(b) == a
        @test DatesPlus.firstdayofweek(c) == a
        @test DatesPlus.firstdayofweek(d) == a
        @test DatesPlus.firstdayofweek(e) == a
        @test DatesPlus.firstdayofweek(f) == a
        @test DatesPlus.firstdayofweek(g) == a
        # Test firstdayofweek over the course of the year
        dt = a
        for i = 0:364
            @test DatesPlus.firstdayofweek(dt) == a + DatesPlus.Week(div(i, 7))
            dt += DatesPlus.Day(1)
        end
        a = DatesPlus.DateTime(2014, 1, 6)
        b = DatesPlus.DateTime(2014, 1, 7)
        c = DatesPlus.DateTime(2014, 1, 8)
        d = DatesPlus.DateTime(2014, 1, 9)
        e = DatesPlus.DateTime(2014, 1, 10)
        f = DatesPlus.DateTime(2014, 1, 11)
        g = DatesPlus.DateTime(2014, 1, 12)
        @test DatesPlus.firstdayofweek(a) == a
        @test DatesPlus.firstdayofweek(b) == a
        @test DatesPlus.firstdayofweek(c) == a
        @test DatesPlus.firstdayofweek(d) == a
        @test DatesPlus.firstdayofweek(e) == a
        @test DatesPlus.firstdayofweek(f) == a
        @test DatesPlus.firstdayofweek(g) == a
        dt = a
        for i = 0:364
            @test DatesPlus.firstdayofweek(dt) == a + DatesPlus.Week(div(i, 7))
            dt += DatesPlus.Day(1)
        end
        @test DatesPlus.firstdayofweek(DatesPlus.DateTime(2013, 12, 24)) == DatesPlus.DateTime(2013, 12, 23)
    end
    @testset "lastdayofweek" begin
        # Sunday = last day of week
        # 2014-01-12 is a Sunday
        a = DatesPlus.Date(2014, 1, 6)
        b = DatesPlus.Date(2014, 1, 7)
        c = DatesPlus.Date(2014, 1, 8)
        d = DatesPlus.Date(2014, 1, 9)
        e = DatesPlus.Date(2014, 1, 10)
        f = DatesPlus.Date(2014, 1, 11)
        g = DatesPlus.Date(2014, 1, 12)
        @test DatesPlus.lastdayofweek(a) == g
        @test DatesPlus.lastdayofweek(b) == g
        @test DatesPlus.lastdayofweek(c) == g
        @test DatesPlus.lastdayofweek(d) == g
        @test DatesPlus.lastdayofweek(e) == g
        @test DatesPlus.lastdayofweek(f) == g
        @test DatesPlus.lastdayofweek(g) == g
        dt = a
        for i = 0:364
            @test DatesPlus.lastdayofweek(dt) == g + DatesPlus.Week(div(i, 7))
            dt += DatesPlus.Day(1)
        end
        a = DatesPlus.DateTime(2014, 1, 6)
        b = DatesPlus.DateTime(2014, 1, 7)
        c = DatesPlus.DateTime(2014, 1, 8)
        d = DatesPlus.DateTime(2014, 1, 9)
        e = DatesPlus.DateTime(2014, 1, 10)
        f = DatesPlus.DateTime(2014, 1, 11)
        g = DatesPlus.DateTime(2014, 1, 12)
        @test DatesPlus.lastdayofweek(a) == g
        @test DatesPlus.lastdayofweek(b) == g
        @test DatesPlus.lastdayofweek(c) == g
        @test DatesPlus.lastdayofweek(d) == g
        @test DatesPlus.lastdayofweek(e) == g
        @test DatesPlus.lastdayofweek(f) == g
        @test DatesPlus.lastdayofweek(g) == g
        dt = a
        for i = 0:364
            @test DatesPlus.lastdayofweek(dt) == g + DatesPlus.Week(div(i, 7))
            dt += DatesPlus.Day(1)
        end
        @test DatesPlus.lastdayofweek(DatesPlus.DateTime(2013, 12, 24)) == DatesPlus.DateTime(2013, 12, 29)
    end
    @testset "first/lastdayofquarter" begin
        @test DatesPlus.firstdayofquarter(DatesPlus.Date(2014, 2, 2)) == DatesPlus.Date(2014, 1, 1)
        @test DatesPlus.firstdayofquarter(DatesPlus.Date(2014, 5, 2)) == DatesPlus.Date(2014, 4, 1)
        @test DatesPlus.firstdayofquarter(DatesPlus.Date(2014, 8, 2)) == DatesPlus.Date(2014, 7, 1)
        @test DatesPlus.firstdayofquarter(DatesPlus.Date(2014, 12, 2)) == DatesPlus.Date(2014, 10, 1)

        @test DatesPlus.firstdayofquarter(DatesPlus.DateTime(2014, 2, 2)) == DatesPlus.DateTime(2014, 1, 1)
        @test DatesPlus.firstdayofquarter(DatesPlus.DateTime(2014, 5, 2)) == DatesPlus.DateTime(2014, 4, 1)
        @test DatesPlus.firstdayofquarter(DatesPlus.DateTime(2014, 8, 2)) == DatesPlus.DateTime(2014, 7, 1)
        @test DatesPlus.firstdayofquarter(DatesPlus.DateTime(2014, 12, 2)) == DatesPlus.DateTime(2014, 10, 1)

        @test DatesPlus.lastdayofquarter(DatesPlus.Date(2014, 2, 2)) == DatesPlus.Date(2014, 3, 31)
        @test DatesPlus.lastdayofquarter(DatesPlus.Date(2014, 5, 2)) == DatesPlus.Date(2014, 6, 30)
        @test DatesPlus.lastdayofquarter(DatesPlus.Date(2014, 8, 2)) == DatesPlus.Date(2014, 9, 30)
        @test DatesPlus.lastdayofquarter(DatesPlus.Date(2014, 12, 2)) == DatesPlus.Date(2014, 12, 31)

        @test DatesPlus.lastdayofquarter(DatesPlus.DateTime(2014, 2, 2)) == DatesPlus.DateTime(2014, 3, 31)
        @test DatesPlus.lastdayofquarter(DatesPlus.DateTime(2014, 5, 2)) == DatesPlus.DateTime(2014, 6, 30)
        @test DatesPlus.lastdayofquarter(DatesPlus.DateTime(2014, 8, 2)) == DatesPlus.DateTime(2014, 9, 30)
        @test DatesPlus.lastdayofquarter(DatesPlus.DateTime(2014, 12, 2)) == DatesPlus.DateTime(2014, 12, 31)
    end
    @testset "first/lastdayofyear" begin
        firstday = DatesPlus.Date(2014, 1, 1)
        lastday = DatesPlus.Date(2014, 12, 31)
        for i = 0:364
            dt = firstday + DatesPlus.Day(i)

            @test DatesPlus.firstdayofyear(dt) == firstday
            @test DatesPlus.firstdayofyear(DateTime(dt)) == DateTime(firstday)

            @test DatesPlus.lastdayofyear(dt) == lastday
            @test DatesPlus.lastdayofyear(DateTime(dt)) == DateTime(lastday)
        end
    end
end
# Adjusters
@testset "Adjuster Constructors" begin
    @test DatesPlus.Date(DatesPlus.ismonday, 2014) == DatesPlus.Date(2014, 1, 6)
    @test DatesPlus.Date(DatesPlus.ismonday, 2014, 5) == DatesPlus.Date(2014, 5, 5)

    @test DatesPlus.DateTime(DatesPlus.ismonday, 2014) == DatesPlus.DateTime(2014, 1, 6)
    @test DatesPlus.DateTime(DatesPlus.ismonday, 2014, 5) == DatesPlus.DateTime(2014, 5, 5)
    @test DatesPlus.DateTime(x->DatesPlus.hour(x)==12, 2014, 5, 21) == DatesPlus.DateTime(2014, 5, 21, 12)
    @test DatesPlus.DateTime(x->DatesPlus.minute(x)==30, 2014, 5, 21, 12) == DatesPlus.DateTime(2014, 5, 21, 12, 30)
    @test DatesPlus.DateTime(x->DatesPlus.second(x)==30, 2014, 5, 21, 12, 30) == DatesPlus.DateTime(2014, 5, 21, 12, 30, 30)
    @test DatesPlus.DateTime(x->DatesPlus.millisecond(x)==500, 2014, 5, 21, 12, 30, 30) == DatesPlus.DateTime(2014, 5, 21, 12, 30, 30, 500)
end
@testset "tonext, toprev, tofirst, tolast" begin
    dt = DatesPlus.Date(2014, 5, 21)
    @test DatesPlus.tonext(dt, DatesPlus.Wed) == DatesPlus.Date(2014, 5, 28)
    @test DatesPlus.tonext(dt, DatesPlus.Wed; same=true) == dt
    @test DatesPlus.tonext(dt, DatesPlus.Thu) == DatesPlus.Date(2014, 5, 22)
    @test DatesPlus.tonext(dt, DatesPlus.Fri) == DatesPlus.Date(2014, 5, 23)
    @test DatesPlus.tonext(dt, DatesPlus.Sat) == DatesPlus.Date(2014, 5, 24)
    @test DatesPlus.tonext(dt, DatesPlus.Sun) == DatesPlus.Date(2014, 5, 25)
    @test DatesPlus.tonext(dt, DatesPlus.Mon) == DatesPlus.Date(2014, 5, 26)
    @test DatesPlus.tonext(dt, DatesPlus.Tue) == DatesPlus.Date(2014, 5, 27)
    # No dayofweek function for out of range values
    @test_throws KeyError DatesPlus.tonext(dt, 8)

    @test DatesPlus.tonext(DatesPlus.Date(0), DatesPlus.Mon) == DatesPlus.Date(0, 1, 3)
    @testset "func, diff steps, same" begin
        @test DatesPlus.tonext(DatesPlus.iswednesday, dt) == DatesPlus.Date(2014, 5, 28)
        @test DatesPlus.tonext(DatesPlus.iswednesday, dt; same=true) == dt
        @test DatesPlus.tonext(DatesPlus.isthursday, dt) == DatesPlus.Date(2014, 5, 22)
        @test DatesPlus.tonext(DatesPlus.isfriday, dt) == DatesPlus.Date(2014, 5, 23)
        @test DatesPlus.tonext(DatesPlus.issaturday, dt) == DatesPlus.Date(2014, 5, 24)
        @test DatesPlus.tonext(DatesPlus.issunday, dt) == DatesPlus.Date(2014, 5, 25)
        @test DatesPlus.tonext(DatesPlus.ismonday, dt) == DatesPlus.Date(2014, 5, 26)
        @test DatesPlus.tonext(DatesPlus.istuesday, dt) == DatesPlus.Date(2014, 5, 27)
        @test DatesPlus.tonext(DatesPlus.ismonday, DatesPlus.Date(0)) == DatesPlus.Date(0, 1, 3)

        @test DatesPlus.tonext(!DatesPlus.iswednesday, dt) == DatesPlus.Date(2014, 5, 22)
        @test DatesPlus.tonext(!DatesPlus.isthursday, dt) == DatesPlus.Date(2014, 5, 23)
    end

    # Reach adjust limit
    @test_throws ArgumentError DatesPlus.tonext(DatesPlus.iswednesday, dt; limit=6)

    @test DatesPlus.tonext(DatesPlus.iswednesday, dt;step=DatesPlus.Day(2)) == DatesPlus.Date(2014, 6, 4)
    @test DatesPlus.tonext(DatesPlus.iswednesday, dt;step=DatesPlus.Day(3)) == DatesPlus.Date(2014, 6, 11)
    @test DatesPlus.tonext(DatesPlus.iswednesday, dt;step=DatesPlus.Day(4)) == DatesPlus.Date(2014, 6, 18)
    @test DatesPlus.tonext(DatesPlus.iswednesday, dt;step=DatesPlus.Day(5)) == DatesPlus.Date(2014, 6, 25)
    @test DatesPlus.tonext(DatesPlus.iswednesday, dt;step=DatesPlus.Day(6)) == DatesPlus.Date(2014, 7, 2)
    @test DatesPlus.tonext(DatesPlus.iswednesday, dt;step=DatesPlus.Day(7)) == DatesPlus.Date(2014, 5, 28)
    @test DatesPlus.tonext(DatesPlus.iswednesday, dt;step=DatesPlus.Week(1)) == DatesPlus.Date(2014, 5, 28)
    @test DatesPlus.tonext(DatesPlus.iswednesday, dt;step=DatesPlus.Week(2)) == DatesPlus.Date(2014, 6, 4)
    @test DatesPlus.tonext(DatesPlus.iswednesday, dt;step=DatesPlus.Week(3)) == DatesPlus.Date(2014, 6, 11)
    @test DatesPlus.tonext(DatesPlus.iswednesday, dt;step=DatesPlus.Week(4)) == DatesPlus.Date(2014, 6, 18)
    @test DatesPlus.tonext(DatesPlus.iswednesday, dt;step=DatesPlus.Week(5)) == DatesPlus.Date(2014, 6, 25)
    @test DatesPlus.tonext(DatesPlus.iswednesday, dt;step=DatesPlus.Week(6)) == DatesPlus.Date(2014, 7, 2)

    @test DatesPlus.tonext(DatesPlus.iswednesday, dt; same=true) == dt
    @test DatesPlus.tonext(DatesPlus.isthursday, dt) == DatesPlus.Date(2014, 5, 22)
end

@testset "toprev" begin
    dt = DatesPlus.Date(2014, 5, 21)
    @test DatesPlus.toprev(dt, DatesPlus.Wed) == DatesPlus.Date(2014, 5, 14)
    @test DatesPlus.toprev(dt, DatesPlus.Wed; same=true) == dt
    @test DatesPlus.toprev(dt, DatesPlus.Thu) == DatesPlus.Date(2014, 5, 15)
    @test DatesPlus.toprev(dt, DatesPlus.Fri) == DatesPlus.Date(2014, 5, 16)
    @test DatesPlus.toprev(dt, DatesPlus.Sat) == DatesPlus.Date(2014, 5, 17)
    @test DatesPlus.toprev(dt, DatesPlus.Sun) == DatesPlus.Date(2014, 5, 18)
    @test DatesPlus.toprev(dt, DatesPlus.Mon) == DatesPlus.Date(2014, 5, 19)
    @test DatesPlus.toprev(dt, DatesPlus.Tue) == DatesPlus.Date(2014, 5, 20)
    # No dayofweek function for out of range values
    @test_throws KeyError DatesPlus.toprev(dt, 8)

    @test DatesPlus.toprev(DatesPlus.Date(0), DatesPlus.Mon) == DatesPlus.Date(-1, 12, 27)
end
@testset "tofirst" begin
    dt = DatesPlus.Date(2014, 5, 21)
    @test DatesPlus.tofirst(dt, DatesPlus.Mon) == DatesPlus.Date(2014, 5, 5)
    @test DatesPlus.tofirst(dt, DatesPlus.Tue) == DatesPlus.Date(2014, 5, 6)
    @test DatesPlus.tofirst(dt, DatesPlus.Wed) == DatesPlus.Date(2014, 5, 7)
    @test DatesPlus.tofirst(dt, DatesPlus.Thu) == DatesPlus.Date(2014, 5, 1)
    @test DatesPlus.tofirst(dt, DatesPlus.Fri) == DatesPlus.Date(2014, 5, 2)
    @test DatesPlus.tofirst(dt, DatesPlus.Sat) == DatesPlus.Date(2014, 5, 3)
    @test DatesPlus.tofirst(dt, DatesPlus.Sun) == DatesPlus.Date(2014, 5, 4)

    @test DatesPlus.tofirst(dt, DatesPlus.Mon, of=DatesPlus.Year) == DatesPlus.Date(2014, 1, 6)
    @test DatesPlus.tofirst(dt, DatesPlus.Tue, of=DatesPlus.Year) == DatesPlus.Date(2014, 1, 7)
    @test DatesPlus.tofirst(dt, DatesPlus.Wed, of=DatesPlus.Year) == DatesPlus.Date(2014, 1, 1)
    @test DatesPlus.tofirst(dt, DatesPlus.Thu, of=DatesPlus.Year) == DatesPlus.Date(2014, 1, 2)
    @test DatesPlus.tofirst(dt, DatesPlus.Fri, of=DatesPlus.Year) == DatesPlus.Date(2014, 1, 3)
    @test DatesPlus.tofirst(dt, DatesPlus.Sat, of=DatesPlus.Year) == DatesPlus.Date(2014, 1, 4)
    @test DatesPlus.tofirst(dt, DatesPlus.Sun, of=DatesPlus.Year) == DatesPlus.Date(2014, 1, 5)

    @test DatesPlus.tofirst(DatesPlus.Date(0), DatesPlus.Mon) == DatesPlus.Date(0, 1, 3)
end
@testset "tolast" begin
    dt = DatesPlus.Date(2014, 5, 21)
    @test DatesPlus.tolast(dt, DatesPlus.Mon) == DatesPlus.Date(2014, 5, 26)
    @test DatesPlus.tolast(dt, DatesPlus.Tue) == DatesPlus.Date(2014, 5, 27)
    @test DatesPlus.tolast(dt, DatesPlus.Wed) == DatesPlus.Date(2014, 5, 28)
    @test DatesPlus.tolast(dt, DatesPlus.Thu) == DatesPlus.Date(2014, 5, 29)
    @test DatesPlus.tolast(dt, DatesPlus.Fri) == DatesPlus.Date(2014, 5, 30)
    @test DatesPlus.tolast(dt, DatesPlus.Sat) == DatesPlus.Date(2014, 5, 31)
    @test DatesPlus.tolast(dt, DatesPlus.Sun) == DatesPlus.Date(2014, 5, 25)

    @test DatesPlus.tolast(dt, DatesPlus.Mon, of=DatesPlus.Year) == DatesPlus.Date(2014, 12, 29)
    @test DatesPlus.tolast(dt, DatesPlus.Tue, of=DatesPlus.Year) == DatesPlus.Date(2014, 12, 30)
    @test DatesPlus.tolast(dt, DatesPlus.Wed, of=DatesPlus.Year) == DatesPlus.Date(2014, 12, 31)
    @test DatesPlus.tolast(dt, DatesPlus.Thu, of=DatesPlus.Year) == DatesPlus.Date(2014, 12, 25)
    @test DatesPlus.tolast(dt, DatesPlus.Fri, of=DatesPlus.Year) == DatesPlus.Date(2014, 12, 26)
    @test DatesPlus.tolast(dt, DatesPlus.Sat, of=DatesPlus.Year) == DatesPlus.Date(2014, 12, 27)
    @test DatesPlus.tolast(dt, DatesPlus.Sun, of=DatesPlus.Year) == DatesPlus.Date(2014, 12, 28)

    @test DatesPlus.tolast(DatesPlus.Date(0), DatesPlus.Mon) == DatesPlus.Date(0, 1, 31)
end
@testset "filter" begin # was recur
    startdate = DatesPlus.Date(2014, 1, 1)
    stopdate = DatesPlus.Date(2014, 2, 1)
    @test length(filter(x->true, startdate:DatesPlus.Day(1):stopdate))  == 32
    @test length(filter(x->true, stopdate:DatesPlus.Day(-1):startdate)) == 32

    Januarymondays2014 = [DatesPlus.Date(2014, 1, 6), DatesPlus.Date(2014, 1, 13), DatesPlus.Date(2014, 1, 20), DatesPlus.Date(2014, 1, 27)]
    @test filter(DatesPlus.ismonday, startdate:DatesPlus.Day(1):stopdate) == Januarymondays2014

    @test_throws MethodError filter((x, y)->x + y, DatesPlus.Date(2013):DatesPlus.Day(1):DatesPlus.Date(2014))
    @test_throws MethodError DatesPlus.DateFunction((x, y)->x + y, Date(0))
    @test_throws ArgumentError DatesPlus.DateFunction((dt)->2, Date(0))
    @test length(filter(x->true, DatesPlus.Date(2013):DatesPlus.Day(1):DatesPlus.Date(2013, 2))) == 32
    @test length(filter(x->true, DatesPlus.Date(2013):DatesPlus.Day(1):DatesPlus.Date(2013, 1, 1))) == 1
    @test length(filter(x->true, DatesPlus.Date(2013):DatesPlus.Day(1):DatesPlus.Date(2013, 1, 2))) == 2
    @test length(filter(x->true, DatesPlus.Date(2013):DatesPlus.Day(1):DatesPlus.Date(2013, 1, 3))) == 3
    @test length(filter(x->true, DatesPlus.Date(2013):DatesPlus.Day(1):DatesPlus.Date(2013, 1, 4))) == 4
    @test length(filter(x->true, DatesPlus.Date(2013):DatesPlus.Day(1):DatesPlus.Date(2013, 1, 5))) == 5
    @test length(filter(x->true, DatesPlus.Date(2013):DatesPlus.Day(1):DatesPlus.Date(2013, 1, 6))) == 6
    @test length(filter(x->true, DatesPlus.Date(2013):DatesPlus.Day(1):DatesPlus.Date(2013, 1, 7))) == 7
    @test length(filter(x->true, DatesPlus.Date(2013):DatesPlus.Day(1):DatesPlus.Date(2013, 1, 8))) == 8
    @test length(filter(x->true, DatesPlus.Date(2013):DatesPlus.Month(1):DatesPlus.Date(2013, 1, 1))) == 1
    @test length(filter(x->true, DatesPlus.Date(2013):DatesPlus.Day(-1):DatesPlus.Date(2012, 1, 1))) == 367

    # Empty range
    @test length(filter(x->true, DatesPlus.Date(2013):DatesPlus.Day(1):DatesPlus.Date(2012, 1, 1))) == 0

    # All leap days in 20th century
    @test length(filter(DatesPlus.Date(1900):DatesPlus.Day(1):DatesPlus.Date(2000)) do x
        DatesPlus.month(x) == DatesPlus.Feb && DatesPlus.day(x) == 29
    end) == 24
end

@testset "Holiday test cases" begin
    dt = DatesPlus.Date(2014, 5, 21)
    dr = DatesPlus.Date(2014):DatesPlus.Day(1):DatesPlus.Date(2015)
    @testset "Thanksgiving" begin
        # 4th Thursday of November
        thanksgiving = x->DatesPlus.dayofweek(x) == DatesPlus.Thu &&
                              DatesPlus.month(x) == DatesPlus.Nov &&
                   DatesPlus.dayofweekofmonth(x) == 4

        d = DatesPlus.Date(2014, 6, 5)

        @test DatesPlus.tonext(d) do x
            thanksgiving(x)
        end == DatesPlus.Date(2014, 11, 27)

        @test DatesPlus.toprev(d) do x
            thanksgiving(x)
        end == DatesPlus.Date(2013, 11, 28)
    end
    @testset "Pittsburgh street cleaning" begin
        @test length(filter(dr) do x
            DatesPlus.dayofweek(x) == DatesPlus.Tue &&
            DatesPlus.April < DatesPlus.month(x) < DatesPlus.Nov &&
            DatesPlus.dayofweekofmonth(x) == 2
        end) == 6
    end
    @testset "U.S. Federal Holidays" begin
        newyears(y) = (y, 1, 1)
        independenceday(y) = (y, 7, 4)
        veteransday(y) = (y, 11, 11)
        christmas(y) = (y, 12, 25)

        isnewyears(dt) = DatesPlus.yearmonthday(dt) == newyears(DatesPlus.year(dt))
        isindependenceday(dt) = DatesPlus.yearmonthday(dt) == independenceday(DatesPlus.year(dt))
        isveteransday(dt) = DatesPlus.yearmonthday(dt) == veteransday(DatesPlus.year(dt))
        ischristmas(dt) = DatesPlus.yearmonthday(dt) == christmas(DatesPlus.year(dt))
        ismartinlutherking(dt) = DatesPlus.dayofweek(dt) == DatesPlus.Mon &&
            DatesPlus.month(dt) == DatesPlus.Jan && DatesPlus.dayofweekofmonth(dt) == 3
        ispresidentsday(dt) = DatesPlus.dayofweek(dt) == DatesPlus.Mon &&
            DatesPlus.month(dt) == DatesPlus.Feb && DatesPlus.dayofweekofmonth(dt) == 3
        # Last Monday of May
        ismemorialday(dt) = DatesPlus.dayofweek(dt) == DatesPlus.Mon &&
                            DatesPlus.month(dt) == DatesPlus.May &&
                            DatesPlus.dayofweekofmonth(dt) == DatesPlus.daysofweekinmonth(dt)
        islaborday(dt) = DatesPlus.dayofweek(dt) == DatesPlus.Mon &&
            DatesPlus.month(dt) == DatesPlus.Sep && DatesPlus.dayofweekofmonth(dt) == 1
        iscolumbusday(dt) = DatesPlus.dayofweek(dt) == DatesPlus.Mon &&
            DatesPlus.month(dt) == DatesPlus.Oct && DatesPlus.dayofweekofmonth(dt) == 2
        isthanksgiving(dt) = DatesPlus.dayofweek(dt) == DatesPlus.Thu &&
            DatesPlus.month(dt) == DatesPlus.Nov && DatesPlus.dayofweekofmonth(dt) == 4

        function easter(y)
            # Butcher's Algorithm: http://www.smart.net/~mmontes/butcher.html
            a = y % 19
            b = div(y, 100)
            c = y % 100
            d = div(b, 4)
            e = b % 4
            f = div(b + 8, 25)
            g = div(b - f + 1, 3)
            h = (19 * a + b - d - g + 15) % 30
            i = div(c, 4)
            k = c % 4
            l = (32 + 2 * e + 2 * i - h - k) % 7
            m = div(a + 11 * h + 22 * l, 451)
            month = div(h + l - 7 * m + 114, 31)
            p = (h + l - 7 * m + 114) % 31
            return (y, month, p + 1)
        end
        iseaster(dt) = DatesPlus.yearmonthday(dt) == easter(DatesPlus.year(dt))

        HOLIDAYS = x->isnewyears(x) || isindependenceday(x) ||
                      isveteransday(x) || ischristmas(x) ||
                      ismartinlutherking(x) || ispresidentsday(x) ||
                      ismemorialday(x) || islaborday(x) ||
                      iscolumbusday(x) || isthanksgiving(x)

        @test length(filter(HOLIDAYS, dr)) == 11

        OBSERVEDHOLIDAYS = x->begin
            # If the holiday is on a weekday
            if HOLIDAYS(x) && DatesPlus.dayofweek(x) < DatesPlus.Saturday
                return true
            # Holiday is observed Monday if falls on Sunday
            elseif DatesPlus.dayofweek(x) == 1 && HOLIDAYS(x - DatesPlus.Day(1))
                return true
            # Holiday is observed Friday if falls on Saturday
            elseif DatesPlus.dayofweek(x) == 5 && HOLIDAYS(x + DatesPlus.Day(1))
                return true
            else
                return false
            end
        end

        observed = filter(OBSERVEDHOLIDAYS, DatesPlus.Date(1999):DatesPlus.Day(1):DatesPlus.Date(2000))
        @test length(observed) == 11
        @test observed[10] == DatesPlus.Date(1999, 12, 24)
        @test observed[11] == DatesPlus.Date(1999, 12, 31)

        # Get all business/working days of 2014
        # Since we have already defined observed holidays,
        # we just look at weekend days and negate the result
        @test length(filter(DatesPlus.Date(2014):DatesPlus.Day(1):DatesPlus.Date(2015)) do x
            !(OBSERVEDHOLIDAYS(x) ||
            DatesPlus.dayofweek(x) > 5)
        end) == 251
    end
end

@testset "proof-of-concepts" begin
    # First day of the next month for each day of 2014
    @test length([DatesPlus.firstdayofmonth(i + DatesPlus.Month(1))
        for i in DatesPlus.Date(2014):DatesPlus.Day(1):DatesPlus.Date(2014, 12, 31)]) == 365

    # From those goofy email forwards claiming a "special, lucky month"
    # that has 5 Fridays, 5 Saturdays, and 5 Sundays and that it only
    # occurs every 823 years
    @test length(filter(Date(2000):DatesPlus.Month(1):Date(2016)) do dt
        sum = 0
        for i = 1:7
            sum += DatesPlus.dayofweek(dt) > 4 ? DatesPlus.daysofweekinmonth(dt) : 0
            dt += DatesPlus.Day(1)
        end
        return sum == 15
    end) == 15 # On average, there's one of those months every year

    r = DatesPlus.Time(x->DatesPlus.second(x) == 5, 1)
    @test r == DatesPlus.Time(1, 0, 5)

    r = filter(x->DatesPlus.second(x) == 5, DatesPlus.Time(0):DatesPlus.Second(1):DatesPlus.Time(10))
    @test length(r) == 600
    @test first(r) == DatesPlus.Time(0, 0, 5)
    @test last(r) == DatesPlus.Time(9, 59, 5)
end
end # module


text:
;; --------------------------------------------------------------------
;;	text_update.jl
;;
;;						Mar/02/2011
;;
;; --------------------------------------------------------------------
(require 'tables)
(load-file "/var/www/data_base/common/rep_common/text_manipulate.jl")
;; --------------------------------------------------------------------
(format t "*** 開始 ***")
(setq file_txt (cadr command-line-args))
(setq id_in (caddr command-line-args))
(setq population_in (cadddr command-line-args))
(format t file_txt)
(format t "%1$s%2$s" id_in population_in)
;
(setq dict_aa (make-table string-hash string=))
(let ((port_in (open-file file_txt 'read)))
	(text_read_proc port_in dict_aa)
)
;
(if (table-bound-p dict_aa id_in)
	(dict_update_proc dict_aa id_in population_in)
)
;
(table-walk dict_display_proc_single dict_aa)
;
(text_write_proc file_txt dict_aa)
;
(format t "*** 終了 ***")

;; --------------------------------------------------------------------

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
function isarmstrong(number)
    d = digits(number)
    return sum(d .^ length(d)) == number
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
"""
$(TYPEDEF)
Encapsulation structure for a DFGNode (Variable or Factor) in Graphs.jl graph.
"""
mutable struct GraphsNode
    index::Int
    dfgNode::DFGNode
end
const FGType = Graphs.GenericIncidenceList{GraphsNode,Graphs.Edge{GraphsNode},Dict{Int,GraphsNode},Dict{Int,Array{Graphs.Edge{GraphsNode},1}}}

mutable struct GraphsDFG{T <: AbstractParams} <: AbstractDFG
    g::FGType
    description::String
    userId::String
    robotId::String
    sessionId::String
    userData::Dict{Symbol, String}
    robotData::Dict{Symbol, String}
    sessionData::Dict{Symbol, String}
    nodeCounter::Int64
    labelDict::Dict{Symbol, Int64}
    addHistory::Vector{Symbol} #TODO: Discuss more - is this an audit trail?
    solverParams::T # Solver parameters
end

GraphsDFG{T}(   g::FGType=Graphs.incdict(GraphsNode,is_directed=false),
                d::String="Graphs.jl implementation",
                n::Int64=0,
                l::Dict{Symbol, Int64}=Dict{Symbol, Int64}(),
                a::Vector{Symbol}=Symbol[];
                userId::String = "UserID",
                robotId::String = "robotID",
                sessionId::String = "sessionID",
                userData::Dict{Symbol, String} = Dict{Symbol, String}(),
                robotData::Dict{Symbol, String} = Dict{Symbol, String}(),
                sessionData::Dict{Symbol, String} = Dict{Symbol, String}(),
                params::T=NoSolverParams()) where T <: AbstractParams = GraphsDFG{T}(g, d, userId, robotId, sessionId, userData, robotData, sessionData, n, l, a, params)


text:
instance(instance::Instance)::Instance = instance
instance(handle::Handle)::Instance = instance(parent(handle))


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct NeuralODE <: DynamicInjection
        name::String
        n_layers::Int
        n_neurons::Int
        n_feedback_states::Int
        nn_parameters::Vector{Float64}
        activation_function::String
        initialize_to_zero::Bool
        base_power::Float64
        states::Vector{Symbol}
        n_states::Int
        ext::Dict{String, Any}
        internal::InfrastructureSystemsInternal
    end

Densely connected NODE with variable number of feedback states.

# Arguments
- `name::String`
- `n_layers::Int`: Number of hidden layers
- `n_neurons::Int`: Neurons per hidden layer
- `n_feedback_states::Int`: Number of feedback states. The total number of states equals n_feedback_states + 2 for real and imaginary output current.
- `nn_parameters::Vector{Float64}`: Weights and biases of the nn which defines the dynamics of the component
- `activation_function::String`: Activation function for input and hidden layers. Output layer has identity activation.
- `initialize_to_zero::Bool`: If true, initialize the output current to zero. Assumes the NODE was trained and will be used with another device at the same bus which accounts for steady state current to match the PF solution.
- `base_power::Float64`: Base power
- `states::Vector{Symbol}`: Real output current, imaginary output current
- `n_states::Int`
- `ext::Dict{String, Any}`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct NeuralODE <: DynamicInjection
    name::String
    "Number of hidden layers"
    n_layers::Int
    "Neurons per hidden layer"
    n_neurons::Int
    "Number of feedback states. The total number of states equals n_feedback_states + 2 for real and imaginary output current."
    n_feedback_states::Int
    "Weights and biases of the nn which defines the dynamics of the component"
    nn_parameters::Vector{Float64}
    "Activation function for input and hidden layers. Output layer has identity activation."
    activation_function::String
    "If true, initialize the output current to zero. Assumes the NODE was trained and will be used with another device at the same bus which accounts for steady state current to match the PF solution."
    initialize_to_zero::Bool
    "Base power"
    base_power::Float64
    "Real output current, imaginary output current"
    states::Vector{Symbol}
    n_states::Int
    ext::Dict{String, Any}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function NeuralODE(name, n_layers, n_neurons, n_feedback_states, nn_parameters, activation_function, initialize_to_zero, base_power=100.0, ext=Dict{String, Any}(), )
    NeuralODE(name, n_layers, n_neurons, n_feedback_states, nn_parameters, activation_function, initialize_to_zero, base_power, ext, [:Ir, :Ii], 2, InfrastructureSystemsInternal(), )
end

function NeuralODE(; name, n_layers, n_neurons, n_feedback_states, nn_parameters, activation_function, initialize_to_zero, base_power=100.0, states=[:Ir, :Ii], n_states=2, ext=Dict{String, Any}(), internal=InfrastructureSystemsInternal(), )
    NeuralODE(name, n_layers, n_neurons, n_feedback_states, nn_parameters, activation_function, initialize_to_zero, base_power, states, n_states, ext, internal, )
end

# Constructor for demo purposes; non-functional.
function NeuralODE(::Nothing)
    NeuralODE(;
        name="init",
        n_layers=0,
        n_neurons=0,
        n_feedback_states=0,
        nn_parameters=Any[0],
        activation_function="0",
        initialize_to_zero=,
        base_power=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`NeuralODE`](@ref) `name`."""
get_name(value::NeuralODE) = value.name
"""Get [`NeuralODE`](@ref) `n_layers`."""
get_n_layers(value::NeuralODE) = value.n_layers
"""Get [`NeuralODE`](@ref) `n_neurons`."""
get_n_neurons(value::NeuralODE) = value.n_neurons
"""Get [`NeuralODE`](@ref) `n_feedback_states`."""
get_n_feedback_states(value::NeuralODE) = value.n_feedback_states
"""Get [`NeuralODE`](@ref) `nn_parameters`."""
get_nn_parameters(value::NeuralODE) = value.nn_parameters
"""Get [`NeuralODE`](@ref) `activation_function`."""
get_activation_function(value::NeuralODE) = value.activation_function
"""Get [`NeuralODE`](@ref) `initialize_to_zero`."""
get_initialize_to_zero(value::NeuralODE) = value.initialize_to_zero
"""Get [`NeuralODE`](@ref) `base_power`."""
get_base_power(value::NeuralODE) = value.base_power
"""Get [`NeuralODE`](@ref) `states`."""
get_states(value::NeuralODE) = value.states
"""Get [`NeuralODE`](@ref) `n_states`."""
get_n_states(value::NeuralODE) = value.n_states
"""Get [`NeuralODE`](@ref) `ext`."""
get_ext(value::NeuralODE) = value.ext
"""Get [`NeuralODE`](@ref) `internal`."""
get_internal(value::NeuralODE) = value.internal

"""Set [`NeuralODE`](@ref) `n_layers`."""
set_n_layers!(value::NeuralODE, val) = value.n_layers = val
"""Set [`NeuralODE`](@ref) `n_neurons`."""
set_n_neurons!(value::NeuralODE, val) = value.n_neurons = val
"""Set [`NeuralODE`](@ref) `n_feedback_states`."""
set_n_feedback_states!(value::NeuralODE, val) = value.n_feedback_states = val
"""Set [`NeuralODE`](@ref) `nn_parameters`."""
set_nn_parameters!(value::NeuralODE, val) = value.nn_parameters = val
"""Set [`NeuralODE`](@ref) `activation_function`."""
set_activation_function!(value::NeuralODE, val) = value.activation_function = val
"""Set [`NeuralODE`](@ref) `initialize_to_zero`."""
set_initialize_to_zero!(value::NeuralODE, val) = value.initialize_to_zero = val
"""Set [`NeuralODE`](@ref) `base_power`."""
set_base_power!(value::NeuralODE, val) = value.base_power = val
"""Set [`NeuralODE`](@ref) `ext`."""
set_ext!(value::NeuralODE, val) = value.ext = val


text:
using NoiseContrastiveEstimation
using Test, Random
using ReverseDiff, LinearAlgebra
using Statistics

@testset "Basic Functional" begin
    R = MersenneTwister(125)
    lϕ(x, θ) = θ[1] * x + θ[2] * x * x / 2
    gϕ(x, ϕ) = [x, x * x / 2]
    data = randn(R, 200) .* 0.3 .+ 1.5
    noised = foldl(hcat, x .+ randn(R, 20) .* 0.4 for x in data)
    J = CNCE(f = lϕ, grad_f = gϕ, data = data, noised = noised)
    @test J((15, -10)).J < J((7.5, -5)).J < J((0, 0)).J
    results = nesterov(J, zeros(2), 0.9, 20.0)
    x = results.sol
    @test J(x).J < J((15, -10)).J
    Jautodiff = CNCE(; f = lϕ, data, noised)
    @test J(x).J == Jautodiff(x).J
    @test J(x).dJ == Jautodiff(x).dJ
end

@testset "Estimate Covariate Model" begin
    R = MersenneTwister(125)
    gϕ(x, θ) = [x[1] * x[1], x[2] * x[2], x[1] * x[2], x[1], x[2]]
    lϕ(x, θ) = θ'gϕ(x, θ)

    data = [[x + 1, y + 1] for (x, y) in zip(randn(R, 200), randn(R, 200))]
    noised = [x + 0.3 * randn(R, 2) for j = 1:50, x in data]
    J = CNCE(f = lϕ, grad_f = gϕ, data = data, noised = noised)
    results = nesterov(J, zeros(5), 0.9, 2.5)
    @test results.sol[1] < 0
    @test results.sol[2] < 0
    @test results.sol[4] > 0
    @test results.sol[5] > 0
    J = StochasticCNCE(
        f = lϕ,
        grad_f = gϕ,
        data = data,
        perturbator = x -> x + 0.3 * randn(R, 2),
        minibatch = 20,
        K = 5,
        rng = R,
    )
    results = nesterov(J, zeros(5), 0.9, 20.0, maxiter = 5000)
    @test results.sol[1] < 0
    @test results.sol[2] < 0
    @test results.sol[4] > 0
    @test results.sol[5] > 0
end

@testset "RBM training" begin
    R = MersenneTwister(125)
    function gϕ(σ, ξ)
        sc = ξ * σ
        g = similar(ξ)
        for i = 1:size(ξ, 1), j = 1:size(ξ, 2)
            g[i, j] = 2 * σ[j] * sc[i]
        end
        g
    end

    function lϕ(σ, ξ)
        sc = ξ * σ
        sc'sc
    end

    data = [sign.(randn(R, 10)) for i = 1:2]
    noised = [sign.(rand(R, 10) .- 1 / 3) .* data[i] for j = 1:20, i in eachindex(data)]
    J = CNCE(lϕ, gϕ, data, noised)
    x = sign.(randn(R, 4, 10)) .* 0.01
    results = nesterov(J, x, 0.95, 0.05)
    @show results
    dx = sign.(results.sol)
    M = [maximum(abs.(dx * data[i] / length(data[i]))) for i = 1:length(data)]
    @test sum(M) == 2
end

@testset "DBM training" begin
    R = MersenneTwister(125)

    function lϕ(σ, ξ)
        L = length(σ)
        t = L ÷ 3
        sc1 = view(ξ, :, 1:t) * view(σ, 1:t)
        sc2 = view(ξ, :, (t+1):2t) * view(σ, (t+1):2t)
        sc3 = view(ξ, :, (2t+1):L) * view(σ, (2t+1):L)
        dot(sc1 + sc2, sc1 + sc2) + dot(sc2 + sc3, sc2 + sc3)
    end

    Nσ = 20
    Nμ = 4
    Ndata = 2
    Nnoise = 20
    p_flip = 1 / 3

    data = [sign.(randn(R, Nσ)) for i = 1:Ndata]
    noised =
        [sign.(rand(R, Nσ) .- p_flip) .* data[i] for j = 1:Nnoise, i in eachindex(data)]
    J = CNCE(; f = lϕ, data, noised)
    x = sign.(randn(R, Nμ, Nσ)) .* 0.01
    results = nesterov(J, x, 0.9, 0.5)
    dx = sign.(results.sol)
    M = [maximum(abs.(dx * data[i] / length(data[i]))) for i = 1:length(data)]
    @test sum(M) == 2
end


text:
"""
A permutation

# Fields:

- `columns`: The columns being indexed as a vector of integers (column numbers)
- `perm`: the permutation - an array or iterator which has the sorted permutation
"""
struct Perm{X}
    columns::Vector{Int}
    perm::X
end

abstract type AbstractIndexedTable end

"""
A tabular data structure that extends [`Columns`](@ref).  Create an `IndexedTable` with the
[`table`](@ref) function.
"""
struct IndexedTable{C<:Columns} <: AbstractIndexedTable
    # `Columns` object which iterates to give an array of rows
    columns::C
    # columns that are primary keys (Vector{Int})
    pkey::Vector{Int}
    # Cache permutations by various subsets of columns
    perms::Vector{Perm}
    # store what percent of the data in each column is unique
    cardinality::Vector{Union{Float64,Missing}}

    columns_buffer::Any
end

Base.@deprecate_binding NextTable IndexedTable

"""
    table(cols; kw...)

Create a table from a (named) tuple of AbstractVectors.

    table(cols::AbstractVector...; names::Vector{Symbol}, kw...)

Create a table from the provided `cols`, optionally with `names`.

    table(cols::Columns; kw...)

Construct a table from a vector of tuples. See [`rows`](@ref) and [`Columns`](@ref).

    table(t::Union{IndexedTable, NDSparse}; kw...)

Copy a Table or NDSparse to create a new table. The same primary keys as the input are used.

    table(x; kw...)

Create an `IndexedTable` from any object `x` that follows the `Tables.jl` interface.

# Keyword Argument Options:

- `pkey`: select columns to sort by and be the primary key.
- `presorted = false`: is the data pre-sorted by primary key columns?
- `copy = true`: creates a copy of the input vectors if `true`. Irrelevant if `chunks` is specified.
- `chunks::Integer`: distribute the table.  Options are:
    - `Int` -- (number of chunks) a safe bet is `nworkers()` after `using Distributed`.
    - `Vector{Int}` -- Number of elements in each of the `length(chunks)` chunks.

# Examples:

    table(rand(10), rand(10), names = [:x, :y], pkey = :x)

    table(rand(Bool, 20), rand(20), rand(20), pkey = [1,2])

    table((x = 1:10, y = randn(10)))

    table([(1,2), (3,4)])
"""
function table end

function table(::Val{:serial}, cols::Tup;
               pkey=Int[],
               chunks=nothing, # unused here
               perms=Perm[],
               presorted=false,
               copy=true,
               cardinality=fill(missing, length(cols)))

    cs = rows(cols)

    if isa(pkey, Union{Int, Symbol})
        pkey = [pkey]
    elseif isa(pkey, Tuple)
        pkey = [pkey...]
    end

    if !presorted && !isempty(pkey)
        pkeys = rows(cs, (pkey...,))
        if !issorted(pkeys)
            perm = sortperm(pkeys)
            if copy
                cs = cs[perm]
            else
                Base.permute!!(refs(cs), perm)
            end
        elseif copy
            cs = copyto!(similar(cs), cs)
        end
    elseif copy
        cs = copyto!(similar(cs), cs)
    end

    intpkey = map(k->colindex(cs, k), pkey)

    IndexedTable{typeof(cs)}(cs,
           intpkey,
           perms,
           cardinality,
           similar(cs, 0))
end

function table(::Val{impl}, cols; kwargs...) where impl
    if impl == :distributed && isa(cols, Tup)
        error("""You requested to create a distributed table.
                 Distributed table is implemented by JuliaDB.
                 run `using JuliaDB` and try again.""")
    else
        error("unknown table implementation invoked")
    end
end

# detect if a distributed table has to be constructed.
_impl(impl::Val) = impl
_impl(impl::Val, x::AbstractArray, z...) = _impl(impl, z...)
_impl(x::AbstractArray...) = _impl(Val{:serial}(), x...)

# table methods go through here
function table(cs::Tup; chunks=nothing, kwargs...)
    impl = chunks !== nothing ? Val(:distributed) : _impl(astuple(cs)...)
    table(impl, cs; chunks=chunks, kwargs...)
end

table(cs::Columns; kwargs...) = table(columns(cs); kwargs...)
table(c::Columns{<:Pair}; kwargs...) = convert(IndexedTable, columns(c).first, columns(c).second; kwargs...)

function table(cols::AbstractArray...; names=nothing, kwargs...)
    if names !== nothing && all(x -> isa(x, Symbol), names)
        cs = namedtuple(names...)(cols)
    else
        cs = cols
    end
    table(cs; kwargs...)
end

# Easy constructor to create a derivative table
function table(t::IndexedTable;
               columns=t.columns,
               chunks=nothing,
               pkey=t.pkey,
               perms=t.perms,
               cardinality=t.cardinality,
               presorted=false,
               copy=true)

    table(columns,
          pkey=pkey,
          perms=perms,
          chunks=chunks,
          cardinality=cardinality,
          presorted=presorted,
          copy=copy)
end

Base.@pure colnames(t::AbstractIndexedTable) = fieldnames(eltype(t))
columns(t::IndexedTable) = columns(t.columns)
# throw a better error message when a custom array
# of different size is used
function column(t::IndexedTable, a::AbstractArray)
    if length(t) != length(a)
        throw(ArgumentError(
                "vector provided must have the same length as table"
             )
        )
    end
    column(rows(t), a)
end

Base.eltype(::Type{IndexedTable{C}}) where {C} = eltype(C)
Base.eltype(t::IndexedTable) = eltype(typeof(t))
Base.copy(t::IndexedTable) = table(t, copy=true)
function Base.empty!(t::IndexedTable)
    empty!(t.perms)
    empty!(rows(t))
    t
end
Base.:(==)(a::IndexedTable, b::IndexedTable) = rows(a) == rows(b)
Base.isequal(a::IndexedTable, b::IndexedTable) = isequal(rows(a), rows(b))

Base.getindex(t::IndexedTable, i::Integer) = getindex(t.columns, i)
Base.getindex(t::IndexedTable, i::Colon) = copy(t)
Base.lastindex(t::IndexedTable) = length(t)

function Base.view(t::IndexedTable, I)
    t.pkey == Int64[] || eltype(I) == Bool || issorted(I) ||
        throw(ArgumentError("`view` called with unsorted index."))
    table(
        view(t.columns, I),
        pkey = t.pkey,
        copy = false,
        presorted = true)
end

Base.length(t::IndexedTable) = length(t.columns)
Base.iterate(t::IndexedTable, i) = iterate(t.columns, i)
Base.iterate(t::IndexedTable) = iterate(t.columns)
function getindex(t::IndexedTable, idxs::AbstractVector{<:Integer})
    if t.pkey == Int64[] || eltype(idxs) == Bool || issorted(idxs)
       #perms = map(t.perms) do p
       #    # TODO: make the ranks continuous
       #    Perm(p.columns, p.perm[idxs])
       #end
        perms = Perm[]
        table(t, columns=t.columns[idxs], perms=perms, copy=false, presorted=true)
    else
        # this is for gracefully allowing this later
        throw(ArgumentError("`getindex` called with unsorted index. This is not allowed at this time."))
    end
end

function Base.getindex(d::ColDict{<:AbstractIndexedTable}, key::Tuple)
    t = d[]
    idx = [colindex(t, k) for k in key]
    pkey = Int[]
    for (i, pk) in enumerate(t.pkey)
        j = something(findfirst(isequal(pk), idx), 0)
        if j > 0
            push!(pkey, j)
        end
    end
    table(d.src, columns=columns(t, key), pkey=pkey)
end

function Base.getindex(d::ColDict{<:AbstractIndexedTable}, key::SpecialSelector)
    getindex(d, lowerselection(d[], key))
end

function ColDict(t::AbstractIndexedTable; copy=nothing)
    ColDict(Base.copy(t.pkey), t,
            convert(Array{Any}, Base.copy(collect(colnames(t)))),
            Any[columns(t)...], copy)
end

function Base.getindex(d::ColDict{<:AbstractIndexedTable})
    table(d.columns...;
          names=d.names,
          copy=d.copy === nothing ? false : d.copy,
          pkey=d.pkey)
end

function subtable(t::Union{Columns, IndexedTable}, idxs; presorted=true)
    table(t, columns=rows(t)[idxs], perms=t.perms, copy=false, presorted=presorted)
end

function primaryperm(t::IndexedTable)
    Perm(t.pkey, Base.OneTo(length(t)))
end

permcache(t::IndexedTable) = vcat(primaryperm(t), t.perms)
cacheperm!(t::IndexedTable, p) = push!(t.perms, p)

"""
    pkeynames(t::Table)

Names of the primary key columns in `t`.

# Examples

    t = table([1,2], [3,4]);
    pkeynames(t)

    t = table([1,2], [3,4], pkey=1);
    pkeynames(t)

    t = table([2,1],[1,3],[4,5], names=[:x,:y,:z], pkey=(1,2));
    pkeynames(t)
"""
function pkeynames(t::AbstractIndexedTable)
    if eltype(t) <: NamedTuple
        (colnames(t)[t.pkey]...,)
    else
        (t.pkey...,)
    end
end

# for a table, selecting the "value" means selecting all fields
valuenames(t::AbstractIndexedTable) = (colnames(t)...,)

ncols(t::AbstractIndexedTable) = length(colnames(t))

"""
    pkeys(itr::IndexedTable)

Primary keys of the table. If Table doesn't have any designated
primary key columns (constructed without `pkey` argument) then
a default key of tuples `(1,):(n,)` is generated.

# Example

    a = table(["a","b"], [3,4]) # no pkey
    pkeys(a)

    a = table(["a","b"], [3,4], pkey=1)
    pkeys(a)
"""
function pkeys(t::IndexedTable)
    if isempty(t.pkey)
        Columns((Base.OneTo(length(t)),))
    else
        rows(t, pkeynames(t))
    end
end

Base.values(t::IndexedTable) = rows(t)

"""
    sort(t    ; select, kw...)
    sort(t, by; select, kw...)

Sort rows by `by`. All of `Base.sort` keyword arguments can be used.

# Examples

    t=table([1,1,1,2,2,2], [1,1,2,2,1,1], [1,2,3,4,5,6],
    sort(t, :z; select = (:y, :z), rev = true)
"""
sort(t::IndexedTable, by...; select = valuenames(t), kwargs...) =
    table(rows(t, select)[sortperm(rows(t, by...); kwargs...)], copy = false)

"""
    sort!(t    ; kw...)
    sort!(t, by; kw...)

Sort rows of `t` by `by` in place. All of `Base.sort` keyword arguments can be used.

# Examples

    t = table([1,1,1,2,2,2], [1,1,2,2,1,1], [1,2,3,4,5,6], names=[:x,:y,:z]);
    sort!(t, :z, rev = true)
    t
"""
function sort!(t::IndexedTable, by...; kwargs...)
    isempty(t.pkey) || error("Tables with primary keys can't be sorted in place")
    Base.permute!!(refs(rows(t)), sortperm(rows(t, by...); kwargs...))
    t
end

"""
    excludecols(itr, cols) -> Tuple of Int

Names of all columns in `itr` except `cols`. `itr` can be any of
`IndexedTable`, `NDSparse`, `StructArrays.StructVector`, or `AbstractVector`

# Examples

    using IndexedTables: excludecols

    t = table([2,1],[1,3],[4,5], names=[:x,:y,:z], pkey=(1,2))

    excludecols(t, (:x,))
    excludecols(t, (2,))
    excludecols(t, pkeynames(t))
    excludecols([1,2,3], (1,))
"""
excludecols(t, cols) = excludecols(t, (cols,))
excludecols(t, cols::SpecialSelector) = excludecols(t, lowerselection(t, cols))
function excludecols(t, cols::Tuple)
    Tuple(setdiff(1:length(colnames(t)), map(x -> colindex(t, x), cols)))
end

"""
    convert(IndexedTable, pkeys, vals; kwargs...)

Construct a table with `pkeys` as primary keys and `vals` as corresponding non-indexed items.
keyword arguments will be forwarded to [`table`](@ref) constructor.

# Example
    convert(IndexedTable, Columns(x=[1,2],y=[3,4]), Columns(z=[1,2]), presorted=true)
"""
function convert(::Type{IndexedTable}, key, val; pkey = (), kwargs...)
    cs = concat_cols(key, val)
    table(cs; pkey=Tuple(1:ncols(key)), kwargs...)
end

convert(T::Type{IndexedTable}, c::Columns{<:Pair}; kwargs...) = convert(T, c.first, c.second; kwargs...)
# showing

global show_compact_when_wide = true
function set_show_compact!(flag=true)
    global show_compact_when_wide
    show_compact_when_wide = flag
end

function showtable(io::IO, t; header=nothing, cnames=colnames(t), divider=nothing, cstyle=[], full=false, ellipsis=:middle, compact=show_compact_when_wide)
    cnames = collect(cnames)
    height, width = displaysize(io)
    showrows = height-5 - (header !== nothing)
    n = length(t)
    header !== nothing && println(io, header)
    if full
        rows = [1:n;]
        showrows = n
    else
        if ellipsis == :middle
            lastfew = div(showrows, 2) - 1
            firstfew = showrows - lastfew - 1
            rows = n > showrows ? [1:firstfew; (n-lastfew+1):n] : [1:n;]
        elseif ellipsis == :end
            lst = n == showrows ?
                showrows : showrows-1 # make space for ellipse
            rows = [1:min(length(t), showrows);]
        else
            error("ellipsis must be either :middle or :end")
        end
    end
    nc = length(columns(t))
    reprs  = [ sprint(io->show(IOContext(io, :compact => true, :limit => true), columns(t)[j][i])) for i in rows, j in 1:nc ]
    strcnames = map(string, cnames)
    widths  = [ max(textwidth(get(strcnames, c, "")), isempty(reprs) ? 0 : maximum(map(textwidth, reprs[:,c]))) for c in 1:nc ]
    if compact && !isempty(widths) && sum(widths) + 2*nc > width
        return showmeta(io, t, cnames)
    end
    for c in 1:nc
        nm = get(strcnames, c, "")
        style = get(cstyle, c, nothing)
        txt = c==nc && divider!=nc ? nm : rpad(nm, widths[c]+(c==divider ? 1 : 2), " ")
        if style == nothing
            print(io, txt)
        else
            Base.with_output_color(print, style, io, txt)
        end
        if c == divider
            print(io, "│")
            length(cnames) > divider && print(io, " ")
        end
    end
    println(io)
    if divider !== nothing
        print(io, "─"^(sum(widths[1:divider])+2*divider-1), "┼")
        if !isempty(widths[divider+1:end])
            print(io, "─"^(sum(widths[divider+1:end])+2*(nc-divider)-1))
        end
    else
        if !isempty(widths)
            print(io, "─"^(sum(widths)+2*nc-2))
        end
    end
    for r in 1:size(reprs,1)
        println(io)
        for c in 1:nc
            print(io, c==nc && nc!=divider ? reprs[r,c] : rpad(reprs[r,c], widths[c]+(c==divider ? 1 : 2), " "))
            if c == divider
                print(io, "│ ")
            end
        end
        if n > showrows && ((ellipsis == :middle && r == firstfew) || (ellipsis == :end && r == size(reprs, 1)))
            if divider === nothing
                println(io)
                print(io, "⋮")
            else
                println(io)
                print(io, " "^(sum(widths[1:divider]) + 2*divider-1), "⋮")
            end
        end
    end
end

function showmeta(io, t, cnames)
    nc = length(columns(t))
    println(io, "Columns:")
    metat = Columns(([1:nc;], [Text(string(get(cnames, i, "<noname>"))) for i in 1:nc],
                     [map(eltype, columns(t))...]))
    showtable(io, metat, cnames=["#", "colname", "type"], cstyle=fill(:bold, nc), full=true, compact=false)
end

function subscriptprint(x::Integer)
    s = string(x)
    cs = Char[]
    lookup = ["₀₁₂₃₄₅₆₇₈₉"...]
    join([lookup[parse(Int, c)+1] for c in s],"")
end

function show(io::IO, t::IndexedTable{T}) where {T}
    header = "Table with $(length(t)) rows, $(length(columns(t))) columns:"
    cstyle = Dict([i=>:bold for i in t.pkey])
    cnames = string.(colnames(t))
    for (i, k) in enumerate(t.pkey)
        cstyle[k] = :bold
        #cnames[k] = cnames[k] * "$(subscriptprint(i))"
    end
    showtable(io, t, header=header, cnames=cnames, cstyle=cstyle)
end


text:
#=
---------------------------------------------------------------------
This function computes the normal vectors and Jacobian of each FACE.

The strategy follows Equation 12.33 in F.X. Giraldo's Introduction to Element-based
Galerkin Methods using Tensor-Product Bases: Analysis, Algorithms, and Applications.

Written by F.X. Giraldo on 5/2019
           Department of Applied Mathematics
           Naval Postgraduate School
           Monterey; CA 93943-5216
---------------------------------------------------------------------
=#

include("map_deriv.jl")

function compute_normals(face,intma,coord,Nface,Np,Nq,ψ,dψ)

    #global arrays
    normals=zeros(DFloat,2,Nq,Nface)
    jac_face=zeros(DFloat,Nq,Nface)

    #local arrays
    x=zeros(DFloat,Np,Np)
    y=zeros(DFloat,Np,Np)

    #loop thru the sides
    for f=1:Nface

        #Store Left & Right Elements
        ilocl=face[1,f]
        ilocr=face[2,f]
        iel=face[3,f]
        ier=face[4,f]

        #Store Element Variables
        for j=1:Np, i=1:Np
            I=intma[i,j,iel]
            x[i,j]=coord[1,I]
            y[i,j]=coord[2,I]
        end #i, j

        #Construct Mapping Derivatives: dx/dksi; dx/deta; dy/dksi
        #dy/deta
        (x_ξ,x_η)=map_deriv(ψ,dψ,x,Np,Nq,DFloat)
        (y_ξ,y_η)=map_deriv(ψ,dψ,y,Np,Nq,DFloat)

        #Compute Normals
        for l=1:Nq
            if (ilocl == 1)
                #Side 1: eta=-1
                i=l
                j=1
                normals[1,l,f]=+y_ξ[i,j]
                normals[2,l,f]=-x_ξ[i,j]
                jac_face[l,f]=sqrt( x_ξ[i,j]^2 + y_ξ[i,j]^2 )
            elseif (ilocl == 2)
                #Side 2: ksi=+1
                i=Nq
                j=l
                normals[1,l,f]=+y_η[i,j]
                normals[2,l,f]=-x_η[i,j]
                jac_face[l,f]=sqrt( x_η[i,j]^2 + y_η[i,j]^2 )
            elseif (ilocl == 3)
                #Side 3: eta=+1
                i=Nq+1-l
                j=Nq
                normals[1,l,f]=-y_ξ[i,j]
                normals[2,l,f]=+x_ξ[i,j]
                jac_face[l,f]=sqrt( x_ξ[i,j]^2 + y_ξ[i,j]^2 )
            elseif (ilocl == 4)
                #Side 4: ksi=-1
                i=1
                j=Nq+1-l
                normals[1,l,f]=-y_η[i,j]
                normals[2,l,f]=+x_η[i,j]
                jac_face[l,f]=sqrt( x_η[i,j]^2 + y_η[i,j]^2 )
            end
        end #l

        #Normalize Normals
        for l=1:Nq
            rnx=sqrt( normals[1,l,f]^2 + normals[2,l,f]^2 )
            normals[1,l,f]=normals[1,l,f]/rnx
            normals[2,l,f]=normals[2,l,f]/rnx
        end #l

    end #s

    return(normals,jac_face)
end


text:
using DataFrames: DataFrame
import CSV

mutable struct Provide{V} <: State{V}
    value::V
end

Provide(; index, init, step, autounit, _value, _type, _...) = begin
    i = value(index)
    t = value(init)
    Δt = value(step)
    df = DataFrame(_value isa String ? CSV.File(_value) : _value)
    df = autounit ? unitfy(df) : df
    v = filter(r -> r[i] >= t && iszero(typeof(Δt)(r[i] - t) % Δt), df)
    v[1, i] != t && error("incompatible index for initial time = $t\n$v")
    !all(isequal(Δt), diff(v[!, i])) && error("incompatible index for time step = $Δt\n$v")
    V = _type
    Provide{V}(v)
end

constructortags(::Val{:Provide}) = (:index, :init, :step, :autounit)

updatetags!(d, ::Val{:Provide}; _...) = begin
    !haskey(d, :index) && (d[:index] = QuoteNode(:index))
    !haskey(d, :init) && (d[:init] = :(context.clock.time))
    !haskey(d, :step) && (d[:step] = :(context.clock.step))
    !haskey(d, :autounit) && (d[:autounit] = true)
end

genvartype(v::VarInfo, ::Val{:Provide}; V, _...) = @q Provide{$V}

gendefault(v::VarInfo, ::Val{:Provide}) = istag(v, :parameter) ? genparameter(v) : genbody(v)

genupdate(v::VarInfo, ::Val{:Provide}, ::MainStep; kw...) = nothing


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module Cloud

using PyCall
using Conda
using Pkg
using Random

"""
    cache_run_config(subscription_id::String, resource_group::String, workspace_name::String, compute_name::String[, path::String])

Writes the AML config information to a file at `path`. If `path` isn't set then the config will be used globally for that OpticSim install.
"""
function cache_run_config(subscription_id::String, resource_group::String, workspace_name::String, compute_name::String, path::String = joinpath(@__DIR__, "amlconf"))
    open(path, "w") do io
        write(io, subscription_id * "\n")
        write(io, resource_group * "\n")
        write(io, workspace_name * "\n")
        write(io, compute_name)
   end
   nothing
end

"""
    get_cached_run_config([path::String])

Reads the AML config information from a file at `path`. If not specified then the global config will be read.
"""
function get_cached_run_config(path::String = joinpath(@__DIR__, "amlconf"))
    open(path, "r") do io
        subscription_id = readline(io)
        resource_group = readline(io)
        workspace_name = readline(io)
        compute_name = readline(io)
        return subscription_id, resource_group, workspace_name, compute_name
    end
end

"""
    submit_run_to_AML(run_name::String, path_to_script::String, script_args::Vector{String} = nothing, sampled_args:Dict{String,Vector{String}} = nothing, config_path::String; hyperdrive_concurrent_runs::Int = 10)
    submit_run_to_AML(run_name::String, path_to_script::String, subscription_id::String, resource_group::String, workspace_name::String, compute_name::String, script_args::Vector{String} = nothing, sampled_args::Dict{String, Vector{String}} = nothing; hyperdrive_concurrent_runs::Int = 10)

Submit a run to AML, `path_to_script` is relative to your local package root (i.e. location of `Project.toml`).
`script_args` are a series of arguments to your script as strings.
`sampled_args` is a dictionary where keys are argument names and values are lists of values (as strings) that that argument will take.
`config_path` is a path to a config file as written by [`cache_run_config`](@ref), if not specified the global config is used. Alternatively this information can be provided directly using the second method above.
`hyperdrive_concurrent_runs` is the maximum number of concurrent runs that will execute on AML (limited by your compute cluster size).
"""
function submit_run_to_AML(run_name::String, path_to_script::String, script_args::Union{Nothing,Vector{String}} = nothing,
                           sampled_args::Union{Nothing,Dict{String,Vector{String}}} = nothing;
                           config_path::String = joinpath(@__DIR__, "amlconf"),
                           hyperdrive_concurrent_runs::Int = 10)
    subscription_id, resource_group, workspace_name, compute_name = get_cached_run_config(config_path)
    submit_run_to_AML(run_name, path_to_script, subscription_id, resource_group, workspace_name, compute_name,
                      script_args, sampled_args, hyperdrive_concurrent_runs=hyperdrive_concurrent_runs)
end

function submit_run_to_AML(run_name::String, path_to_script::String,
                           subscription_id::String, resource_group::String, workspace_name::String, compute_name::String,
                           script_args::Union{Nothing,Vector{String}} = nothing,
                           sampled_args::Union{Nothing,Dict{String,Vector{String}}} = nothing;
                           hyperdrive_concurrent_runs::Int = 10)

    dockerfile = open(joinpath(@__DIR__, "dockerfile")) do file
        read(file, String)
    end

    # add deps to dockerfile
    project_dict = Pkg.TOML.parsefile(Base.active_project())
    packages = collect(keys(project_dict["deps"]))
    sort!(packages)
    pkg_install_cmd = "RUN julia -e \"using Pkg; "
    for package_name in packages
        if "compat" in keys(project_dict) && package_name in keys(project_dict["compat"])
            pkg_install_cmd = pkg_install_cmd * "Pkg.add(name=\\\"" * package_name * "\\\", version=\\\"" * project_dict["compat"][package_name] * "\\\");"
        else
            pkg_install_cmd = pkg_install_cmd * "Pkg.add(\\\"" * package_name * "\\\");"
        end
    end
    if "OpticSim" in packages
        # build OpticSim if it is there
        pkg_install_cmd = pkg_install_cmd * "Pkg.build(\\\"OpticSim\\\");"
    end
    pkg_install_cmd = pkg_install_cmd * "\""

    dockerfile = dockerfile * pkg_install_cmd

    # TODO maybe compile sysimage in docker - would be horribly slow but should speed up import a lot?

    source_directory = joinpath(dirname(Base.active_project()))

    if isfile(joinpath(source_directory, "Manifest.toml")) && !isfile(joinpath(source_directory, ".amlignore"))
        println("No .amlignore file found, creating one")
        open(joinpath(source_directory, ".amlignore"), "w") do io
            write(io, "Manifest.toml\n")
        end
    end

    # set up env for python stuff
    try
        pyimport("azureml.core")
    catch
        # FIXME maybe won't work, might need to restart Julia after this?
        Conda.add("python=3.7")
        Conda.add("pip=20.1.1")
        Conda.pip_interop(true)
        Conda.pip("install", "azureml-sdk")
        Pkg.build("PyCall")
    end

    # copy entry_script from here to source_directory
    entry_script_path = "entry_script_" * randstring() * ".py"
    cp(joinpath(@__DIR__, "entry_script.py"), joinpath(source_directory, entry_script_path))

    py"""
    import os
    import webbrowser
    from azureml.core import Environment, Experiment, Run, Workspace, ScriptRunConfig
    import azureml.train.hyperdrive as hyperdrive
    from azureml.train.hyperdrive.parameter_expressions import choice

    def get_hyperparam_dict(param_dict):
        hyper_param_dict = {}
        num_params = 1
        for key, value in param_dict.items():
            hyper_param_dict[key] = choice(value)
            num_params = num_params * len(value)

        return hyper_param_dict, num_params

    def submit_run(subscription_id, resource_group, workspace_name, compute_name, source_directory, julia_script,
                   script_args, sampled_args, run_name, dockerfile, entry_script_path, hyperdrive_concurrent_runs):
        workspace = Workspace(subscription_id, resource_group, workspace_name)

        compute_target = workspace.compute_targets[compute_name]

        env = Environment("opticsim")
        env.docker.base_image = None
        env.docker.base_dockerfile = dockerfile

        args = [julia_script.replace(os.sep, "/")]
        if script_args is not None:
            args.extend(script_args)

        src = ScriptRunConfig(source_directory=source_directory,
                              script=entry_script_path,
                              arguments=args,
                              compute_target=compute_target,
                              environment=env)
        src.run_config.docker.use_docker = True

        exp_name = os.getlogin() + "-opticsim"
        experiment = Experiment(workspace, exp_name)

        if sampled_args is not None:
            sampling_params, num_params = get_hyperparam_dict(sampled_args)
            param_sampling = hyperdrive.GridParameterSampling(sampling_params)
            hyperdrive_run_config = hyperdrive.HyperDriveConfig(run_config=src,
                                                                hyperparameter_sampling=param_sampling,
                                                                max_concurrent_runs=hyperdrive_concurrent_runs,
                                                                primary_metric_name="",
                                                                primary_metric_goal=hyperdrive.PrimaryMetricGoal.MINIMIZE,
                                                                max_total_runs=num_params)
            run_object = experiment.submit(hyperdrive_run_config, tags={"run_name": run_name})
        else:
            run_object = experiment.submit(src, tags={"run_name": run_name})

        webbrowser.open_new(run_object.get_portal_url())
    """

    py"submit_run"(subscription_id, resource_group, workspace_name, compute_name, source_directory, path_to_script,
                   script_args, sampled_args, run_name, dockerfile, entry_script_path, hyperdrive_concurrent_runs)

    # remove the entry script
    rm(joinpath(source_directory, entry_script_path))
end

export submit_run_to_AML, cache_run_config, get_cached_run_config

end