text:
###############################################################################
#
#  Twists
#
#
#  (C) 2022 Jeroen Hanselman
#
###############################################################################

export quadratic_twist, quadratic_twists, is_twist, twists

@doc Markdown.doc"""
    quadratic_twist(E::EllCrv{T}}, d::T) -> EllCrv{T}

Compute the quadratic twist of $E$ by $d$.
"""
#Needs more testing over characteristic 2. (Magma's twists don't always work over char 2.
#Adapted from Connell's Handbook of elliptic curves.
function quadratic_twist(E::EllCrv{T}, d) where T<: FieldElem

  a1, a2, a3, a4, a6 = a_invars(E)
  K = base_field(E)
  if characteristic(K) != 2
    return EllipticCurve(K, [a1, a2*d + a1^2*(d-1)//4, a3, a4*d^2 + a1*a3*(d^2-1)//2, a6*d^3 + a3^2*(d^3 -1)//4])
  end

  return EllipticCurve(K, [a1, a2+a1^2*d, a3, a4, a6 + a3^2])

end

@doc Markdown.doc"""
    quadratic_twist(E::EllCrv{FinFieldElem}}) -> EllCrv{FinFieldElem}

Compute the unique quadratic twist of $E$.
"""
function quadratic_twist(E::EllCrv{T}) where T<: FieldElem

  K = base_field(E)
  char = characteristic(K)

 if char == 2
   f, h = hyperelliptic_polynomials(E)
   if iseven(degree(K))
     u = normal_basis(GF(Int(char), 1),K)
   else
     u = one(K)
   end

   return EllipticCurve(f + u*h^2, h)
 end

  a = gen(K)
  return quadratic_twist(E, a)

end

#Test if we can't sometimes get two isomorphic curves
@doc Markdown.doc"""
    quadratic_twists(E::EllCrv{FinFieldElem}}) -> Vector{EllCrv{FinFieldElem}}

Compute all quadratic twists of $E$.
"""
function quadratic_twists(E::EllCrv{T}) where T<: FinFieldElem

  return [E, quadratic_twist(E)]

end

function supersingular_twists2(E::EllCrv{T}) where T<: FinFieldElem
#Adapted from Magma code
  K = base_field(E)
  if isodd(degree(K))
    return EllCrv{T}[EllipticCurve(K, [0, 0, 1, 0, 0]), EllipticCurve(K, [0, 0, 1, 1, 0]), EllipticCurve(K, [0, 0, 1, 1, 1]) ]
  end

  u = gen(K);
  c = u
  while absolute_tr(c) == 0
    c = rand(K)
  end
  #First three curves
  part_1 = EllCrv{T}[EllipticCurve(K, [0, 0, 1, 0, 0]), EllipticCurve(K, [0, 0, 1, c, 0]), EllipticCurve(K, [0, 0, 1, 0, c]) ]
  #The other four
  part_2 = EllCrv{T}[EllipticCurve(K, [0, 0, a, 0, a^2*b]) for (a, b) in [[u, 0], [inv(u), 0], [u, c], [inv(u), c]]]
  return vcat(part_1, part_2)
end

function supersingular_twists3(E::EllCrv{T}) where T<: FinFieldElem
#Adapted from Magma code
    K = base_field(E)
    d = degree(K)

    if mod(d, 3) != 0
        c = one(K)
    else
      c = gen(K)
      while absolute_tr(c) == 0
        c = rand(K)
      end
    end

    if isodd(d)
      return EllCrv{T}[EllipticCurve(K, [1,0]), EllipticCurve(K, [-1, 0]), EllipticCurve(K, [-1, c]), EllipticCurve(K, [-1, -c]) ];
    end
    u = gen(K);
    #First four curves
    part_1 = EllCrv{T}[EllipticCurve(K, [-u^i,0]) for i in (0:3)]
    part_2 = EllCrv{T}[EllipticCurve(K, [-1,c]), EllipticCurve(K, [-u^2,(u^3)*c])]
    return vcat(part_1, part_2)

end

@doc Markdown.doc"""
    twists(E::EllCrv{FinFieldElem}}) -> Vector{EllCrv{FinFieldElem}}

Compute all twists of $E$.
"""
function twists(E::EllCrv{T}) where T<: FinFieldElem
#Adapted from Magma code
   K = base_field(E);
   p = characteristic(K)
   j = j_invariant(E)
   if j != 0 && j != 1728
      return EllCrv{T}[E, quadratic_twist(E)]
   elseif p == 2
      return supersingular_twists2(E)
   elseif p == 3
      return supersingular_twists3(E)
   elseif j == 0
      a = gen(K)
      c4, c6 = c_invars(E)
      c = -c6//864
      n = gcd(6, order(K)-1)
      return EllCrv{T}[ EllipticCurve(K, [0,c*a^i]) for i in (0:n-1) ]
   elseif j == 1728
      a = gen(K)
      c4, c6 = c_invars(E)
      c = -c4//48;
      n = gcd(4, order(K)-1)
      return EllCrv{T}[ EllipticCurve(K, [c*a^i,0]) for i in (0:n-1)]
   end
   #This can never happen:
   error("There is a bug in twists.")
end

@doc Markdown.doc"""
    is_twist(EEllCrv, EllCrv) -> Bool

Return true if $F$ is a twist of $E$
"""
function is_twist(E::EllCrv, F::EllCrv)

  return j_invariant(E) == j_invariant(F)
end


text:
"""
    GammaLikelihood(α::Real=1.0, l=exp)

Gamma likelihood with fixed shape `α`.

```math
    p(y|f) = \\operatorname{Gamma}(y | α, l(f))
```
On calling, this returns a Gamma distribution with shape `α` and scale `invlink(f)`.
"""
struct GammaLikelihood{Tl<:AbstractLink,T<:Real} <: AbstractLikelihood
    α::T    # shape parameter
    invlink::Tl
end

GammaLikelihood(l) = GammaLikelihood(1.0, l)
GammaLikelihood(α::Real=1.0, l=exp) = GammaLikelihood(α, link(l))

@functor GammaLikelihood

(l::GammaLikelihood)(f::Real) = Gamma(l.α, l.invlink(f))

(l::GammaLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs))

function expected_loglikelihood(
    ::AnalyticExpectation,
    lik::GammaLikelihood{ExpLink},
    q_f::AbstractVector{<:Normal},
    y::AbstractVector{<:Real},
)
    f_μ = mean.(q_f)
    return sum(
        (lik.α - 1) * log.(y) .- y .* exp.((var.(q_f) / 2) .- f_μ) .- lik.α * f_μ .-
        loggamma(lik.α),
    )
end

default_expectation_method(::GammaLikelihood{ExpLink}) = AnalyticExpectation()

ext: jl


text:
__precompile__()

module PlotlyJS

using Compat; import Compat: String, readstring, view
using JSON
using Blink
using Colors
using DocStringExtensions

# import LaTeXStrings and export the handy macros
using LaTeXStrings
export @L_mstr, @L_str

# export some names from JSON
export json

# globals for this package
const _js_path = joinpath(dirname(dirname(@__FILE__)),
                          "deps", "plotly-latest.min.js")
const _js_cdn_path = "https://cdn.plot.ly/plotly-latest.min.js"
const _mathjax_cdn_path = "https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_SVG"

const _autoresize = [true]
autoresize(b::Bool) = (_autoresize[1] = b; b)
autoresize() = _autoresize[1]

_isijulia() = isdefined(Main, :IJulia) && Main.IJulia.inited

# include these here because they are used below
include("traces_layouts.jl")
include("styles.jl")
abstract AbstractPlotlyDisplay

# core plot object
type Plot{TT<:AbstractTrace}
    data::Vector{TT}
    layout::AbstractLayout
    divid::Base.Random.UUID
    style::Style
end

# include the rest of the core parts of the package
include("display.jl")
include("json.jl")
include("subplots.jl")
include("api.jl")
include("savefig.jl")
include("convenience_api.jl")

# Set some defaults for constructing `Plot`s
function Plot(;style::Style=DEFAULT_STYLE[1])
    Plot(GenericTrace{Dict{Symbol,Any}}[], Layout(), Base.Random.uuid4(), style)
end

function Plot{T<:AbstractTrace}(data::Vector{T}, layout=Layout();
                                style::Style=DEFAULT_STYLE[1])
    Plot(data, layout, Base.Random.uuid4(), style)
end

function Plot(data::AbstractTrace, layout=Layout();
              style::Style=DEFAULT_STYLE[1])
    Plot([data], layout; style=style)
end

function docs()
    schema_path = joinpath(dirname(dirname(@__FILE__)), "deps", "schema.html")
    w = Blink.Window()
    Blink.content!(w, "html", open(readstring, schema_path), fade=false)
end

# NOTE: we export trace constructing types from inside api.jl
# NOTE: we export names of shapes from traces_layouts.jl
export

    # core types
    Plot, GenericTrace, Layout, Shape, ElectronDisplay, JupyterDisplay,
    ElectronPlot, JupyterPlot, AbstractTrace, AbstractLayout,

    # other methods
    savefig, svg_data, png_data, jpeg_data, webp_data, autoresize,

    # plotly.js api methods
    restyle!, relayout!, addtraces!, deletetraces!, movetraces!, redraw!,
    extendtraces!, prependtraces!,

    # non-!-versions (forks, then applies, then returns fork)
    restyle, relayout, addtraces, deletetraces, movetraces, redraw,
    extendtraces, prependtraces,

    # helper methods
    plot, fork, vline, hline, attr,

    # frontend methods
    init_notebook,

    # styles
    use_style!, style, Style

function __init__()
    # --------------------------------------------- #
    # Code to run once when the notebook starts up! #
    # --------------------------------------------- #

    if _isijulia()

        init_notebook()

        @eval begin
            import IJulia
            import IJulia.CommManager: Comm, send_comm
        end

        # set up the comms we will use to send js messages to be executed
        global const _ijulia_eval_comm = Comm(:plotlyjs_eval)
        global const _ijulia_return_comms = ObjectIdDict()

        @eval begin
            function IJulia.display_dict(p::JupyterPlot)
                if p.view.displayed
                    Dict()
                else
                    p.view.displayed = true
                    Dict("text/html" => html_body(p))
                end
            end

            SyncPlot(p::Plot) = SyncPlot(p, JupyterDisplay(p))

            IJulia.display_dict(p::Plot) =
                Dict("text/plain" => sprint(writemime, "text/plain", p))

        end
    else
        @eval SyncPlot(p::Plot) = SyncPlot(p, ElectronDisplay(p))
    end

end

end # module


text:
# shignal_stft.jl
##########################
fs=Int(1000);
Ts=1/fs;
N=Int(12000);
T=N/fs;
#t=(0:Int(N-1))*Ts;

x=zeros(Complex, N)
for m=1:N
    x[m]=sin(2*Ts*(m-1)*(Ts*(m-1)-3)*(Ts*(m-1)-6)*(Ts*(m-1)-9)*(Ts*(m-1)-12))
end

# setting integers
L=256;
M1=128;
M=2*M1;
N0=250;
n0=1+floor(Int, (N-L)/(L-N0));
s=L-N0;

using DSP, FFTW

# window
h=hanning(L);

using LinearAlgebra

# the N-th root of unity
omega=exp(2*pi*im/M);

# sqrt(N) times the Fourier matrix
F=zeros(Complex, M, M);
for m=1:M
    for n=1:M
        F[m,n]=omega^(-(m-1)*(n-1))
    end
end

# sampling and locarization
Y=zeros(Complex, L,n0);
for l=1:L
    for j=1:n0
        Y[l,j]=h[l]*x[s*(j-1)+l]
    end
end

# discrete short-time Fourier transform
Z=F*Y;

# spectrogram
X=zeros(Float64, M1+1,n0);
for k=1:M1+1
    for j=1:n0
        X[k,j]=abs.(Y[M1+2-k,j])
    end
end

using Plots

# plotting
plot(X,
    title="Spectrogram of x(t)=sin(2t(t-3)(t-6)(t-9)(t-12))",
    st=:heatmap,
    color=:cool,
    xaxis="n = 1, 1, 2, ..., n0-1",
    yaxis="k = 0, 1, 2, ..., M/2")
plot!(xticks = ([0 1000 1957;], [0 1000 1957]))
plot!(yticks = ([0 50 100;], [0 50 100]))
savefig("signal_stft.png")


text:
const AffinityPropagation_ = skcl(:AffinityPropagation)
@sk_uns mutable struct AffinityPropagation <: MMI.Unsupervised
    damping::Float64      = 0.5::(0.5 ≤ _ ≤ 1)
    max_iter::Int         = 200::(_ ≥ 1)
    convergence_iter::Int = 15::(_ ≥ 1)
    copy::Bool            = true
    preference::Any       = nothing
    affinity::String      = "euclidean"::(_ in ("euclidean", "precomputed"))
    verbose::Bool         = false
end
@sku_predict AffinityPropagation
function MMI.fitted_params(m::AffinityPropagation, f)
    nc   = length(f.cluster_centers_indices_)
    catv = MMI.categorical(1:nc)
    return (
        cluster_centers_indices = f.cluster_centers_indices_,
        cluster_centers         = f.cluster_centers_,
        labels                  = nc == 0 ? nothing : catv[f.labels_ .+ 1],
        affinity_matrix         = f.affinity_matrix_)
end
meta(AffinityPropagation,
    input   = Table(Continuous),
    target  = AbstractVector{Multiclass},
    # no transform so no output
    weights = false,
    descr   = "Perform Affinity Propagation Clustering of data."
    )

# ============================================================================
const AgglomerativeClustering_ = skcl(:AgglomerativeClustering)
@sk_uns mutable struct AgglomerativeClustering <: MMI.Unsupervised
    n_clusters::Int     = 2::(_ ≥ 1)
    affinity::String    = "euclidean"::(_ in ("euclidean", "l1", "l2", "manhattan", "cosine", "precomputed"))
    memory::Any         = nothing
    connectivity::Any   = nothing
    compute_full_tree::Union{String,Bool} = "auto"::(_ isa Bool || _ == "auto")
    linkage::String     = "ward"::(_ in ("ward", "complete", "average", "single"))
    distance_threshold::Option{Float64}   = nothing::(_ === nothing || _ > 0)
end
function MMI.fitted_params(m::AgglomerativeClustering, f)
    nc   = f.n_clusters_
    catv = MMI.categorical(1:nc)
    return (
        n_clusters = f.n_clusters_,
        labels     = catv[f.labels_ .+ 1],
        n_leaves   = f.n_leaves_,
        n_connected_components = f.n_connected_components_,
        children   = f.children_)
end
meta(AgglomerativeClustering,
    input   = Table(Continuous),
    # no predict nor transform so no target nor output
    weights = false,
    descr   = "Recursively merges the pair of clusters that minimally increases a given linkage distance. Note: there is no `predict` or `transform` method
    associated with it; instead, inspect the `fitted_params`."
    )

# ============================================================================
const Birch_ = skcl(:Birch)
@sk_uns mutable struct Birch <: MMI.Unsupervised
    threshold::Float64    = 0.5::(_ > 0)
    branching_factor::Int = 50::(_ > 0)
    n_clusters::Int       = 3::(_ > 0)
    compute_labels::Bool  = true
    copy::Bool            = true
end
@sku_predict Birch
@sku_transform Birch
function MMI.fitted_params(m::Birch, f)
    nc   = m.n_clusters
    catv = MMI.categorical(1:nc)
    return (
        root               = f.root_,
        dummy_leaf         = f.dummy_leaf_,
        subcluster_centers = f.subcluster_centers_,
        subcluster_labels  = f.subcluster_labels_,
        labels             = catv[f.labels_ .+ 1])
end
meta(Birch,
    input   = Table(Continuous),
    target  = AbstractVector{Multiclass},
    output  = Table(Continuous),
    weights = false,
    descr   = "Memory-efficient, online-learning algorithm provided as an alternative to MiniBatchKMeans. Note: noisy samples are given the label -1."
    )

# ============================================================================
const DBSCAN_ = skcl(:DBSCAN)
@sk_uns mutable struct DBSCAN <: MMI.Unsupervised
    eps::Float64        = 0.5::(_ > 0)
    min_samples::Int    = 5::(_ > 1)
    metric::String      = "euclidean"::(_ in ("euclidean", "precomputed"))
    algorithm::String   = "auto"::(_ in ("auto", "ball_tree", "kd_tree", "brute"))
    leaf_size::Int      = 30::(_ > 1)
    p::Option{Float64}  = nothing
    n_jobs::Option{Int} = nothing
end
function MMI.fitted_params(m::DBSCAN, f)
    nc   = length(f.core_sample_indices_)
    catv = MMI.categorical([-1, (1:nc)...])
    return (
        core_sample_indices = f.core_sample_indices_,
        components          = f.components_,
        labels              = catv[f.labels_ .+ 2])
end
meta(DBSCAN,
    input   = Table(Continuous),
    weights = false,
    descr   = "Density-Based Spatial Clustering of Applications with Noise. Finds core samples of high density and expands clusters from them. Good for data which contains clusters of similar density."
    )

# ============================================================================
const FeatureAgglomeration_ = skcl(:FeatureAgglomeration)
@sk_uns mutable struct FeatureAgglomeration <: MMI.Unsupervised
    n_clusters::Int        = 2::(_ > 0)
    memory::Any            = nothing
    connectivity::Any      = nothing
    # XXX unclear how to pass a proper callable here; just passing mean = nok
    # pooling_func::Function = mean
    affinity::Any          = "euclidean"::(_ isa Function || _ in ("euclidean", "l1", "l2", "manhattan", "cosine",  "precomputed"))
    compute_full_tree::Union{String,Bool} = "auto"::(_ isa Bool || _ == "auto")
    linkage::String        = "ward"::(_ in ("ward", "complete", "average", "single"))
    distance_threshold::Option{Float64}   = nothing
end
@sku_transform FeatureAgglomeration
@sku_inverse_transform FeatureAgglomeration
function MMI.fitted_params(m::FeatureAgglomeration, f)
    nc   = m.n_clusters
    catv = MMI.categorical(1:nc)
    return (
        n_clusters = f.n_clusters_,
        labels     = catv[f.labels_ .+ 1],
        n_leaves   = f.n_leaves_,
        n_connected_components = f.n_connected_components_,
        children   = f.children_,
        distances  = m.distance_threshold === nothing ? nothing : f.distances_)
end
meta(FeatureAgglomeration,
    input   = Table(Continuous),
    output  = Table(Continuous),
    weights = false,
    descr   = "Similar to AgglomerativeClustering, but recursively merges features instead of samples."
    )

# ============================================================================
const KMeans_ = skcl(:KMeans)
@sk_uns mutable struct KMeans <: MMI.Unsupervised
    n_clusters::Int     = 8::(_ ≥ 1)
    n_init::Int         = 10::(_ ≥ 1)
    max_iter::Int       = 300::(_ ≥ 1)
    tol::Float64        = 1e-4::(_ > 0)
    verbose::Int        = 0::(_ ≥ 0)
    random_state::Any   = nothing
    copy_x::Bool        = true
    n_jobs::Option{Int} = nothing
    algorithm::String   = "auto"::(_ in ("auto", "full", "elkane"))
    # long
    precompute_distances::Union{Bool,String} = "auto"::(_ isa Bool || _ == "auto")
    init::Union{AbstractArray,String}        = "k-means++"::(_ isa AbstractArray || _ in ("k-means++", "random"))
end
@sku_transform KMeans
@sku_predict KMeans
function MMI.fitted_params(m::KMeans, f)
    nc   = m.n_clusters
    catv = MMI.categorical(1:nc)
    return (
        cluster_centers = f.cluster_centers_,
        labels          = catv[f.labels_ .+ 1],
        inertia         = f.inertia_)
end
meta(KMeans,
    input   = Table(Continuous),
    target  = AbstractVector{Multiclass},
    output  = Table(Continuous),
    weights = false,
    descr   = "K-Means algorithm: find K centroids corresponding to K clusters in the data."
    )

# ============================================================================
const MiniBatchKMeans_ = skcl(:MiniBatchKMeans)
@sk_uns mutable struct MiniBatchKMeans <: MMI.Unsupervised
    n_clusters::Int         = 8::(_ ≥ 1)
    max_iter::Int           = 100::(_ > 1)
    batch_size::Int         = 100::(_ > 1)
    verbose::Int            = 0
    compute_labels::Bool    = true
    random_state::Any       = nothing
    tol::Float64            = 0.0::(_ ≥ 0)
    max_no_improvement::Int = 10::(_ > 1)
    init_size::Option{Int}  = nothing
    n_init::Int             = 3::(_ > 0)
    init::Union{AbstractArray,String} = "k-means++"::(_ isa AbstractArray || _ in ("k-means++", "random"))
    reassignment_ratio::Float64       = 0.01::(_ > 0)
end
@sku_predict MiniBatchKMeans
@sku_transform MiniBatchKMeans
function MMI.fitted_params(m::MiniBatchKMeans, f)
    nc   = m.n_clusters
    catv = MMI.categorical(1:nc)
    return (
        cluster_centers = f.cluster_centers_,
        labels          = catv[f.labels_ .+ 1],
        inertia         = f.inertia_)
end
meta(MiniBatchKMeans,
    input   = Table(Continuous),
    target  = AbstractVector{Multiclass},
    output  = Table(Continuous),
    weights = false,
    descr   = "Mini-Batch K-Means clustering."
    )

# ============================================================================
const MeanShift_ = skcl(:MeanShift)
@sk_uns mutable struct MeanShift <: MMI.Unsupervised
    bandwidth::Option{Float64}   = nothing
    seeds::Option{AbstractArray} = nothing
    bin_seeding::Bool            = false
    min_bin_freq::Int            = 1::(_ ≥ 1)
    cluster_all::Bool            = true
    n_jobs::Option{Int}          = nothing
    # max_iter::Int                = 300::(_ > 1)
end
@sku_predict MeanShift

function MMI.fitted_params(m::MeanShift, f)
    nc   = size(f.cluster_centers_, 1)
    catv = MMI.categorical(1:nc)
    return (
        cluster_centers = f.cluster_centers_,
        labels          = catv[f.labels_ .+ 1])
end
meta(MeanShift,
    input   = Table(Continuous),
    target  = AbstractVector{Multiclass},
    weights = false,
    descr   = "Mean shift clustering using a flat kernel. Mean shift clustering aims to discover \"blobs\" in a smooth density of samples. It is a centroid-based algorithm, which works by updating candidates for centroids to be the mean of the points within a given region. These candidates are then filtered in a post-processing stage to eliminate near-duplicates to form the final set of centroids."
    )

# ============================================================================
const OPTICS_ = skcl(:OPTICS)
@sk_uns mutable struct OPTICS <: MMI.Unsupervised
    min_samples::Union{Float64,Int} = 5::((_ isa Int && _ > 1) || 0 < _ < 1)
    max_eps::Float64       = Inf
    metric::String         = "minkowski"::(_ in ("precomputed", "cityblock", "cosine", "euclidean", "l1", "l2", "manhattan", "braycurtis", "canberra", "chebyshev", "correlation", "dice", "hamming", "jaccard", "kulsinski", "mahalanobis", "minkowski", "rogerstanimoto", "russellrao", "seuclidean", "sokalmichener", "sokalsneath", "sqeuclidean", "yule"))
    p::Int                 = 2::(_ ≥ 1)
    cluster_method::String = "xi"
    eps::Option{Float64}   = nothing
    xi::Float64            = 0.05::(0 < _ < 1)
    predecessor_correction::Bool = true
    min_cluster_size::Union{Nothing,Float64,Int} = nothing::(_ === nothing || (_ isa Int && _ > 1) || (_ isa Float64 && 0 < _ < 1))
    algorithm::String      = "auto"::(_ in ("auto", "ball_tree", "kd_tree", "brute"))
    leaf_size::Int         = 30::(_ > 1)
    n_jobs::Option{Int}    = nothing
end
function MMI.fitted_params(m::OPTICS, f)
    nc   = size(f.cluster_hierarchy_, 1)
    catv = MMI.categorical([-1, (1:nc)...])
    return (
        labels            = catv[f.labels_ .+ 2],
        reachability      = f.reachability_,
        ordering          = f.ordering_,
        core_distances    = f.core_distances_,
        predecessor       = f.predecessor_,
        cluster_hierarchy = f.cluster_hierarchy_)
end
meta(OPTICS,
    input   = Table(Continuous),
    weights = false,
    descr   = "OPTICS (Ordering Points To Identify the Clustering Structure), closely related to DBSCAN, finds core sample of high density and expands clusters from them. Unlike DBSCAN, keeps cluster hierarchy for a variable neighborhood radius. Better suited for usage on large datasets than the current sklearn implementation of DBSCAN."
    )

# ============================================================================
const SpectralClustering_ = skcl(:SpectralClustering)
@sk_uns mutable struct SpectralClustering <: MMI.Unsupervised
    n_clusters::Int      = 8::(_ ≥ 1)
    eigen_solver::Option{String} = nothing::(_ === nothing || _ in ("arpack", "lobpcg", "amg"))
#    n_components::Option{Int}    = nothing::(_ === nothing || _ ≥ 1)
    random_state::Any     = nothing
    n_init::Int           = 10::(_ ≥ 1)
    gamma::Float64        = 1.0::(_ > 0)
    affinity::String      = "rbf"::(_ in ("nearest_neighbors", "rbf", "precomputed", "precomputed_nearest_neighbors"))
    n_neighbors::Int      = 10::(_ > 0)
    eigen_tol::Float64    = 0.0::(_ ≥ 0)
    assign_labels::String = "kmeans"::(_ in ("kmeans", "discretize"))
    n_jobs::Option{Int}   = nothing
end
function MMI.fitted_params(m::SpectralClustering, f)
    nc   = m.n_clusters
    catv = MMI.categorical(1:nc)
    return (
        labels          = catv[f.labels_ .+ 1],
        affinity_matrix = f.affinity_matrix_)
end
meta(SpectralClustering,
    input   = Table(Continuous),
    weights = false,
    descr   = "Apply clustering to a projection of the normalized Laplacian.
    In practice Spectral Clustering is very useful when the structure of the individual clusters is highly non-convex or more generally when a measure of the center and spread of the cluster is not a suitable description of the complete cluster. For instance when clusters are nested circles on the 2D plane."
    )

# NOTE: the two models below are weird, not bothering with them for now
# # ============================================================================
# SpectralBiclustering_ = skcl(:SpectralBiclustering)
# @sk_uns mutable struct SpectralBiclustering <: MMI.Unsupervised
# end
#
# # ============================================================================
# SpectralCoclustering_ = skcl(:SpectralCoclustering)
# @sk_uns mutable struct SpectralCoclustering <: MMI.Unsupervised
# end


text:
###################################################
## network/tools/speeds.jl
## tools related to speed computations in the network
###################################################

"""
    `typicalSpeeds` : compute an array of typical speeds for each road-type,
    given network data
    Just iterate the iterative constant speeds method method to convergence
    Speeds are in meters per second
"""
function typicalSpeeds(pb::NetworkData;
                       uniqueSpeed=false,
                       maxTrips=100_000, solverargs...)
    speeds = Float64[]
    times = spzeros(Float64, nv(pb.network.graph), nv(pb.network.graph))

    initialTimes = randomTimes(pb.network);

    s = LimitedPaths(pb, initialTimes, pathsPerTrip = 1, maxTrip=maxTrips)

    logError = Inf
    newLogError = nwTripsLogError(s.timings, pb)
    while newLogError < logError
        # get the speeds
        speeds = constantSpeedSolver(s, uniqueSpeed=uniqueSpeed; solverargs...)
        # compute the times
        times = fixedSpeedTimes(pb.network,speeds)
        # update the state
        updateState!(s, times)

        # compute new error
        logError = newLogError
        newLogError = nwTripsLogError(s.timings, pb)
    end
    println("\n==============================================")
    println("Speeds converged, final error = $(newLogError)")
    return speeds
end


text:
"""
    Localization

A single point identified by SMLM.
"""
struct Localization <: DataEntity
    index::Int            # unique identifier within this dataset
    channel::String
    x::Float64
    y::Float64
    z::Float64
    accuracy::Float64     # computed lateral uncertainty based on detected photons
    frame::Int            # first frame observed
    length::Int           # number of adjacent frames observed over
end


text:
"""
    sourceintegration(x::Vector,env::Environment,limits::Vector{T}) where T
    sourceintegration(x::Vector,env::Environment,limits::Vector{Vector{T}}) where T
    sourceintegration(x::FreqArray,env::Environment,limits::Vector{T}) where T
    sourceintegration(x::FreqArray,env::Environment,limits::Vector{Vector{T}}) where T

Source integration of `x` from `limits` given in cartesian coordinates `[xmin,xmax,ymin,ymax]`.
```
src1 = [-0.1,0.1,-0.1,0.1]
src2 = [-0.2,0.0,-0.1,0.1]
limits = [src1,src2]
srcint = sourceintegration(x,env,limits)
```

Optionally supply keyword `int_thres` to set a threshold on source integration. Default is set to `Inf` which sum all values within integration region.
"""
function sourceintegration(x::Vector{T1},env::Environment,limits::Vector{T2};int_thres=Inf) where {T1,T2}
    xi = findall((env.rx.>=limits[1]) .& (env.rx.<=limits[2]))
    yi = findall((env.ry.>=limits[3]) .& (env.ry.<=limits[4]))
    I = LinearIndices((env.Nx,env.Ny))[xi,yi]
    pos_idx = x[I] .> 0
    if iszero(pos_idx)
        return NaN
    else
        fmax,g_idx = findmax(x[I][pos_idx])
        thres_idx = SPL.(x[I][pos_idx]) .- SPL.(fmax) .<= int_thres
        return sum(x[I][pos_idx][thres_idx])
    end
end

function sourceintegration(x::Vector{T1},env::Environment,limits::Vector{Vector{T2}};int_thres=Inf) where {T1,T2}
    out = Float64[]
    for src in limits
        push!(out,sourceintegration(x,env,src;int_thres=int_thres))
    end
    return out
end

function sourceintegration(x::FreqArray,env::Environment,limits::Vector{Vector{T}};int_thres=Inf) where T
    out = Array{Float64}(undef,length(limits),length(x.fc))
    for i in 1:length(x.fc)
        tmp = Float64[]
        for src in limits
            push!(tmp,sourceintegration(x.arr[:,i],env,src;int_thres=int_thres))
        end
        out[:,i] = tmp
    end
    return FreqArray(out,x.fc)
end

function sourceintegration(x::FreqArray,env::Environment,limits::Vector{T};int_thres=Inf) where T
    out = Array{Union{Missing, Float64}}(undef,length(x.fc))
    for i in 1:length(x.fc)
        out[i] = sourceintegration(x.arr[:,i],env,limits;int_thres=int_thres)
    end
    return FreqArray(out,x.fc)
end

"""
    point_to_region(src,dxdy)

Helper function to get integration limits from `src` coordinates and size of area `dxdy`
```
sources = [(0.1,0.1),(-0.1,0.1),(-0.1,-0.1),(0.1,-0.1)]
dxdy = (0.1,0.1)
limits = point_to_region(sources,dxdy)
srcint = sourceintegration(x,env,limits)
```
"""
function point_to_region(src::NTuple{2},dxdy::Tuple{Vararg{T,N}}) where {T,N}
    N == 1 && (dxdy = (dxdy[1],dxdy[1]))
    xmin = src[1]-dxdy[1]/2
    xmax = src[1]+dxdy[1]/2
    ymin = src[2]-dxdy[2]/2
    ymax = src[2]+dxdy[2]/2
    return [xmin,xmax,ymin,ymax]
end
function point_to_region(src::NTuple{2},dxdy::T) where T <: Real
    point_to_region(src::NTuple{2},(dxdy[1],dxdy[1]))
end
function point_to_region(src::T,dxdy) where T <: AbstractArray
    eltype(src) <: NTuple || (src = Tuple.(src)) # Convertion to array of tuples
    out = Array{Array{Float64,1},1}(undef,0)
    for i = 1:length(src)
        push!(out,point_to_region(src[i],dxdy))
    end
    return out
end


text:
module GeneticAlgorithms

    # -------
    using RandomExtensions
    using Distributed

    include("Types.jl")
    using .Types

    import Dates

    import Base, Base.show, Base.isless

    export  runga,
            isless,
            freeze,
            defrost,
            generation_num,
            population

    # -------

    #isless(lhs::Entity, rhs::Entity) = lhs.fitness < rhs.fitness

    fitness!(ent, fitness_score) = ent.fitness = fitness_score

    distance(ind1 , ind2 ) = abs(ind1.fitness - ind2.fitness)

    global _g_model

    # -------

    function freeze(model::GAmodel, entity::EntityData)
        push!(model.freezer, entity)
        println("Freezing: ", entity)
    end

    function freeze(model::GAmodel, entity)
        entitydata = EntityData(entity, model.params.currentGeneration)
        freeze(model, entitydata)
    end

    freeze(entity) = freeze(_g_model, entity)

    function defrost(model::GAmodel, generation::Int)
        filter(model.freezer) do entitydata
            entitydata.generation == generation
        end
    end

    defrost(generation::Int) = defrost(_g_model, generation)

    generation_num(model::GAmodel = _g_model) = model.params.currentGeneration

    population(model::GAmodel = _g_model) = model.population

    function show_simulation(model :: GAmodel, ent)
        printed = model.specific_fitness.simulate_entity(ent, model.instructionsSet)
        @show ent.fitness
        #@show model.specific_fitness.fitness(ent)
        println(printed)
        return "fitness : $(ent.fitness) \n $printed"
    end

    function runga(mdl::Module, fit_mdl :: Module )
        # get the parameters of the specific fitness
        model = GAmodel(fit_mdl.getParams())
        model.specific_fitness = fit_mdl

        # save the ga (mutation and crossover units function)
        model.ga = mdl
        # save once the dictionary of instructions Set
        model.instructionsSet = mdl.getInstructionsSet()
        #run the ga
        runga(model; resume = false)
    end

    #= function runga(mdl::Module)  #deprecated
        model = GAmodel()
        model.instructionsSet = mdl.getInstructionsSet()
        #model.initial_pop_size = initial_pop_size
        model.ga = mdl
        # TODO, this function need a module for the model.specific_fitness !!
        runga(model; resume = false)
    end=#

    function runga(model::GAmodel; resume = false)
        stop = false

        #set _expandAmount & _expandRate, TODO I think they are not placed very well in case we want specifi paramters for each fitness
        _expandAmount = 2
        _expandRate = 4000

        if (!resume)
            #  initialize params.
            model.params.totalFitness = 0;
            model.params.targetFitness = model.specific_fitness.getTargetFitness();
            model.params.targetFitnessCount = 0;
            model.params.currentGeneration = 0;
            #stop = false;

            reset_model(model)
            create_initial_population(model)
            evaluate_population(model)
        end

        if model.params.historyPath != nothing
            #to save params in the history file
            write(model.params.historyPath, "params : " * string(model.params) * " \n")
        end

        while (model.params.currentGeneration < model.params.generations && !stop)
            #CreatenextGeneration();
            crossover_population(model, [])
            evaluate_population(model)

            lastIdx = model.params.populationSize #length(model.population)
            _fitness = model.population[lastIdx].fitness

            if model.params.currentGeneration % 100 == 0
                _log = "$(Dates.now()) , "
                _log *= "Gen : $(model.params.currentGeneration) , "
                l_1 = lastIdx -1
                _log *= "BEST: $_fitness , "
                _log *= "SECOND: $(model.population[l_1].fitness) \n"
                print(_log)
                if model.params.currentGeneration % 1000 ==0

                    _log *= show_simulation(model, model.population[lastIdx]) * "\n"
                end

                if model.params.historyPath != nothing
                    write(model.params.historyPath, _log)
                end

            end

            #println("hehe $_fitness  $(model.params.targetFitness)")
            if (model.params.targetFitness > 0 && _fitness >= model.params.targetFitness)
                #println("/////////////////////// $(model.params.targetFitnessCount)")
                model.params.targetFitnessCount = model.params.targetFitnessCount +1
                if (model.params.targetFitnessCount > 100) # default : set to 1000
                    break;
                end
            else
                # in case we expand and lose the best one
                model.params.targetFitnessCount = 0;
            end

            #=if (OnGenerationFunction != null)
                OnGenerationFunction(this);
            end=#
            if (_expandAmount > 0 && model.params.currentGeneration > 0 && model.params.currentGeneration % _expandRate == 0 && model.params.genomeSize < model.params.maxGenomeSize )
                model.params.genomeSize +=  _expandAmount;
                #_bestStatus.Fitness = 0; # Update display of best program, since genome has changed and we have a better/worse new best fitness.
            end

            model.params.currentGeneration = model.params.currentGeneration + 1
        end

        if model.params.historyPath != nothing
            best = model.population[model.params.populationSize] #[length(model.population)]
            _log = show_simulation(model, best)
            _log *= "\n Generation : $(model.params.currentGeneration) \n"
            _log *= "fitness $(best.fitness) \n"
            _log *= "dna : $(best.dna) \n"
            write(model.params.historyPath, _log)
            close(model.params.historyPath)
        end
        model
    end

    # -------

    function reset_model(model::GAmodel)
        global _g_model = model

        model.params.currentGeneration = 1
        empty!(model.population)
    end

    function create_initial_population(model::GAmodel)
        for i = 1:model.params.populationSize
            entity = model.ga.create_entity(i, model.params.genomeSize)
            push!(model.population, entity)
        end
    end

    function evaluate_population(model::GAmodel)
        #pmap(model.ga.entityToBfInstructions!, model.population)
        #scores = [ model.specific_fitness.fitness(ent, model.instructionsSet) for ent in model.population ]

        # for each entity, translate dna to bf instructions, then perform the specific fitness,
        # it saves its fitness & bonus in itself and return the fitness, see the fitness function
        scores = pmap(
            ent -> (model.ga.entityToBfInstructions!(ent); model.specific_fitness.fitness(ent, model.instructionsSet))
            , model.population)

        model.params.totalFitness = sum(scores)

        #pmap(fitness!, model.population, scores) # TODO decomment for normal use

        # the isless for the entities is defined in bfga
        sort!(model.population; rev = false)
        model.scores = sort!(scores; rev = false)
    end

    function crossover_population(model::GAmodel, groupings)
        thisGeneration = (population(model))
        _length = 1

        model.population = Any[]

        model.params.currentGeneration += 1
        g = model.ga.create_entity(model.params.currentGeneration, model.params.genomeSize)
        g2 = model.ga.create_entity(model.params.currentGeneration, model.params.genomeSize)
        #println("Generation n' $(model.params.currentGeneration) ")
        if model.params.elitism
            l = model.params.populationSize #length(thisGeneration)
            g = thisGeneration[l]
            l_1 = (model.params.populationSize) -1
            g2 = thisGeneration[l-1]
            #push!(model.population, g)
            #push!(model.pop_data, EntityData(g, model.params.currentGeneration))
            #push!(model.population, g2)
            #push!(model.pop_data, EntityData(g2, model.params.currentGeneration))

            _length += 2;
        end
        #println("_______________________________________________________________________________")
        #@show model.scores

        for i in _length:2:model.params.populationSize
            pidx1 = rouletteSelection(model)
            pidx2 = rouletteSelection(model)
            #println("$pidx1    $pidx2") # for debugging
            child1, child2 = nothing, nothing
            parent1 = thisGeneration[pidx1]
            parent2 = thisGeneration[pidx2]

            if (rand() < model.params.crossoverRate)
                child1, child2 = model.ga.crossover(Any[parent1, parent2])
            else
                child1, child2 = parent1, parent2
                parent1.bonus = 0
                parent2.bonus = 0
            end

            push!(model.population, child1)
            push!(model.population, child2)
        end

        if model.params.elitism
            push!(model.population, g)
            #push!(model.pop_data, EntityData(g, model.params.currentGeneration))
            push!(model.population, g2)
            #push!(model.pop_data, EntityData(g2, model.params.currentGeneration))
        end

        model.params.populationSize = length(model.population)
        #pmap(model.ga.clearCode!, model.population)

        # Expand genomes.
        if model.population[1].m_length != model.params.genomeSize
            #println("expannndd")
            newGenomeSize = model.params.genomeSize
            for m in model.population
                if m.m_length != newGenomeSize
                    model.ga.expand(m, newGenomeSize)
                end
            end
        end
    end

    function mutate_population(model::GAmodel) # deprecated, the mutation is performed in the crossover
        pmap(model.ga.mutate, model.population)
        #=for entity in model.population
            model.ga.mutate(entity)
        end=#
    end

    # We can add strangers in the population
    function add_stranger(model :: GAmodel, num)
        #println("-_-")
        for i in 1:num
            entity = model.ga.create_entity(model.params.currentGeneration, model.population[1].m_length)
            push!(model.population, entity)
            #@async push!(model.pop_data, EntityData(entity, model.params.currentGeneration))
        end
    end

    function rouletteSelection(model :: GAmodel)
        #idx = trunc(Int, rand()*length(model.scores))
        #_idx = idx ==0 ? 1 : idx
        #rand(1:trunc(Int, _idx))

        n = model.params.populationSize #length(model.scores)
        randomFitness = rand() * (model.scores[n] == 0 ? 1 : model.scores[n])

        idx = -1
        first = 1
        last = model.params.populationSize
        mid = trunc(Int , (last - first)/2)

        #  ArrayList's BinarySearch is for exact values only
        #  so do this by hand.
        while (idx == -1 && first <= last)
            if (randomFitness < model.scores[mid])
                last = mid;
            elseif (randomFitness > model.scores[mid])
                first = mid;
            end
            mid = trunc(Int, (first + last)/2)
            #  lies between i and i+1
            if ((last - first) == 1)
                idx = last;
            end
        end

        return idx
    end
end


text:
"""
    random_ones(n, m)

returns a Vector{Int} of size `n` with `m` randomly chosen 1s and
the rest 0s.  TODO: explain how it works?

"""
function random_ones(n::Int, m::Int)
    @assert n > 0 && m <= n && m > 0

    # TODO: implement the m > n/2 case with random for zeros!
    output = zeros(Int, n)

    for pivot = (n-m+1):n
        pick = rand(1:n-m+1)
        if output[pick] == 1
            # pick is already chosen, choose the pivot.
            output[pivot] = 1
        else
            # choose pick.
            output[pick] = 1
        end
    end
    return output
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
export GridSpace

struct Region{D}
    mini::NTuple{D,Int}
    maxi::NTuple{D,Int}
end

"""
    Hood{D}
Internal struct for efficiently finding neighboring positions to a given position.
It contains pre-initialized neighbor cartesian indices and delimiters of when the
neighboring indices would exceed the size of the underlying array.
"""
struct Hood{D} # type P stands for Periodic and is a boolean
    whole::Region{D} # allowed values (only useful for non periodic)
    βs::Vector{CartesianIndex{D}} # neighborhood cartesian indices
end

struct GridSpace{D,P} <: DiscreteSpace
    s::Array{Vector{Int},D}
    metric::Symbol
    hoods::Dict{Float64,Hood{D}}
end

"""
    GridSpace(d::NTuple{D, Int}; periodic = true, metric = :chebyshev)
Create a `GridSpace` that has size given by the tuple `d`, having `D ≥ 1` dimensions.
Optionally decide whether the space will be periodic and what will be the distance metric
used, which decides the behavior of e.g. [`nearby_ids`](@ref).
The position type for this space is `NTuple{D, Int}` and valid positions have indices
in the range `1:d[i]` for the `i`th dimension.

`:chebyshev` metric means that the `r`-neighborhood of a position are all
positions within the hypercube having side length of `2*floor(r)` and being centered in
the origin position.

`:euclidean` metric means that the `r`-neighborhood of a position are all positions whose
cartesian indices have Euclidean distance `≤ r` from the cartesian index of the given
position.
"""
function GridSpace(
        d::NTuple{D,Int};
        periodic::Bool = true,
        metric::Symbol = :chebyshev,
        moore = nothing
    ) where {D}
    s = Array{Vector{Int},D}(undef, d)
    if moore ≠ nothing
        @warn "Keyword `moore` is deprecated, use `metric` instead."
        metric = moore == true ? :chebyshev : :euclidean
    end
    for i in eachindex(s)
        s[i] = Int[]
    end
    return GridSpace{D,periodic}(s, metric, Dict{Float64,Hood{D}}())
end

#######################################################################################
# %% Implementation of space API
#######################################################################################
function random_position(model::ABM{<:AbstractAgent,<:GridSpace})
    Tuple(rand(CartesianIndices(model.space.s)))
end

function add_agent_to_space!(a::A, model::ABM{A,<:GridSpace}) where {A<:AbstractAgent}
    push!(model.space.s[a.pos...], a.id)
    return a
end

function remove_agent_from_space!(a::A, model::ABM{A,<:GridSpace}) where {A<:AbstractAgent}
    prev = model.space.s[a.pos...]
    ai = findfirst(i -> i == a.id, prev)
    deleteat!(prev, ai)
    return a
end

function move_agent!(a::A, pos::Tuple, model::ABM{A,<:GridSpace}) where {A<:AbstractAgent}
    remove_agent_from_space!(a, model)
    a.pos = pos
    add_agent_to_space!(a, model)
end

##########################################################################################
# %% Structures for collecting neighbors on a grid
##########################################################################################
# Most of the source code in this section comes from TimeseriesPrediction.jl, specifically
# the file github.com/JuliaDynamics/TimeseriesPrediction.jl/src/spatiotemporalembedding.jl
# It creates a performant envinroment where the cartesian indices of a given neighborhood
# of given radious and metric type are stored and re-used for each search.

Base.length(r::Region{D}) where {D} = prod(r.maxi .- r.mini .+ 1)
function Base.in(idx, r::Region{D}) where {D}
    @inbounds for φ in 1:D
        r.mini[φ] <= idx[φ] <= r.maxi[φ] || return false
    end
    return true
end

# This function initializes the standard cartesian indices that needs to be added to some
# index `α` to obtain its neighborhood
import LinearAlgebra
function initialize_neighborhood!(space::GridSpace{D}, r::Real) where {D}
    d = size(space.s)
    r0 = floor(Int, r)
    if space.metric == :euclidean
        # hypercube of indices
        hypercube = CartesianIndices((repeat([(-r0):r0], D)...,))
        # select subset of hc which is in Hypersphere
        βs = [β for β ∈ hypercube if LinearAlgebra.norm(β.I) ≤ r]
    elseif space.metric == :chebyshev
        βs = vec([CartesianIndex(a) for a in Iterators.product([(-r0):r0 for φ in 1:D]...)])
    else
        error("Unknown metric type")
    end
    whole = Region(map(one, d), d)
    hood = Hood{D}(whole, βs)
    space.hoods[float(r)] = hood
    return hood
end

"""
    grid_space_neighborhood(α::CartesianIndex, space::GridSpace, r::Real)

Return an iterator over all positions within distance `r` from `α` according to the `space`.

The only reason this function is not equivalent with `nearby_positions` is because
`nearby_positions` skips the current position `α`, while `α` is always included in the
returned iterator (because the `0` index is always included in `βs`).

This function re-uses the source code of TimeseriesPrediction.jl, along with the
helper struct `Hood` and generates neighboring cartesian indices on the fly,
reducing the amount of computations necessary (i.e. we don't "find" new indices,
we only add a pre-determined amount of indices to `α`).
"""
function grid_space_neighborhood(α::CartesianIndex, space::GridSpace, r::Real)
    hood = if haskey(space.hoods, r)
        space.hoods[r]
    else
        initialize_neighborhood!(space, r)
    end
    _grid_space_neighborhood(α, space, hood)
end

function _grid_space_neighborhood(
        α::CartesianIndex,
        space::GridSpace{D,true},
        hood,
    ) where {D}
    return ((mod1.(Tuple(α + β), hood.whole.maxi)) for β in hood.βs)
end

function _grid_space_neighborhood(
        α::CartesianIndex,
        space::GridSpace{D,false},
        hood,
    ) where {D}
    return Iterators.filter(x -> x ∈ hood.whole, (Tuple(α + β) for β in hood.βs))
end

grid_space_neighborhood(α, model::ABM, r) = grid_space_neighborhood(α, model.space, r)

##########################################################################################
# %% Extend neighbors API for spaces
##########################################################################################
function nearby_ids(pos::ValidPos, model::ABM{<:AbstractAgent,<:GridSpace}, r = 1)
    nn = grid_space_neighborhood(CartesianIndex(pos), model, r)
    s = model.space.s
    Iterators.flatten((s[i...] for i in nn))
end

function nearby_positions(pos::ValidPos, model::ABM{<:AbstractAgent,<:GridSpace}, r = 1)
    nn = grid_space_neighborhood(CartesianIndex(pos), model, r)
    Iterators.filter(!isequal(pos), nn)
end

#######################################################################################
# %% Further discrete space functions
#######################################################################################
function positions(model::ABM{<:AbstractAgent,<:GridSpace})
    x = CartesianIndices(model.space.s)
    return (Tuple(y) for y in x)
end

function ids_in_position(pos::ValidPos, model::ABM{<:AbstractAgent,<:GridSpace})
    return model.space.s[pos...]
end

###################################################################
# %% pretty printing
###################################################################
Base.size(space::GridSpace) = size(space.s)
function Base.show(io::IO, space::GridSpace{D, P}) where {D, P}
    s = "GridSpace with size $(size(space)), metric=$(space.metric) and periodic=$(P)"
    print(io, s)
end


text:
function installPumi()
  #pumi_version = "e4eabf5"
  #pumi_version = "d1837c5936c28d6ef09abd02c82ea2a4ea9b6f55"
#  pumi_version = "7efdf2319d13ce86994930e867d226ea421c1717"
#  pumi_version = "63fa247358e46352c0d500b1f1735c82447dfddd"
#   pumi_version = "develop"
#  pumi_version = "063ad65"
#  pumi_version = "8b7c260"
  pumi_version="7ba0ce4"
#  pumi_version = "HEAD"

  start_dir = pwd()  # record where we started
  run(`./cleanup.sh`)
#  if isdir("./core")
#    println("deleting existing Core repo in /deps")
#    rm("./core", recursive=true)
#  end

  # download files (if tarballs not present)
  run(`./download.sh`)

  # extract files
  if isdir("core")
    try
      run(`rm -rf ./core`)
    end
  end
  run(`tar xfvz ./core.tar.gz`)

  if isdir("pumi-meshes")
    try
      run(`rm -rf ./pumi-meshes`)
    end
  end
  run(`tar xfvz ./pumi-meshes.tar.gz`)

  # config, make, make install
  cd("./core")
#  run(`git pull`)
  run(`git checkout $pumi_version`)
  mkdir("./build")
  cd("./build")
  mkdir("./install")  # install directory
  run(`../../config.sh`)
  run(`make -j 4`)
  run(`make install`)
  str1 = joinpath( pwd(), "install/lib")
  str3 = joinpath(str1, "pkgconfig")

  # update ENV in preparation for building files in /src
  if haskey(ENV, "LD_LIBRARY_PATH")
    ld_path = ENV["LD_LIBRARY_PATH"]
    ld_path = string(str1, ":", ld_path)
    ENV["LD_LIBRARY_PATH"] = ld_path
  else
    ENV["LD_LIBRARY_PATH"] = str1
  end

  # for backwards compatibility with the pkg-config days
  if haskey(ENV, "PKG_CONFIG_PATH")
    pkg_path = ENV["PKG_CONFIG_PATH"]
    pkg_path = string(str3,":", pkg_path)
    ENV["PKG_CONFIG_PATH"] = pkg_path
  else
    ENV["PKG_CONFIG_PATH"] = str3
  end

  scorec_prefix = joinpath(pwd(), "install")
  println("scorec prefix = ", scorec_prefix)
  ENV["SCOREC_PREFIX"] = scorec_prefix

  cd(start_dir)
  return nothing
end


text:
module TestDatasource

using Base.Test
using SimpleTasksTests.Utils.MockServices
using SimpleTasks.Services.BucketCacheDatasource

import SimpleTasks.Services.Bucket
import SimpleTasks.Services.Cache
import SimpleTasks.Services.Datasource

function test_get_empty_cache()
    key = "somekey"

    bucket_text = "mock contents"
    bucket_file = IOBuffer(bucket_text)
    seekstart(bucket_file)
    bucket_files = Dict()
    bucket_files[key] = bucket_file

    cache_values = Dict()

    bucket = MockBucketService(bucket_files)
    cache = MockCacheService(cache_values)

    datasource = BucketCacheDatasourceService(bucket, cache)

    get_value = Datasource.get(datasource, key)

    @test haskey(cache.mockValues, key)
    new_cache_value = cache.mockValues[key]
    # we are returned the correct value
    @test new_cache_value == get_value
    # cache gets updated with the bucket text
    @test readchomp(new_cache_value) == bucket_text
end

function test_get_multi_empty_cache()
    key = "somekey"
    bucket_text = "mock contents"
    bucket_file = IOBuffer(bucket_text)
    seekstart(bucket_file)

    key2 = "somekey2"
    bucket_text2 = "mock contents2"
    bucket_file2 = IOBuffer(bucket_text2)
    seekstart(bucket_file2)

    bucket_files = Dict()
    bucket_files[key] = bucket_file
    bucket_files[key2] = bucket_file2

    cache_values = Dict()

    bucket = MockBucketService(bucket_files)
    cache = MockCacheService(cache_values)

    datasource = BucketCacheDatasourceService(bucket, cache)

    (get_value1, get_value2) = Datasource.get(datasource, [key, key2])

    @test haskey(cache.mockValues, key)
    new_cache_value = cache.mockValues[key]
    # we are returned the correct value
    @test get_value1 == new_cache_value
    # cache gets updated with the bucket text
    @test readchomp(new_cache_value) == bucket_text

    @test haskey(cache.mockValues, key2)
    new_cache_value2 = cache.mockValues[key2]
    # we are returned the correct value
    @test get_value2 == new_cache_value2
    # cache gets updated with the bucket text
    @test readchomp(new_cache_value2) == bucket_text2
end

function test_get_with_cache()
    key = "somekey"

    bucket_text = "mock contents"
    bucket_file = IOBuffer(bucket_text)
    seekstart(bucket_file)
    bucket_files = Dict()
    bucket_files[key] = bucket_file

    cache_text = "mock contents cached already"
    cache_io = IOBuffer(cache_text)
    seekstart(cache_io)
    cache_values = Dict()
    cache_values[key] = cache_io

    bucket = MockBucketService(bucket_files)
    cache = MockCacheService(cache_values)
    datasource = BucketCacheDatasourceService(bucket, cache)

    get_value = Datasource.get(datasource, key)

    @test haskey(cache.mockValues, key)
    new_cache_io = cache.mockValues[key]
    # we are returned the correct value
    @test get_value == new_cache_io
    # cache does not get updated because it is already in there
    @test readchomp(new_cache_io) == cache_text
end

function test_get_no_bucket()
    key = "somekey"

    bucket_files = Dict()
    cache_values = Dict()

    bucket = MockBucketService(bucket_files)
    cache = MockCacheService(cache_values)
    datasource = BucketCacheDatasourceService(bucket, cache)

    @test_throws Exception Datasource.get(datasource, key)
end

function test_override_cache_get_empty_cache()
    key = "somekey"

    bucket_text = "mock contents"
    bucket_file = IOBuffer(bucket_text)
    seekstart(bucket_file)
    bucket_files = Dict()
    bucket_files[key] = bucket_file

    cache_values = Dict()

    bucket = MockBucketService(bucket_files)
    cache = MockCacheService(cache_values)

    datasource = BucketCacheDatasourceService(bucket, cache)

    get_value = Datasource.get(datasource, key; override_cache=true)

    @test haskey(cache.mockValues, key)
    new_cache_value = cache.mockValues[key]
    # we are returned the correct value
    @test get_value == new_cache_value
    # cache gets forced to update with the new bucket text
    @test readchomp(new_cache_value) == bucket_text
end

function test_override_cache_get_with_cache()
    key = "somekey"

    bucket_text = "mock contents"
    bucket_file = IOBuffer(bucket_text)
    seekstart(bucket_file)
    bucket_files = Dict()
    bucket_files[key] = bucket_file

    cache_text = "mock contents cached already"
    cache_io = IOBuffer(cache_text)
    seekstart(cache_io)
    cache_values = Dict()
    cache_values[key] = cache_io

    bucket = MockBucketService(bucket_files)
    cache = MockCacheService(cache_values)
    datasource = BucketCacheDatasourceService(bucket, cache)

    get_value = Datasource.get(datasource, key; override_cache=true)

    @test haskey(cache.mockValues, key)
    new_cache_io = cache.mockValues[key]
    # we are returned the correct value
    @test get_value == new_cache_io
    # cache gets forced to update with the new bucket text
    @test readchomp(new_cache_io) == bucket_text
end

function test_put_no_value()
    key = "somekey"

    bucket_text = "mock contents"
    bucket_file = IOBuffer(bucket_text)
    seekstart(bucket_file)
    bucket_files = Dict()
    bucket_files[key] = bucket_file

    cache_text = "mock contents cached already"
    cache_io = IOBuffer(cache_text)
    seekstart(cache_io)
    cache_values = Dict()
    cache_values[key] = cache_io

    bucket = MockBucketService(bucket_files)
    cache = MockCacheService(cache_values)
    datasource = BucketCacheDatasourceService(bucket, cache)

    result = Datasource.put!(datasource, key)

    @test result

    @test haskey(bucket.mockFiles, key)
    new_bucket_file = bucket.mockFiles[key]
    # putting back to bucket with new cached file updates bucket
    @test readchomp(new_bucket_file) == cache_text
end

function test_put_multi_no_value()
    key1 = "somekey1"

    bucket_text1 = "mock contents"
    bucket_file1 = IOBuffer(bucket_text1)
    seekstart(bucket_file1)

    key2 = "somekey2"

    bucket_text2 = "mock contents2"
    bucket_file2 = IOBuffer(bucket_text2)
    seekstart(bucket_file2)

    bucket_files = Dict()
    bucket_files[key1] = bucket_file1
    bucket_files[key2] = bucket_file2

    cache_text1 = "mock contents cached already1"
    cache_io1 = IOBuffer(cache_text1)
    seekstart(cache_io1)

    cache_text2 = "mock contents cached already2"
    cache_io2 = IOBuffer(cache_text2)
    seekstart(cache_io2)

    cache_values = Dict()
    cache_values[key1] = cache_io1
    cache_values[key2] = cache_io2

    bucket = MockBucketService(bucket_files)
    cache = MockCacheService(cache_values)
    datasource = BucketCacheDatasourceService(bucket, cache)

    (result1, result2) = Datasource.put!(datasource, [key1, key2])

    @test result1
    @test result2

    @test haskey(bucket.mockFiles, key1)
    new_bucket_file1 = bucket.mockFiles[key1]
    # putting back to bucket with new cached file updates bucket
    @test readchomp(new_bucket_file1) == cache_text1

    @test haskey(bucket.mockFiles, key2)
    new_bucket_file2 = bucket.mockFiles[key2]
    # putting back to bucket with new cached file updates bucket
    @test readchomp(new_bucket_file2) == cache_text2
end

function test_put_not_exist_no_value()
    key = "somekey"

    bucket_text = "mock contents"
    bucket_file = IOBuffer(bucket_text)
    seekstart(bucket_file)
    bucket_files = Dict()
    bucket_files[key] = bucket_file

    cache_values = Dict()

    bucket = MockBucketService(bucket_files)
    cache = MockCacheService(cache_values)
    datasource = BucketCacheDatasourceService(bucket, cache)

    result = Datasource.put!(datasource, key)

    @test !result

    @test haskey(bucket.mockFiles, key)
    new_bucket_file = bucket.mockFiles[key]
    # putting back to bucket when not existing in cache does not modify bucket
    @test readchomp(new_bucket_file) == bucket_text
end

function test_put_new_value()
    key = "somekey"

    bucket_text = "mock contents"
    bucket_file = IOBuffer(bucket_text)
    seekstart(bucket_file)
    bucket_files = Dict()
    bucket_files[key] = bucket_file

    cache_text = "mock contents cached already"
    cache_io = IOBuffer(cache_text)
    seekstart(cache_io)
    cache_values = Dict()
    cache_values[key] = cache_io

    bucket = MockBucketService(bucket_files)
    cache = MockCacheService(cache_values)
    datasource = BucketCacheDatasourceService(bucket, cache)

    new_text = "new value"
    new_value = IOBuffer(new_text)
    seekstart(new_value)
    result = Datasource.put!(datasource, key, new_value)

    @test result == true

    @test haskey(bucket.mockFiles, key)
    new_bucket_file = bucket.mockFiles[key]
    seekstart(new_bucket_file)
    # putting with new value changes both bucket and cache
    @test readchomp(new_bucket_file) == new_text

    @test haskey(cache.mockValues, key)
    new_cache_value = cache.mockValues[key]
    seekstart(new_cache_value)
    # putting with new value changes both bucket and cache
    @test readchomp(new_cache_value) == new_text
end

function test_put_multi_new_value()
    key1 = "somekey1"

    bucket_text1 = "mock contents"
    bucket_file1 = IOBuffer(bucket_text1)
    seekstart(bucket_file1)

    key2 = "somekey2"

    bucket_text2 = "mock contents2"
    bucket_file2 = IOBuffer(bucket_text2)
    seekstart(bucket_file2)

    bucket_files = Dict()
    bucket_files[key1] = bucket_file1
    bucket_files[key2] = bucket_file2

    cache_text1 = "mock contents cached already1"
    cache_io1 = IOBuffer(cache_text1)
    seekstart(cache_io1)

    cache_text2 = "mock contents cached already2"
    cache_io2 = IOBuffer(cache_text2)
    seekstart(cache_io2)

    cache_values = Dict()
    cache_values[key1] = cache_io1
    cache_values[key2] = cache_io2

    bucket = MockBucketService(bucket_files)
    cache = MockCacheService(cache_values)
    datasource = BucketCacheDatasourceService(bucket, cache)

    new_text1 = "new value1"
    new_value1 = IOBuffer(new_text1)
    seekstart(new_value1)

    new_text2 = "new value2"
    new_value2 = IOBuffer(new_text2)
    seekstart(new_value2)

    (result1, result2) = Datasource.put!(datasource, [key1, key2],
        [new_value1, new_value2])

    @test result1
    @test result2

    @test haskey(bucket.mockFiles, key1)
    new_bucket_file1 = bucket.mockFiles[key1]
    seekstart(new_bucket_file1)
    # putting with new value changes both bucket and cache
    @test readchomp(new_bucket_file1) == new_text1

    @test haskey(cache.mockValues, key1)
    new_cache_io1 = cache.mockValues[key1]
    seekstart(new_cache_io1)
    # putting with new value changes both bucket and cache
    @test readchomp(new_cache_io1) == new_text1

    @test haskey(bucket.mockFiles, key2)
    new_bucket_file2 = bucket.mockFiles[key2]
    seekstart(new_bucket_file2)
    # putting with new value changes both bucket and cache
    @test readchomp(new_bucket_file2) == new_text2

    @test haskey(cache.mockValues, key2)
    new_cache_io2 = cache.mockValues[key2]
    seekstart(new_cache_io2)
    # putting with new value changes both bucket and cache
    @test readchomp(new_cache_io2) == new_text2
end

function test_put_not_exist_new_value()
    key = "somekey"

    bucket_text = "mock contents"
    bucket_file = IOBuffer(bucket_text)
    seekstart(bucket_file)
    bucket_files = Dict()
    bucket_files[key] = bucket_file

    cache_values = Dict()

    bucket = MockBucketService(bucket_files)
    cache = MockCacheService(cache_values)
    datasource = BucketCacheDatasourceService(bucket, cache)

    new_text = "new value"
    new_value = IOBuffer(new_text)
    seekstart(new_value)
    result = Datasource.put!(datasource, key, new_value)

    @test result

    @test haskey(bucket.mockFiles, key)
    new_bucket_file = bucket.mockFiles[key]
    # even if the value doesn't exist in cache, the new_text should be used
    @test readchomp(new_bucket_file) == new_text
end

function test_put_multi_new_value_only_cache()
    key1 = "somekey1"

    bucket_text1 = "mock contents"
    bucket_file1 = IOBuffer(bucket_text1)
    seekstart(bucket_file1)

    key2 = "somekey2"

    bucket_text2 = "mock contents2"
    bucket_file2 = IOBuffer(bucket_text2)
    seekstart(bucket_file2)

    bucket_files = Dict()
    bucket_files[key1] = bucket_file1
    bucket_files[key2] = bucket_file2

    cache_text1 = "mock contents cached already1"
    cache_io1 = IOBuffer(cache_text1)
    seekstart(cache_io1)

    cache_text2 = "mock contents cached already2"
    cache_io2 = IOBuffer(cache_text2)
    seekstart(cache_io2)

    cache_values = Dict()
    cache_values[key1] = cache_io1
    cache_values[key2] = cache_io2

    bucket = MockBucketService(bucket_files)
    cache = MockCacheService(cache_values)
    datasource = BucketCacheDatasourceService(bucket, cache)

    new_text1 = "new value1"
    new_value1 = IOBuffer(new_text1)
    seekstart(new_value1)

    new_text2 = "new value2"
    new_value2 = IOBuffer(new_text2)
    seekstart(new_value2)

    (result1, result2) = Datasource.put!(datasource, [key1, key2],
        [new_value1, new_value2]; only_cache=true)

    @test result1
    @test result2

    @test haskey(bucket.mockFiles, key1)
    new_bucket_file1 = bucket.mockFiles[key1]
    new_cache_io1 = cache.mockValues[key1]
    # putting with new value changes both bucket and cache
    @test readchomp(new_bucket_file1) == bucket_text1
    @test readchomp(new_cache_io1) == new_text1

    @test haskey(bucket.mockFiles, key2)
    new_bucket_file2 = bucket.mockFiles[key2]
    new_cache_io2 = cache.mockValues[key2]
    # putting with new value changes both bucket and cache
    @test readchomp(new_bucket_file2) == bucket_text2
    @test readchomp(new_cache_io2) == new_text2

end

function test_delete()
    key = "somekey"

    bucket_text = "mock contents"
    bucket_file = IOBuffer(bucket_text)
    seekstart(bucket_file)
    bucket_files = Dict()
    bucket_files[key] = bucket_file

    cache_text = "mock contents cached already"
    cache_io = IOBuffer(cache_text)
    seekstart(cache_io)
    cache_values = Dict()
    cache_values[key] = cache_io

    bucket = MockBucketService(bucket_files)
    cache = MockCacheService(cache_values)
    datasource = BucketCacheDatasourceService(bucket, cache)

    Datasource.delete!(datasource, key)

    # ensure that the value was deleted from both bucket and cache
    @test !haskey(bucket.mockFiles, key)
    @test !haskey(cache.mockValues, key)
end

function test_delete_multi()
    key1 = "somekey1"

    bucket_text1 = "mock contents"
    bucket_file1 = IOBuffer(bucket_text1)
    seekstart(bucket_file1)

    key2 = "somekey2"

    bucket_text2 = "mock contents2"
    bucket_file2 = IOBuffer(bucket_text2)
    seekstart(bucket_file2)

    bucket_files = Dict()
    bucket_files[key1] = bucket_file1
    bucket_files[key2] = bucket_file2

    cache_text1 = "mock contents cached already1"
    cache_io1 = IOBuffer(cache_text1)
    seekstart(cache_io1)

    cache_text2 = "mock contents cached already2"
    cache_io2 = IOBuffer(cache_text2)
    seekstart(cache_io2)

    cache_values = Dict()
    cache_values[key1] = cache_io1
    cache_values[key2] = cache_io2

    bucket = MockBucketService(bucket_files)
    cache = MockCacheService(cache_values)
    datasource = BucketCacheDatasourceService(bucket, cache)

    Datasource.delete!(datasource, [key1, key2])

    # ensure that all values were deleted from both bucket and cache
    @test !haskey(bucket.mockFiles, key1)
    @test !haskey(cache.mockValues, key1)

    @test !haskey(bucket.mockFiles, key2)
    @test !haskey(cache.mockValues, key2)
end

function test_delete_cache_only()
    key = "somekey"

    bucket_text = "mock contents"
    bucket_file = IOBuffer(bucket_text)
    seekstart(bucket_file)
    bucket_files = Dict()
    bucket_files[key] = bucket_file

    cache_text = "mock contents cached already"
    cache_io = IOBuffer(cache_text)
    seekstart(cache_io)
    cache_values = Dict()
    cache_values[key] = cache_io

    bucket = MockBucketService(bucket_files)
    cache = MockCacheService(cache_values)
    datasource = BucketCacheDatasourceService(bucket, cache)

    Datasource.delete!(datasource, key; only_cache=true)

    # ensure that the value was deleted from both bucket and cache
    @test haskey(bucket.mockFiles, key)
    @test !haskey(cache.mockValues, key)
end

function test_clear_cache()
    key1 = "somekey1"
    key2 = "somekey2"

    bucket_text1 = "mock contents"
    bucket_file1 = IOBuffer(bucket_text1)
    seekstart(bucket_file1)
    bucket_files = Dict()
    bucket_files[key1] = bucket_file1

    cache_text1 = "mock contents cached already"
    cache_io1 = IOBuffer(cache_text1)
    seekstart(cache_io1)

    cache_text2 = "mock contents cached already"
    cache_io2 = IOBuffer(cache_text2)
    seekstart(cache_io2)

    cache_values = Dict()
    cache_values[key1] = cache_io1
    cache_values[key2] = cache_io2

    bucket = MockBucketService(bucket_files)
    cache = MockCacheService(cache_values)
    datasource = BucketCacheDatasourceService(bucket, cache)

    Datasource.clear_cache(datasource)

    # ensure that the value was deleted only from cache
    @test haskey(bucket.mockFiles, key1)

    @test !haskey(cache.mockValues, key1)
    @test !haskey(cache.mockValues, key2)
end

function __init__()
    test_get_empty_cache()
    test_get_multi_empty_cache()
    test_get_with_cache()
    test_get_no_bucket()
    test_override_cache_get_empty_cache()
    test_override_cache_get_with_cache()

    test_put_no_value()
    test_put_multi_no_value()
    test_put_not_exist_no_value()

    test_put_new_value()
    test_put_multi_new_value()
    test_put_not_exist_new_value()

    test_delete()
    test_delete_multi()
    test_delete_cache_only()

    test_clear_cache()
end

end # module TestDatasource

max_forks_repo_name: seung-lab/SimpleTasks


text:
import JuliaDB: select, groupby, filter
function Base.getindex(crps::Union{CorpusRaw,TanzilRaw}, i::Union{Int64,Array{Int64,1},UnitRange{Int64}})
    if i isa Int64
        1 <= i <= length(crps.data) || throw(BoundsError(crps, i))
        return crps.data[i]
    elseif i isa Array{Int64}
        1 <= maximum(i) <= length(crps.data) || throw(BoundsError(crps, maximum(i)))
        return crps.data[i]
    else
        1 <= i.stop <= length(crps.data) || throw(BoundsError(crps, i))
        return crps.data[i.start:i.stop]
    end
end
Base.firstindex(crps::Union{CorpusRaw,TanzilRaw}) = 1
Base.lastindex(crps::Union{CorpusRaw,TanzilRaw}) = length(crps.data)

function Base.getindex(crps::Union{CorpusData,TanzilData}, i::Union{Int64,Array{Int64,1},UnitRange{Int64}})
    if i isa Int64
        1 <= i <= 114 || throw(BoundsError(crps, i))
        tbl = filter(t -> t.chapter === i, crps.data)
        verses = select(tbl, :verse)
        if crps isa TanzilData
            return Chapter(tbl, i, verses, true)
        else
            return Chapter(tbl, i, verses)
        end
    else
        1 <= maximum(i) <= 114 || throw(BoundsError(crps, maximum(i)))
        tbl = filter(t -> in(t.chapter, i), crps.data)
        vrstbl = groupby((verse = :verse => x -> unique(x)[end],), tbl, :chapter)
        verses = select(vrstbl, :verse)
        if crps isa TanzilData
            return Chapter(tbl, i, verses, true)
        else
            return Chapter(tbl, i, verses)
        end
    end
end
Base.firstindex(crps::Union{CorpusData,TanzilData}) = 1
Base.lastindex(crps::Union{CorpusData,TanzilData}) = 114

function Base.getindex(chapter::Chapter, i::Union{Int64,Array{Int64,1},UnitRange{Int64},Tuple})
    if i isa Int64
        1 <= i <= maximum(unique(chapter.verses)) || throw(BoundsError(chapter, i))
        tbl = filter(t -> t.verse === i, chapter.data)
        try
            if chapter.numbers isa Int64
                words = select(tbl, :word)
                return Verse(tbl, i, chapter.numbers, words, chapter.istanzil)
            else
                wrdtbl = groupby((word = :word => x -> unique(x)[end],), tbl, :chapter)
                words = select(wrdtbl, :word)
                return Verse(tbl, i, chapter.numbers, words, chapter.istanzil)
            end
        catch
            words = missing
            return Verse(tbl, i, chapter.numbers, words, chapter.istanzil)
        end
    else
        1 <= maximum(i) <= maximum(unique(chapter.verses)) || throw(BoundsError(chapter, maximum(i)))
        tbl = filter(t -> in(t.verse, i), chapter.data)
        try
            if chapter.numbers isa Int64
                wrdtbl = groupby((word = :word => x -> unique(x)[end],), tbl, :verse)
            else
                wrdtbl = groupby((word = :word => x -> unique(x)[end],), tbl, (:chapter, :verse))
            end
            words = select(wrdtbl, :word)
            return Verse(tbl, i, chapter.numbers, words, chapter.istanzil)
        catch
            words = missing
            return Verse(tbl, i, chapter.numbers, words, chapter.istanzil)
        end
    end
end

function Base.getindex(verse::Verse, i::Union{Int64,Array{Int64,1},UnitRange{Int64},Tuple})
    !verse.istanzil || throw(DomainError("TanzilData: No column word to index."))

    if i isa Int64
        1 <= i <= maximum(unique(verse.words)) || throw(BoundsError(verse, i))
        tbl = filter(t -> t.word === i, verse.data)
        if verse.numbers isa Int64
            if verse.chapters isa Int64
                parts = select(tbl, :part)
                return Word(tbl, i, verse.chapters, verse.numbers, parts)
            else
                prttbl = groupby((part = :part => x -> unique(x)[end],), tbl, :chapter)
                parts = select(prttbl, :part)
                return Word(tbl, i, verse.chapters, verse.numbers, parts)
            end
        else
            if verse.chapters isa Int64
                prttbl = groupby((part = :part => x -> unique(x)[end],), tbl, :verse)
            else
                prttbl = groupby((part = :part => x -> unique(x)[end],), tbl, (:chapter, :verse, :word))
            end
            parts = select(prttbl, :part)
            return Word(tbl, i, verse.chapters, verse.numbers, parts)
        end
    else
        1 <= maximum(i) <= maximum(unique(verse.words)) || throw(BoundsError(verse, maximum(i)))
        tbl = filter(t -> in(t.word, i), verse.data)
        if verse.numbers isa Int64
            if verse.chapters isa Int64
                prttbl = groupby((part = :part => x -> unique(x)[end],), tbl, (:chapter, :verse, :word))
            else
                prttbl = groupby((part = :part => x -> unique(x)[end],), tbl, (:chapter, :verse, :word))
            end
        else
            if verse.chapters isa Int64
                prttbl = groupby((part = :part => x -> unique(x)[end],), tbl, (:verse, :word))
            else
                prttbl = groupby((part = :part => x -> unique(x)[end],), tbl, (:chapter, :verse, :word))
            end
        end
        parts = select(prttbl, :part)
        return Word(tbl, i, verse.chapters, verse.numbers, parts)
    end
end

function Base.getindex(word::Word, i::Union{Int64,Array{Int64,1},UnitRange{Int64},Tuple})
    if i isa Int64
        1 <= i <= unique(word.parts)[end] || throw(BoundsError(word, i))
        tbl = filter(t -> t.part === i, word.data)
        return Part(tbl, i, word.chapters, word.verses, word.numbers)
    else
        1 <= maximum(i) <= maximum(word.parts) || throw(BoundsError(word, maximum(i)))
        tbl = filter(t -> in(t.part, i), word.data)
        return Part(tbl, i, word.chapters, word.verses, word.numbers)
    end
end


text:
using Mocking

associate_vpcwith_hosted_zone_1_patch = @patch function associate_vpcwith_hosted_zone(a...; b...)
    return Dict{String, Any}("ChangeInfo" => Dict{String, Any}("Id" => "/change/C3HC6WDB2UANE2", "Comment" => "", "Status" => "INSYNC", "SubmittedAt" => "2017-01-31T01:36:41.958Z"))
end

change_resource_record_sets_10_patch = @patch function change_resource_record_sets(a...; b...)
    return Dict{String, Any}("ChangeInfo" => Dict{String, Any}("Id" => "/change/C2682N5HXP0BZ4", "Comment" => "Geolocation alias configuration for example.com", "Status" => "PENDING", "SubmittedAt" => "2017-02-10T01:36:41.958Z"))
end

change_resource_record_sets_1_patch = @patch function change_resource_record_sets(a...; b...)
    return Dict{String, Any}("ChangeInfo" => Dict{String, Any}("Id" => "/change/C2682N5HXP0BZ4", "Comment" => "Web server for example.com", "Status" => "PENDING", "SubmittedAt" => "2017-02-10T01:36:41.958Z"))
end

change_resource_record_sets_2_patch = @patch function change_resource_record_sets(a...; b...)
    return Dict{String, Any}("ChangeInfo" => Dict{String, Any}("Id" => "/change/C2682N5HXP0BZ4", "Comment" => "Web servers for example.com", "Status" => "PENDING", "SubmittedAt" => "2017-02-10T01:36:41.958Z"))
end

change_resource_record_sets_3_patch = @patch function change_resource_record_sets(a...; b...)
    return Dict{String, Any}("ChangeInfo" => Dict{String, Any}("Id" => "/change/C2682N5HXP0BZ4", "Comment" => "CloudFront distribution for example.com", "Status" => "PENDING", "SubmittedAt" => "2017-02-10T01:36:41.958Z"))
end

change_resource_record_sets_4_patch = @patch function change_resource_record_sets(a...; b...)
    return Dict{String, Any}("ChangeInfo" => Dict{String, Any}("Id" => "/change/C2682N5HXP0BZ4", "Comment" => "ELB load balancers for example.com", "Status" => "PENDING", "SubmittedAt" => "2017-02-10T01:36:41.958Z"))
end

change_resource_record_sets_5_patch = @patch function change_resource_record_sets(a...; b...)
    return Dict{String, Any}("ChangeInfo" => Dict{String, Any}("Id" => "/change/C2682N5HXP0BZ4", "Comment" => "EC2 instances for example.com", "Status" => "PENDING", "SubmittedAt" => "2017-02-10T01:36:41.958Z"))
end

change_resource_record_sets_6_patch = @patch function change_resource_record_sets(a...; b...)
    return Dict{String, Any}("ChangeInfo" => Dict{String, Any}("Id" => "/change/C2682N5HXP0BZ4", "Comment" => "ELB load balancers for example.com", "Status" => "PENDING", "SubmittedAt" => "2017-02-10T01:36:41.958Z"))
end

change_resource_record_sets_7_patch = @patch function change_resource_record_sets(a...; b...)
    return Dict{String, Any}("ChangeInfo" => Dict{String, Any}("Id" => "/change/C2682N5HXP0BZ4", "Comment" => "Failover configuration for example.com", "Status" => "PENDING", "SubmittedAt" => "2017-02-10T01:36:41.958Z"))
end

change_resource_record_sets_8_patch = @patch function change_resource_record_sets(a...; b...)
    return Dict{String, Any}("ChangeInfo" => Dict{String, Any}("Id" => "/change/C2682N5HXP0BZ4", "Comment" => "Failover alias configuration for example.com", "Status" => "PENDING", "SubmittedAt" => "2017-02-10T01:36:41.958Z"))
end

change_resource_record_sets_9_patch = @patch function change_resource_record_sets(a...; b...)
    return Dict{String, Any}("ChangeInfo" => Dict{String, Any}("Id" => "/change/C2682N5HXP0BZ4", "Comment" => "Geolocation configuration for example.com", "Status" => "PENDING", "SubmittedAt" => "2017-02-10T01:36:41.958Z"))
end

change_tags_for_resource_1_patch = @patch function change_tags_for_resource(a...; b...)
    return Dict{String, Any}()
end

get_hosted_zone_1_patch = @patch function get_hosted_zone(a...; b...)
    return Dict{String, Any}("HostedZone" => Dict{String, Any}("Id" => "/hostedzone/Z3M3LMPEXAMPLE", "ResourceRecordSetCount" => 8, "CallerReference" => "C741617D-04E4-F8DE-B9D7-0D150FC61C2E", "Config" => Dict{String, Any}("PrivateZone" => false), "Name" => "myawsbucket.com."), "DelegationSet" => Dict{String, Any}("NameServers" => Any["ns-2048.awsdns-64.com", "ns-2049.awsdns-65.net", "ns-2050.awsdns-66.org", "ns-2051.awsdns-67.co.uk"]))
end


text:
function normalize_ts(nrm_ts)
    mn=minimum(nrm_ts.price)
    nrm_ts.price=nrm_ts.price.-mn#set minimum to zero
    mx=maximum(nrm_ts.price)
    nrm_ts.price=nrm_ts.price./mx#normalize
    return nrm_ts
end

function rand_set_days(ukbe_data,dn,wn,mn)
    sampled_set=DataFrame()
    selected_months=[]; for j=1:mn; months=[i for i=1:12 if !issubset(i,selected_months)]; push!(selected_months,rand(months));end;sort!(selected_months)
    for mnth in selected_months;#println("month: "*string(mnth))
        hours_in_month=filter(e->month(e.time_stamp)==mnth,ukbe_data)
        days_in_month=round(Int64,length(hours_in_month[!,:time_stamp])/24)
        weeks_in_month=ceil(Int64,days_in_month/7);selected_weeks=[];
        #for w=1:wn; weeks=[i for i=1:weeks_in_month if !issubset(i,selected_weeks)]; push!(selected_weeks,rand(weeks));end;sort!(selected_weeks);selected_days=[];
        selected_weeks=wn;selected_days=[];
        for wek in selected_weeks;#println(wek)
            for d=1:dn; days=[i for i in (wek-1)*7+1:min((wek)*7,days_in_month) if !issubset(i,selected_days)]; if (length(days)>0); push!(selected_days,rand(days));end;end;sort!(selected_days);end
        for dy in selected_days;
            hours_in_day=filter(e->day(e.time_stamp)==dy,hours_in_month)
            sampled_set=vcat(sampled_set,hours_in_day)
        end
    end
    sort!(sampled_set)
    return sampled_set
end


text:
using FillArrays

"""
    DefaultTrajectoryState

The most basic state representation, with a featured graph, the index of the variable to branch on and optionnaly a list of possible values.
"""
struct DefaultTrajectoryState <: GraphTrajectoryState
    fg::GeometricFlux.FeaturedGraph
    variableIdx::Int
    allValuesIdx::Union{AbstractVector{Int}, Nothing}

    DefaultTrajectoryState(fg, variableIdx, allValuesIdx) = new(fg, variableIdx, allValuesIdx)
    DefaultTrajectoryState(fg, variableIdx) = new(fg, variableIdx, nothing)
end

DefaultTrajectoryState(sr::AbstractStateRepresentation) = throw(ErrorException("missing function DefaultTrajectoryState(::$(typeof(sr)))."))

"""
    BatchedFeaturedGraph

A batched representation of the GeometricFlux.FeaturedGraph, to have a closer implementation to GeometricFlux.
"""
struct BatchedFeaturedGraph{T}
    graph::AbstractArray{T, 3}
    nf::AbstractArray{T, 3}
    ef::AbstractArray{T, 3}
    gf::AbstractMatrix{T}

    BatchedFeaturedGraph{T}(graph, nf, ef, gf) where T = new{T}(graph, nf, ef, gf)
    BatchedFeaturedGraph{T}(
        graph;
        nf=Fill(0, (0, size(graph, 1), size(graph, 3))),
        ef=Fill(0, (0, ne(graph[:,:,1]), size(graph, 3))),
        gf=Fill(0, (0, size(graph, 3)))
    ) where T = new{T}(graph, nf, ef, gf)
end
BatchedFeaturedGraph(graph, nf, ef, gf) = BatchedFeaturedGraph{Float32}(graph, nf, ef, gf)

"""
    BatchedDefaultTrajectoryState

The batched version of the `DefaultTrajectoryState`.

It contains all the information that would be stored in a `FeaturedGraph` but reorganised to enable simultaneous
computation on a few graphs.
"""
Base.@kwdef struct BatchedDefaultTrajectoryState{T} <: NonTabularTrajectoryState
    fg::BatchedFeaturedGraph{T}
    variableIdx::AbstractVector{Int}
    allValuesIdx::Union{AbstractMatrix{Int}, Nothing} = nothing
end
BatchedDefaultTrajectoryState(fg, var, val) = BatchedDefaultTrajectoryState{Float32}(fg=fg, variableIdx=var, allValuesIdx=val)

"""
    Flux.functor(::Type{DefaultTrajectoryState}, s)

Utility function used to load data on the working device.

To be noted: this behavior isn't standard, as the function returned creates a `BatchedDefaultTrajectoryState`
rather than a `DefaultTrajectoryState`. This behavior makes it possible to dynamically split the matrices of insterest
from the `FeaturedGraph` wrapper.
"""
function Flux.functor(::Type{DefaultTrajectoryState}, s)
    adj = Flux.unsqueeze(s.fg.graph, 3)
    nf = Flux.unsqueeze(s.fg.nf, 3)
    ef = Flux.unsqueeze(s.fg.ef, 3)
    gf = Flux.unsqueeze(s.fg.gf, 2)
    allValuesIdx = nothing
    if !isnothing(s.allValuesIdx)
        allValuesIdx = Flux.unsqueeze(s.allValuesIdx, 2)
    end
    return (adj, nf, ef, gf), ls -> BatchedDefaultTrajectoryState{Float32}(
        fg = BatchedFeaturedGraph{Float32}(ls[1]; nf=ls[2], ef=ls[3], gf=ls[4]),
        variableIdx = [s.variableIdx],
        allValuesIdx = allValuesIdx
    )
end

"""
    Flux.functor(::Type{Vector{DefaultTrajectoryState}}, s)

Utility function used to load data on the working device.

To be noted: this behavior isn't standard, as the function returned creates a `BatchedDefaultTrajectoryState`
rather than a `Vector{DefaultTrajectoryState}`. This behavior makes it possible to dynamically creates matrices of
the appropriated size to store all the graphs in 3D tensors.
"""
function Flux.functor(::Type{Vector{DefaultTrajectoryState}}, v)
    maxNode = Base.maximum(s -> size(s.fg.nf, 2), v)
    maxEdge = Base.maximum(s -> size(s.fg.ef, 2), v)
    maxGlobal = Base.maximum(s -> length(s.fg.gf), v)
    batchSize = length(v)

    adj = zeros(eltype(v[1].fg.nf), maxNode, maxNode, batchSize)
    nf = zeros(eltype(v[1].fg.nf), size(v[1].fg.nf, 1), maxNode, batchSize)
    ef = zeros(eltype(v[1].fg.ef), size(v[1].fg.ef, 1), 2*maxEdge, batchSize)
    gf = zeros(eltype(v[1].fg.gf), maxGlobal, batchSize)
    variableIdx = ones(Int, batchSize)

    allValuesIdx = nothing
    if !isnothing(v[1].allValuesIdx)
        maxActions = Base.maximum(s -> length(s.allValuesIdx), v)
        allValuesIdx = zeros(Int, maxActions, batchSize)
    end

    Zygote.ignore() do
        # TODO: this could probably be optimized
        foreach(enumerate(v)) do (idx, state)
            adj[1:size(state.fg.graph,1),1:size(state.fg.graph,2),idx] = state.fg.graph
            nf[1:size(state.fg.nf, 1),1:size(state.fg.nf, 2),idx] = state.fg.nf
            ef[1:size(state.fg.ef, 1),1:size(state.fg.ef, 2),idx] = state.fg.ef
            gf[1:size(state.fg.gf, 1),idx] = state.fg.gf
            variableIdx[idx] = state.variableIdx
            if !isnothing(allValuesIdx)
                allValuesIdx[1:length(state.allValuesIdx), idx] = state.allValuesIdx
            end
        end
    end

    return (adj, nf, ef, gf), ls -> BatchedDefaultTrajectoryState{Float32}(
        fg = BatchedFeaturedGraph{Float32}(ls[1]; nf=ls[2], ef=ls[3], gf=ls[4]),
        variableIdx = variableIdx,
        allValuesIdx = allValuesIdx
    )
end
Flux.@functor BatchedFeaturedGraph
Flux.functor(::Type{BatchedDefaultTrajectoryState{T}}, ts) where T = (ts.fg,), ls -> BatchedDefaultTrajectoryState{T}(ls[1], ts.variableIdx, ts.allValuesIdx)


text:
#####
##### 'Tupled closure' implementation.
#####
##### For now, we have hacked an implementation for tuples of two closures.
##### We cannot do the case of arbitrary closure tuples until we figure out how to
##### obtain the length of a closure tuple at compile-time on the GPU.
#####

# Stress divergences

for stress_div in (:∂ⱼ_2ν_Σ₁ⱼ, :∂ⱼ_2ν_Σ₂ⱼ, :∂ⱼ_2ν_Σ₃ⱼ)
    @eval begin
        @inline $stress_div(i, j, k, grid::AbstractGrid, closures::Tuple{C1, C2}, U, Ks) where {C1, C2} = (
              $stress_div(i, j, k, grid, closures[1], U, Ks[1])
            + $stress_div(i, j, k, grid, closures[2], U, Ks[2]))
    end
end

# Tracer flux divergences

@inline ∇_κ_∇c(i, j, k, grid::AbstractGrid, closures::Tuple{C1, C2}, c, iᶜ, Ks, C, buoyancy) where {C1, C2} = (
      ∇_κ_∇c(i, j, k, grid, closures[1], c, iᶜ, Ks[1], C, buoyancy)
    + ∇_κ_∇c(i, j, k, grid, closures[2], c, iᶜ, Ks[2], C, buoyancy))

# Utilities for closures

function calculate_diffusivities!(Ks, arch, grid, closures::Tuple, args...)
    for (α, closure) in enumerate(closures)
        @inbounds K = Ks[α]
        calculate_diffusivities!(K, arch, grid, closure, args...)
    end
    return nothing
end


text:
export AdditiveMonoidalCategory, oplus, ⊕, mzero,
  AdditiveSymmetricMonoidalCategory, FreeAdditiveSymmetricMonoidalCategory,
  MonoidalCategoryWithCodiagonals, CocartesianCategory, FreeCocartesianCategory,
  mmerge, create, copair, incl1, incl2, ∇, □, braid, σ

import Base: collect, ndims

# Monoidal category
###################

""" Doctrine of *monoidal category* in additive notation

The same as `MonoidalCategory` mathematically but with different notation.
"""
@signature Category(Ob,Hom) => AdditiveMonoidalCategory(Ob,Hom) begin
  oplus(A::Ob, B::Ob)::Ob
  oplus(f::Hom(A,B), g::Hom(C,D))::Hom(oplus(A,C),oplus(B,D)) <=
    (A::Ob, B::Ob, C::Ob, D::Ob)
    @op (⊕) := oplus
  mzero()::Ob
end

# Convenience constructors
oplus(xs::Vector{T}) where T = isempty(xs) ? mzero(T) : foldl(oplus, xs)
oplus(x, y, z, xs...) = oplus([x, y, z, xs...])

# Overload `collect` and `ndims` as for multiplicative monoidal categories.
collect(expr::ObExpr{:oplus}) = vcat(map(collect, args(expr))...)
collect(expr::ObExpr{:mzero}) = roottypeof(expr)[]
ndims(expr::ObExpr{:oplus}) = sum(map(ndims, args(expr)))
ndims(expr::ObExpr{:mzero}) = 0

function show_unicode(io::IO, expr::Union{ObExpr{:oplus},HomExpr{:oplus}}; kw...)
  Syntax.show_unicode_infix(io, expr, "⊕"; kw...)
end
show_unicode(io::IO, expr::ObExpr{:mzero}; kw...) = print(io, "O")

function show_latex(io::IO, expr::Union{ObExpr{:oplus},HomExpr{:oplus}}; kw...)
  Syntax.show_latex_infix(io, expr, "\\oplus"; kw...)
end
show_latex(io::IO, expr::ObExpr{:mzero}; kw...) = print(io, "O")

# Symmetric monoidal category
#############################

""" Doctrine of *symmetric monoidal category* in additive notation

The same as `SymmetricMonoidalCategory` mathematically but with different
notation.
"""
@signature AdditiveMonoidalCategory(Ob,Hom) => AdditiveSymmetricMonoidalCategory(Ob,Hom) begin
  braid(A::Ob, B::Ob)::Hom(oplus(A,B),oplus(B,A))
  @op (σ) := braid
end

@syntax FreeAdditiveSymmetricMonoidalCategory(ObExpr,HomExpr) AdditiveSymmetricMonoidalCategory begin
  oplus(A::Ob, B::Ob) = associate_unit(new(A,B), mzero)
  oplus(f::Hom, g::Hom) = associate(new(f,g))
  compose(f::Hom, g::Hom) = associate(new(f,g; strict=true))
end

# Cocartesian category
######################

""" Doctrine of *monoidal category with codiagonals*

A monoidal category with codiagonals is a symmetric monoidal category equipped
with coherent collections of merging and creating morphisms (monoids).
Unlike in a cocartesian category, the naturality axioms need not be satisfied.

For references, see `MonoidalCategoryWithDiagonals`.
"""
@signature AdditiveSymmetricMonoidalCategory(Ob,Hom) => MonoidalCategoryWithCodiagonals(Ob,Hom) begin
  mmerge(A::Ob)::Hom(oplus(A,A),A)
  @op (∇) := mmerge
  create(A::Ob)::Hom(mzero(),A)
  @op (□) := create
end

""" Doctrine of *cocartesian category*

Actually, this is a cocartesian *symmetric monoidal* category but we omit these
qualifiers for brevity.
"""
@signature MonoidalCategoryWithCodiagonals(Ob,Hom) => CocartesianCategory(Ob,Hom) begin
  copair(f::Hom(A,C), g::Hom(B,C))::Hom(oplus(A,B),C) <= (A::Ob, B::Ob, C::Ob)
  incl1(A::Ob, B::Ob)::Hom(A,oplus(A,B))
  incl2(A::Ob, B::Ob)::Hom(B,oplus(A,B))
end

""" Syntax for a free cocartesian category.

In this syntax, the copairing and inclusion operations are defined using
merging and creation, and do not have their own syntactic elements.
Of course, this convention could be reversed.
"""
@syntax FreeCocartesianCategory(ObExpr,HomExpr) CocartesianCategory begin
  oplus(A::Ob, B::Ob) = associate_unit(new(A,B), mzero)
  oplus(f::Hom, g::Hom) = associate(new(f,g))
  compose(f::Hom, g::Hom) = associate(new(f,g; strict=true))

  copair(f::Hom, g::Hom) = compose(oplus(f,g), mmerge(codom(f)))
  incl1(A::Ob, B::Ob) = oplus(id(A), create(B))
  incl2(A::Ob, B::Ob) = oplus(create(A), id(B))
end

function show_latex(io::IO, expr::HomExpr{:mmerge}; kw...)
  Syntax.show_latex_script(io, expr, "\\nabla")
end
function show_latex(io::IO, expr::HomExpr{:create}; kw...)
  Syntax.show_latex_script(io, expr, "\\square")
end


text:
struct TwoPhaseTargetSelector <: AbstractUpdateSelector
    M_target::Manifold
end

function select_update!(selector::TwoPhaseTargetSelector, state, optimizer, pb)
    state.previous_update = state.selected_update

    if state.M == selector.M_target
        state.selected_update = optimizer.manifold_update
    else
        state.selected_update = optimizer.wholespace_update
    end

    return
end

Base.summary(::TwoPhaseTargetSelector) = "TwoPhaseTarget"