text:
using GNSSTools
using Statistics
using Random
using StatsBase
using LinearAlgebra
using ProgressMeter
using Base.Threads
using JLD
using PyPlot
pygui(true)
include("ch4_constellation_hist_structs.jl")

function benchmark_acquisition(sigtype::GNSSTools.SignalType, Ts, CN0s,
                               fd_rates, f_s; iterations=100,
                               include_phase_noise=true,
                               include_thermal_noise=true,
                               include_adc=true,
                               include_carrier=true,
                               h_parms=[1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15]./16.3676e6^2,
                               code_start_idx=missing, phi=missing, f_d=missing,
                               nADC=8, fd_ranges=missing, f_if=1.25e6, prn=1, M=1,
                               Pd_threshold=0.9, p_fa=1e-7,
                               doppler_bin_num=5)
    if sigtype.include_I && sigtype.include_Q
        B = max(sigtype.B_I, sigtype.B_Q)
    elseif sigtype.include_I && !sigtype.include_Q
        B = sigtype.B_I
    elseif !sigtype.include_I && sigtype.include_Q
        B = sigtype.B_Q
    else
        error("No CodeType defined in SigType struct.")
    end
    Pds = Array{Float64}(undef, length(Ts), length(CN0s), length(fd_rates))
    p = Progress(length(Ts)*length(CN0s)*length(fd_rates)*iterations, 1,
                 "Processing")
    @threads for i in 1:length(Ts)
        T = Ts[i]
        t_length = T*M
        replica = definereplica(sigtype, f_s, T)
        signal = definesignal(sigtype, f_s, t_length;
                              prn=prn, f_if=f_if, nADC=nADC,
                              include_phase_noise=include_phase_noise,
                              include_thermal_noise=include_thermal_noise,
                              include_adc=include_adc,
                              include_carrier=include_carrier,
                              receiver_h_parms=h_parms)
        for j in 1:length(CN0s)
            CN0 = CN0s[j]
            for k in 1:length(fd_rates)
                fd_range = fd_ranges[k]
                fd_rate = fd_rates[k]
                Pd = 0
                for iteration in 1:iterations
                    if ismissing(phi)
                        ϕ₀_truth = rand(0:0.0001:2π)
                    else
                        ϕ₀_truth = phi
                    end
                    if ismissing(f_d)
                        f_d_truth = rand(-fd_range:0.1:fd_range)
                    else
                        f_d_truth = f_d
                    end
                    if ismissing(code_start_idx)
                        n0_truth = rand(1:1:replica.sample_num)
                    else
                        n0_truth = code_start_idx
                    end
                    if include_phase_noise
                        new_phase_noise = true
                    else
                        new_phase_noise = false
                    end
                    if include_thermal_noise
                        new_thermal_noise = true
                    else
                        new_thermal_noise = false
                    end
                    definesignal!(signal; CN0=CN0, phi=ϕ₀_truth, f_d=f_d_truth,
                                  code_start_idx=n0_truth, fd_rate=fd_rate,
                                  new_phase_noise=new_phase_noise,
                                  new_thermal_noise=new_thermal_noise)
                    generatesignal!(signal)
                    acq_fd_range = (doppler_bin_num-1)/(2*T)
                    fd_center = round(f_d_truth*T, digits=0)/T
                    fd_course, n0_est, SNR_est, Pd_est, above_threshold,
                    corr_result = courseacquisition(signal,
                                                    replica,
                                                    prn;
                                                    fd_center=fd_center,
                                                    fd_range=acq_fd_range,
                                                    fd_rate=0,
                                                    return_corrresult=true,
                                                    M=M,
                                                    p_fa=p_fa)
                    if above_threshold
                        Pd += 1
                    end
                    next!(p)
                end
                Pds[i,j,k] = Pd / iterations
            end
        end
    end
    return Pds
end

function plot_data_for_T(i, fig, Ts, CN0s, fd_rates, Pds;
                         cmap=get_cmap("viridis"), letter=('a' + (i-1)),
                         plot_title="($(letter))")
    ax = fig.add_subplot(2, 2, i)
    for j in 1:length(Ts)
        T = Ts[j]
        color = cmap(float(j)/length(Ts))
        # 0Hz/s Doppler rate
        if i == 1
            label1 = string("T = $(floor(Int, T*1e3))ms, Low ", L"\dot{f}_d")
            ax.plot(CN0s, Pds[j,:,1], color=color, label=label1, linestyle="-")
        else
            ax.plot(CN0s, Pds[j,:,1], color=color, linestyle="-")
        end
        # Largest minimum Doppler rate
        if i == 2
            label2 = string("T = $(floor(Int, T*1e3))ms, High ", L"\dot{f}_d")
            ax.plot(CN0s, Pds[j,:,2], color=color, label=label2, linestyle=":")
        else
            ax.plot(CN0s, Pds[j,:,2], color=color, linestyle=":")
        end
    end
    xlabel("C/N₀ (dB⋅Hz)")
    ylabel(L"P_d")
    title(plot_title)
    return ax
end

directory = "/home/sjbilardi/Dropbox/Apps/Overleaf/MastersThesis/"
directory = ""
doppler_file_name = string(directory, "ch1_doppler_curves.jld")
doppler_raw_gps, doppler_rate_raw_gps, sig_freq = load(doppler_file_name,
                                                       "doppler_raw_gps",
                                                       "doppler_rate_raw_gps",
                                                       "freq")
doppler_raw_iridium, doppler_rate_raw_iridium = load(doppler_file_name,
                                                     "doppler_raw_iridium",
                                                     "doppler_rate_raw_iridium")
doppler_raw_starlink, doppler_rate_raw_starlink = load(doppler_file_name,
                                                       "doppler_raw_starlink",
                                                       "doppler_rate_raw_starlink")
doppler_raw_oneweb, doppler_rate_raw_oneweb = load(doppler_file_name,
                                                   "doppler_raw_oneweb",
                                                   "doppler_rate_raw_oneweb")

gps_hist = initiate_hist(doppler_raw_gps, doppler_rate_raw_gps; nbins=1000)
iridium_hist = initiate_hist(doppler_raw_iridium, doppler_rate_raw_iridium;
                             nbins=1000)
starlink_hist = initiate_hist(doppler_raw_starlink, doppler_rate_raw_starlink;
                              nbins=1000)
oneweb_hist = initiate_hist(doppler_raw_oneweb, doppler_rate_raw_oneweb;
                            nbins=1000)

iterations = 100
file_names = ["ch4_section3_acquisition_l1ca_constellation_comparison.jld",
              "ch4_section3_acquisition_l5q_constellation_comparison.jld"]
f_ss = [2.5e6, 25e6]  # Hz
f_ifs = [0e6, 0e6]  # Hz
phis = [missing, missing]
code_start_idxs = [missing, missing]
Tss = [[1, 2, 5, 10],
       [1, 20, 40, 80, 100]] .* 1e-3  # ms
CN0s = Array(range(25, 50, step=1))
sigtypes = [define_l1ca_code_type(maximum(Tss[1])),
            define_l5_code_type(;channel="Q")]

# items = [1, 2]
items = [2]
for i in items
    println("Will compare constellations for $(sigtypes[i].name) signal.")
end
for i in items
    f_s = f_ss[i]
    f_if = f_ifs[i]
    phi = phis[i]
    Ts = Tss[i]
    code_start_idx = code_start_idxs[i]
    sigtype = sigtypes[i]
    fd_ranges = [maximum(gps_hist.f_d),      maximum(gps_hist.f_d),
                 maximum(iridium_hist.f_d),  maximum(iridium_hist.f_d),
                 maximum(starlink_hist.f_d), maximum(starlink_hist.f_d),
                 maximum(oneweb_hist.f_d),   maximum(oneweb_hist.f_d)] .* (sigtype.sig_freq/sig_freq)
    fd_rates = [0, minimum(gps_hist.fd_rate),
                0, minimum(iridium_hist.fd_rate),
                0, minimum(starlink_hist.fd_rate),
                0, minimum(oneweb_hist.fd_rate)] .* (sigtype.sig_freq/sig_freq)
    fd_rates_gps = [0, minimum(gps_hist.fd_rate)] .* (sigtype.sig_freq/sig_freq)
    fd_rates_iridium = [0, minimum(iridium_hist.fd_rate)] .* (sigtype.sig_freq/sig_freq)
    fd_rates_starlink = [0, minimum(starlink_hist.fd_rate)] .* (sigtype.sig_freq/sig_freq)
    fd_rates_oneweb = [0, minimum(oneweb_hist.fd_rate)] .* (sigtype.sig_freq/sig_freq)
    Pds = benchmark_acquisition(sigtype, Ts, CN0s, fd_rates, f_s;
                                iterations=iterations, phi=phi,
                                code_start_idx=code_start_idx,
                                f_if=f_if,
                                fd_ranges=fd_ranges)
    Pds_gps = Pds[:,:,1:2]
    Pds_iridium = Pds[:,:,3:4]
    Pds_starlink = Pds[:,:,5:6]
    Pds_oneweb = Pds[:,:,7:8]
    file = string(directory, file_names[i])
    save(file,
         "sig_freq", sig_freq,
         "Ts", Ts,
         "CN0s", CN0s,
         "fd_rates_gps", fd_rates_gps,
         "fd_rates_iridium", fd_rates_iridium,
         "fd_rates_starlink", fd_rates_starlink,
         "fd_rates_oneweb", fd_rates_oneweb,
         "Pds_gps", Pds_gps,
         "Pds_iridium", Pds_iridium,
         "Pds_starlink", Pds_starlink,
         "Pds_oneweb", Pds_oneweb,
         "f_s", f_s,
         "iterations", iterations)
end

directory = "/home/sjbilardi/Dropbox/Apps/Overleaf/MastersThesis/"
file_names = ["ch4_section3_acquisition_l1ca_constellation_comparison",
              "ch4_section3_acquisition_l5q_constellation_comparison"]
items = [1, 2]
for i in items
    file_name = file_names[i]
    file = string(directory, file_name, ".jld")
    results = load(file, "Ts", "CN0s", "fd_rates_gps", "fd_rates_iridium",
                   "fd_rates_starlink", "fd_rates_oneweb", "Pds_gps",
                   "Pds_iridium", "Pds_starlink", "Pds_oneweb")
    Ts, CN0s, fd_rates_gps, fd_rates_iridium, fd_rates_starlink,
    fd_rates_oneweb, Pds_gps, Pds_iridium, Pds_starlink, Pds_oneweb = results
    fig = figure(figsize=(8, 4.5))
    ax1 = plot_data_for_T(1, fig, Ts, CN0s, fd_rates_gps, Pds_gps,
                          plot_title="(a) GPS")
    legend()
    ax2 = plot_data_for_T(2, fig, Ts, CN0s, fd_rates_iridium, Pds_iridium;
                          plot_title="(b) Iridium")
    legend()
    ax3 = plot_data_for_T(3, fig, Ts, CN0s, fd_rates_starlink, Pds_starlink;
                          plot_title="(c) Starlink")
    ax4 = plot_data_for_T(4, fig, Ts, CN0s, fd_rates_oneweb, Pds_oneweb;
                          plot_title="(d) OneWeb")
    subplots_adjust(wspace=0.3, hspace=0.525,
                    bottom=0.1, left=0.1, right=0.9, top=0.93)
    savefig(string(directory, "figures/", file_name, ".svg"), dpi=300)
end


text:
cd("/home/astupidbear/.julia/v0.4/JuLearn/src/ArrayFire/examples/mlp/regression")
ENV["JuLearnBackend"]="ArrayFire"
using JuLearn,ArrayFire#,Plots
@Eval lr=0.1 δ=0. momentum=0. batchsize=50 ▽max=0.1 λ=0.0
include("dataProvider.jl")
include("predict.jl")

inputBase,outputBase=dataProvider();
net=chain([:d,:h,:h,:l],[1 1;1 1000;1000 1;1 1],:regression);
@time err=backprop(net,AFArray(inputBase),AFArray(outputBase),nIter=300)

using PyPlot;
figure();plot(err)
figure();plot(inputBase',Array(predict(net,AFArray(inputBase)))')
hold("on");plot(inputBase',outputBase')

ext: jl


text:
module AsymptoticBoltzmann

using DiffRules
using SpecialFunctions
using QuadGK
using Interpolations
using CSV
using DifferentialEquations
using ODEInterfaceDiffEq
using Roots
using ForwardDiff

include("constants.jl")
include("standard_model/standard_model.jl")
include("thermodynamic_particles.jl")

include("models/abstract_model.jl")
include("models/kinetic_mixing.jl")

include("relic_density.jl")

export KineticMixing
export relic_density

end

max_stars_repo_name: LoganAMorrison/AsymptoticBoltzmann

max_issues_repo_name: LoganAMorrison/AsymptoticBoltzmann

max_forks_repo_name: LoganAMorrison/AsymptoticBoltzmann


text:
## @lc code=start
using LeetCode
count_even(num::Int) = sum(iseven(sum(j-'0' for j in string(i))) for i in 1:num)

## @lc code=end

## @lc test=start
@testset "2180.count-integers-with-even-digit-sum.jl" begin
    @test count_even(4) == 2
    @test count_even(30) == 14
end
## @lc test=end

## @lc info=start
# link: [solution to 2180](https://leetcode-cn.com/problems/count-integers-with-even-digit-sum/solution/python-yi-xing-dai-ma-by-rex-rs-1xty/)
## @lc info=end


text:
using NumericalMethodsforEngineers, Test

a = Float64[16 4 8;4 5 -4;8 -4 22]
b = Float64[4, 2, 5]
n = size(a, 1)
d = zeros(n)

c = a \ b
c |> display
println()

at, dt = ldlt(a, d)

bt = ldlfor(at, b)
bt |> display
println()

for i in 1:n
  at[i, :] = at[i, :] / dt[i]
end

bt = subbac(at, bt)
bt |> display
println()

@test bt == c


text:
# not for export:
const UnivariateFiniteUnion =
    Union{UnivariateFinite, UnivariateFiniteArray}

"""
    classes(d::UnivariateFinite)
    classes(d::UnivariateFiniteArray)

A list of categorial elements in the common pool of classes used to
construct `d`.

    v = categorical(["yes", "maybe", "no", "yes"])
    d = UnivariateFinite(v[1:2], [0.3, 0.7])
    classes(d) # CategoricalArray{String,1,UInt32}["maybe", "no", "yes"]

"""
classes(d::UnivariateFiniteUnion) = d.decoder.classes

"""
    levels(d::UnivariateFinite)

A list of the raw levels in the common pool of classes used to
construct `d`, equal to
`CategoricalArrays.DataAPI.unwrap.(classes(d))`.

    v = categorical(["yes", "maybe", "no", "yes"])
    d = UnivariateFinite(v[1:2], [0.3, 0.7])
    levels(d) # Array{String, 1}["maybe", "no", "yes"]

"""
Missings.levels(d::UnivariateFinite)  =
    CategoricalArrays.DataAPI.unwrap.(classes(d))

function Dist.params(d::UnivariateFinite)
    raw = raw_support(d) # reflects order of pool at instantiation of d
    pairs = tuple([unwrap.(d.decoder(r))=>d.prob_given_ref[r] for r in raw]...)
    levs = unwrap.(classes(d))
    return (levels=levs, probs=pairs)
end

# get the internal integer representations of the support
raw_support(d::UnivariateFiniteUnion) = collect(keys(d.prob_given_ref))

"""
    Dist.support(d::UnivariateFinite)
    Dist.support(d::UnivariateFiniteArray)

Ordered list of classes associated with non-zero probabilities.

    v = categorical(["yes", "maybe", "no", "yes"])
    d = UnivariateFinite(v[1:2], [0.3, 0.7])
    support(d) # CategoricalArray{String,1,UInt32}["maybe", "no"]

"""
Dist.support(d::UnivariateFiniteUnion) =
    map(d.decoder, raw_support(d))

# TODO: If I manually give a class zero probability, it will appear in
# support, which is probably confusing. We may need two versions of
# support - one which are the labels with corresponding keys in the
# dictionary, and one for the true mathematical support (which is the
# one exported)

# not exported:
sample_scitype(d::UnivariateFiniteUnion) = d.scitype

CategoricalArrays.isordered(d::UnivariateFinite) = isordered(classes(d))
CategoricalArrays.isordered(u::UnivariateFiniteArray) = isordered(classes(u))

## DISPLAY

_round_prob(p) = p
_round_prob(p::Union{Float32,Float64}) = round(p, sigdigits=3)

function Base.show(stream::IO, d::UnivariateFinite)
    pairs = Dist.params(d).probs
    arg_str = join(["$(pr[1])=>$(_round_prob(pr[2]))" for pr in pairs], ", ")
    print(stream, "UnivariateFinite{$(d.scitype)}($arg_str)")
end

Base.show(io::IO, mime::MIME"text/plain",
          d::UnivariateFinite) = show(io, d)

# in common case of `Real` probabilities we can do a pretty bar plot:
function Base.show(io::IO, mime::MIME"text/plain",
                   d::UnivariateFinite{<:Finite{K},V,R,P}) where {K,V,R,P<:Real}
    show_bars = false
    if K <= MAX_NUM_LEVELS_TO_SHOW_BARS &&
        all(>=(0), values(d.prob_given_ref))
        show_bars = true
    end
    show_bars || return show(io, d)
    s = support(d)
    x = string.(CategoricalArrays.DataAPI.unwrap.(s))
    y = pdf.(d, s)
    S = d.scitype
    plt = barplot(x, y, title="UnivariateFinite{$S}")
    show(io, mime, plt)
end

show_prefix(u::UnivariateFiniteArray{S,V,R,P,1}) where {S,V,R,P} =
    "$(length(u))-element"
show_prefix(u::UnivariateFiniteArray) = join(size(u),'x')

"""
    isapprox(d1::UnivariateFinite, d2::UnivariateFinite; kwargs...)

Returns `true` if and only if `d1` and `d2` have the same support and
the corresponding probabilities are approximately equal. The key-word
arguments `kwargs` are passed through to each call of `isapprox` on
probability pairs.  Returns `false` otherwise.

"""
function Base.isapprox(d1::UnivariateFinite, d2::UnivariateFinite; kwargs...)
    support1 = Dist.support(d1)
    support2 = Dist.support(d2)
    for c in support1
        c in support2 || return false
        isapprox(pdf(d1, c), pdf(d2, c); kwargs...) ||
            return false # pdf defined below
    end
    return true
end
function Base.isapprox(d1::UnivariateFiniteArray,
                       d2::UnivariateFiniteArray; kwargs...)
    support1 = Dist.support(d1)
    support2 = Dist.support(d2)
    for c in support1
        c in support2 || return false
        isapprox(pdf.(d1, c), pdf.(d2, c); kwargs...) ||
            return false
    end
    return true
end

# TODO: It would be useful to define == as well.

"""
    Dist.pdf(d::UnivariateFinite, x)

Probability of `d` at `x`.

    v = categorical(["yes", "maybe", "no", "yes"])
    d = UnivariateFinite(v[1:2], [0.3, 0.7])
    pdf(d, "yes")     # 0.3
    pdf(d, v[1])      # 0.3
    pdf(d, "no")      # 0.0
    pdf(d, "house")   # throws error

Other similar methods are available too:

    mode(d)    # CategoricalValue{String, UInt32} "maybe"
    rand(d, 5) # CategoricalArray{String,1,UInt32}["maybe", "no", "maybe", "maybe", "no"] or similar
    d = fit(UnivariateFinite, v)
    pdf(d, "maybe") # 0.25
    logpdf(d, "maybe") # log(0.25)

One can also do weighted fits:

    w = [1, 4, 5, 1] # some weights
    d = fit(UnivariateFinite, v, w)
    pdf(d, "maybe") ≈ 4/11 # true

See also `classes`, `support`.
"""
function Dist.pdf(
    d::UnivariateFinite{S,V,R,P},
    cv::CategoricalValue,
) where {S,V,R,P}
    return get(d.prob_given_ref, int(cv), zero(P))
end
Dist.pdf(d::UnivariateFinite{S,V}, c::V) where {S,V} = _pdf(d, c)
Dist.pdf(::UnivariateFinite{S,V}, ::Missing) where {S,V} = missing

# Avoid method ambiguity errors with Dist >= 0.24
Dist.pdf(d::UnivariateFinite{S,V}, c::V) where {S,V<:Real} = _pdf(d, c)

function _pdf(d::UnivariateFinite, c)
    _classes = classes(d)
    c in _classes || throw(DomainError("Value $c not in pool. "))
    pool = CategoricalArrays.pool(_classes)
    class = pool[get(pool, c)]
    return pdf(d, class)
end

Dist.logpdf(d::UnivariateFinite, cv::CategoricalValue) = log(pdf(d,cv))
Dist.logpdf(d::UnivariateFinite{S,V}, c::V) where {S,V} = log(pdf(d, c))
Dist.logpdf(::UnivariateFinite{S,V}, ::Missing) where {S,V} = missing

# Avoid method ambiguity errors with Dist >= 0.24
Dist.logpdf(d::UnivariateFinite{S,V}, c::V) where {S,V<:Real} = log(pdf(d, c))

function Dist.mode(d::UnivariateFinite)
    dic = d.prob_given_ref
    p = values(dic)
    max_prob = maximum(p)
    m = first(first(dic)) # mode, just some ref for now

    # `maximum` of any iterable containing `NaN` would return `NaN`
    # For this case the index `m` won't be updated in the loop below as relations
    # involving NaN as one of it's arguments always returns false
    # (e.g `==(NaN, NaN)` returns false)
    throw_nan_error_if_needed(max_prob)
    for (x, prob) in dic
        if prob == max_prob
            m = x
            break
        end
    end
    return d.decoder(m)
end

function throw_nan_error_if_needed(x)
    if isnan(x)
        throw(
            DomainError(
                NaN,
                "`mode` is invalid for `UnivariateFininite` distribution "*
                "with `pdf` containing `NaN`s"
            )
        )
    end
end

# mode(v::Vector{UnivariateFinite}) = mode.(v)
# mode(u::UnivariateFiniteVector{2}) =
#     [u.support[ifelse(s > 0.5, 2, 1)] for s in u.scores]
# mode(u::UnivariateFiniteVector{C}) where {C} =
#     [u.support[findmax(s)[2]] for s in eachrow(u.scores)]

"""
    _cumulative(d::UnivariateFinite)

Return the cumulative probability vector `C` for the distribution `d`,
using only classes in the support of `d`, ordered according to the
categorical elements used at instantiation of `d`. Used only to
implement random sampling from `d`. We have `C[1] == 0` and `C[end] ==
1`, assuming the probabilities have been normalized.

"""
function _cumulative(d::UnivariateFinite{S,V,R,P}) where {S,V,R,P<:Real}

    # the keys of `d` are in order; see constructor
    p = collect(values(d.prob_given_ref))
    K = length(p)
    p_cumulative = Array{P}(undef, K + 1)
    p_cumulative[1] = zero(P)
    for i in 2:K + 1
        p_cumulative[i] = p_cumulative[i-1] + p[i-1]
    end
    return p_cumulative
end

"""
_rand(rng, p_cumulative, R)

Randomly sample the distribution with discrete support `R(1):R(n)`
which has cumulative probability vector `p_cumulative` (see
[`_cummulative`](@ref)).

"""
function _rand(rng, p_cumulative, R)
    real_sample = rand(rng)*p_cumulative[end]
    K = R(length(p_cumulative))
    index = K
    for i in R(2):R(K)
        if real_sample < p_cumulative[i]
            index = i - R(1)
            break
        end
    end
    return index
end

function Base.rand(rng::AbstractRNG,
                   d::UnivariateFinite{<:Any,<:Any,R}) where R
    p_cumulative = _cumulative(d)
    return Dist.support(d)[_rand(rng, p_cumulative, R)]
end

function Base.rand(rng::AbstractRNG,
                   d::UnivariateFinite{<:Any,<:Any,R},
                   dim1::Int, moredims::Int...) where R # ref type
    p_cumulative = _cumulative(d)
    A = Array{R}(undef, dim1, moredims...)
    for i in eachindex(A)
        @inbounds A[i] = _rand(rng, p_cumulative, R)
    end
    support = Dist.support(d)
    return broadcast(i -> support[i], A)
end

rng(d::UnivariateFinite, args...) = rng(Random.GLOBAL_RNG, d, args...)

function Dist.fit(d::Type{<:UnivariateFinite},
                           v::AbstractVector{C}) where C
    C <: CategoricalValue ||
        error("Can only fit a UnivariateFinite distribution to samples of "*
              "`CategoricalValue type`. ")
    y = skipmissing(v) |> collect
    isempty(y) && error("No non-missing data to fit. ")
    N = length(y)
    count_given_class = Dist.countmap(y)
    classes = Tuple(keys(count_given_class))
    probs = values(count_given_class)./N
    prob_given_class = LittleDict(classes, probs)
    return UnivariateFinite(prob_given_class)
end

# TODO: Implement an MLJ style fit/transform with an option to specify
# Laplacian smoothing.

function Dist.fit(d::Type{<:UnivariateFinite},
                           v::AbstractVector{C},
                           weights::Nothing) where C
    return Dist.fit(d, v)
end

function Dist.fit(d::Type{<:UnivariateFinite},
                           v::AbstractVector{C},
                           weights::AbstractVector{<:Real}) where C
    C <: CategoricalValue ||
        error("Can only fit a UnivariateFinite distribution to samples of "*
              "`CategoricalValue` type. ")
    y = broadcast(identity, skipmissing(v))
    isempty(y) && error("No non-missing data to fit. ")
    classes_seen = filter(in(unique(y)), classes(y[1]))

    # instantiate and initialize prob dictionary:
    prob_given_class = LittleDict{C,Float64}()
    for c in classes_seen
        prob_given_class[c] = 0
    end

    # compute unnormalized  probablilities:
    for i in eachindex(y)
        prob_given_class[y[i]] += weights[i]
    end

    # normalize the probabilities:
    S = sum(values(prob_given_class))
    for c in keys(prob_given_class)
        prob_given_class[c] /=S
    end

    return UnivariateFinite(prob_given_class)
end

# # BROADCASTING OVER SINGLE UNIVARIATE FINITE

# This mirrors behaviour assigned Distributions.Distribution objects,
# which allows `pdf.(d::UnivariateFinite, support(d))` to work.

Broadcast.broadcastable(d::UnivariateFinite) = Ref(d)


text:
function extract_apriori_info(t::RelaxMulEnumInner, x::Vector{MC{N,T}}, Q::Int) where {N,T}
    cv = -Inf
    cc = Inf
    #DEBUG_NL && @show x.box
    for i = 1:Q
        z = x[i]
        cvt = z.cv
        cct = z.cc
        if cvt > cv
            cv = cvt
        end
        if cct < cc
            cc = cct
        end
    end
    return cv, cc
end
function extract_apriori_info(t::RelaxMulEnumInner, x::MCBoxPnt{Q,N,T}, y::MC{N,T}) where {Q,N,T}
    extract_apriori_info(t, x.box, Q)
end

_cut_info(t::RelaxMulEnumInner, v, z, x) = z
_cut_info(t::RelaxMulEnum, v, z, x) = z

relax_info(s::RelaxMulEnumInner, n::Int, t::T) where T = MCBoxPnt{2^n,n,T}
relax_info(s::RelaxMulEnum, n::Int, t::T) where T = MCBoxPnt{2^n,n,T}

const ENUM_OUTER_RND = 1E-9
function f_init!(::RelaxMulEnum, g::DAT, b::RelaxCache{V,N,T}) where {V,N,T}
    b.use_apriori_mul = false
    fprop!(Relax(), g, b)
    xp = copy(b.v.x)
    vlbd = _lbd(b.v)
    vubd = _ubd(b.v)
    xl = copy(vlbd)
    xu = copy(vubd)
    vlbd .-= ENUM_OUTER_RND
    vubd .+= ENUM_OUTER_RND
    for k = 1:_node_count(g)
        if _is_unlocked(b, k) && !_is_num(b, k)
            b._info[k].v = _set(b, k)
        end
    end
    s = _sparsity(g, 1)
    for i = 0:2^N-1
        s = last(bitstring(i), N)
        for (k,j) in enumerate(_sparsity(g, 1))
            b.v.x[j] = s[k] == '1' ? xl[j] : xu[j]
        end
        fprop!(Relax(), g, b)
        for k = _node_count(g):-1:1
            if _is_unlocked(b, k) && !_is_num(b, k)
                b._info[k][i+1] = _set(b, k)
            end
        end
    end
    b.v.x .= xp
    vlbd .= xl
    vubd .= xu
    b.use_apriori_mul = true
    fprop!(RelaxMulEnumInner(), g, b)
    return
end
function fprop!(t::RelaxMulEnum, g::DAT, b::RelaxCache{V,N,T}) where {V,N,T<:RelaxTag}
    vlbd = _lbd(b.v)
    vubd = _ubd(b.v)
    xl = copy(vlbd)
    xu = copy(vubd)
    vlbd .-= ENUM_OUTER_RND
    vubd .+= ENUM_OUTER_RND
    fprop!(RelaxMulEnumInner(false), g, b)
    for k = 1:_node_count(g)
        if _is_unlocked(b, k) && !_is_num(b, k)
            b._info[k].v = _set(b, k)
        end
    end
    vlbd .= xl
    vubd .= xu
    fprop!(RelaxMulEnumInner(true), g, b)
    return
end

fprop!(t::RelaxMulEnumInner, v::Val{EXP}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag} = fprop!(Relax(), v, g, b, k)
fprop!(t::RelaxMulEnumInner, v::Val{EXP10}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag} = fprop!(Relax(), v, g, b, k)
fprop!(t::RelaxMulEnumInner, v::Val{LOG}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag} = fprop!(Relax(), v, g, b, k)
fprop!(t::RelaxMulEnumInner, v::Val{LOG10}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag} = fprop!(Relax(), v, g, b, k)
fprop!(t::RelaxMulEnumInner, v::Val{POW}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag} = fprop!(Relax(), v, g, b, k)
fprop!(t::RelaxMulEnumInner, v::Val{MINUS}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag} = fprop!(Relax(), v, g, b, k)
fprop!(t::RelaxMulEnumInner, v::Val{DIV}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag} = fprop!(Relax(), v, g, b, k)
fprop!(t::RelaxMulEnumInner, v::Val{PLUS}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag} = fprop!(Relax(), v, g, b, k)

function fprop!(t::RelaxMulEnumInner, vt::Variable, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    i = _first_index(g, k)
    x = _val(b, i)
    l = _lbd(b, i)
    u = _ubd(b, i)
    #@show x, l, u
    z = _var_set(MC{N,T}, _rev_sparsity(g, i, k), x, x, l, u)
    #@show z
    if !_first_eval(b)
        z = z ∩ _interval(b, k)
    end
    #@show z
    _store_set!(b, z, k)
    return
end

function fprop_2!(t::RelaxMulEnumInner, v::Val{MULT}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    x = _child(g, 1, k)
    y = _child(g, 2, k)
    x_is_num = _is_num(b, x)
    y_is_num = _is_num(b, y)
    if !x_is_num && y_is_num
        z = _set(b, x)*_num(b, y)
    elseif x_is_num && !y_is_num
        z = _num(b, x)*_set(b, y)
    else
        xs = _set(b, x)
        ys = _set(b, y)
        xinfo = _info(b, x)
        yinfo = _info(b, y)
        xi = t.use_info ? xinfo.v : xs
        yi = t.use_info ? yinfo.v : ys
        xcvU, xccL = extract_apriori_info(t, xinfo, xs)
        ycvU, yccL = extract_apriori_info(t, yinfo, ys)
        xcv = xi.cv; ycv = yi.cv
        xcc = xi.cc; ycc = yi.cc
        xcvg = xi.cv_grad; ycvg = yi.cv_grad
        xccg = xi.cc_grad; yccg = yi.cc_grad
        wIntv = xs.Intv*ys.Intv
        za_l = McCormick.mult_apriori_kernel(xs, ys, wIntv, xcv, ycv, xcvU, ycvU, xcvg, ycvg)
        za_u = McCormick.mult_apriori_kernel(-xs, -ys, wIntv, -xcc, -ycc, -xccL, -yccL, -xccg, -yccg)
        z = (xs*ys) ∩ za_l ∩ za_u
    end
    z = _cut(z, _set(b, k), b.v, b.ϵ_sg, _sparsity(g, k), b.cut, false)
    _store_set!(b, z, k)
    return
end

#TODO: INFO FOR SET BASED ON PASS...
function fprop_n!(t::RelaxMulEnumInner, ::Val{MULT}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    ys = zero(MC{N,T})
    znum = one(Float64)
    b._mult_temp.v = one(MC{N,T})
    b._mult_temp.box .= one(MC{N,T})
    count = 0
    first_set = true
    for i in _children(g, k)
        if _is_num(b, i)
            znum = znum*_num(b, i)
        else
            xs = _set(b, i)
            xinfo = _info(b, i)
            xi = t.use_info ? xinfo.v : xs
            if !first_set
                xcvU, xccL = extract_apriori_info(t, xinfo, xs)
                ycvU, yccL = extract_apriori_info(t, b._mult_temp, ys)
                xcv = xi.cv;       ycv = b._mult_temp.v.cv
                xcc = xi.cc;       ycc = b._mult_temp.v.cc
                xcvg = xi.cv_grad; ycvg = b._mult_temp.v.cv_grad
                xccg = xi.cc_grad; yccg = b._mult_temp.v.cc_grad
                wIntv = xs.Intv*ys.Intv
                za_l = McCormick.mult_apriori_kernel(xs, ys, wIntv, xcv, ycv, xcvU, ycvU, xcvg, ycvg)
                za_u = McCormick.mult_apriori_kernel(-xs, -ys, wIntv, -xcc, -ycc, -xccL, -yccL, -xccg, -yccg)
                ys = (xs*ys) ∩ za_l ∩ za_u
                b._mult_temp.v *= ys
                b._mult_temp.box .*= xinfo.box
            else
                first_set = false
                b._mult_temp.v = xi
                b._mult_temp.box .= xinfo.box
                ys = xs
            end
        end
        count += 1
    end
    z = ys*znum
    z = _cut(z, _set(b, k), b.v, b.ϵ_sg, _sparsity(g, k), b.cut, false)
    _store_set!(b, z, k)
    return
end

#TODO:
function fprop!(t::RelaxMulEnumInner, v::Val{MULT}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    (_arity(g, k) == 2) ? fprop_2!(t, Val(MULT), g, b, k) : fprop_n!(t, Val(MULT), g, b, k)
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function Euler2()

	Fa = 1;
	Fb = 1;

	sum = 0;

	F = 0;

	while F<4e6
		F = Fa+Fb;
		Fa = Fb;
		Fb = F;

		if mod(F,2) == 0
			sum+=F
		end
	end

	return sum

end

@time Euler2()
println(Euler2())

max_stars_repo_name: gribeill/ProjectEuler

max_issues_repo_name: gribeill/ProjectEuler

max_forks_repo_name: gribeill/ProjectEuler


text:
include("../src/CBOR.jl")

using Base.Test
using CBOR
using DataStructures

# Taken (and modified) from Appendix A of RFC 7049

two_way_test_vectors = Dict(
    0 => hex2bytes("00"),
    1 => hex2bytes("01"),
    10 => hex2bytes("0a"),
    23 => hex2bytes("17"),
    24 => hex2bytes("1818"),
    25 => hex2bytes("1819"),
    100 => hex2bytes("1864"),
    1000 => hex2bytes("1903e8"),
    1000000 => hex2bytes("1a000f4240"),
    1000000000000 => hex2bytes("1b000000e8d4a51000"),
    18446744073709551615 => hex2bytes("1bffffffffffffffff"),
    -18446744073709551616 => hex2bytes("3bffffffffffffffff"),
    -1 => hex2bytes("20"),
    -10 => hex2bytes("29"),
    -100 => hex2bytes("3863"),
    -1000 => hex2bytes("3903e7"),

    0.0 => hex2bytes("fa00000000"),
    -0.0 => hex2bytes("fa80000000"),
    1.0 => hex2bytes("fa3f800000"),
    1.1 => hex2bytes("fb3ff199999999999a"),
    1.5 => hex2bytes("fa3fc00000"),
    65504.0 => hex2bytes("fa477fe000"),
    100000.0 => hex2bytes("fa47c35000"),
    3.4028234663852886e+38 => hex2bytes("fa7f7fffff"),
    1.0e+300 => hex2bytes("fb7e37e43c8800759c"),
    5.960464477539063e-8 => hex2bytes("fa33800000"),
    0.00006103515625 => hex2bytes("fa38800000"),
    -4.0 => hex2bytes("fac0800000"),
    -4.1 => hex2bytes("fbc010666666666666"),

    false => hex2bytes("f4"),
    true => hex2bytes("f5"),
    Null() => hex2bytes("f6"),
    Undefined() => hex2bytes("f7"),
    Simple(16) => hex2bytes("f0"),
    Simple(24) => hex2bytes("f818"),
    Simple(255) => hex2bytes("f8ff"),

    Pair(0, "2013-03-21T20:04:00Z") =>
        hex2bytes("c074323031332d30332d32315432303a30343a30305a"),
    Pair(1, 1363896240) => hex2bytes("c11a514b67b0"),
    Pair(1, 1363896240.5) => hex2bytes("c1fb41d452d9ec200000"),
    Pair(23, hex2bytes("01020304")) => hex2bytes("d74401020304"),
    Pair(24, hex2bytes("6449455446")) => hex2bytes("d818456449455446"),
    Pair(32, "http://www.example.com") =>
        hex2bytes("d82076687474703a2f2f7777772e6578616d706c652e636f6d"),

    UInt8[] => hex2bytes("40"),
    hex2bytes("01020304") => hex2bytes("4401020304"),

    "" => hex2bytes("60"),
    "a" => hex2bytes("6161"),
    "IETF" => hex2bytes("6449455446"),
    "\"\\" => hex2bytes("62225c"),
    "\u00fc" => hex2bytes("62c3bc"),
    "\u6c34" => hex2bytes("63e6b0b4"),

    [] => hex2bytes("80"),
    [1, 2, 3] => hex2bytes("83010203"),
    Any[1, [2, 3], [4, 5]] => hex2bytes("8301820203820405"),
    [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25] =>
        hex2bytes("98190102030405060708090a0b0c0d0e0f101112131415161718181819"),

    Dict() => hex2bytes("a0"),
    OrderedDict(1=>2, 3=>4) => hex2bytes("a201020304"),
    OrderedDict("a"=>1, "b"=>[2, 3]) => hex2bytes("a26161016162820203"),
    OrderedDict("a"=>"A", "b"=>"B", "c"=>"C", "d"=>"D", "e"=>"E") =>
        hex2bytes("a56161614161626142616361436164614461656145"),

    ["a", Dict("b"=>"c")] => hex2bytes("826161a161626163")
)

for (data, bytes) in two_way_test_vectors
    @test isequal(bytes, encode(data))
    @test isequal(data, decode(bytes))
end

bytes_to_data_test_vectors = Dict(
    hex2bytes("fa7fc00000") => NaN32,
    hex2bytes("fa7f800000") => Inf32,
    hex2bytes("faff800000") => -Inf32,
    hex2bytes("fb7ff8000000000000") => NaN,
    hex2bytes("fb7ff0000000000000") => Inf,
    hex2bytes("fbfff0000000000000") => -Inf
)

for (bytes, data) in bytes_to_data_test_vectors
    @test isequal(data, decode(bytes))
end

iana_test_vector = Dict(
    BigInt(18446744073709551616) => hex2bytes("c249010000000000000000"),
    BigInt(-18446744073709551617) => hex2bytes("c349010000000000000000")
)

for (data, bytes) in iana_test_vector
    @test isequal(bytes, encode(data))
    @test isequal(data, decode_with_iana(bytes))
end

indef_length_coll_test_vectors = Dict(
    Any[["Hello", " ", "world"], "Hello world", UTF8String] =>
        hex2bytes("7f6548656c6c6f612065776f726c64ff"),
    Any["Hello world".data, "Hello world".data, Array{UInt8, 1}] =>
        hex2bytes("5f18481865186c186c186f18201877186f1872186c1864ff"),
    Any[[1, 2.3, "Twiddle"], [1, 2.3, "Twiddle"], AbstractVector] =>
        hex2bytes("9f01fb40026666666666666754776964646c65ff"),
    Any[OrderedDict(1=>2, 3.2=>"3.2"), OrderedDict(1=>2, 3.2=>"3.2"),
        Associative] => hex2bytes("bf0102fb400999999999999a63332e32ff")
)

coll = []
function list_producer()
    for e in coll
        produce(e)
    end
end
function map_producer()
    for (k, v) in coll
        produce(k)
        produce(v)
    end
end

for (data, bytes) in indef_length_coll_test_vectors
    coll = data[1]
    task =
        if data[3] <: Associative
            Task(map_producer)
        else
            Task(list_producer)
        end
    @test isequal(bytes, encode(Pair(task, data[3])) )
    @test isequal(data[2], decode(bytes))
end


text:
#=
mri_trajectory.jl
based on mri_trajectory.m
Jing Dong
2020-03

using MIRT: image_geom_mri, ndgrid, mri_kspace_spiral, ir_mri_kspace_ga_radial, prompt
=#

export mri_trajectory

using Random: seed!
using Plots; default(markerstrokecolor=:auto, markersize=2)
using Test: @test

"""
    kspace, omega, wi = mri_trajectory( ; ktype, N, fov, arg_wi, kwargs...)

Generate kspace trajectory samples and density compensation functions.

option
* `ktype::Symbol` k-space trajectory type; default `:radial`
* `N::Dims` target image size; default (32,30)
* `fov` field of view in x and y (and z); default [250,250] mm
* `arg_wi` options to pass to `ir_mri_density_comp` - not yet done
* `kwargs` options for the specific trajectory

out
* `kspace [Nk 2|3]` kspace samples in units 1/fov
* `omega [Nk 2|3]` trajectory samples over [-π,π)
* `wi [Nk 1]` (optional) density compensation factors

trajectory types:
`:cartesian`
`:radial`
`:cart_y_2`
`:random`
`:half8`
`:epi_sin`
`:spiral0 :spiral1 :spiral3`
`:rosette3`
`:epi_under`
`:gads` (emulate golden-angle data sharing per winkelmann:07:aor)
"""
function mri_trajectory( ;
		ktype::Symbol = :radial,
		N::Dims = (32,30),
		fov::AbstractVector{<:Number} = 250 * ones(length(N)),
		wi::Union{AbstractArray{<:Real},Nothing} = nothing,
		arg_wi::Any = nothing, # not done
		kwargs...) # optional arguments for specific trajectories

	length(N) ∉ (2,3) && throw("only 2d and 3d done")

	# default density compensation, which works for ordinary Cartesian (only)
	wi = nothing

	# trajectory choices

	if (ktype == :cartesian)
		o1 = ((0:(N[1]-1))/N[1] .- 0.5)*2π
		o2 = ((0:(N[2]-1))/N[2] .- 0.5)*2π
		if (length(N) == 2)
			o1, o2 = ndgrid(o1, o2)
			omega = [o1[:] o2[:]]
		end

		if (length(N) == 3)
			o3 = ((0:(N[3]-1))/N[3] .- 0.5)*2π
			o1, o2, o3 = ndgrid(o1, o2, o3)
			omega = [o1[:] o2[:] o3[:]]
		end

		wi = 1 / prod(fov)
	end

	if (ktype == :epi_under)
		omega, wi = mri_trajectory_epi_under(N, fov ; kwargs...)
	end

	if (ktype == :gads)
		omega, wi = mri_trajectory_gads(N, fov ; kwargs...)
	end

	if (ktype == :radial)
@show N, fov
@show N[1:2], fov[1:2]
		omega, wi = mri_trajectory_radial(N[1:2], fov[1:2] ; kwargs...)

		if (length(N) == 3) # 'stack of radials'
			omega = mri_trajectory_stack(omega, N[3])
			wi = repeat(wi, N[3], 1)
			wi = wi[:] / fov[3]
		end
	end

	if (ktype == :rosette3)
		omega, wi = mri_trajectory_rosette3(N, fov ; kwargs...)
	end

	# echo-planar with sinusoid:
	if (ktype == :epi_sin)
		tmp = ( ; oversample=1) -> oversample
		oversample = tmp( ; kwargs...)
		Npt = round(Int, oversample * prod(N))
		t = (0:(Npt-1))/Npt
		omega = [π*sin.(2π*t*N[2]/2) t*2π .- π]
	end

	# bad spiral:
	if (ktype == :spiral0)
		Nspiral = 400
		omega = LinRange(0, 10*2π, Nspiral)
		omega = π * [cos.(omega) sin.(omega)] .* [omega omega]/maximum(omega)
		wi = abs.(omega[:,1] + 1im * omega[:,2]) # simple |r| weighting
	end

	# crude spiral:
	if (ktype == :spiral1)
		Nspiral = round(Int, prod(N) * π/4)
		omega = LinRange(0, N[1]*2π, Nspiral)
		max_omega = maximum(omega)
		omega = π * [cos.(omega) sin.(omega)].* [omega omega]/max_omega
		wi = abs.(omega[:,1] + 1im * omega[:,2]) # simple |r| weighting
	end

	# 3T spiral:
	if (ktype == :spiral3)
		kspace, omega = mri_kspace_spiral(N = maximum(N[1:2]), fov = maximum(fov[1:2]))

		if length(N) == 3 #stack of these spirals
			omega = mri_trajectory_stack(omega, N[3])
		end

		wi = abs.(omega[:,1] + 1im * omega[:,2]) # simple |r| weighting
		wi = wi / (fov[1]*fov[2]) # approximate scaling
		if length(N) == 3
			length(fov) != 3 && throw("fov not 3D")
			wi = wi / fov[3] # cartesian-style weighting in z
		end
	end

	# random
	if (ktype == :random)
		seed!(0)
		omega = (rand(prod(N), 2) .- 0.5) * 2π
	end

	# 2D half Cartesian + 8 rows
	if (ktype == :half8)
		extra = 8
		o1 = ((0:(N[1]-1))/N[1] .- 0.5)*2π
		o2 = (-N[2]/2:extra)/N[2] * 2π
		o1,o2 = ndgrid(o1, o2)
		omega = [o1[:] o2[:]]
	end

	# 2D FT, undersampled in "y" (phase encode) direction
	if (ktype == :cart_y_2)
		under = 2
		o1 = ((0:(N[1]/1-1))/(N[1]/1) .- 0.5) * 2π
		o2 = ((0:(N[2]/2-1))/(N[2]/under) .- 0.5) * 2π
		o1,o2 = ndgrid(o1, o2)
		omega = [o1[:] o2[:]]
	end

	# convert to physical units
	kspace = zeros(Float32, size(omega))
	for id = 1:length(N)
		dx = fov[id] / N[id]
		kspace[:,id] = omega[:,id] / (2π) / dx
	end

	# wi = ir_mri_density_comp(kspace, arg_wi{:}) # todo

	return kspace, omega, wi
end

# mri_trajectory_stack()
# make 3D "stack of ..." trajectory from a 2D trajectory
function mri_trajectory_stack(omega2, N3::Int)
	o3 = ((0:(N3-1))/N3 .- 0.5)*2π
	o3 = repeat(o3, 1, size(omega2,1)) # [N3,N12]
	o3 = o3'[:] # [N12*N3]
	omega = repeat(omega2, N3) # [N12*N3,2]
	omega = [omega o3] # [N12*N3,3]
	return omega
end

# mri_trajectory_epi_under()
# EPI, with optional under-sampling
function mri_trajectory_epi_under(N, fov ;
		under::Int = 2, # default is every other phase encode
		samp::AbstractVector{<:Bool} = rem.(0:(N[2]-1), under) .== 0,
	)
	nx = N[1]
	ny = N[2]
	omx = (-nx/2:nx/2-1) / nx * 2π # [-π,π) in x
	omy = (-ny/2:ny/2-1) / ny * 2π # [-π,π) in y
	omega = zeros(0,2)
	for iy in findall(samp)
		omega = [omega; [omx omy[iy]*ones(nx)]]
		omx = reverse(omx, dims = 1) # non-flyback
	end
	wi = ones(size(omega,1)) / prod(fov)
	return omega, wi
end

"""
    omega, wi = mri_trajectory_gads(N, fov ; ...)
emulate 2D golden angle radial sampling with data sharing

option:
`Nro` # of samples in each readout/spoke
`shift` shift along read-out due to gradient delays (stress)
`kmax_frac` fractions of maximum krad (0.5) for rings
`under` under-sampling factor for each annulus
"""
function mri_trajectory_gads(N, fov ;
		Nro::Int = maximum(N),
		delta_ro::Real = 1/Nro,
		shift::Real = -0.75,
		kmax_frac::Array = [0.20, 0.35, 0.501],
		nspoke::Array = floor.(Int, π * kmax_frac * Nro),
		under::Array = [1, 1, 0.6],
		start::AbstractArray = (0:2) * π/4
	)

	nspoke = floor.(Int, nspoke .* under)
	if true # make fibonacci for more uniform coverage
		phi = (1 + sqrt(5)) / 2
		n = round.(log.(nspoke * sqrt(5) .+ 0.5) / log(phi))
		nspoke = phi.^n / sqrt(5) .+ 0.5
		nspoke = floor.(Int, nspoke)
	end
	nring = length(nspoke)
	kmax_frac = [0; kmax_frac]
	omega = zeros(0,2)
	for ir=1:nring
		#todo: not sure here
		kspace = ir_mri_kspace_ga_radial(Nspoke = nspoke[ir],
			Nro = Nro, delta_ro = delta_ro, shift = shift, start = start[ir])

		kspace = reshape(kspace, :, 2)
		krad = sqrt.(sum(kspace.^2, dims = 2)[:])
		good = (kmax_frac[ir] .<= krad) .& (krad .< kmax_frac[ir+1])
		kspace = kspace[good, :]
		omega = [omega; 2π*kspace]
	end
	wi = [] # no default DCF
	return omega, wi
end

"""
    mri_trajectory_radial()

option:
* `na_nr` default ensures proper sampling at edge of k-space
* `na` angular spokes; default: na_nr * nr
* `nr` radial samples per spoke
* `ir` default: `0:nr`

todo: generalize to 3D using barger:02:trc
"""
function mri_trajectory_radial(N::Dims, fov ;
		nr::Int = round(Int, maximum(N)/2),
		ir::AbstractArray{Int} = 0:nr,
		na_nr::Real = 2π,
		na::Int = 4*ceil(Int, na_nr * nr/4), # mult of 4
	)

	om = ir/nr * π
	ang = (0:na-1)/na * 2π
	om, ang = ndgrid(om, ang) # [nr+1, na]
	omega = [(om .* cos.(ang))[:] (om .* sin.(ang))[:]]

	# density compensation factors based on "analytical" voronoi
	fov[2] != fov[1] && @warn("fov not square")
	du = 1/fov[1] # assume this radial sample spacing
	wi = π * du^2 / na * 2 * ir # see lauzon:96:eop, joseph:98:sei
	wi = collect(wi)
	wi[ir .== 0] .= π * (du/2)^2 / na

	wi = repeat(wi, na)[:]
	return omega, wi
end

"""
    mri_trajectory_rosette3(N, fov ; ...)
3d rosette, with default parameters from bucholz:08:miw

options:
omax: maximum omega
nt : time samples (65.536 ms for 4 usec dt)
dt : time sample spacing (4 usec)
ti : time samples
"""
function mri_trajectory_rosette3(N, fov ;
		f1::Real = 211,
		f2::Real = 117.13,
		f3::Real = length(N) == 3 ? 73.65 : 0,
		omax::Real = π,
		nshot::Int = 32, # todo: shots vs arms
		nt::Int = length(N) == 3 ? 16385 : 2^9, # time samples (65.536 ms for 4 usec dt)
		dt::Real = 4e-6, # time sample spacing (4 usec)
		ti::AbstractVector = (0:nt-1) * dt, # time samples
	)

	tmp = 2π * ti
	p1 = f1 * tmp
	p2 = f2 * tmp
	p3 = f3 * tmp
	kx = omax * sin.(p1) .* cos.(p2) .* cos.(p3)
	ky = omax * sin.(p1) .* sin.(p2) .* cos.(p3)
	kz = omax * sin.(p1) .* sin.(p3)
	omega = [kx ky kz]
	for is = 1:nshot-1 # n-shot, rotate kx,ky by 2π / N
		ang = is * 2π / nshot
		c = cos(ang)
		s = sin(ang)
		ox = c * kx + s * ky
		oy = -s * kx + c * ky
		omega = [omega; [ox oy kz]]
	end
	wi = omax^3 * abs.(sin.(p1).^2 .* cos.(p1) .* cos.(p3)) # bucholz:08:miw
	if length(N) == 2
		omega = omega[:,1:2]
	end
	return omega, wi
end

"""
    mri_trajectory(:test)
self test
"""
function mri_trajectory(test::Symbol)
	test != :test && throw(DomainError(test, "Not valid"))
	ig = image_geom_mri(nx = 2^5, ny = 2^5-4, fov = 250) # 250 mm FOV
	N = ig.dim

	ktype3 = [:cartesian, :spiral3, :radial] # 3D
	for ktype in ktype3 # 3D tests
		@show ktype
		kspace, omega, wi = mri_trajectory( ; ktype=ktype, N=(N...,4))
		@test size(kspace,2) == size(omega,2) == 3
		@test size(kspace,1) == size(omega,1)
		@test length(wi) ∈ (1, size(kspace,1))
		scatter(omega[:,1]/π, omega[:,2]/π, omega[:,3]/π, label="$ktype")
	end
	prompt()

	ktypes = [:spiral3, :spiral1, :spiral0, :radial, :gads, :rosette3,
		:epi_sin, :epi_under, :cartesian, :half8, :cart_y_2, :random,]
	# arg_wi = [:voronoi]
	plots = Any[]
	for ktype in ktypes
		@show ktype
		args =
			(ktype == :epi_sin) ? Dict(:oversample => 1/1,) :
			(ktype == :radial) ? Dict(:nr => 10, :na_nr => π/2) :
				Dict{Symbol,Nothing}()

		kspace, omega, wi = mri_trajectory( ;
			ktype=ktype, N=N, fov = ig.fovs, args...)
		p = scatter(omega[:,1]/π, omega[:,2]/π, aspect_ratio=1,
			xlim = [-1,1]*1.1, xtick = -1:1, xlabel = "omega1/π",
			ylim = [-1,1]*1.1, ytick = -1:1, ylabel = "omega2/π",
			title = string(ktype), label="")
		push!(plots, p)
	end
	plot(plots...); gui()

	return true
end


text:
using Test, Unitful, PhysicalConstants.CODATA2018
using PlasmaParameters

include("unit/runtests.jl")


text:
#=

    This is a port of the SPPL model ClinicalTrial.pynb

    Note: There is no data in this SPPL model, so I take
    the data from clinical_trial.jl (which is a port of the R2 model ClinicalTrial.cs)

    with the datafiles
    - clinical_trial_control.csv  (size 1000)
    - clinical_trial_treatment.csv (size 1000)

    The SPPL model takes 20 observations and give this exact probability:
       p(isEffective): 0.3658327735005218

    This model:
    Summary Statistics
        parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
            Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

        isEffective    0.3603    0.4801     0.0048    0.0116   1739.9975    1.0002      178.8465
        probControl    9.5703    5.1312     0.0513    0.1206   1744.4818    1.0019      179.3074
      probTreatment   10.3307    5.1155     0.0512    0.1228   1789.7505    1.0001      183.9604
            probAll    9.7317    3.9333     0.0393    0.0986   1621.2317    1.0014      166.6391
           pControl    0.4785    0.2566     0.0026    0.0060   1744.4818    1.0019      179.3074
         pTreatment    0.5165    0.2558     0.0026    0.0061   1789.7505    1.0001      183.9604

    Cf clinical_trial.jl (which use Bernoilli on the probabilities instead of DiscreteUniform.)

=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV
include("jl_utils.jl")

@model function clinical_trial(control,treatment)
    n = length(control)

    isEffective   ~ flip(0.5)
    probControl   ~ DiscreteUniform(0,n)
    probTreatment ~ DiscreteUniform(0,n)
    probAll       ~ DiscreteUniform(0,n)

    if isEffective
        control   ~ filldist(Bernoulli(probControl/n),n)
        treatment ~ filldist(Bernoulli(probTreatment/n),n)
    else
        control   ~ filldist(Bernoulli(probAll/n),n)
        treatment ~ filldist(Bernoulli(probAll/n),n)

    end

    # The probabilities
    pControl ~ Dirac(probControl / n)
    pTreatment ~ Dirac(probTreatment / n)
end

control = parse.(Int,split(readline("clinical_trial_control.csv"),","))
treatment = parse.(Int,split(readline("clinical_trial_treatment.csv"),","))
println("len(control):$(length(control)) len(treatment):$(length(treatment)) mean(control):$(mean(control)) mean(treatment): $(mean(treatment))")
# println("treatment:$treatment")

n = 20
model = clinical_trial(control[1:n], treatment[1:n])

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, HMC(0.1,6), 1_000)
# chns = sample(model, NUTS(), 1_000)

# display(chns[[:isEffective,:probControl,:probTreatment,:probAll,:pControl,:pTreatment]])
display(chns)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

code = create_basis(
    :Hex8,
    "8 node linear hexahedral element",
    (
     (-1.0, -1.0, -1.0), # N1
     ( 1.0, -1.0, -1.0), # N2
     ( 1.0,  1.0, -1.0), # N3
     (-1.0,  1.0, -1.0), # N4
     (-1.0, -1.0,  1.0), # N5
     ( 1.0, -1.0,  1.0), # N6
     ( 1.0,  1.0,  1.0), # N7
     (-1.0,  1.0,  1.0), # N8
    ),
    :(1 + u + v + w + u*v + v*w + w*u + u*v*w),
   )
eval(code)

code = create_basis(
    :Hex20,
    "20 node biquadratic hexahedral element",
    (
     (-1.0, -1.0, -1.0), # N1
     ( 1.0, -1.0, -1.0), # N2
     ( 1.0,  1.0, -1.0), # N3
     (-1.0,  1.0, -1.0), # N4
     (-1.0, -1.0,  1.0), # N5
     ( 1.0, -1.0,  1.0), # N6
     ( 1.0,  1.0,  1.0), # N7
     (-1.0,  1.0,  1.0), # N8
     ( 0.0, -1.0, -1.0), # N9
     ( 1.0,  0.0, -1.0), # N10
     ( 0.0,  1.0, -1.0), # N11
     (-1.0,  0.0, -1.0), # N12
     (-1.0, -1.0,  0.0), # N13
     ( 1.0, -1.0,  0.0), # N14
     ( 1.0,  1.0,  0.0), # N15
     (-1.0,  1.0,  0.0), # N16
     ( 0.0, -1.0,  1.0), # N17
     ( 1.0,  0.0,  1.0), # N18
     ( 0.0,  1.0,  1.0), # N19
     (-1.0,  0.0,  1.0), # N20
    ),
    :(1 + u + v + w + u*v + v*w + u*w + u*v*w + u^2 + v^2 + w^2 + u^2*v + u*v^2 + v^2*w + v*w^2 + u*w^2 + u^2*w + u^2*v*w + u*v^2*w + u*v*w^2),
   )
eval(code)

code = create_basis(
    :Hex27,
    "27 node quadratic hexahedral element",
    (
     (-1.0, -1.0, -1.0), # N1
     ( 1.0, -1.0, -1.0), # N2
     ( 1.0,  1.0, -1.0), # N3
     (-1.0,  1.0, -1.0), # N4
     (-1.0, -1.0,  1.0), # N5
     ( 1.0, -1.0,  1.0), # N6
     ( 1.0,  1.0,  1.0), # N7
     (-1.0,  1.0,  1.0), # N8
     ( 0.0, -1.0, -1.0), # N9
     ( 1.0,  0.0, -1.0), # N10
     ( 0.0,  1.0, -1.0), # N11
     (-1.0,  0.0, -1.0), # N12
     (-1.0, -1.0,  0.0), # N13
     ( 1.0, -1.0,  0.0), # N14
     ( 1.0,  1.0,  0.0), # N15
     (-1.0,  1.0,  0.0), # N16
     ( 0.0, -1.0,  1.0), # N17
     ( 1.0,  0.0,  1.0), # N18
     ( 0.0,  1.0,  1.0), # N19
     (-1.0,  0.0,  1.0), # N20
     ( 0.0,  0.0, -1.0), # N21
     ( 0.0, -1.0,  0.0), # N22
     ( 1.0,  0.0,  0.0), # N23
     ( 0.0,  1.0,  0.0), # N24
     (-1.0,  0.0,  0.0), # N25
     ( 0.0,  0.0,  1.0), # N26
     ( 0.0,  0.0,  0.0), # N27
    ),
    :(1 + u + v + w + u*v + v*w + u*w + u*v*w + u^2 + v^2 + w^2 + u^2*v + u*v^2 + v^2*w + v*w^2 + u*w^2 + u^2*w + u^2*v*w + u*v^2*w + u*v*w^2 + u^2*v^2 + v^2*w^2 + u^2*w^2 + u^2*v^2*w + u*v^2*w^2 + u^2*v*w^2 + u^2*v^2*w^2),
   )
eval(code)


text:
abstract type SDIRKMutableCache <: OrdinaryDiffEqMutableCache end

@cache mutable struct ImplicitEulerCache{uType,rateType,uNoUnitsType,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  fsalfirst::rateType
  atmp::uNoUnitsType
  nlsolver::N
end

function alg_cache(alg::ImplicitEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  γ, c = 1, 1
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  fsalfirst = zero(rate_prototype)

  atmp = similar(u,uEltypeNoUnits)

  ImplicitEulerCache(u,uprev,uprev2,fsalfirst,atmp,nlsolver)
end

mutable struct ImplicitEulerConstantCache{N} <: OrdinaryDiffEqConstantCache
  nlsolver::N
end

function alg_cache(alg::ImplicitEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  γ, c = 1, 1
  J, W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  ImplicitEulerConstantCache(nlsolver)
end

mutable struct ImplicitMidpointConstantCache{N} <: OrdinaryDiffEqConstantCache
  nlsolver::N
end

function alg_cache(alg::ImplicitMidpoint,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  γ, c = 1//2, 1//2
  J, W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  ImplicitMidpointConstantCache(nlsolver)
end

@cache mutable struct ImplicitMidpointCache{uType,rateType,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  nlsolver::N
end

function alg_cache(alg::ImplicitMidpoint,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  γ, c = 1//2, 1//2
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  fsalfirst = zero(rate_prototype)
  ImplicitMidpointCache(u,uprev,fsalfirst,nlsolver)
end

mutable struct TrapezoidConstantCache{uType,tType,N} <: OrdinaryDiffEqConstantCache
  uprev3::uType
  tprev2::tType
  nlsolver::N
end

function alg_cache(alg::Trapezoid,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  γ, c = 1//2, 1
  J, W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)

  uprev3 = u
  tprev2 = t

  TrapezoidConstantCache(uprev3,tprev2,nlsolver)
end

@cache mutable struct TrapezoidCache{uType,rateType,uNoUnitsType,tType,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  fsalfirst::rateType
  atmp::uNoUnitsType
  uprev3::uType
  tprev2::tType
  nlsolver::N
end

function alg_cache(alg::Trapezoid,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  γ, c = 1//2, 1
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  fsalfirst = zero(rate_prototype)

  uprev3 = zero(u)
  tprev2 = t
  atmp = similar(u,uEltypeNoUnits)

  TrapezoidCache(u,uprev,uprev2,fsalfirst,atmp,uprev3,tprev2,nlsolver)
end

mutable struct TRBDF2ConstantCache{Tab,N} <: OrdinaryDiffEqConstantCache
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::TRBDF2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tab = TRBDF2Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.d, tab.γ
  J, W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  TRBDF2ConstantCache(nlsolver,tab)
end

@cache mutable struct TRBDF2Cache{uType,rateType,uNoUnitsType,Tab,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  zprev::uType
  zᵧ::uType
  atmp::uNoUnitsType
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::TRBDF2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = TRBDF2Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.d, tab.γ
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  fsalfirst = zero(rate_prototype)

  atmp = similar(u,uEltypeNoUnits); zprev = similar(u); zᵧ = similar(u)

  TRBDF2Cache(u,uprev,fsalfirst,zprev,zᵧ,atmp,nlsolver,tab)
end

mutable struct SDIRK2ConstantCache{N} <: OrdinaryDiffEqConstantCache
  nlsolver::N
end

function alg_cache(alg::SDIRK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  γ, c = 1, 1
  J, W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  SDIRK2ConstantCache(nlsolver)
end

@cache mutable struct SDIRK2Cache{uType,rateType,uNoUnitsType,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  z₁::uType
  z₂::uType
  atmp::uNoUnitsType
  nlsolver::N
end

function alg_cache(alg::SDIRK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  γ, c = 1, 1
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  fsalfirst = zero(rate_prototype)

  z₁ = similar(u); z₂ = nlsolver.z
  atmp = similar(u,uEltypeNoUnits)

  SDIRK2Cache(u,uprev,fsalfirst,z₁,z₂,atmp,nlsolver)
end

mutable struct SSPSDIRK2ConstantCache{N} <: OrdinaryDiffEqConstantCache
  nlsolver::N
end

function alg_cache(alg::SSPSDIRK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  γ, c = 1//4, 1//1
  J, W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  SSPSDIRK2ConstantCache(nlsolver)
end

@cache mutable struct SSPSDIRK2Cache{uType,rateType,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  z₁::uType
  z₂::uType
  nlsolver::N
end

function alg_cache(alg::SSPSDIRK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  γ, c = 1//4, 1//1
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  fsalfirst = zero(rate_prototype)

  z₁ = similar(u); z₂ = nlsolver.z
  atmp = similar(u,uEltypeNoUnits)

  SSPSDIRK2Cache(u,uprev,fsalfirst,z₁,z₂,nlsolver)
end

mutable struct Kvaerno3ConstantCache{Tab,N} <: OrdinaryDiffEqConstantCache
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::Kvaerno3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tab = Kvaerno3Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ, 2tab.γ
  J, W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  Kvaerno3ConstantCache(nlsolver,tab)
end

@cache mutable struct Kvaerno3Cache{uType,rateType,uNoUnitsType,Tab,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  z₁::uType
  z₂::uType
  z₃::uType
  z₄::uType
  atmp::uNoUnitsType
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::Kvaerno3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = Kvaerno3Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ, 2tab.γ
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  fsalfirst = zero(rate_prototype)

  z₁ = similar(u); z₂ = similar(u); z₃ = similar(u); z₄ = nlsolver.z
  atmp = similar(u,uEltypeNoUnits)

  Kvaerno3Cache(u,uprev,fsalfirst,z₁,z₂,z₃,z₄,atmp,nlsolver,tab)
end

mutable struct Cash4ConstantCache{N,Tab} <: OrdinaryDiffEqConstantCache
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::Cash4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tab = Cash4Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ,tab.γ
  J, W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  Cash4ConstantCache(nlsolver,tab)
end

@cache mutable struct Cash4Cache{uType,rateType,uNoUnitsType,N,Tab} <: SDIRKMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  z₁::uType
  z₂::uType
  z₃::uType
  z₄::uType
  z₅::uType
  atmp::uNoUnitsType
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::Cash4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = Cash4Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ,tab.γ
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  fsalfirst = zero(rate_prototype)

  z₁ = similar(u); z₂ = similar(u); z₃ = similar(u); z₄ = similar(u); z₅ = nlsolver.z
  atmp = similar(u,uEltypeNoUnits)

  Cash4Cache(u,uprev,fsalfirst,z₁,z₂,z₃,z₄,z₅,atmp,nlsolver,tab)
end

mutable struct Hairer4ConstantCache{N,Tab} <: OrdinaryDiffEqConstantCache
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::Union{Hairer4,Hairer42},u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  if alg isa Hairer4
    tab = Hairer4Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  else
    tab = Hairer42Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  end
  γ, c = tab.γ, tab.γ
  J, W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  Hairer4ConstantCache(nlsolver,tab)
end

@cache mutable struct Hairer4Cache{uType,rateType,uNoUnitsType,Tab,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  z₁::uType
  z₂::uType
  z₃::uType
  z₄::uType
  z₅::uType
  atmp::uNoUnitsType
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::Union{Hairer4,Hairer42},u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  if alg isa Hairer4
    tab = Hairer4Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  else # Hairer42
    tab = Hairer42Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  end
  γ, c = tab.γ, tab.γ
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  fsalfirst = zero(rate_prototype)

  z₁ = similar(u); z₂ = similar(u); z₃ = similar(u); z₄ = similar(u); z₅ = nlsolver.z
  atmp = similar(u,uEltypeNoUnits)

  Hairer4Cache(u,uprev,fsalfirst,z₁,z₂,z₃,z₄,z₅,atmp,nlsolver,tab)
end

@cache mutable struct ESDIRK54I8L2SACache{uType,rateType,uNoUnitsType,Tab,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  z₁::uType; z₂::uType; z₃::uType; z₄::uType; z₅::uType; z₆::uType; z₇::uType; z₈::uType
  atmp::uNoUnitsType
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::ESDIRK54I8L2SA,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = ESDIRK54I8L2SATableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ, tab.γ
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  fsalfirst = zero(rate_prototype)

  z₁ = zero(u); z₂ = zero(u); z₃ = zero(u); z₄ = zero(u)
  z₅ = zero(u); z₆ = zero(u); z₇ = zero(u); z₈ = nlsolver.z
  atmp = similar(u,uEltypeNoUnits)

  ESDIRK54I8L2SACache(u,uprev,fsalfirst,z₁,z₂,z₃,z₄,z₅,z₆,z₇,z₈,atmp,nlsolver,tab)
end

mutable struct ESDIRK54I8L2SAConstantCache{N,Tab} <: OrdinaryDiffEqConstantCache
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::ESDIRK54I8L2SA,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tab = ESDIRK54I8L2SATableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ,tab.γ
  J, W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  ESDIRK54I8L2SAConstantCache(nlsolver,tab)
end


text:
#!/usr/bin/env julia

using DipoleChain

using ArgParse

s = ArgParseSettings()
s.autofix_names = true
@add_arg_table s begin
    "-R"
        metavar = "R"
        help = "separation distance"
        arg_type = Float64
        required = true
    "-N"
        metavar = "N"
        help = "number of rotors"
        arg_type = Int
        required = true
    "--l-max"
        metavar = "L"
        help = "local basis truncation"
        arg_type = Int
        required = true
    "--l-total_max"
        metavar = "L"
        help = "many-body basis truncation (default: no truncation)"
        arg_type = Int
    "--beta"
        metavar = "B"
        help = "reciprocal temperature"
        arg_type = Float64
    "--tau"
        metavar = "T"
        help = "reciprocal temperature per link"
        arg_type = Float64
    "-P"
        metavar = "P"
        help = "number of beads"
        arg_type = Int
        required = true
    "--pigs"
        help = "use PIGS"
        action = :store_true
    "--A-start"
        metavar = "jA"
        help = "first rotor of subsystem A (PIGS only)"
        arg_type = Int
        default = 1
    "--A-size"
        metavar = "NA"
        help = "length of subsystem A (PIGS only)"
        arg_type = Int
        default = 1
end
c = parse_args(ARGS, s, as_symbols=true)

R = c[:R]
N = c[:N]
l_max = c[:l_max]
l_total_max = c[:l_total_max]
if isnothing(l_total_max)
    l_total_max = l_max * N
end
beta = c[:beta]
tau = c[:tau]
if isnothing(beta) && isnothing(tau) || !isnothing(beta) && !isnothing(tau)
    println("Exactly one of --beta or --tau must be set.")
    exit()
end
P = c[:P]
pigs = c[:pigs]
A_start = c[:A_start]
A_size = c[:A_size]

if pigs
    if P % 2 == 0
        println("P must be odd for PIGS.")
        exit()
    end
    num_links = P - 1
else
    num_links = P
end

if isnothing(beta)
    beta = tau * num_links
elseif isnothing(tau)
    tau = beta / num_links
end

println("[ ] Constructing basis.")
@time if pigs
    # Always use the ground state symmetry block for PIGS.
    sym_lp = 0
    sym_m = 0
    basis = SingleBlockBasis(N, l_max, l_total_max, sym_lp, sym_m)
else
    basis = MultiBlockBasis(N, l_max, l_total_max)
end
println("[+] Constructed basis: $(basis.size).")

println("[ ] Constructing matrix elements.")
@time K = diag_rot(basis)
@time V = dense_pot(basis)
@time scale_dense_pot(V, 1.0/R^3)
println("[+] Constructed matrix elements.")

println("[ ] Constructing path.")
@time if pigs
    # Trial function is free rotor ground state.
    trial_idx = basis.block.lookup[[0 for _ in 1:2N]]

    half_path = symmetric_trotter_path(div(num_links, 2), tau, K, V)
    path = half_path * half_path
    wf = wf_pigs(half_path, trial_idx)
else
    path = symmetric_trotter_path(num_links, tau, K, V)
end
println("[+] Constructed path.")

println("[ ] Calculating energies.")
@time if pigs
    E0, (avgK, avgV) = energy_mixed(path, trial_idx, K, V)
    # For a spherically symmetric trial function, the kinetic contribution to
    # the mixed estimator should vanish.
    avgK == 0.0 || @warn("avgK == $(avgK)")
else
    Z, E, (avgK, avgV) = energy(path, K, V)
end
println("[+] Calculated energies.")

if pigs
    # Ground state energy.
    println_result("E0 = $(E0)")
else
    # Partition function.
    println_result("Z = $(Z)")
    # Rotational energy.
    println_result("K = $(avgK)")
    # Potential energy.
    println_result("V = $(avgV)")
    # Total energy.
    println_result("E = $(E)")
end

if pigs
    println("[ ] Constructing subsystem bases.")
    @time basis_A = MultiBlockBasis(A_size, l_max, l_total_max)
    @time basis_B = MultiBlockBasis(N-A_size, l_max, l_total_max)
    println("[+] Constructed subsystem bases: $(basis_A.size), $(basis_B.size).")

    println("[ ] Calculating reduced eigenvalues.")
    @time eigvals = reduced_eigenvalues(basis, basis_A, basis_B, A_start, wf)
    println("[+] Calculated reduced eigenvalues.")

    # Von Neumann entropy.
    println_result("SvN = $(S_vn(eigvals))")
    # Order-2 Rényi entropy.
    println_result("S2 = $(S_renyi(eigvals, 2))")
end


text:
""" An embedding holding only its points and no information about the
construction of the embedding. """
mutable struct SimpleEmbedding{D, T} <: AbstractEmbedding{D, T}
    points::AbstractArray{T, 2}
end

export SimpleEmbedding


text:
module TestDiversityMeasure
using Test
using LinearAlgebra

using Diversity
using Diversity.ShortNames
using DataFrames

pop = [3, 3, 4]
pop = pop / sum(pop)
oc = Onecommunity()
sim = [1.0 0 0; 1.0 1.0 0.0; 1.0 1.0 1.0]
ms = GeneralTypes(sim)
meta = Metacommunity(pop, ms, oc)
meta1 = Metacommunity(pop)

ab3 = [1.0 2; 3 0; 0 4]'
ab3 = ab3 / sum(ab3)
sc = Subcommunities(3)
sp = Species(2)
meta2 = Metacommunity(ab3, sp, sc)
nab = NormalisedAlpha(meta2)

@testset "Diversity measures" begin
    diversities = [RawAlpha, NormalisedAlpha, RawBeta, NormalisedBeta,
                   RawRho, NormalisedRho, Gamma]
    shortds = [α, ᾱ, β, β̄, ρ, ρ̄, Γ]
    chars = ["α", "ᾱ", "β", "β̄", "ρ", "ρ̄", "γ"]
    asciis = ["RawAlpha", "NormalisedAlpha",
              "RawBeta", "NormalisedBeta",
              "RawRho", "NormalisedRho", "Gamma"]
    fulls = ["raw alpha diversity", "normalised alpha diversity",
             "distinctiveness", "effective number of subcommunities",
             "redundancy", "representativeness", "gamma diversity"]
    for i in 1:length(diversities)
        @test diversities[i] == shortds[i]
        @test getName(diversities[i](meta)) == chars[i]
        @test getASCIIName(diversities[i](meta2)) == asciis[i]
        @test getFullName(diversities[i](meta1)) == fulls[i]
    end
end

numbers = [1., 2, 4, 8, 16];
numspecies = 100;
fragments = rand(numspecies);
weights = rand(numspecies);
weights /= sum(weights);
Z1 = ones(typeof(weights[1]), (length(weights), length(weights)));
numcommunities = 8;
manyweights = rand(numspecies, numcommunities);
manyweights *= Diagonal(reshape(mapslices(v -> 1. / sum(v), manyweights;
                                          dims=1), numcommunities));

@testset "inddiv / subdiv / metadiv" begin
    @test individualDiversity(nab, 0)[!,:diversity] ≈ inddiv(nab, 0)[!,:diversity]
    @test individualDiversity(nab)(1)[!,:diversity] ≈ inddiv(nab, 1)[!,:diversity]
    idnab = inddiv(nab, [2, 3])
    @test idnab[isapprox.(collect(idnab[!,:q]), 2), :diversity] ≈ inddiv(nab, 2)[!,:diversity]
    allid = inddiv(meta2, Inf)
    @test allid[allid[!,:measure] .== "RawAlpha", :diversity] ≈ inddiv(RawAlpha(meta2), Inf)[!,:diversity]

    @test subcommunityDiversity(nab, Inf)[!,:diversity] ≈ subdiv(nab, Inf)[!,:diversity]
    sdnab = subdiv(nab, [4, 5])
    @test sdnab[isapprox.(collect(sdnab[!,:q]), 4), :diversity] ≈ subdiv(nab, 4)[!,:diversity]
    allsd = subdiv(meta2, Inf)
    @test allsd[allsd[!,:measure] .== "NormalisedAlpha", :diversity] ≈ subdiv(nab, Inf)[!,:diversity]
    allmd = subdiv(meta1, Inf)
    @test allmd[allmd[!,:measure] .== "NormalisedAlpha", :diversity] ≈ metadiv(ᾱ(meta1), Inf)[!,:diversity]

    scg = subcommunityDiversity(Gamma(meta2))
    @test scg(1)[!,:diversity] ≈ scg(1.0)[!,:diversity]

    communities = rand(numspecies, numcommunities);
    communities /= sum(communities);
    @test subdiv(NormalisedAlpha(Metacommunity(communities)), 0)[!,:diversity] ≈ numspecies * ones(size(communities, 2))
    @test subdiv(NormalisedAlpha(Metacommunity(communities)), [0])[!,:diversity] ≈ numspecies * ones(size(communities, 2))
    qs = [0, 1, 2, Inf]
    sna = subdiv(NormalisedAlpha(Metacommunity(communities, Z1)), qs)
    @test nrow(sna) == length(qs) * size(communities, 2)
    for q in qs
        @test sna[isapprox.(sna[!,:q], q), :diversity] ≈ ones(size(communities, 2))
    end
    @test subdiv(RawAlpha(Metacommunity(communities)), 0)[!,:diversity] ≈ numspecies * vec(mapslices(v -> 1. / sum(v), communities; dims=1))

    even = ones((numspecies, numcommunities)) / (numspecies * numcommunities);
    qs = [0, 1, 2, 3, 4, 5, 6, Inf];
    @test metadiv(NormalisedAlpha(Metacommunity(even)), qs)[!,:diversity] ≈ numspecies * ones(length(qs))
    @test metadiv(RawAlpha(Metacommunity(even)), qs)[!,:diversity] ≈ numspecies * numcommunities * ones(length(qs))
    md2 = metadiv(meta2, Inf)
    @test md2[md2[!,:measure] .== "Gamma", :diversity] ≈ metadiv(Gamma(meta2), Inf)[!,:diversity]

    probs = reshape(mapslices(sum, communities; dims=2),
                    size(communities, 1));
    @test metadiv(Gamma(Metacommunity(communities)), qs)[!,:diversity] ≈ qD(probs, qs)
    @test metadiv(Gamma(Metacommunity(communities, Z1)), qs)[!,:diversity] ≈ qDZ(probs, qs, Z1)

    Z = rand(numspecies, numspecies)
    @test metadiv(Gamma(Metacommunity(communities, Z)), qs)[!,:diversity] ≈ qDZ(probs, qs, Z)

    colweights = rand(numcommunities)
    colweights /= sum(colweights)
    allthesame = probs * colweights'
    @test metadiv(RawBeta(Metacommunity(allthesame, Z)), qs)[!,:diversity] ≈ 1.0 ./ qD(colweights, 2 .- qs)
    @test metadiv(NormalisedBeta(Metacommunity(allthesame, Z)), qs)[!,:diversity] ≈ ones(length(qs))
    @test metadiv(NormalisedRho(Metacommunity(allthesame, Z)), qs)[!,:diversity] ≈ ones(length(qs))
    @test metadiv(RawRho(Metacommunity(allthesame, Z)), qs)[!,:diversity] ≈ qD(colweights, qs)

    communitylist = rand(1:numcommunities, numspecies)
    distinct = zeros(Float64, (numspecies, numcommunities))
    for i in 1:numspecies
        distinct[i, communitylist[i]] = weights[i]
    end

    @test metadiv(RawRho(Metacommunity(distinct)), qs)[!,:diversity] ≈ ones(length(qs))
    subnr = subdiv(NormalisedRho(Metacommunity(distinct)), qs)
    for q in qs
        @test subnr[isapprox.(subnr[!,:q], q), :diversity] ≈
              vec(sum(distinct; dims=1))
    end
    @test metadiv(NormalisedBeta(Metacommunity(distinct)), qs)[!,:diversity] ≈
        qD(reshape(sum(distinct; dims=1), numcommunities), qs)
    @test metadiv(RawBeta(Metacommunity(distinct)), qs)[!,:diversity] ≈ ones(length(qs))

    # many (unexported!) diversity levels not yet implemented
    @test_throws ErrorException Diversity.communityDiversity(nab)
end

end


text:
# Stocks/Prices

import Dates: Date
import DataFrames: DataFrame, Not

# std_cols = "date,label,open,high,low,close,volume,uOpen,uHigh,uLow,uClose,uVolume"
const std_cols = "date,open,high,low,close,volume"

"""
    history(
        conn::Connection, symbol::String;
        range::String = "1m", date::String = "",
        filter = std_cols;
        chartCloseOnly::Bool = false
    )

Gets adjusted and unadjusted historical data
for up to 15 years. Useful for building charts.

This endpoint is carried over from the IEX 1.0 API,
and we know it has become difficult to use. We are
working on new endpoints to simplify time series queries.
"""
function history(
    conn::Connection,
    symbol::String;
    range::String = "1m",
    date::String = "",
    filter = std_cols,
    chartCloseOnly::Bool = false
)::DataFrame

    path_params = [:symbol, :range, :date]
    d = get_stocks(
        conn, "chart", path_params;
        symbol = symbol,
        range = range,
        date = date,
        filter = filter,
        chartCloseOnly = chartCloseOnly
    )
    colnames = split("symbol," * filter, ",")
    if length(d) > 0
        df = vcat(map(x -> DataFrame(x), d)...)
        df[!, :date] = Date.(df[!, :date])
        df[!, :symbol] .= symbol
        return df[!, colnames]
    else
        coltypes = [
            String,
            Date,
            String,
            Float64,
            Float64,
            Float64,
            Float64,
            Int64,
            # Float64,
            # Float64,
            # Float64,
            # Float64,
            # Int64,
        ]
        return DataFrame(coltypes, colnames, 0)
    end
end

"""
    batch(conn::Connection, symbol::String, types::String; symbols="",
        range="1m"
    )

Batch various API endpoints for a symbol or list of symbols

Path Parameters

|Option | Description
|:----- | :----------
|symbol | Use 'market' to get_stocks multiple symbols (i.e. .../market/batch?...)

Query String Parameters

|Option  | Details
|:-------|:--------
|types   | Required
|        | Comma delimited list of endpoints to call. The names should match
|        | the individual endpoint names. Limited to 10 endpoints.
|        | should be one of chart,quote,news
|symbols | Optional
|        | Comma delimited list of symbols limited to 100. This parameter is
|        | used only if market option is used.
|        | ***(NOTE: If you do not include a symbols parameter the default
|        | behavior is to return the data for ALL relevant symbols. Please be
|        | cautious when using this behavior as it can drive up message usage
|        | quickly.)***
|range   | Optional
|        | Used to specify a chart range if chart is used in types parameter.
|*       | Optional
|        | Parameters that are sent to individual endpoints can be specified
|        | in batch calls and will be applied to each supporting endpoint.
|        | For example, last can be used for the news endpoint to specify the
|        | number of articles
"""
function batch(
    conn::Connection,
    symbol::String,
    types::String;
    symbols::String = "",
    range = "1m",
    filter = std_cols
)
    path_params = [:symbol, :range, :date]

    if symbol == "market"
        if isempty(symbols)
            throw(ArgumentError("batch requires at least 1 symbol"))
        end
    end

    return get_stocks(
        conn, "batch", path_params;
        symbol = symbol,
        symbols = symbols,
        types = types,
        range = range,
        filter = filter
    )
end

"""

    batched_closes(conn::Connection, symbol::String, symbols::Array{String,1};
        range="1m")

get closes for multiple symbols
"""
function batched_closes(conn::Connection, symbols::Array{String,1}; range="1m")

    hist = batch(conn, "market", "chart"; symbols=join(symbols, ", "),
        range=range, filter="date,close")
    df = []
    for symbol in symbols
        tmp_df = vcat([DataFrame(x) for x in hist[symbol]["chart"]]...)
        tmp_df[!, :symbol] .= symbol
        push!(df, tmp_df)
    end
    vcat(df...)
end

"""
    book(conn::Connection, symbol::String)

Response includes data from deep and quote. Refer to each endpoint for details.
"""
function book(conn::Connection, symbol::String)
    path_params = [:symbol]
    return get_stocks(conn, "book", path_params; symbol = symbol)
end

"""
    symbol_quote(conn::Connection, symbol::String; field::String="latestPrice")
"""
function symbol_quote(conn::Connection, symbol::String; field::String="latestPrice")
    path_params = [:symbol, :field]
    return get_stocks(conn, "quote", path_params; symbol = symbol, field = field)
end

"""
    price(conn::Connection, symbol::String)

Response includes data from deep and quote. Refer to each endpoint for details.
"""
function price(conn::Connection, symbol::String)
    path_params = [:symbol]
    return get_stocks(conn, "price", path_params; symbol = symbol)
end

"""
    ohlc(conn::Connection, symbol::String)

Response includes data from deep and quote. Refer to each endpoint for details.
"""
function ohlc(conn::Connection, symbol::String)
    path_params = [:symbol]
    return get_stocks(conn, "ohlc", path_params; symbol = symbol)
end

"""
    largest_trades(conn::Connection, symbol::String)

Returns 15 minute delayed, last sale eligible trades.
"""
function largest_trades(conn::Connection, symbol::String)
    path_params = [:symbol]
    return get_stocks(conn, "largest-trades", path_params; symbol = symbol)
end

max_stars_repo_name: GodolphinCapitalManagement/Iex

max_issues_repo_name: GodolphinCapitalManagement/Iex

max_forks_repo_name: GodolphinCapitalManagement/Iex


text:
module SubTriangulation

using Test
using Gridap
using STLCutters

import Downloads

function download(id;path="")
  url = "https://www.thingiverse.com/download:$id"
  r = Downloads.request(url)
  if 200 ≤ r.status ≤ 399
    _,ext = splitext(r.url)
    mkpath(path)
    filename = joinpath(path,"$id$ext")
    Downloads.download(url,filename)
  else
    @warn "$id is no longer available on Thingiverse."
    filename = nothing
  end
  filename
end

function main(filename;n=20,δ=0.2,output=nothing)

  println("Running: $(basename(filename))")

  stl = STLGeometry( filename )

  pmin,pmax = get_bounding_box(stl)
  Δ = (pmax-pmin)*δ
  pmin -= Δ
  pmax += Δ
  partition = (n,n,n)
  model = CartesianDiscreteModel(pmin,pmax,partition)

  subcells,subfaces,labels = subtriangulate(model,stl)

  if output !== nothing
    sf = output*"_subfaces"
    sc = output*"_subcells"
    bg = output*"_bgmesh"
    writevtk(subfaces,sf,celldata=["bgcell"=>labels.face_to_bgcell])
    writevtk(subcells,sc,celldata=["inout"=>labels.cell_to_io,"bgcell"=>labels.cell_to_bgcell])
    writevtk(get_grid(model),bg,celldata=["inoutcut"=>labels.bgcell_to_ioc])
  end

end

end # module


text:
using JuMP
using LinearAlgebra

"""
    build_cgcp(optimizer, ...)

Build a conic cut-generating problem to separate a split cut.

Original problem writes
```
min     c'x
s.t.    A x = b
          x ∈ K
```
and, for the split `{π'x ⩽ π0} ∧ {π'x ⩾ π0 +1}`.

The CGCP in full writes (without the normalization)
```
min     α'x - β
s.t.    α =       λ - u0 * π
        α = A'v + μ + v0 * π
        β ⩽         - u0 * π0
        β ⩽ b'v     + v0 * (π0 + 1)
        λ, μ ∈ K*, u0, v0 ⩾ 0
```
and in compact form
```
min     x'λ - u0 * (x'π - π0) + s
s.t.    A'v + (μ - λ) + (u0 + v0) * π       = 0
        b'v + (s - t) + (u0 + v0) * π0 + v0 = 0
        λ, μ ∈ K*, s, t, u0, v0 ⩾ 0
```

# Arguments
* `optimizer`: Solver for the CGCP. Must support conic constraints
* `m::Int`: Number of equality constraints
* `n::Int`: Number of variables
* `A::AbstractMatrix`: the constraint matrix
* `b::Vector`: The right-hand side of equality constraints
* `cones::Vector{Tuple{Vector{Int}, MOI.AbstractSet}}`: Cones
* `x_::Vector`: the fractional point to separate
* `π::Vector, π0::Float64`: Coefficients of the split
* `nrm::Symbol`: The normalization used in CGCP
    * `:Alpha2`: α-normalization `|α| ⩽ 1`
    * `:Conic`: conic normalization `ρ'λ + ρ'μ + u0 + v0 ⩽ 1` with `ρ ∈ int(K)`
"""
function build_cgcp(
    optimizer::MOI.OptimizerWithAttributes,
    m::Int, n::Int, A::SparseMatrixCSC, b::Vector{Float64}, cones,
    x_::Vector{Float64}, π::Vector{Float64}, π0::Float64;
    nrm::Symbol = :Conic,
    compact::Bool=true,
    bridge_type::Union{Nothing, Type}=nothing,
    xref::Vector{Float64}=Float64[]
)
    # Some solvers do not support Rotated Second Order cones directly,
    # so we enable bridges as a (hopefully short-term) work-around
    cgcp = JuMP.direct_model(MOI.instantiate(optimizer, with_bridge_type=bridge_type))

    # Create variables
    if !compact
        @variable(cgcp, α[1:n])
        @variable(cgcp, β)
    end

    @variable(cgcp, v[1:m])

    @variable(cgcp, u0 >= 0)
    @variable(cgcp, v0 >= 0)

    # Conic Farkas multipliers
    λ = Vector{JuMP.VariableRef}(undef, n)
    μ = Vector{JuMP.VariableRef}(undef, n)

    for (kidx, k) in cones
        kd = MOI.dual_set(k)
        d = MOI.dimension(k)
        if isa(kd, MOI.Nonnegatives)
            @assert d == 1
            λ_ = @variable(cgcp)
            μ_ = @variable(cgcp)
            set_lower_bound(λ_, 0.0)
            set_lower_bound(μ_, 0.0)

        elseif isa(kd, MOI.Nonpositives)
            @assert d == 1
            λ_ = @variable(cgcp)
            μ_ = @variable(cgcp)
            set_upper_bound(λ_, 0.0)
            set_upper_bound(μ_, 0.0)

        elseif isa(kd, MOI.Zeros)
            @assert d == 1
            λ_ = @variable(cgcp)
            μ_ = @variable(cgcp)
            set_lower_bound(λ_, 0.0)
            set_lower_bound(μ_, 0.0)
            set_upper_bound(λ_, 0.0)
            set_upper_bound(μ_, 0.0)

        elseif isa(kd, MOI.Reals)
            @assert d == 1
            λ_ = @variable(cgcp)
            μ_ = @variable(cgcp)

        elseif isa(kd, MOI.SecondOrderCone)
            λ_ = @variable(cgcp, [1:d] in SecondOrderCone())
            μ_ = @variable(cgcp, [1:d] in SecondOrderCone())

        elseif isa(kd, MOI.RotatedSecondOrderCone)
            @assert d >= 2
            λ_ = @variable(cgcp, [1:d] in RotatedSecondOrderCone())
            μ_ = @variable(cgcp, [1:d] in RotatedSecondOrderCone())

        else
            error("Primal/dual cones $((typeof(k), typeof(kd))) not supported in CGCP")
        end
        λ[kidx] .= λ_
        μ[kidx] .= μ_
    end
    cgcp[:λ] = λ
    cgcp[:μ] = μ

    # Objective and constraints
    if compact
        @variable(cgcp, η1 >= 0)
        @variable(cgcp, η2 >= 0)

        @constraint(cgcp, A'v .+ ( μ .- λ) .+ (u0 + v0) .* π       .== 0.0)
        @constraint(cgcp, b'v  - (η2 - η1)  + (u0 + v0)  * π0 + v0  == 0.0)

        @objective(cgcp, Min, dot(x_, λ) - u0 * (dot(x_, π) - π0) + η1)

    else
        # Full form
        @constraint(cgcp, f1a, α .==        λ .- (u0 .*  π))
        @constraint(cgcp, f1b, α .== A'v .+ μ .+ (v0 .*  π))
        @constraint(cgcp, f2a, β  <=           -  u0  *  π0)
        @constraint(cgcp, f2b, β  <= b'v       +  v0  * (π0 + 1))

        @objective(cgcp, Min, dot(x_, α) - β)
    end

    # Normalization
    if nrm == :Alpha2
        add_alpha_normalization(m, n, cgcp, cones, π, p=2)

    elseif nrm == :Interior
        length(xref) == length(x_) || error("xref has length $(length(xref))")

        add_interior_normalization(m, n, cgcp, cones, xref - x_, π)

    elseif nrm == :Trivial
        add_trivial_normalization(m, n, cgcp, cones)

    elseif nrm == :Standard
        add_standard_normalization(m, n, cgcp, cones)

    elseif nrm == :Uniform
        add_uniform_normalization(m, n, cgcp, cones)

    elseif nrm == :Combined
        add_combined_normalization(m, n, cgcp, cones, π)

    else
        error("Normalization $nrm is not supported")
    end
    return cgcp
end

function add_alpha_normalization(m::Int, n::Int, cgcp, cones, π; p::Real=2)

    λ, u0 = cgcp[:λ], cgcp[:u0]
    if p == 2
        @constraint(cgcp, normalization, dot(λ + u0 .* π, λ + u0 .* π) <= 1)
    else
        error("Alpha normalization only supported for p={2} (is $p)")
    end

    return cgcp
end

"""
    Interior normalization

"""
function add_interior_normalization(m::Int, n::Int, cgcp, cones, γ, π)
    λ, u0 = cgcp[:λ], cgcp[:u0]
    @constraint(cgcp, normalization, dot(γ, λ) - u0 * dot(γ, π) <= 1)

    return cgcp
end

"""
    add_trivial_normalization

`u0 + v0 <= 1`
"""
function add_trivial_normalization(m::Int, n::Int, cgcp, cones)
    @constraint(cgcp, normalization, cgcp[:u0] + cgcp[:v0] <= 1)
    return cgcp
end

"""
    add_standard_normalization

`ρ'λ + ρ'μ + u0 + v0 ⩽ 1`
"""
function add_standard_normalization(m::Int, n::Int, cgcp, cones)
    ρ = zeros(n)
    for (kidx, k) in cones
        kd = MOI.dual_set(k)

        if isa(kd, MOI.Nonnegatives)
            ρ[kidx] .= 1
        elseif isa(kd, MOI.Nonpositives)
            ρ[kidx] .= -1
        elseif isa(kd, MOI.Zeros)
            # Nothing to do
        elseif isa(kd, MOI.Reals)
            # Should be treated as equality constraints, i.e.,
            # set one of the multipliers to 0 and let the other one free
            @warn "Some variables are fixed. Should treat them as equality constraints" maxlog = 1
        elseif isa(kd, MOI.SecondOrderCone)
            ρ[kidx[1]] = 1
        elseif isa(kd, MOI.RotatedSecondOrderCone)
            ρ[kidx[1:2]] .= 1
        else
            error("Standard normalization for cone $(typeof(k)) is not supported")
        end
    end

    λ, μ = cgcp[:λ], cgcp[:μ]
    @constraint(cgcp, normalization, dot(ρ, λ) + dot(ρ, μ) + cgcp[:u0] + cgcp[:v0] <= 1)

    return cgcp
end

"""
    add_uniform_normalization

`ρ'λ + ρ'μ ⩽ 1`
"""
function add_uniform_normalization(m::Int, n::Int, cgcp, cones)
    ρ = zeros(n)
    for (kidx, k) in cones
        kd = MOI.dual_set(k)

        if isa(kd, MOI.Nonnegatives)
            ρ[kidx] .= 1
        elseif isa(kd, MOI.Nonpositives)
            ρ[kidx] .= -1
        elseif isa(kd, MOI.Zeros)
            # Nothing to do
        elseif isa(kd, MOI.Reals)
            # Should be treated as equality constraints, i.e.,
            # set one of the multipliers to 0 and let the other one free

        elseif isa(kd, MOI.SecondOrderCone)
            ρ[kidx[1]] = 1
        elseif isa(kd, MOI.RotatedSecondOrderCone)
            ρ[kidx[1:2]] .= 1
        else
            error("Uniform normalization for cone $(typeof(k)) is not supported")
        end
    end

    λ, μ = cgcp[:λ], cgcp[:μ]
    @constraint(cgcp, normalization, dot(ρ, λ) + dot(ρ, μ) <= 1)

    return cgcp
end

function add_combined_normalization(m::Int, n::Int, cgcp, cones, π)
    ρ = zeros(n)
    for (kidx, k) in cones
        kd = MOI.dual_set(k)

        if isa(kd, MOI.Nonnegatives)
            ρ[kidx] .= 1
        elseif isa(kd, MOI.Nonpositives)
            ρ[kidx] .= -1
        elseif isa(kd, MOI.Zeros)
            # Nothing to do
        elseif isa(kd, MOI.Reals)
            # Should be treated as equality constraints, i.e.,
            # set one of the multipliers to 0 and let the other one free
        elseif isa(kd, MOI.SecondOrderCone)
            ρ[kidx[1]] = 1
        elseif isa(kd, MOI.RotatedSecondOrderCone)
            ρ[kidx[1:2]] .= 1
        else
            error("Uniform normalization for cone $(typeof(k)) is not supported")
        end
    end

    λ, μ = cgcp[:λ], cgcp[:μ]
    u0, v0 = cgcp[:u0], cgcp[:v0]

    if iszero(dot(ρ, π))
        @constraint(cgcp, normalization, dot(ρ, λ) <= 1)
    else
        @constraint(cgcp, normalization, dot(ρ, λ) + u0 <= 1)
    end

    @constraint(cgcp, normalization2, dot(ρ, λ) + dot(ρ, μ) + u0 + v0 <= 100)

    return cgcp
end

max_stars_repo_name: mtanneau/CLaP

max_issues_repo_name: mtanneau/CLaP

max_forks_repo_name: mtanneau/CLaP


text:
# some mem tests
mod=randn(1000); mod0=2.0
@btime GeoPhyInv.Models.χ!(mod,mod0,1)
@btime GeoPhyInv.Models.χ!(mod,mod0,-1)
@btime GeoPhyInv.Models.χg!(mod,mod0,1)
@btime GeoPhyInv.Models.χg!(mod,mod0,-1)

# some model
mgrid = [range(0.0, stop=10.,step=0.05), range(0.0, stop=10.,step=0.05)];
model = GeoPhyInv.Models.Seismic_zeros(mgrid);
vp0=[2100.,2200.];vs0=[-1., -1.]; ρ0=[2100., 2300.]
GeoPhyInv.Models.adjust_bounds!(model, vp0,vs0,ρ0);

GeoPhyInv.Models.Seismic_addon!(model, randn_perc=1e-3)
nznx=prod(length.(model.mgrid))

# test copyto!
model_new=deepcopy(model)
GeoPhyInv.Models.Seismic_addon!(model_new, randn_perc=1e-3)

@btime copyto!(model_new, model)
@test isequal(model_new, model)

# allocation of model0
model0 = GeoPhyInv.Models.Seismic_zeros(model.mgrid);
GeoPhyInv.Models.adjust_bounds!(model0, model);

for attribvec in [[:χKI,:χρI,:null], [:χKI,:null,:null], [:null,:χρ,:null], [:null,:χρ,:null], [:χvp,:χρ,:null]]
	x=zeros(nznx*count(attribvec.≠:null));
	# reparameterize twice to see conversion formulae are correct
        GeoPhyInv.Models.Seismic_get!(x, model, attribvec);
	@time GeoPhyInv.Models.Seismic_reparameterize!(model0, x, attribvec);
	@test model.χvp ≈ model0.χvp
	#@test model.χvs ≈ model0.χvs
	@test model.χρ ≈ model0.χρ
end

nznx=prod(length.(model.mgrid))
G1 = randn(nznx); G2 = randn(nznx);
g1 = similar(G1); g2 = similar(G2);

gmodel = GeoPhyInv.Models.Seismic_zeros(model.mgrid);

for attribvec in [[:χKI,:χρI,:null], [:χKI,:null,:null], [:χvp,:χρ,:null] , [:null,:χρ,:null], [:χvp,:null,:null]]
	G=randn(nznx*count(attribvec.≠:null));
	g=zeros(nznx*count(attribvec.≠:null));
	# chain rule 1
	@time GeoPhyInv.Models.Seismic_chainrule!(gmodel, model, G, attribvec,1);

	# chain rule 2
	@time GeoPhyInv.Models.Seismic_chainrule!(gmodel, model, g, attribvec, -1);
	@test G ≈ g
end


text:
using Documenter, Gen

makedocs(
    format = :html,
    sitename = "Gen",
    modules = [Gen],
    pages = [
        "Home" => "index.md",
        "Getting Started" => "getting_started.md",
        "Tutorials" => "tutorials.md",
        "Modeling Languages and APIs" => [
            "Generative Functions" => "ref/gfi.md",
            "Probability Distributions" => "ref/distributions.md",
            "Built-in Modeling Language" => "ref/modeling.md",
            "Generative Function Combinators" => "ref/combinators.md",
            "Choice Maps" => "ref/choice_maps.md",
            "Selections" => "ref/selections.md",
            "Optimizing Trainable Parameters" => "ref/parameter_optimization.md",
            "Extending Gen" => "ref/extending.md"
         ],
        "Standard Inference Library" => [
            "Importance Sampling" => "ref/importance.md",
            "Markov chain Monte Carlo" => "ref/mcmc.md",
            "MAP Optimization" => "ref/map.md",
            "Particle Filtering" => "ref/pf.md",
            "Variational Inference" => "ref/vi.md",
            "Learning Generative Functions" => "ref/learning.md"
        ],
        "Internals" => [
            "Optimizing Trainable Parameters" => "ref/internals/parameter_optimization.md"
         ]
    ]
)

deploydocs(
    repo = "github.com/probcomp/Gen.git",
    target = "build"
)