text:
# # data file for gnash1.mod
# # Original AMPL coding by Sven Leyffer, University of Dundee
#
# # An MPEC from F. Facchinei, H. Jiang and L. Qi, A smoothing method for
# # mathematical programs with equilibrium constraints, Universita di Roma
# # Technical report, 03.96. Problem number 8
#
# param	:	c	K	b	:=
# 	1	10	5	1.2
# 	2	 8	5	1.1
# 	3	 6	5	1.0
# 	4	 4	5	0.9
# 	5	 2	5	0.8	;
#
# param L := 20;
# param g := 1.7;
#
# let x := 75;

c = [10, 8, 6, 4, 2]
K = [5, 5, 5, 5, 5]
b = [1.2, 1.1, 1.0, 0.9, 0.8]
L = 20
g = 1.7
# x = 75

# # gnash1m.mod	QQR2-MN-10-8
# # Original AMPL coding by Sven Leyffer, University of Dundee
#
# # Formulation with mixed complementarity of ...
# #
# # An MPEC from F. Facchinei, H. Jiang and L. Qi, A smoothing method for
# # mathematical programs with equilibrium constraints, Universita di Roma
# # Technical report, 03.96. Problem number 8
# #
# # Arises from Gournot Nash equilibrium , 10 instances are available
# # (see gnash1i.dat, i=0,1,...,9)
#
# # Number of variables:   10
# # Number of constraints:  8
#
# # ... parameters for each firm/company
# param c{1..5};			# c_i
# param K{1..5};			# K_i
# param b{1..5};			# \beta_i
#
# # ... parameters for each problem instance
# param L;			# L
# param g;			# \gamma
#
# # ... computed constants
# param gg := 5000^(1/g);
#
# var x >= 0, <= L;
# var y{1..4};
# var l{1..4};   		        # Multipliers
# var Q = x+y[1]+y[2]+y[3]+y[4];	# defined variable Q
#
# minimize f: c[1]*x + b[1]/(b[1]+1)*K[1]^(-1/b[1])*x^((1+b[1])/b[1])
# 		- x*( gg*Q^(-1/g) );
#
# subject to
#
#    F1: 0 = ( c[2] + K[2]^(-1/b[2])*y[1] ) - ( gg*Q^(-1/g) )
# 				- y[1]*( -1/g*gg*Q^(-1-1/g) ) - l[1];
#    F2: 0 = ( c[3] + K[3]^(-1/b[3])*y[2] ) - ( gg*Q^(-1/g) )
# 				- y[2]*( -1/g*gg*Q^(-1-1/g) ) - l[2];
#    F3: 0 = ( c[4] + K[4]^(-1/b[4])*y[3] ) - ( gg*Q^(-1/g) )
# 				- y[3]*( -1/g*gg*Q^(-1-1/g) ) - l[3];
#    F4: 0 = ( c[5] + K[5]^(-1/b[5])*y[4] ) - ( gg*Q^(-1/g) )
# 				- y[4]*( -1/g*gg*Q^(-1-1/g) ) - l[4];
#
#    g1: 0 <= y[1] <= L  complements   l[1];
#    g3: 0 <= y[2] <= L  complements   l[2];
#    g5: 0 <= y[3] <= L  complements   l[3];
#    g7: 0 <= y[4] <= L  complements   l[4];

using JuMP, Complementarity
using Ipopt
using Test

@testset "gnash1m.jl" begin

    gg = 5000^(1/g)

    gnash1m = Model(with_optimizer(Ipopt.Optimizer))

    @variable(gnash1m, 0 <= x <= L)
    @variable(gnash1m, y[1:4])
    @variable(gnash1m, l[1:4])
    @variable(gnash1m, Q >= 0)
    @constraint(gnash1m, Q == x+y[1]+y[2]+y[3]+y[4])
    @NLobjective(gnash1m, Min, c[1]*x + b[1]/(b[1]+1)*K[1]^(-1/b[1])*x^((1+b[1])/b[1])
     		             - x*( gg*Q^(-1/g) ) )

    @NLconstraint(gnash1m, cnstr[i=1:4], 0 == ( c[i+1] + K[i+1]^(-1/b[i+1])*y[i] ) - ( gg*Q^(-1/g) )
                                      - y[i]*( -1/g*gg*Q^(-1-1/g) ) - l[i] )

    for i in 1:4
        @complements(gnash1m, l[i], 0 <= y[i] <= L, smooth)
    end

    # Initial solutions to help reaching the optimality
    set_start_value.(y, [11, 19, 20, 20])
    set_start_value.(l, [0, 0, -3, -5])
    set_start_value(x, 5)
    set_start_value(Q, 70)

    JuMP.optimize!(gnash1m)

    @show JuMP.value.(y)
    @show JuMP.value.(l)
    @show JuMP.value(x)
    @show JuMP.value(Q)

    # JuMP.value.(y) = [11.3327, 19.5003, 20.0, 20.0]
    # JuMP.value.(l) = [-1.35686e-10, -9.7489e-9, -2.51077, -5.18082]
    # JuMP.value(x) = 6.369341692612341
    # JuMP.value(Q) = 77.20235592852553
    # JuMP.objective_value(gnash1m) = -6.116708234438121

    @show JuMP.objective_value(gnash1m)
    @test isapprox(JuMP.objective_value(gnash1m), -6.11671, atol=1e-4)
    @test isapprox( JuMP.value.(l)' * (L .- JuMP.value.(y)), 0, atol=1e-5)
end


text:
using SafeTestsets

module DirectFixture
using Distributions: ContinuousUnivariateDistribution
const CTDist = ContinuousUnivariateDistribution
# Use import so that we can extend the method.
import Fleck: hazards
using Random: AbstractRNG

function hazards(callback::Function, distlist::Array{T, 1}, rng::AbstractRNG) where T <: CTDist
    for dist_idx in 1:length(distlist)
        callback(dist_idx, distlist[dist_idx], true)
    end
end
export CTDist
end

@safetestset markov_direct_initial = "MarkovDirect initial" begin
    using Fleck: MarkovDirect, next
    using Random: MersenneTwister
    using Distributions: Exponential
    using ..DirectFixture
    md = MarkovDirect()
    distributions = fill(Exponential(1.5), 10)
    rng = MersenneTwister(90422342)
    current = 0.0
    when, which = next(md, distributions, current, rng)
    @test when > 0
    @test 1 <= which
    @test which <= 10
end

@safetestset markov_direct_empty = "MarkovDirect empty hazard" begin
    using Fleck: MarkovDirect, next
    using Random: MersenneTwister
    using Distributions: Exponential
    using ..DirectFixture
    md = MarkovDirect()
    distributions = CTDist[]
    rng = MersenneTwister(90497979)
    current = 0.0
    when, which = next(md, distributions, current, rng)
    @test isinf(when)
    @test which === nothing
end

@safetestset markov_direct_prob = "MarkovDirect probabilities correct" begin
    using Fleck: MarkovDirect, next
    using Random: MersenneTwister
    using Distributions: Exponential
    using HypothesisTests: BinomialTest, confint
    using ..DirectFixture

    md = MarkovDirect()
    distributions = vcat(
        fill(Exponential(1.5), 10),
        fill(Exponential(1), 10)
    )
    rng = MersenneTwister(90422342)
    hilo = zeros(Int, 2)
    curtime = 0.0
    for i in 1:10000
        when, which = next(md, distributions, curtime, rng)
        hilo[(which - 1) ÷ 10 + 1] += 1
    end
    ci = confint(BinomialTest(hilo[1], sum(hilo), 3 / 5))
    @test ci[1] < 3 / 5 < ci[2]
end

@safetestset direct_call_prob = "DirectCall probabilities" begin
    using Fleck: DirectCall, set_clock!, next
    using Random: MersenneTwister
    using Distributions: Exponential

    dc = DirectCall(Int)
    rng = MersenneTwister(90422342)
    propensities = [0.3, 0.2, 0.7, 0.001, 0.25]
    for (i, p) in enumerate(propensities)
        set_clock!(dc, i, Exponential(p), :Enabled, rng)
    end
    when, which = next(dc, 100.0, rng)
    @test when > 100
    @test 1 <= which
    @test which <= length(propensities)
end

ext: jl


text:
using EventSimulation
using Distributions
using Statistics
using DataFrames
using PyPlot

mutable struct MMQueue{T<:Distribution} <: AbstractState
    sp::T # service rate
    qlen::Int   # queue length
    busy::Bool  # is server busy?
    arrival_times::Vector{Float64} # when customers arrived
    departure_times::Vector{Float64} # when customers left
end

# customer arrival event
function arrival(s)
    push!(s.state.arrival_times, s.now)
    if s.state.busy # server is working?
        s.state.qlen += 1
    else
        s.state.busy = true
        register!(s, leave, rand(s.state.sp))
    end
end

# customer leave the system event
function leave(s)
    push!(s.state.departure_times, s.now)
    if s.state.qlen > 0 # any customers waiting?
        s.state.qlen -= 1
        register!(s, leave, rand(s.state.sp))
    else
        s.state.busy = false
    end
end

function run(ar::Distribution, sr::Distribution)
    q = MMQueue(sr, 0, false, Float64[], Float64[])
    s = Scheduler(q, Float64)
    repeat_register!(s, arrival, x -> rand(ar))
    go!(s, 1_000_000)
    mean(d-a for (a,d) in zip(s.state.arrival_times, s.state.departure_times))
end

results = DataFrame()
sr = 1.0
sre = Exponential(1/sr)
sru = Uniform(0, 2/sr)
for ar in 0.1:0.1:0.9
    are = Exponential(1/ar)
    aru = Uniform(0, 2/ar)
    push!(results, (ar=ar,
                    sojourn_ee=run(are, sre),
                    sojourn_ue=run(aru, sre),
                    sojourn_eu=run(are, sru),
                    sojourn_uu=run(aru, sru),
                    theoretical=1/(1-ar/sr)/sr))
end

plot(results.ar, Matrix(results[:, 2:end]))
legend(names(results)[2:end])


text:
#######################
# Filtering distribution
# This does not save trajectories, but saves weights and particles at each time step
# Thus each corresponds to the filtering distribution.
export particle_filter, particle_smoother

abstract type ParticleContainer <: AbstractArray{Float64,1}
end

struct FilteringParticleContainer <: ParticleContainer
    X # Particles
    w # weights
    ℓ # Estimate of the log likelihood
end

Base.size(ps::FilteringParticleContainer) = (size(ps.X,2),)
Base.length(ps::FilteringParticleContainer) = prod(size(ps))
Base.getindex(ps::FilteringParticleContainer,i::Int) = ps.X[:,i]
Base.setindex!(ps::FilteringParticleContainer,v,i::Int) = ps.X[:,i] = v
Base.IndexStyle(ps::FilteringParticleContainer) = IndexLinear()

StatsBase.loglikelihood(ps::FilteringParticleContainer) = ps.ℓ[end]

resample!(ps::FilteringParticleContainer,m,t,θ,y,u,sample_function=resample_stratified) = resample!(FilteringParticleContainer,ps.X,ps.w,ps.w[:,t],m,t,θ,y,u,sample_function)

function resample!(::Type{FilteringParticleContainer},X,w,ξs,m,t,θ,y,u,sample_function)
    for i in eachindex(ξs)
        ξs[i] += weight_function(m,X[:,i,t],y,u,t,θ)
    end
    Σ = sum(exp,ξs)
    idx = sample_function(exp.(ξs)./Σ,length(ξs))

    X[:,:,t] .= X[:,idx,t]
    w[:,t] .= w[idx,t]
    ξs[idx],Σ
end

propagate!(ps::FilteringParticleContainer,ξ,Σ,m,t,θ,y,u) = propagate!(FilteringParticleContainer,ps.X,ps.w,ps.ℓ,ξ,Σ,m,t,θ,y,u)

function propagate!(::Type{FilteringParticleContainer},X,w,ℓ,ξ,Σ,m,t,θ,y,u)
    for i in eachindex(ξ)
        #dq = proposal_distribution(m,X[:,i,t],y,u,t,θ)
        #rand!(dq,view(X,:,i,t+1))
        proposal_rand!(view(X,:,i,t+1),m,X[:,i,t],y,u,t,θ)

        w[i,t+1] = w[i,t]
        w[i,t+1] += transition_logpdf(m,X[:,i,t+1],X[:,i,t],u,t,θ)
        w[i,t+1] += observation_logpdf(m,y,X[:,i,t+1],u,t,θ)
        #w[i] *= pdf(f(X[:,i],t,θ),Xn)
        #w[i] *= pdf(g(Xn,t,θ),y)
        w[i,t+1] -= ξ[i]
        w[i,t+1] -= proposal_logpdf(m,X[:,i,t+1],X[:,i,t],y,u,t,θ)

    end
    ℓn = log(sum(exp,w[:,t+1]))
    ℓ[t+1] = ℓ[t] + ℓn - log(length(ξ)) + log(Σ)
    w[:,t+1] .-= ℓn
end

function step(ps::FilteringParticleContainer,m::StateSpaceModel,t,θ,y,u;sample_function=resample_stratified)

    # Resample

    #rs = Resampler(ps,m,y,t,θ,sample_function)

    ξs,Σ = resample!(ps,m,t,θ,y,u,sample_function)

    # Propagate

    #ps2 = ScoreParticleContainer([propose(t,x,w,ξ,α,ps.S,y,θ,m.f,m.g,m.q,m.dp,ps.dg,λ) for (x,w,ξ,α) in rs],ps.ℓ,ps.dg,sum(rs.ξ),λ)

    propagate!(ps,ξs,Σ,m,t,θ,y,u)

    #normalize_weights!(ps2)

    ps
end

function particle_filter(m::StateSpaceModel,θ,Y,N,u=fill([0.0],length(Y));sample_function=resample_stratified,verbose=false)

    T = length(Y)

    if verbose
        p = Progress(T,1)
    end

    ps = initialize_particles(FilteringParticleContainer,m,θ,Y,N)

    for t in 1:T
        ps = step(ps,m,t,θ,Y[t],u[t],sample_function=sample_function)
        if verbose
            next!(p)
        end
    end

    ps
end

function initialize_particles(::Type{FilteringParticleContainer},m::StateSpaceModel,θ,Y,N)

    # We should probably think of μ0 as a distribution for p(x_0)
    # not one for p(x_1)
    # Initialization gives X0, w = 1/N
    #x0 = rand(μ0(θ),N)
    T = length(Y)

    x0 = initial_rand(m,θ,N)
    X0 = zeros(size(x0,1),N,T+1)
    X0[:,:,1] = x0[:,:]

    w = log.(ones(N,T+1)/N)
    ℓ = zeros(T+1)
    #w = [pdf(g(X0[:,end,i],1,θ),y)/N for i in 1:N]
    #ℓ = sum(w)
    #w ./= ℓ

    FilteringParticleContainer(X0,w,ℓ)
end

function particle_smoother(m::StateSpaceModel,θ,Y,pf::FilteringParticleContainer,B,u=fill([0.0],length(Y));sample_function=resample_stratified)
    D = size(pf.X,1)
    N = size(pf.X,2)
    T = size(pf.X,3)-1

    x = zeros(D,T+1,B)
    for i in 1:B
        x[:,end,i] = pf.X[:,sample_function(exp.(pf.w[:,end])./sum(exp,pf.w[:,end]),1),end]
    end

    for t in T-1:-1:0
        wsmooth = pf.w[:,t+1]
        for j in 1:B
            for i in 1:N
            # This is quite slow
                wsmooth[i] += transition_logpdf(m,x[:,t+2,j],pf.X[:,i,t+1],u[t+1],t,θ)
            end

            x[:,t+1,j] = pf.X[:,sample_function(exp.(wsmooth)./sum(exp,wsmooth),1),t+1]
        end
    end
    x
end


text:
n = BigInt(28433)
for i = 1:7830457
  global n
  n *= 2
  n %= BigInt(1e10)
end
n += 1
println(n)

max_stars_repo_name: pdbartlett/misc-stuff

max_issues_repo_name: pdbartlett/misc-stuff

max_forks_repo_name: pdbartlett/misc-stuff


text:
using JSON,JLD,Knet
include("loss.jl")

if !isdefined(Main,:atype)
    global atype = KnetArray{Float32}
end

init(o...)=atype(xavier(Float32,o...))
bias(o...)=atype(zeros(Float32,o...))
elu(x) = relu.(x) .+ (exp.(min.(0,x)) - 1)

function load_resnet(atype;stage=3)
    w,m,meta = ResNetLib.resnet101init(;trained=true,stage=stage)
    global avgimg = meta["normalization"]["averageImage"]
    global descriptions = meta["classes"]["description"]
    return w,m,meta,avgimg
end

function prepocess_image(w,m,imgurl,avgimg;stage=3)
    img = imgdata(imgurl, avgimg)
    return ResNetLib.resnet101(w,m,atype(img);stage=stage);
end

function postpocess_kb(w,x;train=false,pdrop=0.0)
    if train
         x  = dropout(x,pdrop)
    end
    x  = elu(conv4(w[1],x;padding=1,stride=1,mode=1) .+ w[2]) #relu -> elu
    if train
         x  = dropout(x,pdrop)
    end
    x  = elu(conv4(w[3],x;padding=1,stride=1,mode=1) .+ w[4]) #relu -> elu
    h,w,c,b = size(x)
    x  = reshape(x,h*w,c,b)
    return permutedims(x,(2,3,1))
end

function init_postkb(d)
    w = Any[];
    push!(w,init(3,3,1024,d))
    push!(w,bias(1,1,d,1))
    push!(w,init(3,3,d,d))
    push!(w,bias(1,1,d,1))
    return w
end

function process_question(w,r,words,batchSizes;train=false,qdrop=0.0,embdrop=0.0)
    wordemb      = w[1][:,words]

    if train
         wordemb = dropout(wordemb,embdrop)
    end

    B = batchSizes[1]
    eqbatches = all(batchSizes.==B)

    if eqbatches
        wordemb      = reshape(wordemb,size(wordemb,1),1,size(wordemb,2))
        y,hyout,_,rs = rnnforw(r,w[2],wordemb;hy=true,cy=false)
    else
        y,hyout,_,rs = rnnforw(r,w[2],wordemb;batchSizes=batchSizes,hy=true,cy=false)
    end

    q            = vcat(hyout[:,:,1],hyout[:,:,2])

    if train
           q  = dropout(q,qdrop)
    end

    if !eqbatches
        indices      = batchSizes2indices(batchSizes)
        lngths       = length.(indices)
        Tmax         = maximum(lngths)
        td,B         = size(q)
        d            = div(td,2)
        cw           = Any[];

        for i=1:length(indices)
            y1 = y[:,indices[i]]
            df = Tmax-lngths[i]
            if df > 0
                cpad = zeros(Float32,2d,df)
                kpad = atype(cpad) ## look at similar
                ypad = hcat(y1,kpad)
                push!(cw,ypad)
            else
                push!(cw,y1)
            end
        end
        cws_2d =  reshape(vcat(cw...),2d,B*Tmax)
    else
        B      = batchSizes[1]
        Tmax   = length(batchSizes)
        d      = div(size(q,1),2)
        cws_2d = reshape(y,2d,B*Tmax)
    end

    cws_3d =  reshape(w[3]*cws_2d .+ w[4],(d,B,Tmax))

    return q,cws_3d;
end

function batchSizes2indices(batchSizes)
    B = batchSizes[1]
    indices = Any[]
    for i=1:B
        ind = i.+cumsum(filter(x->(x>=i),batchSizes)[1:end-1])
        push!(indices,append!(Int[i],ind))
    end
    return indices
end

function init_rnn(inputSize,embedSize,hiddenSize;outputSize=nothing)
    w = Any[]
    wembed = atype(rand(Float32,embedSize,inputSize))
    push!(w,wembed)
    r,wrnn = rnninit(embedSize,hiddenSize; bidirectional=true, binit=zeros) #change Knet/src/rnn.jl
    #setfgbias!(r,wrnn)
    push!(w,wrnn)
    push!(w,init(hiddenSize,2hiddenSize),bias(hiddenSize,1))
    return w,r
end

function setfgbias!(r,w)
    rnnparam(r,w,1,2,2)[:] = 0.5
    rnnparam(r,w,1,6,2)[:] = 0.5
    rnnparam(r,w,2,2,2)[:] = 0.5
    rnnparam(r,w,2,6,2)[:] = 0.5
end

function control_unit(w,ci₋1,qi,cws,pad;train=false,tap=nothing)
    #cws       : d x B x T
    d,B,T      = size(cws)
    #qi,ci     : d x B
    cqi        = reshape(w[1] * vcat(ci₋1,qi) .+ w[2],(d,B,1)) # eq c1
    #cqi       : d x B x 1
    cvis       = reshape(cqi .* cws,(d,B*T)) #eq c2.1.1
    #cvis      : d x BT
    cvis_2d    = reshape(w[3] * cvis .+ w[4],(B,T)) #eq c2.1.2
    #cvis_2d   : B x T
    if pad != nothing
        cvi    = reshape(softmax(cvis_2d .- pad,2),(1,B,T)) #eq c2.2
    else
        cvi    = reshape(softmax(cvis_2d,2),(1,B,T)) #eq c2.2
    end
    tap!=nothing && get!(tap,"w_attn_$(tap["cnt"])",Array(reshape(cvi,B,T)))
    #cvi       : 1 x B x T
    ci         = reshape(sum(cvi.*cws,3),(d,B)) #eq c2.3
end

function init_control(d)
    w = Any[]
    push!(w,init(d,2d))
    push!(w,bias(d,1))
    push!(w,init(1,d))
    push!(w,bias(1,1))
    return w
end

function read_unit(w,mi₋1,ci,KBhw,KBhw′′;train=false,mdrop=0.0,attdrop=0.0,tap=nothing)
    d,B,N          = size(KBhw)
    d,BN           = size(KBhw′′)
    #KBhw'     : d x B x N  := w[3] * khw + w[4]
    #KBhw''    : d x BN     := reshape(w[5] * khw + w[7]),(d,BN))
    if train
        mi₋1       = dropout(mi₋1,mdrop)
    end

    mi_3d           = reshape(w[1]*mi₋1 .+ w[2],(d,B,1)) #eq r1.1
    #mi_3d     : d x B x 1
    ImKB            = reshape(mi_3d .* KBhw,(d,BN)) # eq r1.2
    #ImKB      : d x BN
    ImKB′           = reshape(w[5] * ImKB .+ KBhw′′,(d,B,N)) #eq r2
    #ImKB'     : d x B x N
    ci_3d           = reshape(ci,(d,B,1))
    #ci_3d     : d x B x 1
    IcmKB_pre       = reshape(ci_3d .* ImKB′,(d,BN)) #eq r3.1.1
    #IcmKB_pre : d x BN
    if train
         IcmKB_pre = dropout(IcmKB_pre,attdrop)#dropout(IcmKB_pre,0.15)
    end
    IcmKB           = reshape(w[6] * IcmKB_pre  .+ w[7],(B,N)) #eq r3.1.2
    #IcmKB     : B x N
    mvi             = reshape(softmax(IcmKB,2),(1,B,N)) #eq r3.2
    #mvi       : 1 x B x N
    tap!=nothing && get!(tap,"KB_attn_$(tap["cnt"])",Array(reshape(mvi,B,N)))
    mnew            = reshape(sum(mvi.*KBhw,3),(d,B)) #eq r3.3
end

function init_read(d)
    w = Any[]
    push!(w,init(d,d))
    push!(w,bias(d,1))
    push!(w,init(d,d))
    push!(w,bias(d,1))
    push!(w,init(d,d))
    push!(w,init(1,d))
    push!(w,bias(1,1))
    return w
end

function write_unit(w,m_new,mi₋1,mj,ci,cj;train=false,selfattn=true,gating=true)
    d,B        = size(m_new)
    #mnew      : d x B
    d,BT       = size(mj)
    T          = div(BT,B)
    #mj        : d x BT
    mi         =  w[1] * vcat(m_new,mi₋1) .+ w[2] #eq w.1 #no such layer in code delete w6-w7
    #mi       : d x B
    !selfattn && return mi
    #iproj     =  w[6] * ci .+ w[7] !!!delete w[6]
    ci_3d      = reshape(ci,d,B,1)
    #ci_3d     : d x B x 1
    cj_3d      = reshape(cj,d,B,T)
    #cj_3d     : d x B x T
    sap        = reshape(ci_3d.*cj_3d,(d,BT)) #eq w2.1.1
    #sap       : d x BT
    sa         = reshape(w[3] * sap .+ w[4],(B,T)) #eq w2.1.2
    #sa        : B x T
    sa′        = reshape(softmax(sa,2),(1,B,T)) #eq w2.1.3
    #sa'       : 1 x B x T
    mj_3d      = reshape(mj,(d,B,T))
    #mj_3d     : d x B x T
    mi_sa      = reshape(sum(sa′ .* mj_3d ,3),(d,B)) #eq w2.2
    #m_sa      : d x B
    mi′′       = w[5] * mi_sa .+ w[6] .+ mi   #eq w2.3
    #mi′′      : d x B

    !gating && return mi′′

    σci′       = sigm.(w[7] * ci .+ w[8])  #eq w3.1
    #σci′      : 1 x B
    mi′′′      = (σci′ .* mi₋1) .+  ((1.-σci′) .* mi′′) #eq w3.2
end

function init_write(d;selfattn=true,gating=true)
    w = Any[]
    push!(w,init(d,2d))
    push!(w,bias(d,1))
    !selfattn && return w
    push!(w,init(1,d))
    push!(w,bias(1,1))
    push!(w,init(d,d))
    push!(w,bias(d,1))
    !gating  && return w
    push!(w,init(1,d))
    push!(w,bias(1,1))
    return w
end

function mac(w,cw,qi,KBhw,KBhw′′,ci₋1,mi₋1,cj,mj,pad;train=false,selfattn=true,gating=true,tap=nothing)
    ci     = control_unit(w[1:4],ci₋1,qi,cw,pad;train=train,tap=tap)
    m_new  = read_unit(w[5:11],mi₋1,ci,KBhw,KBhw′′;train=train,tap=tap)
    mi     = write_unit(w[12:end],m_new,mi₋1,mj,ci,cj;train=train,selfattn=selfattn,gating=gating)
    tap != nothing && (tap["cnt"]+=1)
    return (ci,mi)
end

function init_mac(d;selfattn=true,gating=true)
   w  = Any[];
   append!(w,init_control(d))
   append!(w,init_read(d))
   append!(w,init_write(d;selfattn=selfattn,gating=gating))
   return w
end

function output_unit(w,q,mp;train=false,pdrop=0.0)
  x  = elu(w[1] * vcat(mp,q) .+ w[2])
  if train
      x  = dropout(x,pdrop) #0.15
  end
  y  = w[3] * x .+ w[4]
end

function init_output(d;outsize=28)
    w = Any[];
    push!(w,init(d,3d))
    push!(w,bias(d,1))
    push!(w,init(28,d))
    push!(w,bias(28,1))
    return w;
end

loss_layer(y,answers) = nll(y,answers;average=true)

function forward_net(w,r,qs,KB,batchSizes,pads,xB;answers=nothing,p=12,tap=nothing,selfattn=false,gating=false)
    train         = answers!=nothing

    KBhw          = postpocess_kb(w[1:4],KB;train=train)

    #READ UNIT PRE CALCULATIONS
    d,B,N         = size(KBhw)

    KBhw_2d       = reshape(KBhw,(d,B*N))

#     if train
#         KBhw_2d   = dropout(KBhw_2d,0)
#     end

    #KBhw′_pre     = w[15]*KBhw_2d .+ w[16] # look if it is necessary

    KBhw′′           = w[17]*KBhw_2d   .+  w[18]

    #KBhw′         = reshape(KBhw_2d,(d,B,N))

    #END

    q,cws        = process_question(w[5:8],r,qs,batchSizes;train=train)

    ci           = w[end-1]*xB
    mi           = w[end]*xB
    cj           = ci
    mj           = mi
    qi_c         = w[9]*q .+ w[10]

    for i=1:p
        qi        = qi_c[(i-1)*d+1:i*d,:]
        if train
            ci = dropout(ci,0.15)
            mi = dropout(mi,0.15)
        end
        ci,mi     = mac(w[11:end-6],cws,qi,KBhw,KBhw′′,ci,mi,cj,mj,pads;tap=tap,train=train,selfattn=selfattn,gating=gating)
        if selfattn
            cj        = hcat(cj,ci)
            mj        = hcat(mj,mi)
        end
    end

    y = output_unit(w[end-5:end-2],q,mi;train=train)

    if answers==nothing
        predmat = Array{Float32}(y)
        tap!=nothing && get!(tap,"y",predmat)
        return mapslices(indmax,predmat,1)[1,:]
    else
        return loss_layer(y,answers)
    end
end

function init_network(vocab_size,embed_size,d;p=12,loadresnet=false,selfattn=false,gating=false)
    if loadresnet
        rsnt,m,meta = load_resnet(atype;stage=3);
    else
        rsnt,m,meta = nothing,nothing,nothing;
    end
    w           = Any[];
    wcnn        = init_postkb(d)
    append!(w,wcnn);
    wrnn,r      = init_rnn(vocab_size,embed_size,d)
    append!(w,wrnn);
    #for i=1:p
        push!(w,init(p*d,2d),bias(p*d,1)) #qi embbedding
    #end
    wmac = init_mac(d;selfattn=selfattn,gating=gating) #!!!share weights among cells
    append!(w,wmac);
    wout = init_output(d)
    append!(w,wout);
    push!(w,init_state(d,1;initial=:xavier))
    push!(w,init_state(d,1;initial=:randn)) # m0
    return w,r,rsnt,m,meta;
end

function init_state(d,B;initial=:zero)
    if initial == :zero
        x=zeros(Float32,d,B)
    elseif initial == :randn
        x=randn(Float32,d,B)
    elseif initial == :xavier
        x=xavier(Float32,d,B)
    end
    return atype(x)
end

loss = grad(forward_net)

function savemodel(filename,w,wrun,r,opts,p)
    save(filename,"w",w,"wrun",wrun,"r",r,"opts",opts,"p",p)
end

function loadmodel(filename;onlywrun=false)
    d = load(filename)
    if onlywrun
        wrun=d["wrun"];r=d["r"];opts=nothing;w=nothing;p=d["p"]
    else
        w=d["w"];wrun=d["wrun"];r=d["r"];opts=d["opts"];p=d["p"];
    end
    return w,wrun,r,opts,p;
end

function getQdata(dhome,set)
    JSON.parsefile(dhome*set*".json")
end

function invert(vocab)
       int2tok = Array{String}(length(vocab))
       for (k,v) in vocab; int2tok[v] = k; end
       return int2tok
end

function getDicts(dhome,dicfile)
    dic  = JSON.parsefile(dhome*dicfile*".json")
    qvoc = dic["word_dic"]
    avoc = dic["answer_dic"]
    i2w  = invert(qvoc)
    i2a  = invert(avoc)
    return qvoc,avoc,i2w,i2a
end

function loadFeatures(dhome,set)
    feats    = reinterpret(Float32,read(open(dhome*set*".bin")))
    reshape(feats,(14,14,1024,div(length(feats),200704)))
end

function miniBatch(data;shfl=true,srtd=false,B=32)
    B=32
    L = length(data)
    shfl && shuffle!(data)
    srtd && sort!(data;by=x->length(x[2]))
    batchs = [];
    for i=1:B:L
        b         = min(L-i+1,B)
        questions = Any[]
        answers   = zeros(Int,b)
        images    = Any[]
        families  = zeros(Int,b)

        for j=1:b
            crw = data[i+j-1]
            push!(questions,reverse(Array{Int}(crw[2]).+1))
            push!(images,parse(Int,crw[1][end-9:end-4])+1)
            answers[j]  = crw[3]+1
            families[j] = crw[4]
        end

        lngths     = length.(questions);
        srtindices = sortperm(lngths;rev=true)

        lngths     = lngths[srtindices]
        Tmax       = lngths[1]
        questions  = questions[srtindices]
        answers    = answers[srtindices]
        images     = images[srtindices]
        families   = families[srtindices]

        qs = Int[];
        batchSizes = Int[];
        pads = falses(b,Tmax)

        for k=1:b
           pads[k,lngths[k]+1:Tmax]=true
        end

        if sum(pads)==0
           pads=nothing
        end

        while true
            batch = 0
            for j=1:b
                if length(questions[j]) > 0
                    batch += 1
                    push!(qs,pop!(questions[j]))
                end
            end
            if batch != 0
                push!(batchSizes,batch)
            else
                break;
            end
        end
        push!(batchs,(images,qs,answers,batchSizes,pads,families))
    end
    return batchs
end

function loadTrainingData(dhome="data/")
    info("Loading pretrained features for train&val sets.
    It requires minimum 70GB RAM!!!")
    trnfeats = loadFeatures(dhome,"train")
    valfeats = loadFeatures(dhome,"val")
    info("Loading questions ...")
    trnqstns = getQdata(dhome,"train")
    valqstns = getQdata(dhome,"val")
    info("Loading dictionaries ... ")
    qvoc,avoc,i2w,i2a = getDicts(dhome,"dic")
    return (trnfeats,valfeats),(trnqstns,valqstns),(qvoc,avoc,i2w,i2a)
end

function loadDemoData(dhome="data/demo/")
    info("Loading demo features ...")
    feats = loadFeatures(dhome,"demo")
    info("Loading demo questions ...")
    qstns = getQdata(dhome,"demo")
    info("Loading dictionaries ...")
    dics = getDicts(dhome,"dic")
    return feats,qstns,dics
end

function modelrun(w,r,opts,data,feats;p=12,train=false,wrun=nothing,ema=Float32(0.999),prefix=string(now())[1:11])

    B        = 32 # x 2 = 64 : look at bottom
    KxB      = atype(ones(Float32,1,B))
    #Requires to create 32x2 batchsize
    cumgrads = nothing
    if train
        cumgrads = map(similar,w)
    end

    #Statistics: Add accuracy according to question families
    cnt=total=0.0
    L = length(data)
    println("Timer Starts");flush(STDOUT);tic()

    for i=1:L

        if i % 100 == 0
            toc();tic();
            println(@sprintf("%.2f Completed",100i/L))
            !train && println(@sprintf("%.2f Accuracy",100cnt/total))
            flush(STDOUT)
        end

        if i % 2250 == 0 && train
            savemodel(prefix*"W.jld",w,wrun,r,opts,p); info("model saved"); flush(STDOUT)
            gc(); Knet.gc(); gc();
        end

        filenames,qs,answers,batchSizes,pad,_ = data[i]

        if batchSizes[1] != B
            B    = batchSizes[1]
            KxB  = atype(ones(Float32,1,B))
        end

        X = Any[];
        for k=1:B
               push!(X,view(featdict,:,:,:,filenames[k]))
        end
        xs = cat(4,X...)

        Kxs  = convert(atype,xs)

        Kpad = pad==nothing ? nothing : atype(pad * Float32(1e22))

        if train

            grads = loss(w,r,qs,Kxs,batchSizes,Kpad,KxB;answers=answers,p=p)

            for (cu,gr) in zip(cumgrads,grads)
                axpy!(Float32(0.5),gr,cu);
            end

            if iseven(i)
                Knet.update!(w,cumgrads,opts)
                if wrun != nothing
                    for (wr,wi) in zip(wrun,w)
                        axpy!(1.0-ema,wi-wr,wr)
                    end
                end
                map(x->fill!(x,Float32(0.0)),cumgrads)
            end

        else
            preds  = forward_net(w,r,qs,Kxs,batchSizes,Kpad,KxB)
            cnt   += sum(preds.==answers)
            total += B
        end
    end
end

function train!(w,wrun,r,opts,sets,feats;epochs=10,p=12)
    info("Training Starts....")
    for i=1:epochs
        modelrun(w,r,opts,sets[1],feats[1];wrun=wrun,train=true,p=p)
        if iseven(i)
            modelrun(wrun,r,opts,sets[2],feats[2];train=false,p=p)
        end
    end
    return w,wrun,r,opts;
end

function train(sets,feats;epochs=10,lr=0.0001,mfile=nothing,p=12)
     if mfile==nothing
         w,r,rsnt,_,_ = init_network(90,300,512;loadresnet=false,p=p);
         wrun=deepcopy(w)
         opts = optimizers(w,Adam;lr=lr)
     else
         w,wrun,r,opts,p = loadmodel(mfile)
     end
     train!(w,wrun,r,opts,sets,feats;epochs=10,p=p)
     return w,wrun,r,opts;
end

function train(dhome="data/";mfile=nothing,epochs=10,lr=0.0001,p=12)
     feats,qdata,dics = loadTrainingData(dhome)
     sets = []
     for q in questions; push!(sets,miniBatch(q)); end
     qdata = nothing; gc();
     w,wrun,r,opts = train(sets,feats;epochs=epochs,lr=lr,mfile=mfile,p=p)
     return w,wrun,r,opts,sets,feats,dics;
end

function validate(wrun,r,valset,valfeats)
     modelrun(wrun,r,nothing,valset,valfeats;train=false)
end

function validate(mfile,valset,valfeats)
     _,wrun,r,_ = loadmodel(mfile)
     modelrun(wrun,r,nothing,valset,valfeats;train=false)
     return wrun,r
end

function validate(mfile,dhome)
     _,wrun,r,_ = loadmodel(mfile)
     valfeats   = loadFeatures(dhome,"val")
     qdata      = getQdata(dhome,"val")
     dics       = getDicts(dhome,"dic")
     valset     = minibatch(qdata)
     modelrun(wrun,r,nothing,valset,valfeats;train=false)
     return wrun,r,valset,valfeats
end

function singlerun(w,r,feat,question;p=12)
    KxB        = atype(ones(Float32,1,1))
    pad        = nothing
    batchSizes = ones(Int,length(question))
    Kxs        = convert(atype,feat)
    results    = Dict()
    results["cnt"] = 1
    forward_net(w,r,question,Kxs,batchSizes,pad,KxB;tap=results,p=p)
    prediction = indmax(results["y"])
    return results,prediction
end

function visualize(img,results;p=12)
    s_y,s_x = size(img)./14
    for k=1:p
        α = results["w_attn_$(k)"][:]
        println("step_$(k) most attn. wrds: ",i2w[question[sortperm(α;rev=true)[1:2]]])
        flush(STDOUT)
        display([RGB{N0f8}(α[i],α[i],α[i]) for i=1:length(α)]);
        hsvimg = convert.(HSV,img);
        attn = results["KB_attn_$(k)"]
        for i=1:14,j=1:14
            rngy          = floor(Int,(i-1)*s_y+1):floor(Int,min(i*s_y,320))
            rngx          = floor(Int,(j-1)*s_x+1):floor(Int,min(j*s_x,480))
            hsvimg[rngy,rngx]  = scalepixel.(hsvimg[rngy,rngx],attn[sub2ind((14,14),i,j)])
        end
        display(hsvimg)
    end
end

function scalepixel(pixel,scaler)
     return HSV(pixel.h,pixel.s,pixel.v+2*scaler)
end

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using Base.Test

using PEGParser

function test_simple1()
  @grammar grammar begin
    start = list_content{ _1 }
    list_content = (list_label + list_values){ _2 }
    list_values = list(content, "," + space)
    space = r"[ \t]+"
    list_label = "list:" + -space
    content = r"[^,]+"
  end

  tolist(node, cvalues, ::MatchRule{:default}) = cvalues
  tolist(node, cvalues, ::MatchRule{:content}) = node.value

  data = "list: a, b, c"
  (ast, pos, error) = parse(grammar, data)
  result = transform(tolist, ast)
  @test result == ["a","b","c"]
end

test_simple1()


text:
using ImagesAPI
using Test

@testset "ImagesAPI.jl" begin
    # Write your own tests here.
end


text:
const DNA_SEQ = UInt8(0)
const RNA_SEQ = UInt8(1)
const AA_SEQ = UInt8(2)

const DNA_ALPHABET = ['A', 'T', 'C', 'G', 'N']
# for convinience, include 'C' in RNA alphabet
const RNA_ALPHABET = ['A', 'T', 'C', 'G', 'U', 'N']
const DNA_COMPLEMENT = Dict('A'=>'T', 'T'=>'A', 'C'=>'G', 'G'=>'C', 'N'=>'N', 'n'=>'n', 'a'=>'t', 't'=>'a', 'c'=>'g', 'g'=>'c')

type Sequence
	seq::ASCIIString
	seqtype::UInt8
end

Sequence(s::ASCIIString) = Sequence(s, DNA_SEQ)

Base.length(s::Sequence) = length(s.seq)
Base.getindex(s::Sequence, i::Int64) = getindex(s.seq, i)
Base.getindex(s::Sequence, r::UnitRange{Int64}) = Sequence(getindex(s.seq, r), s.seqtype)
Base.getindex(s::Sequence, indx::AbstractArray{Int64,1}) = Sequence(getindex(s.seq, indx), s.seqtype)
Base.endof(s::Sequence) = length(s)
Base.reverse(s::Sequence) = Sequence(reverse(s.seq), s.seqtype)
Base.display(s::Sequence) = display_sequence(s)
Base.print(s::Sequence) = display_sequence(s)

function dna(str::ASCIIString)
	ret = validate_sequence(str, DNA_ALPHABET)
	if ret > 0
		error("Unsupported DNA nucleotide symbol at pos $ret, value:", str[ret])
	end
	return Sequence(str, DNA_SEQ)
end

function rna(str::ASCIIString)
	ret = validate_sequence(str, RNA_ALPHABET)
	if ret > 0
		error("Unsupported RNA nucleotide symbol at pos $ret, value:", str[ret])
	end
	return Sequence(str, RNA_SEQ)
end

function aa(str::ASCIIString)
	return Sequence(str, AA_SEQ)
end

function validate_sequence(str::ASCIIString, alphabet)
	for i in 1:length(str)
		if !(uppercase(str[i]) in alphabet)
			return i
		end
	end
	return 0
end

function complement(s::Sequence)
	if s.seqtype == AA_SEQ
		error("Amino acid sequence doesn't have reverse complement")
	end
	len = length(s)
	arr = Char['N' for i in 1:len]
	for i in 1:len
		if haskey(DNA_COMPLEMENT, s[i])
			arr[i] = DNA_COMPLEMENT[s[i]]
		end
	end
	return Sequence(ASCIIString(arr), s.seqtype)
end

function reverse_complement(s::Sequence)
	if s.seqtype == AA_SEQ
		error("Amino acid sequence doesn't have reverse complement")
	end
	len = length(s)
	arr = Char['N' for i in 1:len]
	for i in 1:len
		if haskey(DNA_COMPLEMENT, s[len-i+1])
			arr[i] = DNA_COMPLEMENT[s[len-i+1]]
		end
	end
	return Sequence(ASCIIString(arr), s.seqtype)
end

function transcribe(s::Sequence)
	if s.seqtype != DNA_SEQ
		error("Only DNA sequence can be transcribed")
	end
	len = length(s)
	arr = Char[s.seq[i] for i in 1:len]
	for i in 1:len
		if s[i] == 'T'
			arr[i] = 'U'
		end
	end
	return rna(ASCIIString(arr))
end

==(s1::Sequence, s2::Sequence) = (s1.seq == s2.seq) && (s1.seqtype == s2.seqtype)
-(s1::Sequence) = reverse(s1)
!(s1::Sequence) = complement(s1)
~(s1::Sequence) = reverse_complement(s1)

# because we may meet very long sequence
# it's better to display only the head and tail of sequence if it's very long
# if limit is set 0, then no limit is made, the whole sequence will be displayed
function display_sequence(s::Sequence, limit::Int64 = 1000)
	const prefix = Dict(DNA_SEQ=>"dna", RNA_SEQ=>"rna", AA_SEQ=>"aa")
	len = length(s)
	str = prefix[s.seqtype]*":"
	if len <= limit || limit <= 0
		str = str * s.seq
	else
		show = min(500, round(Int, limit/2))
		str = str * s.seq[1:show] * "......"
		str = str * s.seq[len-show:len]
	end
	println(str)
end


text:
###
### Kmer Type definition
###

# Include some basic tuple bitflipping ops - the secret sauce to efficiently
# manipping Kmer's static data.
include("tuple_bitflipping.jl")

"""
    Kmer{A<:Alphabet,K,N} <: BioSequence{A}

A parametric, immutable, bitstype for representing Kmers - short sequences.
Given the number of Kmers generated from raw sequencing reads, avoiding
repetetive memory allocation and triggering of garbage collection is important,
as is the ability to effectively pack Kmers into arrays and similar collections.

In practice that means we an immutable bitstype as the internal representation
of these sequences. Thankfully, this is not much of a limitation - kmers are
rarely manipulated and so by and large don't have to be mutable.

Excepting their immutability, they fulfill the rest of the API and behaviours
expected from a concrete `BioSequence` type, and non-mutating transformations
of the type are still defined.

!!! warning
    Given their immutability, `setindex` and mutating sequence transformations
    are not implemented for Kmers e.g. `reverse_complement!`.
!!! tip
    Note that some sequence transformations that are not mutating are
    available, since they can return a new kmer value as a result e.g.
    `reverse_complement`.
"""
struct Kmer{A<:Alphabet,K,N} <: BioSequence{A}
    data::NTuple{N,UInt64}

    function Kmer{A,K,N}(data::NTuple{N,UInt64}) where {A<:Alphabet,K,N}
        checkmer(Kmer{A,K,N})
        # TODO: Decide on whether this method should always mask the (64N - 2K)
        # MSBs of the input tuple, as we do that in quite a few cases before
        # calling this constructor: see the typemin, typemax, rand, and transformations.jl
        return new(_cliphead(n_unused(Kmer{A,K,N}) * BioSequences.bits_per_symbol(A()), data...))
    end
end

BioSequences.encoded_data(seq::Kmer{A,K,N}) where {A,K,N} = seq.data

# Create a blank ntuple of appropriate length for a given Kmer with N.
@inline blank_ntuple(::Type{Kmer{A,K,N}}) where {A,K,N} = ntuple(x -> zero(UInt64), Val{N}())

###
### Constructors
###

# Create a Mer from a sequence.
function Kmer{A,K,N}(nucs) where {A,K,N}
    seqlen = length(nucs)
    if seqlen != K
        throw(ArgumentError("nucs does not contain the correct number of nucleotides ($seqlen ≠ $K)"))
    end
    data = _build_kmer_data(nucs, Kmer{A,K,N})
    return Kmer{A,K,N}(data)
end

# Convenience version of function above so you don't have to work out correct N.
function Kmer{A,K}(nucs) where {A,K}
    T = kmertype(Kmer{A,K})
    return T(nucs)
end

# Convenience version of function above so you don't have to provide K.
# BUT THIS WILL GET K FROM RUNTIME VALUES - YA CODE IS GONNA BE SLLLOOOOOOWWWW.
Kmer{A}(nucs) where {A,K} = Kmer{A,length(nucs)}(nucs)

@inline function _build_kmer_data(nucs, ::Type{Kmer{A,K,N}}) where {A,K,N}
    checkmer(Kmer{A,K,N})
    # Construct the head.
    elements_per_chunk = div(64, BioSequences.bits_per_symbol(A()))
    bases_in_head = div(64 - (64N - (BioSequences.bits_per_symbol(A()) * K)), BioSequences.bits_per_symbol(A()))
    head = zero(UInt64)
    @inbounds for i in 1:bases_in_head
        nt = convert(eltype(Kmer{A,K,N}), nucs[i])
        if isgap(nt)
            throw(ArgumentError("cannot create a mer with gaps"))
        end
        head = (head << BioSequences.bits_per_symbol(A())) | UInt64(BioSequences.encode(A(), nt))
    end
    # And the rest of the sequence
    idx = Ref(bases_in_head + 1)

    tail = ntuple(Val{N - 1}()) do i
        Base.@_inline_meta
        body = zero(UInt64)
        @inbounds for i in 1:elements_per_chunk
            nt = convert(eltype(Kmer{A,K,N}), nucs[idx[]])
            if isgap(nt)
                throw(ArgumentError("cannot create a mer with gaps"))
            end
            body = (body << BioSequences.bits_per_symbol(A())) | UInt64(BioSequences.encode(A(), nt))
            idx[] += 1
        end
        return body
    end

    return (head, tail...)
end

@inline function _build_kmer_data2(::Type{Kmer{A,K,N}}, nucs) where {A,K,N}
    checkmer(Kmer{A,K,N})
    kmer_data = blank_ntuple(Kmer{A,K,N})
    stop = length(nucs)
    while i ≤ stop
        nt = convert(eltype(A), @inbounds nucs[i])
        bits = BioSequences.encode(A(), nt)
        kmer_data = leftshift_carry(kmer_data, 2, fbits)
        filled += 1
        if filled == K
            return kmer_data
        end
        i += 1
    end
    return nothing
end

@inline function Kmer(nts::Vararg{DNA,K}) where {K}
    T = kmertype(DNAKmer{K})
    data = _build_kmer_data(nts, T)
    return T(data)
end

@inline function Kmer(nts::Vararg{RNA,K}) where {K}
    T = kmertype(RNAKmer{K})
    data = _build_kmer_data(nts, T)
    return T(data)
end

@inline function Kmer(seq::String)
    seq′ = BioSequences.remove_newlines(seq)
    hast = false
    hasu = false
    for c in seq′
        hast |= ((c == 'T') | (c == 't'))
        hasu |= ((c == 'U') | (c == 'u'))
    end
    if (hast & hasu) | (!hast & !hasu)
        throw(ArgumentError("Can't detect alphabet type from string"))
    end
    A = ifelse(hast & !hasu, DNAAlphabet{2}, RNAAlphabet{2})
    return Kmer{A,length(seq′)}(seq′)
end

"""
    kmertype(::Type{Kmer{A,K}}) where {A,K}
Resolve and incomplete kmer typing, computing the N parameter of
`Kmer{A,K,N}`, given only `Kmer{A,K}`.
## Example
```julia
julia> DNAKmer{63}
Kmer{DNAAlphabet{2},63,N} where N
julia> kmertype(DNAKmer{63})
Kmer{DNAAlphabet{2},63,2}
```
"""
@inline function kmertype(::Type{Kmer{A,K}}) where {A,K}
    return Kmer{A,K,BioSequences.seq_data_len(A, K)}
end
@inline kmertype(::Type{Kmer{A,K,N}}) where {A,K,N} = Kmer{A,K,N}

# Aliases
"Shortcut for the type `Kmer{DNAAlphabet{2},K,N}`"
const DNAKmer{K,N} = Kmer{DNAAlphabet{2},K,N}

"Shortcut for the type `DNAKmer{27,1}`"
const DNAKmer27 = DNAKmer{27,1}

"Shortcut for the type `DNAKmer{31,1}`"
const DNAKmer31 = DNAKmer{31,1}

"Shortcut for the type `DNAKmer{63,2}`"
const DNAKmer63 = DNAKmer{63,2}

"Shortcut for the type `Kmer{RNAAlphabet{2},K,N}`"
const RNAKmer{K,N} = Kmer{RNAAlphabet{2},K,N}

"Shortcut for the type `RNAKmer{27,1}`"
const RNAKmer27 = RNAKmer{27,1}

"Shortcut for the type `RNAKmer{31,1}`"
const RNAKmer31 = RNAKmer{31,1}

"Shortcut for the type `RNAKmer{63,2}`"
const RNAKmer63 = RNAKmer{63,2}

"Shortcut for the type `Kmer{AminoAcidAlphabet,K,N}`"
const AAKmer{K,N} = Kmer{AminoAcidAlphabet,K,N}

const DNACodon = DNAKmer{3,1}
const RNACodon = RNAKmer{3,1}

###
### Base Functions
###

@inline ksize(::Type{Kmer{A,K,N}}) where {A,K,N} = K
@inline capacity(::Type{Kmer{A,K,N}}) where {A,K,N} = div(64N, BioSequences.bits_per_symbol(A()))
@inline capacity(seq::Kmer) = capacity(typeof(seq))
@inline n_unused(::Type{Kmer{A,K,N}}) where {A,K,N} = capacity(Kmer{A,K,N}) - K
@inline n_unused(seq::Kmer) = n_unused(typeof(seq))

"""
    checkmer(::Type{Kmer{A,K,N}}) where {A,K,N}

Internal method - enforces good kmer type parameterisation.

For a given Kmer{A,K,N} of length K, the number of words used to
represent it (N) should be the minimum needed to contain all K symbols,
no larger (wasteful) no smaller (just... wrong).

Because it is used on type parameters / variables, these conditions should be
checked at compile time, and the branches / error throws eliminated when the
parameterisation of the Kmer type is good.
"""
@inline function checkmer(::Type{Kmer{A,K,N}}) where {A,K,N}
    if K < 1
        throw(ArgumentError("Bad kmer parameterisation. K must be greater than 0."))
    end
    n = BioSequences.seq_data_len(A, K)
    if n !== N
        # This has been significantly changed conceptually from before. Now we
        # don't just check K, but *enforce* the most appropriate N for K.
        throw(ArgumentError("Bad kmer parameterisation. For K = $K, N should be $n"))
    end
end

@inline Base.length(x::Kmer{A,K,N}) where {A,K,N} = K
@inline Base.summary(x::Kmer{A,K,N}) where {A,K,N} = string(eltype(x), ' ', K, "-mer")

function Base.typemin(::Type{Kmer{A,K,N}}) where {A,K,N}
    return Kmer{A,K,N}(ntuple(i -> zero(UInt64), N))
end

function Base.typemax(::Type{Kmer{A,K,N}}) where {A,K,N}
    return Kmer{A,K,N}((typemax(UInt64), ntuple(i -> typemax(UInt64), N - 1)...))
end

function Base.rand(::Type{Kmer{A,K,N}}) where {A,K,N}
    return Kmer{A,K,N}((rand(UInt64), ntuple(i -> rand(UInt64), N - 1)...))
end

Base.rand(::Type{Kmer{A,K}}) where {A,K} = rand(kmertype(Kmer{A,K}))

function Base.rand(::Type{T}, size::Integer) where {T<:Kmer}
    return [rand(T) for _ in 1:size]
end

###
### Old Mer Base Functions - not transferred to new type.
###
#@inline encoded_data_type(::Type{Mer{A,K}}) where {A,K} = UInt64
#@inline encoded_data_type(::Type{BigMer{A,K}}) where {A,K} = UInt128
#@inline encoded_data_type(x::AbstractMer) = encoded_data_type(typeof(x))
#@inline encoded_data(x::AbstractMer) = reinterpret(encoded_data_type(typeof(x)), x)
#@inline ksize(::Type{T}) where {A,K,T<:AbstractMer{A,K}} = K
#@inline Base.unsigned(x::AbstractMer) = encoded_data(x)
#Base.:-(x::AbstractMer, y::Integer) = typeof(x)(encoded_data(x) - y % encoded_data_type(x))
#Base.:+(x::AbstractMer, y::Integer) = typeof(x)(encoded_data(x) + y % encoded_data_type(x))
#Base.:+(x::AbstractMer, y::AbstractMer) = y + x
#Alphabet(::Type{Mer{A,K} where A<:NucleicAcidAlphabet{2}}) where {K} = Any

include("indexing.jl")

#LongSequence{A}(x::Kmer{A,K,N}) where {A,K,N} = LongSequence{A}([nt for nt in x])
# Convenience method so as don't need to specify A in LongSequence{A}.
LongSequence(x::Kmer{A,K,N}) where {A,K,N} = LongSequence{A}(x)

include("predicates.jl")
include("counting.jl")
include("transformations.jl")

###
### Kmer de-bruijn neighbors
###

# TODO: Decide on this vs. old iterator pattern. I like the terseness of the code vs defining an iterator. Neither allocate.
fw_neighbors(kmer::Kmer{A,K,N}) where {A<:DNAAlphabet,K,N} = ntuple(i -> pushlast(kmer, ACGT[i]), Val{4}())
fw_neighbors(kmer::Kmer{A,K,N}) where {A<:RNAAlphabet,K,N} = ntuple(i -> pushlast(kmer, ACGU[i]), Val{4}())
bw_neighbors(kmer::Kmer{A,K,N}) where {A<:DNAAlphabet,K,N} = ntuple(i -> pushfirst(kmer, ACGT[i]), Val{4}())
bw_neighbors(kmer::Kmer{A,K,N}) where {A<:RNAAlphabet,K,N} = ntuple(i -> pushfirst(kmer, ACGU[i]), Val{4}())

#=
# Neighbors on a de Bruijn graph
struct KmerNeighborIterator{S<:Kmer}
    x::S
end

"""
    neighbors(kmer::S) where {S<:Kmer}

Return an iterator through skip-mers neighboring `skipmer` on a de Bruijn graph.
"""
neighbors(kmer::Kmer) = KmerNeighborIterator{typeof(kmer)}(kmer)

Base.length(::KmerNeighborIterator) = 4
Base.eltype(::Type{KmerNeighborIterator{S}}) where {S<:Kmer} = S

function Base.iterate(it::KmerNeighborIterator{S}, i::UInt64 = 0) where {S<:Kmer}
    if i == 4
        return nothing
    else
        #return S((encoded_data(it.x) << 2) | i), i + 1
        return it.x << 1, i + one(UInt64)
    end
end
=#

###
### String literals
###

macro mer_str(seq, flag)
    seq′ = BioSequences.remove_newlines(seq)
    if flag == "dna" || flag == "d"
        T = kmertype(DNAKmer{length(seq′)})
        return T(seq′)
    elseif flag == "rna" || flag == "r"
        T = kmertype(RNAKmer{length(seq′)})
        return T(seq′)
    elseif flag == "aa" || flag == "a" || flag == "prot" || flag == "p"
        T = kmertype(AAKmer{length(seq′)})
        return T(seq′)
    else
        error("Invalid type flag: '$(flag)'")
    end
end

macro mer_str(seq)
    seq′ = remove_newlines(seq)
    T = kmertype(DNAKmer{length(seq′)})
    return T(seq′)
end


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

"""
    ReliabilityBranching

Branching strategy that uses pseudocosts if there are enough observations
to make an accurate prediction of strong branching scores, or runs the
actual strong branching routine if not enough observations are available.
"""
Base.@kwdef mutable struct ReliabilityBranching <: VariableBranchingRule
    min_samples::Int = 8
    max_sb_calls::Int = 100
    look_ahead::Int = 10
    n_sb_calls::Int = 0
    side_effect::Bool = true
end

function find_branching_var(
    rule::ReliabilityBranching,
    node::Node,
    pool::NodePool,
)::Variable
    nfrac = length(node.fractional_variables)
    pseudocost_scores = [
        _pseudocost_score(
            node,
            pool,
            node.fractional_variables[j],
            node.fractional_values[j],
        ) for j = 1:nfrac
    ]
    σ = sortperm(pseudocost_scores, rev = true)
    sorted_vars = node.fractional_variables[σ]
    _strong_branch_start(node)
    no_improv_count, n_sb_calls = 0, 0
    max_score, max_var = pseudocost_scores[σ[1]], sorted_vars[1]
    for (i, var) in enumerate(sorted_vars)
        use_strong_branch = true
        if n_sb_calls >= rule.max_sb_calls
            use_strong_branch = false
        else
            if var in keys(pool.var_history)
                varhist = pool.var_history[var]
                hlength = min(length(varhist.obj_ratio_up), length(varhist.obj_ratio_down))
                if hlength >= rule.min_samples
                    use_strong_branch = false
                end
            end
        end
        if use_strong_branch
            n_sb_calls += 1
            rule.n_sb_calls += 1
            score = _strong_branch_score(
                node = node,
                pool = pool,
                var = var,
                x = node.fractional_values[σ[i]],
                side_effect = rule.side_effect,
            )
        else
            score = pseudocost_scores[σ[i]]
        end
        if score > max_score
            max_score, max_var = score, var
            no_improv_count = 0
        else
            no_improv_count += 1
        end
        no_improv_count <= rule.look_ahead || break
    end
    _strong_branch_end(node)
    return max_var
end


text:
include("linear_and_integrator.jl")
include("backstepping_position_controller_static_allocator_multicopter.jl")


text:
module LazyArtifacts

# reexport the Artifacts API
using Pkg.Artifacts: Artifacts,
       artifact_exists, artifact_path, artifact_meta, artifact_hash,
       find_artifacts_toml, @artifact_str
export artifact_exists, artifact_path, artifact_meta, artifact_hash,
       find_artifacts_toml, @artifact_str

# define a function for satisfying lazy Artifact downloads
using Pkg.Artifacts: ensure_artifact_installed

end


text:
using FastTransforms, Test

FastTransforms.set_num_threads(ceil(Int, Base.Sys.CPU_THREADS/2))

@testset "libfasttransforms" begin
    n = 64
    α, β, γ, δ, λ, μ = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6
    for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
        x = T(1)./(1:n)
        Id = Matrix{T}(I, n, n)
        for (p1, p2) in ((plan_leg2cheb(Id), plan_cheb2leg(Id)),
                          (plan_ultra2ultra(Id, λ, μ), plan_ultra2ultra(Id, μ, λ)),
                          (plan_jac2jac(Id, α, β, γ, δ), plan_jac2jac(Id, γ, δ, α, β)),
                          (plan_lag2lag(Id, α, β), plan_lag2lag(Id, β, α)),
                          (plan_jac2ultra(Id, α, β, λ), plan_ultra2jac(Id, λ, α, β)),
                          (plan_jac2cheb(Id, α, β), plan_cheb2jac(Id, α, β)),
                          (plan_ultra2cheb(Id, λ), plan_cheb2ultra(Id, λ)))
            y = p1*x
            z = p2*y
            @test z ≈ x
            y = p1*x
            z = p1'y
            y = transpose(p1)*z
            z = transpose(p1)\y
            y = p1'\z
            z = p1\y
            @test z ≈ x
            y = p2*x
            z = p2'y
            y = transpose(p2)*z
            z = transpose(p2)\y
            y = p2'\z
            z = p2\y
            @test z ≈ x
            P = p1*Id
            Q = p2*P
            @test Q ≈ Id
            P = p1*Id
            Q = p1'P
            P = transpose(p1)*Q
            Q = transpose(p1)\P
            P = p1'\Q
            Q = p1\P
            @test Q ≈ Id
            P = p2*Id
            Q = p2'P
            P = transpose(p2)*Q
            Q = transpose(p2)\P
            P = p2'\Q
            Q = p2\P
            @test Q ≈ Id
        end
    end

    for T in (BigFloat, Complex{BigFloat})
        x = T(1)./(1:n)
        Id = Matrix{T}(I, n, n)
        for (p1, p2) in ((plan_leg2cheb(Id), plan_cheb2leg(Id)),
                          (plan_ultra2ultra(Id, λ, μ), plan_ultra2ultra(Id, μ, λ)),
                          (plan_jac2jac(Id, α, β, γ, δ), plan_jac2jac(Id, γ, δ, α, β)),
                          (plan_lag2lag(Id, α, β), plan_lag2lag(Id, β, α)),
                          (plan_jac2ultra(Id, α, β, λ), plan_ultra2jac(Id, λ, α, β)),
                          (plan_jac2cheb(Id, α, β), plan_cheb2jac(Id, α, β)),
                          (plan_ultra2cheb(Id, λ), plan_cheb2ultra(Id, λ)))
            y = p1*x
            z = p2*y
            @test z ≈ x
            y = p1*x
            z = p1'y
            y = transpose(p1)*z
            z = transpose(p1)\y
            y = p1'\z
            z = p1\y
            @test z ≈ x
            y = p2*x
            z = p2'y
            y = transpose(p2)*z
            z = transpose(p2)\y
            y = p2'\z
            z = p2\y
            @test z ≈ x
            P = p1*Id
            Q = p2*P
            @test_skip Q ≈ Id
            P = p1*Id
            Q = p1'P
            P = transpose(p1)*Q
            Q = transpose(p1)\P
            P = p1'\Q
            Q = p1\P
            @test_skip Q ≈ Id
            P = p2*Id
            Q = p2'P
            P = transpose(p2)*Q
            Q = transpose(p2)\P
            P = p2'\Q
            Q = p2\P
            @test_skip Q ≈ Id
        end
    end

    for T in (Float64, Complex{Float64})
        A = T <: Real ? sphones(T, n, 2n-1) : sphones(T, n, 2n-1) + im*sphones(T, n, 2n-1)
        p = plan_sph2fourier(A)
        ps = plan_sph_synthesis(A)
        pa = plan_sph_analysis(A)
        B = copy(A)
        C = ps*(p*A)
        A = p\(pa*C)
        @test A ≈ B

        A = T <: Real ? sphvones(T, n, 2n-1) : sphvones(T, n, 2n-1) + im*sphvones(T, n, 2n-1)
        p = plan_sphv2fourier(A)
        ps = plan_sphv_synthesis(A)
        pa = plan_sphv_analysis(A)
        B = copy(A)
        C = ps*(p*A)
        A = p\(pa*C)
        @test A ≈ B

        A = T <: Real ? diskones(T, n, 4n-3) : diskones(T, n, 4n-3) + im*diskones(T, n, 4n-3)
        p = plan_disk2cxf(A)
        ps = plan_disk_synthesis(A)
        pa = plan_disk_analysis(A)
        B = copy(A)
        C = ps*(p*A)
        A = p\(pa*C)
        @test A ≈ B

        A = T <: Real ? triones(T, n, n) : triones(T, n, n) + im*triones(T, n, n)
        p = plan_tri2cheb(A, α, β, γ)
        ps = plan_tri_synthesis(A)
        pa = plan_tri_analysis(A)
        B = copy(A)
        C = ps*(p*A)
        A = p\(pa*C)
        @test A ≈ B

        A = T <: Real ? tetones(T, n, n, n) : tetones(T, n, n, n) + im*tetones(T, n, n, n)
        p = plan_tet2cheb(A, α, β, γ, δ)
        ps = plan_tet_synthesis(A)
        pa = plan_tet_analysis(A)
        B = copy(A)
        C = ps*(p*A)
        A = p\(pa*C)
        @test A ≈ B
    end
end


text:
# Steps to install Python in Atom
# https://python.plainenglish.io/6-simple-steps-for-a-python-setup-in-atom-ca3100711f62
using PyCall
math = pyimport("math")
math.sin(math.pi / 4)

#Format: YYY-MM-DD
run(`python UC_BAU_postprocess.py '2019-01-01' '2019-01-02'`)

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
using DataFrames
using CSV
using Cairo
using Fontconfig
using Gadfly

df = CSV.read("heat.all.tsv")
p = plot(
    df,
    xgroup = :DataSet,
    x = :Threads,
    y = Col.value(Symbol("Inner(Bad)"), :Inner, :Outer, :Collapse),
    color = Col.index(Symbol("Inner(Bad)"), :Inner, :Outer, :Collapse),
    Geom.subplot_grid(Geom.point, Geom.line),
    Guide.ylabel("Time(s)"),
)

img = PNG("times-plot.png", 20cm, 8cm)
draw(img, p)


text:
using Test
using CliOptions

include("testutils.jl")

@testset "HelpOption()" begin
    @testset "ctor; $(v[1])" for v in [
        ([""], ArgumentError),
        ([""], ArgumentError),
        (["a"], ArgumentError),
        (["-"], ArgumentError),
        (["-a"], ("-a",)),
        (["--"], ArgumentError),
        (["--foo"], ("--foo",)),
    ]
        names, expected = v
        if expected isa Type
            @test_throws expected HelpOption(names...)
        else
            option = HelpOption(names...)
            @test option.names == expected
        end
    end

    @testset "ctor; duplicates, $(v[1])" for v in [
        (["-f", "--foo"], (nothing, nothing)),
        (["-f", "-f"], (ArgumentError, "-f")),
    ]
        names, expected = v
        if expected[1] isa Type
            tr = @test_throws expected[1] HelpOption(names...)
            if tr isa Test.Pass
                msg = stringify(tr.value)
                @test msg == ("ArgumentError: Duplicate names for a HelpOption found: " *
                              expected[2])
            end
        else
            @test HelpOption(names...) !== nothing
        end
    end

    @testset "show(x); $(v[1])" for v in [
        ([], "HelpOption(:h,:help)"),
        (["-a"], "HelpOption(:a)"),
        (["-a", "--foo-bar"], "HelpOption(:a,:foo_bar)"),
    ]
        names, expected_repr = v
        option = HelpOption(names...)
        @test repr(option) == expected_repr
    end

    @testset "show(io, x)" begin
        let option = HelpOption()
            buf = IOBuffer()
            redirect_stdout(buf) do
                show(option)
            end
            @test String(take!(buf)) == "HelpOption(:h,:help)"
        end
    end

    @testset "consume!(); $(v[1])" for v in [
        ([], AssertionError),
        (["-ab"], AssertionError),
        (["-h"], true),
        (["--help"], true),
    ]
        args, expected = v
        option = HelpOption()
        d = Dict{String,Any}()
        ctx = CliOptions.ParseContext()
        if expected isa Type
            @test_throws expected CliOptions.consume!(d, option, args, ctx)
        else
            num_consumed = CliOptions.consume!(d, option, args, ctx)
            @test num_consumed == 1
            @test d["help"] == expected
        end
    end
end


text:
# Generic utilities

module Utils

export
		tau, τ, emptyindices, checkindexseries, isindexseries, deleteatmap,
		splicemap, subsetsequal, substrip, wraptext, ncols, nrows, Repeated

const tau = 2π
const τ = 2π

const emptyindices = Int[]

function checkindexseries(collection, I::AbstractVector{<:Integer})
	if ! isindexseries(I)
		error("expected sorted unique indices")
	end
	checkbounds(collection, I)
end

function isindexseries(I::AbstractVector{<:Integer})
	if isempty(I)
		return true
	end
	lastv = I[1]
	for i in 2:length(I)
		v = I[i]
		if ! (v > lastv)
			return false
		end
		lastv = v
	end
	true
end

function deleteatmap(n::Integer, I::AbstractVector{<:Integer})
	map = collect(1:n)
	for i in I
		map[i] = 0
	end
	offset = 0
	for i in 1:n
		if map[i] == 0
			offset -= 1
		else
			map[i] += offset
		end
	end
	map
end

function splicemap(n::Integer, range::AbstractUnitRange{<:Integer},
		nreplacement::Integer)
	map = collect(1:n)
	for i in range
		map[i] = 0
	end
	offset = nreplacement - length(range)
	for i in last(range)+1:n
		map[i] += offset
	end
	map
end

function substrip(str::AbstractString, i0::Integer, i1::Integer)
	while i0 <= i1 && str[i0] == ' '
		i0 += 1
	end
	while i1 > i0 && str[i1] == ' '
		i1 -= 1
	end
	str[i0:i1]
end

function subsetsequal(A, B, i0, i1)
	for i = i0:i1
		if A[i] != B[i]
			return false
		end
	end
	true
end

function wraptext(text::AbstractString, maxcols::Integer = 80)
	maxcols > 0 || error("maxcols: expected strictly positive value")
	text = strip(text)
	if length(text) <= maxcols
		return [text]
	end
	tokens = split(text)
	lines = Vector{String}[]
	linecols = maxcols
	for token in tokens
		tokencols = length(token)
		if linecols + tokencols + 1 > maxcols
			push!(lines, [token])
			linecols = tokencols
		else
			push!(lines[end], token)
			linecols += tokencols + 1
		end
	end
	[join(line, " ") for line in lines]
end

ncols(M::AbstractMatrix) = size(M, 2)

nrows(M::AbstractMatrix) = size(M, 1)

struct Repeated{T1,T2,N} <: AbstractArray{T1,N}
	val::T1
	size::T2

	Repeated{T1,T2,N}(val::T1, size::T2) where {T1,T2,N} = new(val, size)
end

Repeated(val::T1, size::T2) where {T1,T2} =
		Repeated{T1,T2,length(size)}(val, size)

Repeated(val::T1, n::Integer) where {T1} = Repeated(val, (n,))

Base.size(A::Repeated) = A.size

Base.getindex(A::Repeated, i::Int) = A.val

Base.IndexStyle(::Type{<:Repeated}) = IndexLinear()

end # module


text:
using AugmentedGaussianProcesses; const AGP = AugmentedGaussianProcesses
using Distributions,LinearAlgebra
using Random: seed!; seed!(123)
using Plots; gr();

N = 100; noise = 1e-1
N_test = 10000
X = sort(rand(N))
X_test = collect(range(0,length=N_test,stop=1))
k = RBFKernel(0.1)
K = Symmetric(kernelmatrix(reshape(X,:,1),k)+noise*Diagonal(I,N))
f = rand(MvNormal(zeros(N),K))
y = rand.(Poisson.(10*AGP.logistic.(f)))
function intplot(model,iter)
    p=plot()
    plot!(X,y,t=:scatter);
    pred=(proba_y(model,X_test).-0.5).*2; plot!(X_test,pred);
    # if minimum(pred[2500])
    title!("Iteration : $iter, var : $(AugmentedGaussianProcesses.getvariance(model.kernel[1]))")
    # fstar = model.fstar(X_test,covf=false); plot!(X_test,logit(fstar));
    display(p);
end

plot(X,y,t=:scatter)
poissonmodel= VGP(X,y,k,PoissonLikelihood(),AnalyticVI());
train!(logitmodel,iterations=300)#,callback=intplot)
poissonpred = predict_y(logitmodel,X_test)
plot!(X_test,poissonpred)
negbinmodel = VGP(X,y,k,NegBinomialLikelihood(10),AnalyticVI());
train!(negbinmodel,iterations=50)#,callback=intplot)
negbinpred = predict_y(negbinmodel,X_test)
# svmpred2 = (AugmentedGaussianProcesses.svmpredictproba(svmmodel,X_test).-0.5)*2
# svmfstar,covfstar = svmmodel.fstar(X_test,covf=true)
# y1 = AugmentedGaussianProcesses.svmlikelihood(svmfstar)
# y2 = AugmentedGaussianProcesses.svmlikelihood(-svmfstar)
# plot(X_test,y1+y2)
##
default(legendfontsize=14.0,xtickfontsize=10.0,ytickfontsize=10.0,ylims = (-0.1,20))
p = plot(X,y,t=:scatter,lab="Training Points")
plot!(X_test,poissonpred,lab="Poisson Expectation",lw=7.0)
plot!(X_test,negbinpred,lab="Negative Binomial Expectation",lw=7.0,legend=:topright)

display(p)
savefig(p,String(@__DIR__)*"/Events.png")
findmin(logitpred)
# a,b = logitmodel.fstar(X_test)
# using Expectations, QuadGK, BenchmarkTools

# @btime begin; for i in 1:N_test; d = Normal(a[i],sqrt(max(b[i],1e-8))); expectation(logit,d); end; end;
# @btime begin; for i in 1:N_test; d = Normal(a[i],sqrt(max(b[i],1e-8))); quadgk(x->pdf(d,x)*logit(x),-Inf,Inf); end; end;


text:
@inline function ExaTronKernels.dcopy(n::Int,dx::ROCDeviceArray{Float64,1},incx::Int,
                       dy::ROCDeviceArray{Float64,1},incy::Int)
    tx = workitemIdx().x

    # Ignore incx and incy for now.
    if tx <= n
        @inbounds dy[tx] = dx[tx]
    end
    AMDGPU.sync_workgroup()

    return
end


text:
"""Gaussian process examples."""

using Distributions
using LinearAlgebra
using PyPlot
using Random

μ(X, m) = [m(x) for x in X]
Σ(X, k) = [k(x, xp) for x in X, xp in X]
K(X, Xp, k) = [k(x, xp) for x in X, xp in Xp]

# Mean functions.
mean_const(β) = (x) -> β
mean_lin(β) = (x) -> dot(x, β)
mean_zero() = (x) -> 0.0

# TODO(kykim): More mean functions.
# - Periodic mean: m(x)=a cos(2πx/p)+b sin(2πx/p).
# - Polynomical mean.

# Covariance functions.
cov_const(σ) = (x, xp) -> σ
cov_exp(l) = (x, xp) -> exp(-(x - xp) / l)
cov_gamma_exp(γ, l) = (x, xp) -> exp(-((x - xp) / l)^γ)
cov_sq_exp(l) = (x, xp) -> exp(-(x - xp)^2 / (2 * l^2))

mutable struct GaussianProcess
    m  # mean
    k  # covariance function
    X  # design points
    y  # objective values
    ν  # noise variance
end

function mvnrand(μ, Σ, inflation=1e-6)
    N = MvNormal(μ, Σ + inflation * I)
    return rand(N)
end
Base.rand(gp, X) = mvrand(μ(X, gp.m), Σ(X, gp.k))

function predict(gp, X_pred)
    X_pred = (isa(X_pred, Vector) ? X_pred : [X_pred])
    m, k, ν = gp.m, gp.k, gp.ν
    tmp = K(X_pred, gp.X, k) / (K(gp.X, gp.X, k) + ν * I)
    μ_p = μ(X_pred, m) + tmp * (gp.y - μ(gp.X, m))
    S = K(X_pred, X_pred, k) - tmp * K(gp.X, X_pred, k)
    ν_p = diag(S) .+ eps()
    return (μ_p, ν_p)
end

function plot_gp(gp, lb, ub, true_y::Function=nothing)
    xl = [x for x in lb:0.05:ub]

    # Plot observations.
    plot(gp.X, gp.y, "o", markersize=3, label="obs")

    # Plot the true y if available.
    if true_y != nothing
        plot(xl, [true_y(x) for x in xl], label="true y", lw=0.5)
    end

    # Plot predicted function mean with confidence interval.
    gp_pred = predict(gp, xl)
    gp_pred_mean, gp_pred_var = first(gp_pred), last(gp_pred)
    plot(xl, gp_pred_mean, label="pred mean", lw=0.5)
    gp_pred_c95 = 1.96 * gp_pred_var
    fill_between(xl, gp_pred_mean - gp_pred_c95, gp_pred_mean + gp_pred_c95,
                 alpha=0.1)
    xlabel("x"); ylabel("y"); legend();
    # savefig("myplot.png")
end

#################################
# Below are experiment functions.
#################################

# An example from stor-i.github.io/GaussianProcesses.jl/latest/Regression.
function simple_regression()
    Random.seed!(0)

    n = 10
    x = 2π * rand(n)
    y = sin.(x) + 0.05 * rand(n)

    gp = GaussianProcess(mean_zero(), cov_sq_exp(1.0), x, y, 1.0)
    true_y = (x) -> sin(x)
    figure(1, figsize=(9.0, 6.0))
    plot_gp(gp, -2.0, 8.0, true_y)
end

simple_regression()


text:
"""
Transformers enable you to transform the output generated by a datasets.
Many transformer can be chained together.

```julia
ds = SomeDataset(...)
ds = ds |> Transformer1(...) |> Transformer2(...) |> Transformer3(...)
dl = Dataloader(ds)
```
"""
abstract type Transformer <: Dataset end
Base.length(t::Transformer) = length(t.ds)
function (t::Transformer)(ds::Dataset)
    t.ds = ds
    return t
end

"""
Add noise to samples in a dataset.

# Usage

```julia
ds = TestDataset((28,28,1),(10,),100)
ds = ds |> NoiseTransformer(noiselevel=0.05)
```
"""
mutable struct NoisingTransfomer <: Transformer
    ds::Union{Dataset, Nothing}
	noiselevel
	axis

    NoisingTransfomer(noiselevel=0.01; axis=[:X]) = new(nothing,noiselevel,axis)
end

genNoise(arr::AbstractArray, level::Number) = randn(eltype(arr), size(arr)) .* level

function Base.getindex(t::NoisingTransfomer, idx)
	X,Y = t.ds[idx]
	if :X in t.axis
		X += genNoise(X, t.noiselevel)
	end

	if :Y in t.axis
		Y += genNoise(Y, t.noiselevel)
	end

	return (X,Y)
end

"""
Normalize data before feeding it to a model. A typical usecase is to normalize
an image.

# Usage

```julia
ds = TestDataset((28,28,1),(10,),100) |> Normalize()
```
"""
mutable struct Normalizer <: Transformer
    ds::Union{Dataset, Nothing}
	dims::Array{Int}
	means
	stds
	axis

    Normalizer(means=0, stds=1, dims=[1,2]; axis=[:X]) =
		new(nothing,dims, means, stds, axis)
end

function normalize(x, means, stds)
	(x .- means) ./ stds
end

function Base.getindex(t::Normalizer, idx)
	X,Y = t.ds[idx]
	if :X in t.axis
		X += normalize(X, t.means, t.stds)
	end

	if :Y in t.axis
		Y += normalize(Y, t.means, t.stds)
	end

	return (X,Y)
end

"""
Randomly crop an image in the shape of WxHxC to a certain size.

# Usage

```julia
# Only crop X
ds = TestDataset((100,100,3),(10,),100)
ds = ds |> ImageCrop((50,50))

# Crop both X and Y
ds = TestDataset((100,100,3),(50,50,3),100)
ds = ds |> ImageCrop((50,50), (20,20))
```

"""
mutable struct ImageCrop <: Transformer
    ds::Union{Dataset, Nothing}
	shapeX
	shapeY

    ImageCrop(shapeX , shapeY=nothing) = new(nothing,shapeX, shapeY)
end

function imageCrop(arr::AbstractArray,shape)
	maxX = size(arr,1) - shape[1]
	maxY = size(arr,2) - shape[2]
	offX = rand(1:maxX)
	offY = rand(1:maxY)
	arr[offX:offX+shape[1]-1,offY:offY+shape[2]-1, :]
end

function Base.getindex(t::ImageCrop, idx)
	X,Y = t.ds[idx]
	if t.shapeX !== nothing
		X = imageCrop(X, t.shapeX)
	end

	if t.shapeY !== nothing
		Y = imageCrop(Y, t.shapeY)
	end

	return (X,Y)
end

function onehot(x, labels, dtype::Type)
	result = zeros(dtype, length(labels))
    result[findfirst(x .== labels)] = 1
    result
end

function onehot(X::AbstractArray, labels, dtype::Type)
	result = zeros(dtype, length(labels))
	for x in X
		result[findfirst(x .== labels)] = 1
	end
	result
end

"""
OneHot transformer for a single sample. Uses
the onehot function to perform the actual transoformation.

Supports both single and multi class encoding:

	5 => [0,0,0,0,0,1,0,0,0,0]
	[5,7] => [0,0,0,0,0,1,0,1,0,0]

# Usage

```julia
ds = mnistds() |> OneHotEncoder(0:9)
```
"""
mutable struct OneHotEncoder <: Transformer
    ds::Union{Dataset, Nothing}
	labels
	axis
	dtype::Type

    OneHotEncoder(labels; axis=[:Y], dtype=getContext().dtype) =
		new(nothing, labels, axis, dtype)
end

function Base.getindex(t::OneHotEncoder, idx)
	X,Y = t.ds[idx]
	if :X in t.axis
		X = onehot(X, t.labels, t.dtype)
	end

	if :Y in t.axis
		Y = onehot(Y, t.labels, t.dtype)
	end

	return (X,Y)
end

"""
Return a subset of samples of the underlying dataset. Optionally
you can specify a start index (default = 1).

# Usage

data = MyDataset() |> Subset(1:100)

"""
mutable struct Subset <: Transformer
    ds::Union{Dataset, Nothing}
	idxs::AbstractArray{Int}

    Subset(idxs::AbstractArray{Int}) = new(nothing, idxs)
end

Base.length(t::Subset) = length(t.idxs)
Base.getindex(t::Subset, idx) = t.ds[t.idxs[idx]]

"""
Split one dataset into two subsets, one for training and one for validation.
By default Split will shuffle the indexes first in order to get a representative
validaiton set. The default split is 80% for training and 20% for validation.

# Usage

```julia
data_train, data_valid = ImageDataset(...) |> Split(0.25, shuffle=false)
```
"""
mutable struct Split <: Transformer
    ds::Union{Dataset, Nothing}
	valid_perc::Float64
	shuffle::Bool

    Split(valid_perc=0.2; shuffle=true) =
		new(nothing, valid_perc, shuffle)
end

function (t::Split)(ds::Dataset)
    t.ds = ds
	maxl = length(ds)

	idxs = t.shuffle ? Random.shuffle(1:maxl) : 1:maxl
	cut = round(Int, maxl*(1.0-t.valid_perc))

    a, b = (Subset(idxs[1:cut]), Subset(idxs[cut+1:end]))
	return (a(ds), b(ds))
end

#=
# MiniBatch helper functions
=#
function update_mb!(arr::AbstractArray, elem::AbstractArray, idx::Int)
	@assert size(arr)[1:end-1] == size(elem) "$(size(arr)) $(size(elem))"
	idxs = Base.OneTo.(size(elem))
	arr[idxs..., idx] = elem
end

function update_mb!(t::Tuple, elems::Tuple, idx::Int)
	@assert length(t) == length(elems)
	for (arr,elem) in zip(t,elems)
		update_mb!(arr, elem, idx)
	end
end

# Create a minibatch that has a similar size to the structure returned
# from the dataset.
create_mb(arr::AbstractArray, batchsize) = similar(arr, size(arr)..., batchsize)
create_mb(t::Tuple, batchsize)= Tuple(collect(create_mb(elem, batchsize) for elem in t))

"""
A MiniBatch is a transformer that will collect multipel samples from the underlying
dataset and put it in a minibatch. So it will call a dataset to get a single sample
and combines them into a minibatch. The MiniBatch itself is an iterator.

It will do this using threading, so when a dataset will read some samples from disk
this won't become a botttleneck.

Example:

	data = ImageDataset(filenames,labels) |> MiniBatch(16, shuffle=False)
	fit!(workout, data)

"""
mutable struct MiniBatch <: Transformer
    ds::Union{Dataset, Nothing}
    batchsize::Int
    shuffle::Bool

    MiniBatch(batchsize=8; shuffle=true) =
        new(nothing, batchsize, shuffle)
end

Base.length(dl::MiniBatch) = length(dl.ds) ÷ dl.batchsize

function Base.iterate(dl::MiniBatch, state=undef)
    maxl = length(dl.ds)
    bs = dl.batchsize

    if state == undef
        idxs = dl.shuffle ? Random.shuffle(1:maxl) : 1:maxl
        state = (idxs,1)
    end
    idxs, count = state

    if count > (maxl-bs) return nothing end

	l = Threads.SpinLock()
	minibatch = nothing

    Threads.@threads for i in 1:bs

		idx = i + count - 1
		sample = dl.ds[idx]
		@assert sample isa Tuple "Datasets should return Tuples, not $(typeof(sample))"

		if minibatch === nothing
			Threads.lock(l)
			if minibatch === nothing
				minibatch = create_mb(sample, bs)
			end
			Threads.unlock(l)
		end

		update_mb!(minibatch, sample, i)
    end
	Threads.unlock(l)
    return ((minibatch), (idxs, count + bs))
end


text:
using CairoMakie
using GeoData
using NCDatasets

function squeeze(A)
    singleton_dims = tuple((d for d in 1:ndims(A) if size(A, d) == 1)...)
    return dropdims(A, dims=singleton_dims)
end

function plot_initial_state(sose_profiles, sose_grid, interpolated_profiles, grid, lat, lon, start_date; z_bottom=-grid.Lz, filepath)

    zc_sose = sose_grid.zC
    zf_sose = sose_grid.zF

    U₀ = sose_profiles.U[1, :]
    V₀ = sose_profiles.V[1, :]
    Θ₀ = sose_profiles.Θ[1, :]
    S₀ = sose_profiles.S[1, :]

    Ugeo₀ = sose_profiles.Ugeo[1, :]
    Vgeo₀ = sose_profiles.Vgeo[1, :]

    zc = grid.zC[1:grid.Nz]
    zf = grid.zF[1:grid.Nz+1]

    ℑU₀ = interpolated_profiles.U.(zc, 0)
    ℑV₀ = interpolated_profiles.V.(zc, 0)
    ℑΘ₀ = interpolated_profiles.Θ.(zc, 0)
    ℑS₀ = interpolated_profiles.S.(zc, 0)

    ℑUgeo₀ = interpolated_profiles.Ugeo.(zc, 0)
    ℑVgeo₀ = interpolated_profiles.Vgeo.(zc, 0)

    colors = ["dodgerblue2", "crimson", "forestgreen"]

    fig = Figure(resolution=(2500, 1080))

    ax_U = fig[1, 1] = Axis(fig, xlabel="m/s", ylabel="z (m)")
    lines!(ax_U, U₀, zc_sose, label="U (SOSE)", linewidth=3, color=colors[1])
    lines!(ax_U, ℑU₀, zc, label="U (SOSE, interpolated)", linewidth=3, color=colors[1], linestyle=:dash)
    lines!(ax_U, Ugeo₀, zc_sose, label="U (geo)", linewidth=3, color=colors[2])
    lines!(ax_U, ℑUgeo₀, zc, label="U (geo, interpolated)", linewidth=3, color=colors[2], linestyle=:dash)
    axislegend(ax_U, position=:rb, framevisible=false)
    ylims!(ax_U, (z_bottom, 0))

    ax_V = fig[1, 2] = Axis(fig, xlabel="m/s", ylabel="z (m)")
    lines!(ax_V, V₀, zc_sose, label="V (SOSE)", linewidth=3, color=colors[1])
    lines!(ax_V, ℑV₀, zc, label="V (SOSE, interpolated)", linewidth=3, color=colors[1], linestyle=:dash)
    lines!(ax_V, Vgeo₀, zc_sose, label="V (geo)", linewidth=3, color=colors[2])
    lines!(ax_V, ℑVgeo₀, zc, label="V (geo, interpolated)", linewidth=3, color=colors[2], linestyle=:dash)
    axislegend(ax_V, position=:rb, framevisible=false)
    ylims!(ax_V, (z_bottom, 0))

    ax_Θ = fig[1, 3] = Axis(fig, xlabel="°C", ylabel="z (m)")
    lines!(ax_Θ, Θ₀, zc_sose, label="Θ (SOSE)", linewidth=3, color=colors[3])
    lines!(ax_Θ, ℑΘ₀, zc, label="Θ (SOSE, interpolated)", linewidth=3, color=colors[3], linestyle=:dash)
    axislegend(ax_Θ, position=:rb, framevisible=false)
    ylims!(ax_Θ, (z_bottom, 0))
    xlims!(ax_Θ, extrema(z_bottom < -grid.Lz ? filter(!iszero, Θ₀) : ℑΘ₀))

    ax_S = fig[1, 4] = Axis(fig, xlabel="psu", ylabel="z (m)")
    lines!(ax_S, S₀, zc_sose, label="S (SOSE)", linewidth=3, color=colors[3])
    lines!(ax_S, ℑS₀, zc, label="S (SOSE, interpolated)", linewidth=3, color=colors[3], linestyle=:dash)
    axislegend(ax_S, position=:rb, framevisible=false)
    ylims!(ax_S, (z_bottom, 0))
    xlims!(ax_S, extrema(z_bottom < -grid.Lz ? filter(!iszero, S₀) : ℑS₀))

    plot_title = @sprintf("SOSE LESbrary.jl (%.2f°N, %.2f°E) initial state at %s", lat, lon, start_date)
    supertitle = fig[0, :] = Label(fig, plot_title, textsize=30)

    @info "Saving $filepath..."
    save(filepath, fig)

    return nothing
end

function plot_surface_forcings(filepath_prefix)

    ds_b = NCDstack(filepath_prefix * "_large_scale.nc")

    times = ds_b[:time] / day
    Nt = length(times)

    colors = ["dodgerblue2", "crimson", "forestgreen"]

    # Makie doesn't support DateTime plotting yet :(
    # date_times = [start_date + Millisecond(round(Int, 1000t)) for t in times]

    τx_SOSE = ds_b[:τx].data
    τy_SOSE = ds_b[:τy].data
    τ_SOSE = @. √(τx_SOSE^2 + τy_SOSE^2)
    QΘ_SOSE = ds_b[:QT].data
    QS_SOSE = ds_b[:QS].data
    mld_SOSE = ds_b[:mld_SOSE].data
    mld_LES = ds_b[:mld_LES].data

    fig = Figure(resolution=(1920, 1080))
    plot_title = @sprintf("Realistic SOSE LESbrary.jl (%.2f°N, %.2f°E) surface forcings", lat, lon)

    ax_τ = fig[1, 1] = Axis(fig, ylabel="N/m²")
    line_τx = lines!(ax_τ, times, τx_SOSE, label="τx", linewidth=3, color=colors[1])
    line_τy = lines!(ax_τ, times, τy_SOSE, label="τy", linewidth=3, color=colors[2])
    line_τ  = lines!(ax_τ, times, τ_SOSE, label="√(τx² + τy²)", linewidth=3, color=colors[3])
    axislegend(ax_τ, position=:rb, framevisible=false)
    xlims!(ax_τ, extrema(times))
    ylims!(ax_τ, extrema([extrema(τx_SOSE)..., extrema(τy_SOSE)..., extrema(τ_SOSE)...]))
    hidexdecorations!(ax_τ, grid=false)

    ax_QΘ = fig[2, 1] = Axis(fig, ylabel="QΘ (W/m²)")
    line_QΘ = lines!(ax_QΘ, times, QΘ_SOSE, linewidth=3, color=colors[3])
    xlims!(ax_QΘ, extrema(times))
    ylims!(ax_QΘ, extrema(QΘ_SOSE))
    hidexdecorations!(ax_τ, grid=false)

    ax_QS = fig[3, 1] = Axis(fig, xlabel="time (days)", ylabel="QS (kg/m²/s)")
    line_QS = lines!(ax_QS, times, QS_SOSE, linewidth=3, color=colors[3])
    xlims!(ax_QS, extrema(times))
    ylims!(ax_QS, extrema(QS_SOSE))

    ax_mld = fig[4, 1] = Axis(fig, xlabel="time (days)", ylabel="Mixed layer depth (m)")
    line_mld_SOSE = lines!(ax_mld, times, mld_SOSE, label="SOSE", linewidth=3, color=colors[3], linestyle=:dash)
    line_mld_LES = lines!(ax_mld, times, mld_LES, label="LES", linewidth=3, color=colors[3])
    axislegend(ax_mld, position=:rb, framevisible=false)
    xlims!(ax_mld, extrema(times))
    ylims!(ax_mld, extrema([extrema(mld_SOSE)..., extrema(mld_LES)...]))

    supertitle = fig[0, :] = Label(fig, plot_title, textsize=30)

    filepath = filepath_prefix * "_surface_forcings.png"
    save(filepath, fig)

    @info "Figure saved: $filepath"

    return nothing
end

function animate_fields(lat, lon, start_date, filepath_prefix; frame_skip=2, fps=30)

    ds_xy = NCDstack(filepath_prefix * "_surface.nc")
    ds_yz = NCDstack(filepath_prefix * "_slice.nc")

    times = ds_xy[:time]
    Nt = length(times)

    xc = ds_xy[:xC].data
    xf = ds_xy[:xF].data
    yc = ds_xy[:yC].data
    yf = ds_xy[:yF].data
    zc = ds_yz[:zC].data
    zf = ds_yz[:zF].data

    fig = Figure(resolution=(1920, 1080))

    u_max = max(maximum(abs, ds_xy[:u]), maximum(abs, ds_yz[:u]))
    v_max = max(maximum(abs, ds_xy[:v]), maximum(abs, ds_yz[:v]))
    w_max = max(maximum(abs, ds_xy[:w]), maximum(abs, ds_yz[:w]))
    U_max = max(u_max, v_max, w_max)
    U_lims = 0.5 .* (-U_max, +U_max)

    frame = Node(1)

    plot_title = @lift @sprintf("Realistic SOSE LESbrary.jl (%.2f°N, %.2f°E): time = %s", lat, lon, DateTime(start_date) + Millisecond(round(Int, 1000 * times[$frame])))

    u_xy = @lift ds_xy[:u][Ti=$frame].data |> squeeze
    v_xy = @lift ds_xy[:v][Ti=$frame].data |> squeeze
    w_xy = @lift ds_xy[:w][Ti=$frame].data |> squeeze
    T_xy = @lift ds_xy[:T][Ti=$frame].data |> squeeze
    S_xy = @lift ds_xy[:S][Ti=$frame].data |> squeeze

    u_yz = @lift ds_yz[:u][Ti=$frame].data |> squeeze
    v_yz = @lift ds_yz[:v][Ti=$frame].data |> squeeze
    w_yz = @lift ds_yz[:w][Ti=$frame].data |> squeeze
    T_yz = @lift ds_yz[:T][Ti=$frame].data |> squeeze
    S_yz = @lift ds_yz[:S][Ti=$frame].data |> squeeze

    ax_u_xy = fig[1, 1] = Axis(fig, title="u-velocity")
    hm_u_xy = heatmap!(ax_u_xy, xf, yc, u_xy, colormap=:balance, colorrange=U_lims)
    hidedecorations!(ax_u_xy)

    ax_v_xy = fig[1, 2] = Axis(fig, title="v-velocity")
    hm_v_xy = heatmap!(ax_v_xy, xc, yf, v_xy, colormap=:balance, colorrange=U_lims)
    hidedecorations!(ax_v_xy)

    ax_w_xy = fig[1, 3] = Axis(fig, title="w-velocity")
    hm_w_xy = heatmap!(ax_w_xy, xc, yc, w_xy, colormap=:balance, colorrange=U_lims)
    hidedecorations!(ax_w_xy)

    ax_T_xy = fig[1, 4] = Axis(fig, title="conservative temperature")
    hm_T_xy = heatmap!(ax_T_xy, xc, yc, T_xy, colormap=:thermal, colorrange=extrema(ds_xy[:T]))
    hidedecorations!(ax_T_xy)

    ax_S_xy = fig[1, 5] = Axis(fig, title="absolute salinity")
    hm_S_xy = heatmap!(ax_S_xy, xc, yc, S_xy, colormap=:haline, colorrange=extrema(ds_xy[:S]))
    hidedecorations!(ax_S_xy)

    ax_u_yz = fig[2, 1] = Axis(fig)
    hm_u_yz = heatmap!(ax_u_yz, xf, zc, u_yz, colormap=:balance, colorrange=U_lims)
    hidedecorations!(ax_u_yz)

    ax_v_yz = fig[2, 2] = Axis(fig)
    hm_v_yz = heatmap!(ax_v_yz, xc, zc, v_yz, colormap=:balance, colorrange=U_lims)
    hidedecorations!(ax_v_yz)

    ax_w_yz = fig[2, 3] = Axis(fig)
    hm_w_yz = heatmap!(ax_w_yz, xc, zf, w_yz, colormap=:balance, colorrange=U_lims)
    hidedecorations!(ax_w_yz)

    ax_T_yz = fig[2, 4] = Axis(fig)
    hm_T_yz = heatmap!(ax_T_yz, xc, zc, T_yz, colormap=:thermal, colorrange=extrema(ds_yz[:T]))
    hidedecorations!(ax_T_yz)

    ax_S_yz = fig[2, 5] = Axis(fig)
    hm_S_yz = heatmap!(ax_S_yz, xc, zc, S_yz, colormap=:haline, colorrange=extrema(ds_yz[:S]))
    hidedecorations!(ax_S_yz)

    supertitle = fig[0, :] = Label(fig, plot_title, textsize=30)

    filepath = filepath_prefix * "_surface_slice_movie.mp4"
    record(fig, filepath, 1:frame_skip:Nt, framerate=fps) do n
        @info "Animating surface and slice movie frame $n/$Nt..."
        frame[] = n
    end

    @info "Movie saved: $filepath"

    return nothing
end

function animate_profiles(lat, lon, start_date, filepath_prefix; frame_skip=2, fps=30)

    ds_p = NCDstack(filepath_prefix * "_profiles.nc")
    ds_b = NCDstack(filepath_prefix * "_large_scale.nc")

    times = ds_p[:time]
    Nt = length(times)

    zc = ds_p[:zC].data
    zf = ds_p[:zF].data

    fig = Figure(resolution=(2500, 1080))

    frame = Node(1)

    plot_title = @lift @sprintf("Realistic SOSE LESbrary.jl (%.2f°N, %.2f°E): time = %s", lat, lon, DateTime(start_date) + Millisecond(round(Int, 1000 * times[$frame])))

    U_LES = @lift ds_p[:U][Ti=$frame].data
    V_LES = @lift ds_p[:V][Ti=$frame].data
    T_LES = @lift ds_p[:T][Ti=$frame].data
    S_LES = @lift ds_p[:S][Ti=$frame].data
    B_LES = @lift ds_p[:B][Ti=$frame].data

    U_SOSE = @lift ds_b[:u][Ti=$frame].data
    V_SOSE = @lift ds_b[:v][Ti=$frame].data
    T_SOSE = @lift ds_b[:T][Ti=$frame].data
    S_SOSE = @lift ds_b[:S][Ti=$frame].data
    # ∂ρ∂z_SOSE = @lift ds_b[:∂ρ∂z][Ti=$frame].data

    U_geo = @lift ds_b[:Ugeo][Ti=$frame].data
    V_geo = @lift ds_b[:Vgeo][Ti=$frame].data

    # time_so_far = @lift ds_b[:time][1:$frame].data
    # τx_SOSE = @lift ds_b[:τx][Ti=1:$frame].data
    # τy_SOSE = @lift ds_b[:τy][Ti=1:$frame].data
    # QΘ_SOSE = @lift ds_b[:QT][Ti=1:$frame].data
    # QS_SOSE = @lift ds_b[:QS][Ti=1:$frame].data

    colors = ["dodgerblue2", "crimson", "forestgreen"]

    ax_U = fig[1, 1] = Axis(fig, xlabel="m/s", ylabel="z (m)")
    line_U_SOSE = lines!(ax_U, U_SOSE, zc, label="U (SOSE)", linewidth=3, color=colors[1], linestyle=:dash)
    line_U_geo  = lines!(ax_U, U_geo, zc, label="U (geo)", linewidth=3, color=colors[1], linestyle=:dot)
    line_U_LES  = lines!(ax_U, U_LES, zc, label="U (LES)", linewidth=3, color=colors[1])
    axislegend(ax_U, position=:rb, framevisible=false)
    xlims!(ax_U, extrema([extrema(ds_p[:U])..., extrema(ds_b[:u])..., extrema(ds_b[:Ugeo])...]))
    ylims!(ax_U, extrema(zf))

    ax_V = fig[1, 2] = Axis(fig, xlabel="m/s", ylabel="z (m)")
    line_V_SOSE = lines!(ax_V, V_SOSE, zc, label="V (SOSE)", linewidth=3, color=colors[2], linestyle=:dash)
    line_V_geo  = lines!(ax_V, V_geo, zc, label="V (geo)", linewidth=3, color=colors[2], linestyle=:dot)
    line_V_LES  = lines!(ax_V, V_LES, zc, label="V (LES)", linewidth=3, color=colors[2])
    axislegend(ax_V, position=:rb, framevisible=false)
    xlims!(ax_V, extrema([extrema(ds_p[:U])..., extrema(ds_b[:v])..., extrema(ds_b[:Vgeo])...]))
    ylims!(ax_V, extrema(zf))

    ax_T = fig[1, 3] = Axis(fig, xlabel="°C", ylabel="z (m)")
    line_T_SOSE = lines!(ax_T, T_SOSE, zc, label="Θ (SOSE)", linewidth=3, color=colors[3], linestyle=:dash)
    line_T_LES  = lines!(ax_T, T_LES, zc, label="Θ (LES)", linewidth=3, color=colors[3])
    axislegend(ax_T, position=:rb, framevisible=false)
    xlims!(ax_T, extrema([extrema(ds_p[:T])..., extrema(ds_b[:T])...]))
    ylims!(ax_T, extrema(zf))

    ax_S = fig[1, 4] = Axis(fig, xlabel="psu", ylabel="z (m)")
    line_S_SOSE = lines!(ax_S, S_SOSE, zc, label="S (SOSE)", linewidth=3, color=colors[3], linestyle=:dash)
    line_S_LES  = lines!(ax_S, S_LES, zc, label="S (LES)", linewidth=3, color=colors[3])
    axislegend(ax_S, position=:rb, framevisible=false)
    xlims!(ax_S, extrema([extrema(ds_p[:S])..., extrema(ds_b[:S])...]))
    ylims!(ax_S, extrema(zf))

    ax_B = fig[1, 5] = Axis(fig, xlabel="m/s²", ylabel="z (m)")
    line_B_LES  = lines!(ax_B, B_LES, zc, label="B (LES)", linewidth=3, color=colors[3])
    axislegend(ax_B, position=:rb, framevisible=false)
    xlims!(ax_B, extrema(ds_p[:B]))
    ylims!(ax_B, extrema(zf))

    # ax_N = fig[1, 6] = Axis(fig, xlabel="kg/m⁴", ylabel="z (m)")
    # line_N_SOSE  = lines!(ax_N, ∂ρ∂z_SOSE, zf, label="∂ρ∂z (SOSE)", linewidth=3, color=colors[3], linestyle=:dash)
    # axislegend(ax_N, position=:rb, framevisible=false)
    # xlims!(ax_N, extrema(ds_b[:∂ρ∂z]))
    # ylims!(ax_N, extrema(zf))

    # ax_τ = fig[2, :] = Axis(fig, ylabel="N/m²")
    # line_τx = lines!(ax_τ, time_so_far, τx_SOSE, label="τx", linewidth=3, color=colors[1])
    # line_τy = lines!(ax_τ, time_so_far, τy_SOSE, label="τy", linewidth=3, color=colors[2])
    # axislegend(ax_τ, position=:rb, framevisible=false)
    # xlims!(ax_τ, extrema(ds_b[:time]))
    # ylims!(ax_τ, extrema([extrema(ds_b[:τx])..., extrema(ds_b[:τy])...]))
    # ax_τ.height = Relative(0.15)

    # ax_QΘ = fig[3, :] = Axis(fig, ylabel="QΘ (W/m²)")
    # line_QΘ = lines!(ax_QΘ, time_so_far, QΘ_SOSE, linewidth=3, color=colors[1])
    # xlims!(ax_QΘ, extrema(ds_b[:time]))
    # ylims!(ax_QΘ, extrema(ds_b[:QT]))
    # ax_QΘ.height = Relative(0.15)

    # ax_QS = fig[4, :] = Axis(fig, ylabel="QS (kg/m²/s))")
    # line_QS = lines!(ax_QS, time_so_far, QS_SOSE, linewidth=3, color=colors[1])
    # xlims!(ax_QS, extrema(ds_b[:time]))
    # ylims!(ax_QS, extrema(ds_b[:QS]))
    # ax_QS.height = Relative(0.15)

    supertitle = fig[0, :] = Label(fig, plot_title, textsize=30)

    filepath = filepath_prefix * "_first_order_statistics.mp4"
    record(fig, filepath, 1:frame_skip:Nt, framerate=fps) do n
        @info "Animating first-order statistics movie frame $n/$Nt..."
        frame[] = n
    end

    @info "Movie saved: $filepath"

    return nothing
end


text:
using BenchmarkTools, HePPCAT
using Random, LinearAlgebra

const SUITE = BenchmarkGroup()

## todo: rework to benchmark individual updates instead

# SAGE
SUITE["sage"] = BenchmarkGroup()
rng = MersenneTwister(123)
nfull, vfull = (10, 40), (1, 4)
for d = 5:10:15, k = 1:3:4, L = 1:2
    n, v = nfull[1:L], vfull[1:L]
    F, Z = randn(rng, d, k), [randn(rng, k, nl) for nl in n]
    Y = [F * Zl + sqrt(vl) * randn(rng, d, nl) for (Zl, vl, nl) in zip(Z, v, n)]

    Yflat = collect.(eachcol(hcat(Y...)))
    F0 = randn(rng, d, k)

    for T = 0:10:20
        SUITE["sage"]["block: d=$d, k=$k, L=$L, T=$T"] = @benchmarkable HePPCAT.ppca($Y, $k, $T, $F0, $(Val(:sage)))
        SUITE["sage"]["flat:  d=$d, k=$k, L=$L, T=$T"] = @benchmarkable HePPCAT.ppca($Yflat, $k, $T, $F0, $(Val(:sage)))
    end
end

# MM
SUITE["mm"] = BenchmarkGroup()
rng = MersenneTwister(123)
nfull, vfull = (40, 10), (4, 1)
for d = 50:25:100, k = 1:2:3, L = 1:2
    n, v = nfull[1:L], vfull[1:L]
    F, Z = randn(rng, d, k), [randn(rng, k, nl) for nl in n]
    Y = [F * Zl + sqrt(vl) * randn(rng, d, nl) for (Zl, vl, nl) in zip(Z, v, n)]

    Yflat = collect.(eachcol(hcat(Y...)))
    F0 = randn(rng, d, k)

    for T = 0:10:10
        SUITE["mm"]["flat:  d=$d, k=$k, L=$L, T=$T"] = @benchmarkable HePPCAT.ppca($Yflat, $k, $T, $F0, $(Val(:mm)))
    end
end

# PGD
SUITE["pgd"] = BenchmarkGroup()
rng = MersenneTwister(123)
nfull, vfull = (40, 10), (4, 1)
for d = 50:25:100, k = 1:2:3, L = 1:2
    n, v = nfull[1:L], vfull[1:L]
    F, Z = randn(rng, d, k), [randn(rng, k, nl) for nl in n]
    Y = [F * Zl + sqrt(vl) * randn(rng, d, nl) for (Zl, vl, nl) in zip(Z, v, n)]

    Yflat = collect.(eachcol(hcat(Y...)))
    F0 = randn(rng, d, k)

    for T = 0:10:10
        SUITE["pgd"]["flat:  d=$d, k=$k, L=$L, T=$T"] = @benchmarkable HePPCAT.ppca($Yflat, $k, $T, $F0, $(Val(:pgd)))
    end
end

# SGD
SUITE["sgd"] = BenchmarkGroup()
rng = MersenneTwister(123)
nfull, vfull = (40, 10), (4, 1)
for d = 50:25:100, k = 1:2:3, L = 1:2
    n, v = nfull[1:L], vfull[1:L]
    F, Z = randn(rng, d, k), [randn(rng, k, nl) for nl in n]
    Y = [F * Zl + sqrt(vl) * randn(rng, d, nl) for (Zl, vl, nl) in zip(Z, v, n)]

    Yflat = collect.(eachcol(hcat(Y...)))
    F0 = randn(rng, d, k)

    for T = 0:10:10
        SUITE["sgd"]["flat:  d=$d, k=$k, L=$L, T=$T"] = @benchmarkable HePPCAT.ppca($Yflat, $k, $T, $F0, $(Val(:sgd)))
    end
end


text:
abstract type AbstractObject end
abstract type AbstractWorld end
abstract type AbstractLandmark <: AbstractObject end
abstract type AbstractSensor end
abstract type AbstractCamera <: AbstractSensor end
abstract type AbstractRobot <: AbstractObject end
abstract type AbstractAgent end
abstract type AbstractEstimator end
abstract type AbstractParticle end
abstract type AbstractMcl <: AbstractEstimator end
abstract type AbstractResetMcl <: AbstractMcl end
abstract type AbstractMDPAgent <: AbstractAgent end
# KalmanFilter isa AbstractEstimator
# Mcl, KldMcl, FastSlam1, FastSlam2 isa AbstractMcl
# ResetMcl, AMcl isa AbstractResetMcl

max_stars_repo_name: soblin/JuliaProbo

max_issues_repo_name: soblin/JuliaProbo

max_forks_repo_name: soblin/JuliaProbo


text:
__precompile__() # this module is safe to precompile
"""
   DifferentialExtinction
   Original Python code (solar_diffex.py) by Andrea Lin
   Ported to Julia by Eric Ford
   Based on Collier-Cameron et al. (2019) MNRAS 487, 1082.  doi:10.1093/mnras/stz1215
   See also Astronomical Algorithms (Meeus, 1991) p151 & p178
"""
module DifferentialExtinction
using ..SunAsAStar

using Dates
using DataFrames, CSV

import ..SunAsAStar: AU, R_sun, valid_obs

f_z(alt) = π/2-alt

function f_airmass(z)
	secz = sec(z)
	return secz * (1.0 - 0.0012*(secz^2 - 1.0))
end

default_k_extinction = 0.157482  # Fit for WIYN, revisit at some point?
"""fractional brightness diff (from sun center)"""
function f_epsilon(z_center, dist_sun; k = default_k_extinction ) #
        z_upper = z_center - asin(R_sun/dist_sun)
        z_lower = z_center + asin(R_sun/dist_sun)
        delta_airmass = f_airmass(z_upper) - f_airmass(z_lower)
        #=
		 k_ex = 0.157482  # WARNING: Hardwired for WIYN.
		 X = [ones(length(eph.columns["airmass"])) eph.columns["airmass"]]
		 y =  eph.columns["magextinct"]
		 coeff = X \y
		 k_ex = coeff[2]
		 =#
		return k * delta_airmass/2.0 * log(10.)/2.5
end

"parallactic angle"
function f_q(H_sun, z_center, lat)
	#lat = deg2rad(wiyn["lat"]) #* u.deg # observer's latitude # no longer hard-wired
	return asin(sin(H_sun) * cos(lat)/sin(z_center))
end

f_phi(pole_angle, q) = (pole_angle - q)

function f_pole_angle_ecl(alpha_sun)
	e_tilt = deg2rad(23.4) # ecliptic tilt
    asin(sin(e_tilt) * cos(alpha_sun))
end

f_psi(pole_angle_ecl, q) = (pole_angle_ecl - q)

 """
 f_B0(JD)
  see Astronomical Algorithms (Meeus, 1991) p151 & p178
 """
 function f_B0(JD)

    T = (JD - 2451545.) / 36525.
    L0 = deg2rad(280.46645 + 36000.76983*T + 0.0003032*T^2)
    MA = deg2rad(357.52910 + 35999.05030*T - 0.0001559*T^2 - 0.00000048*T^3)
    sun_center = (1.914600 - 0.004817*T - 0.000014*T^2) * sin(MA) +
    			 (0.019993 - 0.000101*T) * sin(2*MA) +
    			 0.000290 * sin(3*MA)
	sun_center = deg2rad(sun_center)
    sun_true_lon = L0 + sun_center
    lam = sun_true_lon - deg2rad(0.00569) - deg2rad(0.00478 * sin(deg2rad(125.04 - 1934.136*T)))
    kappa = deg2rad(73.6667 + 1.3958333 * (JD - 2396758.) / 36525.)
    return asin(sin(lam - kappa) * sin(deg2rad(7.25)))
end

function diff_extinct_corr(JD, alt, HA, RA, pole_angle, dist_sun, lat)
    z = f_z(deg2rad(alt))
    epsilon = f_epsilon(z, dist_sun)
    q = f_q(deg2rad(HA*15), z, deg2rad(lat))

    phi = (f_phi(deg2rad(pole_angle), q))
    psi = f_psi(f_pole_angle_ecl(deg2rad(RA)), q)

    inc_sol_axis = f_B0(JD)

	ld_coeff = 0.6 # limb darkening coeff (linear)

	# sidereal diff rotation (Snodgrass & Ulrich 1990)
	sr_A = 2.972e-6 #* u.rad / u.s
	sr_B = -0.484e-6 #* u.rad / u.s
	sr_C = -0.361e-6 #* u.rad / u.s

	omega_orb = (2*π) / (365.25 * 24. * 3600.) #* u.rad / u.s # omega of Earth around Sun

    diff_rot_fac = sr_A*(7*ld_coeff - 15.) / (20*(3-ld_coeff)) +
      sr_B*(19 *ld_coeff - 35 + cos(inc_sol_axis)^2 * (3*ld_coeff - 35)) / (280*(3-ld_coeff)) +
      sr_C*(187*ld_coeff - 315 - cos(inc_sol_axis)^2 * (630 - 118*ld_coeff) -
      105 *cos(inc_sol_axis)^4 * (ld_coeff-1.)) / (6720*(3-ld_coeff))

    delta_vr = (epsilon * sin(phi) * R_sun * sin(inc_sol_axis) * diff_rot_fac -
       (7*ld_coeff - 15) * epsilon * sin(psi) / (20*(3-ld_coeff)) * R_sun * omega_orb)
	delta_vr *= 1e3 # km/s -> m/s.  WARNING: Andrea's python code returns cm/s

    return delta_vr
end

"""
   calc_Δv_diff_extinction(bjd; obs, ...)
   calc_Δv_diff_extinction(start, stop; obs, dt, dt_unit, ...)
Estimate apparent Δv (m/s) in solar observations due to differential extinction.
"""
function calc_Δv_diff_extinction end

function calc_Δv_diff_extinction(epoch::Float64; obs::Symbol, max_airmass::Float64 = 10.0, verbose::Bool = false)
	calc_Δv_diff_extinction_astropy(epoch; obs=obs, max_airmass=max_airmass, verbose=verbose)
end

function calc_Δv_diff_extinction_horizons(;epochs::Union{Vector{Float64},Dict}, obs::Symbol, max_airmass::Float64 = 10.0, verbose::Bool = false)
	loc = get_obs_loc(obs)
	eph = Horizons.query_horizons(epochs=epochs,obs=obs,max_airmass=max_airmass, verbose=verbose)
	delta_vr = diff_extinct_corr.(eph.columns["datetime_jd"], eph.columns["EL"],
	 	eph.columns["hour_angle"], eph.columns["RA"],
		eph.columns["NPole_ang"], eph.columns["delta"].*AU, loc["lat"])
	return (;jd=eph.columns["datetime_jd"], Δv=delta_vr)
end

function calc_Δv_diff_extinction_astropy(epoch::Float64; obs::Symbol, max_airmass::Float64 = 10.0, verbose::Bool = false)
	@assert 1 < max_airmass <= 10.0
	sol = get_solar_info(epoch, obs=obs)
	ha = sol[:lst] .- sol[:ra]./15
	idx_visible = 0 < sol[:airmass] <= max_airmass
	#if !idx_visible return (jd=epoch, Δv = 0.0) end
	if !idx_visible return Δv = 0.0 end
	horizons = Horizons.get_horizons_output(epoch,obs=obs)
	#delta = SolarRotation.get_earth_sun_dist(epoch, obs=obs)
	loc = get_obs_loc(obs)
	delta_vr = diff_extinct_corr(epoch, sol[:alt], ha, sol[:ra], horizons[:sol_npole_ang], sol[:sol_dist_au].*AU, loc["lat"]) #eph["NPole_ang"], eph["delta"])
	return delta_vr
end

function calc_Δv_diff_extinction_astropy(epochs::Vector{Float64}, obs::Symbol, max_airmass::Float64 = 10.0, verbose::Bool = false)
	@assert 1 <= length(epochs) <= 10000
	@assert 1 < max_airmass <= 10.0
	sol = get_solar_info(epochs, obs=obs)
	ha = sol.lst .- sol.ra./15
	idx_visible = 0 .< sol.airmass .<= max_airmass
	horizons = Horizons.get_horizons_output.(epochs[idx_visible],obs=obs)
	delta = get_earth_sun_dist(epochs[idx_visible], obs=obs)  # CHECK:  OK TO REMOVE DOT?
	#delta = sun.sol_dist_au
	loc = get_obs_loc(obs)
	#println(" jd = ", epochs, " delta_h = ", horizons.sol_dist_au, " delta_ap = ", delta)
	#println(" diff = ", horizons.sol_dist_au.-delta)
	delta_vr = diff_extinct_corr.(epochs[idx_visible], sol.alt[idx_visible], ha[idx_visible], sol.ra[idx_visible], horizons.sol_npole_ang, horizons.sol_dist_au.*AU, loc["lat"]) #eph["NPole_ang"], eph["delta"])
	return delta_vr
end

function compare_Horizons_Astropy(epochs::Vector{Float64}, obs::Symbol; max_airmass::Float64 = 10.0, verbose::Bool = false)
	if verbose println("# compare_Horizons_Astropy: ") end
	loc = get_obs_loc(obs)
	eph = Horizons.query_horizons(epochs=epochs,obs=obs,max_airmass=max_airmass, verbose=verbose)
	delta_vr = diff_extinct_corr.(eph.columns["datetime_jd"], eph.columns["EL"],
	 	eph.columns["hour_angle"], eph.columns["RA"],
		eph.columns["NPole_ang"], eph.columns["delta"].*AU, loc["lat"])

	sol = get_solar_info(epochs, obs=obs)
	ha = sol.lst .- sol.ra./15
	idx_visible = 0 .< sol.airmass .<= max_airmass
	delta = get_earth_sun_dist.(epochs[idx_visible], obs=obs)

	if verbose
		println("Horizons airmass =", eph.columns["airmass"])
		println("Astropy  airmass =", sol.airmass[idx_visible])
		#println("Δjd = ", eph.columns["datetime_jd"], " vs ", epochs[idx_visible])
		println("Δjd = ", eph.columns["datetime_jd"].-epochs[idx_visible])
		println("Δalt = ", eph.columns["EL"].- sol.alt[idx_visible])
		println("Δha = ", eph.columns["hour_angle"].- ha[idx_visible])
		println("ΔRA = ", eph.columns["RA"].- sol.ra[idx_visible])
	end
	horizons = Horizons.get_horizons_output.(epochs[idx_visible],obs=obs)
	delta_vr_new = diff_extinct_corr.(epochs[idx_visible], sol.alt[idx_visible], ha[idx_visible], sol.ra[idx_visible], zeros(length(epochs[idx_visible])), fill(AU,length(epochs[idx_visible])), loc["lat"]) #eph["NPole_ang"], eph["delta"])

	if verbose
		println("delta_vr_old = ", delta_vr, "\ndelta_vr_new = ", delta_vr_new)
		println("delta vr = ", delta_vr.- delta_vr_new)
	end
	return (eph=eph, sol=sol, delta_vr_old=delta_vr, delta_vr_new=delta_vr_new)
end

#=
function calc_Δv_diff_extinction(bjd::Real; obs::Symbol, max_airmass::Float64 = 10.0, verbose::Bool = false)
	epochs = [bjd]
	(bjd_out, Δv) = calc_Δv_diff_extinction(epochs=epochs, obs=obs, max_airmass=max_airmass, verbose=verbose)
	return first(Δv)
end
=#

#=
function calc_Δv_diff_extinction(start_date::String, stop_date::String; obs::Symbol, dt::Integer = 1, dt_unit::String = "m", max_airmass::Float64 = 10.0, verbose::Bool = false)
	@assert 1 <= dt <= 24*60
	@assert dt_unit ∈ ["m", "h", "d"]
	@assert 1 < max_airmass <= 10.0
	epochs = Dict("start"=>start_date, "stop"=>stop_date, "step"=>string(dt) * dt_unit)
	#(bjd_out, Δv) = calc_Δv_diff_extinction_horizons(epochs=epochs, obs=obs, max_airmass=max_airmass, verbose=verbose)
	# Doesn't work yet
	(bjd_out, Δv) = calc_Δv_diff_extinction_astropy(epochs, obs=obs, max_airmass=max_airmass, verbose=verbose)

	return (bjd=bjd_out, Δv=Δv)
end
=#

#=
function calc_Δv_diff_extinction(start_bjd::Real, stop_bjd::Real; obs::Symbol, dt::Integer = 1, dt_unit::String = "m", max_airmass::Float64 = 10.0, verbose::Bool=false)
	start_date = string(Date(Dates.julian2datetime(start_bjd)))
	stop_date = string(Date(Dates.julian2datetime(stop_bjd)))
	calc_Δv_diff_extinction(start_date, stop_date, obs=obs, dt=dt, dt_unit=dt_unit, max_airmass=max_airmass, verbose=verbose)
end

function calc_Δv_diff_extinction_demo(; verbose::Bool = false)
	start_date = "2021-01-01"
	stop_date = "2021-02-11"
	calc_Δv_diff_extinction(start_date, stop_date, obs=:WIYN, verbose=verbose)
end
=#

export calc_Δv_diff_extinction

end # module DifferentialExtinction