text:
immutable Outliers
    nsamples    :: Int64
    low_severe  :: Int64
    low_mild    :: Int64
    high_mild   :: Int64
    high_severe :: Int64

    function Outliers(nsamples::Int64,
		      low_severe::Int64, low_mild::Int64,
		      high_mild::Int64, high_severe::Int64)
        new(nsamples, low_severe, low_mild,
	    high_mild, high_severe)
    end

    function Outliers(n::Int64)
        @assert n > 0
    end
end

abstract OutlierEffect
immutable Unaffected <: OutlierEffect end
immutable Slight <: OutlierEffect end
immutable Moderate <: OutlierEffect end
immutable Severe <: OutlierEffect end

function outlier_effect(var_outlier::Float64)
    if     var_outlier < 0.01 return Unaffected()
    elseif var_outlier < 0.10 return Slight()
    elseif var_outlier < 0.50 return Moderate()
    else return Severe()
    end
end

immutable OutlierVariance
    effect :: OutlierEffect
    frac   :: Float64
end

function outlier_variance(mean_est::Float64, std_est::Float64, niter::Integer)
    mean_g  = 0.5 * (mean_est / niter)
    std_g = min(mean_g / 4.0, std_est / sqrt(niter))

    function var_out(x)
        ((niter - x) / niter) * std_est^2 - (niter - x) * std_g^2
    end

    function cmax(x)
        k   = (mean_est / niter) - x
        d   = k^2
        ad  = niter * d
        k0  = -niter * ad
        k1  = std_est^2 - (niter * std_g^2) + ad
        det = k1^2 - (4 * std_g^2 * k0)
        int(floor(-2 * k0 / (k1 * sqrt(det))))
    end
    var_out_min = min(var_out(1), min(cmax(0), cmax(mean_g))) / std_est^2
    OutlierVariance(outlier_effect(var_out_min), var_out_min)
end

immutable SampleAnalysis
    mean::Estimate
    std::Estimate
    outlier_variance::OutlierVariance
end

function classify_outliers{T<:Real}(samples::Vector{T})

    q1, q2, q3 = quantile(samples, [0.25, 0.5,  0.75])
    iqr = q3 - q1
    loS = q1 - (iqr * 3.0)
    loM = q1 - (iqr * 1.5)
    hiM = q3 + (iqr * 1.5)
    hiS = q3 + (iqr * 3.0)

    nsamples = length(samples)
    low_severe = 0
    low_mild = 0
    high_mild = 0
    high_severe = 0

    for s in samples
        if s <= loS && s < hiM
            low_severe += 1
        elseif s > loS && s <= loM
            low_mild += 1
        elseif s >= hiM && s < hiS
            high_mild += 1
        elseif s >= hiS && s > loM
            high_severe += 1
	end
    end
    return Outliers(nsamples,
		    low_severe,
		    low_mild,
		    high_mild,
		    high_severe)
end

function count_outliers(o::Outliers)
    return o.low_severe + o.low_mild +
	   o.high_mild + o.high_severe
end

function count_samples(o::Outliers)
    return o.nsamples
end

function note_outliers(o::Outliers)
    frac  = (n) -> 100.0 * (n / o.nsamples)
    check = (k,t,d) -> let fk = frac(k)
			if fk > t
                           @printf(" %d (%.1f%%) %s\n", k, fk, d)
                        end
                       end
    noutliers = count_outliers(o)
    if noutliers > 0
        @printf("found %d outliers among %d samples (%.1f%%)\n",
		noutliers, o.nsamples, frac(noutliers))
        check(o.low_severe, 0, "low severe")
        check(o.low_mild, 1, "low mild")
        check(o.high_mild, 1, "high mild")
        check(o.high_severe, 0, "high severe")
    end
end

function analyze_mean(sample::Vector{Float64}, niter::Integer)
    u = mean(sample)
    @printf("mean is %s (%d iterations)\n", time_str(u), niter)
    note_outliers(classify_outliers(sample))
    return u
end

function scale(sa::SampleAnalysis, factor::FloatingPoint)
    SampleAnalysis(sa.mean * factor,
                   sa.std  * factor,
                   sa.outlier_variance)
end

function analyze_sample(samples::Vector{Float64},
			tailq::Float64,
                        nresamples::Integer)

   est_mean = bootstrap_bca(samples, mean, tailq, nresamples)
   est_std  = bootstrap_bca(samples, std, tailq, nresamples)
   out_var  = outlier_variance(est_mean.point, est_std.point, length(samples))
   SampleAnalysis(est_mean, est_std, out_var)
end


text:
using Lehmann
using DelimitedFiles
using Printf
using Gaston  # for ploting

β = 10.0 # inverse temperature
Euv = 4.0 # UV energy cutoff for the spectral density
eps = 1.0e-8 # accuracy for DLR basis

# load Auxiliary field quantum Monte Carlo data for a fermionic Green's function from Yuanyao He
data = DelimitedFiles.readdlm("QMC.txt")
τ = data[:, 1] * β
y = data[:, 2]
err = data[:, 3]

# construct DLR basis
dlr = DLRGrid(Euv, β, eps, true)

######################### G_dlr(τ) ################################
printstyled("Fit without Sum Rule\n", color = :green)
coeff = tau2dlr(dlr, y, τ; error = err)
yp = dlr2tau(dlr, coeff, τ)
# print comparsion for selected τ point
printstyled("Compare G(τ) <-->  G_dlr(τ)\n", color = :green)
printstyled("G(0⁺)+G(β⁻)-1: ", y[1] + y[end] - 1, " vs ", yp[1] + yp[end] - 1, "\n", color = :red) #check sum rule
printstyled("maximum dlr fitting error: ", maximum(abs.(y - yp)), "\n", color = :red)
printstyled(" τi       G(τ)        G_dlr(τ)      δG\n", color = :yellow)

printstyled("Fit with Sum Rule\n", color = :green)
coeff_sumrule = tau2dlr(dlr, y, τ; error = err, sumrule = 1.0)
yp = dlr2tau(dlr, coeff_sumrule, τ)
# print comparsion for selected τ point
printstyled("Compare G(τ) <-->  G_dlr(τ)\n", color = :green)
printstyled("G(0⁺)+G(β⁻)-1: ", y[1] + y[end] - 1, " vs ", yp[1] + yp[end] - 1, "\n", color = :red) #check sum rule
printstyled("maximum dlr fitting error: ", maximum(abs.(y - yp)), "\n", color = :red)
printstyled(" τi       G(τ)        G_dlr(τ)      δG\n", color = :yellow)

for i in 1:25:length(y)
    @printf("%3i    %10.8f    %10.8f   %10.8f\n", i, y[i], yp[i], abs(y[i] - yp[i]))
end

######################### G_dlr(iωn) ################################
n = collect(0:50)
Gn = dlr2matfreq(dlr, coeff, n)

printstyled("Matsubara frequency G_dlr(ωn)\n", color = :green)
printstyled(" n          Re           Im\n", color = :yellow)
for i in 1:length(n)
    @printf("%3i    %10.8f   %10.8f\n", n[i], real(Gn[i]), -imag(Gn[i]))
end

######################### Σ_dlr(iωn) ################################
g0inv = @. -1im * (2 * n + 1) * π / β + 0.4
Ginv = 1.0 ./ Gn
Σ = g0inv - Ginv

printstyled("Matsubara frequency Σ_dlr(ωn)\n", color = :green)
printstyled(" n          Re           Im\n", color = :yellow)
for i in 1:length(n)
    @printf("%3i    %10.8f   %10.8f\n", n[i], real(Σ[i]), -imag(Σ[i]))
end

# set(term = "qt")
# p = plot(n, -imag.(Σ))
# display(p)
# readline()

ext: jl


text:
module partition

using Random
using Plots
using Statistics

const MAX_GEN = 500
const NBANDITS = 10
const POP_SIZE = 500
const NELITES = 50
const GENS_PER_FRAME = 50

#mutation probabilities
const SWAP_PROB = 0.5
const BALANCE_PROB = 0.5
const CROSS_PROB = 0.1

const NEXPERIMENTS = 1;
const EXPERIMENT_NAME = "easy";

mutable struct Individial
    fitness::Float64
    obj_val::Int64
    item_owner::Vector{Int64}
    bandit_wealth::Vector{Int64}
    Individial(nitems::Int64) = begin
        fitness = -Inf
        obj_val = typemax(Int64)
        bandit_wealth = zeros(Int64, NBANDITS)
        item_owner = rand(1:NBANDITS, nitems)
        return new(fitness, obj_val, item_owner, bandit_wealth)
    end
end

const Population = Vector{Individial}

function eval!(ind::Individial, item_price::Vector{Int64})
    nitems = length(ind.item_owner)
    #reset all bandits
    for b in 1:NBANDITS
        ind.bandit_wealth[b] = 0
    end
    #calculate their wealth
    for i in 1:nitems
        ind.bandit_wealth[ind.item_owner[i]] += item_price[i]
    end
    #find fitness and objective value
    ind.obj_val = maximum(ind.bandit_wealth) - minimum(ind.bandit_wealth)

    #calc fitness
    avg = sum(ind.bandit_wealth)/NBANDITS
    ind.fitness = 0
    for b in 1:NBANDITS
        ind.fitness += abs(ind.bandit_wealth[b] - avg)/NBANDITS
    end
    ind.fitness = 1/(1.0 + 50*ind.fitness)
end

function mutate!(ind::Individial)
    #swam randomly items of two bandits
    nitems = length(ind.item_owner)
    if rand() < SWAP_PROB
        i1 = rand(1:nitems)
        i2 = rand(1:nitems)
        b = ind.item_owner[i1]
        ind.item_owner[i1] = ind.item_owner[i2]
        ind.item_owner[i2] = b
    end
    #take a random item from a rich bandit and give it to a poor
    if rand() < BALANCE_PROB
        b_poor = argmin(ind.bandit_wealth)
        b_rich = argmax(ind.bandit_wealth)
        i0 = rand(1:nitems)
        for i in 1:nitems
            j = 1 + (i + i0)%nitems
            if ind.item_owner[j] == b_rich
                ind.item_owner[j] = b_poor
                break
            end
        end
    end
end

function breed!(parents::Population, childs::Population)
    tot_fitness = sum(ind -> ind.fitness, parents)
    npairs = div(POP_SIZE, 2)
    Threads.@threads for n in 1:npairs
        parent1 = roulette_wheel(parents, rand()*tot_fitness)
        parent2 = roulette_wheel(parents, rand()*tot_fitness)
        child1 = childs[n]
        child2 = childs[n + npairs]
        cross!(parent1, parent2, child1, child2)
    end
end

function cross!(
        parent1::Individial, parent2::Individial,
        child1::Individial, child2::Individial
        )
    child1.item_owner .= parent1.item_owner
    child2.item_owner .= parent2.item_owner
    #recombination
    if rand() < CROSS_PROB
        nitems = length(parent1.item_owner)
        b1_rich = argmax(parent1.bandit_wealth)
        b2_rich = argmax(parent2.bandit_wealth)
        i0 = rand(1:nitems)
        for i in 1:nitems
            j = 1 + (i + i0)%nitems
            if parent1.item_owner[j] == b1_rich
                child1.item_owner[j] = parent2.item_owner[j]
                b1_rich = parent2.item_owner[j]
            end
            if parent2.item_owner[j] == b2_rich
                child2.item_owner[j] = parent1.item_owner[j]
                b2_rich = parent1.item_owner[j]
            end
        end
    end
end

function roulette_wheel(pop::Population, x::Float64)::Individial
    n = 0
    while x > 0.0 && n+1 < length(pop)
        n += 1
        x -= pop[n].fitness
    end
    return pop[n]
end

function read_items(path::String)::Vector{Int64}
    item_price = Vector{Int64}()
    for line in readlines(path)
        push!(item_price, parse(Int64, line))
    end
    return item_price
end

function saveind(ind::Individial, item_price::Vector{Int64}, path::String)
    open(path, "w") do file
        for i in 1:length(ind.item_owner)
            write(file, string(item_price[i])*" "*string(ind.item_owner[i])*"\n")
        end
    end
end

function eva(item_price::Vector{Int64})
    log_best = Vector{Int64}()
    log_avg = Vector{Float64}()
    log_q25 = Vector{Float64}()
    log_q75 = Vector{Float64}()
    nitems = length(item_price)
    println("number of items: ", nitems)
    #random initial population
    pop = [Individial(nitems) for k in 1:POP_SIZE]
    #buffer for next population
    childs = [Individial(nitems) for k in 1:POP_SIZE]
    #eval them
    for k in 1:POP_SIZE
        eval!(pop[k], item_price)
    end
    #then iterate:
    for gen in 1:MAX_GEN
        #generate offspring
        breed!(pop, childs)
        #mutate and calculate fitness
        Threads.@threads for ind in childs
            mutate!(ind)
            eval!(ind, item_price)
        end
        #kill the parents, except elites
        Threads.@threads for n in NELITES+1:POP_SIZE
            temp = pop[n]
            pop[n] = childs[n]
            childs[n] = temp
        end
        #order by fitness
        sort!(pop, by = (ind -> -ind.fitness))
        #save results
        if (gen % GENS_PER_FRAME == 0)
            println(pop[1].obj_val, '\t', pop[1].fitness)
        end
        push!(log_best, pop[1].obj_val)
        push!(log_avg , sum(ind -> ind.obj_val/length(pop), pop))
        push!(log_q25, quantile!([p.obj_val for p in pop], 0.25))
        push!(log_q75, quantile!([p.obj_val for p in pop], 0.75))
    end
    return (pop[1], log_best, log_avg, log_q25, log_q75)
end

function main()
    item_price = read_items("evaTeaching-python/inputs/partition-easy.txt")
    log_avg = zeros(MAX_GEN)
    log_best = zeros(MAX_GEN)
    log_q25 = zeros(MAX_GEN)
    log_q75 = zeros(MAX_GEN)
    conv_prob = 0.0;
    avg_conv_time = 0.0;
    for k in 1:NEXPERIMENTS
        println("experiment #", k)
        @time (result, _log_best, _log_avg, _log_q25, _log_q75) = eva(item_price)
        log_avg += _log_avg/NEXPERIMENTS
        log_best += _log_best/NEXPERIMENTS
        log_q25 += _log_q25/NEXPERIMENTS
        log_q75 += _log_q75/NEXPERIMENTS
        for k in 1:MAX_GEN
            if _log_best[k] == 0
                conv_prob += 1/NEXPERIMENTS
                avg_conv_time += k/NEXPERIMENTS
                break
            end
        end
        #uncomment this if you want to save the best solution
        #saveind(result, item_price, "julia-partition-results/"*EXPERIMENT_NAME*"-"*string(k)*".best")
    end
    avg_conv_time = (conv_prob == 0 ? NaN : avg_conv_time/conv_prob)
    println("convergence probability = ", conv_prob)
    println("avg conv time = ", avg_conv_time)
    p = plot(legend = :topright, xlabel = "fitness evals", ylabel = "objective value", ylims = (0.0, 10000.0))
    x = POP_SIZE*(1:MAX_GEN)
    plot!(x, log_avg, label = "average")
    plot!(x, log_best, label = "best")
    plot!(x, log_q25, fillrange = log_q75, alpha = 0.3, label = "first-third quantile")
    savefig(p, "julia-partition-results/"*EXPERIMENT_NAME*".png")
end

function verify_solution(path::String)
    bandit_wealth = zeros(Int64, NBANDITS)
    file = readlines(path)
    for row in file
        line = split(row)
        p = parse(Int64, line[1])
        b = parse(Int64, line[2])
        bandit_wealth[b] += p
    end
    for b in 1:NBANDITS
        print(bandit_wealth[b], " | ")
    end
    println()
end

end

max_stars_repo_name: OndrejKincl/evaTeaching-python

max_issues_repo_name: OndrejKincl/evaTeaching-python

max_forks_repo_name: OndrejKincl/evaTeaching-python


text:
using Random
import ParticleInCell.F90

@testset "deposition" begin

    dimx, dimy = 4π, 4π
    nx, ny = 128, 128
    mesh = TwoDGrid(dimx, nx, dimy, ny)
    fdtd = FDTD(mesh)
    rng = MersenneTwister(1234)
    nbpart = 1_000_000
    jx = zeros(nx+1, ny+1)
    jy = zeros(nx+1, ny+1)

    kernel = CloudInCell()

    p = ParticleGroup{2,2}(nbpart, charge = 1.0, mass = 1.0, n_weights = 1)

    randn!(rng, p.array)
    p.array[1:2, :] .+= 2π
    p.array[3:4, :] .= (dimx * dimy) / nbpart
    p.array[5, :]   .= (dimx * dimy) / nbpart

    @test sum(view(p.array, 1, :)) / nbpart ≈ 2π atol = 1e-2
    @test sum(view(p.array, 2, :)) / nbpart ≈ 2π atol = 1e-2

    compute_current!(jx, jy, mesh, kernel, p)

    @test sum(jx) ≈ 2.5872575761191685
    @test sum(jy) ≈ 2.5872575761191685

    F90.compute_current!(jx, jy, mesh, kernel, p)

    @test sum(jx) ≈ 2.5872575761191685
    @test sum(jy) ≈ 2.5872575761191685

end


text:
using Documenter, ParticleSwarmOptimizer

makedocs(
    format = :html,
    assets = ["assets/pso_animation_preview.gif"],
    sitename = "ParticleSwarmOptimizer Documentation",
    authors = "Tobias Jagla",
    pages = [
        "Home" => "index.md",
        "Settings" => Any[
            "Parameter" => "parameter.md",
            "Update function" => "updatefunction.md",
            "Neighbourhood" => "neighbourhood.md"
        ],
        "Extension" => Any[
            "Update function" => "own_updatefunction.md",
            "Neighbourhood" => "own_neighbourhood.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/control13/ParticleSwarmOptimizer.jl.git",
    target = "build",
    julia = "0.6",
    deps   = nothing,
    make   = nothing
)


text:
__precompile__()

module AstroTime

using EarthOrientation
using Reexport
import RemoteFiles: path, isfile

export @timescale

include("utils.jl")

include("LeapSeconds.jl")
include("TimeScales.jl")
include("Periods.jl")
include("Epochs.jl")

@reexport using .TimeScales
@reexport using .Periods
@reexport using .LeapSeconds
@reexport using .Epochs

function __init__()
    isfile(LSK_FILE) && push!(LSK_DATA, path(LSK_FILE))
end

"""
    @timescale scale

Define a new timescale and the corresponding `Epoch` type alias.

# Example

```jldoctest
julia> @timescale Custom

julia> Custom <: TimeScale
true
julia> CustomEpoch == Epoch{Custom, T} where T <: Number
true
```
"""
macro timescale(scale)
    if !(scale isa Symbol)
        error("Invalid time scale name.")
    end
    epoch = Symbol(scale, "Epoch")
    return quote
        struct $(esc(scale)) <: TimeScale end
        const $(esc(epoch)) = Epoch{$(esc(scale))}
        nothing
    end
end

function update()
    EarthOrientation.update()
    download(LSK_FILE)
    push!(LSK_DATA, path(LSK_FILE))
    nothing
end

end # module


text:
function url(pkg)
    path = joinpath(Pkg.dir(), "METADATA", pkg, "url")
    ispath(path) || error("No METADATA for $pkg.")
    url_git = strip(readstring(path))
end

isrecipe(x) = false
function isrecipe(ex::Expr)
    (ex.head == :macrocall) &&
    (ex.args[1] == Symbol("@recipe"))
end

isline(ex) = false
isline(ex::Expr) = ex.head == :line


text:
using DiffEqFlux, DifferentialEquations, Optim, Flux, Zygote, Test, BenchmarkTools

u0 = Float32[2.; 0.]
datasize = 30
tspan = (0.0f0,1.5f0)

function trueODEfunc(du,u,p,t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end
t = range(tspan[1],tspan[2],length=datasize)
prob = ODEProblem(trueODEfunc,u0,tspan)
ode_data = Array(solve(prob,ROCK4(),saveat=t))

fastdudt2 = FastChain((x,p) -> x.^3,
            FastDense(2,50,tanh),
            FastDense(50, 2))

fast_n_ode = NeuralODE(fastdudt2,tspan,ROCK4(),saveat=t)

function fast_predict_n_ode(p)
  fast_n_ode(u0,p)
end

function fast_loss_n_ode(p)
    pred = fast_predict_n_ode(p)
    loss = sum(abs2,ode_data .- pred)
    loss,pred
end

dudt2 = Chain(x -> x.^3,
            Dense(2,50,tanh),
            Dense(50, 2))

n_ode = NeuralODE(dudt2,tspan,ROCK4(),saveat=t)

function predict_n_ode(p)
  n_ode(u0,p)
end

function loss_n_ode(p)
    pred = predict_n_ode(p)
    loss = sum(abs2,ode_data .- pred)
    loss,pred
end

p = fast_n_ode.p
_p,re = Flux.destructure(dudt2)

@test fastdudt2(ones(2),_p) ≈ dudt2(ones(2))
@test fast_loss_n_ode(p)[1] ≈ loss_n_ode(p)[1]
@test Zygote.gradient((p)->fast_loss_n_ode(p)[1], p)[1] ≈ Zygote.gradient((p)->loss_n_ode(p)[1], p)[1]

@btime Zygote.gradient((p)->fast_loss_n_ode(p)[1], p)
@btime Zygote.gradient((p)->fast_loss_n_ode(p)[1], p)
@btime Zygote.gradient((p)->loss_n_ode(p)[1], p)
@btime Zygote.gradient((p)->loss_n_ode(p)[1], p)

@benchmark Zygote.gradient((p)->fast_loss_n_ode(p)[1], p)
@benchmark Zygote.gradient((p)->loss_n_ode(p)[1], p)


text:
#===============================================================================
                                The main routine
===============================================================================#
"""
    mcmc(setup_mcmc::MCMCSetup, schedule::MCMCSchedule, setup)  <: ModelSetup

function for running the mcmc. receives as imput `MCMCSetup`, `MCMCschedule` and `setup`.
See [MCMCSetup](@ref), [MCMCSchedule](@ref). setup typically is DiffusionSetup making
additional setup choices when using the MCMC infrastructure to sample diffusion processes.
See [DiffusionSetup](@ref)
"""
function mcmc(setup_mcmc::MCMCSetup, schedule::MCMCSchedule, setup::T) where T <: ModelSetup
    ws, ll, θ = create_workspace(setup)
    ws_mcmc = create_workspace(setup_mcmc, schedule, θ)
    adpt = adaptation_object(setup, ws)

    aux = nothing
    for step in schedule
        step.save && save_imputed!(ws)
        for i in step.idx
            if step.param_updt || typeof(ws_mcmc.updates[i]) <: Imputation
                ws = next(ws, ws_mcmc.updates[i])
                ll, acc, θ = update!(ws_mcmc.updates[i], ws, θ, ll, step, aux)
                aux = aux_params(ws_mcmc.updates[i], aux)
                update!(ws_mcmc, acc, θ, i)
                step.verbose && print("\n")
            end
        end
        step.verbose && print("-----------------------------------------------",
                              "------\n")
        step.readjust && readjust!(ws_mcmc, step.iter)
        step.fuse && fuse!(ws_mcmc, schedule)
        ll = adaptation!(ws, adpt, step.iter, ll)
    end
    display_summary(ws, ws_mcmc)
    ws, ws_mcmc
end

display_summary(::Any, ::Any) = print("nothing to display...\n")
adaptation_object(::Any, ::Any) = nothing
adaptation!(::Any, ::Nothing, ::Any, ll) = ll


text:
using Mimi
using MimiDICE2013

DICE = MimiDICE2013.get_model()
run(DICE)

explore(DICE)


text:
abstract type AbstractVanillaWidget{T} <: AbstractObservable{T} end
Observables.observe(a::AbstractVanillaWidget) = observe(value(a))
Base.getproperty(a::AbstractVanillaWidget, s::Symbol) = getproperty(props(a), s)
Base.propertynames(a::AbstractVanillaWidget) = propertynames(props(a))

_get_value(x) = x
_get_value(x::AbstractObservable) = x[]
_get_value(x::AbstractDict) = Dict(key => _get_value(val) for (key, val) in x)
_get_value(x::NamedTuple) = map(_get_value, x)

to_observable(o::Observable) = o
to_observable(o::AbstractObservable) = Observables.observe(o)
to_observable(o) = Observable(o)

const specialprops = (:style, :attributes, :events)

observify((key, val)::Pair) = key in specialprops ? (key => to_observable(val)) : (key => val)

function addreactivity!(scp::Scope, tag, children...; addevents = Dict(),
    id = string("reactivenode-", uuid4()), properties...)

    n = node(tag, children...; id = id, map(_get_value, values(properties))...)

    onmount(scp, js"""
    function() {
        var node = _webIOScope.dom.querySelector('#'+$id)
        var events = $addevents;
        for (var key of Object.keys(events)) {
            node.addEventListener(key, events[key]);
        }
    }
    """)

    for (key, val) in properties
        key in specialprops && continue
        val isa AbstractObservable || continue
        skey = string(key)
        setobservable!(scp, skey, observe(val))
        onjs(scp[skey], js"""
        function (val) {
            var node = _webIOScope.dom.querySelector('#'+$id)
            node[$skey] = val;
        }
        """)
    end

    for propkey in specialprops
        dict = get(properties, propkey, Dict())
        for (key, val) in pairs(dict)
            val isa AbstractObservable || continue
            skey, spropkey = string(key), string(propkey)
            obskey = spropkey*"."*skey
            setobservable!(scp, obskey, observe(val))
            onjs(scp[obskey], js"""
            function (val) {
                var node = _webIOScope.dom.querySelector('#'+$id)
                node[$spropkey][$skey] = val;
            }
            """)
        end
    end
    scp.dom = n
    return
end


text:
# This file was generated, do not modify it. # hide
using MLJ
using PyPlot
ioff() # hide
import Statistics


text:
### PlotlyJS set up
import PlotlyJS
import WebIO

# WebIO.install_jupyter_nbextension()

function set_seriescolor(seriescolor::Array, gens::Array)
    colors = []
    for i in 1:length(gens)
        count = i % length(seriescolor)
        if count == 0
            count = length(seriescolor)
        end
        colors = vcat(colors, seriescolor[count])
    end
    return colors
end

function plotly_stack_gen(stacked_gen::StackedGeneration, seriescolor::Array; kwargs...)
    set_display = get(kwargs, :display, true)
    save_fig = get(kwargs, :save, nothing)
    traces = PlotlyJS.GenericTrace{Dict{Symbol, Any}}[]
    gens = stacked_gen.labels
    seriescolor = set_seriescolor(seriescolor, gens)
    for gen in 1:length(gens)
        push!(
            traces,
            PlotlyJS.scatter(;
                name = gens[gen],
                x = stacked_gen.time_range,
                y = stacked_gen.data_matrix[:, gen],
                stackgroup = "one",
                mode = "lines",
                fill = "tonexty",
                line_color = seriescolor[gen],
                fillcolor = seriescolor[gen],
            ),
        )
    end
    p = PlotlyJS.plot(
        traces,
        PlotlyJS.Layout(title = "Variables", yaxis_title = "Generation (MW)"),
    )
    set_display && PlotlyJS.display(p)
    if !isnothing(save_fig)
        PlotlyJS.savefig(p, joinpath(save_fig, "Stack_Generation.png"))
    end
end

function plotly_stack_gen(stacks::Array{StackedGeneration}, seriescolor::Array; kwargs...)
    set_display = get(kwargs, :display, true)
    save_fig = get(kwargs, :save, nothing)
    plots = []
    for stack in 1:length(stacks)
        trace = PlotlyJS.GenericTrace{Dict{Symbol, Any}}[]
        gens = stacks[stack].labels
        seriescolor = set_seriescolor(seriescolor, gens)
        for gen in 1:length(gens)
            if stack == 1
                leg = true
            else
                leg = false
            end
            push!(
                trace,
                PlotlyJS.scatter(;
                    name = gens[gen],
                    showlegend = leg,
                    x = stacks[stack].time_range,
                    y = stacks[stack].data_matrix[:, gen],
                    stackgroup = "one",
                    mode = "lines",
                    fill = "tonexty",
                    line_color = seriescolor[gen],
                    fillcolor = seriescolor[gen],
                ),
            )
        end
        p = PlotlyJS.plot(
            trace,
            PlotlyJS.Layout(title = "Variables", yaxis_title = "Generation (MW)"),
        )
        plots = vcat(plots, p)
    end
    plots = vcat(plots...)
    set_display && PlotlyJS.display(plots)
    if !isnothing(save_fig)
        PlotlyJS.savefig(plots, joinpath(save_fig, "Bar_Generation.png"))
    end
end

function _check_matching_variables(results::Array)
    variables = DataFrames.DataFrame()
    first_keys = collect(keys(results[1].variables))
    for res in 2:length(results)
        var = collect(keys(results[res].variables))
        if var != first_keys
            throw(IS.ConflictingInputsError("The given results do not have matching variable lists."))
        end
    end
end

function plotly_stack_plots(results::IS.Results, seriescolor::Array; kwargs...)
    set_display = get(kwargs, :display, true)
    save_fig = get(kwargs, :save, nothing)
    for (key, var) in results.variables
        traces = PlotlyJS.GenericTrace{Dict{Symbol, Any}}[]
        var = results.variables[key]
        gens = collect(names(var))
        seriescolor = set_seriescolor(seriescolor, gens)
        for gen in 1:length(gens)
            push!(
                traces,
                PlotlyJS.scatter(;
                    name = gens[gen],
                    x = results.time_stamp[:, 1],
                    y = convert(Matrix, var)[:, gen],
                    stackgroup = "one",
                    mode = "lines",
                    fill = "tonexty",
                    line_color = seriescolor[gen],
                    fillcolor = seriescolor[gen],
                ),
            )
        end
        p = PlotlyJS.plot(
            traces,
            PlotlyJS.Layout(title = "$key", yaxis_title = "Generation (MW)"),
        )
        set_display && PlotlyJS.display(p)
        if !isnothing(save_fig)
            Plots.savefig(p, joinpath(save_fig, "$(key)_Stack.png"))
        end
    end
end

function plotly_stack_plots(results::Array, seriescolor::Array; kwargs...)
    set_display = get(kwargs, :display, true)
    save_fig = get(kwargs, :save, nothing)
    _check_matching_variables(results)
    for key in collect(keys(results[1, 1].variables))
        plots = []
        for res in 1:size(results, 2)
            traces = PlotlyJS.GenericTrace{Dict{Symbol, Any}}[]
            var = results[1, res].variables[key]
            gens = collect(names(var))
            seriescolor = set_seriescolor(seriescolor, gens)
            for gen in 1:length(gens)
                if res == 1
                    leg = true
                else
                    leg = false
                end
                push!(
                    traces,
                    PlotlyJS.scatter(;
                        name = gens[gen],
                        showlegend = leg,
                        x = results[1, res].time_stamp[:, 1],
                        y = convert(Matrix, var)[:, gen],
                        stackgroup = "one",
                        mode = "lines",
                        fill = "tonexty",
                        line_color = seriescolor[gen],
                        fillcolor = seriescolor[gen],
                    ),
                )
            end
            p = PlotlyJS.plot(
                traces,
                PlotlyJS.Layout(
                    title = "$key",
                    yaxis_title = "Generation (MW)",
                    grid = (rows = 3, columns = 1, pattern = "independent"),
                ),
            )
            plots = vcat(plots, p)
        end
        plots = vcat(plots...)
        set_display && PlotlyJS.display(plots)
        if !isnothing(save_fig)
            PlotlyJS.savefig(plots, joinpath(save_fig, "Bar_Generation.png"))
        end
    end
end

function plotly_bar_gen(bar_gen::BarGeneration, seriescolor::Array; kwargs...)
    time_range = bar_gen.time_range
    set_display = get(kwargs, :display, true)
    save_fig = get(kwargs, :save, nothing)
    time_span = IS.convert_compound_period(
        convert(Dates.TimePeriod, time_range[2] - time_range[1]) * length(time_range),
    )
    traces = PlotlyJS.GenericTrace{Dict{Symbol, Any}}[]
    gens = bar_gen.labels
    seriescolor = set_seriescolor(seriescolor, gens)
    for gen in 1:length(gens)
        push!(
            traces,
            PlotlyJS.scatter(;
                name = gens[gen],
                x = ["$time_span, $(time_range[1])"],
                y = bar_gen.bar_data[:, gen],
                type = "bar",
                marker_color = seriescolor[gen],
            ),
        )
    end
    p = PlotlyJS.plot(
        traces,
        PlotlyJS.Layout(
            title = "Variables",
            yaxis_title = "Generation (MW)",
            color = seriescolor,
            barmode = "stack",
        ),
    )
    set_display && PlotlyJS.display(p)
    if !isnothing(save_fig)
        PlotlyJS.savefig(p, joinpath(save_fig, "Bar_Generation.png"))
    end
end

function plotly_bar_gen(bar_gen::Array{BarGeneration}, seriescolor::Array; kwargs...)
    time_range = bar_gen[1].time_range
    set_display = get(kwargs, :display, true)
    save_fig = get(kwargs, :save, nothing)
    time_span = IS.convert_compound_period(
        convert(Dates.TimePeriod, time_range[2] - time_range[1]) * length(time_range),
    )
    seriescolor = set_seriescolor(seriescolor, bar_gen[1].labels)
    plots = []
    for bar in 1:length(bar_gen)
        traces = PlotlyJS.GenericTrace{Dict{Symbol, Any}}[]
        gens = bar_gen[bar].labels
        for gen in 1:length(gens)
            if bar == 1
                leg = true
            else
                leg = false
            end
            push!(
                traces,
                PlotlyJS.scatter(;
                    name = gens[gen],
                    showlegend = leg,
                    x = ["$time_span, $(time_range[1])"],
                    y = bar_gen[bar].bar_data[:, gen],
                    type = "bar",
                    barmode = "stack",
                    marker_color = seriescolor[gen],
                ),
            )
        end
        p = PlotlyJS.plot(
            traces,
            PlotlyJS.Layout(
                title = "Variables",
                yaxis_title = "Generation (MW)",
                color = seriescolor,
                barmode = "stack",
                stackgroup = "one",
            ),
        )
        plots = vcat(plots, p)
    end
    plots = vcat(plots...)
    set_display && PlotlyJS.display(plots)
    if !isnothing(save_fig)
        PlotlyJS.savefig(plots, joinpath(save_fig, "Bar_Generation.png"))
    end
end

function plotly_bar_plots(results::Array, seriescolor::Array; kwargs...)
    set_display = get(kwargs, :display, true)
    save_fig = get(kwargs, :save, nothing)
    time_range = results[1].time_stamp
    time_span = IS.convert_compound_period(
        convert(Dates.TimePeriod, time_range[2, 1] - time_range[1, 1]) *
            size(time_range, 1),
    )
    for key in collect(keys(results[1].variables))
        plots = []
        for res in 1:length(results)
            var = results[res].variables[key]
            traces = PlotlyJS.GenericTrace{Dict{Symbol, Any}}[]
            gens = collect(names(var))
            seriescolor = set_seriescolor(seriescolor, gens)
            for gen in 1:length(gens)
                if res == 1
                    leg = true
                else
                    leg = false
                end
                push!(
                    traces,
                    PlotlyJS.scatter(;
                        name = gens[gen],
                        showlegend = leg,
                        x = ["$time_span, $(time_range[1, 1])"],
                        y = sum(convert(Matrix, var)[:, gen], dims = 1),
                        type = "bar",
                        barmode = "stack",
                        stackgroup = "one",
                        marker_color = seriescolor[gen],
                    ),
                )
            end
            p = PlotlyJS.plot(
                traces,
                PlotlyJS.Layout(
                    title = "$key",
                    yaxis_title = "Generation (MW)",
                    barmode = "stack",
                ),
            )
            plots = vcat(plots, p)
        end
        plots = vcat(plots...)
        set_display && PlotlyJS.display(plots)
        if !isnothing(save_fig)
            PlotlyJS.savefig(plots, joinpath(save_fig, "$(key)_Bar.png"))
        end
    end
end

function plotly_bar_plots(res::IS.Results, seriescolor::Array; kwargs...)
    set_display = get(kwargs, :display, true)
    save_fig = get(kwargs, :save, nothing)
    time_range = res.time_stamp
    time_span = IS.convert_compound_period(
        convert(Dates.TimePeriod, time_range[2, 1] - time_range[1, 1]) *
            size(time_range, 1),
    )
    for (key, var) in res.variables
        traces = PlotlyJS.GenericTrace{Dict{Symbol, Any}}[]
        gens = collect(names(var))
        seriescolor = set_seriescolor(seriescolor, gens)
        for gen in 1:length(gens)
            push!(
                traces,
                PlotlyJS.scatter(;
                    name = gens[gen],
                    x = ["$time_span, $(time_range[1, 1])"],
                    y = sum(convert(Matrix, var)[:, gen], dims = 1),
                    type = "bar",
                    barmode = "stack",
                    marker_color = seriescolor[gen],
                ),
            )
        end
        p = PlotlyJS.plot(
            traces,
            PlotlyJS.Layout(
                barmode = "stack",
                title = "$key",
                yaxis_title = "Generation (MW)",
            ),
        )
        set_display && PlotlyJS.display(p)
        if !isnothing(save_fig)
            PlotlyJS.savefig(p, joinpath(save_fig, "$(key)_Bar.png"))
        end
    end
end

RecipesBase.@recipe function StackedPlot(
    results::StackedArea,
    variable::String,
    seriescolor::Array,
)
    seriescolor := seriescolor
    time = convert.(Dates.DateTime, results.time_range)
    n = length(time)
    data = results.data_matrix
    z = cumsum(data, dims = 2)
    # Plot attributes
    grid := false
    title := variable
    if size(results.labels, 2) == 1 # work-around for a weird glitch in recipes
        label := results.labels[1]
    else
        label := results.labels
    end
    legend := :outerright
    interval = time[2] - time[1]
    time_interval = convert(Dates.Hour, interval * n)
    xlabel := "$time_interval"
    ylabel := "Generation (MW)"
    xtick := [time[1], time[n]]
    # create filled polygon
    sy = vcat(z[:, 1], zeros(n))
    sx = [time[1:n]; reverse(time[1:n])]
    for c in 1:size(z, 2)
        if c !== 1
            sy = hcat(sy, vcat(z[:, c], reverse(z[:, c - 1])))
        end
    end
    RecipesBase.@series begin
        seriestype := :shape
        sx, sy
    end
end

RecipesBase.@recipe function StackedPlot(
    results::Array{StackedArea},
    variable::String,
    seriescolor::Array,
)
    layout := (length(results), 1)
    for i in 1:length(results)
        res = results[i]
        time = convert.(Dates.DateTime, res.time_range)
        n = length(time)
        data = res.data_matrix
        z = cumsum(data, dims = 2)
        # Plot attributes
        grid := false
        title := variable
        if size(res.labels, 2) == 1 # work-around for a weird glitch in recipes
            label := res.labels[1]
        else
            label := res.labels
        end
        legend := :outerright
        interval = time[2] - time[1]
        time_interval = convert(Dates.Hour, interval * n)
        xlabel := "$time_interval"
        ylabel := "Generation (MW)"
        xtick := [time[1], time[n]]
        # create filled polygon
        sy = vcat(z[:, 1], zeros(n))
        sx = [time[1:n]; reverse(time[1:n])]
        for c in 1:size(z, 2)
            if c !== 1
                sy = hcat(sy, vcat(z[:, c], reverse(z[:, c - 1])))
            end
        end
        RecipesBase.@series begin
            subplot := i
            seriestype := :shape
            seriescolor := seriescolor[:, 1:size(sy, 2)]
            sx, sy
        end
    end
end

RecipesBase.@recipe function StackedGeneration(res::StackedGeneration, seriescolor::Array)
    seriescolor := seriescolor
    time = convert.(Dates.DateTime, res.time_range)
    n = length(time)
    data = res.data_matrix
    z = cumsum(data, dims = 2)
    # Plot Attributes
    grid := false
    title := "Generator"
    if size(res.labels, 2) == 1 # work-around for a weird glitch in recipes
        label := res.labels[1]
    else
        label := res.labels
    end
    legend := :outerright
    interval = time[2] - time[1]
    time_interval = convert(Dates.Hour, interval * n)
    xlabel := "$time_interval"
    ylabel := "Generation (MW)"
    xtick := [time[1], time[n]]
    # Create filled polygon
    sy = vcat(z[:, 1], zeros(n))
    sx = [time[1:n]; reverse(time[1:n])]
    for c in 2:size(z, 2)
        if c !== 1
            sy = hcat(sy, vcat(z[:, c], reverse(z[:, c - 1])))
        end
    end
    RecipesBase.@series begin
        seriestype := :shape
        sx, sy
    end

end

RecipesBase.@recipe function StackedGeneration(
    results::Array{StackedGeneration},
    seriescolor::Array,
)
    layout := (length(results), 1)
    for i in 1:length(results)
        res = results[i]
        time = convert.(Dates.DateTime, res.time_range)
        n = length(time)
        data = res.data_matrix
        z = cumsum(data, dims = 2)
        # Plot Attributes
        grid := false
        title := "Generator"
        if size(res.labels, 2) == 1 # work-around for a weird glitch in recipes
            label := res.labels[1]
        else
            label := res.labels
        end
        legend := :outerright
        interval = time[2] - time[1]
        time_interval = convert(Dates.Hour, interval * n)
        xlabel := "$time_interval"
        ylabel := "Generation (MW)"
        xtick := [time[1], time[n]]
        # Create filled polygon
        sy = vcat(z[:, 1], zeros(n))
        sx = [time[1:n]; reverse(time[1:n])]
        for c in 2:size(z, 2)
            if c !== 1
                sy = hcat(sy, vcat(z[:, c], reverse(z[:, c - 1])))
            end
        end
        RecipesBase.@series begin
            subplot := i
            seriestype := :shape
            seriescolor := seriescolor[:, 1:size(sy, 2)]
            sx, sy
        end
    end
end

RecipesBase.@recipe function BarPlot(res::BarPlot, variable::String, seriescolor::Array)
    seriescolor := seriescolor
    time = convert.(Dates.DateTime, res.time_range)
    n = length(time)
    data_point = res.bar_data
    data = [data_point; data_point]
    z = cumsum(data, dims = 2)
    # Plot Attributes
    grid := false
    title := variable
    seriestype := :shape
    if size(res.labels, 2) == 1 # work-around for a weird glitch in recipes
        label := res.labels[1]
    else
        label := res.labels
    end

    legend := :outerright
    interval = time[2] - time[1]
    time_interval = convert(Dates.Hour, interval * n)
    xlabel := "$time_interval, $(time[1])"
    ylabel := "Generation(MW)"
    xlims := (1, 8)
    xticks := false
    n = 2
    # Create filled polygon
    for c in 1:size(z, 2)
        sx = [[4, 5]; [5, 4]]
        sy = vcat(z[:, c], c == 1 ? zeros(n) : reverse(z[:, c - 1]))
        RecipesBase.@series sx, sy
    end

end

RecipesBase.@recipe function BarPlot(
    results::Array{BarPlot},
    variable::String,
    seriescolor::Array,
)
    layout := (length(results), 1)
    for i in 1:length(results)
        res = results[i]
        time = convert.(Dates.DateTime, res.time_range)
        n = length(time)
        data_point = res.bar_data
        data = [data_point; data_point]
        z = cumsum(data, dims = 2)
        # Plot Attributes
        grid := false
        title := variable
        seriestype := :shape
        if size(res.labels, 2) == 1 # work-around for a weird glitch in recipes
            label := res.labels[1]
        else
            label := res.labels
        end

        legend := :outerright
        interval = time[2] - time[1]
        time_interval = convert(Dates.Hour, interval * n)
        xlabel := "$time_interval, $(time[1])"
        ylabel := "Generation(MW)"
        xlims := (1, 8)
        xticks := false
        n = 2
        # Create filled polygon
        seriescolor := seriescolor[:, 1:size(data_point, 2)]
        for c in 1:size(z, 2)
            subplot := i
            sx = [[4, 5]; [5, 4]]
            sy = vcat(z[:, c], c == 1 ? zeros(n) : reverse(z[:, c - 1]))
            RecipesBase.@series sx, sy
        end
    end
end

RecipesBase.@recipe function BarGen(res::BarGeneration, seriescolor::Array)
    seriescolor := seriescolor
    time = convert.(Dates.DateTime, res.time_range)
    n = 2
    data_point = res.bar_data
    data = [data_point; data_point]
    z = cumsum(data, dims = 2)
    # Plot Attributes
    grid := false
    title := "Generator"
    start_time = time[1]
    interval = time[2] - time[1]
    time_interval = convert(Dates.Hour, interval * length(time))
    xlabel := "$time_interval, $(time[1])"
    ylabel := "Generation(MW)"
    seriestype := :shape
    if size(res.labels, 2) == 1 # work-around for a weird glitch in recipes
        label := res.labels[1]
    else
        label := res.labels
    end
    legend := :outerright
    xlims := (1, 8)
    xticks := false
    for c in 1:size(z, 2)
        sx = [[4, 5]; [5, 4]]
        sy = vcat(z[:, c], c == 1 ? zeros(n) : reverse(z[:, c - 1]))
        RecipesBase.@series sx, sy
    end
end

RecipesBase.@recipe function BarGen(results::Array{BarGeneration}, seriescolor::Array)
    layout := (length(results), 1)
    for i in 1:length(results)
        seriescolor := seriescolor
        res = results[i]
        time = convert.(Dates.DateTime, res.time_range)
        n = 2
        data_point = res.bar_data
        data = [data_point; data_point]
        z = cumsum(data, dims = 2)
        # Plot Attributes
        grid := false
        title := "Generator"
        start_time = time[1]
        interval = time[2] - time[1]
        time_interval = convert(Dates.Hour, interval * length(time))
        xlabel := "$time_interval, $(time[1])"
        ylabel := "Generation(MW)"
        seriestype := :shape
        if size(res.labels, 2) == 1 # work-around for a weird glitch in recipes
            label := res.labels[1]
        else
            label := res.labels
        end
        legend := :outerright
        xlims := (1, 8)
        xticks := false
        seriescolor := seriescolor[:, 1:size(data_point, 2)]
        for c in 1:size(z, 2)
            subplot := i
            sx = [[4, 5]; [5, 4]]
            sy = vcat(z[:, c], c == 1 ? zeros(n) : reverse(z[:, c - 1]))
            RecipesBase.@series sx, sy
        end
    end
end


text:
using Documenter
using FileType

makedocs(
    sitename = "FileType",
    format = Documenter.HTML(),
    modules = [FileType]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "https://github.com/Agent-Hellboy/FileType.jl.git"
)


text:
import JUDI.split, JUDI.subsample

export H5SeisCon

mutable struct H5SeisCon
  seis::PyCall.PyObject
  pkey::String
  # 'pkeyvals' must be Float64 as roundoff may cause problems
  # when retrieving data from 'seis' (getPKeyTraceSize for example)
  pkeyvals::Array{Float64,1}
end

function H5SeisCon(;
  seis::PyCall.PyObject=nothing,
  pkey::String="SP",
  pkeyvals=nothing)
  if isnothing(seis)
    @error "Seis is Nothing\n"
    return
  end

  if isnothing(pkey)
    @error "pkey is Nothing\n"
    return
  end

  if !seis.hasPKeySort(pkey)
    @info "Updating trace header limits and adding pkey sort: $pkey..."
    seis.updateTraceHeaderLimits(Int(1e7))
    seis.addPKeySort(pkey)
  end

  # must go after 'addPKeySort'
  if isnothing(pkeyvals)
    pkeyvals = seis.getPKeyValues(pkey)
  end

  return H5SeisCon(seis, pkey, pkeyvals)
end

size(con::H5SeisCon) = size(con.pkeyvals)
length(con::H5SeisCon) = length(con.pkeyvals)

function getindex(con::H5SeisCon, a::TA) where {TA<:Union{Array{<:Integer,1}, AbstractRange, Integer}}
  read_H5SeisCon(con, a)
end

function getindex(con::H5SeisCon, a::Colon)
  read_H5SeisCon(con, 1:length(con))
end

function split(s::H5SeisCon, inds::Union{Vector{Ti}, AbstractRange{Ti}}) where {Ti<:Integer}
  return H5SeisCon(s.seis, s.pkey, view(s.pkeyvals, inds))
end

split(s::H5SeisCon, inds::Integer) = split(s, [inds])

subsample(a::Array{H5SeisCon,1}, i::Int64) = a[i]

iterate(S::H5SeisCon, state::Integer=1) = state > length(S) ? nothing : (S[state], state+1)


text:
#
# Copyright (c) 2008, Maxim Likhachev
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
#     * Redistributions of source code must retain the above copyright
#       notice, this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above copyright
#       notice, this list of conditions and the following disclaimer in the
#       documentation and/or other materials provided with the distribution.
#     * Neither the name of the Carnegie Mellon University nor the names of its
#       contributors may be used to endorse or promote products derived from
#       this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES LOSS OF USE, DATA, OR PROFITS OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
#

# Julia translation.

using Plots
using Printf
using LinearAlgebra
default(legend = false)

# Play with these based on your needs
# ------------ From here ------------
resolution = 0.05
numberofangles = 16 # preferably a power of 2, definitely multiple of 8
numberofprimsperangle = 5

# Multipliers (multiplier is used as costmult*cost)
forwardcostmult = 1
backwardcostmult = 50
forwardandturncostmult = 2
sidestepcostmult = 50
turninplacecostmult = 50

# ------------- To here -------------

function genmprim_unicycle(outfilename)
    # 0 degreees
    basemprimendpts0_c = zeros(numberofprimsperangle, 4) # x,y,theta,costmult
    # x aligned with the heading of the robot, angles are positive
    # counterclockwise
    # 0 theta change
    basemprimendpts0_c[1,:] = [1,0,0,forwardcostmult]
    basemprimendpts0_c[2,:] = [8,0,0,forwardcostmult]
    basemprimendpts0_c[3,:] = [-1,0,0,backwardcostmult]
    # 1/16 theta change
    basemprimendpts0_c[4,:] = [8,1,1,forwardandturncostmult]
    basemprimendpts0_c[5,:] = [8,-1,-1,forwardandturncostmult]
    # turn in place
    # basemprimendpts0_c[6,:] = [0,0,1,turninplacecostmult]
    # basemprimendpts0_c[7,:] = [0,0,-1,turninplacecostmult]

    # 45 degrees
    basemprimendpts45_c = zeros(numberofprimsperangle, 4) # x,y,theta,costmult (multiplier is used as costmult*cost)
    # x aligned with the heading of the robot, angles are positive
    # counterclockwise
    # 0 theta change
    basemprimendpts45_c[1,:] = [1,1,0,forwardcostmult]
    basemprimendpts45_c[2,:] = [6,6,0,forwardcostmult]
    basemprimendpts45_c[3,:] = [-1,-1,0,backwardcostmult]
    # 1/16 theta change
    basemprimendpts45_c[4,:] = [5,7,1,forwardandturncostmult]
    basemprimendpts45_c[5,:] = [7,5,-1,forwardandturncostmult]
    # turn in place
    # basemprimendpts45_c[6,:] = [0,0,1,turninplacecostmult]
    # basemprimendpts45_c[7,:] = [0,0,-1,turninplacecostmult]

    # 22.5 degrees
    basemprimendpts22p5_c = zeros(numberofprimsperangle, 4) # x,y,theta,costmult (multiplier is used as costmult*cost)
    # x aligned with the heading of the robot, angles are positive
    # counterclockwise
    # 0 theta change
    basemprimendpts22p5_c[1,:] = [2,1,0,forwardcostmult]
    basemprimendpts22p5_c[2,:] = [6,3,0,forwardcostmult]
    basemprimendpts22p5_c[3,:] = [-2,-1,0,backwardcostmult]
    # 1/16 theta change
    basemprimendpts22p5_c[4,:] = [5,4,1,forwardandturncostmult]
    basemprimendpts22p5_c[5,:] = [7,2,-1,forwardandturncostmult]
    # turn in place
    # basemprimendpts22p5_c[6,:] = [0,0,1,turninplacecostmult]
    # basemprimendpts22p5_c[7,:] = [0,0,-1,turninplacecostmult]

    open(outfilename, "w") do f

        @printf(f, "resolution_m: %f\n", resolution)
        @printf(f, "numberofangles: %d\n", numberofangles)
        @printf(f, "totalnumberofprimitives: %d\n", numberofprimsperangle*numberofangles)

        for angleind in 1:numberofangles
            #text(----) TODO
            display(plot())

            for primind in 1:numberofprimsperangle
                @printf(f, "primID: %d\n", primind-1)
                @printf(f, "startangle_c: %d\n", angleind-1)

                # current angle
                currentangle = (angleind-1)*2*pi/numberofangles
                currentangle_36000int = round((angleind-1)*36000/numberofangles)

                # compute which template to use
                if (currentangle_36000int %  9000 == 0)
                    basemprimendpts_c = basemprimendpts0_c[primind,:]
                    angle = currentangle
                elseif (currentangle_36000int % 4500 == 0)
                    basemprimendpts_c = basemprimendpts45_c[primind,:]
                    angle = currentangle - 45*pi/180
                elseif ((currentangle_36000int-7875) % 9000 == 0)
                    basemprimendpts_c = basemprimendpts33p75_c[primind,:]
                    basemprimendpts_c[1] = basemprimendpts33p75_c[primind, 2] # reverse x and y
                    basemprimendpts_c[2] = basemprimendpts33p75_c[primind, 1]
                    basemprimendpts_c[3] = -basemprimendpts33p75_c[primind, 3] # reverse the angle as well
                    angle = currentangle - 78.75*pi/180
                    println("78p75")
                elseif ((currentangle_36000int-6750) % 9000 == 0)
                    basemprimendpts_c = basemprimendpts22p5_c[primind,:]
                    basemprimendpts_c[1] = basemprimendpts22p5_c[primind, 2] # reverse x and y
                    basemprimendpts_c[2] = basemprimendpts22p5_c[primind, 1]
                    basemprimendpts_c[3] = -basemprimendpts22p5_c[primind, 3] # reverse the angle as well
                    angle = currentangle - 67.5*pi/180
                    println("67p5")
                elseif ((currentangle_36000int-5625) % 9000 == 0)
                    basemprimendpts_c = basemprimendpts11p25_c[primind,:]
                    basemprimendpts_c[1] = basemprimendpts11p25_c[primind, 2] # reverse x and y
                    basemprimendpts_c[2] = basemprimendpts11p25_c[primind, 1]
                    basemprimendpts_c[3] = -basemprimendpts11p25_c[primind, 3] # reverse the angle as well
                    angle = currentangle - 56.25*pi/180
                    println("56p25")
                elseif ((currentangle_36000int-3375) % 9000 == 0)
                    basemprimendpts_c = basemprimendpts33p75_c[primind,:]
                    angle = currentangle - 33.75*pi/180
                    println("33p75")
                elseif ((currentangle_36000int-2250) % 9000 == 0)
                    basemprimendpts_c = basemprimendpts22p5_c[primind,:]
                    angle = currentangle - 22.5*pi/180
                    println("22p5")
                elseif ((currentangle_36000int-1125) % 9000 == 0)
                    basemprimendpts_c = basemprimendpts11p25_c[primind,:]
                    angle = currentangle - 11.25*pi/180
                    println("11p25")
                else
                    println("ERROR: invalid angular resolution. angle = $(currentangle_36000int)")
                    return
                end

                # now figure out what action will be
                baseendpose_c = basemprimendpts_c[1:3]
                additionalactioncostmult = basemprimendpts_c[4]
                endx_c = round(baseendpose_c[1]*cos(angle) - baseendpose_c[2]*sin(angle))
                endy_c = round(baseendpose_c[1]*sin(angle) + baseendpose_c[2]*cos(angle))
                endtheta_c = (angleind - 1 + baseendpose_c[3]) % numberofangles
                endpose_c = [endx_c endy_c endtheta_c]

                println("rotation angle=$(angle*180/pi)")

                if baseendpose_c[2] == 0 && baseendpose_c[3] == 0
                    # println("endpose=$(endpose_c[1]) $(endpose_c[2]) $(endpose_c[3])")
                end

                # generate intermediate poses (remember they are w.r.t 0,0 (and not
                # centers of the cells)
                numofsamples = 10
                intermcells_m = zeros(numofsamples,3)
                startpt = [0,0,currentangle]
                endpt = [endpose_c[1]*resolution,endpose_c[2]*resolution,((angleind - 1 + baseendpose_c[3]) % numberofangles)*2*pi/numberofangles]
                intermcells_m = zeros(numofsamples,3)
                if ((endx_c == 0 && endy_c == 0) || baseendpose_c[3] == 0) # turn in place or move forward
                    for iind in 1:numofsamples
                        intermcells_m[iind,:] = [startpt[1] + (endpt[1] - startpt[1])*(iind-1)/(numofsamples-1),startpt[2] + (endpt[2] - startpt[2])*(iind-1)/(numofsamples-1),0]
                        rotation_angle = baseendpose_c[3] * (2*pi/numberofangles)
                        intermcells_m[iind,3] = (startpt[3] + (rotation_angle)*(iind-1)/(numofsamples-1)) % (2*pi)
                    end
                else # unicycle-based move forward or backward
                    R = [cos(startpt[3]) sin(endpt[3])-sin(startpt[3]);sin(startpt[3]) -(cos(endpt[3]) -cos(startpt[3]))]
                    S = pinv(R)*[endpt[1] - startpt[1],endpt[2] - startpt[2]]
                    l = S[1]
                    tvoverrv = S[2]
                    rv = (baseendpose_c[3]*2*pi/numberofangles + l/tvoverrv)
                    tv = tvoverrv*rv

                    if l < 0
                        println("WARNING: l = $l < 0 -> bad action start/end points")
                        l = 0
                    end
                    # compute rv
                    # rv = baseendpose_c[3]*2*pi/numberofangles
                    # compute tv
                    # tvx = (endpt[1] - startpt[1])*rv/(sin(endpt[3]) - sin(startpt[3]))
                    # tvy = -(endpt[2] - startpt[2])*rv/(cos(endpt[3]) - cos(startpt[3]))
                    # tv = (tvx + tvy)/2.0
                    # generate samples
                    for iind = 1:numofsamples
                        dt = (iind-1)/(numofsamples-1)

                        if dt*tv < l
                            intermcells_m[iind,:] = [startpt[1] + dt*tv*cos(startpt[3]),startpt[2] + dt*tv*sin(startpt[3]),startpt[3]]
                        else
                            dtheta = rv*(dt - l/tv) + startpt[3]
                            intermcells_m[iind,:] = [startpt[1] + l*cos(startpt[3]) + tvoverrv*(sin(dtheta) - sin(startpt[3])),startpt[2] + l*sin(startpt[3]) - tvoverrv*(cos(dtheta) - cos(startpt[3])),dtheta]
                        end
                    end
                    # correct
                    errorxy = [endpt[1] - intermcells_m[numofsamples,1]
                               endpt[2] - intermcells_m[numofsamples,2]]
                    println("l=$l errx=$(errorxy[1]) erry=$(errorxy[2])")
                    interpfactor = collect(0:1/(numofsamples-1):1)
                    intermcells_m[:,1] = intermcells_m[:,1] + errorxy[1]*interpfactor
                    intermcells_m[:,2] = intermcells_m[:,2] + errorxy[2]*interpfactor
                end
                # write out
                @printf(f, "endpose_c: %d %d %d\n", endpose_c[1], endpose_c[2], endpose_c[3])
                @printf(f, "additionalactioncostmult: %d\n", additionalactioncostmult)
                @printf(f, "intermediateposes: %d\n", size(intermcells_m,1))
                for interind in 1:size(intermcells_m, 1)
                    @printf(f, "%.4f %.4f %.4f\n", intermcells_m[interind,1], intermcells_m[interind,2], intermcells_m[interind,3]);
                end
                display(plot!(intermcells_m[:,1], intermcells_m[:,2], xlim=(-0.3,0.3), ylim=(-0.3,0.3)))
            end
            readline()
        end

    end
end

if length(ARGS) < 1
    println("Ooooops!")
    exit()
end

genmprim_unicycle(ARGS[1])


text:
module TestResult

import Wordle
import Test: @test, @test_throws

@test Wordle.Result("⬜🟨🟩⬜🟨") === Wordle.Result(
    (
        Wordle.PRESENT_NOWHERE,
        Wordle.PRESENT_ELSEWHERE,
        Wordle.PRESENT_HERE,
        Wordle.PRESENT_NOWHERE,
        Wordle.PRESENT_ELSEWHERE,
    )
)

@test Wordle.Result("⬜🟨⬜🟨🟩") === Wordle.Result(
    (
        Wordle.PRESENT_NOWHERE,
        Wordle.PRESENT_ELSEWHERE,
        Wordle.PRESENT_NOWHERE,
        Wordle.PRESENT_ELSEWHERE,
        Wordle.PRESENT_HERE,
    )
)

end # module


text:
# General math functions of quaternions

"""
    abs2(q)

Sum the squares of the components of the quaternion

# Examples
```jldoctest
julia> abs2(Quaternion(1,2,4,10))
121
```
"""
Base.abs2(q::AbstractQuaternion) = sum(q.components.^2)
Base.abs2(q::Rotor{T}) where {T<:Real} = one(T)

"""
    abs(q)

Square-root of the sum the squares of the components of the quaternion

# Examples
```jldoctest
julia> abs(Quaternion(1,2,4,10))
11.0
```
"""
Base.abs(q::AbstractQuaternion) = sqrt(abs2(q))
Base.abs(q::Rotor{T}) where {T<:Real} = one(T)

"""
    abs2vec(q)

Sum the squares of the "vector" components of the quaternion

# Examples
```jldoctest
julia> abs2vec(Quaternion(1,2,3,6))
49
```
"""
abs2vec(q::AbstractQuaternion) = @inbounds q.components[2]^2 + q.components[3]^2 + q.components[4]^2

"""
    absvec(q)

Square-root of the sum of the squares of the "vector" components of the quaternion

# Examples
```jldoctest
julia> absvec(Quaternion(1,2,3,6))
7.0
```
"""
absvec(q::AbstractQuaternion) = sqrt(abs2vec(q))

# norm(q::Quaternion) = Base.abs2(q)  ## This might just be confusing

Base.inv(q::AbstractQuaternion) = conj(q) / abs2(q)
Base.inv(q::Rotor) = conj(q)  # Specialize to ensure output is also a Rotor

"""
    log(q)

Logarithm of a quaternion.

As with the usual complex logarithm, the quaternion logarithm has multiple
branches, though the quaternion branches are three-dimensional: for any unit
"vector" quaternion q̂, you could add any integer multiple of 2πq̂ to the result
of this function and still get the same result after exponentiating (within
numerical accuracy).  This function is the principal logarithm.

This function has discontinuous (and fairly arbitrary) behavior along the
negative real axis: if the "vector" components of the quaternion are precisely
zero *and* the scalar component is negative, the returned quaternion will have
scalar component `log(-q.w)`, but will also have a `z` component of π.  The
choice of the `z` direction is arbitrary; the "vector" component of the
returned quaternion could be π times any unit vector.

Note that this function is not specialized to unit-quaternion inputs, so the
scalar component of the returned value will be nonzero unless the input has
*precisely* unit magnitude.

# Examples
```jldoctest
julia> log(exp(1.2imy))
0.0 + 0.0𝐢 + 1.2𝐣 + 0.0𝐤

julia> log(Quaternion(exp(7)))
7.0 + 0.0𝐢 + 0.0𝐣 + 0.0𝐤

julia> log(Quaternion(-exp(7)))
7.0 + 0.0𝐢 + 0.0𝐣 + 3.141592653589793𝐤
```

"""
function Base.log(q::Quaternion{T}) where {T}
    q = float(q)
    absolute2vec = abs2vec(q)
    if iszero(absolute2vec)
        if q.w < 0
            return Quaternion(log(-q.w), 0, 0, π)
        end
        return Quaternion(log(q.w), 0, 0, 0)
    end
    absolutevec = sqrt(absolute2vec)
    f = atan(absolutevec, q.w) / absolutevec  # acos((w^2-absolutevec^2) / (w^2+absolutevec^2)) / 2absolutevec
    Quaternion(log(abs2(q))/2, f*q.x, f*q.y, f*q.z)
end
function Base.log(q::Rotor{T}) where {T}
    q = float(q)
    absolute2vec = abs2vec(q)
    if iszero(absolute2vec)
        if q.w < 0
            return QuatVec{float(T)}(0, 0, 0, π)
        end
        return QuatVec{float(T)}(0, 0, 0, 0)
    end
    absolutevec = sqrt(absolute2vec)
    f = atan(absolutevec, q.w) / absolutevec  # acos(q.w) / absolutevec
    QuatVec(0, f*q.x, f*q.y, f*q.z)
end

"""
    exp(q)

Exponential of a quaternion

# Examples
```jldoctest
julia> exp(imx*π/4)  # Rotation through π/2 (note the extra 1/2) about the x axis
0.7071067811865476 + 0.7071067811865475𝐢 + 0.0𝐣 + 0.0𝐤
```
"""
function Base.exp(q::Quaternion{T}) where {T}
    q = float(q)
    absolute2vec = abs2vec(q)
    if iszero(absolute2vec)
        return Quaternion(exp(q.w), 0, 0, 0)
    end
    absolutevec = sqrt(absolute2vec)
    e = exp(q.w)
    s, c = sincos(absolutevec)
    esinc = e * s / absolutevec
    Quaternion(e*c, esinc*q.x, esinc*q.y, esinc*q.z)
end
function Base.exp(q::QuatVec{T}) where {T}
    q = float(q)
    absolute2vec = abs2vec(q)
    if iszero(absolute2vec)
        return Rotor{float(T)}(1, 0, 0, 0)
    end
    absolutevec = sqrt(absolute2vec)
    s, c = sincos(absolutevec)
    sinc = s / absolutevec
    Rotor(c, sinc*q.x, sinc*q.y, sinc*q.z)
end

@doc raw"""
    sqrt(q)

Square-root of a quaternion.

The general formula whenever the denominator is nonzero is

```math
\sqrt{q} = \frac{|q| + q} {\sqrt{2|q| + 2q.w}}
```

This can be proven by expanding `q` as `q.w + q.vec` and multiplying the
expression above by itself.

When the denominator is zero, this function has discontinuous (and fairly
arbitrary) behavior, just as with the quaternion [`log`](@ref) function.  In
this case, either all components are zero — in which case the result is simply
the zero quaternion — or the "vector" components of the quaternion are
precisely zero and the scalar component is negative.  If the latter is true,
the denominator above will be a pure-imaginary number.  Because the quaternions
come with infinitely many elements that square to -1, it is not clear *which*
imaginary should be used, so we arbitrarily choose to set the result
proportional to the `z` quaternion.  The choice of the `z` direction is
arbitrary; the "vector" component of the returned quaternion could be in any
direction.

# Examples
```jldoctest
julia> q = Quaternion(1.2, 3.4, 5.6, 7.8);

julia> sqrtq = √q;

julia> sqrtq^2 ≈ q
true

julia> √Quaternion(4)
2.0 + 0.0𝐢 + 0.0𝐣 + 0.0𝐤

julia> √Quaternion(-4)
0.0 + 0.0𝐢 + 0.0𝐣 + 2.0𝐤
```
"""
function Base.sqrt(q::Quaternion{T}) where {T}
    q = float(q)
    absolute2vec = abs2vec(q)
    if iszero(absolute2vec)
        if q.w < 0
            return Quaternion(0, 0, 0, sqrt(-q.w))
        end
        return Quaternion(sqrt(q.w), 0, 0, 0)
    end
    absolute2 = absolute2vec + q.w^2
    c1 = sqrt(absolute2) + q.w
    c2 = sqrt(inv(2*c1))
    Quaternion(c1*c2, q.x*c2, q.y*c2, q.z*c2)
end
function Base.sqrt(q::Rotor{T}) where {T}
    q = float(q)
    absolute2vec = abs2vec(q)
    if iszero(absolute2vec)
        if q.w < 0
            return Rotor{float(T)}(0, 0, 0, 1)
        end
        return Rotor{float(T)}(1, 0, 0, 0)
    end
    c1 = 1 + q.w
    c2 = sqrt(inv(2*c1))
    Rotor(c1*c2, q.x*c2, q.y*c2, q.z*c2)
end

"""
    angle(q)

Phase angle in radians of the rotation represented by this quaternion.

Note that this may be different from your interpretation of the angle of a
complex number in an important way.  Because quaternions act on vectors by
conjugation — as in `q*v*conj(q)` — there are *two* copies of `q` involved in
that expression; in some sense, a quaternion acts "twice".  Therefore, this
angle may be twice what you expect from an analogy with complex numbers —
dpending on how you interpret the correspondence between complex numbers and
quaternions.  Also, while rotations in the complex plane have a natural choice
of axis (the positive `z` direction), that is not the case for quaternions,
which means that the sign of this angle is arbitrary, and we always choose it
to be positive.

# Examples
```jldoctest
julia> θ=1.2;

julia> R=exp(θ * imz / 2);

julia> angle(R)
1.2

```
"""
Base.angle(q::Quaternion{T}) where T = 2 * absvec(log(q))
Base.angle(q::Rotor{T}) where T = 2 * absvec(log(q))

function Base.:^(q::Quaternion, s::Real)
    exp(s * log(q))
end
function Base.:^(q::Rotor, s::Real)
    q = float(q)
    absolutevec = absvec(q)
    if absolutevec ≤ eps(typeof(absolutevec))
        if q.w < 0
            # log(q) ≈ π𝐤
            sin_πs, cos_πs = sincospi(oftype(absolutevec, s))
            return Rotor{eltype(q)}([cos_πs, 0, 0, sin_πs])
        end
        # log(q) ≈ 0
        return one(q)
    end
    f1 = s * atan(absolutevec, q.w)
    sin_f1, cos_f1 = sincos(f1)
    f2 = sin_f1 / absolutevec
    Rotor{typeof(cos_f1)}([cos_f1, f2*q.x, f2*q.y, f2*q.z])
end
Base.:^(q::Quaternion, s::Integer) = (s ≥ 0 ? Base.power_by_squaring(q, s) : inv(Base.power_by_squaring(q, -s)))
Base.:^(q::QuatVec, s::Integer) = (s ≥ 0 ? Base.power_by_squaring(q, s) : inv(Base.power_by_squaring(q, -s)))
Base.:^(q::Rotor, s::Integer) = (s ≥ 0 ? Base.power_by_squaring(q, s) : inv(Base.power_by_squaring(q, -s)))


text:
# Slope Limiters
# note that r = (u - ul)/(ur - u)
# and we set a = u - ul, b = ur - u to define the limiters

function minmod(a,b,c)
  if (a > 0 && b > 0 && c > 0)
    min(a,b,c)
  elseif (a < 0 && b < 0 && c < 0)
    max(a,b,c)
  else
    zero(a)
  end
end

struct GeneralizedMinmodLimiter <: AbstractSlopeLimiter
    θ::Float64
end

GeneralizedMinmodLimiter(;θ=1.0) = GeneralizedMinmodLimiter(θ)
(limiter::GeneralizedMinmodLimiter)(a,b) =minmod(limiter.θ*a,0.5*(a+b),limiter.θ*b)

function minmod(a,b)
  0.5*(sign(a)+sign(b))*min(abs(a),abs(b))
end

struct MinmodLimiter <: AbstractSlopeLimiter end
(::MinmodLimiter)(a,b) = minmod(a,b)

struct OsherLimiter <: AbstractSlopeLimiter
    β::Float64
end
OsherLimiter(;β=1.0) = OsherLimiter(β)
(limiter::OsherLimiter)(a,b) = max(zero(a),min(a,limiter.β*b))

struct SuperbeeLimiter <: AbstractSlopeLimiter end
(::SuperbeeLimiter)(a,b) = max(zero(a),min(2*a,b),min(a,2*b))

"""
function compute_slopes(u, mesh, θ, M, ::Type{Val{true}})
  Estimate slopes of the discretization of function u,
      using a generalized minmod limiter
  inputs:
  `u` discrete approx of function u
  'nonscalar' bool
  `θ` parameter of generalized minmod limiter
  `mesh` problem mesh
  `Type{Val}` bool to choose threaded version
"""
function compute_slopes(u::AbstractArray, mesh::AbstractFVMesh, slopeLimiter::AbstractSlopeLimiter, nonscalar::Bool, ::Type{Val{true}})
  ∇u = similar(u)
  Threads.@threads for j in cell_indices(mesh)
      inner_slopes_loop!(∇u,j,u,mesh,slopeLimiter,nonscalar)
  end
  ∇u
end

function compute_slopes(u::AbstractArray, mesh::AbstractFVMesh, slopeLimiter::AbstractSlopeLimiter, nonscalar::Bool, ::Type{Val{false}})
  ∇u = similar(u)
  for j in cell_indices(mesh)
      inner_slopes_loop!(∇u,j,u,mesh,slopeLimiter,nonscalar)
  end
  ∇u
end

function inner_slopes_loop!(∇u,j,u,mesh,slopeLimiter::AbstractSlopeLimiter,nonscalar)
  ul = cellval_at_left(j,u,mesh)
  ur = cellval_at_right(j+1,u,mesh)
  if nonscalar
    @inbounds for i = 1:size(u,1)
      ∇u[i,j] = slopeLimiter(u[i,j]-ul[i], ur[i]-u[i,j])
    end
  else
    ∇u[j] = slopeLimiter(u[j]-ul, ur-u[j])
  end
end


text:
using Pkg

Pkg.add("PyPlot")
Pkg.add("DelimitedFiles")
Pkg.add("Random")
Pkg.add("StatsBase")
Pkg.add("CSV")
Pkg.add("DataFrames")
Pkg.add("SpecialFunctions")
Pkg.add("StatsFuns")
Pkg.add("Distributions")


text:
# Smile Tutorial 1: Creating a Bayesian Network
# this is the Julia version for the tutorial at
# https://dslpitt.org/genie/wiki/SMILE_Tutorial_1:_Creating_a_Bayesian_Network

using Smile

net = Network()
success = add_node(net, DSL_CPT, "success")
forecast = add_node(net, DSL_CPT, "forecast")

nodeA = get_node(net, success)
nodeDefA = definition(nodeA)
nodeValA = value(nodeA)

somenames = IdArray()
add(somenames, "Success")
add(somenames, "Failure")
set_number_of_outcomes( nodeDefA, somenames )

flush(somenames)
add(somenames, "Good")
add(somenames, "Moderate")
add(somenames, "Poor")
set_number_of_outcomes( definition(get_node(net, forecast)), somenames )

add_arc(net, success, forecast)

theprobs = DoubleArray()
set_size(theprobs, 2)
set_at(theprobs, 0, 0.2)
set_at(theprobs, 1, 1.8)
set_definition(definition(get_node(net, success)), theprobs)

thecoordinates = SysCoordinates(definition(get_node(net, forecast)))

write_file(net, "output.xdsl")


text:
macro show(ex)
    :(println(($string(ex)), " = ", $ex) )
end

const doublecolon = @eval (:(::Int)).head

quoted_expr(ex) = expr(:quote, ex)
quoted_exprs(exprs) = { quoted_expr(ex) for ex in exprs }

quoted_tuple(t) = expr(:tuple, {t...})

is_expr(ex, head::Symbol) = (isa(ex, Expr) && (ex.head == head))
function is_expr(ex, head::Symbol, nargs::Int)
    is_expr(ex, head) && length(ex.args) == nargs
end

# == @expect ==================================================================

fail_expect(predexpr) = error("Expected: ", string(predexpr))

default_code_checkexpect(ex) = :($ex ? nothing : fail_expect($quoted_expr(ex)))
code_checkexpect(ex)         = default_code_checkexpect(ex)

macro expect(args...)
    if !(1 <= length(args) <= 2)
        error("\n@expect: expected one or two arguments")
    end
    predexpr = args[1]
    if length(args) >= 2
        # explicit error message given
        message = args[2]
        msg_parts = is_expr(message, :tuple) ? message.args : {message}
        return :(($predexpr) ? nothing : ($expr(:call, :error, msg_parts...)))
    else # no explicit message
        return code_checkexpect(predexpr)
    end
end

# == @setdefault ==============================================================

# macro setdefault(args...)
#     # allow @setdefault(refexpr, default)
#     if (length(args)==1) && is_expr(args[1], :tuple)
#         args = args[1].args
#     end
#     refexpr, default = tuple(args...)
macro setdefault(refexpr, default)
    @expect is_expr(refexpr, :ref)
    dict_expr, key_expr = tuple(refexpr.args...)
    @gensym dict key #defval
    quote
        ($dict)::Associative = ($dict_expr)
        ($key) = ($key_expr)
        if has(($dict), ($key))
            ($dict)[($key)]
        else
            ($dict)[($key)] = ($default) # returns the newly inserted value
        end
    end
end

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
struct GibbsWithGradients{T3, V, B3, CS} <: SpinSampler{T3}
    A::T3           # q × N × M
    A₁::T3          # 1 × N × M
    α::T3           # 1 × 1 × M
    I::V            # M element vector
    batchshift::V   # M element vector
    𝕡::T3           # q × N × M
    𝕀::B3           # q × N × M with Bool entries
    q::Int          # number of spin states
    N::Int          # number of spins
    M::Int          # sample size
    catsampler::CS
    function GibbsWithGradients(X::AbstractArray{<:AbstractFloat,3})
        q, N, M = size(X)
        _eltype = promote_type(eltype(X), Float32)
        A = similar(X, _eltype, q, N, M)
        A₁= similar(X, _eltype, 1, N, M)
        α = similar(X, _eltype, 1, 1, M)
        I = similar(X, Int, M)
        batchshift = oftype(I, range(0; step = q*N, length = M))
        𝕡 = similar(X, _eltype, q, N, M)
        𝕀 = similar(X, Bool, q, N, M)
        catsampler = InverseCDF(𝕡, q*N, M)
        return new{typeof(A), typeof(I), typeof(𝕀), typeof(catsampler)}(
            A, A₁, α, I, batchshift, 𝕡, 𝕀, q, N, M, catsampler
        )
    end
end

(f::GibbsWithGradients)(X, energy) = gibbswithgradients!(X, energy, f)

function gibbswithgradients!(X::AbstractArray{T,3}, energy, buffer = GibbsWithGradients(X)) where T
    (; A, A₁, α, I, batchshift, 𝕡, 𝕀, q, N, M, catsampler) = buffer
    _onehalf = convert(T, 1//2)
    E, ∇E = energy(X)

    _gwg_proposal!(𝕡, X, ∇E, A₁)
    unitrange_rand!(I, 1:q*N, 𝕡, catsampler)
                                    # flipped spin and proposed spin state for each sample
    I .+= batchshift                # shift index with sample index

    α .= reshape(E, 1,1,M)
    @view(α[:]) .-= log.(@view(𝕡[I]))   # α = E(x) - ln(π(x′|x))

    fill!(A, zero(T))               # initialize
    fill!(@view(A[I]), one(T))      # onehot encode flipped spin and proposed spin state
    sum!(A₁, A)                     # onehot encode site of flipped spin
    @. A = X + (A - A₁ * X)         # proposed state
    @. A = A > _onehalf              # ensure the entries are in (0,1)
    @. 𝕀 = A₁ == X == 1             # onehot encode flipped spin and original spin state

    Ẽ, ∇Ẽ = energy(A)
    _gwg_proposal!(𝕡, A, ∇Ẽ, A₁)
    α .-= reshape(Ẽ, 1,1,M)         # α = E(x) - ln(π(x′|x)) - E(x′)

    @view(α[:]) .+= log.(𝕡[𝕀])      # α = E(x) - ln(π(x′|x)) - E(x′) + ln(π(x|x′))

    @. α = exp(-relu(-α))       # α = min{ 1, p(x′)π(x|x′) / p(x)π(x′|x) }
    @. α = α ≥ rand()           # accept proposal

    @. X = (1-α) * X + α * A    # update X
    @. X = X > _onehalf          # ensure the entries are in (0,1)
    return X
end

function _gwg_proposal!(𝕡::AbstractArray{T,3}, X, ∇E, A₁) where T
    _typemax = typemax(T)
    _onehalf = convert(T, 1//2)
    @. 𝕡 = X * ∇E
    sum!(A₁, 𝕡)
    @. 𝕡 = -_onehalf * (∇E - A₁)
    @. 𝕡-= (X > _onehalf) * _typemax # disallow transition to current state
    softmax!(𝕡, dims = 1:2)
end


text:
using PooksoftTextMagicAPI
using Test

function default_test()
    return true
end

@testset "pooksoft_text_magic_api_test_set" begin
    @test default_test() == true
end


text:
models = (
    AdaBoostClassifier,
    BaggingClassifier,
    GradientBoostingClassifier,
    RandomForestClassifier,
    ExtraTreesClassifier
)

fparams = (
    AdaBoostClassifier=(:estimators, :estimator_weights, :estimator_errors, :classes, :n_classes),
    BaggingClassifier=(:base_estimator, :estimators, :estimators_samples, :estimators_features, :classes, :n_classes, :oob_score, :oob_decision_function),
    GradientBoostingClassifier=(:n_estimators, :feature_importances, :train_score, :loss, :init, :estimators, :oob_improvement),
    RandomForestClassifier=(:estimators, :classes, :n_classes, :n_features, :n_outputs, :feature_importances, :oob_score, :oob_decision_function),
    ExtraTreesClassifier=(:estimators, :classes, :n_classes, :feature_importances, :n_features, :n_outputs, :oob_score, :oob_decision_function)
)

@testset "Fit/Predict" begin
    X, y = simple_binaryclf()
    for mod in models
        m = mod()
        acc, f = test_clf(m, X, y)
        @test acc >= 0.75
        @test keys(MB.fitted_params(m, f)) == getproperty(fparams, Symbol(mod))
    end
end

@testset "Docstrings" begin
    for mod in models
        m = mod()
        @test !isempty(MB.docstring(m))
    end
end

models = (
    AdaBoostRegressor,
    BaggingRegressor,
    GradientBoostingRegressor,
    RandomForestRegressor,
    ExtraTreesRegressor
)

fparams = (
    AdaBoostRegressor=(:estimators, :estimator_weights, :estimator_errors, :feature_importances),
    BaggingRegressor=(:estimators, :estimators_samples, :estimators_features, :oob_score, :oob_prediction),
    GradientBoostingRegressor=(:feature_importances, :train_score, :loss, :init, :estimators, :oob_improvement),
    RandomForestRegressor=(:estimators, :feature_importances, :n_features, :n_outputs, :oob_score, :oob_prediction),
    ExtraTreesRegressor=(:estimators, :feature_importances, :n_features, :n_outputs, :oob_score, :oob_prediction)
)

@testset "Fit/Predict" begin
    X, y, ls = simple_regression()
    for mod in models
        m = mod()
        r, f = test_regression(m, X, y, ls)
        @test r < 5
        @test keys(MB.fitted_params(m, f)) == getproperty(fparams, Symbol(mod))
    end
end


text:
########## Storage API

abstract type AbstractSalsaStorage end

# The functions for a Salsa Storage are declared below, and should all be implemented for
# concrete instances of AbstractSalsaStorage.
#
# NOTE that most users of Salsa can use the `DefaultStorage` included with Salsa.jl, but
# these are the methods to implement for a custom storage backend.

########## Lookups

# Methods here have a counterpart in the runtime API.

"""
    _memoized_lookup_internal(runtime::Runtime, key::InputKey{F,TT})
    _memoized_lookup_internal(runtime::Runtime, key::DerivedKey{F,TT})

Storage-specific implementation of [`Salsa.memoized_lookup`](@ref). Implementations should
look up `key` in the storage associated with that `runtime`, and return the (cached) value
if it exists.

For derived keys, if no valid cached value is found, compute the value from scratch by
running the derived function `F` indicated in the `DerivedKey{F}` type with the function
arguments extracted from `key.args`. The result is then cached again.

This is the main interface that Salsa uses to access values from the Salsa Storage, when
requested by users.

Providers of a Salsa Storage layer should implement this function for both types of Salsa
keys, by adding a method that dispatches on their specific type of `Runtime`, e.g.:

```julia
function Salsa._memoized_lookup_internal(
    rt::Salsa.RuntimeWithStorage{MyStorage}, key::DerivedKey{F,TT}
) where {F,TT}
    # ...
    # Lookup key in `storage(rt)`, and if not present, execute `F(rt, key.args)` to compute
    # it and cache it.
    # ...
end
```
"""
function _memoized_lookup_internal end

"""
    _unmemoized_input_lookup_internal(runtime::Runtime, input_key::InputKey)

Storage-specific implementation of [`Salsa.unmemoized_input_lookup`](@ref). Implementations
should look up the value of `input_key` without tracking it as a dependency.
"""
function _unmemoized_input_lookup_internal end

"""
    _previous_output_internal(runtime::Runtime)

Storage-specific implementation of [`Salsa.previous_output`](@ref). Implementations can
assume that they're being called from within a `@derived` function and should return the
currently cached value (which at this point must be invalid) for that function and its
arguments.
"""
function _previous_output_internal end

"""
    _unwrap_salsa_value(runtime::Runtime, v)

Some storage engines might need to store values in a wrapper type with additional metadata
(such as validity). This function is the last thing called before returning the value to the
user and gives the storage implementation a chance to unwrap the cached value.
"""
function _unwrap_salsa_value end

########## Writes

"""
    set_input!(runtime::Runtime, input_key::InputKey, value)

Set the value of an input in the storage associated with `runtime`. This function is called
by Salsa's higher-level API such as that generated by the `@declare_input` macro.
"""
function set_input! end

"""
    delete_input!(runtime::Runtime, input_key::InputKey)

Delte an input in the storage associated with `runtime`. This function is called by Salsa's
higher-level API such as that generated by the `@declare_input` macro.
"""
function delete_input! end


text:
@doc "Calculate the distance between point A (latA,lonA) and Point B(latB,lonB)"
function distanceInKm(latA::AbstractFloat,lonA::AbstractFloat,latB::AbstractFloat,lonB::AbstractFloat)
  return 2 * 6372.8 * asin(sqrt(sind((latB-latA)/2)^2 + cosd(latA) * cosd(latB) * sind((lonB - lonA)/2)^2))
end

@doc "Calculate the distance between point A(latA,lonA) and a vector of Points(latB,lonB)"
function distanceInKm(latA::AbstractFloat,lonA::AbstractFloat,latB::Vector{AbstractFloat},lonB::Vector{AbstractFloat})
    n = size(latB,1)
    resp::Vector{AbstractFloat}=Vector{AbstractFloat}(n)
    if (size(lonB,1)!=n)
     error("latB and lonB have different sizes")
    end
    @simd for i in 1:n  # @simd vectoraze for loop
      @inbounds resp[i] = distanceInKm(latA,lonA,latB[i],lonB[i]) # @inbounds to not check bounds
    end
end

@doc "Calculate the distance between a vector of points"
function distanceInKm(latA::Vector{AbstractFloat},lonA::Vector{AbstractFloat},latB::Vector{AbstractFloat},lonB::Vector{AbstractFloat})
    n = size(latA,1)
    isEqualSize::Bool = n == size(lonA,1) == size(latB,1)==size(lonB,1)

    if (!isEqualSize)
     error("the vectors have different sizes")
    end
    resp::Vector{AbstractFloat}=Vector{AbstractFloat}(n)
    @simd for i in 1:n  # @simd vectoraze for loop
      @inbounds resp[i] = distanceInKm(latA[i],lonA[i],latA[i],lonA[i])
    end
end

@doc "Calculate the Azimute between A and B"
function azimuthLineAtoB(latA::AbstractFloat,lonA::AbstractFloat,latB::AbstractFloat,lonB::AbstractFloat)
     numer = cosd(latB) * sind(lonB - lonA)
     denom = sind(latB) * cosd(latA) - cosd(latB) * sind(latA) * cosd(lonB - lonA)
     return rad2deg((atan2(numer, denom) + 2π) % 2π)
end


text:
using AMQPClient

channel!(conn) =
    channel(conn, AMQPClient.UNUSED_CHANNEL, true)

function close_channel!(chan)
    if isopen(chan)
        close(chan)
        AMQPClient.wait_for_state(chan, AMQPClient.CONN_STATE_CLOSED)
    end
end

function amqp_channel!(f::Function, conn)
    chan = channel!(conn)
    try
        f(chan)
    finally
        close_channel!(chan)
    end
end


text:
using Requires

write_patchwork_prelude(io::IO) =
    write(io, "<script>", Patchwork.js_runtime(), "</script>")

write_canvas_prelude(io::IO) = begin
    write(io, Canvas.custom_elements())
    write_patchwork_prelude(io)
end

@require Morsel begin
    # Allow route handlers to return Patchwork nodes

    import Meddle: MeddleRequest, Response

    function Morsel.prepare_response{ns, tag}(
            data::Elem{ns, tag}, req::MeddleRequest, res::Response)
        io = IOBuffer()
        Patchwork.write_patchwork_prelude(io)

        writemime(io, MIME"text/html"(), data)
        prepare_response(takebuf_string(io), req, res)
    end

    function Morsel.prepare_response(
            data::Tile, req::MeddleRequest, res::Response)

        io = IOBuffer()

        Canvas.write_canvas_prelude(io)

        writemime(io, MIME"text/html"(), render(data))
        prepare_response(takebuf_string(io), req, res)
    end
end

@require IJulia begin
    # Load custom element definitions

    using IJulia.CommManager
    import Base.Random: UUID, uuid4

    function setup_transport(sig::Input)
        id = makeid(sig)
        comm = Comm(:CanvasSignal, data=[:signalId => id])
        comm.on_msg = (msg) -> push!(sig, decodeJSON(sig, msg.content["data"]["value"]))
        return id
    end
end

@require Blink begin
    # This is still defunct though
    import BlinkDisplay, Graphics

    Blink.windowinit() do w
        Blink.head(w, custom_elements())
    end

    Graphics.media(Tile, Graphics.Media.Graphical)
end

@require Gadfly begin
    import Gadfly: Compose

    function convert(::Type{Tile}, p::Gadfly.Plot)
        backend = Compose.Patchable(
                     Compose.default_graphic_width,
                     Compose.default_graphic_height)
        convert(Tile, Compose.draw(backend, p))
    end
end

@require DataFrames begin
    include(Pkg.dir("Canvas", "src", "library", "table.jl"))
end


text:
function initswarmuniform(f, Opts::PSOOptions)

    # Initialize Swarm
    swarm = Swarm(Opts.NDims, Opts.SwarmSize)

    # Create Matrix of Initial Positions wth Uniform Distribution
    InitSpan = MMatrix{Opts.NDims,2}(zeros(Opts.NDims,2))
    @inbounds for i in 1:Opts.NDims
        if -Opts.InitialSwarmSpan[i]/2 > Opts.LowerBounds[i]
            InitSpan[i,1] = -Opts.InitialSwarmSpan[i]/2
        else
            InitSpan[i,1] = -Opts.LowerBounds[i]
        end
        if Opts.InitialSwarmSpan[i]/2 < Opts.UpperBounds[i]
            InitSpan[i,2] = Opts.InitialSwarmSpan[i]/2
        else
            InitSpan[i,2] = Opts.UpperBounds[i]
        end
    end

    if Opts.InitialSwarmMatrix[1] === NaN
        @views Pos = SizedMatrix{Opts.NDims, Opts.SwarmSize}(
               InitSpan[:,1] .+ (InitSpan[:,2] .- InitSpan[:,1]).*rand(Opts.NDims, Opts.SwarmSize)
        )
    elseif size(Opts.InitialSwarmMatrix)[2] >= Opts.SwarmSize
        Pos = SizedMatrix{Opts.NDims, Opts.SwarmSize}(
            view(Opts.InitialSwarmMatrix,:,1:Opts.SwarmSize)
        )
    else
        @views Pos = SizedMatrix{Opts.NDims, Opts.SwarmSize}(
            [
            Opts.InitialSwarmMatrix
            InitSpan[:,1] .+ (InitSpan[:,2] .- InitSpan[:,1]) .*
                rand(Opts.NDims, Opts.SwarmSize - size(Opts.InitialSwarmMatrix)[2])
            ]
        )
    end

    # Create Matrix of Initial Velocities with Uniform Distribution
    Vel = SizedMatrix{Opts.NDims, Opts.SwarmSize}(
        2 .* Opts.InitialSwarmSpan .* (rand(Opts.NDims, Opts.SwarmSize) .- 0.5)
    )

    # Set Each Particles Position and Velocity
    @inbounds for i in 1:Opts.SwarmSize
        swarm[i].x = SVector{Opts.NDims}(view(Pos,:,i))
        swarm[i].p = SVector{Opts.NDims}(view(Pos,:,i))
        swarm[i].v = SVector{Opts.NDims}(view(Vel,:,i))
    end

    # Compute Objective Function for all Particles in Swarm
    if Opts.UseParallel
        Threads.@threads for i in 1:Opts.SwarmSize
            swarm[i].fx = f(swarm[i].x)
            swarm[i].fp = deepcopy(swarm[i].fx)
        end
    else
        for i in 1:Opts.SwarmSize
            swarm[i].fx = f(swarm[i].x)
            swarm[i].fp = deepcopy(swarm[i].fx)
        end
    end

    # Check Objective Function Values if Desired
    if Opts.FunValCheck == "on"
        fvalcheck(swarm)
    end

    # Set Global Best
    setgbest!(swarm)

    # Initialize Neighborhood Size
    swarm.N = max(2, floor(Opts.SwarmSize*Opts.MinNeighborsFraction))

    # Initialize Inertia
    if Opts.InertiaRange[2] > 0
        if Opts.InertiaRange[2] > Opts.InertiaRange[1]
            swarm.W = Opts.InertiaRange[2]
        else
            swarm.W = Opts.InertiaRange[1]
        end
    else
        if Opts.InertiaRange[2] < Opts.InertiaRange[1]
            swarm.W = Opts.InertiaRange[2]
        else
            swarm.W = Opts.InertiaRange[1]
        end
    end

    # Initialize Self and Social Adjustment Weights
    swarm.y1 = Opts.SelfAdjustmentWeight
    swarm.y2 = Opts.SocialAdjustmentWeight

    # Print Status
    if Opts.Display == "iter"
        printstatus(swarm, 0.0, 0, 0)
    end

    return swarm
end

max_stars_repo_name: GrantHecht/JuSwarm

max_forks_repo_name: GrantHecht/JuSwarm


text:
# Utility functions

function isprobvec(p::Vector{Float64})
    s = 0.
    for i = 1:length(p)
        pi = p[i]
        s += pi
        if pi < 0
            return false
        end
    end
    return abs(s - 1.0) <= 1.0e-12
end

function _randu(Ku::Uint, U::Uint)   # ~ U[0:Ku-1]
	x = rand(Uint)
	while x > U
		x = rand(Uint)
	end
	rem(x, Ku)
end

function randi(K::Int)  # Fast method to draw a random integer from 1:K
	Ku = uint(K)
	U = div(typemax(Uint), Ku) * Ku
	int(_randu(Ku, U)) + 1
end

function randi(a::Int, b::Int)  # ~ U[a:b]
	Ku = uint(b - a + 1)
	U = div(typemax(Uint), Ku) * Ku
	int(_randu(Ku, U)) + a
end

immutable RandIntSampler
	a::Int
	Ku::Uint
	U::Uint

	function RandIntSampler(K::Int) # 1:K
		Ku = uint(K)
		U = div(typemax(Uint), Ku) * Ku
		new(1, Ku, U)
	end

	function RandIntSampler(a::Int, b::Int)  # a:b
		Ku = uint(b - a + 1)
		U = div(typemax(Uint), Ku) * Ku
		new(a, Ku, U)
	end
end

rand(s::RandIntSampler) = int(_randu(s.Ku, s.U)) + s.a


text:
export ContinuousWave
export GaussianPulse, DGaussianPulse
export BlackmanHarrisPulse, DBlackmanHarrisPulse

function fourier(t,signal)
	δt = step(t) ; N = length(t)
	signal_fft_mag = abs.(fft(signal))* δt
	freqs = fftfreq(N,1/δt)
	return freqs, signal_fft_mag
end

abstract type AbstractWaveform end

include("./Waveforms/ContinuousWaveForm.jl")
include("./Waveforms/GaussianWaveForm.jl")
include("./Waveforms/BlackmanHarrisWaveForm.jl")


text:
@testset "convolutional constructors" begin
    s = (8, 6, 2, 6)
    x = randn(Float32, s) |> gpu
    xsize = s[1:3]
    ldim = 2
    kernelsizes = (3, 5)
    nchannels = (4, 8)
    scalings = (2, 1)
    densedims = [128]

    # encoder
    encoder = GMExtensions.conv_encoder(xsize, ldim, kernelsizes, nchannels, scalings) |> gpu
    y = encoder(x)
    @test length(encoder.layers) == 2*length(kernelsizes) + 1 + 1 # 2*(conv + maxpool) + reshape + dense
    @test size(y) == (ldim, s[end])
    @test eltype(y) == Float32

    encoder = GMExtensions.conv_encoder(xsize, ldim, kernelsizes, nchannels, scalings;
        densedims = densedims) |> gpu
    @test length(encoder.layers) == 2*length(kernelsizes) + 1 + 2 # 2*(conv + maxpool) + reshape + 2*dense

    # decoder
    decoder = GMExtensions.conv_decoder(xsize, ldim, reverse(kernelsizes),
        reverse(nchannels), reverse(scalings)) |> gpu
    z = decoder(y)
    @test length(decoder.layers) == 1 + 1 + 2*length(kernelsizes) # dense + reshape + 2*(conv + maxpool)
    @test size(z) == size(x)
    @test eltype(z) == Float32

    decoder = GMExtensions.conv_decoder(xsize, ldim, reverse(kernelsizes),
        reverse(nchannels), reverse(scalings); densedims = densedims) |> gpu
    @test length(decoder.layers) == 2 + 1 + 2*length(kernelsizes) # dense + reshape + 2*(conv + maxpool)

    # also test trainability
    enc_params = get_params(encoder)
    dec_params = get_params(decoder)
    loss(x) = Flux.mse(decoder(encoder(x)), x)
    opt = ADAM()
    data = (x,)
    GenerativeModels.update_params!(encoder, data, loss, opt)
    @test all(param_change(enc_params, encoder))
    GenerativeModels.update_params!(decoder, data, loss, opt)
    @test all(param_change(dec_params, decoder))

    # test vectorization of output
    decoder = GMExtensions.conv_decoder(xsize, ldim, reverse(kernelsizes),
        reverse(nchannels), reverse(scalings); densedims = densedims, vec_output = true) |> gpu
    z = decoder(y)
    @test size(z) == (reduce(*, xsize), 6)
    @test size(reshape(z, s...)) == s

    decoder = GMExtensions.conv_decoder(xsize, ldim, reverse(kernelsizes),
        reverse(nchannels), reverse(scalings); densedims = densedims, vec_output = true,
        vec_output_dim = reduce(*, xsize) + 1) |> gpu
    z = decoder(y)
    @test size(z) == (reduce(*, xsize) + 1, 6)

end


text:
using Comonicon, PlutoCLI
Comonicon.install(PlutoCLI)