text:
"""An `AbstractSpec` object contains all the immutable specifications in a
wavelet filter bank"""
abstract AbstractSpec{
    T<:Number,
    D<:AbstractDomain,
    G<:AbstractPointGroup,
    W<:RedundantWaveletClass}

"""Enforces properties of the wavelets to satisfy null mean, limited spatial
support, and Littlewood-Paley inequality.

* If the wavelet is dyadic, the maximum required quality factor
`max_qualityfactor` must be equal to `1.0`.

* The truncation threshold `ɛ` must be `in [0.0, 1.0[`.

* The signal length must be a power of two above 4, i.e. `log2_size > 1`.

* The maximum required quality factor `max_qualityfactor` must be between `1.0`
and the number of filters per octaves.

* The maximum scale must `max_scale` must be at least about 5 times above the
maximum quality factor.

* The mother frequency `motherfrequency` must be in `]0.0, 0.5]`.

* The lowest center frequency must be greater or equal than the number of
per octaves, i.e. `(log2_size-n_octaves) >= 1 + log2(n_filters_per_octave)`."""
function checkspec_super(spec::AbstractSpec)
    if isdyadic(spec.class) && (spec.max_qualityfactor != 1.0)
        error("Maximum quality factor must be equal to `1.0` ",
        "for dyadic wavelets.\n",
        "`class` = ", spec.class, "\n",
        "`max_qualityfactor` = ", spec.max_qualityfactor, ".")
    end
    if (spec.ɛ >= 1.0) || (spec.ɛ < 0.0) || (spec.ɛ === -0.0)
        error("`ɛ` must be in `[0.0, 1.0[`. A typical value is `1e-4`.")
    end
    if spec.log2_size .< 2
        error("Too short signal length.\n",
        "`log2_size =` ", spec.log2_size, "\n",
        "`log2_size` must be `≧2`")
    end
    if spec.max_qualityfactor < 1.0
        error("Too small maximum quality factor.\n",
        "`max_qualityfactor =` ", spec.max_qualityfactor, "\n",
        "`max_qualityfactor` must be `≧1.0`.")
    end
    if spec.motherfrequency <= 0.0 || spec.motherfrequency > 0.5
        error("`motherfrequency` must be in `]0.0, 0.5]`.")
    end
    if spec.n_filters_per_octave < 1
        error("Too few filters per octave.\n",
        "`n_filters_per_octave = `", spec.n_filters_per_octave, "\n",
        "`n_filters_per_octave` must be `≧1`.")
    end
    if spec.n_octaves < 1
        error("Too few octaves.\n",
        "`n_octaves = `", spec.n_octaves, "\n",
        "`n_octaves` must be `≧1`")
    end
    if spec.max_qualityfactor > spec.n_filters_per_octave
        error("Too few filters per octave for the given quality factor.\n",
        "`max_qualityfactor = `", spec.max_qualityfactor, "\n",
        "`n_filters_per_octave = `", spec.n_filters_per_octave, "\n",
        """The inequality `n_filters_per_octave ≧ max_qualityfactor` must be
        satisfied.""")
    end
    if spec.n_octaves >= spec.log2_size
        error("Too many octaves.\n",
        "`log2_size = `", spec.log2_size, "\n",
        "`n_octaves = `", spec.n_octaves, "\n",
        """The inequality `minimum(log2_size) > n_octaves` must be satisfied.
        Either increase `log2_size` or decrease `n_octaves`.""")
    end
    if spec.log2_size-spec.n_octaves < 1+log2(spec.n_filters_per_octave)
        error("Too many filters per octave for the given length.\n",
        "`log2_size = `", spec.log2_size, "\n",
        "`log2(n_filters_per_octave) = `", log2(spec.n_filters_per_octave), "\n",
        "`n_octaves = `", spec.n_octaves, "\n",
        """The inequality
        `minimum(log2_size)-n_octaves ≧ 1 + log2(n_filters_per_octave)`
        must be satisfied. Either increase `log2_size`, decrease `n_octaves`,
        or decrease `n_filters_per_octave`.""")
    end
    empirical_max_ψscale = mapreduce(get_scale, max, spec.ψmetas)
    empirical_max_scale = max(empirical_max_ψscale, spec.ϕmeta.scale)
    if empirical_max_scale > (spec.max_scale + 1e-3)
        error("Required time-frequency localization is too tight.\n",
        "`max_qualityfactor = `", spec.max_qualityfactor, "\n",
        "`max_scale = `", spec.max_scale, "\n",
        "`motherfrequency = `", spec.motherfrequency, "\n",
        "The wavelet ", typeof(spec), "cannot have both a bandwidth `< ``",
        spec.motherfrequency / spec.max_qualityfactor,
        "and a scale `<` ", spec.max_scale, ".\n",
        "Either decrease `max_qualityfactor` or decrease `max_scale`.")
    end
    if empirical_max_scale > (2.0^spec.log2_size + 1e-3)
        min_resolution = 2.0^(-spec.n_octaves/spec.n_filters_per_octave)
        min_centerfrequency = spec.motherfrequency * min_resolution
        max_bandwidth = min_centerfrequency / spec.max_qualityfactor
        size = 1 .<< spec.log2_size
        error("Spatial localization is coarser than signal length.\n",
        "`log2_size = ``", spec.log2_size,
        "`max_qualityfactor = ``", spec.max_qualityfactor,
        "`max_scale = ``", spec.max_scale,
        "`motherfrequency = ``", spec.motherfrequency, "\n",
        "`n_octaves = ``", spec.n_octaves, "\n",
        "The wavelet ", typeof(spec), "cannot have both a bandwidth ``< ",
        "motherfrequency*2^(-n_octaves)/qualityfactor = ``", max_bandwidth,
        "and a scale `< 2^(log2_size) = `", size, ".\n",
        """Either increase log2_size, decrease max_qualityfactor,
        set `max_scale <= log2_size`, or decrease `n_octaves`.""")
    end
    return true
end

"""Given an signal input type `T`, returns a conservative default value for
the truncation threshold `ɛ of the wavelets.
The chosen values are powers of ten within the same order of magnitude as
floating-point machine precision `eps(T)`."""
default_ɛ(T::Type{Float16}) = 1e-3
default_ɛ(T::Type{Float32}) = 1e-7
default_ɛ(T::Type{Float64}) = 1e-15
default_ɛ{RealT}(T::Type{Complex{RealT}}) = default_ɛ(RealT)

"""Given a maximum quality factor and a number of filter per octaves (both of
which may be `Void`), returns the maximum quality factor in a wavelet filter
bank."""
default_max_qualityfactor(max_q::Real, nfo::Any) = Float64(max_q)
default_max_qualityfactor(max_q::Void, nfo::Integer) = Float64(nfo)
default_max_qualityfactor(max_q::Void, nfo::Void) = 1.0

"""The dimensionless mother center frequency `ξ` (corresponding to a log-period
`γ=0`) is computed as the midpoint between the center frequency of the second
wavelet `ξ*2^(-1/n_filters_per_octave)` (corresponding to `γ=1`) and the
negative mother center frequency `(1-ξ)`. Hence the equation
`2ξ = ξ*2^(-1/n_filters_per_octave) + (1-ξ)`, of which we
derive `ξ = 1 / (3 - 2^(-1/n_filters_per_octave))`. This formula is valid
only when the wavelet is a symmetric bump in the Fourier domain."""
default_motherfrequency(class::RedundantWaveletClass, n_filters_per_octave) =
    inv(3.0 - exp2(-inv(n_filters_per_octave)))

"""Given a maximum quality factor and a number of filter per octaves (both of
which may be `Void`), returns the default number of filters per octave in a
wavelet filter bank."""
default_n_filters_per_octave(nfo::Integer, max_q::Any) = Int(nfo)
default_n_filters_per_octave(nfo::Void, max_q::Real) = ceil(Int, max_q)
default_n_filters_per_octave(nfo::Void, max_q::Void) = 1

"""Returns the maximal number octaves in a filter bank such that all scales are
below `2^(log2_size)`."""
default_n_octaves(n_octaves::Int, args...) = n_octaves
function default_n_octaves(n_octaves::Void, class::RedundantWaveletClass,
        log2_size::Int, max_qualityfactor::Float64, max_scale::Real,
        motherfrequency::Float64, n_filters_per_octave::Int, args...)
    siglength = 1 << minimum(log2_size)
    if max_scale > siglength
        min_centerfrequency = uncertainty(class) / siglength * max_qualityfactor
    else
        min_centerfrequency = uncertainty(class) / max_scale * 1.0
    end
    n_octaves_a = floor(Int, log2(motherfrequency / min_centerfrequency))
    n_octaves_b = minimum(log2_size) - 1 - ceil(Int, log2(n_filters_per_octave))
    return min(n_octaves_a, n_octaves_b)
end

"""Fallback of the uncertainty constant from the spec to its class. The RHS
method must be specifically implemented by `AbstractSpec` concrete subtypes."""
uncertainty(spec::AbstractSpec) = uncertainty(spec.class)


text:
# m5.3s.jl

using Pkg, DrWatson

using MonteCarloMeasurements
using StanSample
using StatisticalRethinking

# ### snippet 5.1

df = CSV.read(sr_datadir("WaffleDivorce.csv"), DataFrame);
scale!(df, [:Marriage, :MedianAgeMarriage, :Divorce])

# Define the Stan language model

stan5_3 = "
data {
    int N;
    vector[N] divorce_s;
    vector[N] marriage_s;
    vector[N] medianagemarriage_s;
}
parameters {
    real a;
    real bA;
    real bM;
    real<lower=0> sigma;
}
model {
    vector[N] mu = a + + bA * medianagemarriage_s + bM * marriage_s;
    a ~ normal( 0 , 0.2 );
    bA ~ normal( 0 , 0.5 );
    bM ~ normal( 0 , 0.5 );
    sigma ~ exponential( 1 );
    divorce_s ~ normal( mu , sigma );
}
";

# Define the SampleModel
m5_3s = SampleModel("m5.3", stan5_3);

# Input data

m5_3_data = Dict(
    "N" => size(df, 1),
    "divorce_s" => df[:, :Divorce_s],
    "marriage_s" => df[:, :Marriage_s],
    "medianagemarriage_s" => df[:, :MedianAgeMarriage_s]
);

# Sample using cmdstan

rc5_3s = stan_sample(m5_3s, data=m5_3_data);

if success(rc5_3s)

    # Rethinking results

    rethinking_results = "
           mean   sd  5.5% 94.5%
    a      0.00 0.10 -0.16  0.16
    bM    -0.07 0.15 -0.31  0.18
    bA    -0.61 0.15 -0.85 -0.37
    sigma  0.79 0.08  0.66  0.91
    ";

    part5_3s = read_samples(m5_3s, :particles)
    part5_3s |> display

    read_samples(m5_3s) |> display

end

# End m5.3s.jl

ext: jl


text:
using OptimalEstimationProject
using LinearAlgebra
using AstroTime
using MATLAB

function main()
    # Spacecraft simulator
    t0      = TAIEpoch("2021-08-08T12:00:00")
    scsim   = SpacecraftSim(t0.second)

    # Initial state error variance
    σr      = 0.1   # [km]
    σv      = 0.01  # [km / s]
    σm      = 1e-6  # [kg]
    sqrtP0  = Diagonal([σr, σr, σr, σv, σv, σv, σm])
    P0      = sqrtP0.^2

    # Initial state estimate
    (y0true, us) = OptimalEstimationProject.GetStateAndControl(scsim, t0.second)
    xhat0   = y0true[1:7] + sqrtP0*randn(7)

    # Measurement time stamps
    Δt      = 30        # [sec]
    gpsΔt   = 15*60     # [sec]
    ts      = range(t0.second + Δt; step = Δt, stop = t0.second + scsim.ts[2]*86400)

    # Measurement statistics
    σρ      = 10.0e-3 # [km]     Pseudorange noise standard deviation
    σr      = 5.0e-3 # [km]     GPS broadcast ephemeris standard deviation
    σa      = 1.0e-6 # [km/s^2] Accelerometer noise standard deviation

    # Process noise covariance
    R       = Diagonal((σρ^2 + 3*σr^2)*ones(32))
    Q       = Diagonal([0.0, 0.0, 0.0, 1e-12, 1e-12, 1e-12, 1e-3])

    # GPS Simulation Span
    startWeek   = 2170
    startDay    = 0
    endWeek     = 2174
    endDay      = 6

    # Create GPS simulation object
    gpssim = GPSSim(startWeek, startDay, endWeek, endDay; σρ = σρ, σr = σr)

    # Create IMU simulation object
    imusim = IMUSim(σa)

    # Create EKF
    ekf = EKF(xhat0, P0, Q, (σρ^2 + 3*σr^2), σa^2, ts, gpsΔt, gpssim, imusim, scsim; steps2save = 2, lunaPerts = true);

    # Run filter
    runFilter!(ekf)

    # Get true trajectory for plotting
    xtrue = zeros(length(ekf.txp), 3)
    for i in 1:length(ekf.txp)
        (y0l, us) = OptimalEstimationProject.GetStateAndControl(scsim, ekf.txp[i])
        xtrue[i, :] .= y0l[1:3]
    end

    # Plotting
    plot(ekf, xtrue, length(ekf.txp))
end

main()

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
@dict_readable struct WorkspaceFoldersChangeEvent
    added::Vector{WorkspaceFolder}
    removed::Vector{WorkspaceFolder}
end

@dict_readable struct DidChangeWorkspaceFoldersParams
    event::WorkspaceFoldersChangeEvent
end

@dict_readable struct DidChangeConfigurationParams
    settings::Any
end

struct ConfigurationItem <: Outbound
    scopeUri::Union{DocumentUri,Missing}
    section::Union{String,Missing}
end

struct ConfigurationParams <: Outbound
    items::Vector{ConfigurationItem}
end

##############################################################################
# File watching
const FileChangeType = Int
const FileChangeTypes = Dict("Created" => 1, "Changed" => 2, "Deleted" => 3)
@dict_readable struct FileEvent
    uri::String
    type::FileChangeType
end

@dict_readable struct DidChangeWatchedFilesParams
    changes::Vector{FileEvent}
end

const WatchKind = Int
const WatchKinds = Dict("Create" => 1,
                       "Change" => 2,
                       "Delete" => 4)
struct FileSystemWatcher <: Outbound
    globPattern::String
    kind::Union{WatchKind,Missing}
end

struct DidChangeWatchedFilesRegistrationOptions <: Outbound
    watchers::Vector{FileSystemWatcher}
end

##############################################################################

struct CreateFileOptions <: Outbound
    overwrite::Union{Bool,Missing}
    ignoreIfExists::Union{Bool,Missing}
end

struct CreateFile <: Outbound
    kind::String
    uri::DocumentUri
    options::Union{CreateFileOptions,Missing}
    CreateFile(uri, options=missing) = new("create", uri, options)
end

struct RenameFileOptions <: Outbound
    overwrite::Union{Bool,Missing}
    ignoreIfExists::Union{Bool,Missing}
end

struct RenameFile <: Outbound
    kind::String
    oldUri::DocumentUri
    newUri::DocumentUri
    options::Union{RenameFileOptions,Missing}
    RenameFile(uri, options=missing) = new("rename", uri, options)
end

struct DeleteFileOptions <: Outbound
    overwrite::Union{Bool,Missing}
    ignoreIfNotExists::Union{Bool,Missing}
end

struct DeleteFile <: Outbound
    kind::String
    uri::DocumentUri
    options::Union{DeleteFileOptions,Missing}
    DeleteFile(uri, options=missing) = new("delete", uri, options)
end

##############################################################################
# Registration

struct Registration <: Outbound
    id::String
    method::String
    registerOptions::Union{Any,Missing}
end

struct RegistrationParams <: Outbound
    registrations::Vector{Registration}
end

struct TextDocumentRegistrationOptions <: Outbound
    documentSelector::Union{DocumentSelector,Nothing}
end

struct TextDocumentChangeRegistrationOptions <: Outbound
    documentSelector::DocumentSelector
    syncKind::Int
end

struct TextDocumentSaveRegistrationOptions <: Outbound
    documentSelector::DocumentSelector
    includeText::Union{Bool,Missing}
end

struct Unregistration <: Outbound
    id::String
    method::String
end

struct UnregistrationParams <: Outbound
    unregistrations::Vector{Unregistration}
end

struct StaticRegistrationOptions <: Outbound
    id::Union{String,Missing}
end

##############################################################################


text:
initial = open("conway_cubes.txt") do f
    BitArray(permutedims(hcat(map(r -> map(x ->  x == '#' ? true : false, collect(r)), readlines(f))...), (2,1)))
end

function expand(m)
    if sum(m[1, :, :, :]) > 0 # roof
        # println("stretching roof")
        m = vcat(falses(1, size(m, 2), size(m, 3), size(m, 4)), m)
    end
    if sum(m[end, :, :, :]) > 0 # floor
        # println("stretching floor")
        m = vcat(m, falses(1, size(m, 2), size(m, 3), size(m, 4)))
    end
    if sum(m[:, 1, :, :]) > 0 # left
        # println("stretching left")
        m = hcat(falses(size(m, 1), 1, size(m, 3), size(m, 4)), m)
    end
    if sum(m[:, end, :, :]) > 0 # right
        # println("stretching right")
        m = hcat(m, falses(size(m, 1), 1, size(m, 3), size(m, 4)))
    end
    if sum(m[:, :, 1, :]) > 0 # back
        # println("stretching back")
        m = cat(falses(size(m, 1), size(m, 2), 1, size(m, 4)), m, dims = 3)
    end
    if sum(m[:, :, end, :]) > 0 # front
        # println("stretching front")
        m = cat(m, falses(size(m, 1), size(m, 2), 1, size(m, 4)), dims = 3)
    end
    if sum(m[:, :, :, 1]) > 0 # hyperback
        # println("stretching hyperback")
        m = cat(falses(size(m, 1), size(m, 2), size(m, 3), 1), m, dims = 4)
    end
    if sum(m[:, :, :, end]) > 0 # front
        # println("stretching hyperfront")
        m = cat(m, falses(size(m, 1), size(m, 2), size(m, 3), 1), dims = 4)
    end
    return m
end

function index2coord(i, s)
    w = (i - 1) ÷ (s[1] * s[2] * s[3])
    z = (i - 1 - w * s[1] * s[2] * s[3]) ÷ (s[1] * s[2])
    x = (i - 1 - z * s[1] * s[2] - w * s[1] * s[2] * s[3]) ÷ s[1]
    y = i - 1 - x * s[1] - z * s[1] * s[2] - w * s[1] * s[2] * s[3]
    return (x+1, y+1, z+1, w+1)
end

function create_range(x, max)
    xmin = x <= 1 ? 1 : x - 1
    xmax = x >= max ? max : x + 1
    return xmin:xmax
end

function print_map(x)
    for w in 1:size(x, 4)
        for i in 1:size(x, 3)
            println("w = $(w), z = $(i)")
            for j in 1:size(x, 1)
                println(join(map(x -> x ? "#" : ".", x[j, :, i, w])))
            end
        end
    end
end

n = reshape(initial, (size(initial)..., 1, 1))
# print_map(n)

for j in 1:6
    global n
    n = expand(n)
    nc = deepcopy(n)

    for i in eachindex(n)
        x, y, z, w = index2coord(i, size(n))
        xr = create_range(x, size(n, 2))
        yr = create_range(y, size(n, 1))
        zr = create_range(z, size(n, 3))
        wr = create_range(w, size(n, 4))
        s = sum(view(n, yr, xr, zr, wr)) - Int(n[i])
        # println("$(i): $(index2coord(i, size(n))): $(n[i]): $(s)")
        if n[i] && s ∉ 2:3
            nc[i] = false
        end
        if !n[i] && s == 3
            nc[i] = true
        end
    end
    n = nc
    # println("i = $(j)")
    # print_map(n)
    # println()
end

println(sum(n))


text:
int main() {
  boolean[] bools = new boolean[10];
  double[] doubles = new double[10];
  int j=0;
  double x = 0.0;
  while (j<bools.length) {
    doubles[j] = x*2.5;
    bools[j] = doubles[j] > 15.0;
    j++;
    x = x + 1.0;
  }

  j = 0;
  for(boolean b : bools) {
    if(b) {
      printString("yes");
    } else {
      printDouble(doubles[j]);
    }
    j++;
  }
  return 0;
}

max_stars_repo_name: fredefox/javalette

max_issues_repo_name: fredefox/javalette

max_forks_repo_name: fredefox/javalette


text:
function capturing_analysis(expr, out, is_literal)

    @match expr begin
        ::QuoteNode =>
             capturing_analysis(expr.value, out, true)

        if is_literal end && Expr(:$, args...) =>

            foreach(x -> capturing_analysis(x, out, false), args)

        if is_literal end && Expr(_, args...) =>

            foreach(x -> capturing_analysis(x, out, true), args)

        if is_literal end && _ => nothing

        # not literal
        ::Symbol => (push!(out, expr); nothing)

        Expr(:quote, args...) =>
            foreach(x -> capturing_analysis(x, out, true), args)

        :(Do($(args...))) =>
            foreach(args) do arg
                @match arg begin
                    Expr(:kw, key :: Symbol, value) =>
                        begin
                            push!(out, key)
                        end
                    _ => nothing
                end
            end
        :($a || $b) =>
        let out1 = Set(Symbol[]),
            out2 = Set(Symbol[])

            capturing_analysis(a, out1, false)
            capturing_analysis(b, out2, false)

            union!(out, intersect(out1, out2))
            nothing
        end
        # type pattern
        Expr(:(::), a, _) => capturing_analysis(a, out, false)
        # dict pattern
        :(Dict($(args...))) =>
            foreach(args) do arg
                @match arg begin
                    :($_ => $v) => capturing_analysis(v, out, false)
                    _ => nothing
                end
            end
        # app pattern
        :($_($(args...))) => foreach(x -> capturing_analysis(x, out, false), args)
        # other expr
        Expr(_, args...) => foreach(x -> capturing_analysis(x, out, false), args)
        # ref pattern
        Expr(:&, _)              ||
        # predicate
        Expr(:function, _...)    ||
        Expr(:if, _...)          ||
        x => nothing
    end
end


text:
@testset "Univariate functions" begin
    # bspline0 from Daisy benchmarks
    f = x -> (1 - x)^3 / 6.0
    dom = Interval(-4.5, -0.3)
    xref = Interval(0.36616666666666675, 27.729166666666668)

    for solver in available_solvers
        x = enclose(f, dom, solver)
        rleft, rright = relative_precision(x, xref)
        if solver isa TaylorModelsEnclosure
            # for this example, TaylorModels cannot tightly approximate the lhs
            @test rleft ≤ 30 && rright ≤ 1e-10
        elseif solver isa SumOfSquaresEnclosure
            # solver returns negative -1.06259e-06
            @test rleft ≥ -1e-5 && rright ≤ 1e-5
        else
            @test rleft ≤ 1e-5 && rright ≤ 1e-5
        end
    end
end

@testset "Univariate example the Quickstart Guide" begin
    f(x) = -x^3/6 + 5x
    dom = Interval(1, 4)

    x = enclose(f, dom, NaturalEnclosure())
    xref = Interval(-5.66667, 19.8334)
    rleft, rright = relative_precision(x, xref)
    @test rleft ≤ 1e-5 && rright ≤ 1e-5

    x = enclose(f, dom, TaylorModelsEnclosure(order=4))
    xref = Interval(-4.2783, 12.7084)
    rleft, rright = relative_precision(x, xref)
    @test rleft ≤ 1e-5 && rright ≤ 1e-5

    x = enclose(f, dom, MooreSkelboeEnclosure())
    xref = Interval(4.83299, 10.5448)
    rleft, rright = relative_precision(x, xref)
    @test rleft ≤ 1e-5 && rright ≤ 1e-5

    x = enclose(f, dom, SumOfSquaresEnclosure(backend=SDPA.Optimizer))
    xref = Interval(4.8333, 10.541)
    rleft, rright = relative_precision(x, xref)
    @test rleft ≤ 1e-5 && rright ≤ 1e-5
end

@testset "Test univariate polynomial input" begin
    @polyvar x
    p = -x^3/6 + 5x
    dom = Interval(1, 4)

    x = enclose(p, dom)
    xref = Interval(-5.66667, 19.8334)
    rleft, rright = relative_precision(x, xref)
    @test rleft ≤ 1e-5 && rright ≤ 1e-5
end


text:
using SafeTestsets

@safetestset "Analog Circuits" begin include("analog.jl") end
#@safetestset "Digital Circuits" begin include("digital.jl") end
@safetestset "RC Circuit Demo" begin include("demo.jl") end
@safetestset "Thermal Circuits" begin include("thermal.jl") end


text:
# Glyph has to be a Bokehjs type as it's defined directly in Bokeh's JSON, but these methods are
# more porcelain than plumbing, so are defined here.
function Bokehjs.Glyph(glyphtype::Symbol,
                       linecolor::NullString,
                       linewidth::NullInt,
                       linealpha::NullFloat,
                       fillcolor::NullString,
                       fillalpha::NullFloat,
                       size::NullInt,
                       dash::Union{Void, Vector{Int64}},
                       fields::Union{Void, Dict{Symbol, Symbol}})
    props = Dict(
        :linecolor => linecolor == nothing ? omit : Dict{Symbol, BkAny}(:value => linecolor),
        :linewidth => linewidth == nothing ? omit : Dict{Symbol, BkAny}(:units => :data, :value =>linewidth),
        :linealpha => linealpha == nothing ? omit : Dict{Symbol, BkAny}(:units => :data, :value =>linealpha),
        :fillalpha => fillalpha == nothing ? omit : Dict{Symbol, BkAny}(:units => :data, :value =>fillalpha),
        :size => size == nothing ? omit : Dict{Symbol, BkAny}(:units => :screen, :value =>size),
        :fillcolor => fillcolor == nothing ? omit : Dict{Symbol, BkAny}(:value =>fillcolor),
    )
    if fields != nothing
        for (field, val) in fields
            if !haskey(props, field)
                error("unknown field $(field) passed to Glyph")
            end
            props[field] = Dict{Symbol,BkAny}(:field => val, :units => :data)
        end
    end
    Glyph(Bokehjs.uuid4(),
          glyphtype,
          props[:linecolor],
          props[:linewidth],
          props[:linealpha],
          props[:fillcolor],
          props[:fillalpha],
          props[:size],
          dash == nothing ? omit : dash,
          Dict(:units =>:data, :field => :x),
          Dict(:units =>:data, :field => :y))
end

function Bokehjs.Glyph(;glyphtype=nothing,
                        linecolor=nothing,
                        linewidth=nothing,
                        linealpha=nothing,
                        fillcolor=nothing,
                        fillalpha=nothing,
                        size=nothing,
                        dash=nothing,
                        fields=nothing)
    glyphtype = glyphtype == nothing ? :Line : glyphtype
    Glyph(glyphtype, linecolor, linewidth, linealpha, fillcolor, fillalpha,
          size, dash, fields)
end

function Bokehjs.Glyph(glyphtype::Symbol; kwargs...)
    Glyph(glyphtype=glyphtype; kwargs...)
end

function Base.show(io::IO, g::Bokehjs.Glyph)
    names = Glyph.names
    features = AbstractString[]
    for field in fieldnames(Glyph)
        showname = field == :_type_name ? :type : field
        g.(field) != nothing && push!(features, "$showname: $(g.(field))")
    end
    print(io, "Glyph(", join(features, ", "), ")")
end

type BokehDataSet
    data::Dict{Symbol, Vector}
    glyph::Glyph
    legend::NullString

    function BokehDataSet(data::Dict{Symbol, Vector}, glyph::Glyph,
                          legend::NullString=nothing)
        new(data, glyph, legend)
    end
end

function BokehDataSet(xdata::RealVect, ydata::RealVect, args...)
    data = Dict{Symbol, Vector}(:x => xdata, :y => ydata)
    BokehDataSet(data, args...)
end

type Plot
    datacolumns::Vector{BokehDataSet}
    tools::Vector{Symbol}
    filename::AbstractString
    title::AbstractString
    width::Int
    height::Int
    x_axis_type::NullSymbol
    y_axis_type::NullSymbol
    legendsgo::NullSymbol
end

# both the string and the reversed string will be tried, eg. "ox" is equivilent
# to "xo"
const STRINGTOKENS = let
    temp = Dict{ASCIIString,Dict{Symbol,Any}}(
            "--" => Dict(:dash      => [4, 4]),
            "-." => Dict(:dash      => [1, 4, 2]),
            "ox" => Dict(:glyphtype => :CircleX),
            "o+" => Dict(:glyphtype => :CircleCross),
            "sx" => Dict(:glyphtype => :SquareX),
            "s+" => Dict(:glyphtype => :SquareCross),
        )

    for (k, v) in Dict(temp)
        temp[reverse(k)] = v
    end
    temp
end

# heavily borrowed from Winston, thanks Winston!
const CHARTOKENS = Dict(
    '-'=> Dict(:dash=>nothing),
    ':'=> Dict(:dash=>[1, 4]),
    ';'=> Dict(:dash=>[1, 4, 2]),
    '+'=> Dict(:glyphtype=>:Cross),
    'o'=> Dict(:glyphtype=>:Circle),
    '*'=> Dict(:glyphtype=>:Asterisk),
    '.'=> Dict(:glyphtype=>:Circle, :size=>2),
    'x'=> Dict(:glyphtype=>:X),
    's'=> Dict(:glyphtype=>:Square),
    'd'=> Dict(:glyphtype=>:Diamond),
    '^'=> Dict(:glyphtype=>:Triangle),
    'v'=> Dict(:glyphtype=>:InvertedTriangle),
    'y'=> Dict(:linecolor=>"yellow"),
    'm'=> Dict(:linecolor=>"magenta"),
    'c'=> Dict(:linecolor=>"cyan"),
    'r'=> Dict(:linecolor=>"red"),
    'g'=> Dict(:linecolor=>"green"),
    'b'=> Dict(:linecolor=>"blue"),
    'w'=> Dict(:linecolor=>"white"),
    'k'=> Dict(:linecolor=>"black"),
)

Base.convert(::Type{Array{Glyph, 1}}, glyph::Glyph) = [glyph]

function Base.convert(::Type{Array{Glyph, 1}}, styles::AbstractString)
    map(style -> convert(Glyph, style), split(styles, '|'))
end

function Base.convert(::Type{Glyph}, style::AbstractString)
    # start wil default style
    styd = Dict(:glyphtype=>:Line, :linecolor=>"blue", :linewidth=>1,
                :linealpha=>1.0)

    # loop over all string tokens
    for key in keys(STRINGTOKENS)
        splitstyle = split(style, key)

        # if we found a match, apply it
        if length(splitstyle) > 1
            for (k, v) in STRINGTOKENS[key]
                styd[k] = v
            end

            # then remove it from `style` to start next loop
            style = join(splitstyle)
        end
    end

    # if we still have any char tokens left, process them
    for char in style
        if haskey(CHARTOKENS, char)
            for (k, v) in CHARTOKENS[char]
                styd[k] = v
            end
        else
            warn("unrecognized char '$char'")
        end
    end

    # set sane defaults for glyphs that are supposed to be filled or empty
    filledglyphs = [:Circle, :Square, :Diamond, :Triangle, :InvertedTriangle]
    if in(styd[:glyphtype], filledglyphs)
        styd[:fillcolor] = styd[:linecolor]
        styd[:fillalpha] = DEFAULT_FILL_ALPHA
        styd[:linealpha] = DEFAULT_FILL_ALPHA
        if !haskey(styd, :size)
            styd[:size] = DEFAULT_SIZE
        end
    end

    emptyglyphs = [:CircleX, :CircleCross, :SquareX, :SquareCross]
    if in(styd[:glyphtype], emptyglyphs)
        styd[:fillcolor] = "transparent"
        styd[:size] = DEFAULT_SIZE
    end

     # set sane defaults for glyphs that are supposed to be filled or empty
    Glyph(;styd...)
end


text:
abstract type LTISystem{TE<:TimeEvolution} <: AbstractSystem end
+(sys1::LTISystem, sys2::LTISystem) = +(promote(sys1, sys2)...)
-(sys1::LTISystem, sys2::LTISystem) = -(promote(sys1, sys2)...)
*(sys1::LTISystem, sys2::LTISystem) = *(promote(sys1, sys2)...)
/(sys1::LTISystem, sys2::LTISystem) = /(promote(sys1, sys2)...)

# Fallback number
+(sys1::LTISystem, sys2::Number) = +(promote(sys1, sys2)...)
-(sys1::LTISystem, sys2::Number) = -(promote(sys1, sys2)...)
*(sys1::LTISystem, sys2::Number) = *(promote(sys1, sys2)...)
/(sys1::LTISystem, sys2::Number) = /(promote(sys1, sys2)...)

+(sys1::Number, sys2::LTISystem) = +(promote(sys1, sys2)...)
-(sys1::Number, sys2::LTISystem) = -(promote(sys1, sys2)...)
*(sys1::Number, sys2::LTISystem) = *(promote(sys1, sys2)...)
/(sys1::Number, sys2::LTISystem) = /(promote(sys1, sys2)...)

# Fallback Matrix
+(sys1::LTISystem, sys2::AbstractMatrix) = +(promote(sys1, sys2)...)
-(sys1::LTISystem, sys2::AbstractMatrix) = -(promote(sys1, sys2)...)
*(sys1::LTISystem, sys2::AbstractMatrix) = *(promote(sys1, sys2)...)
/(sys1::LTISystem, sys2::AbstractMatrix) = /(promote(sys1, sys2)...)

+(sys1::AbstractMatrix, sys2::LTISystem) = +(promote(sys1, sys2)...)
-(sys1::AbstractMatrix, sys2::LTISystem) = -(promote(sys1, sys2)...)
*(sys1::AbstractMatrix, sys2::LTISystem) = *(promote(sys1, sys2)...)
/(sys1::AbstractMatrix, sys2::LTISystem) = /(promote(sys1, sys2)...)

Base.:^(sys::LTISystem, p::Integer) = Base.power_by_squaring(sys, p)

# TODO We should have proper fallbacks for matrices too

feedback(sys1::Union{LTISystem,Number,AbstractMatrix{<:Number}},
         sys2::Union{LTISystem,Number,AbstractMatrix{<:Number}}; kwargs...) = feedback(promote(sys1, sys2)...; kwargs...)

Base.inv(G::LTISystem) = 1/G

"""`issiso(sys)`

Returns `true` if `sys` is SISO, else returns `false`."""
function issiso(sys::LTISystem)
    return ninputs(sys) == 1 && noutputs(sys) == 1
end

"""`timeevol(sys::LTISystem)`
Get the timeevol::TimeEvolution from system `sys`, usually sys.timeevol """
timeevol(sys::LTISystem) = sys.timeevol

"""`iscontinuous(sys)`

Returns `true` for a continuous-time system `sys`, else returns `false`."""
iscontinuous(sys::LTISystem) = timeevol(sys) isa Continuous
"""`isdiscrete(sys)`

Returns `true` for a discrete-time system `sys`, else returns `false`."""
isdiscrete(sys::LTISystem) = timeevol(sys) isa Discrete

function Base.getproperty(sys::LTISystem, s::Symbol)
    if s === :Ts
        # if !isdiscrete(sys) # NOTE this line seems to be breaking inference of isdiscrete (is there a test for this?)
        if isdiscrete(sys)
            return timeevol(sys).Ts
        else
            @warn "Getting time 0.0 for non-discrete systems is deprecated. Check `isdiscrete` before trying to access time."
            return 0.0
        end
    else
        return getfield(sys, s)
    end
end

Base.propertynames(sys::LTISystem, private::Bool=false) =
    (fieldnames(typeof(sys))..., (isdiscrete(sys) ? (:Ts,) : ())...)

common_timeevol(systems::LTISystem...) = common_timeevol(timeevol(sys) for sys in systems)

"""`isstable(sys)`

Returns `true` if `sys` is stable, else returns `false`."""
function isstable(sys::LTISystem)
    if iscontinuous(sys)
        if any(real.(poles(sys)).>=0)
            return false
        end
    else
        if any(abs.(poles(sys)).>=1)
            return false
        end
    end
    return true
end

# Fallback since LTISystem not AbstractArray
Base.size(sys::LTISystem, i::Integer) = size(sys)[i]


text:
function print_param_info(p::Param)

    #Grid information
    tmax        = p.tmax
    t_save      = p.out_every_t
    deltat      = p.deltat

    xmin        = p.xmin
    xmax        = p.xmax
    Nx          = p.xnodes

    ymin        = p.ymin
    ymax        = p.ymax
    Ny          = p.ynodes

    dx      = (xmax - xmin)/(Nx-1)
    dy      = (ymax - ymin)/(Ny-1)

    if(p.Boundaries == :periodic )
        dx      = (xmax - xmin)/(Nx)
        dy      = (ymax - ymin)/(Ny)

    end

    #Print information
    println("---------------------------------------------")
    println("Parameter Reader:")
    println("---------------------------------------------")
    println("Grid Information:")
    println("x  = [" , xmin , "," , xmax ,"]" )
    println("y  = [" , ymin , "," , ymax ,"]" )
    println("dx = "  , dx)
    println("dy = "  , dy)
    println("---------------------------------------------")
    println("Time Information:")
    println("t  = [" , 0.0 , "," , tmax ,"]" )
    println("dt = "  , deltat)
    println("---------------------------------------------")
    println("Initial Conditions:")
    println("A0 = " , p.A0)
    println("σ  = " , p.σ)
    println("r0 = " , p.r0)
    println("ω  = " , p.ω)
    println("m  = " , p.m)
    println("---------------------------------------------")
    println("Field and System:")
    println("μ           = " , p.μ)
    println("dissipation = " , p.dissipation)

    if( p.dissipation)

    println("Orbit r     = " , p.R_orbit)
    println("BH1 r       = " , p.R1)
    println("BH1 α       = " , p.alpha1)
    println("BH2 r       = " , p.R2)
    println("BH2 α       = " , p.alpha2)

    println("Ω           = " , p.Ω)

    end

    println("Boundaries  = " , p.Boundaries)
    println("---------------------------------------------")
    println("saving to file  = " ,  string(p.folder ,"/", p.fname , ".h5"))
    println("---------------------------------------------")

end


text:
#=

This file contains helper functions related to downloading/using artifacts

=#

""" Shorthand for @artifact_str """
function artifact_helper(name::AbstractString)

    joinpath(ensure_artifact_installed(name, find_artifacts_toml(@__DIR__),
                                       quiet_download = true), name) * ".par"

end

"""
    $(FUNCTIONNAME)(molecule::AbstractString, database::AbstractString = "hitran")

Given a molecule name and a database name, retrieve the transition file
"""
function artifact(molecule::AbstractString,
                  database::AbstractString = "hitran")
    @assert (database == "hitran") "This package currently only supports the hitran database!"
    @assert molecule in Absorption.mol_names "Not a supported molecule!"
    mol_id = Absorption.mol_number(molecule)
    return artifact_helper("hitran_molec_id_$(mol_id)_$(molecule)")
end


text:
let sys = read("data/dense.json", DenseSystem)

    @test isapprox(sys, DenseSystem(test_params_flat...))

    buf = IOBuffer()
    write(buf, sys)
    seekstart(buf)
    @test read(buf, DenseSystem) == sys
end

let sys = read("data/diagonal.json", DenseSystem)

    @test isapprox(sys, DenseSystem(test_params_flat_diag...))

    buf = IOBuffer()
    write(buf, sys)
    seekstart(buf)
    @test read(buf, DenseSystem) == sys
end

let sys = read("data/diagonal.json", DiagonalSystem)

    @test isapprox(sys, DiagonalSystem(test_params_flat_diag...))

    buf = IOBuffer()
    write(buf, sys)
    seekstart(buf)
    @test read(buf, DiagonalSystem) == sys
end

@test_throws SurfaceCouplingException read("data/dense.json", DiagonalSystem)


text:
"""
Stimulus

A wrapper class and API for handling stimuli without needing to worry about the
specific implementation for that experiment. Basically, I want to be able to
call `compute_STRFs(raster, stimulus)` regardless of the data source.
"""
module GrayScaleStimuli

using JLD
import Base: show, isequal, hash
include("metadata.jl")

################################################################################
#### Abstract type to use as a placeholder/to get all docstrings
################################################################################
""""
    AbstractStimulus

Abstract type used as placeholder.
"""
abstract type AbstractStimulus end

################################################################################
#### Miscellaneous functions and constants
################################################################################

_NI(m) = error("Not Implemented: $m")

"""
    frame_size(S)

The size of the image, in pixels. Returns a tuple `(w,h)`; each frame of the image is
displayed on screen as an image that many pixels in size.
"""
frame_size(S::AbstractStimulus) = _NI("frame_size($(typeof(S)))")

"""
    frame_time(S)

The duration of a single frame, in seconds
"""
frame_time(S::AbstractStimulus) = _NI("frame_time($(typeof(S)))")

"""
    frame_rate(S::AbstractStimulus)

Convenience method returns `1/frame_time(S)`. Each concrete subtype of
`AbstractStimulus` should implement `frame_time`.
"""
frame_rate(S::AbstractStimulus) = 1.0 / frame_time(S)

"""
    n_frames(S)

The number of frames
"""
n_frames(S::AbstractStimulus) = _NI("n_frames($(typeof(S)))")

"""
    resolution(S::AbstractStimulus)

The resolution of the individual frames
"""
resolution(S::AbstractStimulus) = _NI("resolution($(typeof(S)))")

"""
    patch_size(S::AbstractStimulus)

The size, in pixels, of the individual patches of constant value. For any stimulus `S`, we
should have `frame_size(S) == resolution(S) .* patch_size(S)`
"""
patch_size(S::AbstractStimulus) = _NI("patch_size($(typeof(S)))")

"""
    matrix_form([Float64], S::AbstractStimulus, frames=1:n_frames(S))

A minimal representation of `S` as a `Matrix` of the specified type. The `i,j`th entry is the `i`th pixel of the `j`th frame.
"""
matrix_form(::Type{T}, S::AbstractStimulus, frames=1:n_frames(S)) where T = _NI("matrix_form($T, $(tpyoef(S)), $frames)")
matrix_form(S::AbstractStimulus, frames=1:n_frames(S)) = matrix_form(Float64, S, frames)

"""
    _compute_STRFs(spike_hist, stim; kwargs...)

Generic implementation of STRF computation. Uses `matrix_form(stim)` to perform
computations. Returns a matrix.

Life will be easy if we assume that the bins of the histogram are 1-to-1 with
the frames of the stimulus. So, implement concrete subtypes of
`AbstractStimulus` with that in mind.

For keyword arguments: specifying `window_length_frames` overrides any value
computed from `window_length_s`
"""
function _compute_STRFs(spike_hist::Matrix{Float64}, stim::AbstractStimulus; kwargs...)
    # preprocessing stuff: get the window length
    dkwargs = Dict(kwargs)
    window_length_s = get(dkwargs, :window_length_s, 0.5)
    window_length_frames = round(Int, window_length_s / frame_time(stim))
    if haskey(dkwargs, :window_length_frames)
        window_length_frames = dkwargs[:window_length_frames]
    end

    # get the matrix form of the stimulus
    stimulus = matrix_form(stim)

    N_bins, N_cells = size(spike_hist)
    N_bits_per_frame, N_frames = size(stimulus)
    if N_bins != N_frames
        error("_compute_STRFs: mismatched histogram and matrix form (N_bins = $N_bins, N_frames = $N_frames)")
    end

    avg_spikes_per_frame = zeros(spike_hist)
    for j = 1:N_cells
        avg_spikes_per_frame[window_length_frames:end,j] = spike_hist[window_length_frames:end,j] / sum(spike_hist[window_length_frames:end-1,j])
    end

    RFs = zeros(N_bits_per_frame, N_cells, window_length_frames)
    for k = 1:window_length_frames
        RFs[:, :, window_length_frames - k + 1] = stimulus * avg_spikes_per_frame
        avg_spikes_per_frame = circshift(avg_spikes_per_frame, [-1,0])
    end
    return RFs
end

################################################################################
#### Include definitions of concrete types
################################################################################
include("GrayScaleStimulus.jl")
include("optimizers.jl")
include("STRF_response.jl")
include("animation.jl")
include("fileio.jl")
# include("CRCNS/CRCNS_Stimulus.jl")
# include("animation.jl")
# include("Ganmor/Ganmor_stimulus.jl")

export show, isequal, hash, savestimulus, loadstimulus,
       frame_size, frame_time, frame_rate, n_frames, resolution, patch_size, matrix_form, frame_image,
       time_to_index, index_to_time,
       GrayScaleStimulus, #CRCNS_Stimulus, Ganmor_Stimulus,
       compute_STRFs, STRF_response, scale_response,
       imagemean, matrixcovariance,
       animated_gif

end


text:
"""
Recursively merge two nested dictionaries
"""
function merge_rec(a::AbstractDict, b::AbstractDict)
    d = Dict()
    for (k, v) in a
        if haskey(b, k)
            d[k] = merge_rec(v, b[k])
        else
            d[k] = v
        end
    end

    for (k, v) in b
        if !haskey(a, k)
            d[k] = v
        end
    end

    return d
end

merge_rec(a, b) = b

"""
Read in a yaml or folder of yamls into a dictionary.
"""
function build_yaml(path::AbstractString)
    yml = Dict{String,Any}()
    if isdir(path)
        for file in readdir(path; join = true)
            if endswith(file, r"\.ya?ml")
                this_yml = YAML.load_file(file)
                merge!(yml, this_yml)
            end
        end
    else
        @assert endswith(path, r"\.ya?ml") "Must provide a yaml file"
        yml = YAML.load_file(path)
    end

    return yml
end

"""
Get the head and tail tuples from args
"""
headtail(a, b...) = (a, b)

"""
Compare the original params to what was parsed and print warnings for any original
params that are unused in the final parsed parasms.
"""
function warn_unused_params(parsed::AbstractDict, params::AbstractDict; key_path = "")
    count = 0

    for (k, v) in params
        kp = "$key_path.$k"
        if haskey(parsed, k)
            count += warn_unused_params(parsed[k], params[k]; key_path = kp)
        else
            @warn "[$kp] is unused"
            count += 1
        end
    end

    return count
end

# neither is a dict
warn_unused_params(parsed, params; key_path = "") = 0

# one is a dict
function warn_unused_params(parsed, params::AbstractDict; key_path = "")
    @warn "[$key_path] has unused params: $params"
    return 1
end
function warn_unused_params(parsed::AbstractDict, params; key_path = "")
    @warn "[$key_path] has sub-keys, got unused params: $params"
    return 1
end

"""
Get the keys of a dict, or the items in a list.
"""
keys_or_vals(d::AbstractDict) = collect(keys(d))
keys_or_vals(a::AbstractArray) = a


text:
"""
  function calculate_D(sMMProblem::SMMProblem, theta0::Array{Float64,1}; method::Symbol = :central)

Function to calculate the jacobian of the simulated moments.
If the simulation is long enough, this provideds a good approximation for
the expected value of the jacobian. The output is the "D" matrix in the terminology
of Gouriéroux and Monfort (1996).
"""
function calculate_D(sMMProblem::SMMProblem, theta0::Array{Float64,1}; method::Symbol = :central)

  Calculus.jacobian(x -> sMMProblem.simulate_empirical_moments_array(x), theta0, method)

end

"""
  calculate_Avar!(sMMProblem::SMMProblem, theta0::Array{Float64,1}; method::Symbol = :central)

Calculate the asymptotic variance of the SMM estimator using the sandwich formula.
This function assume that your model is simulating time series.
As a result, the asymptotic variance of the SMM estimator is (1+tau) times
the asymptotic variance of the GMM one, where tau is the ratio of the sample size
to the size of the simulated sample. See Duffie and Singleton (1993) and Gouriéroux and Monfort (1996).
"""
function calculate_Avar!(sMMProblem::SMMProblem, theta0::Array{Float64,1}, tData::Int64, tSimulation::Int64; method::Symbol = :central)

  # Safety Checks
  if sMMProblem.Sigma0 == Array{Float64}(0,0)
    Base.error("Please initialize sMMProblem.Sigma0 using the function set_Sigma0!.")
  end

  tau = tData/tSimulation

  # Jacobian matrix
  D = calculate_D(sMMProblem, theta0, method = method)

  # Weigthing matrix (only diagonal matrix for the moment)
  W = diagm([1/sMMProblem.empiricalMoments[k][2] for k in keys(sMMProblem.empiricalMoments)])

  Sigma1 = transpose(D)*W*D

  Sigma2 = transpose(D)*W*sMMProblem.Sigma0*W*D

  inv_Sigma1 = inv(Sigma1)

  # sandwich formula, where tau captures the extra noise introduced by simulation,
  # compared to the usual GMM estimator.
  # I use the function Symmetric, because otherwise a test issymmetric()
  # on the asymptotic variance will return "false" because of numerical approximations
  sMMProblem.Avar = Symmetric((1+tau)*inv_Sigma1*Sigma2*inv_Sigma1)

end

"""
  calculate_se(sMMProblem::SMMProblem, tData::Int64, tSimulation::Int64)

Function to calculate the standard error associated to the the ith parameter,
respecting the ordering given by the ordered dictionary sMMProblem.priors
"""
function calculate_se(sMMProblem::SMMProblem, tData::Int64, i::Int64)

  # Safety Checks
  if sMMProblem.Avar == Array{Float64}(0,0)
    Base.error("Please caclulate the asymptotic variance using the function calculate_Avar!.")
  end

  sqrt((1/tData)*sMMProblem.Avar[i,i])

end

"""
  calculate_t(sMMProblem::SMMProblem, tData::Int64, tSimulation::Int64)

Function to calculate the t-statistic associated to following test:
H0: theta_i = 0
H1: theta_i != 0
The ordering of parameters is the one given by the ordered dictionary sMMProblem.priors
"""
function calculate_t(sMMProblem::SMMProblem, theta0::Array{Float64,1}, tData::Int64, i::Int64)

  # Safety Checks
  if sMMProblem.Avar == Array{Float64}(0,0)
    Base.error("Please caclulate the asymptotic variance using the function calculate_Avar!.")
  end

  theta0[i]/calculate_se(sMMProblem, tData, i)

end

"""
  calculate_pvalue(sMMProblem::SMMProblem, tData::Int64, tSimulation::Int64)

Function to calculate the p-value associated to following test:
H0: theta_i = 0
H1: theta_i != 0
The ordering of parameters is the one given by the ordered dictionary sMMProblem.priors
"""
function calculate_pvalue(sMMProblem::SMMProblem, theta0::Array{Float64,1}, tData::Int64, i::Int64)

  # Safety Checks
  if sMMProblem.Avar == Array{Float64}(0,0)
    Base.error("Please caclulate the asymptotic variance using the function calculate_Avar!.")
  end

  t =  calculate_t(sMMProblem, theta0, tData, i)

  # asymptotically normally distributed N(0,1)
  d = Normal(0,1)

  # p-value  = 2 times area to the right of t
  pvalue = 2*(1.0 - cdf(d, t))

  return pvalue

end

"""
  calculate_CI(sMMProblem::SMMProblem, theta0::Array{Float64,1}, tData::Int64, i::Int64)

Function to calculate an alpha confidence interval for the ith parameter.
"""
function calculate_CI(sMMProblem::SMMProblem, theta0::Array{Float64,1}, tData::Int64, i::Int64, alpha::Float64)

  # Safety Checks
  if sMMProblem.Avar == Array{Float64}(0,0)
    Base.error("Please caclulate the asymptotic variance using the function calculate_Avar!.")
  end

  # standard error
  se = calculate_se(sMMProblem, tData, i)

  # asymptotically normally distributed N(0,1)
  d = Normal(0,1)

  critical_value = - quantile(d, alpha)

  #return lower and upper bound of the confidence interval
  return theta0[i] - se*critical_value, theta0[i] + se*critical_value

end

function summary_table(sMMProblem::SMMProblem, theta0::Array{Float64,1}, tData::Int64, alpha::Float64)

  # Safety Checks
  if sMMProblem.Avar == Array{Float64}(0,0)
    Base.error("Please caclulate the asymptotic variance using the function calculate_Avar!.")
  end

  df = DataFrame(Estimate = Float64[], StdError = Float64[], tValue = Float64[], pValue = Float64[], ConfIntervalLower = Float64[], ConfIntervalUpper = Float64[])

  for i = 1:length(theta0)

    se = calculate_se(sMMProblem, tData, i)
    t = calculate_t(sMMProblem, theta0, tData, i)
    p = calculate_pvalue(sMMProblem, theta0, tData, i)
    CI_lower, CI_upper = calculate_CI(sMMProblem, theta0, tData, i, alpha)

    push!(df, [theta0[i], se, t, p, CI_lower, CI_upper])

  end

  return df

end


text:
W, H, x, y, r = parse.(Int, split(readline()))
if (0 <= x-r && x+r <= W) && (0 <= y-r && y+r <= H)
    println("Yes")
else
    println("No")
end


text:
using SweepOperator

# necessary to support 0.4 and 0.5
if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

# setup
n, p = 1000, 10
x = randn(n, p)
xtx = x'x

@testset "Sweep One By One" begin
    A = deepcopy(xtx)
    B = deepcopy(xtx)
    for j in 1:p
        sweep!(A, j)
        sweep!(A, j, true)
    end
    @test A ≈ B

    A = deepcopy(xtx)
    B = deepcopy(xtx)
    for j in 1:p
        sweep!(A, j, true)
    end
    for j in 1:p
        sweep!(A, j, true)
    end
    @test A ≈ B
end

@testset "Sweep All" begin
    A = deepcopy(xtx)
    B = deepcopy(xtx)
    sweep!(A, 1:p)
    sweep!(A, 1:p, true)
    @test A ≈ B
end

@testset "Linear Regression" begin
    x = randn(n, p)
    y = x * collect(1.:p) + randn(n)
    xy = [x y]
    xytxy = xy'xy
    sweep!(xytxy, 1:p)
    @test xytxy[1:p, end] ≈ x\y
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ c0fb9740-1f5f-11eb-36f0-f98f57070f9d
begin
	using NeXLCore
	using PlutoUI, DataFrames
	md"""
	## Mass Fraction to Atomic Fraction

	Compute various representations of a material's composition including uncertainty components.
	"""
end

# ╔═╡ df8df9ee-1f5f-11eb-2b7b-3f998003064d
md"""
Enter the composition of the material in this text field using the format:

    element1, mass_fraction1[, dmass_fraction1 [, atomic_weight1]]
    element2, mass_fraction2[, dmass_fraction2 [, atomic_weight2]]
    ...
    elementn, mass_fractionn[, dmass_fractionn [, atomic_weightn]]

where items in [] are optional.

$(@bind matstr TextField( (80, 10) ))

Material name: $(@bind namstr TextField(default="Unknown"))
"""

# ╔═╡ d7ec2bc0-1f61-11eb-000f-ad121d90310a
begin
	function parse_matstr(nm, str)
		res = Dict{Label, UncertainValue}()
		for line in split(str, "\n")
			items = split(line, ",")
			elm, qty, dqty = missing, 0.0, 0.0
			try
				elm = length(items)>=1 ? element(strip(items[1])) : missing
			catch
				elm = missing
			end
			try
				qty = length(items)>=2 ? parse(Float64, strip(items[2])) : 0.0
			catch
				qty = 0.0
			end
			try
			 	dqty = length(items)>=3 ? parse(Float64, strip(items[3])) : 0.0
			catch
				dqty = 0.0
			end
			if (!ismissing(elm)) && ((qty ≠ 0.0) || (dqty ≠ 0.0))
				res[MassFractionLabel(nm,elm)] = uv(qty, dqty)
				aw = a(elm)
				try
					aw = length(items)>=4 ? parse(Float64, strip(items[4])) : a(elm)
				catch
					aw = a(elm)
				end
				res[AtomicWeightLabel(nm,elm)] = aw
			end
		end
		return length(res) > 0 ? mf2comp(nm, uvs(res)) : nothing
	end
	nothing
end

# ╔═╡ d1faad00-1f60-11eb-2f22-d197787a36ec
begin
	mat = parse_matstr(namstr, matstr)
	!isnothing(mat) ? sort(asa(DataFrame, mat, false),:Variable) : nothing
end

# ╔═╡ Cell order:
# ╟─c0fb9740-1f5f-11eb-36f0-f98f57070f9d
# ╟─df8df9ee-1f5f-11eb-2b7b-3f998003064d
# ╟─d7ec2bc0-1f61-11eb-000f-ad121d90310a
# ╟─d1faad00-1f60-11eb-2f22-d197787a36ec


text:
### Standard Metropolis-Hastings Monte Carlo runner
immutable SMHRunner <: AbstractMHRunner
    numburnin::Int
    numiterations::Int
    policy::MHRuntimePolicy #policies determining runtime behaviour of the MCMC

    function SMHRunner(nburnin::Int,niterations::Int,policy::MHRuntimePolicy)
        @assert 0 <= nburnin && 0 < niterations
        new(nburnin,niterations,policy)
    end
end

###Factory functions
_runner(::Type{Val{:standard}},policy::MHRuntimePolicy,niterations::Int;numburnin =0) = SMHRunner(numburnin,niterations,policy)

function run!(runner_::SMHRunner,model_::AbstractModel,sampler_::AbstractSampler,tuner_::AbstractTuner)
    samplerstate_ = samplerstate(sampler_,1,runner_.policy.sampletype,runner_.policy.calculationtype)
    indicator_,tunerstate_,chain_ = createcommon(runner_,tuner_,numparas(model_),1,1)
    tic()
    burnin!(runner_,model_,sampler_,samplerstate_,tuner_,tunerstate_,indicator_,chain_)
    main!(runner_,model_,sampler_,samplerstate_,tuner_,tunerstate_,indicator_,chain_)
    chain_.runtime = toq()
    chain_
end

function burnin!(runner_::SMHRunner,model_::AbstractModel,sampler_::AbstractSampler,samplerstate_::AbstractSamplerState,
                 tuner_::AbstractTuner,tunerstate_::AbstractTunerState,indicator_::AbstractIndicatorMatrix,chain_::AbstractChain)
    storeduring = _storeduring(:burnin,runner_.policy)
    initialize!(runner_,model_,samplerstate_,chain_,storeduring)
    for i=1:runner_.numburnin
        iterate!(runner_,model_,samplerstate_,indicator_,chain_,storeduring)
        accepted!(tunerstate_,indicator_)
        needstuning(tuner_,i)?tune!(runner_,samplerstate_,tuner_,tunerstate_):nothing
        preparenext!(runner_,samplerstate_,indicator_) #prepare for the next iteration
    end
end

function main!(runner_::SMHRunner,model_::AbstractModel,sampler_::AbstractSampler,samplerstate_::AbstractSamplerState,
               tuner_::AbstractTuner,tunerstate_::AbstractTunerState,indicator_::AbstractIndicatorMatrix,chain_::AbstractChain)
    storeduring = _storeduring(:main,runner_.policy)
    for i=1:runner_.numiterations
        iterate!(runner_,model_,samplerstate_,indicator_,chain_,storeduring)
        needstuning(tuner_,i)?println("Iteration $i/$(runner_.numiterations)"):nothing
        preparenext!(runner_,samplerstate_,indicator_)
    end
end

function iterate!(runner_::SMHRunner,model_::AbstractModel,samplerstate_::AbstractSamplerState,
                  indicator_::AbstractIndicatorMatrix,chain_::AbstractChain,storeduring::Bool)
    propose!(samplerstate_)
    geometry!(model_,proposals(samplerstate_))
    a = acceptance!(samplerstate_)
    transitionprobability!(indicator_,a)
    sampleindicator!(indicator_)
    storeduring?store!(runner_,samplerstate_,indicator_,chain_):nothing
end

function store!(runner_::SMHRunner,samplerstate_::AbstractSamplerState,
                indicator_::AbstractIndicatorMatrix,chain_::AbstractChain)
    if indicatorsamples(indicator_)[2] == 2
        store!(chain_,from(samplerstate_),1)
    else
        store!(chain_,proposals(samplerstate_),1)
    end
    accepted!(chain_,indicator_)
end

preparenext!(runner_::SMHRunner,samplerstate_::AbstractSamplerState,indicator_::AbstractIndicatorMatrix) = prepare!(samplerstate_,indicatorsamples(indicator_)[2]!=2)

function tune!(runner_::SMHRunner,samplerstate_::AbstractSamplerState,tuner_::AbstractTuner,tunerstate_::AbstractTunerState)
    tvals = tune(tuner_,tunerstate_)
    showstep(tuner_,tunerstate_)
    tune!(samplerstate_,tvals...)
    nextindex!(tunerstate_)
end

function show(io::IO,r::SMHRunner)
    println(io,"Standard Metropolis-Hastings runner with:")
    println(io," numburnin: ",r.numburnin)
    println(io," numiterations: ",r.numiterations)
    print(io," policy: ")
    show(io,r.policy)
    println(io)
    nothing
end


text:
export little_group_elements
export little_group
export little_symmetry
export little_symmetry_iso

"""
    little_group_elements(tsym, psym)

Return elements of `psym` compatible with the translation symmetry `tsym`.
"""
function little_group_elements(tsym::FiniteTranslationSymmetry, psym::PointSymmetry)::Vector{Int}
    return Int[
        i_elem
            for (i_elem, elem) in enumerate(elements(psym))
                if iscompatible(tsym, elem)
    ]
end

"""
    little_group_elements(tsym, tsym_irrep_index, psym)

Return elements of `psym` compatible with the translation symmetry `tsym` at
irrep `tsym_irrep_index`
"""
function little_group_elements(
    tsym::FiniteTranslationSymmetry,
    tsym_irrep_index::Integer,
    psym::PointSymmetry,
)::Vector{Int}
    return Int[
        i_elem
            for (i_elem, pop) in enumerate(elements(psym))
                if iscompatible(tsym, tsym_irrep_index, pop)
    ]
end

"""
    little_group(tsym, tsym_irrep_index, psym)
"""
function little_group(
    tsym::FiniteTranslationSymmetry,
    tsym_irrep_index::Integer,
    psym::PointSymmetry,
)::FiniteGroup
    lg_elems = little_group_elements(tsym, tsym_irrep_index, psym)
    return little_group(tsym, psym, lg_elems)
end

"""
    little_group(tsym, psym, elements)

Generate a little group with given elements.
The elements of the little group, which may be sparse, are compressed into consecutive integers.
"""
function little_group(
    tsym::FiniteTranslationSymmetry,
    psym::PointSymmetry,
    elements::AbstractVector{<:Integer},
)
    element_lookup = Dict(x=>i for (i, x) in enumerate(elements))
    ord_group = length(elements)
    mtab = zeros(Int, (ord_group, ord_group))
    for i in 1:ord_group, j in 1:ord_group
        mtab[i,j] = element_lookup[group_product(psym.group, elements[i], elements[j])]
    end
    return FiniteGroup(mtab)
end

"""
    little_symmetry(tsym, psym, lg_elements)

Generate a little symmetry.
First generate a little group with `lg_elements`, look up `PointSymmetryDatabase` using
the names of the elements, and then checks for group isomorphism.

# Arguments
* `tsym::FiniteTranslationSymmetry`: translation symmetry
* `psym::PointSymmetry`: point symmetry
* `lg_elements::AbstractVector{<:Integer}`: indices of the elements of the little group
"""
function little_symmetry(
    tsym::FiniteTranslationSymmetry,
    psym::PointSymmetry,
    lg_elements::AbstractVector{<:Integer},
)
    (lg_raw, lg_matrep_raw, lg_element_names_raw) = let
        lg_raw = little_group(tsym, psym, lg_elements)
        lg_matrep_raw = psym.matrix_representations[lg_elements]
        lg_element_names_raw = psym.element_names[lg_elements]
        (lg_raw, lg_matrep_raw, lg_element_names_raw)
    end

    psym2 = PointSymmetryDatabase.find(lg_element_names_raw)
    ϕ = group_isomorphism(lg_raw, psym2.group)
    isnothing(ϕ) && error("Group not isomorphic")

    lg_matrep = lg_matrep_raw[ϕ]
    lg_element_names = lg_element_names_raw[ϕ]

    return PointSymmetry(
        psym2.group,
        psym2.generators,
        psym2.conjugacy_classes,
        psym2.character_table,
        psym2.irreps,
        lg_element_names,
        lg_matrep,
        psym2.hermann_mauguin,
        psym2.schoenflies,
    )
end

"""
    little_symmetry(tsym, psym)
"""
function little_symmetry(tsym::FiniteTranslationSymmetry, psym::PointSymmetry)
    return little_symmetry(tsym, psym, little_group_elements(tsym, psym))
end

"""
    little_symmetry(tsym, tsym_irrep_index, psym)
"""
function little_symmetry(tsym::FiniteTranslationSymmetry, tsym_irrep::Integer, psym::PointSymmetry)
    if !iscompatible(tsym, psym)
        throw(ArgumentError("translation and point symmetries are not compatible"))
    end

    tsym_irrep == 1 && return little_symmetry(tsym, psym)

    (lg_raw, lg_matrep_raw, lg_element_names_raw) = let
        lg_elements = little_group_elements(tsym, tsym_irrep, psym)
        lg_raw = little_group(tsym, psym, lg_elements)
        lg_matrep_raw = psym.matrix_representations[lg_elements]
        lg_element_names_raw = psym.element_names[lg_elements]
        (lg_raw, lg_matrep_raw, lg_element_names_raw)
    end

    psym2 = PointSymmetryDatabase.find(lg_element_names_raw)
    ϕ = group_isomorphism(lg_raw, psym2.group)
    isnothing(ϕ) && error("Group not isomorphic")

    lg_matrep = lg_matrep_raw[ϕ]
    lg_element_names = lg_element_names_raw[ϕ]

    return PointSymmetry(
        psym2.group,
        psym2.generators,
        psym2.conjugacy_classes,
        psym2.character_table,
        psym2.irreps,
        lg_element_names,
        lg_matrep,
        psym2.hermann_mauguin,
        psym2.schoenflies,
    )
end

# COV_EXCL_START

"""
    little_symmetry_iso(tsym, tsym_irrep_index, psym)

Find little symmetry using group isomorphism
"""
function little_symmetry_iso(
    tsym::FiniteTranslationSymmetry,
    tsym_irrep_index::Integer,
    psym::PointSymmetry,
)
    @warn "Function little_symmetry_iso deprecated. Correctness not guaranteed." maxlog=1
    if !iscompatible(tsym, psym)
        throw(ArgumentError("translation and point symmetries are not compatible"))
    end
    tsym_irrep_index == 1 && return little_symmetry(tsym, psym)
    (lg_matrep, lg_element_names, lg_irrep) = let
        lg_elements = little_group_elements(tsym, tsym_irrep_index, psym)

        lg_raw = little_group(tsym, psym, lg_elements)
        lg_matrep_raw = psym.matrix_representations[lg_elements]
        lg_element_names_raw = psym.element_names[lg_elements]

        (lg_irrep, ϕ) = IrrepDatabase.find(lg_raw)

        lg_matrep = lg_matrep_raw[ϕ]
        lg_element_names = lg_element_names_raw[ϕ]
        (lg_matrep, lg_element_names, lg_irrep)
    end

    generators = minimal_generating_set(lg_irrep.group)
    hermann_mauguin = join(lg_element_names[generators])
    schoenflies = "unknown"

    simple_element_names = sort(simplify_name.(lg_element_names))
    for i in 1:32
        psym = PointSymmetryDatabase.get(i)
        if ( sort(simplify_name.(psym.element_names)) == simple_element_names
             && !isnothing(group_isomorphism(lg_irrep.group, psym.group)) )
            hermann_mauguin = psym.hermann_mauguin
            break
        end
    end

    return PointSymmetry(
        lg_irrep.group,
        generators,
        lg_irrep.conjugacy_classes,
        lg_irrep.character_table,
        lg_irrep.irreps,
        lg_element_names,
        lg_matrep,
        hermann_mauguin,
        schoenflies,
    )
end
# COV_EXCL_STOP


text:
module EnsembleSampler

using Distributed

using Printf

using ..Acor

## Store ensembles as (ndim, nwalkers) arrays.

"""
    propose(ps, qs)

Return the proposed stretch-move updates to the points `ps` treating
the `qs` as fixed.
"""
function propose(ps::Array{Float64, 2}, qs::Array{Float64, 2})
    n = size(ps, 2)

    zs = exp.(log(0.5) .+ rand(Float64, (n,)).*(log(2.0) - log(0.5)))
    inds = rand(1:n, (n,))

    ps_out = zeros(size(ps))
    for i in 1:n
        ps_out[:,i] = qs[:,inds[i]] + zs[i]*(ps[:,i] - qs[:,inds[i]])
    end
    ps_out, zs
end

"""
    lnprobs(xs, lnprobfn)

Return an array of log-probabilities of the points `xs` under the
log-probability function `lnprobfn`.
"""
function lnprobs(xs::Array{Float64, 2}, lnprobfn)
    if length(workers()) > 1
        xseq = Any[xs[:,i] for i=1:size(xs,2)]
        Array{Float64,1}(pmap(lnprobfn, xseq, batch_size=div(size(xs,2),(4*length(workers())))))
    else
        Float64[lnprobfn(xs[:,i]) for i in 1:size(xs, 2)]
    end
end

"""
    update_half(ensemble, lnprobs, lnprobfn, half)

Return the updated `(ensemble, lnprobs)` from updating either the odd
(`half % 2 == 0`) or even (`half % 2 == 0`) indices of the `(ndim,
nwalkers)` array of ensemble.
"""
function update_half(ensemble::Array{Float64, 2}, lnprob::Array{Float64, 1}, lnprobfn, half)
    n = size(ensemble, 2)
    nd = size(ensemble, 1)

    if n % 2 != 0
        error("update_half requires even number of walkers")
    end

    n_half = div(n, 2)

    if half % 2 == 0
        ps = ensemble[:, 1:n_half]
        qs = ensemble[:, n_half+1:end]
        lnps = lnprob[1:n_half]
    else
        ps = ensemble[:, n_half+1:end]
        qs = ensemble[:, 1:n_half]
        lnps = lnprob[n_half+1:end]
    end

    ps_new, zs = propose(ps, qs)
    lnp_new = lnprobs(ps_new, lnprobfn)

    lnpacc = lnp_new - lnps + nd*log.(zs)

    outhalf = zeros((nd, n_half))
    lnpouthalf = zeros(n_half)
    for i in 1:n_half
        if lnpacc[i] > 0 || log(rand()) < lnpacc[i]
            outhalf[:,i] = ps_new[:,i]
            lnpouthalf[i] = lnp_new[i]
        else
            outhalf[:,i] = ps[:,i]
            lnpouthalf[i] = lnps[i]
        end
    end

    out = zeros((nd, n))
    lnpout = zeros(n)
    if half % 2 == 0
        out[:,1:n_half] = outhalf
        out[:,n_half+1:end] = ensemble[:,n_half+1:end]

        lnpout[1:n_half] = lnpouthalf
        lnpout[n_half+1:end] = lnprob[n_half+1:end]
    else
        out[:,1:n_half] = ensemble[:,1:n_half]
        out[:,n_half+1:end] = outhalf

        lnpout[1:n_half] = lnprob[1:n_half]
        lnpout[n_half+1:end] = lnpouthalf
    end

    out, lnpout
end

"""
   update(ensemble, lnprob, lnprobfn)

Return the updated `(ensemble, lnprob)` from a MCMC step.
"""
function update(ensemble::Array{Float64, 2}, lnprob::Array{Float64,1}, lnprobfn)
    e, lp = update_half(ensemble, lnprob, lnprobfn, 0)
    update_half(e, lp, lnprobfn, 1)
end

"""
    run_mcmc(ensemble, lnprob, lnprobfn, steps; thin=1)

Return the pair `(chain, chainlnprob)` representing `steps` of MCMC
updates (thinned by `thin`).

`chain` will have shape `(ndim, nwalkers, nsave)` where `nsave =
div(steps, thin)`.  `chainlnprob` will have shape `(nwalkers, nsave)`.
"""
function run_mcmc(ensemble, lnprob, lnprobfn, steps; thin=1)
    nsave = div(steps, thin)

    chain = zeros((size(ensemble, 1), size(ensemble, 2), nsave))
    chainlnprob = zeros((size(ensemble, 2), nsave))
    for i in 1:steps
        ensemble, lnprob = update(ensemble, lnprob, lnprobfn)

        if i % thin == 0
            chain[:,:,div(i, thin)] = ensemble
            chainlnprob[:,div(i,thin)] = lnprob
        end
    end

    chain, chainlnprob
end

"""
    basic_callback(ensemble, lnprob, n, thin)

Print basic summary of the `run_to_neff` step just completed as a
callback.

Callbacks to `run_to_neff` will be passed the ensemble, the
corresponding lnprobs, the number of MCMC steps taken and the thin
parameter.

"""
function basic_callback(pts, lnprobs, n, thin)
    println(@sprintf("Advanced for %d steps (thin = %d); (max, mean, var)(log(pi)) = (%.1f, %.1f, %.1f)", n, thin, maximum(lnprobs), mean(lnprobs), var(lnprobs)))
    println("    ACLs are [" * join([@sprintf("%.1f", ac) for ac in Acor.acl(pts)], ",") * "]")
    flush(STDOUT)
end

"""
    run_to_neff(ensemble, lnprob, lnprobfn, neff; callback=nothing)

Return `(chain, chainlnprobs)` that results from running the MCMC
until accumulating at least `neff` independent ensembles worth of
samples.

Runs the MCMC in successively-doubled number of steps (with
corresponding thinning) until the desired number of effective
ensembles worth of samples are accumulated.  Resets the positions of
the walkers and the doubling steps every time the highest posterior
value increases by a significant amount (based on the dimension of the
sampling problem).

`chain` will have shape `(ndim, nwalkers, nensembles)` and
`chainlnprobs` will have shape `(nwalkers, nensembles)`.

"""
function run_to_neff(ensemble, lnprob, lnprobfn, neff; callback=nothing)
    n = 8*neff
    n0 = n
    nd = size(ensemble, 1)
    nw = size(ensemble, 2)
    thin = 1

    ps = reshape(ensemble, (nd, nw, 1))
    lnps = reshape(lnprob, (nw, 1))

    lnpmax = maximum(lnps)

    while true
        ps, lnps = run_mcmc(ps[:,:,end], lnps[:,end], lnprobfn, n, thin=thin)

        if !(callback == nothing)
            callback(ps, lnps, n, thin)
        end

        acls = Acor.acl(ps)
        amax = maximum(acls)
        nee = size(ps, 3)/amax

        if nee > neff
            break
        end

        # Now check whether we need to re-centre around the maximum...
        lpm = maximum(lnps)

        if lpm > lnpmax + 2.0*sqrt(nd/2.0) # If we have increased max log(L) by 2-sigma
            ne = size(ps, 3)

            # Count number of samples we have above mean - 3*sigma
            lnpthresh = lpm - nd/2.0 - 3.0*sqrt(nd/2.0)

            pslnps_above = Set()
            for i in 1:nw
                for j in 1:ne
                    p = ps[:,i,j]
                    l = lnps[i,j]

                    if l > lnpthresh
                        push!(pslnps_above, (p,l))
                    end
                end
            end
            na = length(pslnps_above)

            if na > 2*nd
                # Then we have at least 2*nd samples above threshold,
                # and we can re-start the sampling from general
                # position
                lnpmax = lpm
                n = n0
                thin = 1

                ps_above = zeros(nd, na)
                lnps_above = zeros(na)
                for (i, (p, lp)) in enumerate(pslnps_above)
                    ps_above[:,i] = p
                    lnps_above[i] = lp
                end

                new_ps = zeros(nd, nw, 1)
                new_lnps = zeros(nw, 1)

                for i in 1:nw
                    p = nothing
                    l = -Inf
                    while l < lnpthresh
                        j = rand(1:na)
                        p = ps_above[:,j]
                        l = lnps_above[j]
                    end
                    new_ps[:,i,1] = p
                    new_lnps[i,1] = l
                end

                ps = new_ps
                lnps = new_lnps
            end
        else
            n *= 2
            thin *= 2
        end
    end

    ps, lnps
end

end


text:
if VERSION > v"1.6"
    exit() # Don't want to build the system image, since we will use LBT
end

using PackageCompiler
using MKL_jll

# if no environment variable ENV["USE_BLAS64"] is set install.jl
# tries to change USE_BLAS64 = false
const USEBLAS64 = parse(Bool,get(ENV, "USE_BLAS64","false"))

include("../src/install.jl")
enable_mkl_startup()


text:
using ServerMetrics
using Sockets
using Dates
using ServerMetrics: AbstractMetric, NumericMetric, MetricGroup
using ServerMetrics: register!, unregister!, clear_registry!
using ServerMetrics: get_metric
using ServerMetrics: get_cells
using DataStructures: SortedDict
using Test: @test, @test_throws, @test_logs

get_cell_values(m) = Dict(cell.labels => cell.value[] for cell in get_cells(m))
metric_values(m::NumericMetric) = m.value[]
metric_values(m::MetricGroup) = Dict(k => v.value[] for (k, v) in m.cells)
metric_values(m::AbstractMetric) = metric_values(m.content)
is_dummy(x) = isa(x, ServerMetrics.DummyCell)

@testcase "get_cell! with invalid labels" begin
    cnt = Counter(;action=String, code=Int64)
    @test is_dummy(ServerMetrics.get_cell!(cnt; action="get"))
    @test is_dummy(ServerMetrics.get_cell!(cnt; action="get", code=200, unknown="foo"))
    @test is_dummy(ServerMetrics.get_cell!(cnt; action="get", code=nothing))
end

@testcase "get_cell_if_exists does not create new cells" begin
    cnt = Counter(;day=String)
    @test ServerMetrics.get_cell_if_exists(cnt; day="Monday") == nothing
    @test metric_values(cnt) == Dict()

    @test ServerMetrics.get_cell!(cnt; day="Monday") isa NumericMetric
    @test ServerMetrics.get_cell_if_exists(cnt; day="Monday") isa NumericMetric
    @test metric_values(cnt) ==Dict(Dict(:day => "Monday") => 0.0)
end

@testcase "Invalid access is ignored" begin
    cnt = Counter(;action=String)
    inc!(cnt; action="get")
    inc!(cnt; unknown=nothing) # This is invalid
    inc!(cnt; action="put", mystery=1) # This is invalid
    inc!(cnt) # This is also invalid (missing action label)
    @test metric_values(cnt) == Dict(
        Dict(:action => "get") => 1.0
    )
end

@testcase "get_cell! for singular counter" begin
    cnt = Counter()
    @test is_dummy(ServerMetrics.get_cell!(cnt; label=10))
    @test ServerMetrics.get_cell_if_exists(cnt; label=10) == nothing
    @test ServerMetrics.get_cell_if_exists(cnt) == cnt.content
    @test ServerMetrics.get_cell!(cnt) == cnt.content
end

@testcase "get_cell with Counter" begin
    cnt = Counter(;action=String, code=Int64)
    @test metric_values(cnt) == Dict()

    @test ServerMetrics.get_cell!(cnt; action="get", code=200) isa NumericMetric
    @test metric_values(cnt) == Dict(
        Dict(:action => "get", :code => 200) => 0.0,
    )

    ServerMetrics.get_cell!(cnt; action="get", code=500)
    ServerMetrics.get_cell!(cnt; action="put", code=200)
    @test metric_values(cnt) == Dict(
        Dict(:action => "get", :code => 200) => 0.0,
        Dict(:action => "get", :code => 500) => 0.0,
        Dict(:action => "put", :code => 200) => 0.0,
    )
end

@testcase "Counter with labels " begin
    c = Counter(;action=String, response_code=Int64)
    inc!(c; action="get", response_code=404)
    inc!(c; action="put", response_code=200)
    inc!(c, 2.0; action="get", response_code=404)

    @test metric_values(c) == Dict(
        Dict(:action => "get", :response_code => 404) => 3.0,
        Dict(:action => "put", :response_code => 200) => 1.0,
    )
end

@testcase "Gauge with labels" begin
    g = Gauge(;resource=String)
    set!(g, 20.0; resource="disk")
    inc!(g; resource="disk")
    dec!(g, 5.0; resource="water_level")
    inc!(g, 2.0; resource="water_level")
    @test metric_values(g) == Dict(
        Dict(:resource => "disk") => 21.0,
        Dict(:resource => "water_level") => -3.0,
    )
end

@testcase "Counter with no labels has one cell" begin
    c = Counter()
    @test length(get_cells(c)) == 1
    inc!(c, 2.5)
    @test metric_values(c) == 2.5
end

@testcase "Simple counter manipulations" begin
    c = ServerMetrics.Counter()
    @test metric_values(c) == 0.0
    inc!(c)
    @test metric_values(c) == 1.0
    inc!(c, 2)
    @test metric_values(c) == 3.0

    # Incrementing by negative value is logged and ignored
    inc!(c, -1.0)
    @test metric_values(c) == 3.0
end

@testcase "Gauge constructors" begin
    g1 = ServerMetrics.Gauge()
    @test metric_values(g1) == 0.0

    g2 = ServerMetrics.Gauge(10.0)
    @test metric_values(g2) == 10.0
end

@testcase "Simple gauge manipulations" begin
    g = ServerMetrics.Gauge()
    @test metric_values(g) == 0.0
    inc!(g)
    @test metric_values(g) == 1.0
    inc!(g, 2)
    @test metric_values(g) == 3.0

    inc!(g, -1.0)
    @test metric_values(g) == 3.0

    dec!(g, -1.0)
    @test metric_values(g) == 3.0
    # TODO(janrous): Perhaps we might consider chg!(g, v) which will
    # inc! or dec! based on the sign?
end

@testcase "Registry enforces unique names" begin
    r = ServerMetrics.MetricRegistry()
    c1 = ServerMetrics.Counter()
    register!(r, c1, "first_name")
    c2 = ServerMetrics.Counter()
    register!(r, c2, "second_name")
    @test ServerMetrics.name(c1) == "first_name"
    @test ServerMetrics.name(c2) == "second_name"
    c3 = ServerMetrics.Counter()
    @test_throws KeyError register!(r, c3, "first_name")
    @test ServerMetrics.name(c3) == nothing
end

@testcase "Once registered, name remains fixed" begin
    r = ServerMetrics.MetricRegistry()
    c = ServerMetrics.Counter()
    @test ServerMetrics.name(c) === nothing
    register!(r, c, "my_name")
    @test ServerMetrics.name(c) == "my_name"
    unregister!(r, "my_name")
    @test ServerMetrics.name(c) == "my_name"
end

@testcase "overwriting existing metrics" begin
    r = ServerMetrics.MetricRegistry()
    c1 = ServerMetrics.Counter()
    c2 = ServerMetrics.Counter()
    register!(r, c1, "first_name")
    @test_throws KeyError register!(r, c2, "first_name")
    @test get_metric(r, "first_name") == c1
    @test_logs (:warn,) register!(r, c2, "first_name"; overwrite=true)
    @test get_metric(r, "first_name") == c2
end

# Returns sorted Array of metric names registered within given registry.
list_metrics(r::ServerMetrics.MetricRegistry) = sort(collect(keys(r.metrics)))

# Uniqueness of metric name is enforced within a single registry. Two registries
# can have metric with the same name.
@testcase "Metrics unique within registry" begin
    r1 = ServerMetrics.MetricRegistry()
    r2 = ServerMetrics.MetricRegistry()
    g1 = ServerMetrics.Gauge(1.0)
    g2 = ServerMetrics.Gauge(3.0)
    register!(r1, g1, "unique")
    register!(r2, g2, "unique")
    @test list_metrics(r1) == ["unique"]
    @test list_metrics(r2) == ["unique"]
    @test metric_values(get_metric(r1, "unique")) == 1.0
    @test metric_values(get_metric(r2, "unique")) == 3.0
end

# Single metric can only be registered once.
@testcase "Name constraints on metrics are enforced" begin
    r1 = ServerMetrics.MetricRegistry()
    r2 = ServerMetrics.MetricRegistry()
    c = ServerMetrics.Counter()
    register!(r1, c, "my_counter")

    # Registering under the same name twice is okay
    @test_logs register!(r2, c, "my_counter")

    # Registering under different name is not okay
    @test_throws AssertionError register!(r2, c, "my_counter_2")
end

@testcase "Multiple registration with same name okay" begin
    r1 = ServerMetrics.MetricRegistry()
    r2 = ServerMetrics.MetricRegistry()
    c = ServerMetrics.Counter()
    register!(r1, c, "my_counter")
    @test_logs register!(r2, c, "my_counter")
end

# Metric could be registered sequentially in two registries if it's unregistered from one.
@testcase "Sequential multiple registration" begin
    r1 = ServerMetrics.MetricRegistry()
    r2 = ServerMetrics.MetricRegistry()
    c = ServerMetrics.Counter()
    register!(r1, c, "my_counter")
    unregister!(r1, "my_counter")
    register!(r2, c, "my_counter")
    @test list_metrics(r1) == []
    @test list_metrics(r2) == ["my_counter"]
end

@testcase "Metric name length enforced" begin
    @test_throws ArgumentError register!(
        MetricRegistry(),
        ServerMetrics.Counter(),
        "c"^201
    )
end

@testcase "Metric name validation" begin
    r = MetricRegistry()
    c = ServerMetrics.Counter()
    @test_throws ArgumentError register!(r, c, "hyphens-not-permitted")
    @test_throws ArgumentError register!(r, c, ".tralala")
    @test_throws ArgumentError register!(r, c, "2tralala")
    @test_throws ArgumentError register!(r, c, "tra\\lala")
    @test_throws ArgumentError register!(r, c, "t,r()a22ala")
    @test_throws ArgumentError register!(r, c, "prometheus.dislikes.this")
end

# Turn Dict keys into Return sorted array of keys from Dict
key_set(d::SortedDict) = Set([first(kv) for kv in d])

Base.@kwdef struct ThreeMetrics <: AbstractMetricCollection
    a = Counter()
    b = Counter()
    c = Gauge(5.0)
end

Base.@kwdef struct CounterAndTwoGauges <: AbstractMetricCollection
    c = Counter()
    g1 = Gauge()
    g2 = Gauge()
end

Base.@kwdef struct CounterAndGauge <: AbstractMetricCollection
    c = Counter()
    g = Gauge()
end

@testcase "Int64 counter manipulations" begin
    c = Counter()
    inc!(c, 2)
    @test metric_values(c) == 2.0
end

@testcase "Int64 gauge manipulations" begin
    g = Gauge(0)

    inc!(g, 1)
    @test metric_values(g) == 1.0

    dec!(g, 2)
    @test metric_values(g) == -1.0

    set!(g, 5)
    @test metric_values(g) == 5.0
end

@testcase "Unregister from registry" begin
    r = MetricRegistry()
    @test list_metrics(r) == []
    register!(r, Counter(), "aaa")
    register!(r, Counter(), "bbb")
    @test list_metrics(r) == ["aaa", "bbb"]
    @test_throws KeyError unregister!(r, "nonexistent")
    unregister!(r, "aaa")
    @test list_metrics(r) == ["bbb"]
    unregister!(r, "bbb")
    @test list_metrics(r) == []
end

@testcase "Clear registry" begin
    r = MetricRegistry()
    @test list_metrics(r) == []
    clear_registry!(r)
    @test list_metrics(r) == []
    register!(r, Counter(), "aaa")
    register!(r, Counter(), "bbb")
    @test list_metrics(r) == ["aaa", "bbb"]
    clear_registry!(r)
    @test list_metrics(r) == []
    register!(r, Counter(), "bbb")
    @test list_metrics(r) == ["bbb"]
end

Base.@kwdef struct MixedMetrics <: AbstractMetricCollection
    my_counter = Counter()
    my_gauge = Gauge()

    other_stuff::String = "default_value"
    anything::Any = nothing
end

@testcase "Collection with nonmetric fields" begin
    r = MetricRegistry(MixedMetrics())
    @test list_metrics(r) == ["my_counter", "my_gauge"]
end

Base.@kwdef struct ArbitraryStruct
    my_counter_1 = Counter()
    name::String = "String thingy"
end

@testcase "Collection of any type" begin
    r = MetricRegistry(ArbitraryStruct())
    @test list_metrics(r) == ["my_counter_1"]
end

Base.@kwdef struct FirstGroup <: AbstractMetricCollection
    first_counter = Counter()
    first_gauge = Gauge()
end

Base.@kwdef struct SecondGroup <: AbstractMetricCollection
    second_counter = Counter()
    second_gauge = Gauge()
end

@testcase "Register non-overlapping collections" begin
    first_group = FirstGroup()
    second_group = SecondGroup()
    r = MetricRegistry()

    @test list_metrics(r) == []

    register_collection!(r, first_group)
    @test list_metrics(r) == ["first_counter", "first_gauge"]

    register_collection!(r, second_group)
    @test list_metrics(r) == ["first_counter", "first_gauge", "second_counter", "second_gauge"]
end

@testcase "MAX_CELLS_PER_METRIC triggers" begin
    cnt = Counter(;order=Int64)
    for i in 1:205
        inc!(cnt, i; order=i)
    end
    # only the last 200 cells should be present. Timestmaps may not change fast enough
    # to ensure that LRU will behave in a stable manner.
    @test length(metric_values(cnt)) == 200
end

@testcase "Zero all metrics" begin
    r = MetricRegistry()
    c = register!(r, ServerMetrics.Counter(), "my_counter")
    g = register!(r, ServerMetrics.Gauge(), "my_gauge")
    c_free = ServerMetrics.Counter()

    inc!(c, 2.0)
    set!(g, 5.0)
    inc!(c_free)

    @test metric_values(c) == 2.0
    @test metric_values(g) == 5.0
    @test metric_values(c_free) == 1.0

    ServerMetrics.zero_all_metrics(r)

    @test metric_values(c) == 0.0
    @test metric_values(g) == 0.0
    @test metric_values(c_free) == 1.0
end

@testcase "value_of for simple metrics" begin
    r = MetricRegistry()
    c = register!(r, ServerMetrics.Counter(), "my_counter")
    g = register!(r, ServerMetrics.Gauge(), "my_gauge")
    inc!(c)
    set!(g, 50)
    @test ServerMetrics.value_of(r, "my_counter") == 1.0
    @test ServerMetrics.value_of(r, "my_gauge") == 50.0
    @test ServerMetrics.value_of(r, "unknown") == nothing
end

@testcase "value_of for labelled metrics" begin
    r = MetricRegistry()
    c = register!(r, ServerMetrics.Counter(;response_code=Int64, action=String), "my_counter")
    inc!(c; action="GET", response_code=200)
    inc!(c; action="PUT")  # This inc! is invalid due to bad labels and should be silently ignored.

    # The following doesn't work because all labels are missing.
    @test ServerMetrics.value_of(r, "my_counter") == nothing

    # The following doesn't work because unknown "badlabel" is present.
    @test ServerMetrics.value_of(r, "my_counter"; action="GET", response_code=200, badlabel=1) == nothing

    # The following doesn't work because "response_code" label is missing.
    @test ServerMetrics.value_of(r, "my_counter"; action="PUT") == nothing

    # The following works and refers to the cell created above.
    @test ServerMetrics.value_of(r, "my_counter"; action="GET", response_code=200) == 1.0
end


text:
using UnicodePlots

# a simple example for DNN_CTCLoss_With_Softmax usage
# by the way, it's for a single input case
S = 5
T = 100
E = 1000
x = 0.5*rand(S,T);
P = Variable(x,keepsgrad=true);
L = zeros(1, E)
for i = 1:E
    L[i] = DNN_CTC_With_Softmax(P, [2 3 4 5 2 3 4 5])
    backward()
    update!(P, 1E-3)
end

lineplot(vec(L))

# a simple example for DNN_Batch_CTCLoss_With_Softmax usage
# by the way, it's for a batch of inputs case
E = 100
S = 5
T1 = 10
T2 = 20
x1 = 0.5*rand(S,T1);  # first input
x2 = 0.3*rand(S,T2);  # second input
seq1 = [2 3 4];       # first input's label
seq2 = [3 2];         # second input's label

x   = [x1 x2];        # a batch of inputs
seq = [seq1 seq2];    # a batch of labels

P = Variable(x, keepsgrad=true);
L = zeros(1, E)
for i = 1:E
    L[i] = DNN_Batch_CTC_With_Softmax(P, seq, [T1 T2], [length(seq1) length(seq2)])
    backward()
    update!(P, 1E-3)
end

lineplot(vec(L))

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta