text:
using Stheno: ZeroKernel, OneKernel, ConstKernel, CustomMean, pw, Stretched, Scaled
using Stheno: EQ, Exp, Linear, Noise, PerEQ, Matern32, Matern52, RQ, Sum, Product, stretch,
    Poly, GammaExp, Wiener, WienerVelocity
using LinearAlgebra

@timedtestset "kernel" begin

    @timedtestset "base kernels" begin
        rng, N, N′, D = MersenneTwister(123456), 5, 6, 2
        x0 = collect(range(-2.0, 2.0; length=N)) .+ 1e-3 .* randn(rng, N)
        x1 = collect(range(-1.7, 2.3; length=N)) .+ 1e-3 .* randn(rng, N)
        x2 = collect(range(-1.7, 3.3; length=N′)) .+ 1e-3 .* randn(rng, N′)
        x0_r, x1_r = range(-5.0, step=1, length=N), range(-4.0, step=1, length=N)
        x2_r, x3_r = range(-5.0, step=2, length=N), range(-3.0, step=1, length=N′)
        x4_r = range(-2.0, step=2, length=N′)

        X0 = ColVecs(randn(rng, D, N))
        X1 = ColVecs(randn(rng, D, N))
        X2 = ColVecs(randn(rng, D, N′))

        ȳ, Ȳ, Ȳ_sq = randn(rng, N), randn(rng, N, N′), randn(rng, N, N)

        @timedtestset "ZeroKernel" begin
            differentiable_kernel_tests(ZeroKernel(), ȳ, Ȳ, Ȳ_sq, x0, x1, x2)
            differentiable_kernel_tests(ZeroKernel(), ȳ, Ȳ, Ȳ_sq, X0, X1, X2)
        end

        @timedtestset "OneKernel" begin
            differentiable_kernel_tests(OneKernel(), ȳ, Ȳ, Ȳ_sq, x0, x1, x2)
            differentiable_kernel_tests(OneKernel(), ȳ, Ȳ, Ȳ_sq, X0, X1, X2)
        end

        @timedtestset "ConstKernel" begin
            @test ew(ConstKernel(5), x0) == 5 .* ones(length(x0))
            differentiable_kernel_tests(ConstKernel(5.0), ȳ, Ȳ, Ȳ_sq, x0, x1, x2)
            differentiable_kernel_tests(ConstKernel(5.0), ȳ, Ȳ, Ȳ_sq, X0, X1, X2)
        end

        @timedtestset "EQ" begin
            differentiable_kernel_tests(EQ(), ȳ, Ȳ, Ȳ_sq, x0, x1, x2)
            differentiable_kernel_tests(EQ(), ȳ, Ȳ, Ȳ_sq, X0, X1, X2)
        end

        @timedtestset "PerEQ" begin
            differentiable_kernel_tests(PerEQ(), ȳ, Ȳ, Ȳ_sq, x0, x1, x2; atol=1e-6)
        end

        @timedtestset "Exp" begin
            differentiable_kernel_tests(Exp(), ȳ, Ȳ, Ȳ_sq, x0 .+ 1, x1, x2)
            differentiable_kernel_tests(Exp(), ȳ, Ȳ, Ȳ_sq, X0, X1, X2)
        end

        @timedtestset "Matern32" begin
            differentiable_kernel_tests(Matern32(), ȳ, Ȳ, Ȳ_sq, x0, x1, x2)
            differentiable_kernel_tests(Matern32(), ȳ, Ȳ, Ȳ_sq, X0, X1, X2)
        end

        @timedtestset "Matern52" begin
            differentiable_kernel_tests(Matern52(), ȳ, Ȳ, Ȳ_sq, x0, x1, x2)
            differentiable_kernel_tests(Matern52(), ȳ, Ȳ, Ȳ_sq, X0, X1, X2)
        end

        @timedtestset "RQ" begin
            @timedtestset "α=1.0" begin
                differentiable_kernel_tests(RQ(1.0), ȳ, Ȳ, Ȳ_sq, x0, x1, x2)
                differentiable_kernel_tests(RQ(1.0), ȳ, Ȳ, Ȳ_sq, X0, X1, X2)
            end
            @timedtestset "α=1.5" begin
                differentiable_kernel_tests(RQ(1.5), ȳ, Ȳ, Ȳ_sq, x0, x1, x2)
                differentiable_kernel_tests(RQ(1.5), ȳ, Ȳ, Ȳ_sq, X0, X1, X2)
            end
            @timedtestset "α=100.0" begin
                differentiable_kernel_tests(RQ(100.0), ȳ, Ȳ, Ȳ_sq, x0, x1, x2)
                differentiable_kernel_tests(RQ(100.0), ȳ, Ȳ, Ȳ_sq, X0, X1, X2)
            end
            @timedtestset "single-input" begin
                adjoint_test((α, x, x′)->ew(RQ(α), x, x′), ȳ, 1.5, x0, x1)
                adjoint_test((α, x, x′)->pw(RQ(α), x, x′), Ȳ, 1.5, x0, x2)
                adjoint_test((α, x)->ew(RQ(α), x), ȳ, 1.5, x0)
                adjoint_test((α, x)->pw(RQ(α), x), Ȳ_sq, 1.5, x0)
            end
            @timedtestset "multi-input" begin
                adjoint_test((α, x, x′)->ew(RQ(α), x, x′), ȳ, 1.5, X0, X1)
                adjoint_test((α, x, x′)->pw(RQ(α), x, x′), Ȳ, 1.5, X0, X2)
                adjoint_test((α, x)->ew(RQ(α), x), ȳ, 1.5, X0)
                adjoint_test((α, x)->pw(RQ(α), x), Ȳ_sq, 1.5, X0)
            end
        end

        @timedtestset "Linear" begin
            differentiable_kernel_tests(Linear(), ȳ, Ȳ, Ȳ_sq, x0, x1, x2)
            differentiable_kernel_tests(Linear(), ȳ, Ȳ, Ȳ_sq, X0, X1, X2)
        end

        @timedtestset "Poly" begin
            @test pw(Poly(1, 0.0), x0) ≈ pw(Linear(), x0)
            @testset "Poly{$p}" for p in [1, 2, 3]
                differentiable_kernel_tests(Poly(p, 0.5), ȳ, Ȳ, Ȳ_sq, x0, x1, x2)
                differentiable_kernel_tests(Poly(p, 0.5), ȳ, Ȳ, Ȳ_sq, X0, X1, X2)
            end
        end

        @timedtestset "GammaExp" begin
            @test pw(GammaExp(2.0), x0) ≈ pw(stretch(EQ(), sqrt(2)), x0)
            @testset "γ=$γ" for γ in [1.0, 1.5, 2.0]
                differentiable_kernel_tests(GammaExp(γ), ȳ, Ȳ, Ȳ_sq, x0, x1, x2)
                differentiable_kernel_tests(GammaExp(γ), ȳ, Ȳ, Ȳ_sq, X0, X1, X2)
            end
        end

        @timedtestset "Wiener" begin
            x0 = collect(range(0.1, 2.0; length=N)) .+ 1e-3 .* randn(rng, N)
            x1 = collect(range(0.1, 2.3; length=N)) .+ 1e-3 .* randn(rng, N)
            x2 = collect(range(0.1, 3.3; length=N′)) .+ 1e-3 .* randn(rng, N′)
            kernel_tests(Wiener(), x0, x1, x2)
        end

        @timedtestset "WienerVelocity" begin
            x0 = collect(range(0.1, 2.0; length=N)) .+ 1e-3 .* randn(rng, N)
            x1 = collect(range(0.1, 2.3; length=N)) .+ 1e-3 .* randn(rng, N)
            x2 = collect(range(0.1, 3.3; length=N′)) .+ 1e-3 .* randn(rng, N′)
            kernel_tests(WienerVelocity(), x0, x1, x2)
        end

        @timedtestset "Noise" begin
            @test ew(Noise(), x0, x0) == zeros(length(x0))
            @test pw(Noise(), x0, x0) == zeros(length(x0), length(x0))
            @test ew(Noise(), x0) == ones(length(x0))
            @test pw(Noise(), x0) == Diagonal(ones(length(x0)))
        end

        @timedtestset "Sum" begin
            differentiable_kernel_tests(Sum(EQ(), Exp()), ȳ, Ȳ, Ȳ_sq, x0, x1, x2)
            differentiable_kernel_tests(Sum(EQ(), Exp()), ȳ, Ȳ, Ȳ_sq, X0, X1, X2)
            @test EQ() + Exp() isa Sum
        end

        @timedtestset "Product" begin
            differentiable_kernel_tests(Product(EQ(), Exp()), ȳ, Ȳ, Ȳ_sq, x0, x1, x2)
            differentiable_kernel_tests(Product(EQ(), Exp()), ȳ, Ȳ, Ȳ_sq, X0, X1, X2)
            @test EQ() * Exp() isa Product
        end

        @timedtestset "Scaled" begin
            differentiable_kernel_tests(Scaled(0.5, EQ()), ȳ, Ȳ, Ȳ_sq, x0, x1, x2)
            differentiable_kernel_tests(Scaled(0.5, EQ()), ȳ, Ȳ, Ȳ_sq, X0, X1, X2)
            adjoint_test(σ²->pw(Scaled(σ², EQ()), X0), Ȳ_sq, 0.5)
            @test 0.5 * EQ() isa Scaled
            @test EQ() * 0.5 isa Scaled
        end

        @timedtestset "Stretched" begin
            @timedtestset "Scalar a" begin
                k = stretch(EQ(), 0.1)
                differentiable_kernel_tests(k, ȳ, Ȳ, Ȳ_sq, x0, x1, x2; atol=1e-7, rtol=1e-7)
                differentiable_kernel_tests(k, ȳ, Ȳ, Ȳ_sq, X0, X1, X2; atol=1e-7, rtol=1e-7)
            end
            @timedtestset "Vector a" begin
                k = stretch(EQ(), randn(rng, D))
                differentiable_kernel_tests(k, ȳ, Ȳ, Ȳ_sq, X0, X1, X2)
            end
            @timedtestset "Matrix a" begin
                k = stretch(EQ(), randn(rng, D, D))
                differentiable_kernel_tests(k, ȳ, Ȳ, Ȳ_sq, X0, X1, X2)
            end
            @timedtestset "Convenience Constructors" begin
                unstretched = [
                    (eq(), EQ()),
                    (rq(1.5), RQ(1.5)),
                    (γ_exponential(2.0), GammaExp(2.0)),
                    (poly(2, 0.6), Poly(2, 0.6)),
                ]
                stretched = [
                    (eq(0.7), stretch(EQ(), 0.7)),
                    (rq(1.5, 0.5), stretch(RQ(1.5), 0.5)),
                    (γ_exponential(2.0, 0.53), stretch(GammaExp(2.0), 0.53)),
                    (poly(2, 0.6, 0.4), stretch(Poly(2, 0.6), 0.4)),
                ]
                @testset "$k_conv" for (k_conv, k) in unstretched
                    @test pw(k_conv, x0) ≈ pw(k, x0)
                    differentiable_kernel_tests(k_conv, ȳ, Ȳ, Ȳ_sq, x0, x1, x2)
                end
                @testset "$k_conv" for (k_conv, k) in stretched
                    @test pw(k_conv, x0) ≈ pw(k, x0)
                    differentiable_kernel_tests(k_conv, ȳ, Ȳ, Ȳ_sq, x0, x1, x2)
                end

                # Wiener kernels aren't differentiable.
                not_diff_unstretched = [
                    (wiener(), Wiener()),
                    (wiener_velocity(), WienerVelocity()),
                ]
                not_diff_stretched = [
                    (wiener(0.67), stretch(Wiener(), 0.67)),
                    (wiener_velocity(0.65), stretch(WienerVelocity(), 0.65)),
                ]
                @testset "$k_conv" for (k_conv, k) in not_diff_unstretched
                    @test pw(k_conv, x0) ≈ pw(k, x0)
                    kernel_tests(k_conv, abs.(x0), abs.(x1), abs.(x2))
                end
                @testset "$k_conv" for (k_conv, k) in not_diff_stretched
                    @test pw(k_conv, x0) ≈ pw(k, x0)
                    kernel_tests(k_conv, abs.(x0), abs.(x1), abs.(x2))
                end
            end
        end
    end

    @timedtestset "(is)zero" begin
        @test zero(ZeroKernel()) == ZeroKernel()
        @test zero(EQ()) == ZeroKernel()
        @test iszero(ZeroKernel()) == true
        @test iszero(EQ()) == false
    end
end


text:
# # Two Dimensions
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/random_flow_field.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/random_flow_field.ipynb)
#
# Simulate trajectories of a particle cloud in a two-dimensional flow field.
# A doubly-periodic domain and simple flow fields are used.
# For additional documentation e.g. see :
# [1](https://JuliaClimate.github.io/IndividualDisplacements.jl/dev/),
# [2](https://JuliaClimate.github.io/MeshArrays.jl/dev/),
# [3](https://docs.juliadiffeq.org/latest/solvers/ode_solve.html),
# [4](https://en.wikipedia.org/wiki/Displacement_(vector))
#
# Exercises:
# - change the initial distribution of particles
# - increase the duration of the trajectories simulation
# - treat the non-periodic domain case by padding `u,v` with zeros
# - replace `u,v` with your own two-dimensional flow fields
#
# ![particles in random flow](https://github.com/JuliaClimate/IndividualDisplacements.jl/raw/master/examples/figs/RandomFlow.gif)

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 1. Import Software

using IndividualDisplacements, DataFrames
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/jupyter/flow_fields.jl"));

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 2. Flow Fields
#
# Let's start with a simple, doubly periodic flow field defined by a streamfunction
# and create the `FlowFields` data structure which will then be drive the
# individual displacement and trajectory computations

nx=16
dx= π/nx
XC = dx*(collect(1:2*nx) .- 0.5)
YC = dx*(collect(1:nx) .- 0.5)

fac=0.1
f(x, y) = sin(x) + cos(y) #streamfunction
ϕ = fac*[f(x, y) for x in XC,y in YC] #streamfunction
uC = -fac*[sin(y) for x in XC,y in YC] #dphi/dy at cell center
vC = -fac*[cos(x) for x in XC,y in YC]; #-dphi/dx at cell center

# When starting as done here with collocated velocities in m/s units (`uC,vC` at
# the grid cell center), one easily obtains the normalized, staggered C-grid
# velocities (`u,v`) as shown below.
# `u` (resp `v`) is staggered by `-0.5` grid point in direction `1` for `u` (`2` for `v`)
# from the grid cell center `uC,vC`. This staggering of variables as well as the
# normalization of velocities by the grid scale are key conventions
# adopted in `IndividualDisplacements.jl`.

u=0.5*(circshift(uC, (1,0))+uC) /dx #staggered u converted to grid point units (m/s -> 1/s)
v=0.5*(circshift(vC, (0,1))+vC) /dx #staggered v converted to grid point units (m/s -> 1/s)

𝐹=FlowFields(u,u,v,v,[0.,10.])

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 3. Initialize Individuals
#
# For example, we can initialize 100 particles within a central subdomain as follows.

np,nq=size(u)
x=np*(0.4 .+ 0.2*rand(100))
y=nq*(0.4 .+ 0.2*rand(100));

# The following constructor function wraps everything in the `Individuals` data structure.

𝐼=Individuals(𝐹,x,y)

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 4. Compute Trajectories
#
# The time period is `𝐼.𝑃.𝑇` by default, unless `∫!(𝐼,𝑇)` is called instead.

∫!(𝐼)

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 5. Plot Results

p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/jupyter/recipes_plots.jl"));

🔴_by_t = groupby(𝐼.🔴, :t)
phi_scatter(ϕ,🔴_by_t[end])

# Or to generate a simple animation:
#
# ```
# anim = @animate for t in eachindex(🔴_by_t)
#    phi_scatter(ϕ,🔴_by_t[t])
# end
#
# pth=tempdir()*"/"
# gif(anim, pth*"RandomFlow.gif", fps = 10)
# ```

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 6. Alternatives (optional)

# Making the flow field time variable simply involves providing two different
# flow fields. Linear interpolation will then be applied in the specified time
# range. For example :
#
# ```
# 𝐹=FlowFields(-u,u,-v,v,[0.,10.])
# ```

# It is often convenient to set up `FlowFields` using the MeshArrays.jl package, which
# handles such staggered grids in general fashion, instead of plain Arrays. Given a doubly
# periodic flow field defined as before this is readily done via the `convert_to_FlowFields` function.
#
# ```
# u,v,ϕ=random_flow_field()
# #𝐹=FlowFields(u,u,v,v,[0.,10.])
# 𝐹=convert_to_FlowFields(u,v,10.0);
# ```

# We should note that `u,v` derive from streamfunction `ϕ` which is defined at the corner point in `random_flow_field`.
# Doing this ensures that the resulting `u,v` is non-divergent / purely-rotational over the C-grid domain. In brief:
#
# ```
# u=-(circshift(ϕ, (0,-1))-ϕ)
# v=(circshift(ϕ, (-1,0))-ϕ)
# ```

# In contrast uC, vC computed as earlier in this example may contain contain both
# [rotational and divergent](https://en.wikipedia.org/wiki/Helmholtz_decomposition)
# components -- if this were an important consideration, `MeshArrays.jl` provides
# tools that can be used to e.g. ensure non-divergence.

# When using `MeshArray` flow fields instead of plain `Array` flow fields, initial
# positions also include a subdomain array index (`a`;  all ones in our example).
#
# ```
# a=ones(size(x))
# #isa(𝐹,𝐹_MeshArray2D)
# 𝐼=Individuals(𝐹,x,y,a)
# ```

ext: jl


text:
module Utils
const ⊏ = occursin
end  # module


text:
module MdSlicing
export ProblemSpec, ConvexMesh, fillfps, slice1d, slice2d, sliceplorer, hypersliceplorer, plot

include("types.jl")
include("meshes.jl")
include("hyperslicesets.jl")
include("intersecttest.jl")
include("hyperslice.jl")
include("sliceplorer.jl")
include("plots.jl")

# FIXME: put types back in
# function fillfps(ps::ProblemSpec, fp::Array{Union{Missing,Float64},1}, n::Int = 50)
function fillfps(ps, fp, n=50)
    if !any(map(ismissing, fp)) # special case
        return [fp]
    end
    fpvals =
    [ismissing(x[1]) ? (x[2][1] .+ (x[2][2] - x[2][1]) .* rand(n)) : fill(x[1], n) for x in zip(fp, ps.vals)]
  # transpose everything
    [[fpvals[c][r] for c = 1:length(fp)] for r = 1:n]
end

end # module


text:
using JLD
using DiffPrecTest
using DataFrames, DataFramesMeta, CSV
using Statistics, StatsBase, LinearAlgebra
using SparseArrays
using Random, Distributions
using LaTeXStrings, LaTeXTabulars

cd("exper/sim3")

pArr = [100, 200]
elemArr = [(5,5), (8, 7), (50, 25), (21, 20), (30, 30)]
methodArr = ["Sym-N", "Asym-N", "YinXia", "Sym-B", "Asym-B", "O-Sym-N", "O-Asym-N", "O-Sym-B", "O-Asym-B"]

df = DataFrame(p = Int[], row = Int[], col = Int[], trueValue=Float64[], method=String[], bias=Float64[], coverage=Float64[], lenCI=Float64[])

trueValue = zeros(Float64, 2, 5)

for ip=1:2
  for iElem=1:5
    global trueValue

    p = pArr[ip]
    Random.seed!(1234)

    # generate model
    Ωx = Matrix{Float64}(I, p, p)
    for l=1:p-1
      Ωx[l  , l+1] = 0.6
      Ωx[l+1, l  ] = 0.6
    end
    for l=1:p-2
      Ωx[l  , l+2] = 0.3
      Ωx[l+2, l  ] = 0.3
    end
    Δ = zeros(p, p)
    # generate Delta
    for j=1:5
      i = 5 + (j-1)
      Δ[i, i] = rand(Uniform(0.1, 0.2))
    end
    for j=1:5
      i = 5 + (j-1)
      v = rand(Uniform(0.2, 0.5))
      Δ[i  , i+1] = v
      Δ[i+1, i  ] = v
    end
    Ωy = Ωx - Δ
    d = Vector{Float64}(undef, p)
    rand!(Uniform(0.5, 2.5), d)
    d .= sqrt.(d)
    D = Diagonal(d)
    Σx = inv(Symmetric(D * Ωx * D))
    Σy = inv(Symmetric(D * Ωy * D))

    tΔ = D * Δ * D

    ri, ci = elemArr[iElem]
    indE = (ci - 1) * p + ri

    trueValue[ip, iElem] = tΔ[ri, ci]

    res = load("results/sim3_res_$(ip)_$(iElem).jld", "results")

    for j in 1:9
      sr = computeSimulationResult([res[j, i] for i=1:1000], tΔ[indE])
      push!(df, [pArr[ip], elemArr[iElem][1], elemArr[iElem][2], tΔ[indE], methodArr[j], sr.bias, sr.coverage, sr.lenCoverage])
    end
  end
end

CSV.write("sim3.csv", df)

### create a table

latex_tabular("sim3_table.tex",
  Tabular("lllccc"),
  [
    Rule(:top),
    ["", "", "", "Coverage", "Length", L"{\rm Bias} \times 10^3"],
    Rule(:mid),
    vcat([
      hcat(
        (a = fill("", 15, 1);
        b = latexstring("p = $(pArr[ip])");
        a[1] = "\\multirow{15}{*}{$(b)}";
        a
        ),
        vcat([
          hcat(
            (a = fill("", 3, 1);
            b = latexstring("\\Delta_{$(elemArr[iElem][1]),$(elemArr[iElem][2])} = $(trueValue[ip, iElem])");
            a[1] = "\\multirow{3}{*}{$(b)}";
            a
            ),
            Matrix{Any}(
              @linq df |>
              where( (:method .== "Sym-N") .| (:method .== "YinXia") .| (:method .== "O-Sym-N"),
                    :p .== pArr[ip],
                    :row .== elemArr[iElem][1], :col .== elemArr[iElem][2]) |>  select(:method,
                    coverage = convert.(Int, :coverage * 1000),
                    lenCI = round.(:lenCI, digits=3),
                    bias = round.(:bias * 1000, digits=1))
            )
          )
          for iElem=1:5
          ]...
        )
      )
      for ip=1:2
    ]...),
    Rule(:bottom)
  ]
)


text:
"""
    valid_actions(env)

Return a collection of actions that are appropriate for the current state. This should be a subset of actions(env).
"""
function valid_actions end

"""
    valid_action_mask(env)

Return a mask (any `AbstractArray{Bool}`) indicating which actions are valid.

This function only applies to environments with finite action spaces. If both `valid_actions` and `valid_action_mask` are provided, `valid_actions(env)` should return the same set as `actions(env)[valid_action_mask(env)]`.
"""
function valid_action_mask end

"""
    observations(env)

Return a collection of all observations that might be returned by `observe(env)`.

This function is a *static property* of the environment; the value it returns should not change based on the state.
"""
function observations end


text:
"""
    ReplicatedSample(Z::AbstractVector)

This type represents `K` iid samples ``Z_{i1},\\dotsc, Z_{iK}`` drawn from
the same distribution ``F_i``. In the setting, for which Aurora was developed,
the distribution ``F_i``` is parameterized by its mean ``\\mu_i``, i.e.,

```math
\\mu_i = \\mathbb E_{F_i}[ Z_{ij}],
```

as well as a nuisance parameter ``\alpha_i``,
so that ``F_i = F(\\cdot \\mid \\mu_i, \\alpha_i)`` and

```math
Z_{i1},\\dotsc, Z_{iK} \\mid  \\mid \\mu_i, \\alpha_i \\; \\sim \\; F(\\cdot \\mid \\mu_i, \\alpha_i).
```
"""
struct ReplicatedSample{TS<:AbstractVector,SS,RS} #<: EBayesSample{TS}
    Z::TS
    sorted_Z::SS
    rank_idx::RS
end

function ReplicatedSample(Z)
    sorted_Z = sort(Z)
    rank_idx = StatsBase.ordinalrank(Z)
    ReplicatedSample(Z, sorted_Z, rank_idx)
end

function Base.broadcasted(::typeof(ReplicatedSample), Zs::AbstractMatrix)
    ReplicatedSample.(copy.(eachrow(Zs)))
end

Statistics.mean(Z::ReplicatedSample) = mean(Z.Z)
Statistics.std(Z::ReplicatedSample) = std(Z.Z)

StatsBase.nobs(Z::ReplicatedSample) = length(Z.Z)

function StatsBase.response(Z::ReplicatedSample, j::Integer)
    Z.Z[j]
end

function design_row!(vec, Z::ReplicatedSample, j)
    K = nobs(Z)
    rank_j = Z.rank_idx[j]
    if rank_j > 1
        vec[1:(rank_j-1)] = Z.sorted_Z[1:(rank_j-1)]
    end
    if rank_j < K
        vec[(rank_j):(K-1)] = Z.sorted_Z[(rank_j + 1):K]
    end
    vec
end

function design_response_row!(vec, Z::ReplicatedSample, j)
    vec[1] = 1.0
    K = nobs(Z)
    design_row!(view(vec, 2:K), Z, j)
    vec[K+1] = StatsBase.response(Z, j)
    vec
end

function design_matrix!(matrix, Zs::AbstractVector{<:ReplicatedSample}, j)
    for (i,Z) in enumerate(Zs)
        design_row!(view(matrix, :, i), Z, j)
    end
    matrix
end

function design_matrix(Zs::AbstractVector{<:ReplicatedSample}, j)
    n = length(Zs)
    K = nobs(Zs[1])
    design_matrix!(zeros(K-1,n), Zs, j)
end


text:
# knot.jl
# File for creating knots in the physical space in 1D, 2D or 3D.
@doc """
### calcKnot

Computes a uniform knot vector for the mapping object. It is assumed that the
number of control points is known from before. It is user specified. The knot
vector is computed in the parametric space [0,1] in all the 3 directions

**Arguments**

*  `map` : Object of mapping type

"""->

function calcKnot(map::AbstractMappingType)

  # Assertions
  for i = 1:3
    @assert length(map.edge_knot[i]) >= 2*map.order[i] "Length of knot vector
    not sufficient to ensure nonperiodic/clamped/open knot vector in $i
    direction"
  end

  for di = 1:3
    order = map.order[di]
    nctl = map.nctl[di]
    # Zero out the knot vector
    map.edge_knot[di][:] = 0.0
    step = 1 / (nctl-order+1)
    low = 0.0
    for i = (order+1):nctl
      map.edge_knot[di][i] = low + step
      low += step
    end
    # Fill the end of the knot vector
    map.edge_knot[di][nctl+1:nctl+order] = 1.0
  end

  return nothing
end

@doc """
### calcKnotLMS

Computes the knot vector in the parametric space [0,1]. This is for a least
squares curve approximation.

REFERENCE: Equation 9.68 & eqation 9.69, Page 412 The NURBS book, 2nd Edition

**Inputs**

*  `order` : Order of B-spline
*  `nctl`  : Number of control points for that knot vector
*  `X`     : ordinates used to make the knot vector in parametric space
*  `U`     : Knot vector to be populated

**Outputs**

*  None

"""->

function calcKnotLMS(order, nctl, X, U)
  # REFERENCE : Gaetan's pySpline knots.f90, knots_lms

  @assert length(U) = nctl + order "Length of specified knot vector is not equal
  to the sum of number of control points and order"

  lengthX = length(X)

  # Populate the extremas of the knot vector
  for i = 1:order
    U[i] = 0.0
    U[nctl+i] = 1.0
  end  # End for i = 1:order

  # Populate the middle of the knot vector. Recall the structure of knot vectors
  if isodd(lengthX)
    d = lengthX/(nctl-order+1)
    for j = 1:nctl-order
      i = convert(Int, floor(j*d)) # Conversion for integer based indexing
      alpha = j*d-i
      U[order+j] = (1-alpha)*X[i] + alpha*X[i+2]
    end
  else
    d = lengthX/(nctl-order+1)
    for j = 1:nctl-order
      i = convert(Int, floor(j*d)) # Conversion for integer based indexing
      alpha = j*d-i+0.5
      U[order+j] = (1-alpha)*X[i] + alpha*X[i]
    end
  end  # End if isodd(lengthX)

  return nothing
end

@doc """
### calcKnotInterp

Create a knot vector for an interpolating spline. This is again done in the
parametric space between [0,1]

**Inputs**

*  `X`  : Ordinates used to make the knot vector in the parametric space
*  `nd` : Number of derivatives specified
*  `order` : Order of the B-spline curve
*  `U`  : Knot vector to be computed

**Outputs**

*  None

"""->
function calcKnotInterp(X, nd, order, U)

  lengthX = length(X)
  nctl = lengthX + nd

  @assert nctl+order == length(U) # Sanity check

  # Populate the ends of the knot vector in the parametric space
  for j = 1:order
    U[j] = 0.0
    U[nctl+j] = 1.0
  end

  if nd == lengthX
    if order == 3
      for i = 1:(lengthX-2)
        U[order+2*i-1] = 0.5*(X[i] + X[i+1])
        U[order+2*i] = X[i+1]
      end # end for i = 1:(lengthX-2)
      U[nctl] = 0.5*(X[lengthX-1] + 1)
    elseif k == 4
      U[5] = 0.5 * X[2]
      U[nctl] = 0.5 * (X[lengthX-1] + 1)
      for i = 1:(lengthX-3)
        U[order+2*i] = (1/3)*(2*X[i+1] + X[i+2])
        U[order+2*i+1] = (1/3)*(X[i+1] + 2*X[i+2])
      end  # End for i = 1:(lengthX-3)
    end  # end if order == 3
  elseif nd == 0
    if order == 2
      for j = 1:lengthX-order
        U[order+j] = X[j+1]
      end
    elseif order == 3
      for j = 1:lengthX-order
        U[order+j] = 0.5*(X[j+1] + X[j+2])
      end
    elseif order == 4
      for j = 1:lengthX-order
        U[order+j] = X[j+2]
      end
    else
      error("Interpolation is only available for k = 2, 3 or 4")
    end
  else
    error("More options for nd not available yet")
  end  # End if nd == lengthX

  return nothing
end


text:
# Consts

# Variables

# methods
function foo()
    tag::String = "Hello there!";
    println(tag);

    i::Int = 1;
    j::Int = 2;
    sum::Float16 = i+j;
    println(sum);
end;

foo();

max_stars_repo_name: fahimfarhan/legendary-coding-odyssey

max_issues_repo_name: fahimfarhan/legendary-coding-odyssey

max_forks_repo_name: fahimfarhan/legendary-coding-odyssey


text:
using GPUArrays, Test

include("testsuite.jl")

@testset "JLArray" begin
    include("jlarray.jl")
    using .JLArrays

    jl([1])

    JLArrays.allowscalar(false)
    TestSuite.test(JLArray)
end

@testset "Array" begin
    TestSuite.test(Array)
end


text:
module Benchmark
using GerryChain
using LightGraphs
using PyPlot
using LinearAlgebra
using StatsBase
using ArgParse

export main

include("test.jl")
include("main.jl")

end # module

max_stars_repo_name: mggg/benchmarks

max_issues_repo_name: mggg/benchmarks

max_forks_repo_name: mggg/benchmarks


text:
using Compat
using FilePaths
import Playground

# Get our current working path
deps_dir = parent(Path(@__FILE__))

# Setup the build directory
build_dir = join(deps_dir, "usr", "build")
mkdir(build_dir; recursive=true, exist_ok=true)

bin_exec = haskey(ENV, "PLAYGROUND_BIN_EXEC") ? parse(ENV["PLAYGROUND_BIN_EXEC"]) : false

julia_bin = joinpath(Base.JULIA_HOME, "julia")
script =
    """
    #!$julia_bin --optimize=0

    import Playground
    Playground.main(ARGS)
    """

if bin_exec
    # Actually build the playground executable
    info("Trying to build playground executable in $build_dir ...")
    using BuildExecutable
    # We avoid passing a Path to build_executable as it can cause problems.
    build_executable("playground", "script.jl", String(build_dir), "generic"; force=true)
else
    info("Writing a playground script to $build_dir:\n```\n$script```\n")
    write(join(build_dir, "playground"), script)
    chmod(join(build_dir, "playground"), "+x")
end

config_file = join(build_dir, "config.yml")

info("Writing default config to $config_file.")
if exists(config_file)
    backup_file = join(build_dir, ".config.yml_$(Dates.today()).bak")
    info("Backing up existing config file to $backup_file")
    copy(config_file, backup_file; exist_ok=true, overwrite=true)
end

write(config_file, Playground.DEFAULT_CONFIG)

if is_unix()
    copy(
        join(deps_dir, "usr", "bin", "INSTALL.sh"),
        join(build_dir, "INSTALL.sh");
        exist_ok=true,
        overwrite=true
    )
end

copy(join(deps_dir, "..", "LICENSE"), join(build_dir, "LICENSE"); exist_ok=true, overwrite=true)
copy(join(deps_dir, "..", "README.md"), join(build_dir, "README.md"); exist_ok=true, overwrite=true)

# Only install the config and executable if the
# PLAYGROUND_INSTALL env variable has been set.
# This is just cause there isn't a `Pkg.install` or
# `Pkg.build("Pkg", install=true)`
install = haskey(ENV, "PLAYGROUND_INSTALL") ? parse(ENV["PLAYGROUND_INSTALL"]) : false

# Store our install paths
install_dir = join(home(), ".playground")
config_installed = join(install_dir, "config.yml")
require_installed = join(install_dir, "REQUIRE")
playground_installed = join(install_dir, "bin", p"playground")
playground_compiled = join(build_dir, "playground")

if install
    # Set up the user level playground directory
    info("Setting up user playground directory...")
    mkdir(install_dir; recursive=true, exist_ok=true)
    mkdir(join(install_dir, p"bin"); recursive=true, exist_ok=true)

    info("Linking playground config to $config_installed.")

    if exists(config_installed)
        info("~/.playground/config.yml already exists. Skipping.")
        info("Please see $config_file if you have any problems with your existing config.yml file.")
    else
        symlink(config_file, config_installed; exist_ok=true, overwrite=true)
    end

    info("Linking playground executable to $playground_installed")

    if exists(playground_installed)
        backup_file = join(install_dir, "bin", ".playground_$(Dates.today()).bak")
        info("Backing up existing playground executable to $backup_file")
        copy(playground_installed, backup_file; exist_ok=true, overwrite=true)
    end

    symlink(playground_compiled, playground_installed; exist_ok=true, overwrite=true)

    if !exists(require_installed)
        info("Creating an empty REQUIRE file at $require_installed.")
        touch(require_installed)
    end

    info(
        "Adding $(join(install_dir, "bin")) to your PATH " *
        "variable will make `playground` and any julia versions installed via " *
        "Playground.jl available on your search path."
    )
else
    warn(
        "Compiled playground executable $playground_compiled " *
        "not installed to $playground_installed"
    )
end


text:
using Test
using HTensors
using BenchmarkTools

include("make.jl")

u1 = rand(10,4)
u2 = rand(9,6)
u3 = rand(8,3)
u4 = rand(11,4)
u12 = rand(4,6,7)
u34 = rand(3,4,5)
u1234 = rand(7,5,1)

t1234 = make4(u1,u2,u3,u4,u12,u34,u1234)
@test nmode(t1234) == 4
@test shape(t1234) == [10,9,8,11]

a1234 = array(t1234)
@test size(a1234) == (10,9,8,11,1)

te1234 = HTensorEvaluator(t1234)
te1234c = HTensorEvaluator(t1234,1e6)

idxs = [2, 3, 5, 8]
@test a1234[idxs..., :] == getelts(t1234, idxs)
@test a1234[idxs..., :] ≈ getelts(te1234, idxs)
@test a1234[idxs..., :] ≈ getelts(te1234c, idxs)

bm = @benchmark getelts($te1234, $idxs)
@test memory(bm) == 0
bm = @benchmark getelts($te1234c, $idxs)
@test memory(bm) == 0

idxs = [3:8, 2, :, [3,5,7]]
@test a1234[idxs..., :] == getelts(t1234, idxs)


text:
if !isdefined(Main, :CustomStencil)
    include("../module/CustomStencil.jl")
    using Main.CustomStencil
end
include("misc.jl")
radius = 4
l = 2*radius + 1
coefs = zeros(Float64, (l,l,l))

for i = 1:radius
    global radius
    posi = radius + i + 1
    negi = radius -i + 1
    mid = radius + 1
    coefs[posi,mid,mid] = 1
    coefs[negi,mid,mid] = 1
    coefs[mid, posi, mid] = 1
    coefs[mid, negi, mid] = 1
    coefs[mid, mid, posi] = 1
    coefs[mid, mid, negi] = 1

end
mid = radius + 1
coefs[mid,mid,mid] = 1
#
# star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1, 4, c[i]*(
#             D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
#             D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))

st_inst = NewStencilInstance(coefs)

nx = 8
ny = 8
nz = 8
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

t_steps = 5
## Actual time steps
data = CreateData(dx,dy,dz,radius)
dout = ApplyStencil(st_inst, data, t_steps)
#heatmap(dout1[:,:,Int(round(dz/2))])
o = Array(dout)


text:
using CategoricalArrays, DataFrames, GLM, CSV, Statistics

function recode_raw_survey_data!(df)
    recode!(df.Q3, 1=> 10.0, 2=> 12.0, 3=> 13.0, 4=> 14.0, 5=> 16.0, 6=> 18.0, 7=> 20.0)
    recode!(df.perturb1, 2=> 1.0, 1=> -1.0)
end

function confusion(var1, var2)
    println("11: ", sum(var1.*var2)/length(var1))
    println("10: ", sum(var1.*(1 .-var2)/length(var1)))
    println("01: ", sum((1 .-var1).*var2/length(var1)))
    println("00: ", sum((1 .-var1).*(1 .-var2)/length(var1)))
end

function clean_data(datapath)
    data = DataFrame!(CSV.File(datapath))
    recode_raw_survey_data!(data)
    recode!(data.Q11, missing => -100)
    data.correct = (data.Q11 .≈ 3) .+ (data.Q11 .≈ -0.5)
    println(nrow(data))
    return head(data, 316)
end

function summary_stats(datapath)
    data = clean_data(datapath)
    println("Conditional prob I like math: ", mean((data.Q10 .>3).*data.correct)/mean(data.correct))
    println("Cond prob incorrect: ", mean((data.Q10 .>3).*(1 .- data.correct))/mean(1 .- data.correct))
    for var in [data.Q10, data.Q12, data.Q3, data.correct]
        println(mean(var), " (", std(var), ") ")
    end
    print(lm(@formula(Q10 ~ correct ), data))
    #println("Education and Skill: ", cor(data.Q3, data.Q12))
    #println("Correct and Outcome: ", cor(data.correct, data.Q10))
end

function linReg(y, x)
    x = [ones(length(y)) x ]
    β = inv(x'*x)*x'*y
    return β
end

function calculate_gradients(datapath, β, psize)
    data = clean_data(datapath)
    yhat = β[1] .+ (data.perturb1.*psize .+ β[2]).*data.correct
    π = -1 .*(yhat .- data.Q10).^2
    gradient = linReg(π, data.perturb1 .* psize)[2]
    println("Gradient: ", gradient)
    return gradient
end

function get_b0_gradient(datapath, β, psize)
    data = clean_data(datapath)
    net = mean(data.Q10) - β[1] - mean(data.correct.*β[2])
    println(net)
    return net
end

function get_intercept(datapath, β)
    data = clean_data(datapath)
    β0  = mean(data.Q10) - mean(β.*data.correct)
    return β0
end

function o_mse_calc(files, beta)
    mses = zeros(length(files))
    n  = zeros(length(files))
    for (i, f) in zip(1:length(files), files)
        data = clean_data(f)
        x = Matrix([ones(size(data)[1]) data[!, ["correct"]]])
        y = data.Q10
        yhat = x*beta
        mses[i] =  1 - mean((yhat .- y).^2)/var(y)
        println("MSE: ", mean((yhat .- y).^2))
        n[i] = length(yhat)
    end
    println(mses)
    return sum(mses.*n)/sum(n)
end

max_stars_repo_name: evanmunro/personalized-policy

max_issues_repo_name: evanmunro/personalized-policy

max_forks_repo_name: evanmunro/personalized-policy


text:
# Uncomment the depwarns when we drop 0.4 support
# See also ngenerate.jl.

macro Dict(pairs...)
    esc(Expr(:block, # :(Base.depwarn("@Dict is deprecated, use Dict instead", Symbol("@Dict"))),
                       Expr(:call, :Dict, pairs...)))
end
macro AnyDict(pairs...)
    esc(Expr(:block, # :(Base.depwarn("@AnyDict is deprecated, use Dict{Any,Any} instead",
                     #               Symbol("@AnyDict"))),
             Expr(:call, :(Base.AnyDict), pairs...)))
end

if VERSION >= v"0.4.0-dev+3184"
    include("ngenerate.jl")
    using .CompatCartesian
    export @ngenerate, @nsplat
end


text:
struct Multivector{V<:MVector}<:GeometricAlgebraType
    coefs::V
end

Multivector(coefs::AbstractVector) = Multivector(MVector{length(coefs)}(coefs))

vectype(::Type{<:Multivector{V}}) where {V} = V

eltype(T::Type{<:Multivector}) = eltype(vectype(T))
eltype(::Type{Multivector}) = Any
eltype(mv::Multivector) = eltype(typeof(mv))

length(T::Type{<:Multivector}) = length(vectype(T))
length(mv::Multivector) = length(typeof(mv))

zero(::Type{<:Multivector{V}}) where {V} = Multivector(zeros(V))

grades(mv::Multivector) = unique(grade.(nonzero_blades(mv)))

setindex!(mv::Multivector, args...) = setindex!(mv.coefs, args...)

store!(mv::Multivector, b::Blade) = setindex!(mv, b.coef, b.index)

add!(mv::Multivector, b::Blade) = setindex!(mv, mv[b.index] + b.coef, b.index)

getindex(mv::Multivector, args...) = getindex(mv.coefs, args...)

index(::Multivector, i) = i
index(::Multivector, b::Blade) = b.index

broadcastable(mv::Multivector) = mv.coefs

(==)(x::Multivector, y::Multivector) = all(x.coefs .== y.coefs)
@commutative (==)(x::Blade, y::GeometricAlgebraType) = x.coef == value(y, x)

(≈)(x::Multivector, y::Multivector; kwargs...) = all(.≈(x.coefs, y.coefs; kwargs...))
(≈)(x::Blade, y::GeometricAlgebraType; kwargs...) = ≈(x.coef, value(y, x); kwargs...)
(≈)(y::GeometricAlgebraType, x::Blade; kwargs...) = ≈(x.coef, value(y, x); kwargs...)

function (==)(x::GeometricAlgebraType, y::GeometricAlgebraType)
    bx, by = collect(nonzero_blades(x)), collect(nonzero_blades(y))
    length(bx) == length(by) && all(bx .== by)
end

function (≈)(x::GeometricAlgebraType, y::GeometricAlgebraType; kwargs...)
    bx, by = collect(nonzero_blades(x)), collect(nonzero_blades(y))
    length(bx) == length(by) && all(.≈(bx, by; kwargs...))
end

indices(mv::Multivector) = 1:length(mv)

is_homogeneous(mv::Multivector) = length(grades(mv)) ≤ 1

function convert(T::Type{<:Multivector}, kvec::KVector)
    mv = zero(Multivector, eltype(T) == Any ? eltype(kvec) : promote_type(eltype(T), eltype(kvec)))
    mv.coefs[indices(kvec)] .= kvec.coefs
    mv
end

function kvectors end

function kvector end

function show(io::IO, @nospecialize mv::Multivector)
    bs = nonzero_blades(mv)
    if isempty(bs)
        print(io, string(scalar(zero(eltype(mv)))))
    else
        print(io, join(string.(collect(bs)), " + "))
    end
end


text:
export WLS, RandomCascade
using Wavelets, Statistics

"""
    WLS(shufflemethod::Surrogate = IAAFT();
        f::Union{Nothing, Function} = Statistics.cor,
        rescale::Bool = true,
        wt::Wavelets.WT.OrthoWaveletClass = Wavelets.WT.Daubechies{16}())

A wavelet surrogate generated by the following procedure:

1. Compute the wavelet transform of the signal. This results in a set of
    detail coefficients over a set of dyadic scales. As in Keylock (2006),
    we here use the maximal overlap discrete wavelet transform, or MODWT,
    so that the number of coefficients at each scale are the same.
2. Shuffle the detail coefficients at each dyadic scale using the
    provided `shufflemethod`. See "Shuffling methods" below for alternatives.
3. Apply the inverse wavelet transform to the shuffled detail coefficients
    to obtain a surrogate time series.

## Shuffling methods

You may choose to use any surrogate from this package to perform the
randomization of the detail coefficients at each dyadic scale.

The following methods have been discussed in the literature (more may exist):

- Random permutations of wavelet coefficients within each scale (Breakspear et al., 2003). To get this behaviour, use `WLS(x, RandomShuffle(), rescale = false, f = nothing)`.
- Cyclic rotation of wavelet coefficients within each scale (Breakspear et al., 2003).  To get this behaviour, use `WLS(x, Circshift(1:length(x)), rescale = false, f = nothing)`.
- Block resampling of wavelet coefficients within each scale (Breakspear et al., 2003). To get this behaviour, use `WLS(x, BlockShuffle(nblocks, randomize = true), rescale = false, f = nothing)`.
- IAAFT resampling of wavelet coefficients within each scale (Keylock, 2006). To get this behaviour, use `WLS(x, IAAFT(), rescale = true, f = Statistics.cor)`.
    This method preserves the local mean and variance structure of the signal, but
    randomises nonlinear properties of the signal (i.e. Hurst exponents)[^Keylock2006].
    These surrogates can therefore be used to test for changes in nonlinear properties
    of the original signal. In contrast to IAAFT surrogates, the IAAFT-wavelet surrogates
    also preserves nonstationarity. Using other `shufflemethod`s does not necessarily
    preserve nonstationarity. To deal with nonstationary signals, Keylock (2006) recommends
    using a wavelet with a high number of vanishing moments. Thus, our default is to
    use a Daubechies wavelet with 16 vanishing moments. *Note: The iterative procedure after
    the rank ordering step (step [v] in [^Keylock2006]) is not performed in
    this implementation.*

The default method and parameters replicate the behaviour of Keylock (2006)'s IAAFT
wavelet surrogates.

## Error minimization

For the [`IAAFT`](@ref) approach introduced in Keylock (2006), detail coefficients
at each level are circularly rotated to minimize an error function. The methods
introduced in Breakspear et al. (2003) do not apply this error minimization.

In our implementation, you can turn this option on/off using the `f` parameter of
the `WLS` constructor. If `f = nothing` turns off error minization. If `f` is set
to a two-argument function that computes some statistic, for example
`f = Statistics.cor`, then detail coefficients at each scale are circularly
rotated until that function is maximized (and hence the "error" minimized).
If you want to *minimize* some error function, then instead provide an appropriate
transform of your function. For example, if using the root mean squared deviation,
define `rmsd_inv(x, y) = 1 - StatsBase.rmsd(x, y)` and set `f = rmsd_inv`.

## Rescaling

If `rescale == true`, then surrogate values are mapped onto the
values of the original time series, as in the [`AAFT`](@ref) algorithm.
If `rescale == false`, surrogate values are not constrained to the
original time series values. If [`AAFT`](@ref) or [`IAAFT`](@ref) shuffling
is used, `rescale` should be set to `true`. For other methods, it does not
necessarily need to be.

[^Breakspear2003]: Breakspear, M., Brammer, M., & Robinson, P. A. (2003). Construction of multivariate surrogate sets from nonlinear data using the wavelet transform. Physica D: Nonlinear Phenomena, 182(1-2), 1-22.
[^Keylock2006]: C.J. Keylock (2006). "Constrained surrogate time series with preservation of the mean and variance structure". Phys. Rev. E. 73: 036707. doi:10.1103/PhysRevE.73.036707.
"""
struct WLS{WT <: Wavelets.WT.OrthoWaveletClass, S <: Surrogate, E <: Union{Nothing, Function}} <: Surrogate
    shufflemethod::Surrogate # should preserve values of the original series
    rescale::Bool
    f::E
    wt::WT

    function WLS(method::S = IAAFT(); rescale::Bool = true, wt::WT = Wavelets.WT.Daubechies{16}(), f::E = Statistics.cor) where {S <: Surrogate, WT <: Wavelets.WT.OrthoWaveletClass, E <: Union{Nothing, Function}}
        new{WT, S, E}(method, rescale, f, wt)
    end
end

# Initialize without error minimization.
function get_init_noerrorminimize(x, method, rng = Random.default_rng())
    wl = wavelet(method.wt)
    L = length(x)
    x_sorted = sort(x)

    # Wavelet coefficients (step [i] in Keylock)
    W = modwt(x, wl)
    T = eltype(W)
    R = zeros(size(W))

    Nscales = ndyadicscales(L)

    # Will contain surrogate realizations of the wavelet coefficients
    # at each scale (step [ii] in Keylock).
    sW = zeros(T, size(W))

    # Surrogate generators for each set of coefficients
    sgs = [surrogenerator(W[:, i], method.shufflemethod, rng) for i = 1:Nscales]

    init = (
        wl = wl,
        W = W,
        Nscales = Nscales,
        sW = sW,
        sgs = sgs,
        x_sorted = x_sorted,
        R = R,
    )

    return init
end

# Initialize error minimization, as in Keylock (2006).
function get_init_errorminimize(x, method::WLS, rng = Random.default_rng())
    wl = wavelet(method.wt)
    L = length(x)
    x_sorted = sort(x)

    # Wavelet coefficients (step [i] in Keylock)
    W = modwt(x, wl)
    T = eltype(W)
    Nscales = ndyadicscales(L)

    # Will contain surrogate realizations of the wavelet coefficients
    # at each scale (step [ii] in Keylock).
    sW = zeros(T, size(W))

    # We will also need a matrix to store the mirror images of the
    # surrogates (last part of step [ii])
    sWmirr = zeros(T, size(W))

    # Surrogate generators for each set of coefficients
    sgs = [surrogenerator(W[:, i], method.shufflemethod, rng) for i = 1:Nscales]

    # Temporary array for the circular shift error minimizing step
    circshifted_s = zeros(T, size(W))
    circshifted_smirr = zeros(T, size(W))
    R = zeros(size(W))

    s = similar(x)

    init = (wl = wl, W = W, Nscales = Nscales, L = L,
            sW = sW, sgs = sgs, sWmirr = sWmirr,
            circshifted_s = circshifted_s,
            circshifted_smirr = circshifted_smirr,
            x_sorted = x_sorted, R = R,
    )

    return init
end

function surrogenerator(x::AbstractVector{T}, method::WLS, rng = Random.default_rng()) where T
    init = isnothing(method.f) ?
        get_init_noerrorminimize(x, method, rng) :
        get_init_errorminimize(x, method, rng)

    s = similar(x)

    return SurrogateGenerator(method, x, s, init, rng)
end

function (sg::SurrogateGenerator{<:WLS})()
    if isnothing(sg.method.f)
        wls_noerrorminimize(sg)
    else
        wls_errorminimize(sg)
    end
end

function wls_noerrorminimize(sg::SurrogateGenerator{<:WLS})
    s = sg.s

    fds = (:wl, :W, :Nscales, :sW, :sgs, :x_sorted, :R)
    wl, W, Nscales, sW, sgs, x_sorted, R = getfield.(Ref(sg.init), fds)

    for λ in 1:Nscales
        sW[:, λ] .= sgs[λ]()
    end

    s .= imodwt(sW, wl)

    if sg.method.rescale
        s[sortperm(s)] .= x_sorted
    end

    return s
end

function wls_errorminimize(sg::SurrogateGenerator{<:WLS})
    s = sg.s

    # Error minimization function (if we reached this function,
    # f is never `Nothing`), so we can use it safely.
    f = sg.method.f

    fds = (:wl, :W, :Nscales, :L, :sW, :sgs, :sWmirr,
        :circshifted_s, :circshifted_smirr,
        :x_sorted, :R)

    wl, W, Nscales, L, sW, sgs, sWmirr,
        circshifted_s, circshifted_smirr,
        x_sorted, R = getfield.(Ref(sg.init), fds)

    # Create surrogate versions of detail coefficients at each dyadic scale [first part of step (ii) in Keylock]
    for λ in 1:Nscales
        sW[:, λ] .= sgs[λ]()
    end

    # Mirror the surrogate coefficients [last part of step (ii) in Keylock]
    sWmirr .= reverse(sW, dims = 1)

    # In the original paper, surrogates and mirror images are matched to original
    # detail coefficients in a circular manner until some error criterion is
    # minimized. Then, the surrogate or its mirror image, depending on which provides
    # the best fit to the original coefficients, is chosen as the representative
    # for a particular dyadic scale. Here, we instead use maximal correlation as
    # the criterion for matching.
    optimal_shifts = zeros(Int, Nscales)
    optimal_shifts_mirr = zeros(Int, Nscales)
    maxcorrs = zeros(Nscales)
    maxcorrs_mirr = zeros(Nscales)

    for i in 0:L-1
        circshift!(circshifted_s, sW, (i, 0))
        circshift!(circshifted_smirr, sWmirr, (i, 0))

        for λ in 1:Nscales
            origW = W[:, λ]
            c = f(origW, circshifted_s[:, λ])
            if c > maxcorrs[λ]
                maxcorrs[λ] = c
                optimal_shifts[λ] = i
            end

            c_mirr = f(origW, circshifted_smirr[:, λ])
            if c_mirr > maxcorrs_mirr[λ]
                maxcorrs_mirr[λ] = c_mirr
                optimal_shifts_mirr[λ] = i
            end
        end
    end

    # Decide which coefficients are retained (either surrogate or mirror surrogate coefficients)
    for λ in 1:Nscales
        if maxcorrs[λ] >= maxcorrs_mirr[λ]
            R[:, λ] .= circshift(sW[:, λ], optimal_shifts[λ])
        else
            R[:, λ] .= circshift(sWmirr[:, λ], optimal_shifts_mirr[λ])
        end
    end

    s .= imodwt(R, wl)

    if sg.method.rescale
        s[sortperm(s)] .= x_sorted
    end

    return s
end

"""
    RandomCascade(paddingmode::String = "zeros")

A random cascade multifractal wavelet surrogate (Paluš, 2008)[^Paluš2008].

If the input signal length is not a power of 2, the signal must be
padded before the surrogate is constructed. `paddingmode` determines
how the signal is padded. Currently supported padding modes: `"zeros"`.
The final surrogate (constructed from the padded signal) is subset
to match the length of the original signal.

Random cascade surrogate preserve multifractal properties of the input
time series, that is, interactions among dyadic scales and nonlinear
dependencies[^Paluš2008].

[^Paluš2008]: Paluš, Milan (2008). Bootstrapping Multifractals: Surrogate Data from Random Cascades on Wavelet Dyadic Trees. Physical Review Letters, 101(13), 134101–. doi:10.1103/PhysRevLett.101.134101
"""
struct RandomCascade{WT <: Wavelets.WT.OrthoWaveletClass} <: Surrogate
    wt::WT
    paddingmode::String

    function RandomCascade(; wt::WT = Wavelets.WT.Daubechies{16}(), paddingmode::String = "zeros") where {WT <: Wavelets.WT.OrthoWaveletClass}
        new{WT}(wt, paddingmode)
    end
end

function surrogenerator(x::AbstractVector{T}, method::RandomCascade, rng = Random.default_rng()) where T
    nlevels = ndyadicscales(length(x))
    mode = method.paddingmode

    # Pad input so that input to discrete wavelet transform has length which is a power of 2
    x̃ = zeros(2^(nlevels + 1))
    if mode == "zeros"
        copyto!(x̃, x)
    else
        throw(ArgumentError("""`paddingmode` must be one of ["zeros"]"""))
    end

    wl = wavelet(method.wt)

    # Wavelet coefficients (step [i] in Keylock)
    c = dwt(x̃, wl, nlevels)

    # Surrogate coefficients will be partly identical to original coefficients,
    # so we simply copy them and replace the necessary coefficients later.
    cₛ = copy(c)

    # Multiplication factors and index vectors can be pre-allocated for
    # levels 2:nlevels-1; they are overwritten for each new surrogate.
    Ms = [zeros(dyadicdetailn(j-1)) for j = 2:nlevels-1]
    ixs = [zeros(Int, dyadicdetailn(j-1)) for j = 2:nlevels-1]

    init = (
        wl = wl,
        c = c,
        cₛ = cₛ,
        nlevels = nlevels,
        s̃ = similar(x̃),
        Ms = Ms,
        ixs = ixs,
    )

    return SurrogateGenerator(method, x, similar(x), init, rng)
end

function (sg::SurrogateGenerator{<:RandomCascade})()
    s, rng = sg.s, sg.rng
    c, cₛ, s̃, wl, nlevels, Ms, ixs =
        sg.init.c, sg.init.cₛ, sg.init.s̃, sg.init.wl, sg.init.nlevels, sg.init.Ms, sg.init.ixs

    cₛ[dyadicdetailrange(0)] = @view c[dyadicdetailrange(0)]
    cₛ[dyadicdetailrange(1)] = @view c[dyadicdetailrange(1)]

    for (l, j) = enumerate(2:nlevels-1)
        cⱼ₋₁ = @view c[dyadicdetailrange(j - 1)]
        cⱼ = @view c[dyadicdetailrange(j)]

        M = Ms[l]
        ct = 1
        @inbounds for k = 1:length(cⱼ₋₁)
            if k % 2 == 0
                M[ct] = cⱼ[2*k] / cⱼ₋₁[k]
            else
                M[ct] = cⱼ[2*(k+1)] / cⱼ₋₁[k+1]
            end
            ct += 1
        end

        shuffle!(rng, M)
        new_coeffs!(M, cⱼ₋₁)
        ix = ixs[l]
        sortperm!(ix, M)
        cₛ[dyadicdetailrange(j-1)] .= @view cⱼ₋₁[ix]
    end
    s̃ .= idwt(cₛ, wl, nlevels)

    # Surrogate length must match length of original signal.
    s .= @view s̃[1:length(s)]
    return s
end

function new_coeffs!(M, cⱼ₋₁)
    @inbounds for k = 1:length(cⱼ₋₁)
        M[k] = M[k] * cⱼ₋₁[k]
    end
end


text:
module Getters
using PhaseSpaceIO: IAEAParticle, EGSParticle

FIELDNAMES = union(Set(fieldnames(IAEAParticle{0,0})),
               # Set(fieldnames(EGSParticle{Float32}))
               Set([:latch, :zlast])
              )

for field in FIELDNAMES
    @eval ($field)(p) = p.$field
    @eval export $field
end

export energy
energy(p) = E(p)*weight(p)

end


text:
using RoundingEmulator
using Test

include("special_values.jl")
include("setrounding_raw.jl")


text:
using Base.Test
using PDFExtract

contents = readpdf("sample.pdf")

images = saveimages("sample.pdf", options=["-dpi","100"])