text:
@testset "Utilities" begin
  @test_throws ArgumentError readgeotable("doesnotexist.csv")
  @test_throws AssertionError readgeotable(joinpath(datadir,"data3D.tsv"), delim='\t', coordnames=[:a])
end


text:
"""
    dense_bottleneck(inplanes, growth_rate)

Create a Densenet bottleneck layer
([reference](https://arxiv.org/abs/1608.06993)).

# Arguments
- `inplanes`: number of input feature maps
- `outplanes`: number of output feature maps on bottleneck branch
               (and scaling factor for inner feature maps; see ref)
"""
function dense_bottleneck(inplanes, outplanes)
  inner_channels = 4 * outplanes
  m = Chain(conv_bn((1, 1), inplanes, inner_channels; bias = false, rev = true)...,
            conv_bn((3, 3), inner_channels, outplanes; pad = 1, bias = false, rev = true)...)

  SkipConnection(m, (mx, x) -> cat(x, mx; dims = 3))
end

"""
    transition(inplanes, outplanes)

Create a DenseNet transition sequence
([reference](https://arxiv.org/abs/1608.06993)).

# Arguments
- `inplanes`: number of input feature maps
- `outplanes`: number of output feature maps
"""
transition(inplanes, outplanes) =
  [conv_bn((1, 1), inplanes, outplanes; bias = false, rev = true)...,
   MeanPool((2, 2))]

"""
    dense_block(inplanes, growth_rates)

Create a sequence of DenseNet bottlenecks increasing
the number of output feature maps by `growth_rates` with each block
([reference](https://arxiv.org/abs/1608.06993)).

# Arguments
- `inplanes`: number of input feature maps to the full sequence
- `growth_rates`: the growth (additive) rates of output feature maps
                  after each block (a vector of `k`s from the ref)
"""
dense_block(inplanes, growth_rates) = [dense_bottleneck(i, o)
  for (i, o) in zip(inplanes .+ cumsum([0, growth_rates[1:(end - 1)]...]), growth_rates)]

"""
    densenet(inplanes, growth_rates; reduction = 0.5, nclasses = 1000)

Create a DenseNet model
([reference](https://arxiv.org/abs/1608.06993)).

# Arguments
- `inplanes`: the number of input feature maps to the first dense block
- `growth_rates`: the growth rates of output feature maps within each
                  [`dense_block`](#) (a vector of vectors)
- `reduction`: the factor by which the number of feature maps is scaled across each transition
- `nclasses`: the number of output classes
"""
function densenet(inplanes, growth_rates; reduction = 0.5, nclasses = 1000)
  layers = conv_bn((7, 7), 3, inplanes; stride = 2, pad = (3, 3), bias = false)
  push!(layers, MaxPool((3, 3), stride = 2, pad = (1, 1)))

  outplanes = 0
  for (i, rates) in enumerate(growth_rates)
    outplanes = inplanes + sum(rates)
    append!(layers, dense_block(inplanes, rates))
    (i != length(growth_rates)) &&
      append!(layers, transition(outplanes, floor(Int, outplanes * reduction)))
    inplanes = floor(Int, outplanes * reduction)
  end
  push!(layers, BatchNorm(outplanes, relu))

  return Chain(Chain(layers...),
               Chain(AdaptiveMeanPool((1, 1)),
                     MLUtils.flatten,
                     Dense(outplanes, nclasses)))
end

"""
    densenet(nblocks; growth_rate = 32, reduction = 0.5, nclasses = 1000)

Create a DenseNet model
([reference](https://arxiv.org/abs/1608.06993)).

# Arguments
- `nblocks`: number of dense blocks between transitions
- `growth_rate`: the output feature map growth rate of dense blocks (i.e. `k` in the ref)
- `reduction`: the factor by which the number of feature maps is scaled across each transition
- `nclasses`: the number of output classes
"""
densenet(nblocks; growth_rate = 32, reduction = 0.5, nclasses = 1000) =
  densenet(2 * growth_rate, [fill(growth_rate, n) for n in nblocks];
           reduction = reduction, nclasses = nclasses)

"""
    DenseNet(nblocks::NTuple{N, <:Integer};
             growth_rate = 32, reduction = 0.5, nclasses = 1000)

Create a DenseNet model
([reference](https://arxiv.org/abs/1608.06993)).
See also [`densenet`](#).

# Arguments
- `nblocks`: number of dense blocks between transitions
- `growth_rate`: the output feature map growth rate of dense blocks (i.e. `k` in the paper)
- `reduction`: the factor by which the number of feature maps is scaled across each transition
- `nclasses`: the number of output classes
"""
struct DenseNet
  layers
end

function DenseNet(nblocks::NTuple{N, <:Integer};
                  growth_rate = 32, reduction = 0.5, nclasses = 1000) where N
  layers = densenet(nblocks; growth_rate = growth_rate,
                             reduction = reduction,
                             nclasses = nclasses)

  DenseNet(layers)
end

@functor DenseNet

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

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

const densenet_config = Dict(121 => (6, 12, 24, 16),
                             161 => (6, 12, 36, 24),
                             169 => (6, 12, 32, 32),
                             201 => (6, 12, 48, 32))

"""
    DenseNet(config::Int = 121; pretrain = false)
    DenseNet(transition_config::NTuple{N,Int})

Create a DenseNet model with specified configuration. Currently supported values are (121, 161, 169, 201)
([reference](https://arxiv.org/abs/1608.06993)).
Set `pretrain = true` to load the model with pre-trained weights for ImageNet.

!!! warning
    `DenseNet` does not currently support pretrained weights.

See also [`Metalhead.densenet`](#).
"""
function DenseNet(config::Int = 121; pretrain = false)
  @assert config in keys(densenet_config) "`config` must be one out of $(sort(collect(keys(densenet_config))))."
  model = DenseNet(densenet_config[config])

  pretrain && loadpretrain!(model, string("DenseNet", config))
  return model
end

# deprecations
@deprecate DenseNet121(; kw...) DenseNet(121; kw...)
@deprecate DenseNet161(; kw...) DenseNet(161; kw...)
@deprecate DenseNet169(; kw...) DenseNet(169; kw...)
@deprecate DenseNet201(; kw...) DenseNet(201; kw...)

ext: jl


text:
using .BaseAgent
using Params
using Traffic.Simple
using Channel.Gilbert
using Util
using Params.QParams

export CooperativeQ

type CooperativeQ <: Agent
    s :: AgentState
    Q :: Array{Float64, 3}
    visit :: Array{Float64, 3}
    a :: Int
    state :: (Int, Int)
    P_tx :: Float64
    bitrate :: Float64
    status :: Status
    expertness :: Float64
    buf_interval :: Int
    beta_idle :: Float64
end

const n_p_levels = length(Params.P_levels)
const idle_action = Params.n_channel * n_p_levels + 1

function CooperativeQ(i, P, pos)
    Q = rand(int(Params.n_channel), P.buf_levels + 1, idle_action)
    Q *= Params.P_tx * Params.t_slot # a good initial randomization
    visit = zeros(Params.n_channel, P.buf_levels + 1, idle_action)
    CooperativeQ(AgentState(i, P, pos), Q, visit, 0, (0, 0), 0, 0, Initialized, 0, div(Params.B + 1, P.buf_levels), P.beta_idle)
end

function policy!(a :: CooperativeQ)
    # check whether we have packets first
    if rand() < epsilon
        a.a = rand(1:idle_action)
    else
        a.a = indmax(a.Q[a.s.chan, div(a.s.B_empty, a.buf_interval) + 1, :])
    end
    nothing
end

function BaseAgent.act(a :: CooperativeQ, env, t)
    if a.status == Initialized
        if a.s.B_max - a.s.B_empty == 0
            a.a = -1
            return idle(a)
        end

        policy!(a)
        a.state = (a.s.chan, div(a.s.B_empty, a.buf_interval) + 1)
        a.visit[a.s.chan, a.state[2], a.a] += 1

        if a.a == idle_action
            return idle(a)
        end

        chan = int8(fld(a.a - 1, n_p_levels) + 1)
        return switch(a, chan)
    elseif a.status == Switched
        return sense(a, env, detect_traffic)
    elseif a.status == Sensed
        if detect_traffic(a.s, env.traffics[a.s.chan], a.s.t_remaining)
            return idle(a)
        else
            a.P_tx = Params.P_levels[(a.a - 1) % n_p_levels + 1]
            a.bitrate = capacity(env.channels[a.s.chan], a.P_tx, a.s.pos)
            pkgs_to_send = min(a.s.B_max - a.s.B_empty, floor(a.s.t_remaining * a.bitrate / Params.pkt_size))

            if pkgs_to_send == 0
                return idle(a)
            end
            return transmit!(a, a.P_tx, env, pkgs_to_send)
        end
    elseif a.status == Transmitted
        return idle(a)
    end
end

function alpha(a :: CooperativeQ)
    return 0.2 + 0.8 / (1 + a.visit[a.state[1], a.state[2], a.a])
end

function BaseAgent.feedback(a :: CooperativeQ, res :: Result, idle :: Bool = false, n_pkt :: Int = 0)
    feedback(a.s, res, idle, n_pkt)
    # if we didn't take any action that's worth learning, return immediately
    if a.a == -1
        return nothing
    end
    r :: Float64 = 0
    if idle
        r = - a.beta_idle * a.bitrate * Params.t_slot / a.s.E_slot
    elseif res == Success
        K = 1 # a.P_tx ^ 2 * Params.t_slot / a.bitrate
        r = K * Params.pkt_size * n_pkt / a.s.E_slot
    elseif res == Collision
        r = - beta_md * a.bitrate * Params.t_slot / a.s.E_slot
    elseif res == LostInChannel
        r = - beta_loss * a.bitrate * Params.t_slot / a.s.E_slot
    else
        return nothing
    end
    # UPDATE Q
    max_Q_next = maximum(a.Q[a.s.chan, floor(a.s.B_empty / a.buf_interval) + 1])
    Q_now = a.Q[a.state[1], a.state[2], a.a]
    a.Q[a.state[1], a.state[2], a.a] += alpha(a) * (r + discount * max_Q_next - Q_now)
    # Update expertness
    if r > 0
        a.expertness += r # sum of reinforcements for now
    end
    nothing
end

# Coordinator object for Cooperative Q
type Coordinator
    US :: Array{Float64, 3} # U*S where U*S*V' = svd(Q) for each agent
    Vt :: Array{Float64, 3} # V' where U*S*V' = svd(Q) for each agent
    function Coordinator(P :: ParamT)
        c = new()
        c.US = zeros(P.n_agent, int(Params.n_channel) * (P.buf_levels + 1), P.d_svd)
        c.Vt = zeros(P.n_agent, P.d_svd, Params.n_channel * length(Params.P_levels) + 1)
    end
end

function BaseAgent.initcoordinator(:: Type{CooperativeQ}, P :: ParamT)
    return Coordinator(P)
end

function BaseAgent.cooperate(agents :: Vector{CooperativeQ}, P :: ParamT, coordinator :: Coordinator, t)
    const n_agent = int(P.n_agent)
    const sharingperiod = P.sharingperiod
    # size of US + size of Vt, used for data sharing time and energy computation
    # const sizeQ = P.d_svd * 64 * (n_channel * (P.buf_levels + 1) + (n_channel * length(P_levels) + 1))
    const sizeQ = 64 * n_channel * (P.buf_levels + 1) * (n_channel * length(P_levels) + 1)
    # time slots required for an agent to send/receive Q matrix by sending/receiving US & Vt
    const timeQ = int(ceil(sizeQ ./ Params.controlcapacity ./ t_slot))
    # time required for an agent to send/receive Q matrix by sending/receiving US & Vt
    const rawtimeQ0 = sizeQ ./ Params.controlcapacity ./ t_slot
    # Shape of the Q matrix
    const shapeQ = (int(Params.n_channel) * (P.buf_levels + 1), Params.n_channel * length(Params.P_levels) + 1)
    Q = zeros(n_agent, int(Params.n_channel), P.buf_levels + 1, Params.n_channel * length(Params.P_levels) + 1)
    if n_agent > 1
        US = coordinator.US
        Vt = coordinator.Vt
        # join Cooperative Q learners
        if t > Params.t_saturation && t % sharingperiod < 2 * n_agent * timeQ
            tt = t % sharingperiod
            i = fld(tt, timeQ) + 1
            # sending phase (from CR)
            if tt < n_agent * timeQ
                # receiving complete, update combined Q matrix
                # using only positive expertness ones
                if tt - (i-1) * timeQ == timeQ - 1
                    #= u, s, v = svd(reshape(agents[i].expertness .* agents[i].Q, shapeQ))
                    US[i, :, :] = u[:,1:P.d_svd]*diagm(s[1:P.d_svd])
                    Vt[i, :, :] = v'[1:P.d_svd, :] =#
                    Q[i,:,:,:] = agents[i].expertness .* agents[i].Q
                    expertness[i] = agents[i].expertness
                end
                if tt == (i - 1) * timeQ
                    # initialize control channel usage overhead
                    rawtimeQ = rawtimeQ0
                end
                # Calculate control channel energy overhead
                if rawtimeQ > 1
                    agents[i].s.E_slot += Params.P_tx * t_slot # add control channel energy overhead
                else
                    agents[i].s.E_slot += Params.P_tx * t_slot * rawtimeQ # add control channel energy overhead
                end
                rawtimeQ -= 1 # doesn't matter when rawtimeQ < 1, it'll be reset anyways
            else # receiving phase (to CR)
                i -= n_agent
                # sending complete, update agent's Q matrix
                if tt - (i+n_agent-1) * timeQ == timeQ - 1
                    #fill!(agents[i].Q, 0)
                    # use up-to-date Q function on CR's side
                    weight = 1 - P.trustQ
                    agents[i].Q *= weight #* reshape(Q[j,:,:,:], size(agents[i].Q))
                    # learning from experts (LE)
                    # choose some experts randomly
                    experts = zeros(Bool, n_agent)
                    for expert in shuffle!([[x for x=1:i-1],[x for x=i+1:n_agent]])[1:fld(n_agent, 2)]
                        experts[expert] = true
                    end
                    for j=1:n_agent
                        if i == j
                            continue
                        elseif !experts[j]
                            continue
                        else
                            if expertness[j] > expertness[i]
                                weight = P.trustQ * (expertness[j] - expertness[i])
                                # normalize
                                weight ./= sum((expertness - expertness[i]) .* (expertness .> expertness[i]))
                                #= us = slice(US, (j, 1:size(US)[2], 1:size(US)[3]))
                                vt = slice(Vt, (j, 1:size(Vt)[2], 1:size(Vt)[3]))
                                agents[i].Q += weight * reshape(us * vt, size(agents[i].Q)) =#
                                agents[i].Q += weight * reshape(Q[j,:,:,:], size(agents[i].Q))
                            end
                        end
                    end
                    agents[i].expertness = 0 # *= 1 - P.trustQ
                end
            end
        end
    end
end


text:
# This is experimental, as few arches support it, and I can't think of many uses other than floating point RNGs.

@inline __ifmalo(v1, v2, v3) = ((((v1 % UInt64)) * ((v2 % UInt64))) & 0x000fffffffffffff) + (v3 % UInt64)
@inline _ifmalo(v1, v2, v3) = __ifmalo(v1, v2, v3)
function ifmahi_quote(W)
    mask = W > 1 ? llvmconst(W, "i64 4503599627370495") : "4503599627370495"
    shift = W > 1 ? llvmconst(W, "i128 52") : "52"
    t64 = W > 1 ? "<$W x i64>" : "i64"
    t128 = W > 1 ? "<$W x i128>" : "i128"
    instrs = """
        %a52 = and $t64 %0, $mask
        %b52 = and $t64 %1, $mask
        %a128 = zext $t64 %a52 to $t128
        %b128 = zext $t64 %b52 to $t128
        %c128 = mul $t128 %a128, %b128
        %csr = lshr $t128 %c128, $shift
        %c64 = trunc $t128 %csr to $t64
        %res = add $t64 %c64, %2
        ret $t64 %res
    """
    jt = W > 1 ? :(_Vec{$W,UInt64}) : :UInt64
    call = :($LLVMCALL($instrs, $jt, Tuple{$jt,$jt,$jt}, data(v1), data(v2), data(v3)))
    W > 1 && (call = Expr(:call, :Vec, call))
    Expr(:block, Expr(:meta,:inline), call)
end
@generated _ifmahi(v1::UInt64, v2::UInt64, v3::UInt64) = ifmahi_quote(1)
@generated __ifmahi(v1::Vec{W,UInt64}, v2::Vec{W,UInt64}, v3::Vec{W,UInt64}) where {W} = ifmahi_quote(W)

function ifmaquote(W::Int, lo::Bool)
    op = lo ? "@llvm.x86.avx512.vpmadd52l.uq.$(64W)" : "@llvm.x86.avx512.vpmadd52h.uq.$(64W)"
    decl = "declare <$W x i64> $op(<$W x i64>, <$W x i64>, <$W x i64>)"
    instrs = "%res = call <$W x i64> $op(<$W x i64> %0, <$W x i64> %1, <$W x i64> %2)\n ret <$W x i64> %res"
    llvmcall_expr(decl, instrs, :(_Vec{$W,UInt64}), :(Tuple{_Vec{$W,UInt64},_Vec{$W,UInt64},_Vec{$W,UInt64}}), "<$W x i64>", ["<$W x i64>","<$W x i64>","<$W x i64>"], [:(data(v3)), :(data(v1)), :(data(v2))])
end
@generated _ifmalo(v1::Vec{W,UInt64},v2::Vec{W,UInt64},v3::Vec{W,UInt64}, ::True) where {W} = ifmaquote(W, true)
@generated _ifmahi(v1::Vec{W,UInt64},v2::Vec{W,UInt64},v3::Vec{W,UInt64}, ::True) where {W} = ifmaquote(W, false)

@inline _ifmalo(v1::Vec{W,UInt64},v2::Vec{W,UInt64},v3::Vec{W,UInt64}, ::False) where {W} = __ifmalo(v1, v2, v3)
@inline _ifmahi(v1::Vec{W,UInt64},v2::Vec{W,UInt64},v3::Vec{W,UInt64}, ::False) where {W} = __ifmahi(v1, v2, v3)

@inline function _ifmalo(v1::Vec{W,UInt64}, v2::Vec{W,UInt64}, v3::Vec{W,UInt64}) where {W}
    use_ifma = (has_feature(Val(:x86_64_avx512ifma)) & _ispow2(StaticInt{W}())) &
        (gt(StaticInt{W}(), StaticInt{8}()) & le(StaticInt{W}() * StaticInt{8}(), register_size()))
    _ifmalo(v1, v2, v3, use_ifma)
end
@inline function _ifmahi(v1::Vec{W,UInt64}, v2::Vec{W,UInt64}, v3::Vec{W,UInt64}) where {W}
    use_ifma = (has_feature(Val(:x86_64_avx512ifma)) & _ispow2(StaticInt{W}())) &
        (gt(StaticInt{W}(), StaticInt{8}()) & le(StaticInt{W}() * StaticInt{8}(), register_size()))
    _ifmahi(v1, v2, v3, use_ifma)
end

"""
    ifmalo(v1, v2, v3)

Multiply unsigned integers `v1` and `v2`, adding the lower 52 bits to `v3`.

Requires `has_feature(Val(:x86_64_avx512ifma))` to be fast.
"""
@inline ifmalo(v1, v2, v3) = _ifmalo(v1 % UInt64, v2 % UInt64, v3 % UInt64)
"""
    ifmalo(v1, v2, v3)

Multiply unsigned integers `v1` and `v2`, adding the upper 52 bits to `v3`.

Requires `has_feature(Val(:x86_64_avx512ifma))` to be fast.
"""
@inline ifmahi(v1, v2, v3) = ((a,b,c) = promote(v1 % UInt64, v2 % UInt64, v3 % UInt64); _ifmahi(a, b, c))

@inline function _vfmadd_fast_uint64(a::Vec{W,UInt64}, b::Vec{W,UInt64}, c::Vec{W,UInt64}, ::True) where {W}
    ifmalo(a, b, c)
end
@inline function _vfmadd_fast_uint64(a::Vec{W,UInt64}, b::Vec{W,UInt64}, c::Vec{W,UInt64}, ::False) where {W}
    add_fast(mul_fast(a, b), c)
end

@inline function vfmadd_fast(a::Vec{W,UInt64},b::Vec{W,UInt64},c::Vec{W,UInt64}) where {W}
    _vfmadd_fast_uint64(a, b, c, has_feature(Val(:x86_64_avx512ifma)))
end


text:
using DiffEqFlux, Flux
using LinearAlgebra, Distributions
using GalacticOptim
using Test
using DataInterpolations

function run_test(f, layer, atol)

    data_train_vals = rand(500)
    data_train_fn = f.(data_train_vals)

    function loss_function(θ)
        data_pred = [layer(x, θ) for x in data_train_vals]
        loss = sum(abs.(data_pred.-data_train_fn))/length(data_train_fn)
        return loss
    end

    function cb(p,l)
        @show l
        return false
    end

    optfunc = GalacticOptim.OptimizationFunction((x, p) -> loss_function(x), GalacticOptim.AutoZygote())
    optprob = GalacticOptim.OptimizationProblem(optfunc, layer.saved_points)
    res = GalacticOptim.solve(optprob, ADAM(0.1), cb=cb, maxiters = 100)

    optprob = GalacticOptim.OptimizationProblem(optfunc, res.minimizer)
    res = GalacticOptim.solve(optprob, ADAM(0.1), cb=cb, maxiters = 100)
    opt = res.minimizer

    data_validate_vals = rand(100)
    data_validate_fn = f.(data_validate_vals)

    data_validate_pred = [layer(x,opt) for x in data_validate_vals]

    return sum(abs.(data_validate_pred.-data_validate_fn))/length(data_validate_fn) < atol
end

##test 01: affine function, Linear Interpolation
a, b = rand(2)
f = x -> a*x + b
layer = SplineLayer((0.0,1.0),0.01,LinearInterpolation)
@test run_test(f, layer, 0.05)

##test 02: non-linear function, Quadratic Interpolation
a, b, c = rand(3)
f = x -> a*x^2+ b*x + x
layer = SplineLayer((0.0,1.0),0.01,QuadraticInterpolation)
@test run_test(f, layer, 0.05)

##test 03: non-linear function, Quadratic Spline
a, b, c = rand(3)
f = x -> a*sin(b*x+c)
layer = SplineLayer((0.0,1.0),0.1,QuadraticSpline)
@test run_test(f, layer, 0.05)

##test 04: non-linear function, Cubic Spline
f = x -> exp(x)*x^2
layer = SplineLayer((0.0,1.0),0.1,CubicSpline)
@test run_test(f, layer, 0.05)


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

#trunc
dt = Dates.Date(2012, 12, 21)
@test trunc(dt, Dates.Year) == Dates.Date(2012)
@test trunc(dt, Dates.Month) == Dates.Date(2012, 12)
@test trunc(dt, Dates.Day) == Dates.Date(2012, 12, 21)
dt = Dates.DateTime(2012, 12, 21, 16, 30, 20, 200)
@test trunc(dt, Dates.Year) == Dates.DateTime(2012)
@test trunc(dt, Dates.Month) == Dates.DateTime(2012, 12)
@test trunc(dt, Dates.Day) == Dates.DateTime(2012, 12, 21)
@test trunc(dt, Dates.Hour) == Dates.DateTime(2012, 12, 21, 16)
@test trunc(dt, Dates.Minute) == Dates.DateTime(2012, 12, 21, 16, 30)
@test trunc(dt, Dates.Second) == Dates.DateTime(2012, 12, 21, 16, 30, 20)
@test trunc(dt, Dates.Millisecond) == Dates.DateTime(2012, 12, 21, 16, 30, 20, 200)
t = Dates.Time(1, 2, 3, 4, 5, 6)
@test trunc(t, Dates.Hour) == Dates.Time(1)
@test trunc(t, Dates.Minute) == Dates.Time(1, 2)
@test trunc(t, Dates.Second) == Dates.Time(1, 2, 3)
@test trunc(t, Dates.Millisecond) == Dates.Time(1, 2, 3, 4)
@test trunc(t, Dates.Microsecond) == Dates.Time(1, 2, 3, 4, 5)
@test trunc(t, Dates.Nanosecond) == Dates.Time(1, 2, 3, 4, 5, 6)

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

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

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

@test Dates.firstdayofmonth(Jan) == Dates.DateTime(2013, 1, 1)
@test Dates.firstdayofmonth(Feb) == Dates.DateTime(2013, 2, 1)
@test Dates.firstdayofmonth(Mar) == Dates.DateTime(2013, 3, 1)
@test Dates.firstdayofmonth(Apr) == Dates.DateTime(2013, 4, 1)
@test Dates.firstdayofmonth(May) == Dates.DateTime(2013, 5, 1)
@test Dates.firstdayofmonth(Jun) == Dates.DateTime(2013, 6, 1)
@test Dates.firstdayofmonth(Jul) == Dates.DateTime(2013, 7, 1)
@test Dates.firstdayofmonth(Aug) == Dates.DateTime(2013, 8, 1)
@test Dates.firstdayofmonth(Sep) == Dates.DateTime(2013, 9, 1)
@test Dates.firstdayofmonth(Oct) == Dates.DateTime(2013, 10, 1)
@test Dates.firstdayofmonth(Nov) == Dates.DateTime(2013, 11, 1)
@test Dates.firstdayofmonth(Dec) == Dates.DateTime(2013, 12, 1)

@test Dates.firstdayofmonth(Date(Jan)) == Dates.Date(2013, 1, 1)
@test Dates.firstdayofmonth(Date(Feb)) == Dates.Date(2013, 2, 1)
@test Dates.firstdayofmonth(Date(Mar)) == Dates.Date(2013, 3, 1)
@test Dates.firstdayofmonth(Date(Apr)) == Dates.Date(2013, 4, 1)
@test Dates.firstdayofmonth(Date(May)) == Dates.Date(2013, 5, 1)
@test Dates.firstdayofmonth(Date(Jun)) == Dates.Date(2013, 6, 1)
@test Dates.firstdayofmonth(Date(Jul)) == Dates.Date(2013, 7, 1)
@test Dates.firstdayofmonth(Date(Aug)) == Dates.Date(2013, 8, 1)
@test Dates.firstdayofmonth(Date(Sep)) == Dates.Date(2013, 9, 1)
@test Dates.firstdayofmonth(Date(Oct)) == Dates.Date(2013, 10, 1)
@test Dates.firstdayofmonth(Date(Nov)) == Dates.Date(2013, 11, 1)
@test Dates.firstdayofmonth(Date(Dec)) == Dates.Date(2013, 12, 1)

# Test first day of week; 2014-01-06 is a Monday = 1st day of week
a = Dates.Date(2014, 1, 6)
b = Dates.Date(2014, 1, 7)
c = Dates.Date(2014, 1, 8)
d = Dates.Date(2014, 1, 9)
e = Dates.Date(2014, 1, 10)
f = Dates.Date(2014, 1, 11)
g = Dates.Date(2014, 1, 12)
@test Dates.firstdayofweek(a) == a
@test Dates.firstdayofweek(b) == a
@test Dates.firstdayofweek(c) == a
@test Dates.firstdayofweek(d) == a
@test Dates.firstdayofweek(e) == a
@test Dates.firstdayofweek(f) == a
@test Dates.firstdayofweek(g) == a
# Test firstdayofweek over the course of the year
dt = a
for i = 0:364
    @test Dates.firstdayofweek(dt) == a + Dates.Week(div(i, 7))
    dt += Dates.Day(1)
end
a = Dates.DateTime(2014, 1, 6)
b = Dates.DateTime(2014, 1, 7)
c = Dates.DateTime(2014, 1, 8)
d = Dates.DateTime(2014, 1, 9)
e = Dates.DateTime(2014, 1, 10)
f = Dates.DateTime(2014, 1, 11)
g = Dates.DateTime(2014, 1, 12)
@test Dates.firstdayofweek(a) == a
@test Dates.firstdayofweek(b) == a
@test Dates.firstdayofweek(c) == a
@test Dates.firstdayofweek(d) == a
@test Dates.firstdayofweek(e) == a
@test Dates.firstdayofweek(f) == a
@test Dates.firstdayofweek(g) == a
dt = a
for i = 0:364
    @test Dates.firstdayofweek(dt) == a + Dates.Week(div(i, 7))
    dt += Dates.Day(1)
end
@test Dates.firstdayofweek(Dates.DateTime(2013, 12, 24)) == Dates.DateTime(2013, 12, 23)
# Test last day of week; Sunday = last day of week
# 2014-01-12 is a Sunday
a = Dates.Date(2014, 1, 6)
b = Dates.Date(2014, 1, 7)
c = Dates.Date(2014, 1, 8)
d = Dates.Date(2014, 1, 9)
e = Dates.Date(2014, 1, 10)
f = Dates.Date(2014, 1, 11)
g = Dates.Date(2014, 1, 12)
@test Dates.lastdayofweek(a) == g
@test Dates.lastdayofweek(b) == g
@test Dates.lastdayofweek(c) == g
@test Dates.lastdayofweek(d) == g
@test Dates.lastdayofweek(e) == g
@test Dates.lastdayofweek(f) == g
@test Dates.lastdayofweek(g) == g
dt = a
for i = 0:364
    @test Dates.lastdayofweek(dt) == g + Dates.Week(div(i, 7))
    dt += Dates.Day(1)
end
a = Dates.DateTime(2014, 1, 6)
b = Dates.DateTime(2014, 1, 7)
c = Dates.DateTime(2014, 1, 8)
d = Dates.DateTime(2014, 1, 9)
e = Dates.DateTime(2014, 1, 10)
f = Dates.DateTime(2014, 1, 11)
g = Dates.DateTime(2014, 1, 12)
@test Dates.lastdayofweek(a) == g
@test Dates.lastdayofweek(b) == g
@test Dates.lastdayofweek(c) == g
@test Dates.lastdayofweek(d) == g
@test Dates.lastdayofweek(e) == g
@test Dates.lastdayofweek(f) == g
@test Dates.lastdayofweek(g) == g
dt = a
for i = 0:364
    @test Dates.lastdayofweek(dt) == g + Dates.Week(div(i, 7))
    dt += Dates.Day(1)
end
@test Dates.lastdayofweek(Dates.DateTime(2013, 12, 24)) == Dates.DateTime(2013, 12, 29)

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

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

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

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

firstday = Dates.Date(2014, 1, 1)
lastday = Dates.Date(2014, 12, 31)
for i = 0:364
    dt = firstday + Dates.Day(i)

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

    @test Dates.lastdayofyear(dt) == lastday
    @test Dates.lastdayofyear(DateTime(dt)) == DateTime(lastday)
end

# Adjusters
# Adjuster Constructors
@test Dates.Date(Dates.ismonday, 2014) == Dates.Date(2014, 1, 6)
@test Dates.Date(Dates.ismonday, 2014, 5) == Dates.Date(2014, 5, 5)

@test Dates.DateTime(Dates.ismonday, 2014) == Dates.DateTime(2014, 1, 6)
@test Dates.DateTime(Dates.ismonday, 2014, 5) == Dates.DateTime(2014, 5, 5)
@test Dates.DateTime(x->Dates.hour(x)==12, 2014, 5, 21) == Dates.DateTime(2014, 5, 21, 12)
@test Dates.DateTime(x->Dates.minute(x)==30, 2014, 5, 21, 12) == Dates.DateTime(2014, 5, 21, 12, 30)
@test Dates.DateTime(x->Dates.second(x)==30, 2014, 5, 21, 12, 30) == Dates.DateTime(2014, 5, 21, 12, 30, 30)
@test Dates.DateTime(x->Dates.millisecond(x)==500, 2014, 5, 21, 12, 30, 30) == Dates.DateTime(2014, 5, 21, 12, 30, 30, 500)

# tonext, toprev, tofirst, tolast
dt = Dates.Date(2014, 5, 21)
@test Dates.tonext(dt, Dates.Wed) == Dates.Date(2014, 5, 28)
@test Dates.tonext(dt, Dates.Wed; same=true) == dt
@test Dates.tonext(dt, Dates.Thu) == Dates.Date(2014, 5, 22)
@test Dates.tonext(dt, Dates.Fri) == Dates.Date(2014, 5, 23)
@test Dates.tonext(dt, Dates.Sat) == Dates.Date(2014, 5, 24)
@test Dates.tonext(dt, Dates.Sun) == Dates.Date(2014, 5, 25)
@test Dates.tonext(dt, Dates.Mon) == Dates.Date(2014, 5, 26)
@test Dates.tonext(dt, Dates.Tue) == Dates.Date(2014, 5, 27)
# No dayofweek function for out of range values
@test_throws KeyError Dates.tonext(dt, 8)

@test Dates.tonext(Dates.Date(0), Dates.Mon) == Dates.Date(0, 1, 3)

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

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

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

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

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

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

@test Dates.toprev(Dates.Date(0), Dates.Mon) == Dates.Date(-1, 12, 27)

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

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

@test Dates.tofirst(Dates.Date(0), Dates.Mon) == Dates.Date(0, 1, 3)

#tolast
@test Dates.tolast(dt, Dates.Mon) == Dates.Date(2014, 5, 26)
@test Dates.tolast(dt, Dates.Tue) == Dates.Date(2014, 5, 27)
@test Dates.tolast(dt, Dates.Wed) == Dates.Date(2014, 5, 28)
@test Dates.tolast(dt, Dates.Thu) == Dates.Date(2014, 5, 29)
@test Dates.tolast(dt, Dates.Fri) == Dates.Date(2014, 5, 30)
@test Dates.tolast(dt, Dates.Sat) == Dates.Date(2014, 5, 31)
@test Dates.tolast(dt, Dates.Sun) == Dates.Date(2014, 5, 25)

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

@test Dates.tolast(Dates.Date(0), Dates.Mon) == Dates.Date(0, 1, 31)

# filter (was recur)
startdate = Dates.Date(2014, 1, 1)
stopdate = Dates.Date(2014, 2, 1)
@test length(filter(x->true, startdate:stopdate)) == 32
@test length(filter(x->true, stopdate:Dates.Day(-1):startdate)) == 32

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

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

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

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

# Thanksgiving: 4th Thursday of November
thanksgiving = x->Dates.dayofweek(x) == Dates.Thu &&
                      Dates.month(x) == Dates.Nov &&
           Dates.dayofweekofmonth(x) == 4

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

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

@test Dates.toprev(d) do x
    thanksgiving(x)
end == Dates.Date(2013, 11, 28)

# Pittsburgh street cleaning
dr = Dates.Date(2014):Dates.Date(2015)
@test length(filter(dr) do x
    Dates.dayofweek(x) == Dates.Tue &&
    Dates.April < Dates.month(x) < Dates.Nov &&
    Dates.dayofweekofmonth(x) == 2
end) == 6

# U.S. Federal Holidays
newyears(y) = (y, 1, 1)
independenceday(y) = (y, 7, 4)
veteransday(y) = (y, 11, 11)
christmas(y) = (y, 12, 25)

isnewyears(dt) = Dates.yearmonthday(dt) == newyears(Dates.year(dt))
isindependenceday(dt) = Dates.yearmonthday(dt) == independenceday(Dates.year(dt))
isveteransday(dt) = Dates.yearmonthday(dt) == veteransday(Dates.year(dt))
ischristmas(dt) = Dates.yearmonthday(dt) == christmas(Dates.year(dt))
ismartinlutherking(dt) = Dates.dayofweek(dt) == Dates.Mon &&
    Dates.month(dt) == Dates.Jan && Dates.dayofweekofmonth(dt) == 3
ispresidentsday(dt) = Dates.dayofweek(dt) == Dates.Mon &&
    Dates.month(dt) == Dates.Feb && Dates.dayofweekofmonth(dt) == 3
# Last Monday of May
ismemorialday(dt) = Dates.dayofweek(dt) == Dates.Mon &&
                    Dates.month(dt) == Dates.May &&
                    Dates.dayofweekofmonth(dt) == Dates.daysofweekinmonth(dt)
islaborday(dt) = Dates.dayofweek(dt) == Dates.Mon &&
    Dates.month(dt) == Dates.Sep && Dates.dayofweekofmonth(dt) == 1
iscolumbusday(dt) = Dates.dayofweek(dt) == Dates.Mon &&
    Dates.month(dt) == Dates.Oct && Dates.dayofweekofmonth(dt) == 2
isthanksgiving(dt) = Dates.dayofweek(dt) == Dates.Thu &&
    Dates.month(dt) == Dates.Nov && Dates.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) = Dates.yearmonthday(dt) == easter(Dates.year(dt))

const 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

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

observed = filter(OBSERVEDHOLIDAYS, Dates.Date(1999):Dates.Date(2000))
@test length(observed) == 11
@test observed[10] == Dates.Date(1999, 12, 24)
@test observed[11] == Dates.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(Dates.Date(2014):Dates.Date(2015)) do x
    !(OBSERVEDHOLIDAYS(x) ||
    Dates.dayofweek(x) > 5)
end) == 251

# First day of the next month for each day of 2014
@test length([Dates.firstdayofmonth(i + Dates.Month(1))
    for i in Dates.Date(2014):Dates.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):Dates.Month(1):Date(2016)) do dt
    sum = 0
    for i = 1:7
        sum += Dates.dayofweek(dt) > 4 ? Dates.daysofweekinmonth(dt) : 0
        dt += Dates.Day(1)
    end
    return sum == 15
end) == 15 # On average, there's one of those months every year

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

r = filter(x->Dates.second(x) == 5, Dates.Time(0):Dates.Time(10))
@test length(r) == 600
@test first(r) == Dates.Time(0, 0, 5)
@test last(r) == Dates.Time(9, 59, 5)

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
# Types and functions
# -------------------

type HiddenLayer
    # The hidden layer in ELM
    #
    # Properties
    # ----------
    # `n_hidden_neurons` is an Integer representing the total hidden neurons
    # `weight_matrix` is the input to hidden layer weight matrix which is randomly chosen
    # `bias_vector` is the randomly chosen input to hidden layer bias vector
    # `act_func` is the activation function for hidden neurons

    n_hidden_neurons::Integer
    weight_matrix::Matrix{Float64}
    bias_vector::Vector{Float64}
    act_func::Function
end

type ExtremeLearningMachine
    # Extreme Learning Machine
    #
    # Properties
    # ----------
    # `hidden_layer` is the hidden layer inside this ELM
    # `output_weights` is the weight matrix calculated during training
    # `C` is the regularisation parameter that improves generalisation (not implemented as of now)

    hidden_layer::HiddenLayer
    output_weights::Matrix{Float64}
    C::Integer

    function ExtremeLearningMachine(n_hidden_neurons::Integer;
                                    C = 100)
	this = new()
	this.C = C

	weight_matrix = rand(0, 0)
	bias_vector = rand(n_hidden_neurons)

	this.hidden_layer = HiddenLayer(n_hidden_neurons, weight_matrix, bias_vector, sigmoid)

	this
    end
end

function sigmoid(x)
    # Sigmoid activation

    1 ./ (1 + exp(-x))
end

function find_activations(layer::HiddenLayer,
                          x::Matrix{Float64})
    # Calculates the activations of the hidden layer neurons
    #
    # Parameters
    # ----------
    # `layer` is the HiddenLayer with neurons
    # `x` is the input matrix
    #
    # Returns
    # -------
    # Activation matrix after passing through hidden layer

    n_observations = size(x)[1]

    act_matrix = zeros(layer.n_hidden_neurons, n_observations)

    for i = 1:n_observations
	act_matrix[:, i] = layer.act_func(layer.weight_matrix * x[i, :]' + layer.bias_vector)
    end

    act_matrix
end

function fit!(elm::ExtremeLearningMachine,
              x::Union(Matrix{Float64}, DataFrame),
              y::Union(Vector{Float64}, DataArray{Float64, 1}, DataArray{Int64, 1}))
    # Trains the elm using the given training data
    #
    # Parameters
    # ----------
    # `elm` the ELM to train
    # `x` input data
    # `y` output data

    if typeof(x) == DataFrame
        x_mat = float(array(x))
    else
        x_mat = x
    end

    if typeof(y) != Vector{Float64}
        y_vec = float(array(y))
    else
        y_vec = y
    end

    n_observations, n_inputs = size(x_mat)
    weight_matrix = rand(elm.hidden_layer.n_hidden_neurons, n_inputs) * 2 - 1
    elm.hidden_layer.weight_matrix = weight_matrix

    act_matrix = find_activations(elm.hidden_layer, x_mat)

    output_weights = y_vec' * pinv(act_matrix)
    elm.output_weights = output_weights
end

function predict(elm::ExtremeLearningMachine,
                 x::Union(Matrix{Float64}, DataFrame))
    # Predicts the output
    #
    # Parameters
    # ----------
    # `elm` the trained ELM
    # `x` input data to predict (Matrix)

    if typeof(x) == DataFrame
        x_mat = float(array(x))
    else
        x_mat = x
    end

    act_matrix = find_activations(elm.hidden_layer, x_mat)
    vec(elm.output_weights * act_matrix)
end


text:
import Oracle
import Oracle.OraNumbers

using Test

@testset "OraNumbers" begin
    @testset "basic numbers" begin
        z = zero(OraNumbers.OraNumber)
        o = one(OraNumbers.OraNumber)

        @test string(z) == "0"
        @test string(o) == "1"
        @test string(OraNumbers.OraNumber(nothing)) == "nothing"
    end

    @testset "comparison" begin
        z = zero(OraNumbers.OraNumber)
        o = one(OraNumbers.OraNumber)

        @test z == z
        @test o == o
        @test z != o
    end

    @testset "normalize" begin
        normalized_one = one(OraNumbers.OraNumber)
        non_normalized_one = OraNumbers.inc_exponent(normalized_one, 2)
        @test OraNumbers.isnormalized(normalized_one) == true
        @test OraNumbers.isnormalized(non_normalized_one) == false
        new_one = OraNumbers.normalize(non_normalized_one)
        @test OraNumbers.isnormalized(new_one)
        @test new_one == normalized_one
    end

    @testset "arithmetic" begin
        z = zero(OraNumbers.OraNumber)
        o = one(OraNumbers.OraNumber)

        @test z == -z
        @test z == -(-z)
        @test o == -(-o)

        @test string(-z) == "0"
        @test string(-o) == "-1"
    end
end


text:
module Nonconvex

const debugging = Ref(false)
const show_residuals = Ref(false)

export  Model,
        addvar!,
        add_ineq_constraint!,
        add_eq_constraint!,
        getmin,
        getmax,
        setmin!,
        setmax!,
        optimize,
        Workspace,
        MMA87,
        MMA02,
        MMALag,
        IpoptAlg,
        NLoptAlg,
        AugLag,
        PercivalAlg,
        JuniperIpoptAlg,
        KKTCriteria,
        IpoptCriteria,
        FunctionWrapper,
        MMAOptions,
        IpoptOptions,
        NLoptOptions,
        AugLagOptions,
        PercivalOptions,
        JuniperIpoptOptions,
        Tolerance

using Parameters, Zygote, ChainRulesCore, ForwardDiff, MathOptInterface
using Ipopt, NLopt, ADNLPModels, Percival, NLPModelsModifiers, JuMP
using LinearAlgebra, Setfield, Requires, SparseArrays, Reexport
using Juniper
using Optim: Optim, AbstractOptimizer

@reexport using LinearAlgebra

abstract type Workspace end

# General

include("utilities/params.jl")
include("functions/functions.jl")
include("functions/value_jacobian.jl")
include("functions/function_docs.jl")
include("functions/counting_function.jl")
include("functions/aggregations.jl")

# Models

include("models/model.jl")
include("models/model_docs.jl")
include("mma_approximation/mma_approx.jl")
include("mma_approximation/xmma_approx.jl")
include("mma_approximation/mma_approx_docs.jl")
include("models/mma_model.jl")
include("models/dual_model.jl")
include("models/mmalag_model.jl")

# MMA

include("utilities/convergence.jl")
include("utilities/callbacks.jl")
include("utilities/options.jl")
include("algorithms/mma_algorithm.jl")

# MMA-AugLag2

include("algorithms/stoch_optimizers.jl")
include("algorithms/nonstoch_optimizers.jl")
include("algorithms/ammal.jl")

# Augmented Lagrangian

include("models/auglag_model.jl")
include("algorithms/auglag_algorithm.jl")

# Wrappers

include("wrappers/ipopt.jl")
include("wrappers/nlopt.jl")
include("wrappers/percival.jl")
include("wrappers/moi.jl")
include("wrappers/juniper.jl")

end


text:
"""
    struct JSONModel <: MetabolicModel
        json::Dict{String,Any}
    end

A struct used to store the contents of a JSON model, i.e. a model read from a
file ending with `.json`. These model files typically store all the model
parameters in arrays of JSON objects (i.e. Julia dictionaries).

Usually, not all of the fields of the input JSON can be easily represented when
converting to other models, care should be taken to avoid losing information.

Direct work on this precise model type is not very efficient, as the accessor
functions need to repeatedly find the information in the JSON tree. This gets
very slow especially if calling many accessor functions sequentially. To avoid
that, convert to e.g. [`StandardModel`](@ref) as soon as possible.

# Example
````
model = load_json_model("some_model.json")
model.json # see the actual underlying JSON
reactions(model) # see the list of reactions
````
"""
struct JSONModel <: MetabolicModel
    json::Dict{String,Any}
end

# helper access macros, see examples below for usage
macro _json_sectionkey(namesConstant::Symbol, error)
    esc(:(
        begin
            _key = _guesskey(keys(model.json), _constants.keynames.$namesConstant)
            if isnothing(_key)
                $error
            end
            _key
        end
    ))
end

macro _json_section(namesConstant::Symbol, error)
    esc(:(
        begin
            _key = @_json_sectionkey $namesConstant ($error)
            model.json[_key]
        end
    ))
end

macro _json_section_firstid(namesConstant::Symbol, id::Symbol, error)
    return esc(:(
        begin
            _s = @_json_section $namesConstant ($error)
            _s[findfirst(x -> x["id"] == $id, _s)]
        end
    ))
end

function _parse_annotations(x)::Annotations
    Annotations([k => if typeof(v) == String
        [v]
    else
        convert(Vector{String}, v)
    end for (k, v) in x])
end

_parse_notes(x)::Notes = _parse_annotations(x)

"""
    reactions(model::JSONModel)

Extract reaction names (stored as `.id`) from JSON model.
"""
function reactions(model::JSONModel)
    rs = @_json_section rxns throw(
        DomainError(keys(model.json), "JSON model has no reaction keys"),
    )

    return [string(get(rs[i], "id", "rxn$i")) for i in eachindex(rs)]
end

"""
    metabolites(model::JSONModel)

Extract metabolite names (stored as `.id`) from JSON model.
"""
function metabolites(model::JSONModel)
    ms = @_json_section mets throw(
        DomainError(keys(model.json), "JSON model has no metabolite keys"),
    )

    return [string(get(ms[i], "id", "met$i")) for i in eachindex(ms)]
end

"""
    genes(model::JSONModel)

Extract gene names from a JSON model.
"""
function genes(model::JSONModel)
    gs = @_json_section genes return []

    return [string(get(gs[i], "id", "gene$i")) for i in eachindex(gs)]
end

"""
    stoichiometry(model::JSONModel)

Get the stoichiometry. Assuming the information is stored in reaction object
under key `.metabolites`.
"""
function stoichiometry(model::JSONModel)
    rxn_ids = reactions(model)
    met_ids = metabolites(model)

    rs = @_json_section rxns throw(
        DomainError(keys(model.json), "JSON model has no reaction keys"),
    )

    S = SparseArrays.spzeros(length(met_ids), length(rxn_ids))
    for (i, rid) in enumerate(rxn_ids)
        r = rs[findfirst(x -> x["id"] == rid, rs)]
        for (met_id, coeff) in r["metabolites"]
            j = findfirst(==(met_id), met_ids)
            if isnothing(j)
                throw(
                    DomainError(
                        met_id,
                        "Unknown metabolite found in stoichiometry of $(rxn_ids[i])",
                    ),
                )
            end
            S[j, i] = coeff
        end
    end
    return S
end

"""
    bounds(model::JSONModel)

Get the bounds for reactions, assuming the information is stored in
`.lower_bound` and `.upper_bound`.
"""
function bounds(model::JSONModel)
    rs = @_json_section rxns return (
        sparse(fill(-_constants.default_reaction_bound, n_reactions(model))),
        sparse(fill(_constants.default_reaction_bound, n_reactions(model))),
    )
    return (
        sparse([get(rxn, "lower_bound", -_constants.default_reaction_bound) for rxn in rs]),
        sparse([get(rxn, "upper_bound", _constants.default_reaction_bound) for rxn in rs]),
    )
end

"""
    objective(model::JSONModel)

Collect `.objective_coefficient` keys from model reactions.
"""
function objective(model::JSONModel)
    rs = @_json_section rxns return spzeros(n_reactions(model))

    return sparse([float(get(rxn, "objective_coefficient", 0.0)) for rxn in rs])
end

"""
    reaction_gene_associaton(model::JSONModel, rid::String)

Parses the `.gene_reaction_rule` from reactions.
"""
function reaction_gene_association(model::JSONModel, rid::String)
    rxn = @_json_section_firstid rxns rid return nothing
    return _maybemap(_parse_grr, get(rxn, "gene_reaction_rule", nothing))
end

"""
    reaction_subsystem(model::JSONModel, rid::String)

Parses the `.subsystem` out from reactions.
"""
function reaction_subsystem(model::JSONModel, rid::String)
    rxn = @_json_section_firstid rxns rid return nothing
    return get(rxn, "subsystem", nothing)
end

"""
    metabolite_formula(model::JSONModel, mid::String)

Parse and return the metabolite `.formula`
"""
function metabolite_formula(model::JSONModel, mid::String)
    met = @_json_section_firstid mets mid return nothing
    return _maybemap(_parse_formula, get(met, "formula", nothing))
end

"""
    metabolite_charge(model::JSONModel, mid::String)

Return the metabolite `.charge`
"""
function metabolite_charge(model::JSONModel, mid::String)
    met = @_json_section_firstid mets mid return nothing
    return get(met, "charge", 0)
end

"""
    metabolite_compartment(model::JSONModel, mid::String)

Return the metabolite `.compartment`
"""
function metabolite_compartment(model::JSONModel, mid::String)
    met = @_json_section_firstid mets mid return nothing
    return get(met, "compartment", nothing)
end

"""
    gene_annotations(model::JSONModel, gid::String)::Annotations

Gene annotations from the [`JSONModel`](@ref).
"""
function gene_annotations(model::JSONModel, gid::String)::Annotations
    gene = @_json_section_firstid genes gid return Dict()
    _maybemap(_parse_annotations, get(gene, "annotation", nothing))
end

"""
    gene_notes(model::JSONModel, gid::String)::Notes

Gene notes from the [`JSONModel`](@ref).
"""
function gene_notes(model::JSONModel, gid::String)::Notes
    gene = @_json_section_firstid genes gid return Dict()
    _maybemap(_parse_notes, get(gene, "notes", nothing))
end

"""
    reaction_annotations(model::JSONModel, rid::String)::Annotations

Reaction annotations from the [`JSONModel`](@ref).
"""
function reaction_annotations(model::JSONModel, rid::String)::Annotations
    rxn = @_json_section_firstid rxns rid return Dict()
    _maybemap(_parse_annotations, get(rxn, "annotation", nothing))
end

"""
    reaction_notes(model::JSONModel, rid::String)::Notes

Reaction notes from the [`JSONModel`](@ref).
"""
function reaction_notes(model::JSONModel, rid::String)::Notes
    rxn = @_json_section_firstid rxns rid return Dict()
    _maybemap(_parse_notes, get(rxn, "notes", nothing))
end

"""
    metabolite_annotations(model::JSONModel, mid::String)::Annotations

Metabolite annotations from the [`JSONModel`](@ref).
"""
function metabolite_annotations(model::JSONModel, mid::String)::Annotations
    met = @_json_section_firstid mets mid return Dict()
    _maybemap(_parse_annotations, get(met, "annotation", nothing))
end

"""
    metabolite_notes(model::JSONModel, mid::String)::Notes

Metabolite notes from the [`JSONModel`](@ref).
"""
function metabolite_notes(model::JSONModel, mid::String)::Notes
    met = @_json_section_firstid mets mid return Dict()
    _maybemap(_parse_notes, get(met, "notes", nothing))
end

"""
    reaction_stoichiometry(model::JSONModel, rxn_id::String)::Dict{String, Float64}

Return the reaction equation of reaction with id `rxn_id` in model. The reaction
equation maps metabolite ids to their stoichiometric coefficients.
"""
function reaction_stoichiometry(m::JSONModel, rxn_id::String)::Dict{String, Float64}
    ind = findfirst(x -> x["id"] == rxn_id, m.json["reactions"])
    m.json["reactions"][ind]["metabolites"]
end

"""
    Base.convert(::Type{JSONModel}, mm::MetabolicModel)

Convert any [`MetabolicModel`](@ref) to [`JSONModel`](@ref).
"""
function Base.convert(::Type{JSONModel}, mm::MetabolicModel)
    if typeof(mm) == JSONModel
        return mm
    end

    rxn_ids = reactions(mm)
    met_ids = metabolites(mm)
    gene_ids = genes(mm)
    S = stoichiometry(mm)
    lbs, ubs = bounds(mm)
    ocs = objective(mm)

    json_model = JSONModel(Dict{String,Any}())
    json = Dict{String,Any}()
    json["id"] = "model" # default

    #TODO: add notes, names and similar fun stuff when they are available

    json[first(_constants.keynames.genes)] = [
        Dict([
            "id" => gid,
            "annotation" => gene_annotations(mm, gid),
            "notes" => gene_notes(mm, gid),
        ],) for gid in gene_ids
    ]

    json[first(_constants.keynames.mets)] = [
        Dict([
            "id" => mid,
            "formula" => _maybemap(_unparse_formula, metabolite_formula(mm, mid)),
            "charge" => metabolite_charge(mm, mid),
            "compartment" => metabolite_compartment(mm, mid),
            "annotation" => metabolite_annotations(mm, mid),
            "notes" => metabolite_notes(mm, mid),
        ]) for mid in met_ids
    ]

    json[first(_constants.keynames.rxns)] = [
        begin
            res = Dict{String,Any}()
            res["id"] = rid
            res["subsystem"] = reaction_subsystem(mm, rid)
            res["annotation"] = reaction_annotations(mm, rid)
            res["notes"] = reaction_notes(mm, rid)

            grr = reaction_gene_association(mm, rid)
            if !isnothing(grr)
                res["gene_reaction_rule"] = _unparse_grr(String, grr)
            end

            res["lower_bound"] = lbs[ri]
            res["upper_bound"] = ubs[ri]
            res["objective_coefficient"] = ocs[ri]
            I, V = findnz(S[:, ri])
            res["metabolites"] =
                Dict{String,Float64}([met_ids[ii] => vv for (ii, vv) in zip(I, V)])
            res
        end for (ri, rid) in enumerate(rxn_ids)
    ]

    return JSONModel(json)
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -1.263105417837706 <= q <= 2.4766737340949954)

                     add_NL_constraint(m, :(swish(0.7635571507341847 + -0.35250618093605723*swish(0.5051094873509578 + -0.17934163301535788*swish(-0.331542222824047 + 0.7450661875335634*swish(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + -0.10668129539960924*swish(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))) + -0.9610191438924418*swish(0.6781688636512797 + -0.22580923517577522*swish(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + 0.2965219706242044*swish(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4])))) + 0.1709251594821839*swish(0.058205314222445015 + 0.7954955842799771*swish(-0.331542222824047 + 0.7450661875335634*swish(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + -0.10668129539960924*swish(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))) + -0.9762293687340509*swish(0.6781688636512797 + -0.22580923517577522*swish(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + 0.2965219706242044*swish(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))))) + swish(-0.29663194049169395 + 0.013502946151282558*swish(0.5051094873509578 + -0.17934163301535788*swish(-0.331542222824047 + 0.7450661875335634*swish(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + -0.10668129539960924*swish(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))) + -0.9610191438924418*swish(0.6781688636512797 + -0.22580923517577522*swish(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + 0.2965219706242044*swish(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4])))) + -0.2459643040595232*swish(0.058205314222445015 + 0.7954955842799771*swish(-0.331542222824047 + 0.7450661875335634*swish(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + -0.10668129539960924*swish(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))) + -0.9762293687340509*swish(0.6781688636512797 + -0.22580923517577522*swish(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + 0.2965219706242044*swish(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# Also needs packages Parameters and VAWTools

plotyes = true
using PyPlot
using WhereTheWaterFlows
if !(@isdefined WWF)
    const WWF = WhereTheWaterFlows
end

"A artificial DEM"
function ele(x, y; withpit=false, randfac=0.0)
    out = - (x^2 - 1)^2 - (x^2*y - x - 1)^2 + 6 + 0.1*x + 3*y
    if withpit
        out -= 2*exp(-(x^2 + y^2)*50)
    end
    out += randfac*randn(size(out))
    return out<0 ? 0.0 : out
end
dx = 0.01
xs = -1.5:dx:1
ys = -0.5:dx:3.0
dem = ele.(xs, ys', randfac=0.1, withpit=true);
plotyes && WWF.heatmap(xs,ys,dem)
area, slen, dir, nout, nin, pits  = WWF.waterflows(dem, fillpits=true);
@assert size(dem)==(length(xs), length(ys))
plotyes && WWF.plotit(xs,ys,dem)
plotyes && WWF.plotarea(xs, ys, area, pits)

c, bnds = WWF.catchments(dir, pits);
plotyes && WWF.heatmap(xs, ys, c)

demf = WWF.fill_dem(dem, pits, dir) #, small=1e-6)
plotyes && WWF.heatmap(xs,ys, demf.-dem)
# WWF.plotarea(xs, ys, demf) # will show some artefacts in filled-in area unless `small` is set

## Sorry, this only works on my machine...
import VAWTools
using Parameters
gorner = VAWTools.read_agr("/home/mauro/model_runs/glads/gorner/input/topo_gorner/data/dhm_gorezg2007.grid")
@unpack x,y,v = gorner;
#v = v + rand(size(v)...)*1e-2;
#vv = VAWTools.boxcar(v, 4); #  + rand(size(v)...)*1e-1; # add some noise to kill lakes
dem = v;

area, slen, dir, nout, nin, pits  = WWF.waterflows(dem);
c, bnds = WWF.catchments(dir, pits);
plotyes && WWF.plotarea(x,y,area, pits)

demf = WWF.fill_dem(dem, pits, dir)
plotyes && WWF.heatmap(x,y, demf.-dem)

dir_, nin_, nout_, pits_, c_, bnds_ = WWF.drainpits(dem, dir, nin, nout, pits);
area_, slen_ = WWF._waterflows(dir_, nout_, nin_, pits_);
WWF.plotarea(x, y, area_, pits_)

using BenchmarkTools
@btime WWF.d8dir_feature(v) # 51.695 ms (388933 allocations: 83.48 MiB)
#@btime WWF.d8dir_feature1(v)

if false # small
    @unpack x,y,v = VAWTools.downsample(gorner, 40);
    WWF.plotarea(x,y,v, prefn=x->x);
    dir, nout, nin, pits = WWF.d8dir_feature(v);
    area, slen, dir, nout, nin, pits  = WWF.waterflows(v);
    c, bnds = WWF.catchments(dir, pits);

    figure(); WWF.heatmap(x,y,slen)
    figure(); WWF.heatmap(x,y,nin)
    figure(); WWF.heatmap(x,y,c)
    figure(); WWF.plotdir(x,y,dir)

    dir_, nin_, nout_, pits_, c_, bnds_ = WWF.drainpits(v, dir, nin, nout, pits);
    area_, slen_ = WWF.waterflows(dir_, nout_, nin_, pits_);
    #figure(); plotdir(x,y,dir_); figure(); plotdir(x,y,dir)
    WWF.plotarea(x,y,area_, pits_)
end


text:
## reductions ##

###### Functors ######

# Note that functors are merely used as internal machinery to enhance code reuse.
# They are not exported.
# When function arguments can be inlined, the use of functors can be removed.

abstract Func{N}

type IdFun <: Func{1} end
type AbsFun <: Func{1} end
type Abs2Fun <: Func{1} end
type ExpFun <: Func{1} end
type LogFun <: Func{1} end

type AddFun <: Func{2} end
type MulFun <: Func{2} end
type AndFun <: Func{2} end
type OrFun <: Func{2} end

evaluate(::IdFun, x) = x
evaluate(::AbsFun, x) = abs(x)
evaluate(::Abs2Fun, x) = abs2(x)
evaluate(::ExpFun, x) = exp(x)
evaluate(::LogFun, x) = log(x)
evaluate(f::Callable, x) = f(x)

evaluate(::AddFun, x, y) = x + y
evaluate(::MulFun, x, y) = x * y
evaluate(::AndFun, x, y) = x & y
evaluate(::OrFun, x, y) = x | y
evaluate(f::Callable, x, y) = f(x, y)

###### Generic reduction functions ######

# Note that getting type-stable results from reduction functions,
# or at least having type-stable loops, is nontrivial (#6069).

## foldl

function _foldl(op, v0, itr, i)
    if done(itr, i)
        return v0
    else
        (x, i) = next(itr, i)
        v = evaluate(op, v0, x)
        while !done(itr, i)
            (x, i) = next(itr, i)
            v = evaluate(op, v, x)
        end
        return v
    end
end

function foldl(op::Callable, v0, itr, i)
    is(op, +) && return _foldl(AddFun(), v0, itr, i)
    is(op, *) && return _foldl(MulFun(), v0, itr, i)
    is(op, &) && return _foldl(AndFun(), v0, itr, i)
    is(op, |) && return _foldl(OrFun(), v0, itr, i)
    return _foldl(op, v0, itr, i)
end

foldl(op::Callable, v0, itr) = foldl(op, v0, itr, start(itr))

function foldl(op::Callable, itr)
    i = start(itr)
    done(itr, i) && error("Argument is empty.")
    (v0, i) = next(itr, i)
    return foldl(op, v0, itr, i)
end

## foldr

function foldr(op::Callable, v0, itr, i=endof(itr))
    # use type stable procedure
    if i == 0
        return v0
    else
        v = op(itr[i], v0)
        while i > 1
            x = itr[i -= 1]
            v = op(x, v)
        end
        return v
    end
end

foldr(op::Callable, itr) = (i = endof(itr); foldr(op, itr[i], itr, i-1))

## reduce

reduce(op::Callable, v, itr) = foldl(op, v, itr)

function reduce(op::Callable, itr) # this is a left fold
    if is(op, +)
        return sum(itr)
    elseif is(op, *)
        return prod(itr)
    elseif is(op, |)
        return any(itr)
    elseif is(op, &)
        return all(itr)
    end
    return foldl(op, itr)
end

# pairwise reduction, requires n > 1 (to allow type-stable loop)
function r_pairwise(op::Callable, A::AbstractArray, i1,n)
    if n < 128
        @inbounds v = op(A[i1], A[i1+1])
        for i = i1+2:i1+n-1
            @inbounds v = op(v,A[i])
        end
        return v
    else
        n2 = div(n,2)
        return op(r_pairwise(op,A, i1,n2), r_pairwise(op,A, i1+n2,n-n2))
    end
end

function reduce(op::Callable, A::AbstractArray)
    n = length(A)
    n == 0 ? error("argument is empty") : n == 1 ? A[1] : r_pairwise(op,A, 1,n)
end

function reduce(op::Callable, v0, A::AbstractArray)
    n = length(A)
    n == 0 ? v0 : n == 1 ? op(v0, A[1]) : op(v0, r_pairwise(op,A, 1,n))
end

###### Specific reduction functions ######

## in & contains

function in(x, itr)
    for y in itr
        if y == x
            return true
        end
    end
    return false
end
const ∈ = in
∉(x, itr)=!∈(x, itr)
∋(itr, x)= ∈(x, itr)
∌(itr, x)=!∋(itr, x)

function contains(itr, x)
    depwarn("contains(collection, item) is deprecated, use in(item, collection) instead", :contains)
    in(x, itr)
end

function contains(eq::Function, itr, x)
    for y in itr
        if eq(y, x)
            return true
        end
    end
    return false
end

## countnz & count

function countnz(itr)
    n = 0
    for x in itr
        if x != 0
            n += 1
        end
    end
    return n
end

function countnz(a::AbstractArray)
    n = 0
    for i = 1:length(a)
        @inbounds x = a[i]
        if x != 0
            n += 1
        end
    end
    return n
end

function count(pred::Function, itr)
    n = 0
    for x in itr
        if pred(x)
            n += 1
        end
    end
    return n
end

## sum

# result type inference for sum

sumtype{T}(::Type{T}) = typeof(zero(T) + zero(T))
sumzero{T}(::Type{T}) = zero(T) + zero(T)
addzero(x) = x + zero(x)

typealias SumResultNumber Union(Uint,Uint64,Uint128,Int,Int64,Int128,Float32,Float64,Complex64,Complex128)

sumtype{T<:SumResultNumber}(::Type{T}) = T
sumzero{T<:SumResultNumber}(::Type{T}) = zero(T)
addzero(x::SumResultNumber) = x

sumzero{T<:AbstractArray}(::Type{T}) = error("Summing over an empty collection of arrays is not allowed.")
addzero(a::AbstractArray) = a

# general sum over iterables

function _sum(f, itr, s)  # deal with non-empty cases
    # pre-condition: s = start(itr) && !done(itr, s)
    (v, s) = next(itr, s)
    done(itr, s) && return addzero(evaluate(f, v)) # adding zero for type stability
    # specialize for length > 1 to have type-stable loop
    (x, s) = next(itr, s)
    result = evaluate(f, v) + evaluate(f, x)
    while !done(itr, s)
        (x, s) = next(itr, s)
        result += evaluate(f, x)
    end
    return result
end

function sum(itr)
    s = start(itr)
    if done(itr, s)
        if applicable(eltype, itr)
            return sumzero(eltype(itr))
        else
            throw(ArgumentError("sum(itr) is undefined for empty collections; instead, do isempty(itr) ? z : sum(itr), where z is the correct type of zero for your sum"))
        end
    end
    _sum(IdFun(), itr, s)
end

function sum(f::Union(Function,Func{1}), itr)
    s = start(itr)
    done(itr, s) && error("Argument is empty.")
    _sum(f, itr, s)
end

sum(x::Number) = x
sum(A::AbstractArray{Bool}) = countnz(A)

sumabs(itr) = sum(AbsFun(), itr)
sumabs2(itr) = sum(Abs2Fun(), itr)

sumabs(x::Number) = abs(x)
sumabs2(x::Number) = abs2(x)

# Note: sum_seq uses four accumulators, so each accumulator gets at most 256 numbers
sum_pairwise_blocksize(f) = 1024

# a fast implementation of sum in sequential order (from left to right).
# to allow type-stable loops, requires length > 1
function sum_seq(f, a::AbstractArray, ifirst::Int, ilast::Int)

    @inbounds if ifirst + 6 >= ilast  # length(a) < 8
        i = ifirst
        s = evaluate(f, a[i]) + evaluate(f, a[i+1])
        i = i+1
        while i < ilast
            s += evaluate(f, a[i+=1])
        end
        return s

    else # length(a) >= 8

        # more effective utilization of the instruction
        # pipeline through manually unrolling the sum
        # into four-way accumulation. Benchmark shows
        # that this results in about 2x speed-up.

        s1 = evaluate(f, a[ifirst]) + evaluate(f, a[ifirst + 4])
        s2 = evaluate(f, a[ifirst + 1]) + evaluate(f, a[ifirst + 5])
        s3 = evaluate(f, a[ifirst + 2]) + evaluate(f, a[ifirst + 6])
        s4 = evaluate(f, a[ifirst + 3]) + evaluate(f, a[ifirst + 7])

        i = ifirst + 8
        il = ilast - 3
        while i <= il
            s1 += evaluate(f, a[i])
            s2 += evaluate(f, a[i+1])
            s3 += evaluate(f, a[i+2])
            s4 += evaluate(f, a[i+3])
            i += 4
        end

        while i <= ilast
            s1 += evaluate(f, a[i])
            i += 1
        end

        return s1 + s2 + s3 + s4
    end
end

# Pairwise (cascade) summation of A[i1:i1+n-1], which has O(log n) error growth
# [vs O(n) for a simple loop] with negligible performance cost if
# the base case is large enough.  See, e.g.:
#        http://en.wikipedia.org/wiki/Pairwise_summation
#        Higham, Nicholas J. (1993), "The accuracy of floating point
#        summation", SIAM Journal on Scientific Computing 14 (4): 783–799.
# In fact, the root-mean-square error growth, assuming random roundoff
# errors, is only O(sqrt(log n)), which is nearly indistinguishable from O(1)
# in practice.  See:
#        Manfred Tasche and Hansmartin Zeuner, Handbook of
#        Analytic-Computational Methods in Applied Mathematics (2000).
#
# sum_impl requires length(a) > 1
#
function sum_impl(f, a::AbstractArray, ifirst::Int, ilast::Int)
    if ifirst + sum_pairwise_blocksize(f) >= ilast
        sum_seq(f, a, ifirst, ilast)
    else
        imid = (ifirst + ilast) >>> 1
        sum_impl(f, a, ifirst, imid) + sum_impl(f, a, imid+1, ilast)
    end
end
sum_impl{T<:Integer}(f::Union(IdFun,AbsFun,Abs2Fun), a::AbstractArray{T}, ifirst::Int, ilast::Int) =
    sum_seq(f, a, ifirst, ilast)

function sum(f::Union(Function,Func{1}), a::AbstractArray)
    n = length(a)
    n == 0 && error("Argument is empty.")
    n == 1 && return addzero(evaluate(f, a[1]))
    sum_impl(f, a, 1, n)
end

for (fname, func, cutoff) in ((:sum, :IdFun, 16), (:sumabs, :AbsFun, 32), (:sumabs2, :Abs2Fun, 32))
    @eval function $fname{T}(a::AbstractArray{T})
        n = length(a)
        n == 0 && return addzero(evaluate($func(), zero(T)))
        n == 1 && return addzero(evaluate($func(), a[1]))
        if n < $cutoff
            # It is important that this is inlined to provide good
            # performance for small inputs
            @inbounds s = evaluate($func(), a[1]) + evaluate($func(), a[2])
            for i = 3:length(a)
                @inbounds s += evaluate($func(), a[i])
            end
            return s
        end
        sum_impl($func(), a, 1, n)
    end
end

# Kahan (compensated) summation: O(1) error growth, at the expense
# of a considerable increase in computational expense.
function sum_kbn{T<:FloatingPoint}(A::AbstractArray{T})
    n = length(A)
    if n == 0
        return sumzero(T)
    end
    s = addzero(A[1])
    c = sumzero(T)
    for i in 2:n
        @inbounds Ai = A[i]
        t = s + Ai
        if abs(s) >= abs(Ai)
            c += ((s-t) + Ai)
        else
            c += ((Ai-t) + s)
        end
        s = t
    end
    s + c
end

## prod

prodtype{T}(::Type{T}) = typeof(zero(T) * zero(T))
prodone{T}(::Type{T}) = one(T) * one(T)
multone(x) = x * one(x)

function _prod(f, itr, s)
    (x, s) = next(itr, s)
    v = evaluate(f, x)
    done(itr, s) && return multone(v) # multiplying by one for type stability
    # specialize for length > 1 to have type-stable loop
    (x, s) = next(itr, s)
    result = v * evaluate(f, x)
    while !done(itr, s)
        (x, s) = next(itr, s)
        result *= evaluate(f, x)
    end
    return result
end

function prod(itr)
    s = start(itr)
    if done(itr, s)
        if applicable(eltype, itr)
            T = eltype(itr)
            return prodone(T)
        else
            throw(ArgumentError("prod(itr) is undefined for empty collections; instead, do isempty(itr) ? o : prod(itr), where o is the correct type of identity for your product"))
        end
    end
    _prod(IdFun(), itr, s)
end

function prod(f::Function, itr)
    s = start(itr)
    done(itr, s) && error("Argument is empty.")
    _prod(f, itr, s)
end

prod(x::Number) = x

prod(A::AbstractArray{Bool}) =
    error("use all() instead of prod() for boolean arrays")

function prod_impl{T}(f, A::AbstractArray{T}, first::Int, last::Int)
    # pre-condition: last > first
    i = first
    @inbounds v = evaluate(f, A[i])
    @inbounds result = v * evaluate(f, A[i+=1])
    while i < last
        @inbounds result *= evaluate(f, A[i+=1])
    end
    return result
end

function prod{T}(A::AbstractArray{T})
    n = length(A)
    n == 0 && return prodone(T)
    n == 1 && return multone(A[1])
    prod_impl(IdFun(), A, 1, n)
end

function prod(f::Function, A::AbstractArray)
    n = length(A)
    n == 0 && error("Argument is empty.")
    n == 1 && return multone(evaluate(f, A[1]))
    prod_impl(f, A, 1, n)
end

## maximum & minimum

function maximum(f::Union(Function,Func{1}), itr)
    s = start(itr)
    if done(itr, s)
        error("argument is empty")
    end
    (x, s) = next(itr, s)
    v = evaluate(f, x)
    while !done(itr, s)
        (x, s) = next(itr, s)
        v = scalarmax(v, evaluate(f, x))
    end
    return v
end

function minimum(f::Union(Function,Func{1}), itr)
    s = start(itr)
    if done(itr, s)
        error("argument is empty")
    end
    (x, s) = next(itr, s)
    v = evaluate(f, x)
    while !done(itr, s)
        (x, s) = next(itr, s)
        v = scalarmin(v, evaluate(f, x))
    end
    return v
end

maximum(itr) = maximum(IdFun(), itr)
minimum(itr) = minimum(IdFun(), itr)

function maximum_impl(f, A::AbstractArray, first::Int, last::Int)
    # locate the first non NaN number
    v = evaluate(f, A[first])
    i = first + 1
    while v != v && i <= last
        @inbounds v = evaluate(f, A[i])
        i += 1
    end

    while i <= last
        @inbounds x = evaluate(f, A[i])
        if x > v
            v = x
        end
        i += 1
    end
    v
end

function minimum_impl(f, A::AbstractArray, first::Int, last::Int)
    # locate the first non NaN number
    v = evaluate(f, A[first])
    i = first + 1
    while v != v && i <= last
        @inbounds v = evaluate(f, A[i])
        i += 1
    end

    while i <= last
        @inbounds x = evaluate(f, A[i])
        if x < v
            v = x
        end
        i += 1
    end
    v
end

function maximum(f::Union(Function,Func{1}), A::AbstractArray)
    n = length(A)
    n == 0 && error("Argument is empty.")
    n == 1 && return evaluate(f, A[1])
    maximum_impl(f, A, 1, n)
end

function minimum(f::Union(Function,Func{1}), A::AbstractArray)
    n = length(A)
    n == 0 && error("Argument is empty.")
    n == 1 && return evaluate(f, A[1])
    minimum_impl(f, A, 1, n)
end

maximum(A::AbstractArray) = maximum(IdFun(), A)
minimum(A::AbstractArray) = minimum(IdFun(), A)

minabs(A::AbstractArray) = minimum(AbsFun(), A)

# maxabs accepts empty array
function maxabs(A::AbstractArray)
    n = length(A)
    n == 0 && return abs(zero(T))
    n == 1 && return abs(A[1])
    maximum_impl(AbsFun(), A, 1, n)
end

maximum(x::Real) = x
minimum(x::Real) = x
maxabs(x::Number) = abs(x)
minabs(x::Number) = abs(x)

## extrema

extrema(r::Range) = (minimum(r), maximum(r))

function extrema(itr)
    s = start(itr)
    if done(itr, s)
        error("argument is empty")
    end
    (v, s) = next(itr, s)
    vmin = v
    vmax = v
    while !done(itr, s)
        (x, s) = next(itr, s)
        if x == x
            if x > vmax
                vmax = x
            elseif x < vmin
                vmin = x
            end
        end
    end
    return (vmin, vmax)
end

function extrema{T<:Real}(A::AbstractArray{T})
    if isempty(A); error("argument must not be empty"); end
    n = length(A)

    # locate the first non NaN number
    v = A[1]
    i = 2
    while v != v && i <= n
        @inbounds v = A[i]
        i += 1
    end

    vmin = v
    vmax = v
    while i <= n
        @inbounds v = A[i]
        if v > vmax
            vmax = v
        elseif v < vmin
            vmin = v
        end
        i += 1
    end

    return (vmin, vmax)
end

## all & any

function all(itr)
    for x in itr
        if !x
            return false
        end
    end
    return true
end

function any(itr)
    for x in itr
        if x
            return true
        end
    end
    return false
end

function any(pred::Union(Function,Func{1}), itr)
    for x in itr
        if pred(x)
            return true
        end
    end
    return false
end

function all(pred::Union(Function,Func{1}), itr)
    for x in itr
        if !pred(x)
            return false
        end
    end
    return true
end

###### mapreduce ######

function mapreduce(f::Callable, op::Callable, itr)
    s = start(itr)
    if done(itr, s)
        error("argument is empty")
    end
    (x, s) = next(itr, s)
    v = f(x)
    if done(itr, s)
        return v
    else # specialize for length > 1 to have a hopefully type-stable loop
        (x, s) = next(itr, s)
        result = op(v, f(x))
        while !done(itr, s)
            (x, s) = next(itr, s)
            result = op(result, f(x))
        end
        return result
    end
end

function mapreduce(f::Callable, op::Callable, v0, itr)
    v = v0
    for x in itr
        v = op(v,f(x))
    end
    return v
end

# pairwise reduction, requires n > 1 (to allow type-stable loop)
function mr_pairwise(f::Callable, op::Callable, A::AbstractArray, i1,n)
    if n < 128
        @inbounds v = op(f(A[i1]), f(A[i1+1]))
        for i = i1+2:i1+n-1
            @inbounds v = op(v,f(A[i]))
        end
        return v
    else
        n2 = div(n,2)
        return op(mr_pairwise(f,op,A, i1,n2), mr_pairwise(f,op,A, i1+n2,n-n2))
    end
end
function mapreduce(f::Callable, op::Callable, A::AbstractArray)
    n = length(A)
    n == 0 ? error("argument is empty") : n == 1 ? f(A[1]) : mr_pairwise(f,op,A, 1,n)
end
function mapreduce(f::Callable, op::Callable, v0, A::AbstractArray)
    n = length(A)
    n == 0 ? v0 : n == 1 ? op(v0, f(A[1])) : op(v0, mr_pairwise(f,op,A, 1,n))
end

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7, x8, x9)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7, x8, x9]);F = [361*x1^2*x4^2+1748*x1^2*x4*x9+2209*x1^2*x7^2+2116*x1^2*x8^2+2116*x1^2*x9^2-7896*x1*x3*x4*x7+1026*x1*x3*x4*x8-7050*x1*x3*x6*x7-1034*x1*x3*x7^2-2392*x1*x3*x7*x8+2484*x1*x3*x8*x9+1330*x1*x4^3+3220*x1*x4^2*x9+1254*x1*x4*x6*x8-5704*x1*x5*x8^2+3036*x1*x6*x8*x9-5546*x1*x7^2*x9+920*x1*x8^2*x9+625*x2^4+1000*x2^3*x3+400*x2^2*x3^2+2050*x2^2*x4*x7+4450*x2^2*x4*x9-300*x2^2*x5*x7+1640*x2*x3*x4*x7+3560*x2*x3*x4*x9-240*x2*x3*x5*x7+7056*x3^2*x4^2+12600*x3^2*x4*x6+1848*x3^2*x4*x7+6781*x3^2*x6^2+3486*x3^2*x6*x7+1526*x3^2*x7^2+729*x3^2*x8^2+1890*x3*x4^2*x8+9912*x3*x4*x7*x9-1700*x3*x5*x6^2-1350*x3*x5*x6*x7-340*x3*x5*x6*x9+3224*x3*x5*x7*x8-270*x3*x5*x7*x9+8850*x3*x6*x7*x9+1782*x3*x6*x8^2-5712*x3*x6*x8*x9+1298*x3*x7^2*x9-5056*x3*x7*x8*x9+1225*x4^4+2310*x4^2*x6*x8+1681*x4^2*x7^2+7298*x4^2*x7*x9+7921*x4^2*x9^2-492*x4*x5*x7^2-1068*x4*x5*x7*x9+625*x5^2*x6^2+250*x5^2*x6*x9+36*x5^2*x7^2+3844*x5^2*x8^2+25*x5^2*x9^2+4200*x5*x6*x8*x9-1240*x5*x8^2*x9+840*x5*x8*x9^2+1089*x6^2*x8^2+3481*x7^2*x9^2+7156*x8^2*x9^2-1786*x1*x3*x7+5980*x1*x4*x8-3040*x1*x4*x9-7360*x1*x9^2+3192*x3^2*x4+2850*x3^2*x6+418*x3^2*x7-3380*x3*x4*x7-6528*x3*x6*x7-5184*x3*x7^2+2242*x3*x7*x9-4320*x3*x8*x9-5600*x4^2*x9-8060*x4*x5*x8+1300*x4*x8*x9+4800*x5*x6*x7+960*x5*x7*x9-5280*x6*x8*x9+16128*x7*x8*x9+361*x3^2+4225*x4^2+9216*x7^2+6400*x9^2, 2500*x2^3+3000*x2^2*x3+800*x2*x3^2+4100*x2*x4*x7+8900*x2*x4*x9-600*x2*x5*x7+1640*x3*x4*x7+3560*x3*x4*x9-240*x3*x5*x7, -7896*x1*x4*x7+1026*x1*x4*x8-7050*x1*x6*x7-1034*x1*x7^2-2392*x1*x7*x8+2484*x1*x8*x9+1000*x2^3+800*x2^2*x3+1640*x2*x4*x7+3560*x2*x4*x9-240*x2*x5*x7+14112*x3*x4^2+25200*x3*x4*x6+3696*x3*x4*x7+13562*x3*x6^2+6972*x3*x6*x7+3052*x3*x7^2+1458*x3*x8^2+1890*x4^2*x8+9912*x4*x7*x9-1700*x5*x6^2-1350*x5*x6*x7-340*x5*x6*x9+3224*x5*x7*x8-270*x5*x7*x9+8850*x6*x7*x9+1782*x6*x8^2-5712*x6*x8*x9+1298*x7^2*x9-5056*x7*x8*x9-1786*x1*x7+6384*x3*x4+5700*x3*x6+836*x3*x7-3380*x4*x7-6528*x6*x7-5184*x7^2+2242*x7*x9-4320*x8*x9+722*x3, 722*x1^2*x4+1748*x1^2*x9-7896*x1*x3*x7+1026*x1*x3*x8+3990*x1*x4^2+6440*x1*x4*x9+1254*x1*x6*x8+2050*x2^2*x7+4450*x2^2*x9+1640*x2*x3*x7+3560*x2*x3*x9+14112*x3^2*x4+12600*x3^2*x6+1848*x3^2*x7+3780*x3*x4*x8+9912*x3*x7*x9+4900*x4^3+4620*x4*x6*x8+3362*x4*x7^2+14596*x4*x7*x9+15842*x4*x9^2-492*x5*x7^2-1068*x5*x7*x9+5980*x1*x8-3040*x1*x9+3192*x3^2-3380*x3*x7-11200*x4*x9-8060*x5*x8+1300*x8*x9+8450*x4, -5704*x1*x8^2-300*x2^2*x7-240*x2*x3*x7-1700*x3*x6^2-1350*x3*x6*x7-340*x3*x6*x9+3224*x3*x7*x8-270*x3*x7*x9-492*x4*x7^2-1068*x4*x7*x9+1250*x5*x6^2+500*x5*x6*x9+72*x5*x7^2+7688*x5*x8^2+50*x5*x9^2+4200*x6*x8*x9-1240*x8^2*x9+840*x8*x9^2-8060*x4*x8+4800*x6*x7+960*x7*x9, -7050*x1*x3*x7+1254*x1*x4*x8+3036*x1*x8*x9+12600*x3^2*x4+13562*x3^2*x6+3486*x3^2*x7-3400*x3*x5*x6-1350*x3*x5*x7-340*x3*x5*x9+8850*x3*x7*x9+1782*x3*x8^2-5712*x3*x8*x9+2310*x4^2*x8+1250*x5^2*x6+250*x5^2*x9+4200*x5*x8*x9+2178*x6*x8^2+2850*x3^2-6528*x3*x7+4800*x5*x7-5280*x8*x9, 4418*x1^2*x7-7896*x1*x3*x4-7050*x1*x3*x6-2068*x1*x3*x7-2392*x1*x3*x8-11092*x1*x7*x9+2050*x2^2*x4-300*x2^2*x5+1640*x2*x3*x4-240*x2*x3*x5+1848*x3^2*x4+3486*x3^2*x6+3052*x3^2*x7+9912*x3*x4*x9-1350*x3*x5*x6+3224*x3*x5*x8-270*x3*x5*x9+8850*x3*x6*x9+2596*x3*x7*x9-5056*x3*x8*x9+3362*x4^2*x7+7298*x4^2*x9-984*x4*x5*x7-1068*x4*x5*x9+72*x5^2*x7+6962*x7*x9^2-1786*x1*x3+418*x3^2-3380*x3*x4-6528*x3*x6-10368*x3*x7+2242*x3*x9+4800*x5*x6+960*x5*x9+16128*x8*x9+18432*x7, 4232*x1^2*x8+1026*x1*x3*x4-2392*x1*x3*x7+2484*x1*x3*x9+1254*x1*x4*x6-11408*x1*x5*x8+3036*x1*x6*x9+1840*x1*x8*x9+1458*x3^2*x8+1890*x3*x4^2+3224*x3*x5*x7+3564*x3*x6*x8-5712*x3*x6*x9-5056*x3*x7*x9+2310*x4^2*x6+7688*x5^2*x8+4200*x5*x6*x9-2480*x5*x8*x9+840*x5*x9^2+2178*x6^2*x8+14312*x8*x9^2+5980*x1*x4-4320*x3*x9-8060*x4*x5+1300*x4*x9-5280*x6*x9+16128*x7*x9, 1748*x1^2*x4+4232*x1^2*x9+2484*x1*x3*x8+3220*x1*x4^2+3036*x1*x6*x8-5546*x1*x7^2+920*x1*x8^2+4450*x2^2*x4+3560*x2*x3*x4+9912*x3*x4*x7-340*x3*x5*x6-270*x3*x5*x7+8850*x3*x6*x7-5712*x3*x6*x8+1298*x3*x7^2-5056*x3*x7*x8+7298*x4^2*x7+15842*x4^2*x9-1068*x4*x5*x7+250*x5^2*x6+50*x5^2*x9+4200*x5*x6*x8-1240*x5*x8^2+1680*x5*x8*x9+6962*x7^2*x9+14312*x8^2*x9-3040*x1*x4-14720*x1*x9+2242*x3*x7-4320*x3*x8-5600*x4^2+1300*x4*x8+960*x5*x7-5280*x6*x8+16128*x7*x8+12800*x9];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
export  transitionFluxReport

"""
```
transitionFluxReport(
    mpSim::MPsim,
    timeGrid::Vector{Float64},
    transitions::TransitionType... )
```
This function generates reports for the flux of the unique transitions in `transitions` of the manpower simulation `mpSim`. The reports are generated on the grid of time points `timeGrid`.

The transitions are given either as:
* a `String`: the name of the transition, where `"attrition"` (personnel attrition) and `"retirement"` (default retirement scheme) are valid names,
* a `Tuple{String, String}`: the source and target nodes of the transition, where `""`, `"out"`, and `"external"` refer to "not in the system", or
* a `Tuple{String, String, String}`: the name and the source and target nodes of the transition.

If multiple transitions have the same name, the fluxes will be summed together.

This function will issue a warning and not generate any report in the following two cases:
1. There are no positive time points in the time grid;
2. None of the entered transitions exist in the simulation.

The function returns a `DataFrame` holding the flux report. The first column are the time points, and each other column corresponds to one of the valid transitions. In case the function issues a warning, its return value will be an empty `DataFrame`.
"""
function transitionFluxReport( mpSim::MPsim, timeGrid::Vector{Float64},
    transitions::TransitionType... )::DataFrame

    result = DataFrame()

    if now( mpSim ) == 0
        @warn "Simulation hasn't started yet, can't make report."
        return result
    end  # if now( mpSim ) == 0

    timeGrid = timeGrid[0.0 .<= timeGrid .<= now( mpSim )]
    timeGrid = unique( sort( timeGrid, rev = true ) )

    if isempty( timeGrid )
        @warn "No valid time points in time grid, cannot generate report."
        return result
    end  # if isempty( timeGrid )

    if timeGrid[end] > 0.0
        push!( timeGrid, 0.0 )
    end  # if timeGrid[end] > 0.0

    reverse!( timeGrid )

    transitions = filter( transition -> validateTransition( mpSim, transition ),
        collect( transitions ) )

    if isempty( transitions )
        @warn "No valid transitions in list, cannot generate report."
        return result
    end  # if isempty( transitions )

    unique!( transitions )
    results = zeros( Int, length( timeGrid ), length( transitions ) )

    for ii in eachindex( transitions )
        results[:, ii] = createTransitionFluxReport( mpSim, timeGrid,
            transitions[ii] )
    end  # for transition in transitions

    transitionNames = vcat( "timeStart", "timePoint",
        generateTransitionName.( transitions ) )

    return DataFrame( hcat( vcat( 0, timeGrid[1:(end - 1)] ), timeGrid,
        results ), Symbol.( transitionNames ) )

end  # transitionFluxReport( mpSim, timeGrid, transitions )

"""
```
transitionFluxReport(
    mpSim::MPsim,
    timeRes::Real,
    transitions::TransitionType... )
```
This function generates reports for the flux of the unique transitions in `transitions` of the manpower simulation `mpSim`. The reports are generated on a grid of time points with resolution `timeRes`.

The transitions are given either as:
* a `String`: the name of the transition, where `"attrition"` (personnel attrition) and `"retirement"` (default retirement scheme) are valid names,
* a `Tuple{String, String}`: the source and target nodes of the transition, where `""`, `"out"`, and `"external"` refer to "not in the system", or
* a `Tuple{String, String, String}`: the name and the source and target nodes of the transition.

If multiple transitions have the same name, the fluxes will be summed together.

This function will issue a warning and not generate any report in the following two cases:
1. The resolution of the time grid is ⩽ 0;
2. None of the entered transitions exist in the simulation.

The function returns a `DataFrame` holding the flux report. The first column are the time points, and each other column corresponds to one of the valid transitions. In case the function issues a warning, its return value will be an empty `DataFrame`.
"""
transitionFluxReport( mpSim::MPsim, timeRes::Real,
    transitions::TransitionType... )::DataFrame = transitionFluxReport( mpSim,
    generateTimeGrid( mpSim, timeRes ), transitions... )

include( joinpath( repPrivPath, "transitionflux.jl" ) )


text:
export testAbstractMeganetElement

"""
simple tests that every AbstractMeganetElement must pass
"""
function testAbstractMeganetElement(L::AbstractMeganetElement{T};out::Bool=false,nex=8) where {T <: Number}

    @testset "features immutable" begin
        theta = initTheta(L)
        theta .+= .1 # To test if Y changes for affineScalingLayer
        Y     = randn(T,nFeatIn(L),nex)
        Yo    = copy(Y)
        # THE FOLLOWING HACK (copy(Y)) EFFECTIVELY CANCELS THIS TEST...COMMENTING IT OUT
        Zd,Z,tmp  = apply(L,theta,copy(Y),[],true)
#         @test norm(Y-Yo)/norm(Yo) < 1e4*eps(T)

        dY    = randn(T,nFeatIn(L),nex)
        Z1    =  JYmv(L,dY,theta,Y,tmp)[2]
        @test norm(Y-Yo)/norm(Yo) < 1e4*eps(T)
        Z12    =  JYmv(L,dY,theta,Y,tmp)[2]
        @test norm(Z1-Z12)/(norm(Z1)+(norm(Z1)==0)) < 1e4*eps(T)

        dth    = randn(T,nTheta(L))
        Z1    =  Jthetamv(L,dth,theta,Y,tmp)[2]
        @test norm(Y-Yo)/norm(Yo) < 1e4*eps(T)
        Z12    =  Jthetamv(L,dth,theta,Y,tmp)[2]
        @test norm(Z1-Z12)/(norm(Z1)+(norm(Z1)==0)) < 1e4*eps(T)

        Zout = randn(T,nFeatOut(L),nex)
        Z1 =  JthetaTmv(L,copy(Zout),(T)[],theta,Y,tmp)
        @test norm(Y-Yo)/norm(Yo) < 1e4*eps(T)
        Z12 =  JthetaTmv(L,copy(Zout),(T)[],theta,Y,tmp)
        @test norm(Z1-Z12)/(norm(Z1)+(norm(Z1)==0)) < 1e4*eps(T)

        Zout = randn(T,nFeatOut(L),nex)
        Z1 =  JYTmv(L,copy(Zout),(T)[],theta,Y,tmp)
        @test norm(Y-Yo)/norm(Yo) < 1e4*eps(T)
        Z12 =  JYTmv(L,copy(Zout),(T)[],theta,Y,tmp)
        @test norm(Z1-Z12)/(norm(Z1)+(norm(Z1)==0)) < 1e4*eps(T)

    end

    @testset "apply without derivatives" begin
        theta = initTheta(L)
        Y     = randn(T,nFeatIn(L),nex)
        Z     = apply(L,theta,Y,[],false)
    end

    @testset "derivative w.r.t Y test" begin

        theta = initTheta(L)
       Y     = randn(T,nFeatIn(L),nex)

       function testFun(x,v=[])
           if !(isempty(v))
               Z = apply(L,theta,copy(x),[],true)
               return Z[2], reshape(JYmv(L,v,theta,x,Z[3])[2],size(Z[2]))
           else
               return apply(L,theta,copy(x),[])[2]
           end
       end
       chkDer, = checkDerivative(testFun,copy(Y),out=out)
       @test chkDer
    end

    @testset "adjoint w.r.t. Y test" begin
    theta = initTheta(L)
       Y     = randn(T,nFeatIn(L),nex)
       dY    = randn(T,nFeatIn(L),nex)
       Z     = randn(T,nFeatOut(L),nex)

       tmp = apply(L,theta,copy(Y),[],true)
       Z1 =  JYmv(L,copy(dY),theta,copy(Y),tmp[3])[2]
       Z2 =  JYTmv(L,copy(Z),(T)[],theta,copy(Y),tmp[3])

       t1 = dot(Z1,Z)
       t2 = dot(Z2,dY)
       @test abs(t1-t2)/abs(t2) < 1e4*eps(T)

    end

    if nTheta(L)>0
        @testset "derivative w.r.t theta test" begin
        theta = initTheta(L)
          Y     = randn(T,nFeatIn(L),nex)

          function testFunTh(x,v=[])
              if !(isempty(v))
                  Z = apply(L,x,copy(Y),[],true)
                  return Z[2], reshape(Jthetamv(L,v,x,copy(Y),Z[3])[2],size(Z[2]))
              else
                  return apply(L,x,copy(Y),[])[2]
              end
          end
          chkDer, = checkDerivative(testFunTh,copy(theta),out=out)
          @test chkDer
        end

        @testset "adjoint w.r.t. theta test" begin
        theta = initTheta(L)
          Y     = randn(T,nFeatIn(L),nex)
          dtheta    = randn(T,nTheta(L))
          Z     = randn(T,nFeatOut(L),nex)

          tmp = apply(L,theta,copy(Y),[],true)
          Z1 =  Jthetamv(L,copy(dtheta),copy(theta),copy(Y),copy(tmp[3]))[2]
          Z2 =  JthetaTmv(L,copy(Z),(T)[],theta,copy(Y),tmp[3])

          t1 = dot(Z1,Z)
          t2 = dot(Z2,dtheta)
          if out
              println("t1=$t1 \t t2=$t2 \t err = $(t1-t2)")
          end
          @test abs(t1-t2)/abs(t2) < 1e4*eps(T)
        end
    end
end


text:
using CUDA
using KernelAbstractions

NMax = 370;

an = (zeros(Complex{Float32},NMax))
an .+= 1.3-0.01im
bn = (zeros(Complex{Float32},NMax))
bn .+= 1.2-0.1im
wignerA = (zeros(Complex{Float32},2NMax+1,NMax,NMax)); fill!(wignerA,1.0-0im);
wignerB = (zeros(Complex{Float32},2NMax+1,NMax,NMax)); fill!(wignerB,1.0-0im);
Sνν = zeros(2NMax)

@kernel function compute_Sl_νν!(@Const(wignerA),@Const(wignerB), an, bn, Sνν, @Const(lMax))
    # Get length of an:
    san = size(an)[1];
    # Indices over n and m
    n, m = @index(Global, NTuple)

    # Outer loop over l
    for l = 1:lMax
        if max(l-n,n) <= m <= min(san,n+l)
            #println((an[n]' + bn[n]') * (an[m] + bn[m]) * wignerA[l,n,m]^2)
            Sνν[l] += real((an[n]' + bn[n]') * (an[m] + bn[m]) * wignerA[l,n,m]^2)
            Sνν[l] += 1/2 *  real(abs2(an[n] + bn[n]) * wignerA[l,n,m]^2)
        end
    end

end

kernel! = compute_Sl_νν!(CPU(),4);
kernelGPU! = compute_Sl_νν!(CUDADevice());
@time event = kernel!(wignerA, wignerB, an, bn, Sνν, 2NMax, ndrange=(NMax,NMax)); wait(event)

wignerA_CU = CuArray(wignerA);
wignerB_CU = CuArray(wignerB);
an_CU      = CuArray(an);
bn_CU      = CuArray(bn);
Sνν_CU     = CuArray(Sνν);
@time event = kernelGPU!(wignerA_CU, wignerB_CU, an_CU, bn_CU, Sνν_CU, 2NMax, ndrange=(NMax,NMax)); wait(event)