text:
"""
    turing_model(formula, data; model=Normal, priors=DefaultPrior(), standardize=false)

Create a Turing model using `formula` syntax and a `data` source.

# `formula`

`formula` is the the same friendly interface to specify used to specify statistical models by
[`brms`](https://paul-buerkner.github.io/brms/),
[`rstarnarm`](https://mc-stan.org/rstanarm/index.html),
[`bambi`](https://bambinos.github.io/bambi),
[`StatsModels.jl`](https://juliastats.org/StatsModels.jl/latest/) and
[`MixedModels.jl`](https://juliastats.org/MixedModels.jl/dev/).
The syntax is done by using the `@formula` macro and then specifying the dependent variable
followed by a tilde `~` then the independent variables separated by a plus sign `+`.

Example: `@formula(y ~ x1 + x2 + x3)`.

Moderations/interactions can be specified with the asterisk sign `*`, e.g. `x1 * x2`.
This will be expanded to `x1 + x2 + x1:x2`, which, following the principle of hierarchy,
the main effects must also be added along with the interaction effects. Here `x1:x2`
means that the values of `x1` will be multiplied (interacted) with the values of `x2`.

Random-effects (a.k.a. group-level effects) can be specified with the `(term | group)` inside
the `@formula`, where `term` is the independent variable and `group` is the **categorical**
representation (i.e., either a column of `String`s or a `CategoricalArray` in `data`).
You can specify a random-intercept with `(1 | group)`.

Example: `@formula(y ~ (1 | group) + x1)`.

**Notice: random-effects are currently only implemented for a single group-level intercept.
Future versions of `TuringGLM.jl` will support slope random-effects and multiple group-level
effets.**

# `data`

`data` can be any `Tables.jl`-compatible data interface.
The most popular ones are `DataFrame`s and `NamedTuple`s.

# `model`

`model` represents the likelihood function which you want to condition your data on.
It has to be a subtype of `Distributions.UnivariateDistribution`.
Currently, `TuringGLM.jl` supports:

* `Normal` (the default if not specified): linear regression
* `TDist`: robust linear regression
* `Bernoulli`: logistic regression
* `Poisson`: Poisson count data regression
* `NegativeBinomial`: negative binomial robust count data regression

# `priors`

`TuringGLM.jl` comes with state-of-the-art default priors, based on the literature and the
Stan community.
By default, `turing_model` will use `DefaultPrior`. But you can specify your own with
`priors=CustomPrior(predictors, intercept, auxiliary)`. All models take a `predictors` and
`intercept` priors.

In robust models, e.g. Linear Regression with Student-t likelihood or Count Regression
with Negative Binomial likelihood, often there is an extra auxiliary parameter that is
needed to parametrize to model to overcome under- or over-dispersion. If you are specifying
a custom prior for one of these type of models, then you should also specify a prior for
the auxiliary parameter.

Non-robust models do not need an auxiliary parameter and you can pass `nothing` as the
auxiliary argument.

Example for a non-robust model: `@formula(y, ...), data; priors=CustomPrior(Normal(0, 2.5), Normal(10, 5), nothing)`

Example for a robust model: `@formula(y, ...), data; priors=CustomPrior(Normal(0, 2.5), Normal(10, 5), Exponential(1))`

# `standardize`

Whether `true` or `false` to standardize your data to mean 0 and standard deviation 1
before inference. Some science fields prefer to analyze and report effects in terms of
standard devations. Also, whenever measurement scales differs, it is often suggested to
standardize the effects for better comparison. By default, `turing_model` sets `standardize=false`.
"""
function turing_model(
    formula::FormulaTerm,
    data;
    model::Type{T}=Normal,
    priors::Prior=DefaultPrior(),
    standardize::Bool=false,
) where {T<:UnivariateDistribution}
    return _turing_model(formula, data, T; priors, standardize)
end

function _turing_model(
    formula::FormulaTerm,
    data,
    ::Type{T};
    priors::Prior=DefaultPrior(),
    standardize::Bool=false,
) where {T<:UnivariateDistribution}
    # extract y, X and Z
    y = data_response(formula, data)
    X = data_fixed_effects(formula, data)
    Z = data_random_effects(formula, data)

    # μ and σ identities
    μ_X = 0
    σ_X = 1
    μ_y = 0
    σ_y = 1

    if standardize
        μ_X, σ_X, X = standardize_predictors(X)
        μ_y, σ_y, y = standardize_predictors(y)
        if !isnothing(Z)
            #TODO: implement random-effects slope
            throw(
                ArgumentError(
                    "TuringGLM currently does not support random-effects for slope terms"
                ),
            )
        end
    end

    # Random-Effects Conditionals
    prior = _prior(priors, y, T)
    ranef = TuringGLM.ranef(formula)
    model = if ranef === nothing
        _model(μ_X, σ_X, prior, T)
    else
        intercept_ranef = intercept_per_ranef(ranef)
        group_var = first(ranef).rhs
        idx = get_idx(term(group_var), data)
        # print for the user the idx
        println("The idx are $(last(idx))\n")
        _model(μ_X, σ_X, prior, intercept_ranef, idx, T)
    end
    return model(y, X)
end

# Default priors
_prior(prior::Prior, y, ::Type{<:UnivariateDistribution}) = prior
function _prior(::DefaultPrior, y, ::Type{Normal})
    return CustomPrior(TDist(3), median(y) + mad(y) * TDist(3), nothing)
end
function _prior(::DefaultPrior, y, ::Type{TDist})
    return CustomPrior(TDist(3), median(y) + mad(y) * TDist(3), Gamma(2, 0.1))
end
function _prior(::DefaultPrior, y, ::Type{Bernoulli})
    return CustomPrior(TDist(3), 2.5 * TDist(3), nothing)
end
function _prior(::DefaultPrior, y, ::Type{Poisson})
    return CustomPrior(TDist(3), 2.5 * TDist(3), nothing)
end
function _prior(::DefaultPrior, y, ::Type{NegativeBinomial})
    return CustomPrior(TDist(3), 2.5 * TDist(3), Gamma(0.01, 0.01))
end
function _prior(::DefaultPrior, y, T::Type{<:UnivariateDistribution})
    return throw(
        ArgumentError(
            "No default prior implemented for likelihood type $T. Please check the documentation for supported likelihoods.",
        ),
    )
end

# Models with Normal likelihood
function _model(μ_X, σ_X, prior, intercept_ranef, idx, ::Type{Normal})
    idxs = first(idx)
    n_gr = length(unique(first(idx)))
    @model function normal_model_ranef(
        y,
        X;
        predictors=size(X, 2),
        idxs=idxs,
        n_gr=n_gr,
        intercept_ranef=intercept_ranef,
        μ_X=μ_X,
        σ_X=σ_X,
        prior=prior,
        residual=1 / std(y),
    )
        α ~ prior.intercept
        β ~ filldist(prior.predictors, predictors)
        σ ~ Exponential(residual)
        μ = α .+ X * β
        if !isempty(intercept_ranef)
            τ ~ mad(y) * truncated(TDist(3), 0, Inf)
            zⱼ ~ filldist(Normal(), n_gr)
            αⱼ = zⱼ .* τ
            μ .+= αⱼ[idxs]
        end
        #TODO: implement random-effects slope
        y ~ MvNormal(μ, σ^2 * I)
        return (; α, β, σ, τ, zⱼ, αⱼ, y)
    end
end
function _model(μ_X, σ_X, prior, ::Type{Normal})
    @model function normal_model(
        y, X; predictors=size(X, 2), μ_X=μ_X, σ_X=σ_X, prior=prior, residual=1 / std(y)
    )
        α ~ prior.intercept
        β ~ filldist(prior.predictors, predictors)
        σ ~ Exponential(residual)
        y ~ MvNormal(α .+ X * β, σ^2 * I)
        return (; α, β, σ, y)
    end
end

# Models with Student-t likelihood
function _model(μ_X, σ_X, prior, intercept_ranef, idx, ::Type{TDist})
    @model function student_model_ranef(
        y,
        X;
        predictors=size(X, 2),
        idxs=first(idx),
        n_gr=length(unique(first(idx))),
        intercept_ranef=intercept_ranef,
        μ_X=μ_X,
        σ_X=σ_X,
        prior=prior,
        residual=1 / std(y),
    )
        α ~ prior.intercept
        β ~ filldist(prior.predictors, predictors)
        σ ~ Exponential(residual)
        ν ~ prior.auxiliary
        μ = α .+ X * β
        if !isempty(intercept_ranef)
            τ ~ 0 + mad(y) * truncated(TDist(3), 0, Inf)
            zⱼ ~ filldist(Normal(), n_gr)
            αⱼ = zⱼ .* τ
            μ .+= αⱼ[idxs]
        end
        #TODO: implement random-effects slope
        y ~ arraydist(μ + σ * TDist.(ν))
        return (; α, β, σ, ν, τ, zⱼ, αⱼ, y)
    end
end
function _model(μ_X, σ_X, prior, ::Type{TDist})
    @model function student_model(
        y, X; predictors=size(X, 2), μ_X=μ_X, σ_X=σ_X, prior=prior, residual=1 / std(y)
    )
        α ~ prior.intercept
        β ~ filldist(prior.predictors, predictors)
        σ ~ Exponential(residual)
        ν ~ prior.auxiliary
        y ~ arraydist((α .+ X * β) .+ σ .* TDist.(ν))
        return (; α, β, σ, ν, y)
    end
end

# Models with Bernoulli likelihood
function _model(μ_X, σ_X, prior, intercept_ranef, idx, ::Type{Bernoulli})
    @model function bernoulli_model_ranef(
        y,
        X;
        predictors=size(X, 2),
        idxs=first(idx),
        n_gr=length(unique(first(idx))),
        intercept_ranef=intercept_ranef,
        μ_X=μ_X,
        σ_X=σ_X,
        prior=prior,
    )
        α ~ prior.intercept
        β ~ filldist(prior.predictors, predictors)
        μ = α .+ X * β
        if !isempty(intercept_ranef)
            τ ~ mad(y) * truncated(TDist(3), 0, Inf)
            zⱼ ~ filldist(Normal(), n_gr)
            αⱼ = zⱼ .* τ
            μ .+= αⱼ[idxs]
        end
        #TODO: implement random-effects slope
        y ~ arraydist(LazyArray(@~ BernoulliLogit.(μ)))
        return (; α, β, τ, zⱼ, αⱼ, y)
    end
end
function _model(μ_X, σ_X, prior, ::Type{Bernoulli})
    @model function bernoulli_model(
        y, X; predictors=size(X, 2), μ_X=μ_X, σ_X=σ_X, prior=prior
    )
        α ~ prior.intercept
        β ~ filldist(prior.predictors, predictors)
        y ~ arraydist(LazyArray(@~ BernoulliLogit.(α .+ X * β)))
        return (; α, β, y)
    end
end

# Models with Poisson likelihood
function _model(μ_X, σ_X, prior, intercept_ranef, idx, ::Type{Poisson})
    @model function poisson_model_ranef(
        y,
        X;
        predictors=size(X, 2),
        idxs=first(idx),
        n_gr=length(unique(first(idx))),
        intercept_ranef=intercept_ranef,
        μ_X=μ_X,
        σ_X=σ_X,
        prior=prior,
    )
        α ~ prior.intercept
        β ~ filldist(prior.predictors, predictors)
        μ = α .+ X * β
        if !isempty(intercept_ranef)
            τ ~ mad(y) * truncated(TDist(3), 0, Inf)
            zⱼ ~ filldist(Normal(), n_gr)
            αⱼ = zⱼ .* τ
            μ .+= αⱼ[idxs]
        end
        #TODO: implement random-effects slope
        y ~ arraydist(LazyArray(@~ LogPoisson.(μ)))
        return (; α, β, τ, zⱼ, αⱼ, y)
    end
end
function _model(μ_X, σ_X, prior, ::Type{Poisson})
    @model function poisson_model(
        y, X; predictors=size(X, 2), μ_X=μ_X, σ_X=σ_X, prior=prior
    )
        α ~ prior.intercept
        β ~ filldist(prior.predictors, predictors)
        y ~ arraydist(LazyArray(@~ LogPoisson.(α .+ X * β)))
        return (; α, β, y)
    end
end

# Models with NegativeBinomial likelihood
function _model(μ_X, σ_X, prior, intercept_ranef, idx, ::Type{NegativeBinomial})
    @model function negbin_model_ranef(
        y,
        X;
        predictors=size(X, 2),
        idxs=first(idx),
        n_gr=length(unique(first(idx))),
        intercept_ranef=intercept_ranef,
        μ_X=μ_X,
        σ_X=σ_X,
        prior=prior,
    )
        α ~ prior.intercept
        β ~ filldist(prior.predictors, predictors)
        ϕ⁻ ~ prior.auxiliary
        ϕ = 1 / ϕ⁻
        μ = α .+ X * β
        if !isempty(intercept_ranef)
            τ ~ mad(y) * truncated(TDist(3), 0, Inf)
            zⱼ ~ filldist(Normal(), n_gr)
            αⱼ = zⱼ .* τ
            μ .+= αⱼ[idxs]
        end
        #TODO: implement random-effects slope
        y ~ arraydist(LazyArray(@~ NegativeBinomial2.(exp.(μ), ϕ)))
        return (; α, β, ϕ, τ, zⱼ, αⱼ, y)
    end
end
function _model(μ_X, σ_X, prior, ::Type{NegativeBinomial})
    @model function negbin_model(y, X; predictors=size(X, 2), μ_X=μ_X, σ_X=σ_X, prior=prior)
        α ~ prior.intercept
        β ~ filldist(prior.predictors, predictors)
        ϕ⁻ ~ prior.auxiliary
        ϕ = 1 / ϕ⁻
        y ~ arraydist(LazyArray(@~ NegativeBinomial2.(exp.(α .+ X * β), ϕ)))
        return (; α, β, ϕ, y)
    end
end

# Unsupported models
function _model(μ_X, σ_X, prior, intercept_ranef, idx, T::Type{<:UnivariateDistribution})
    return throw(
        ArgumentError(
            "No Turing model implemented for likelihood type $T. Please check the documentation for supported likelihoods.",
        ),
    )
end
function _model(μ_X, σ_X, prior, T::Type{<:UnivariateDistribution})
    return throw(
        ArgumentError(
            "No Turing model implemented for likelihood type $T. Please check the documentation for supported likelihoods.",
        ),
    )
end

"""
    NegativeBinomial2(μ, ϕ)

An alternative parameterization of the Negative Binomial distribution:

```math
\\text{Negative-Binomial}(n \\mid \\mu, \\phi) \\sim \\binom{n + \\phi - 1}{n} \\left( \\frac{\\mu}{\\mu + \\phi} \\right)^{n!} \\left( \\frac{\\phi}{\\mu + \\phi} \\right)^{\\phi!}
```

where the expectation is μ and variance is (μ + μ²/ϕ).

The alternative parameterization is inspired by the [Stan's `neg_binomial_2` function](https://mc-stan.org/docs/functions-reference/nbalt.html).
"""
function NegativeBinomial2(μ::T, ϕ::T) where {T<:Real}
    p = max(1 / (1 + μ / ϕ), 1e-6) # numerical stability
    r = ϕ
    return NegativeBinomial(r, p)
end


text:
using BenchmarkTools
using FFTW
using LinearAlgebra
using MPI
using PencilFFTs
using Random

function fft_inplace!(a, b, fft_plan)
    mul!(b, fft_plan, a)
end

function runtests()
    suite = BenchmarkGroup()

    if ~MPI.Initialized()
        MPI.Init()
    end

    n = 128
    dims = (n, n, n)

    suite["PencilFFTs"] = BenchmarkGroup()
    for (flag_name, flags) in [
                           ("FFTW.ESTIMATE", FFTW.ESTIMATE),
                           ("FFTW.MEASURE", FFTW.MEASURE),
                           ("FFTW.PATIENT", FFTW.PATIENT),
                          ]

        # Distribute 12 processes on a 1 × 1 grid.
        comm = MPI.COMM_WORLD
        Nproc = MPI.Comm_size(comm)

        # Let MPI_Dims_create choose the decomposition.
        proc_dims = let pdims = zeros(Int, 2)
            MPI.Dims_create!(Nproc, pdims)
            pdims[1], pdims[2]
        end

        transform = Transforms.FFT()

        fft_plan = PencilFFTPlan(dims, transform, proc_dims, comm; fftw_flags=flags)

        # Allocate and initialise input data, and apply transform.
        a_pen = allocate_input(fft_plan)
        b_pen = allocate_output(fft_plan)
        randn!(a_pen)
        randn!(b_pen)

        suite["PencilFFTs"][flag_name] = @benchmarkable fft_inplace!($a_pen, $b_pen, $fft_plan)
    end

    suite["FFTW"] = BenchmarkGroup()
    for (flag_name, flags) in [
                           ("FFTW.ESTIMATE", FFTW.ESTIMATE),
                           ("FFTW.MEASURE", FFTW.MEASURE),
                           ("FFTW.PATIENT", FFTW.PATIENT),
                          ]
        a_arr = Array{ComplexF64}(undef, dims)
        b_arr = Array{ComplexF64}(undef, dims)
        rand!(a_arr)
        rand!(b_arr)

        fft_plan = plan_fft(a_arr; flags=flags)

        suite["FFTW"][flag_name] = @benchmarkable fft_inplace!($a_arr, $b_arr, $fft_plan)
    end

    @show tune!(suite)

    results = run(suite)

    @show res_min = minimum(results)
    @show res_med = median(results)

    @show judge(res_med["FFTW"], res_med["PencilFFTs"])
end

runtests()

ext: jl


text:
Pkg.add("ODBC")
Pkg.add("Dates")
using ODBC
ODBC.connect("pubsODBC")
	# ODBC Connection Object
	# ----------------------
	# Connection Data Source: pubsODBC
	# pubsODBC Connection Number: 1
	# Contains resultset? No
results = query("select * from titles")
# elapsed time: 0.189238158 seconds
# 18x10 DataFrame
# | Row | title_id |
# |-----|----------|
# | 1   | "BU1032" |
# | 2   | "BU1111" |
# | 3   | "BU2075" |
# | 4   | "BU7832" |
# | 5   | "MC2222" |
# | 6   | "MC3021" |
# ⋮
# | 12  | "PS2091" |
# | 13  | "PS2106" |
# | 14  | "PS3333" |
# | 15  | "PS7777" |
# | 16  | "TC3218" |
# | 17  | "TC4203" |
# | 18  | "TC7777" |

# | Row | title                                                             |
# |-----|-------------------------------------------------------------------|
# | 1   | "The Busy Executive's Database Guide"                             |
# | 2   | "Cooking with Computers: Surreptitious Balance Sheets"            |
# | 3   | "You Can Combat Computer Stress!"                                 |
# | 4   | "Straight Talk About Computers"                                   |
# | 5   | "Silicon Valley Gastronomic Treats"                               |
# | 6   | "The Gourmet Microwave"                                           |
# ⋮
#  Row | _type          | pub_id | price | advance | royalty | ytd_sales |
# -----|----------------|--------|-------|---------|---------|-----------|
#  1   | "business    " | "1389" | 19.99 | 5000.0  | 10      | 4095      |
#  2   | "business    " | "1389" | 11.95 | 5000.0  | 10      | 3876      |
#  3   | "business    " | "0736" | 2.99  | 10125.0 | 24      | 18722     |
#  4   | "business    " | "1389" | 19.99 | 5000.0  | 10      | 4095      |
#  5   | "mod_cook    " | "0877" | 19.99 | 0.0     | 12      | 2032      |
#  6   | "mod_cook    " | "0877" | 2.99  | 15000.0 | 24      | 22246     |
# fname = "books.dat"
# query("select * from titles", file=fname, delim = '\t')

# updates:
updsql = "update titles set type = 'psychology' where title_id='BU1032'"
query(updsql)

listdrivers()
listdsns()

max_stars_repo_path: Module 1/Chapter08/odbc.jl

max_issues_repo_path: Module 1/Chapter08/odbc.jl

max_forks_repo_path: Module 1/Chapter08/odbc.jl


text:
using SeisIO, JLD2, Distributed, DelimitedFiles, Sockets,Logging

export segy2jld2
#

# Converting PASSCAL segy data to JLD2 format. There are options to save to JLD2 format
# as single file or separate files for each segy file.
"""
Key features:
1. Can be used to convert individual SEGY (PASSCAL VARIATION FORMAT ONLY FOR NOW) file or a group/list of SEGY files.
2. Handles the situation when multiple data files were download for the same station
    and channel. Those files will be merged into one single SeisData object with
    multiple SeisChannel entries.
3. All SEGY headers are kept and saved in :misc filed of the SeisData object.

"""

"""
Naming convention for segyfiles when need to extract/update station information.
    YY.DDD.HH.MM.SS.SITE.CHAN
    02.089.00.01.50.7460.1
"""

"""

    segy2jld2(infile)

Saves individual segy file (in) to JLD2 file.
The output file name will be the same as the input, with new extension *.jld2.

"""

function segy2jld2(infile::String;siteinfodict::Dict=Dict())
    print("Same file, new extension of jld2.\n")
    print("Converting PASSCAL SEGY file: [ ",infile," ] ... \n")

    if occursin("segy",infile) || occursin("SEGY",infile)
        filebase=infile[1:end-4]
    else
        filebase=infile
    end
    outfile = join([filebase,".jld2"]);

    datin = SeisIO.read_data("passcal",infile, swap=true, full=true);
    if length(siteinfodict) > 0
        net,sta,loc,chan=split(datin.id[1],".")
        if occursin("/",filebase)
            filebasetemp=split(filebase,"/")
            parsename=split(filebasetemp[end],".")
        else
            parsename=split(filebase,".")
        end
        if length(net)==0 #no net information
            net=siteinfodict[collect(keys(siteinfodict))[1]][1]
            @warn "Using the first netcode $net in the siteinfodict as the netcode. This should be changed for a more generic usage."
        end
        if length(sta)==0 #no station information
            sta=parsename[6] #get the station name from the infile name, the 6th element
            @warn "Using the station name $sta from the 6th element of the filename."
        end
        if length(chan)==0 #no channel information
            chan=parsename[7] #get the channel name from the infile name, the 7th element
            @warn "Using the channel name $chan from the 7th element of the filename."
        end
        siteid=join([net,sta,loc,chan],".")
        update_siteinfo!(datin,siteid,siteinfodict)
    end
    datin.misc[1]["dlerror"] = 0

    stationlist = String[];
    push!(stationlist,datin.id[1])  #pass the 1st element of datin.id[] to the station list
    # stationinfo = Dict(["stationlist" => stationlist,"stationmethod" => "SEGY2jld2","stationsrc" => "converted"])
    stimetemp = datin.t[1][3]*1e-6; #convert mseconds to seconds;
    datlength = datin.t[1][2]/datin.fs; #convert from number of samples to time duration.
    starttime = u2d(stimetemp[1]); #convert union time to calendar time;
    endtime = u2d(stimetemp[1] + datlength[1]);
    DLtimestamplist = String[];
    push!(DLtimestamplist,replace(string(starttime),"-" => "."));

    #save to jld2
    stemp = joinpath(DLtimestamplist[1],stationlist[1]);
    file = jldopen(outfile, "w");
    file["info/DLtimestamplist"] = DLtimestamplist;
    file["info/stationlist"] = stationlist;
    file["info/starttime"]   = string(starttime)
    file["info/endtime"]     = string(endtime)
    file["info/DL_time_unit"]= datlength[1]; #length of the data in time (duration), not number of samples

    #save the waveform data as SeisData format.
    file[stemp] = SeisData(datin);

    JLD2.close(file);

    print("Saved JLD2 file to: ", outfile, "\n");
end

"""

    segy2jld2(infile,outfile)

Saves individual segy file (in) to JLD2 file (out), which is explicitly specified.

"""

function segy2jld2(infile::String,outfile::String;siteinfodict::Dict=Dict())
    print("Converting PASSCAL SEGY file: [ ",infile," ] ... \n")
    if occursin("segy",infile) || occursin("SEGY",infile)
        filebase=infile[1:end-4]
    else
        filebase=infile
    end

    datin = SeisIO.read_data("passcal",infile, swap=true, full=true);
    if length(siteinfodict) > 0
        net,sta,loc,chan=split(datin.id[1],".")
        if occursin("/",filebase)
            filebasetemp=split(filebase,"/")
            parsename=split(filebasetemp[end],".")
        else
            parsename=split(filebase,".")
        end
        if length(net)==0 #no net information
            net=siteinfodict[collect(keys(siteinfodict))[1]][1]
            @warn "Using the first netcode $net in the siteinfodict as the netcode. This should be changed for a more generic usage."
        end
        if length(sta)==0 #no station information
            sta=parsename[6] #get the station name from the infile name, the 6th element
            @warn "Using the station name $sta from the 6th element of the filename."
        end
        if length(chan)==0 #no channel information
            chan=parsename[7] #get the channel name from the infile name, the 7th element
            @warn "Using the channel name $chan from the 7th element of the filename."
        end
        siteid=join([net,sta,loc,chan],".")
        update_siteinfo!(datin,siteid,siteinfodict)
    end

    datin.misc[1]["dlerror"] = 0

    stationlist = String[];
    push!(stationlist,datin.id[1])  #pass the 1st element of datin.id[] to the station list
    # stationinfo = Dict(["stationlist" => stationlist,"stationmethod" => "SEGY2jld2","stationsrc" => "converted"])
    stimetemp = datin.t[1][3]*1e-6; #convert mseconds to seconds;
    datlength = datin.t[1][2]/datin.fs; #convert from number of samples to time duration.
    starttime = u2d(stimetemp[1]); #convert union time to calendar time;
    endtime = u2d(stimetemp[1] + datlength[1]);
    DLtimestamplist = String[];
    push!(DLtimestamplist,replace(string(starttime),"-" => "."));

    #save to jld2
    stemp = joinpath(DLtimestamplist[1],stationlist[1]);
    file = jldopen(outfile, "w");
    file["info/DLtimestamplist"] = DLtimestamplist;
    file["info/stationlist"] = stationlist;
    file["info/starttime"]   = string(starttime)
    file["info/endtime"]     = string(endtime)
    file["info/DL_time_unit"]= datlength[1]; #length of the data in time (duration), not number of samples

    #save the waveform data as SeisData format.
    file[stemp] = SeisData(datin);

    JLD2.close(file);

    print("Saved JLD2 file to: ", outfile, "\n");
end

"""

    segy2jld2(filelist,timestamp,outfile;verbose=false)

Saves segy files based on a filelist generated using ls() in Julia to
JLD2 file (out), which is explicitly specified.

    Arguments:
    filelist: An String Array including all segyfiles with the same tiemstamp, e.g., same day
    tiemstamp: this must be specified for the segy files, assuming they all share
                the same time stamp. Suggested format: 2019.06.27T00.00.00
    outfile: output file name for the JLD2 file. All data fiels will be readin and
                saved to a single JLD2 file.
    verbose: if true, the code will print out more messages. Default is false.
"""

function segy2jld2(filelist::Array{String,1},timestamp::String,outfile::String;verbose::Bool=false,siteinfodict::Dict=Dict())
    # print("I convert and pack all SEGY files to a JLD2 file.\n")

    #read the filenames from the infilelist into a strinng array
    # filelist = readfilelist(infilelist)

    stationlist = String[];
    DLtimestamplist = String[];
    push!(DLtimestamplist,timestamp); #use the specified timestamp for the group.

    #get some intial information from the first SEGY file.
    datin0 = SeisIO.read_data("passcal",filelist[1], swap=true, full=true);
    stimemin = datin0.t[1][3]*1e-6; #initial starttime in seconds;
    ftemp = datin0.t[1][2]/datin0.fs; #convert from number of samples to time duration.
    datlengthmax = ftemp[1]
    etimemax = stimemin + datlengthmax;

    file = jldopen(outfile, "w");
    file["info/DLtimestamplist"] = DLtimestamplist;
    stemp_pre = "-"
    count = 0
    # stationlisttemp = String[];
    #To deal with multiple channel issue, I am doing the loop twice.
    # 1. The first loop with collect multiple channels into one SeisData object.
    # 2. The second loop will work through the actual number of SeisData objects, counted by `count`
    Sall = SeisData[];
    Str_all = String[];

    #in the loop for each SEGY file. we will find the minimum starttime, maximum endtime, and maximum SEGYlength.
    for infile = filelist
        println("Converting SEGY file: [ ",infile," ]")

        if occursin("segy",infile) || occursin("SEGY",infile)
            filebase=infile[1:end-4]
        else
            filebase=infile
        end

        datin = SeisIO.read_data("passcal",infile, swap=true, full=true);
        if length(siteinfodict) > 0
            net,sta,loc,chan=split(datin.id[1],".")
            if occursin("/",filebase)
                filebasetemp=split(filebase,"/")
                parsename=split(filebasetemp[end],".")
            else
                parsename=split(filebase,".")
            end
            if length(net)==0 #no net information
                net=siteinfodict[collect(keys(siteinfodict))[1]][1]
                @warn "Using the first netcode $net in the siteinfodict as the netcode. This should be changed for a more generic usage."
            end
            if length(sta)==0 #no station information
                sta=parsename[6] #get the station name from the infile name, the 6th element
                @warn "Using the station name $sta from the 6th element of the filename."
            end
            if length(chan)==0 #no channel information
                chan=parsename[7] #get the channel name from the infile name, the 7th element
                @warn "Using the channel name $chan from the 7th element of the filename."
            end
            siteid=join([net,sta,loc,chan],".")
            update_siteinfo!(datin,siteid,siteinfodict)
        end

        datin.misc[1]["dlerror"] = 0
        push!(stationlist,datin.id[1])  #pass the 1st element of datin.id[] to the station list
        stimetemp = datin.t[1][3]*1e-6; #convert mseconds to seconds;
        if stimetemp < stimemin
            stimemin = stimetemp
        end
        SEGYlength = datin.t[1][2]/datin.fs; #convert from number of samples to time duration.
        if SEGYlength[1] > datlengthmax
            datlengthmax = SEGYlength[1]
        end
        etimetemp = stimetemp + SEGYlength[1]
        if etimetemp > etimemax
            etimemax = etimetemp
        end

        #save to jld2
        stemp = joinpath(timestamp,datin.id[1]);
        if stemp == stemp_pre
            println(stemp," exists. Appending to form multiple channels.")
            append!(Sall[count],SeisData(datin))
        else
            push!(Sall,SeisData(datin))
            push!(Str_all,stemp)
            count += 1
            # file[stemp] = SeisData(datin);
        end
        stemp_pre = stemp
    end
    starttime = u2d(stimemin); #convert union time to calendar time;
    endtime = u2d(etimemax);
    #the following info dat is saved after running through the whole group.
    file["info/stationlist"] = unique(stationlist);
    file["info/starttime"]   = string(starttime)
    file["info/endtime"]     = string(endtime)
    file["info/DL_time_unit"]= datlengthmax; #length of the data in time (duration), not number of samples

    #2. start to save the data to jld2
    for i = 1:count
        file[Str_all[i]] = Sall[i]
    end #end of loop for all SEGY files within one group/directory

    JLD2.close(file);

    print("Saved JLD2 file to: ", outfile, "\n");
end

"""

    segy2jld2(datdirlist,timestamplist,outfile;verbose=false)

Searches for all segy files scanning throgh the datdirlist (Array{String,1}). The segy files
    in each segy directory are assigned to one time stamp based on timestamplist (Array{String,1}).
    The two Arrays must be equal in size. The list can be generated from a file list (using ls in shell),
    using a utility function in this package: readfilelist(infilelist::String).

All segy files will be saved into one single JLD2 file (out), which is explicitly specified as outfile.

    Arguments:
    datdirlist: An array of strings storing the directories grouped by time stamp, e.g., one day.
    timestamplist: An array of strings storing the time stamps assigning to each SEGY file group.
    outfile: output file name for the JLD2 file. All data fiels will be readin and
                saved to a single JLD2 file.
    verbose: if true, the code will print out more messages. Default is false.
"""

function segy2jld2(datdirlist::Array{String,1},timestamplist::Array{String,1},outfile::String;verbose::Bool=false,siteinfodict::Dict=Dict())
    if length(datdirlist) != length(timestamplist)
        error("datdirlist and timestamplist must be the same length!")
    end

    stationlist = String[];

    # datlengthmax = 0; #initiate the length as zero, will find the maximum value after looping through the files.

    file = jldopen(outfile, "w");
    file["info/DLtimestamplist"] = timestamplist;
    file["info/starttime"] = timestamplist[1];
    #here we use the first value in timestamplist as the starttime for a group of multiple timestamp data.
    file["info/endtime"] = timestamplist[end]; #similar to starttime, here we use the last value in timestamplist as the endtime.

    #in the loop for each segy file in each SEGY directory. we will find the maximum datlength.
    # countdir=1;
    t1all=0
    t2all=0
    for (datdir,ts) in zip(datdirlist,timestamplist)

        t0=time()
        filelist=ls(joinpath(datdir,"*"));

        print("Converting for directory: [ ",datdir," ] ... \n")

        count = 0
        # stationlisttemp = String[];
        #To deal with multiple channel issue, I am doing the loop twice.
        # 1. The first loop with collect multiple channels into one SeisData object.
        # 2. The second loop will work through the actual number of SeisData objects, counted by `count`
        Sall = SeisData[];
        Str_all = String[];
        stemp_pre = "-"
        for infile = filelist
            datin = SeisIO.read_data("passcal",infile,swap=true,full=true);
            if length(siteinfodict) > 0
                net,sta,loc,chan=split(datin.id[1],".")
                if occursin("/",filebase)
                    filebasetemp=split(filebase,"/")
                    parsename=split(filebasetemp[end],".")
                else
                    parsename=split(filebase,".")
                end
                if length(net)==0 #no net information
                    net=siteinfodict[collect(keys(siteinfodict))[1]][1]
                    @warn "Using the first netcode $net in the siteinfodict as the netcode. This should be changed for a more generic usage."
                end
                if length(sta)==0 #no station information
                    sta=parsename[6] #get the station name from the infile name, the 6th element
                    @warn "Using the station name $sta from the 6th element of the filename."
                end
                if length(chan)==0 #no channel information
                    chan=parsename[7] #get the channel name from the infile name, the 7th element
                    @warn "Using the channel name $chan from the 7th element of the filename."
                end
                siteid=join([net,sta,loc,chan],".")
                update_siteinfo!(datin,siteid,siteinfodict)
            end

            datin.misc[1]["dlerror"] = 0

            if datdir == datdirlist[1] && infile == filelist[1]
                # println(datdir)
                file["info/DL_time_unit"]= datin.t[1][2]/datin.fs; #length of the data in time (duration), not number of samples;
                #convert from number of samples to time duration.
            end

            stemp = joinpath(ts,datin.id[1]);

            if stemp == stemp_pre
                if verbose == true
                    println(stemp," exists. Append to form multiple channels.")
                end
                append!(Sall[count],SeisData(datin))
            else
                push!(Sall,SeisData(datin))
                push!(Str_all,stemp)
                count += 1
                # file[stemp] = SeisData(datin);
            end

            stemp_pre = stemp;

        end #end of loop for all segy files within one group/directory

        #2. start to save the data to jld2
        for i = 1:count
            file[Str_all[i]] = Sall[i]
        end #end of loop for all segy files within one group/directory

        stationlist = unique(vcat(stationlist,keys(file[ts])))
        if verbose == true
            print("                ----------------> ",length(filelist)," segy files, time used: ",time() - t0,"\n")
        end

    end #end of loop for all SEGY directories.

    # println([t1all, t2all])
    #the following info is saved after running through the whole group.
    file["info/stationlist"] = stationlist;

    JLD2.close(file);

    print("Saved JLD2 file to: ", outfile, "\n");
end


text:
include("../model/SensorModel.jl")
include("../model/VelocityControlledDifferentialDrive.jl")
using Base.Test

let robot_pose = initial_state(0, 0, 0) # X
    marker_position = [1.5 0]
    ε = 0.001

    sensed_position = detect_marker(robot_pose, marker_position)

    @test sensed_position[1] ≈ 0 atol=ε
    @test sensed_position[2] ≈ 1.5 atol=ε
end

let robot_pose = initial_state(0, 0, π/2) # B
    marker_position = [1.5 0]
    ε = 0.001

    sensed_position = detect_marker(robot_pose, marker_position)

    @test sensed_position[1] ≈ 1.5 atol=ε
    @test sensed_position[2] ≈ 0 atol=ε
end

let robot_pose = initial_state(-1, 0, 0)
    marker_position = [1.5 0]
    ε = 0.001

    sensed_position = detect_marker(robot_pose, marker_position)

    @test sensed_position[1] ≈ 0 atol=ε
    @test sensed_position[2] ≈ 2.5 atol=ε
end

let robot_pose = initial_state(0, -1, π/2) # C
    marker_position = [1.5 0]
    ε = 0.001

    sensed_position = detect_marker(robot_pose, marker_position)

    @test sensed_position[1] ≈ 1.5 atol=ε
    @test sensed_position[2] ≈ 1 atol=ε
end

let robot_pose = initial_state(-1, -1, 0) # D
    marker_position = [1.5 0]
    ε = 0.001

    sensed_position = detect_marker(robot_pose, marker_position)

    @test sensed_position[1] ≈ -1 atol=ε
    @test sensed_position[2] ≈ 2.5 atol=ε
end

let robot_pose = initial_state(1, 2, π/2) # F
    marker_position = [1.5 0]
    ε = 0.001

    sensed_position = detect_marker(robot_pose, marker_position)

    @test sensed_position[1] ≈ 0.5 atol=ε
    @test sensed_position[2] ≈ -2 atol=ε
end

let robot_pose = initial_state(-1, -0.5, π/2)
    marker_position = [1.5 0.5]
    ε = 0.001

    sensed_position = detect_marker(robot_pose, marker_position)

    @test sensed_position[1] ≈ 2.5 atol=ε
    @test sensed_position[2] ≈ 1 atol=ε
end

let robot_pose = initial_state(2, 1, π/2 + π/4)
    marker_position = [1 2]
    ε = 0.001

    sensed_position = detect_marker(robot_pose, marker_position)

    @test sensed_position[1] ≈ 0 atol=ε
    @test sensed_position[2] ≈ 1.414214 atol=ε
end

let robot_pose = initial_state(-1, 1, π)
    marker_position = [1.5 1]
    ε = 0.001

    sensed_position = detect_marker(robot_pose, marker_position)

    @test sensed_position[1] ≈ 0 atol=ε
    @test sensed_position[2] ≈ -2.5 atol=ε
end

let robot_pose = initial_state(-1, 1, -π)
    marker_position = [1.5 1]
    ε = 0.001

    sensed_position = detect_marker(robot_pose, marker_position)

    @test sensed_position[1] ≈ 0 atol=ε
    @test sensed_position[2] ≈ -2.5 atol=ε
end

max_stars_repo_name: grgomrton/robot-model

max_issues_repo_name: grgomrton/robot-model

max_forks_repo_name: grgomrton/robot-model


text:
using StochasticGene
using Test

function fit_rna_test(root=".")
    gene = "CENPL"
    cell = "HCT116"
    fish = false
    nalleles = 2
    r = [0.01, 0.1, 1.0, 0.01006327034802035, 0.05]
    decayrate = 0.01006327034802035
    model = StochasticGene.model_rna(r,2,nalleles,1,.01,[1,2,3],(),decayrate,1.0,.05,fish)
    data = StochasticGene.data_rna(gene,"MOCK","data/HCT116_testdata",fish,"scRNA_test",root)
    options = StochasticGene.MHOptions(10000,2000,0,120.,1.,100.)
    fit,stats,waic = StochasticGene.run_mh(data,model,options,1);
    return stats.meanparam, fit.llml, model
end

function teststeadystatemodel(r,G,nhist,nalleles)
    g1 = StochasticGene.steady_state(r[1:2*G],G-1,nhist,nalleles)
    g2 = StochasticGene.simulatorGM(r[1:2*G],G-1,nhist,nalleles)
    return g1, g2
end

@testset "StochasticGene" begin

        p,ll,model =  fit_rna_test()
        @test isapprox(p,[0.0097; 0.10; 0.99],rtol=0.05)
        @test isapprox(ll,4450,rtol=0.05)

        r = [0.01, 0.1, 1.0, 0.01]
        h1,h2 = teststeadystatemodel(r,2,60,2)

        @test isapprox(h1,h2,rtol=0.05)

end


text:
# Autogenerated wrapper script for HarfBuzz_jll for x86_64-apple-darwin
export libharfbuzz, libharfbuzz_gobject, libharfbuzz_subset

using Cairo_jll
using Fontconfig_jll
using FreeType2_jll
using Glib_jll
using Graphite2_jll
using Libffi_jll
JLLWrappers.@generate_wrapper_header("HarfBuzz")
JLLWrappers.@declare_library_product(libharfbuzz, "@rpath/libharfbuzz.0.dylib")
JLLWrappers.@declare_library_product(libharfbuzz_gobject, "@rpath/libharfbuzz-gobject.0.dylib")
JLLWrappers.@declare_library_product(libharfbuzz_subset, "@rpath/libharfbuzz-subset.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Cairo_jll, Fontconfig_jll, FreeType2_jll, Glib_jll, Graphite2_jll, Libffi_jll)
    JLLWrappers.@init_library_product(
        libharfbuzz,
        "lib/libharfbuzz.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libharfbuzz_gobject,
        "lib/libharfbuzz-gobject.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libharfbuzz_subset,
        "lib/libharfbuzz-subset.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module BenchConstruct
using BenchmarkTools
using CompressHashDisplace

suite = BenchmarkGroup()

DICTIONARY = joinpath(@__DIR__, "..", "assets", "american-english.txt")

dict = Dict{String, Int}()

for (line, word) in enumerate(readlines(DICTIONARY))
    dict[word] = line
end

suite["fd_construct"] = @benchmarkable FrozenDict($dict)
suite["fud_construct"] = @benchmarkable FrozenUnsafeDict($dict)

end # module

BenchConstruct.suite


text:
using BinDeps

@BinDeps.setup

libgumbo = library_dependency("libgumbo", aliases = ["libgumbo-1"])

provides(Sources,
         URI("http://jamesporter.me/julia/gumbo-1.0.tar.gz"),
         libgumbo,
         unpacked_dir="gumbo-1.0")

provides(Binaries, URI("https://cache.julialang.org/https://bintray.com/artifact/download/tkelman/generic/gumbo.7z"),
         libgumbo, unpacked_dir="usr$WORD_SIZE/bin", os = :Windows)

provides(BuildProcess,
         Autotools(libtarget="libgumbo.la"),
         libgumbo, os = :Unix)

@BinDeps.install Dict(:libgumbo => :libgumbo)


text:
"""
    CircularArrayBuffer{E, T, N}

Using a `N` dimension Array to simulate a circular buffer of `N-1` dimensional elements.
Here `E` is the type of element and `T` is same with the `eltype` of `E`.
Call `eltype(b::CircularArrayBuffer{E,T,N})` will return `E`.

# Examples

```julia-repl
julia> b = CircularArrayBuffer{Float64}(2)
0-element CircularArrayBuffer{Float64,Float64,1}

julia> push!(b, rand())
1-element CircularArrayBuffer{Float64,Float64,1}:
 0.9709012378596478

julia> push!(b, rand())
2-element CircularArrayBuffer{Float64,Float64,1}:
 0.9709012378596478
 0.4510778027035365

julia> push!(b, rand())
2-element CircularArrayBuffer{Float64,Float64,1}:
 0.4510778027035365
 0.6774251288208646

julia> b = CircularArrayBuffer{Array{Float64,2}}(3, (3,3))
3×3×0 CircularArrayBuffer{Array{Float64,2},Float64,3}

julia> push!(b, randn(3,3))
3×3×1 CircularArrayBuffer{Array{Float64,2},Float64,3}:
[:, :, 1] =
 -0.548592   0.926179  -1.40998
 -0.0888621  0.177208   0.342665
  0.0925987  1.18531    0.962738
```
"""
mutable struct CircularArrayBuffer{E, T, N} <: AbstractArray{T, N}
    buffer::Array{T, N}
    first::Int
    length::Int
    stepsize::Int
    CircularArrayBuffer{T}(capacity::Int) where T = new{T, T, 1}(Vector{T}(undef, capacity), 1, 0, 1)
    function CircularArrayBuffer{T}(capacity::Int, element_size::Vararg{Int, N}) where {T<:AbstractArray,N}
        ndims(T) == N || throw(DimensionMismatch("the ndims of the specified type $T doesn't math the length of element_size $element_size"))
        new{T, eltype(T), N+1}(Array{eltype(T), N+1}(undef, element_size..., capacity), 1, 0, *(element_size...))
    end
    CircularArrayBuffer{T}(capacity::Int, element_size::Tuple{Vararg{Int, N}}) where {T<:AbstractArray,N} = CircularArrayBuffer{T}(capacity, element_size...)
end

eltype(cb::CircularArrayBuffer{E, T, N}) where {E, T, N} = E
size(cb::CircularArrayBuffer{E, T, N}) where {E, T, N} = (size(cb.buffer)[1:N-1]..., cb.length)

for func in [:view, :getindex]
     @eval @__MODULE__() begin
        $func(cb::CircularArrayBuffer{E, T, 1}, i::Int) where {E, T} = $func(cb.buffer, _buffer_index(cb, i))
        $func(cb::CircularArrayBuffer{E, T, 2}, i::Int) where {E, T} = $func(cb.buffer, :, _buffer_index(cb, i))
        $func(cb::CircularArrayBuffer{E, T, 3}, i::Int) where {E, T} = $func(cb.buffer, :, :, _buffer_index(cb, i))
        $func(cb::CircularArrayBuffer{E, T, 4}, i::Int) where {E, T} = $func(cb.buffer, :, :, :, _buffer_index(cb, i))
        $func(cb::CircularArrayBuffer{E, T, N}, i::Int) where {E, T, N} = $func(cb.buffer, [(:) for _ in 1 : N-1]...,  _buffer_index(cb, i))

        $func(cb::CircularArrayBuffer{E, T, 1}, I::Vector{Int}) where {E, T} = $func(cb.buffer, map(i -> _buffer_index(cb, i), I))
        $func(cb::CircularArrayBuffer{E, T, 2}, I::Vector{Int}) where {E, T} = $func(cb.buffer, :, map(i -> _buffer_index(cb, i), I))
        $func(cb::CircularArrayBuffer{E, T, 3}, I::Vector{Int}) where {E, T} = $func(cb.buffer, :, :, map(i -> _buffer_index(cb, i), I))
        $func(cb::CircularArrayBuffer{E, T, 4}, I::Vector{Int}) where {E, T} = $func(cb.buffer, :, :, :, map(i -> _buffer_index(cb, i), I))
        $func(cb::CircularArrayBuffer{E, T, N}, I::Vector{Int}) where {E, T, N} = $func(cb.buffer, [(:) for _ in 1 : N-1]..., map(i -> _buffer_index(cb, i), I))

        $func(cb::CircularArrayBuffer{E, T, 1}, i::UnitRange{Int}) where {E, T} = $func(cb.buffer, _buffer_index(cb, i))
        $func(cb::CircularArrayBuffer{E, T, 2}, i::UnitRange{Int}) where {E, T} = $func(cb.buffer, :,  _buffer_index(cb, i))
        $func(cb::CircularArrayBuffer{E, T, 3}, i::UnitRange{Int}) where {E, T} = $func(cb.buffer, :, :,  _buffer_index(cb, i))
        $func(cb::CircularArrayBuffer{E, T, 4}, i::UnitRange{Int}) where {E, T} = $func(cb.buffer, :, :, :,  _buffer_index(cb, i))
        $func(cb::CircularArrayBuffer{E, T, N}, i::UnitRange{Int}) where {E, T, N} = $func(cb.buffer, [(:) for _ in 1 : N-1]...,  _buffer_index(cb, i))
    end
end

setindex!(cb::CircularArrayBuffer{E, T, 1}, data, i) where {E, T} = cb.buffer[_buffer_index(cb, i)] = data
setindex!(cb::CircularArrayBuffer{E, T, 2}, data, i) where {E, T} = cb.buffer[:, _buffer_index(cb, i)] = data
setindex!(cb::CircularArrayBuffer{E, T, 3}, data, i) where {E, T} = cb.buffer[:, :, _buffer_index(cb, i)] = data
setindex!(cb::CircularArrayBuffer{E, T, 4}, data, i) where {E, T} = cb.buffer[:, :, :, _buffer_index(cb, i)] = data
setindex!(cb::CircularArrayBuffer{E, T, N}, data, i) where {E, T, N} = cb.buffer[[(:) for _ in 1 : N-1]..., _buffer_index(cb, i)] = data

## kept only to show elements
## never manually get/set the inner elements because it is very slow
## see https://discourse.julialang.org/t/varargs-performance/13578
getindex(cb::CircularArrayBuffer{E, T, N}, I::Vararg{Int, N}) where {E, T, N} = getindex(cb.buffer, I[1:N-1]...,  _buffer_index(cb, I[end]))

capacity(cb::CircularArrayBuffer) = size(cb.buffer)[end]
length(cb::CircularArrayBuffer) = cb.length
isfull(cb::CircularArrayBuffer) = length(cb) == capacity(cb)
isempty(cb::CircularArrayBuffer) = length(cb) == 0
empty!(cb::CircularArrayBuffer) = (cb.length = 0; cb)

"""
    push!(cb::CircularArrayBuffer{E, T, N}, data::E) where {E, T, N}

Add an element to the back and overwrite front if full.
Make sure that `length(data) == cb.stepsize`
"""
@inline function push!(cb::CircularArrayBuffer{E, T, N}, data::AbstractArray{T}) where {E, T, N}
    # length(data) == cb.stepsize || throw(DimensionMismatch("the length of buffer's stepsize doesn't match the length of data, $(cb.stepsize) != $(length(data))"))
    # if full, increment and overwrite, otherwise push
    if cb.length == capacity(cb)
        cb.first = (cb.first == capacity(cb) ? 1 : cb.first + 1)
    else
        cb.length += 1
    end
    nxt_idx = _buffer_index(cb, cb.length)
    cb.buffer[cb.stepsize * (nxt_idx - 1) + 1: cb.stepsize * nxt_idx] = data
    cb
end

@inline function push!(cb::CircularArrayBuffer{E, T, 1}, data::T) where {E, T}
    if cb.length == capacity(cb)
        cb.first = (cb.first == capacity(cb) ? 1 : cb.first + 1)
    else
        cb.length += 1
    end
    cb.buffer[_buffer_index(cb, cb.length)] = data
    cb
end

@inline function _buffer_index(cb::CircularArrayBuffer, i::Int)
    n = capacity(cb)
    idx = cb.first + i - 1
    if idx > n
        idx - n
    else
        idx
    end
end

@inline function _buffer_index(cb::CircularArrayBuffer, i::UnitRange{Int})
    start = _buffer_index(cb, i.start)
    stop = _buffer_index(cb, i.stop)
    start ≤ stop ? (start:stop) : vcat(start:capacity(cb), 1:stop)
end


text:
using DataFrames, CSV

for (csv_filename, asset_class, append_arg) in (
    ("list_currencies.csv", "@asset Cash", false),
    ("list_bonds.csv", "@asset Bond", true),
    ("list_commodities.csv", "@asset Commodity", true)
)
    df = CSV.read(joinpath("tools", csv_filename), DataFrame, header=4)

    dropmissing!(df, "Alphabetic Code")
    unique!(df, :Name)

    nrows = size(df)[1]
    df = hcat(DataFrame(Macro = fill(asset_class, nrows)), combine(df, "Alphabetic Code" => :Code,
            :Name => x -> replace.(titlecase.(x; strict=false), r" |’" => ""),
            renamecols=false)
    )
    CSV.write(joinpath("src", "pkgdefaults.jl"), df, delim="  ", header = false, append=append_arg)
end


text:
using Test, MultiKDE, Distributions

# Hard-code test cases and expected results
hardcode_data = [[0.75762386,6.0,8.0,],
[0.40709311,2.0,8.0,],
[0.73019489,1.0,7.0,],
[0.9350646,1.0,8.0,],
[0.20622478,7.0,3.0,],
[0.82868765,4.0,9.0,],
[0.83275865,6.0,5.0,],
[0.06151981,6.0,5.0,],
[0.47064985,2.0,0.0,],
[0.81312999,5.0,4.0,],
[0.60976384,5.0,0.0,],
[0.58597128,0.0,9.0,],
[0.97580821,5.0,9.0,],
[0.21986135,2.0,4.0,],
[0.33835384,6.0,6.0,],
[0.52131305,9.0,8.0,],
[0.0651185,6.0,9.0,],
[0.75884609,9.0,3.0,],
[0.59935732,5.0,4.0,],
[0.88377118,6.0,3.0,],]
hardcode_data_pred = [[0.11249392,4.0,6.0,],
[0.28704927,5.0,3.0,],
[0.06442488,1.0,6.0,],
[0.88971842,6.0,9.0,],
[0.58269063,5.0,2.0,],
[0.04070221,3.0,4.0,],
[0.33534097,2.0,3.0,],
[0.25188701,3.0,9.0,],
[0.54533024,5.0,9.0,],
[0.51731902,2.0,2.0,],
[0.40469525,5.0,1.0,],
[0.34114321,8.0,4.0,],
[0.39484506,7.0,9.0,],
[0.45614159,6.0,3.0,],
[0.01384118,2.0,2.0,],
[0.09847895,3.0,5.0,],
[0.95260732,8.0,3.0,],
[0.92190954,2.0,1.0,],
[0.10851023,7.0,7.0,],
[0.41038077,2.0,3.0,],]

# Expected result of above input
hardcode_expected_result = [
    [4.014556252605015e-07,4.93587849306923e-07,5.017127433753062e-05,4.297840888944418,1.4103686336015966,0.0004239071556898233,0.0002450830674101923,0.001520097546290904,1.6695504130319276e-05,9.585990125997658e-05,0.0024919490954259764,0.024667731573921547,1.2143707817645687e-05,8.952814734339356e-05,3.470292679120249e-09,9.825400944075402e-06,0.00017389979047319838,0.10795426189791625,2.5029651369084925e-05,4.859412084689039,],
    [0.006731618223323871,0.034852108047552656,0.0008987511262355018,0.26494502607810944,0.002093798643830388,0.04134463379287447,0.001734435606041325,0.034237993437472795,0.34814274656951194,0.0014581092899801078,0.002014698187056216,0.002549489235012434,0.03711991361931681,0.32252782622899556,0.0013659916119672653,0.0024460021110439176,0.06131450492717148,0.0012846078414541185,0.000750247341719304,0.0017450935407583167,],
    [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,],
    [0.07012047583564995,0.16682618301528643,0.04595196932831633,0.26088694468743184,0.15452574048888829,0.09336064500373596,0.11310195315610574,0.09485512130103643,0.3361451517102332,0.10902096724118503,0.13788547837516862,0.05885101643100109,0.1329887267194227,0.3387753430435409,0.0787994864890616,0.052985136882964974,0.10988250614008524,0.07668631403187032,0.06941533271304429,0.11576170881780576,],
]
hardcode_dims = [ContinuousDim(), CategoricalDim(10), UnorderedCategoricalDim(10)]
hardcode_bws = [[0.01, 0.1, 1], [1, 0.2, 0.01], [0.02, 1, 1], [0.5, 0.5, 0.5]]

# Arguments for sanity tests
continuous_range = (0, 200)
categorical_level = 10
unordered_extra_level = 10
# zero_val_tolerance = 1e-2
s1_tol = 1
univ_sanity_trial, multi_sanity_trial = 5, 1
num_s2_test = 100
num_obs = 20
dim_types = [ContinuousDim, CategoricalDim, UnorderedCategoricalDim]
bw_s1 = 0.1
bw_s2 = 0.1
diff_delta = 2e-1
multi_test_num_dims = [1, 2, 3, 4, 5]
min_range = 1
# Objects for sanity tests
unordered_objs = Vector{Any}([Vector{Int},  Vector{Float64}, Vector{Array}, Vector{Vector}, Vector{Pair}, Vector{ErrorException}, Vector{Tuple}])
append!(unordered_objs, 1:unordered_extra_level)
unordered_level = length(unordered_objs)

function randint(n_mx)
    ceil(Int,rand()*n_mx)
end

function randint(N, n_mx)
    [ceil(Int,rand()*n_mx) for _ in 1:N]
end

function gen_batch_data(dim_type::DimensionType, n)
    dt = []
    for _ in 1:n
        _dt = gen_data(dim_type)
        if !(_dt isa Real)
            _dt = findall(x->x===_dt, unordered_objs)[1]
        end
        push!(dt, _dt)
    end
    dt
end

function gen_batch_data(dim_types::Vector, n)
    dt = []
    for _ in 1:n
        push!(dt, [gen_data(dim_types[i]) for i in 1:length(dim_types)])
    end
    dt
end

function gen_data(dim_type)
    local _dt
    if dim_type isa CategoricalDim
        _dt = rand(1:categorical_level)
    elseif dim_type isa UnorderedCategoricalDim
        _dt = unordered_objs[rand(1:length(unordered_objs))]
    elseif dim_type isa ContinuousDim
        _dt = rand(Uniform(continuous_range[1], continuous_range[2]))
    else
        error("Nonsupported dim_type")
    end
    _dt
end

# Check if the point is local maximum(gradient=0) using difference differential
function is_maximum(multi_kde::KDEMulti, x; delta=1e-10)
    x1 = [_kde.type isa ContinuousDim ? _x-delta :
                            (_kde.type isa UnorderedCategoricalDim ? _x : _x-delta)
                            for (_kde, _x) in zip(multi_kde.KDEs, x)]
    x2 = [_kde.type isa ContinuousDim ? _x+delta :
                            (_kde.type isa UnorderedCategoricalDim ? _x : _x+delta)
                            for (_kde, _x) in zip(multi_kde.KDEs, x)]
    val, val_x1, val_x2 = MultiKDE.pdf(multi_kde, x), MultiKDE.pdf(multi_kde, x1), MultiKDE.pdf(multi_kde, x2)
    (val-val_x1 >= 0) && (val-val_x2 >= 0)
end

function is_maximum(kde::KDEUniv, x; delta=1e-10)
    x1 = kde.type isa ContinuousDim ? x-delta : (kde.type isa UnorderedCategoricalDim ? x-delta : x-delta)
    x2 = kde.type isa ContinuousDim ? x+delta : (kde.type isa UnorderedCategoricalDim ? x+delta : x+delta)
    val, val_x1, val_x2 = MultiKDE.pdf(kde, x, keep_all=false), MultiKDE.pdf(kde, x1, keep_all=false), MultiKDE.pdf(kde, x2, keep_all=false)
    (val-val_x1 >= 0) && (val-val_x2 >= 0)
end

# Same as above, but calculating the gradient with direction that x points to y
function kde_gradient_to_point(multi_kde::KDEMulti, x, y; delta=1e-10)
    vec = [_kde.type isa ContinuousDim ? _x-_y :
                            (_kde.type isa UnorderedCategoricalDim ? nothing : _x-_y)
                            for (_kde, _x, _y) in zip(multi_kde.KDEs, x, y)]
    x1 = [_vec===nothing ? _x : _x-(_vec*delta) for (_x, _vec) in zip(x, vec)]
    x2 = [_vec===nothing ? _x : _x+(_vec*delta) for (_x, _vec) in zip(x, vec)]
    (MultiKDE.pdf(multi_kde, x2) - MultiKDE.pdf(multi_kde, x1)) / delta
end

function kde_gradient_to_point(kde::KDEUniv, x, y; delta=1e-10)
    vec = kde.type isa ContinuousDim ? x-y : (kde.type isa UnorderedCategoricalDim ? nothing : x-y)
    x1 = vec===nothing ? x : x-(vec*delta)
    x2 = vec===nothing ? x : x+(vec*delta)
    (MultiKDE.pdf(kde, x2, keep_all=false) - MultiKDE.pdf(kde, x1, keep_all=false)) / delta
end

@testset "KDEUniv" begin
    println("Test set for KDEUniv")
    for _ in 1:univ_sanity_trial
        _dim = dim_types[randint(length(dim_types))]
        _dim = _dim==ContinuousDim ? _dim() : (_dim==CategoricalDim ? _dim(categorical_level) : _dim(unordered_level))
        dt = gen_batch_data(_dim, num_obs)
        # Filter observations that are too close
        local dt_s1
        if _dim isa ContinuousDim
            dt_s1 = filter(e->all([abs(e-_dt)>2 for _dt in dt if _dt≠e]), dt)
        else
            dt_s1 = dt
        end
        kde_s1 = KDEUniv(_dim, bw_s1, dt_s1)
        kde_s2 = KDEUniv(_dim, bw_s2, dt)
        # 1. Sanity-check 1: When bandwidth is small enough, the gradient of PDF of observations shoule be close to zero.
        println("Performing sanity check 1...", string(_dim))
        ## Gradient function
        is_maximum_s1 = [is_maximum(kde_s2, _dt, delta=diff_delta) for _dt in kde_s1.data]
        # println("------------------------------")
        # println(sort(kde_s1.data))
        # println([MultiKDE.pdf(kde_s1, _dt, keep_all=false) for _dt in sort(kde_s1.data)])
        # println(is_maximum_s1)
        # println("------------------------------")
        # Debug plotting
        # if sum(is_maximum_s1) < (length(is_maximum_s1)-s1_tol)
        #     x = LinRange(minimum(kde_s1.data), maximum(kde_s1.data), 20000)
        #     y = [MultiKDE.pdf(kde_s1, _x, keep_all=false) for _x in x]
        #     highest = maximum(y)
        #     plot(x, y)
        #     plot!(kde_s1.data, [highest+0.05 for _ in 1:length(y)], seriestype=:scatter, size=(900, 450))
        #     println(sort(kde_s1.data))
        #     println([MultiKDE.pdf(kde_s1, _dt, keep_all=false) for _dt in sort(kde_s1.data)])
        #     println(kde_s1.kernel)
        #     println(sort(kde_s1.data))
        #     println("Plotted")
        #     savefig(string(sum(is_maximum_s1), ".png"))
        # end
        @test sum(is_maximum_s1) >= (length(is_maximum_s1)-s1_tol)
        # 2. Sanity-check 2: Using a hyperbox to include all observations, for every vector out of box that points opposite to box center, the gradient should be negative.
        println("Performing sanity check 2...")
        mn, mx = minimum(kde_s2.data), maximum(kde_s2.data)
        middle = (mx-mn)/2
        range = min(mx-mn, min_range)
        grad_s2 = []
        for i in 1:num_s2_test
            point_out = rand([mn-range, mx+range])
            _grad = kde_gradient_to_point(kde_s2, point_out, middle, delta=diff_delta)
            push!(grad_s2, _grad)
        end
        println(grad_s2)
        @test sum((grad_s2 .≈ 0) .| (grad_s2 .< 0)) == length(grad_s2)
    end
end

@testset "KDEMulti" begin
    println("Test set for KDEMulti")
    # 1. Hard-code cases check
    # result = RealVectorVector()
    result = Vector()
    for bandwidth in hardcode_bws
        multi_kde = KDEMulti(hardcode_dims, bandwidth, hardcode_data)
        _result = [MultiKDE.pdf(multi_kde, _data) for _data in hardcode_data_pred]
        push!(result, _result)
    end
    @test sum(result .≈ hardcode_expected_result) == length(hardcode_expected_result)
    # Variables needed for Sanity-check
    # Try different combinations
    for num_dim in multi_test_num_dims
        for _ in 1:multi_sanity_trial
            println("Number of dimensions = ", num_dim)
            # Prepare Random KDEMulti
            _dims = dim_types[randint(num_dim, length(dim_types))]
            _dims = [_dim==ContinuousDim ? _dim() : (_dim==CategoricalDim ? _dim(categorical_level) : _dim(unordered_level)) for _dim in _dims]
            _bw_s1 = [bw_s1 for _ in 1:num_dim]
            _bw_s2 = [bw_s2 for _ in 1:num_dim]
            dt = gen_batch_data(_dims, num_obs)
            candidates = Dict{Int, Vector}()
            for (i, __dims) in zip(1:length(_dims), _dims)
                if __dims isa UnorderedCategoricalDim
                    candidates[i] = unordered_objs
                end
            end
            local dt_s1
            if any([_dim isa ContinuousDim for _dim in _dims])
                valid_idx = []
                for (_k, _dim) in zip(1:length(_dims), _dims)
                    if _dim isa ContinuousDim
                        _dim_dt = [_dt[_k] for _dt in dt]
                        _valid_idx = filter(_i->all([abs(_dim_dt[_i]-_dim_dt[_j])>2 for _j in 1:length(_dim_dt) if _dim_dt[_i]≠_dim_dt[_j]]), 1:length(_dim_dt))
                        push!(valid_idx, _valid_idx)
                    end
                end
                valid_idx = intersect(valid_idx...)
                dt_s1 = dt[valid_idx]
            else
                dt_s1 = dt
            end
            # println("------------------------------")
            # println(_dims)
            # println(_bw_s1)
            # println(_bw_s2)
            # println(dt_s1)
            # println(candidates)
            # println("------------------------------")
            multi_kde_s1 = KDEMulti(_dims, _bw_s1, dt_s1, candidates)
            multi_kde_s2 = KDEMulti(_dims, _bw_s2, dt, candidates)
            # 2. Sanity-check 1: When bandwidth is small enough, the gradient of PDF of observations shoule be close to zero.
            println("Performing sanity check 1...")
            ## Gradient function
            is_maximum_s1 = [is_maximum(multi_kde_s1, _dt, delta=diff_delta) for _dt in dt_s1]
            println(is_maximum_s1)
            # @test sum(grad_s1 .< zero_val_tolerance) >= (length(grad_s1)-zero_num_tolerance)
            @test sum(is_maximum_s1) >= (length(is_maximum_s1))
            # 3. Sanity-check 2: Using a hyperbox to include all observations, for every vector out of box that points opposite to box center, the gradient should be negative.
            println("Performing sanity check 2...")
            mins = [minimum(kde.data) for kde in multi_kde_s2.KDEs]
            maxs = [maximum(kde.data) for kde in multi_kde_s2.KDEs]
            middle = [_dim isa ContinuousDim ? (_max-_min) / 2 :
                        (_dim isa UnorderedCategoricalDim ? (multi_kde_s2.index_to_unordered[_kde][round(Int, (_max+_min)/2)]) : round(Int, (_max+_min)/2))
                        for (i, _dim, _kde, _min, _max) in zip(1:length(multi_kde_s2.dims), multi_kde_s2.dims, multi_kde_s2.KDEs, mins, maxs)]
            ranges = [min((maxs[i]-mins[i]), min_range) for i in 1:length(mins)]
            grad_s2 = []
            for i in 1:num_s2_test
                point_out = [_dim isa ContinuousDim ? rand([mins[i]-rand(Uniform(0, ranges[i])), maxs[i]+rand(Uniform(0, ranges[i]))]) :
                                (_dim isa UnorderedCategoricalDim ? rand([multi_kde_s2.index_to_unordered[_kde][1], multi_kde_s2.index_to_unordered[_kde][length(multi_kde_s2.index_to_unordered[_kde])]]) :
                                                                            rand([mins[i]-rand(1:ranges[i]), maxs[i]+rand(1:ranges[i])]))
                                for (i, _dim, _kde) in zip(1:length(multi_kde_s2.dims), multi_kde_s2.dims, multi_kde_s2.KDEs)]
                _grad = kde_gradient_to_point(multi_kde_s2, point_out, middle, delta=diff_delta)
                push!(grad_s2, _grad)
            end
            println(grad_s2)
            @test sum((grad_s2 .≈ 0) .| (grad_s2 .< 0)) == length(grad_s2)
        end
    end
end


text:
module MPI

using Libdl, Serialization
using Requires
using DocStringExtensions

export mpiexec, UBuffer, VBuffer

function serialize(x)
    s = IOBuffer()
    Serialization.serialize(s, x)
    take!(s)
end

function deserialize(x)
    s = IOBuffer(x)
    Serialization.deserialize(s)
end

primitive type SentinelPtr Sys.WORD_SIZE
end

primitive type MPIPtr Sys.WORD_SIZE
end
@assert sizeof(MPIPtr) == sizeof(Ptr{Cvoid})
Base.cconvert(::Type{MPIPtr}, x::SentinelPtr) = x
Base.unsafe_convert(::Type{MPIPtr}, x::SentinelPtr) = reinterpret(MPIPtr, x)

function _doc_external(fname)
"""
- `$fname` man page: [OpenMPI](https://www.open-mpi.org/doc/current/man3/$fname.3.php), [MPICH](https://www.mpich.org/static/docs/latest/www3/$fname.html)
"""
end

try
    include(joinpath(dirname(@__DIR__), "deps", "deps.jl"))
    include(joinpath(dirname(@__DIR__), "deps", "compile_time_mpi_constants.jl"))
catch e
    error("MPI.jl not properly configured, please run `Pkg.build(\"MPI\")`.")
end
include("define_load_time_mpi_constants.jl")
include("read_load_time_mpi_constants.jl")
include("implementations.jl")
include("error.jl")
include("info.jl")
include("group.jl")
include("comm.jl")
include("environment.jl")
include("datatypes.jl")
include("buffers.jl")
include("operators.jl")
include("nonblocking.jl")
include("pointtopoint.jl")
include("collective.jl")
include("topology.jl")
include("onesided.jl")
include("io.jl")
include("errhandler.jl")
include("mpiexec_wrapper.jl")

include("deprecated.jl")

function __init__()

    @static if Sys.isunix()
        # dlopen the MPI library before any ccall:
        # - RTLD_GLOBAL is required for Open MPI
        #   <https://www.open-mpi.org/community/lists/users/2010/04/12803.php>
        # - also allows us to ccall global symbols, which enables
        #   profilers which use LD_PRELOAD
        # - don't use RTLD_DEEPBIND; this leads to segfaults at least
        #   on Ubuntu 15.10
        #   <https://github.com/JuliaParallel/MPI.jl/pull/109>
        Libdl.dlopen(libmpi, Libdl.RTLD_LAZY | Libdl.RTLD_GLOBAL)
    end

    __init__deps()

    read_load_time_mpi_constants()

    # disable UCX memory cache, since it doesn't work correctly
    # https://github.com/openucx/ucx/issues/5061
    if !haskey(ENV, "UCX_MEMTYPE_CACHE")
        ENV["UCX_MEMTYPE_CACHE"] = "no"
    end

    # Julia multithreading uses SIGSEGV to sync threads
    # https://docs.julialang.org/en/v1/devdocs/debuggingtips/#Dealing-with-signals-1
    # By default, UCX will error if this occurs (issue #337)
    if !haskey(ENV, "UCX_ERROR_SIGNALS")
        # default is "SIGILL,SIGSEGV,SIGBUS,SIGFPE"
        ENV["UCX_ERROR_SIGNALS"] = "SIGILL,SIGBUS,SIGFPE"
    end

    @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" include("cuda.jl")
end

end


text:
#= get a hires esri shaded relief image.

We are going to download a static, hires shaded relief image from the ESRI servers.

By default basemap in python will call
http://server.arcgisonline.com/ArcGIS/rest/services/World_Shaded_Relief/MapServer/export?bbox=-2045969.4173062348,-2567274.0527822566,2546397.796126721,1048917.582910654&bboxSR=2163&imageSR=2163&size=5000,3937&dpi=96&format=png32&f=image
for our iamge.

But this only gets 4096 pixels in xpixels because of esri limits.  But we can tile ourselves...

We are going to tile for a 10k by Y image, where Y is computed so that
Y = floor(Int,10000*0.7874352088211344)

Then generate calls to download with chrome.

These will be named

ersi-shaded-relief-veryhires-1
ersi-shaded-relief-veryhires-2
...
ersi-shaded-relief-veryhires-6

as there are 6 of them.

To do this, we'll need to adjust the xmin and ymin

The standard call
#http://server.arcgisonline.com/ArcGIS/rest/services/World_Shaded_Relief/MapServer/export?bbox=-2045969.4173062348,-2567274.0527822566,2546397.796126721,1048917.582910654&bboxSR=2163&imageSR=2163&size=2600,4094&dpi=96&format=png32&f=image
=#

xmin,ymin,xmax,ymax = -2045969.4173062348,-2567274.0527822566,2546397.796126721,1048917.582910654
println("$xmin,$ymin,$(xmin + (xmax-xmin)/2),$ymax")
println("$(xmin + (xmax-xmin)/2),$ymin,$xmax,$ymax")

## Montage them
# http://www.imagemagick.org/Usage/montage/#montage
cmd = `montage esri-shaded-relief-veryhires-1.png esri-shaded-relief-veryhires-2.png -geometry +0+0 -tile 2x1 out.png`
run(cmd)

## This shows that everything works for a 2x1 tile. Now let's try our general tile.
aspect = 0.7874352088211344
xsize = 10005 # chosen so we get xsize/3 and ysize/2 are integers
ysize = floor(Int,xsize*aspect)
# we want stuff that's easy
@show xsize/3, ysize/2

xvals = range(xmin,xmax,length=4) # length = +1 because we want all the intermediate
yvals = range(ymin,ymax,length=3)
xtile = Int(xsize/3)
ytile = Int(ysize/2)

curimage = 1
for j=1:length(yvals)-1
  for i=1:length(xvals)-1
    global curimage
    url="http://server.arcgisonline.com/ArcGIS/rest/services/World_Shaded_Relief/MapServer/export?bbox=$(xvals[i]),$(yvals[j]),$(xvals[i+1]),$(yvals[j+1])&bboxSR=2163&imageSR=2163&size=$xtile,$ytile&dpi=96&format=png32&f=image"
    display(url)
    download(url, "ersi-shaded-relief-veryhires-$curimage.png")
    curimage += 1
  end
end

##
cmd = `montage ersi-shaded-relief-veryhires-4.png ersi-shaded-relief-veryhires-5.png ersi-shaded-relief-veryhires-6.png ersi-shaded-relief-veryhires-1.png ersi-shaded-relief-veryhires-2.png ersi-shaded-relief-veryhires-3.png  -geometry +0+0 -tile 3x2 esri-shaded-relief-veryhires.png`
run(cmd)

#= Old initial notes

We are going to setup code to cache the arcgis shaded relief output as this
can take a while to download and fails sometime.

When I used to use the bm.arcgisimage(...) call that is now commented out,
it would download the following picture.

http://server.arcgisonline.com/ArcGIS/rest/services/World_Shaded_Relief/MapServer/export?bbox=-2045969.4173062348,-2567274.0527822566,2546397.796126721,1048917.582910654&bboxSR=2163&imageSR=2163&size=5000,3937&dpi=96&format=png32&f=image

based on the code

"""
if not hasattr(self,'epsg'):
            msg = dedent("""
            Basemap instance must be creating using an EPSG code
            (http://spatialreference.org) in order to use the wmsmap method""")
            raise ValueError(msg)
        ax = kwargs.pop('ax', None) or self._check_ax()
        # find the x,y values at the corner points.
        p = pyproj.Proj(init="epsg:%s" % self.epsg, preserve_units=True)
        xmin,ymin = p(self.llcrnrlon,self.llcrnrlat)
        xmax,ymax = p(self.urcrnrlon,self.urcrnrlat)
        if self.projection in _cylproj:
            Dateline =\
            _geoslib.Point(self(180.,0.5*(self.llcrnrlat+self.urcrnrlat)))
            hasDateline = Dateline.within(self._boundarypolyxy)
            if hasDateline:
                msg=dedent("""
                arcgisimage cannot handle images that cross
                the dateline for cylindrical projections.""")
                raise ValueError(msg)
        if self.projection == 'cyl':
            xmin = (180./np.pi)*xmin; xmax = (180./np.pi)*xmax
            ymin = (180./np.pi)*ymin; ymax = (180./np.pi)*ymax
        # ypixels not given, find by scaling xpixels by the map aspect ratio.
        if ypixels is None:
            ypixels = int(self.aspect*xpixels)
        # construct a URL using the ArcGIS Server REST API.
        basemap_url = \
"%s/rest/services/%s/MapServer/export?\
bbox=%s,%s,%s,%s&\
bboxSR=%s&\
imageSR=%s&\
size=%s,%s&\
dpi=%s&\
format=png32&\
transparent=true&\
f=image" %\
(server,service,xmin,ymin,xmax,ymax,self.epsg,self.epsg,xpixels,ypixels,dpi)
        # print URL?
        if verbose: print(basemap_url)
        # return AxesImage instance.
        return self.imshow(imread(urlopen(basemap_url)),ax=ax,
                           origin='upper')
"""

http://server.arcgisonline.com/ArcGIS/rest/services/World_Shaded_Relief/MapServer/export?bbox=-2045969.4173062348,-2567274.0527822566,2546397.796126721,1048917.582910654&bboxSR=2163&imageSR=2163&size=5000,3937&dpi=96&format=png32&f=image

I wanted a 7500 pixel version, so I need the aspect ratio to adjust ypixels.

This involves finding the aspect to compute the ypixels, so I got the aspect

0.7874352088211344

http://server.arcgisonline.com/ArcGIS/rest/services/World_Shaded_Relief/MapServer/export?bbox=-2045969.4173062348,-2567274.0527822566,2546397.796126721,1048917.582910654&bboxSR=2163&imageSR=2163&size=7500,5905&dpi=96&format=png32&f=image

=#

max_stars_repo_name: dgleich/flowpaper-code

max_issues_repo_name: dgleich/flowpaper-code

max_forks_repo_name: dgleich/flowpaper-code


text:
__precompile__()

module InteractBulma

using Reexport
using WebIO
@reexport using InteractBase

import InteractBase: dropdown

export Bulma

struct Bulma<:InteractBase.WidgetTheme; end

const custom_css = joinpath(@__DIR__, "..", "assets", "custom.css")
const main_css = joinpath(@__DIR__, "..", "assets", "main.css")
const bulma_min_css = joinpath(@__DIR__, "..", "assets", "bulma.min.css")
const bulma_slider_min_css = joinpath(@__DIR__, "..", "assets", "bulma-slider.min.css")
const bulma_switch_min_css = joinpath(@__DIR__, "..", "assets", "bulma-switch.min.css")
const bulma_accordion_min_css = joinpath(@__DIR__, "..", "assets", "bulma-accordion.min.css")
const bulma_checkradio_min_css = joinpath(@__DIR__, "..", "assets", "bulma-checkradio.min.css")
const bulma_tooltip_min_css = joinpath(@__DIR__, "..", "assets", "bulma-tooltip.min.css")
const font_awesome = InteractBase.font_awesome

function InteractBase.libraries(::Bulma)
    bulmalibs = InteractBase.isijulia() ?
        [main_css] :
        [
            custom_css,
            bulma_min_css,
            bulma_slider_min_css,
            bulma_switch_min_css,
            bulma_accordion_min_css,
            bulma_checkradio_min_css,
            bulma_tooltip_min_css,
        ]
    vcat(font_awesome, bulmalibs)
end

function __init__()
    InteractBase.settheme!(Bulma())
    nothing
end

include("widgets.jl")

end # module


text:
Base.broadcastable(model::EoSModel) = Ref(model)

"""
    eos(model::EoSModel, V, T, z=SA[1.0])

basic OpenSAFT function, returns the total Helmholtz free energy.

# Inputs:
- `model::EoSModel` Thermodynamic model to evaluate
- `V` Total volume, in [m³]
- `T` Temperature, in [K]
- `z` mole amounts, in [mol], by default is `@SVector [1.0]`

# Outputs:

- Total Helmholtz free energy, in [J]

by default, it calls `R̄*T*∑(z)*(a_ideal(ideal_model,V,T,z) + a_res(model,V,T,z))` where `ideal_model == idealmodel(model)`, where `a_res` is the reduced residual Helmholtz energy and `a_ideal` is the reduced ideal Helmholtz energy.
You can mix and match ideal models if you provide:
- `idealmodel(model)`: extracts the ideal model from your Thermodynamic model
- `a_res(model,V,T,z)`: residual reduced Helmholtz free energy

"""
function eos(model::EoSModel, V, T, z=SA[1.0])
    return N_A*k_B*∑(z)*T * (a_ideal(idealmodel(model),V,T,z)+a_res(model,V,T,z))
end
"""
    idealmodel(model::EoSModel)

retrieves the ideal model from the input's model.

# Examples:

```julia-repl
julia> pr = PR(["water"],idealmodel=IAPWS95Ideal)
PR{IAPWS95Ideal} with 1 component:
 "water"
Contains parameters: a, b, acentricfactor, Tc, Mw
julia> OpenSAFT.idealmodel(pr)
IAPWS95Ideal()
```

"""
idealmodel(model::EoSModel) = model.idealmodel

function eos(model::IdealModel, V, T, z=SA[1.0])
    return N_A*k_B*sum(z)*T * a_ideal(model,V,T,z)
end

# function eos(model::CubicModel, V, T, z)
#     return N_A*k_B*sum(z)*T * a_tot(model,V,T,z)
# end

"""
    component_names(model)::Vector{tring}

returns a vector of strings of each component.
"""
component_names(model) = model.components

function v_rackett(model,T)
    tc = only(paramvals(model.params.Tc))
    pc = only(paramvals(model.params.Pc))
    vc = only(paramvals(model.params.Vc))
    TT = promote_type(typeof(tc),typeof(T))
    zc = pc*vc/(R̄*tc)
    _2_7 = TT(2//7)
    _1 = TT(1.0)
    R̄*tc/pc*zc^(_1 + (_1 - T/tc)^(_2_7))
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    bat_default(f::Function, argname::Symbol, objective...)
    bat_default(f::Function, argname::Val, objective...)

Get the default value for argument `argname` of function `f` to use
for `objective`(s).

`objective`(s) are mandatory arguments of function `f` that semantically
constitute it's main objective(s), and that that a good default choice of
optional arguments (e.g. choice of algorithm(s), etc.) may depend on.
Which arguments are considered to be objectives is function-specific.

For example:

```julia
bat_default(bat_sample, :algorithm, density::PosteriorDensity) == MetropolisHastings()
bat_default(bat_sample, Val(:algorithm), samples::DensitySampleVector) == OrderedResampling()
```
"""
function bat_default end
export bat_default

@inline bat_default(f::Function, argname::Symbol, objective...) = bat_default(f, Val{argname}(), objective...)

"""
    argchoice_msg(f::Function, argname::Val, x)

*BAT-internal, not part of stable public API.*

Generates an information message regarding the choice of value `x` for
argument `argname` of function `f`.

The value `x` will often be the result of [`bat_default`](@ref).
"""
function argchoice_msg end

function bat_default_withinfo(f::Function, argname::Val, objective...)
    default = bat_default(f::Function, argname::Val, objective...)
    @info argchoice_msg(f, argname::Val, default)
    default
end

function bat_default_withdebug(f::Function, argname::Val, objective...)
    default = bat_default(f::Function, argname::Val, objective...)
    @debug argchoice_msg(f, argname::Val, default)
    default
end

result_with_args(r::NamedTuple) = merge(r, (optargs = NamedTuple(), kwargs = NamedTuple()))

result_with_args(r::NamedTuple, optargs::NamedTuple) = merge(r, (optargs = optargs, kwargs = NamedTuple()))

result_with_args(r::NamedTuple, optargs::NamedTuple, kwargs::NamedTuple) = merge(r, (optargs = optargs, kwargs = kwargs))

result_with_args(r::NamedTuple, optargs::NamedTuple, kwargs::Base.Iterators.Pairs) = result_with_args(r, optargs, values(kwargs))


text:
using Documenter, SpikingNNFunctions

makedocs(;
    modules=[SpikingNNFunctions],
    format=Documenter.HTML(
            canonical = "https://darsnack.github.io/SpikingNNFunctions.jl/stable"
        ),
    pages=[
        "Home" => "index.md",
        "API" => [
            "Neurons" => "neurons.md",
            "Synapses" => "synapses.md",
            "Thresholds" => "thresholds.md",
        ],
    ],
    repo="https://github.com/darsnack/SpikingNNFunctions.jl/blob/{commit}{path}#L{line}",
    sitename="SpikingNNFunctions.jl",
    authors="Kyle Daruwalla",
    assets=String[],
)

deploydocs(;
    repo="github.com/darsnack/SpikingNNFunctions.jl",
)


text:
using tictactoe
using Test

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

max_stars_repo_name: HariharanJayashankar/tictactoe

max_issues_repo_name: HariharanJayashankar/tictactoe

max_forks_repo_name: HariharanJayashankar/tictactoe


text:
using SparseArrays
using LinearAlgebra
using MatrixDepot
using BenchmarkTools
using Metis
using UnicodePlots
using Profile
using PrettyTables
using ChainPartitioners
using ChainPartitioners: SparseHint

for mtx in [
            #"DIMACS10/chesapeake",
            "Boeing/pwtk",
            #"Schmid/thermal1",
            #"Rothberg/3dtube",
           ]
    A = permutedims(1.0 * sparse(mdopen(mtx).A))

    work_model = AffineWorkModel(0, 10, 1)
    net_model = AffineConnectivityModel(0, 10, 1, 100)
    sym_model = AffineMonotonizedSymmetricConnectivityModel(0, 0, 1, 100, 90)
    comm_model = AffinePrimaryConnectivityModel(0, 10, 1, 0, 100)
    loc_model = AffineSecondaryConnectivityModel(0, 10, 1, 0, 100)
    eps=0.1

    println()
    println()
    println(spy(A, maxwidth=50, maxheight=50, title="$mtx"))
    rows = []
    for (method, key, rng) in [
        #(DisjointPartitioner(EquiSplitter(), EquiSplitter()), "split_equal__split_equal", false),
        #(SymmetricPartitioner(EquiSplitter()), "split_equal__symmetric", false),
        #(DisjointPartitioner(BisectIndexBottleneckSplitter(work_model), EquiSplitter()), "split_work__split_equal", false),
        #(AlternatingPartitioner(BisectIndexBottleneckSplitter(work_model), MagneticPartitioner()), "split_work__map_local", true),
        #(AlternatingNetPartitioner(SparseHint(), BisectIndexBottleneckSplitter(net_model), FlipBisectIndexBottleneckSplitter(loc_model)), "split_nets__split_comm", false),
        #(AlternatingNetPartitioner(SparseHint(), BisectIndexBottleneckSplitter(net_model), FlipBisectIndexBottleneckSplitter(loc_model), BisectIndexBottleneckSplitter(comm_model), FlipBisectIndexBottleneckSplitter(loc_model)), "split_nets__split_comm__split_comm__split_comm", false),
        #(AlternatingNetPartitioner(SparseHint(), BisectIndexBottleneckSplitter(net_model), MagneticPartitioner()), "split_nets__map_local", true),
        #(AlternatingNetPartitioner(SparseHint(), BisectIndexBottleneckSplitter(net_model), MagneticPartitioner(), BisectIndexBottleneckSplitter(comm_model)), "split_nets__map_local__split_comm", true),
        #(AlternatingNetPartitioner(SparseHint(), BisectIndexBottleneckSplitter(net_model), GreedyBottleneckPartitioner(loc_model)), "split_nets__map_greedy", true),
        #(AlternatingNetPartitioner(SparseHint(), BisectIndexBottleneckSplitter(net_model), GreedyBottleneckPartitioner(loc_model), BisectIndexBottleneckSplitter(comm_model)), "split_nets__map_greedy__split_comm", true),
        #(AlternatingNetPartitioner(SparseHint(), BisectCostBottleneckSplitter(net_model, eps), FlipBisectCostBottleneckSplitter(loc_model, eps)), "split_nets__split_comm__approx", false),
        #(AlternatingNetPartitioner(SparseHint(), BisectCostBottleneckSplitter(net_model, eps), FlipBisectCostBottleneckSplitter(loc_model, eps), BisectCostBottleneckSplitter(comm_model, eps), FlipBisectCostBottleneckSplitter(loc_model, eps)), "split_nets__split_comm__split_comm__split_comm__approx", false),
        #(AlternatingNetPartitioner(SparseHint(), BisectCostBottleneckSplitter(net_model, eps), MagneticPartitioner()), "split_nets__map_local__approx", true),
        #(AlternatingNetPartitioner(SparseHint(), BisectCostBottleneckSplitter(net_model, eps), MagneticPartitioner(), BisectCostBottleneckSplitter(comm_model, eps)), "split_nets__map_local__split_comm__approx", true),
        #(AlternatingNetPartitioner(SparseHint(), BisectCostBottleneckSplitter(net_model, eps), GreedyBottleneckPartitioner(loc_model)), "split_nets__map_greedy__approx", true),
        #(AlternatingNetPartitioner(SparseHint(), BisectCostBottleneckSplitter(net_model, eps), GreedyBottleneckPartitioner(loc_model), BisectCostBottleneckSplitter(comm_model, eps)), "split_nets__map_greedy__split_comm__approx", true),
        #(AlternatingPartitioner(LazyBisectCostBottleneckSplitter(net_model, eps), LazyFlipBisectCostBottleneckSplitter(loc_model, eps)), "split_nets__split_comm__lazy__approx", false),
        #(PermutingPartitioner(CuthillMcKeePermuter(), SymmetricPartitioner(LazyBisectCostBottleneckSplitter(sym_model, eps))), "rcm_order__split_sym__symmetric__lazy__approx", false),
        (SymmetricPartitioner(LazyBisectCostBottleneckSplitter(sym_model, eps)), "split_sym__symmetric__lazy__approx", false),
        #(MetisPartitioner(weight=work_model), "metis", false)
        #(AlternatingPartitioner(LazyBisectCostBottleneckSplitter(net_model, eps), FlipBisectCostBottleneckSplitter(loc_model, eps), LazyBisectCostBottleneckSplitter(comm_model, eps), FlipBisectCostBottleneckSplitter(loc_model, eps)), "split_nets__split_comm__split_comm__split_comm__lazy__approx", false),
        #(AlternatingPartitioner(LazyBisectCostBottleneckSplitter(net_model, eps), MagneticPartitioner()), "split_nets__map_local__lazy__approx", true),
        #(AlternatingPartitioner(LazyBisectCostBottleneckSplitter(net_model, eps), MagneticPartitioner(), LazyBisectCostBottleneckSplitter(comm_model, eps)), "split_nets__map_local__split_comm__lazy__approx", true),
        #(AlternatingPartitioner(LazyBisectCostBottleneckSplitter(net_model, eps), GreedyBottleneckPartitioner(loc_model)), "split_nets__map_greedy__lazy__approx", true),
        #(AlternatingPartitioner(LazyBisectCostBottleneckSplitter(net_model, eps), GreedyBottleneckPartitioner(loc_model), LazyBisectCostBottleneckSplitter(comm_model, eps)), "split_nets__map_greedy__split_comm__lazy__approx", true),
    ]
        for K in [16]
            setup_time = time(@benchmark (Π, Φ) = partition_plaid($A, $K, $method))
            @profile begin
            (Π, Φ) = partition_plaid(A, K, method)
            end
            Profile.print()
            push!(rows, [K key setup_time])
            println(method)
        end
    end
    pretty_table(vcat(rows...), ["K", "method", "setuptime"])
end


text:
import Base: fetch, produce, consume
import Base: map, put!, push!, take!, isopen, close

include("constants.jl")
include("crc32.jl")
include("protocol.jl")
include("io.jl")
include("channels.jl")
include("client.jl")
include("requests.jl")