text:
using GAFramework

# Solves for function a + 2b + 3c + 4d + 5e = 42
# by minimizing abs(a + 2b + 3c + 4d + 5e - 42)
m = CoordinateModel(x -> abs(x[1] + 2x[2] + 3x[3] + 4x[4] + 5x[5] - 42),
                  [-100.,-100.,-100.,-100.,-100], # minimum corner
                  [100.,100,100.,100.,100]) # maximum corner in rectangle
# Prints fitness every iteration
state = GAState(m, ngen=100, npop=3_000, elite_fraction=0.01,
                mutation_params=Dict(:rate=>0.1), print_fitness_iter=10)
println(ga(state))

# Solves for function a + 2b + 3c + 4d + 5e = 42
# by minimizing abs(a + 2b + 3c + 4d + 5e - 42)
# without restricting the result to the [-100,100] rectangle
# but still initializing the population in that rectangle
m = CoordinateModel(x -> abs(x[1] + 2x[2] + 3x[3] + 4x[4] + 5x[5] - 42),
                  [-100.,-100.,-100.,-100.,-100], # minimum corner
                  [100.,100,100.,100.,100], false) # maximum corner in rectangle
# Prints fitness every iteration
state = GAState(m, ngen=100, npop=3_000, elite_fraction=0.01,
                mutation_params=Dict(:rate=>0.1), print_fitness_iter=10)
println(ga(state))

# Minimizes exp(|x|)
m = CoordinateModel(x -> exp(dot(x,x)),
                  [-1.,-1.,-1.,-1.,-1], # minimum corner
                  [1.,1.,1.,1.,1]) # maximum corner in rectangle
# Saves creature to file every 10 iterations
state = GAState(m, ngen=100, npop=3_000, elite_fraction=0.01,
                mutation_params=Dict(:rate=>0.1), print_fitness_iter=10,
                save_creature_iter=10, file_name_prefix="minexp_3000")
println(ga(state))

# Change maximum number of generations and then continue
state.ngen = 200
ga(state)

# It's possible to also change the other fields in state and then continue on
# but that is usually only when you made a mistake with a filename or something

if false
    # The following doesn't work b/c JLD doesn't support types w/ functions
    # Minimizes |tan(x + 1000)|
    m = CoordinateModel(x -> abs(tanh(x[1]+1000)), [-2000.], [2000.])
    # Saves full state to file every 10 iterations
    state = GAState(m, ngen=100, npop=3_000, elite_fraction=0.01,
                    mutation_params=Dict(:rate=>0.1), print_fitness_iter=10,
                    save_state_iter=25, file_name_prefix="tan1000_3000")
    println(ga(state))

    # Reload state from file and restart
    state = loadgastate("tan1000_3000_state_25.jld")

    ga(state)

end


text:
using LinearAlgebra
export sample_urn_without_replacement, adj_mat_to_vec, vec_to_adj_mat
export rand_perturb, get_num_vertices

function sample_urn_without_replacement(
    n::Int, bins::Vector{Int}
    )
    var_bins = copy(bins)
    n_tot = sum(bins)
    @assert n_tot ≥ n "Sum of bins must be ≥ n."
    out = Int[]
    for i in 1:n
        counter = 0
        z = rand(1:n_tot)
        for (obs,bin_size) in enumerate(var_bins)
            counter += bin_size
            if counter ≥ z
                var_bins[obs] -= 1
                n_tot -= 1
                push!(out, obs)
                break
            end
        end
    end
    return out
end

function adj_mat_to_vec(
        A::Matrix;
        directed::Bool=true,
        self_loops::Bool=true
    )

    V = size(A,1)
    if directed & self_loops
        return A[:]
    elseif directed & !self_loops
        return A[[i!=j ? true : false for i=1:V,j=1:V]]
    elseif !directed & self_loops
        return A[[i ≤ j ? true : false for i=1:V,j=1:V]]
    elseif !directed & !self_loops
        return A[[i < j ? true : false for i=1:V,j=1:V]]
    end
end

function vec_to_adj_mat(
    x::Vector;
    directed::Bool=true,
    self_loops::Bool=true
    )
    M = length(x)
    if directed & self_loops
        V = round(Int, √(M))
        A = zeros(eltype(x), V, V)
        A[:] = copy(x)
        return A
    elseif directed & !self_loops
        V = round(Int, √(M+1/4) + 1/2)
        A = zeros(eltype(x), V, V)
        A[[i!=j ? true : false for i=1:V,j=1:V]] = copy(x)
        return A
    elseif !directed & self_loops
        V = round(Int, √(2M+1/4) - 1/2)
        A = zeros(eltype(x), V, V)
        A[[i ≤ j ? true : false for i=1:V,j=1:V]] = copy(x)
        A[[i > j ? true : false for i=1:V,j=1:V]] = A'[[i > j ? true : false for i=1:V,j=1:V]]
        @assert issymmetric(A) "Output must be symmetric for un-directed case."
        return A
    elseif !directed & !self_loops
        V = round(Int, √(2M+1/4) + 1/2)
        A = zeros(eltype(x), V, V)
        A[[i < j ? true : false for i=1:V,j=1:V]] = copy(x)
        A += A'
        @assert issymmetric(A) "Output must be symmetric for un-directed case."
        return A
    end

end

rand_perturb(A::Matrix{Bool}, τ::Float64) = [rand() < τ ? !x : x for x in A]

function get_num_vertices(
    M::Int, # Number of edges
    directed::Bool, self_loops::Bool
    )
    if directed & self_loops
        return round(Int, √(M))
    elseif directed & !self_loops
        return round(Int, √(M+1/4) + 1/2)
    elseif !directed & self_loops
        return round(Int, √(2M+1/4) - 1/2)
    elseif !directed & !self_loops
        return round(Int, √(2M+1/4) + 1/2)
    end
end

ext: jl

max_stars_repo_path: src/Graph Models/utils.jl

max_issues_repo_path: src/Graph Models/utils.jl

max_forks_repo_path: src/Graph Models/utils.jl


text:
# ## Scaling

function scaleprojection(ratios::NTuple{N}, T = Float32) where N
    a = zeros(Float32, N, N)
    a[I(N)] = SVector{N}(ratios)
    return LinearMap(SArray{Tuple{N, N}}(a))
end

struct ScaleRatio{N} <: ProjectiveTransform
    ratios::NTuple{N}
end

function getprojection(scale::ScaleRatio, bounds; randstate = nothing)
    return scaleprojection(scale.ratios)
end

"""
    ScaleKeepAspect(minlengths)

Projective transformation that scales the shortest side of `item`
to `minlengths`, keeping the original aspect ratio.
"""
struct ScaleKeepAspect{N} <: ProjectiveTransform
    minlengths::NTuple{N, Int}
end

function getprojection(scale::ScaleKeepAspect{N}, bounds; randstate = nothing) where N
    ratio = maximum(scale.minlengths ./ boundssize(bounds))
    return scaleprojection(Tuple(ratio for _ in 1:N))
end

"""
    ScaleFixed(sizes)

Projective transformation that scales sides to `sizes`, disregarding
aspect ratio.

See also [`ScaleKeepAspect`](#).
"""
struct ScaleFixed{N} <: ProjectiveTransform
    sizes::NTuple{N, Int}
end

function getprojection(scale::ScaleFixed{N}, bounds; randstate = nothing) where N
    ratios = scale.sizes ./ boundssize(bounds)
    return scaleprojection(ratios)
end

"""
    Rotate(γ)

Rotate 2D spatial data by angle γ around the center.
"""
struct Rotate <: ProjectiveTransform
    γ
end

function getprojection(rotate::Rotate, bounds; randstate = nothing)
    middlepoint = sum(bounds) ./ length(bounds)
    r = rotate.γ / 360 * 2pi
    return recenter(RotMatrix(r), middlepoint)
end

"""
    Reflect(γ)

Reflect 2D spatial data by angle γ around the center.
"""
struct Reflect <: ProjectiveTransform
    γ
end

function getprojection(reflect::Reflect, bounds; randstate = nothing)
    middlepoint = sum(bounds) ./ length(bounds)
    r = reflect.γ / 360 * 2pi
    return recenter(reflectionmatrix(r), middlepoint)
end

FlipX() = Reflect(180)
FlipY() = Reflect(90)

reflectionmatrix(r) = SMatrix{2, 2, Float32}(cos(2r), sin(2r), sin(2r), -cos(2r))

"""
    PinOrigin()

Projective transformation that translates the data so that
the upper left bounding corner is at the origin `(0, 0)` (or
the multidimensional equivalent).

Projective transformations on images return `OffsetArray`s,
but not on keypoints. Hardware like GPUs do not support OffsetArrays,
so they will be unwrapped and no longer match up with the keypoints.

Pinning the data to the origin makes sure that the resulting
`OffsetArray` has the same indices as a regular array, starting
at one.
"""
struct PinOrigin <: ProjectiveTransform end

function getprojection(::PinOrigin, bounds; randstate = nothing)
    return Translation(-bounds[1])
end

function apply(::PinOrigin, item::Union{Image, MaskMulti, MaskBinary}; randstate = nothing)
    item = @set item.data = parent(itemdata(item))
    item = @set item.bounds = makebounds(size(itemdata(item)))
    return item
end

# `PinOrigin` should not compose with a cropped transform otherwise the pinning won't work.

compose(cropped::CroppedProjectiveTransform, pin::PinOrigin) = Sequence(cropped, pin)

# ## Resize crops

RandomResizeCrop(sz) = ScaleKeepAspect(sz) |> RandomCrop(sz) |> PinOrigin()
CenterResizeCrop(sz) = ScaleKeepAspect(sz) |> CenterCrop(sz) |> PinOrigin()

ResizePadDivisible(sz, by) = ScaleKeepAspect(sz) |> PadDivisible(by) |> PinOrigin()


text:
# ========================== f_diff_expansion_coeffs_adams_moulton(k) ===========

@doc raw"""
    f_diff_expansion_coeffs_adams_moulton(k [, T=Int])

``k^{th}``-order Adams-Moulton expansion coefficients,

```math
-\frac{\nabla}{ln(1-\nabla)} = \sum_{p=0}^{\infty}b_p\nabla^p= 1 - \frac{1}{2}\nabla - \frac{1}{12}\nabla^2 - \frac{1}{24}\nabla^3 +\cdots.
```
The weights are stored in *forward* order: ``[b_0^k,\ \cdots,\ b_k^k]`` -
order of use in summation.
#### Examples:
```
k = 5
b = f_diff_expansion_coeffs_adams_moulton(k::Int); println(b)
 Rational[1//1, -1//2, -1//12, -1//24, -19//720, -3//160]

D = denominator(gcd(b)); println(D)
 1440

o = convert(Vector{Int},(b .* D)); println(o)
 [1440, -720, -120, -60, -38, -27]
```
"""
function f_diff_expansion_coeffs_adams_moulton(k::Int; T=Int)
# =====================================================================================
#   Adams-Moulton expansion coefficients
# =====================================================================================
    o = Base.zeros(Rational{T},k+1)
    s::Rational{T} = 0//1

    o[1] = 1//1

    b = Base.copy(o)
    a = Base.append!([s],[1//(i+1) for i=1:k])

    for p=1:k
        b = CamiXon.polynom_product_expansion(a, b, k)
        Base.isodd(p) ? o = o .- b : o = o .+ b
    end

    return o  # Note that D = devisor(gcd(o))

end
function f_diff_expansion_coeffs_adams_moulton(k::Int) # short argument: better performance
# =====================================================================================
#   Adams-Moulton expansion coefficients
# =====================================================================================
    o = Base.zeros(Rational{Int},k+1)
    s = 0//1

    o[1] = 1//1

    b = Base.copy(o)
    a = Base.append!([s],[1//(i+1) for i=1:k])

    for p=1:k
        b = CamiXon.polynom_product_expansion(a, b, k)
        Base.isodd(p) ? o = o .- b : o = o .+ b
    end

    return o  # Note that D = devisor(gcd(o))

end

# ========================== create_adams_moulton_weights(k)====================

@doc raw"""
    create_adams_moulton_weights(k::Int [; rationalize=false [, devisor=false [, T=Int]]])

``k^{th}``-order Adams-Moulton weights vector,
```math
y[n+1] = y[n] + \frac{1}{D}\sum_{j=0}^{k}a^k[j]f[n+1-k+j]
```
The weights are stored in the vector ``a^k \equiv[a_k^k/D,\ \cdots,\ a_0^k/D]`` under the convention
``a^k[j] \equiv a_{k-j}^k/D``, where ``a_j^k`` are the Adams-Moulton weight coefficients and ``D`` the
corresponding Adams-Moulton divisor. By default the output is in Float64,
optionally the output is rational, with or without specification of the gcd devisor.
#### Example:
```
[create_adams_moulton_weights(k::Int; rationalize=true, devisor=true, T=Int) for k=1:8]
8-element Vector{Tuple{Int64, Int64, Vector{Int64}}}:
 (1, 2, [1, 1])
 (2, 12, [-1, 8, 5])
 (3, 24, [1, -5, 19, 9])
 (4, 720, [-19, 106, -264, 646, 251])
 (5, 1440, [27, -173, 482, -798, 1427, 475])
 (6, 60480, [-863, 6312, -20211, 37504, -46461, 65112, 19087])
 (7, 120960, [1375, -11351, 41499, -88547, 123133, -121797, 139849, 36799])
 (8, 3628800, [-33953, 312874, -1291214, 3146338, -5033120, 5595358, -4604594, 4467094, 1070017])
```
"""
function create_adams_moulton_weights(k::Int; rationalize=false, devisor=false, T=Int)

    ∇ = CamiXon.f_diff_weights_array(k)
    coeffs = CamiXon.f_diff_expansion_coeffs_adams_moulton(k; T)

    o = CamiXon.f_diff_expansion_weights(coeffs,∇)

    if rationalize
        D = Base.denominator(Base.gcd(o))       # Adams-Moulton devisor
        o = devisor ? (k, D, Base.round.(T, o * D)) : o
    else
        o = Base.convert(Vector{Float64},o)
    end

    return o

end

# ======================== f_diff_expansion_coeffs_adams_bashford(k) ===========

@doc raw"""
    f_diff_expansion_coeffs_adams_bashford(k [, T=Int])

``(k+1)``-point Adams-Bashford expansion coefficients ``B_p``.

```math
-\frac{\nabla}{(1-\nabla)ln(1-\nabla)}=\sum_{p=0}^{\infty}B_p\nabla^p=1+\ \frac{1}{2}∇+\ \frac{5}{12}∇^2+\ \cdots.
```
The weights are stored in *forward* order: ``[B_0^k,\ \cdots,\ B_k^k]`` -
order of use in summation.
#### Examples:
```
k = 5
o = f_diff_expansion_coeffs_adams_bashford(k::Int); println(o)
 Rational{Int64}[1//1, 1//2, 5//12, 3//8, 251//720, 95//288]
```
"""
function f_diff_expansion_coeffs_adams_bashford(k::Int; T=Int)
# =====================================================================================
#   Adams-Bashford expansion coefficients
# =====================================================================================
    a = Base.ones(Rational{T},k+1)

    b = CamiXon.f_diff_expansion_coeffs_adams_moulton(k; T)
    o = CamiXon.polynom_product_expansion(a, b, k)

    return o  # Note that D = denominator(gcd(o))

end
function f_diff_expansion_coeffs_adams_bashford(k::Int) # short argument: better performance

    a = Base.ones(Rational{Int},k+1)

    b = CamiXon.f_diff_expansion_coeffs_adams_moulton(k)
    o = CamiXon.polynom_product_expansion(a, b, k)

    return o  # Note that D = denominator(gcd(o))

end


text:
# Iterator
# --------

function Base.eltype(::HTSReadWriter)
    HTSRecord
end

function Base.IteratorSize(::Type{T}) where T<: HTSReadWriter
    Base.SizeUnknown()
end

function Base.iterate(hf::HTSReadWriter, state::Missing)
    record = HTSRecord()
    isnothing(tryread!(hf, record)) && return nothing
    return (record, state)
end

function Base.iterate(hf::HTSReadWriter)
    state = missing
    iterate(hf, state)
end

#function eachrecord(hf::HTSReadWriter)
#end

# MultiRegion Iterator
# --------------------

export HTSRegionsIterator

"""
Iterate alignment records over the specified regions. If a read
overlaps more than one region, it will only be returned once.

The region string can take one of the following forms:

| region          | Outputs                                                     |
|:--------------- |:----------------------------------------------------------- |
| REF             | All reads with RNAME REF                                    |
| REF:            | All reads with RNAME REF                                    |
| REF:START       | Reads with RNAME REF overlapping START to end of REF        |
| REF:-END        | Reads with RNAME REF overlapping start of REF to END        |
| REF:START-END   | Reads with RNAME REF overlapping START to END               |
| .               | All reads from the start of the file                        |
| *               | Unmapped reads at the end of the file (RNAME '*' in SAM)    |

"""
HTSRegionsIterator

mutable struct HTSRegionsIterator{T}
    htsreader::HTSReadWriter{T}
    itr_ptr::Ptr{htslib.hts_itr_t}
    function HTSRegionsIterator(hf::HTSReadWriter{T}, regions::AbstractVector) where T
        if !isa(regions, Vector{String})
            return HTSRegionsIterator(hf, map(r -> region_string(hf, r), regions))
        end
        # also need to preserve pointer array from gc
        @assert isa(regions, Vector{String})
        regarray = pointer.(regions)
        GC.@preserve regions regarray begin
            idxptr = unsafe_load(pointer(hf)).idx
            hdrptr = pointer(header(hf))
            regarrayptr = convert(Ptr{Ptr{Char}}, pointer(regarray))
            itrptr = htslib.sam_itr_regarray(idxptr, hdrptr, regarrayptr, Cuint(length(regions)))
            ans = new{T}(hf, itrptr)
            Base.finalizer(ans) do ans
                ptr = pointer(ans)
                ptr != C_NULL && htslib.sam_itr_destroy(ptr)
                return nothing
            end
            ans
        end
    end
end

function HTSRegionsIterator(hf::HTSReadWriter, chr, start, stop)
    regions = String[]
    if isa(chr, AbstractString)
        chr = Ref(chr)
    end
    for (c, s, e) in zip(chr, start, stop)
        region = region_string(hf, Tuple{typeof(c),typeof(s),typeof(e)}((c, s, e)))
        push!(regions, region)
    end
    HTSRegionsIterator(hf, regions)
end

function HTSRegionsIterator(hf::HTSReadWriter, region::Union{Tuple, AbstractString})
    regions = Vector{typeof(region)}(undef, 1)
    regions[1] = region
    return HTSRegionsIterator(hf, regions)
end

# bad solution, this should be a feature request to htslib
function region_string(hf::HTSReadWriter, regionspec::Tuple{<:AbstractString,<:Real,<:Real})::String
    (chr, start, stop) = regionspec
    @assert start >= 0
    @assert stop >= 0
    start = Int(start)
    stop = Int(stop)
    return "$chr:$start-$stop"
end

function region_string(hf::HTSReadWriter, regionspec::Tuple{<:Real,<:Real,<:Real})::String
    (chr, start, stop) = regionspec
    @assert chr >= 1
    @assert start >= 0
    @assert stop >= 0
    chr_str = refname(header(hf), chr)
    start = Int(start)
    stop = Int(stop)
    return "$chr_str:$start-$stop"
end

function region_string(hf::HTSReadWriter, regionspec::AbstractString)::String
    String(regionspec)
end

function BioGenerics.IO.tryread!(itr::HTSRegionsIterator, record::HTSRecord)::Union{HTSRecord,Nothing}
    GC.@preserve itr begin
        res = htslib.sam_itr_next(pointer(parent(itr)), pointer(itr), pointer(record))
        res < -1 && error("error reading record")
        res == -1 && return nothing
        return record
    end
end

function Base.pointer(x::HTSRegionsIterator)
    x.itr_ptr
end

function Base.parent(x::HTSRegionsIterator)
    x.htsreader
end

function Base.eltype(::HTSRegionsIterator)
    HTSRecord
end

function Base.IteratorSize(::Type{T}) where T <: HTSRegionsIterator
    Base.SizeUnknown()
end

function Base.iterate(itr::HTSRegionsIterator, state::Missing)
    res = tryread!(itr, HTSRecord())
    isnothing(res) && return nothing
    return (res, missing)
end

function Base.iterate(itr::HTSRegionsIterator)
    iterate(itr, missing)
end

## implementation of read! and read are common between HTSRegionsIterator and HTSReadWriter
function Base.read!(hf::HTSRegionsIterator, record::HTSRecord)
    isnothing(tryread!(hf, record)) && throw(EOFError())
    record
end

function Base.read(hf::HTSRegionsIterator)
    read!(hf, HTSRecord())
end


text:
## test data available for use in PopGen
export dataset, @nancycats, @gulfsharks

"""
    dataset(::String)
Load an example dataset from either `"gulfsharks"` (SNP) or `"nancycats"` (microsatellite).
Can also use `"sharks"` and `"cats"` as shorthands. Use `?@nancycats` and
`?@gulfsharks` to learn more about these datasets.

### Example
```
ncats = dataset("nancycats")
gsharks = dataset("sharks")
```
"""
function dataset(name::String)
    if lowercase(name) in ["nancycats", "cats"]
        filename = normpath(joinpath(@__DIR__,"../..","data", "nancycats.gen"))
    elseif lowercase(name) in ["gulfsharks", "sharks"]
        filename = normpath(joinpath(@__DIR__,"../..","data", "gulfsharks.csv"))
    else
        throw(ArgumentError("Please choose either the \"nancycats\" or \"gulfsharks\" datasets"))
    end
    read_from(filename, silent = true)
end

"""
    @nancycats
Returns `PopData` of corresponding "nancycats" dataset as featured in
the R package `adegenet`. This dataset is composed of 9 microsatellite
loci in 237 individuals across 17 populations.

Example:
```
ncats = @nancycats
```
"""
macro nancycats()
    return esc(quote
       dataset("nancycats")
    end)
end

"""
    @gulfsharks
Returns `PopData` corresponding the Blacknose shark dataset as used in
Dimens et al. 2019. This is a mid-sized SNP dataset of 2209 haplotypes
in 212 individuals, across 7 populations.

Example:

```
sharks = @gulfsharks
```
"""
macro gulfsharks()
    return esc(quote
       dataset("gulfsharks")
    end)
end


text:
# Test layers and data/model movements on and off the GPU
# Add tests for layers and their gradients on the GPU
# Most of the forward passes should be fine being applied
# to bitstype objects, but this gives higher coverage for our use-cases
# Check that getting the gradients does not throw

# generic movement tests
@testset "Basic GPU Movement" begin
  @test gradient(x -> sum(gpu(x)), rand(3,3)) isa Tuple
  @test gradient(x -> sum(cpu(x)), gpu(rand(3,3))) isa Tuple
end

# TODO: These layers get into scalar indexing issues.
const BROKEN_LAYERS = Union{DepthwiseConv}

const ACTIVATIONS = [identity, relu, tanh,
                     sigmoid, exp, softplus,
                     elu, selu]

function gpu_gradtest(name::String, layers::Vector, x_cpu = nothing, args...; test_cpu = true)
  isnothing(x_cpu) && error("Missing input to test the layers against.")
  @testset "$name GPU grad tests" begin
    for layer in layers
      @testset "$layer Layer GPU grad test" begin

        # compute output and grad of parameters
        l_cpu = layer(args...)
        ps_cpu = Flux.params(l_cpu)
        y_cpu, back_cpu = pullback(() -> sum(l_cpu(x_cpu)), ps_cpu)
        gs_cpu = back_cpu(1f0)

        x_gpu = gpu(x_cpu)
        l_gpu = l_cpu |> gpu
        ps_gpu = Flux.params(l_gpu)

        if typeof(l_gpu) <: BROKEN_LAYERS
          @test_broken gradient(() -> sum(l_gpu(x_gpu)), ps_gpu) isa Flux.Zygote.Grads
        else
          y_gpu, back_gpu = pullback(() -> sum(l_gpu(x_gpu)), ps_gpu)
          gs_gpu = back_gpu(1f0) # TODO many layers error out when backprop int 1, should fix

          # compute grad of input
          xg_cpu = gradient(x -> sum(l_cpu(x)), x_cpu)[1]
          xg_gpu = gradient(x -> sum(l_gpu(x)), x_gpu)[1]

          # test
          if test_cpu
            if VERSION >= v"1.7" && layer === GroupedConvTranspose && args[end] == selu
              @test_broken y_gpu ≈ y_cpu rtol=1f-3 atol=1f-3
            else
              @test y_gpu ≈ y_cpu rtol=1f-3 atol=1f-3
            end
            if isnothing(xg_cpu)
              @test isnothing(xg_gpu)
            else
              if layer === GroupedConvTranspose
                @test Array(xg_gpu) ≈ xg_cpu rtol = 2f-2 atol = 1f-3
              else
                @test Array(xg_gpu) ≈ xg_cpu rtol = 1f-3 atol = 1f-3
              end
            end
          end
          @test gs_gpu isa Flux.Zygote.Grads
          for (p_cpu, p_gpu) in zip(ps_cpu, ps_gpu)
            if isnothing(gs_cpu[p_cpu])
              @test isnothing(gs_gpu[p_gpu])
            else
              @test gs_gpu[p_gpu] isa Flux.CUDA.CuArray
              if test_cpu
                @test Array(gs_gpu[p_gpu]) ≈ gs_cpu[p_cpu] rtol=1f-3 atol=1f-3
              end
            end
          end
        end
      end
    end
  end
end

# Just to give testset in gpu_gradtest meaningful labels
ConvNoBias(args...) = Conv(args...; bias = false)
ConvTransposeNoBias(args...) = ConvTranspose(args...; bias = false)
CrossCorNoBias(args...) = CrossCor(args...; bias = false)
DepthwiseConvNoBias(args...) = DepthwiseConv(args...; bias = false)
GroupedConv(args...) = Conv(args..., groups = 5)
GroupedConvTranspose(args...) = ConvTranspose(args..., groups = 5)

for act in ACTIVATIONS
  r = rand(Float32, 28, 28, 1, 1)
  conv_layers = [Conv, ConvNoBias,
                 ConvTranspose, ConvTransposeNoBias,
                 CrossCor, CrossCorNoBias,
                 DepthwiseConv, DepthwiseConvNoBias]
  gpu_gradtest("Convolution with $act", conv_layers, r, (2,2), 1=>3, act, test_cpu = false)

  groupedconv = [GroupedConv, GroupedConvTranspose]
  gpu_gradtest("GroupedConvolution with $act", groupedconv, rand(Float32, 28, 28, 100, 2), (3,3), 100 => 25, act, test_cpu = true)

  batch_norm = [BatchNorm]
  gpu_gradtest("BatchNorm 1 with $act", batch_norm, rand(Float32, 28,28,3,4), 3, act, test_cpu = false) #TODO fix errors
  gpu_gradtest("BatchNorm 2 with $act", batch_norm, rand(Float32, 5,4), 5, act, test_cpu = false)

  instancenorm = [InstanceNorm]
  gpu_gradtest("InstanceNorm with $act", instancenorm, r, 1, act, test_cpu = false)

  groupnorm = [GroupNorm]
  gpu_gradtest("GroupNorm with $act", groupnorm, rand(Float32, 28,28,3,1), 3, 1, act, test_cpu = false)
end

r = rand(Float32, 28, 28, 1, 1)

pooling_layers = [MaxPool, MeanPool]
gpu_gradtest("Pooling", pooling_layers, r, (2,2))

adaptive_pooling_layers = [AdaptiveMaxPool, AdaptiveMeanPool]
gpu_gradtest("AdaptivePooling", adaptive_pooling_layers, r, (7,7), test_cpu = false)

dropout_layers = [Dropout, AlphaDropout]
gpu_gradtest("Dropout", dropout_layers, r, 0.5f0; test_cpu = false) # dropout is not deterministic

layer_norm = [LayerNorm]
gpu_gradtest("LayerNorm 1", layer_norm, rand(Float32, 28,28,3,4), 1, test_cpu = false) #TODO fix errors
gpu_gradtest("LayerNorm 2", layer_norm, rand(Float32, 5,4), 5)

upsample = [x -> Upsample(scale=x)]
gpu_gradtest("Upsample 2d", upsample, rand(Float32, 3, 4, 2, 3), (2,2))
gpu_gradtest("Upsample 1d", upsample, rand(Float32, 3, 4, 2, 3), (2,))

pixelshuffle = [PixelShuffle]
gpu_gradtest("PixelShuffle 2d", pixelshuffle, rand(Float32, 3, 4, 18, 3), 3)
gpu_gradtest("PixelShuffle 1d", pixelshuffle, rand(Float32, 3, 18, 3), 3)

embedding = [Flux.Embedding]
gpu_gradtest("Embedding", embedding, [1,3,5], 5, 2)
gpu_gradtest("Embedding repeated indices", embedding, [1,3,5,3], 5, 2)
gpu_gradtest("Embedding integer index", embedding, 1, 5, 2)
gpu_gradtest("Embedding 2d index", embedding, [1 2; 3 4], 5, 2)
gpu_gradtest("Embedding OneHotVec index", embedding, OneHotVector(1, 5), 5, 2)
gpu_gradtest("Embedding OneHotMatrix index", embedding,  OneHotMatrix([1,2,3], 5), 5, 2)
gpu_gradtest("Embedding OneHotMatrix repeated indices", embedding, OneHotMatrix([1,2,2], 5), 5, 2)

@testset "function layers" begin
  x = rand(Float32, 3,3)
  gpu_autodiff_test(x -> sum(Flux.normalise(x; dims=1)), x)
  gpu_autodiff_test(x -> sum(Flux.normalise(x; dims=2)), x)
  gpu_autodiff_test(x -> sum(Flux.normalise(x)), x)
end

@testset "Zeros mapped for $cl" for cl in (Conv, ConvTranspose, CrossCor, DepthwiseConv)
  l = cl((2,2), 1=>3, bias = false) |> gpu
  ip = zeros(Float32, 28,28,1,1) |> gpu
  if typeof(l) <: BROKEN_LAYERS
    @test_broken sum(l(ip)) ≈ 0.f0
    @test_broken gradient(() -> sum(l(ip)), Flux.params(l)) isa Flux.Zygote.Grads
  else
    @test sum(l(ip)) ≈ 0.f0
    gs = gradient(() -> sum(l(ip)), Flux.params(l))
    @test l.bias ∉ gs.params
  end
end

@testset "Dense with Zeros bias" begin
  l = Dense(ones(Float32, 4, 3), Flux.Zeros()) |> gpu
  ip = zeros(Float32, 3, 7) |> gpu

  @test sum(l(ip)) ≈ 0.f0
  gs = gradient(() -> sum(l(ip)), Flux.params(l))
  @test l.bias ∉ gs.params
end

@testset "Extended BatchNorm" begin
  m_cpu = BatchNorm(2)
  m_gpu = m_cpu |> gpu
  x_cpu = rand(Float32, 3, 2, 2)
  x_gpu = x_cpu |> gpu

  ## In :auto mode, track statistics only in gradient contest
  μ_cpu = copy(m_cpu.μ)
  m_cpu(x_cpu)
  @test m_cpu.μ ≈ μ_cpu
  gradient(() -> sum(m_cpu(x_cpu)), Flux.params(m_cpu))
  @test !(m_cpu.μ ≈ μ_cpu)

  μ_gpu = copy(m_gpu.μ)
  m_gpu(x_gpu)
  @test m_gpu.μ ≈ μ_gpu
  gradient(() -> sum(m_gpu(x_gpu)), Flux.params(m_gpu))
  @test !(m_gpu.μ ≈ μ_gpu)

  @test Array(m_gpu.μ) ≈ m_cpu.μ

  ## In testmode, never track statistics
  testmode!(m_cpu)
  μ_cpu = copy(m_cpu.μ)
  m_cpu(x_cpu)
  @test m_cpu.μ ≈ μ_cpu
  gradient(() -> sum(m_cpu(x_cpu)), Flux.params(m_cpu))
  @test m_cpu.μ ≈ μ_cpu

  testmode!(m_gpu)
  μ_gpu = copy(m_gpu.μ)
  m_gpu(x_gpu)
  @test m_gpu.μ ≈ μ_gpu
  gradient(() -> sum(m_gpu(x_gpu)), Flux.params(m_gpu))
  @test m_gpu.μ ≈ μ_gpu

  ## In trainmode, always track statistics
  trainmode!(m_cpu)
  μ_cpu = copy(m_cpu.μ)
  m_cpu(x_cpu)
  @test !(m_cpu.μ ≈ μ_cpu)
  μ_cpu = copy(m_cpu.μ)
  gradient(() -> sum(m_cpu(x_cpu)), Flux.params(m_cpu))
  @test !(m_cpu.μ ≈ μ_cpu)

  trainmode!(m_gpu)
  μ_gpu = copy(m_gpu.μ)
  m_gpu(x_gpu)
  @test !(m_gpu.μ ≈ μ_gpu)
  μ_gpu = copy(m_gpu.μ)
  gradient(() -> sum(m_gpu(x_gpu)), Flux.params(m_gpu))
  @test !(m_gpu.μ ≈ μ_gpu)

  ## No errors if input type mistmatch
  # x_cpu = rand(Float64, 3, 2, 2)
  # x_gpu = x_cpu |> gpu
  # m_cpu(x_cpu)
  # gradient(() -> sum(m_cpu(x_cpu)), Flux.params(m_cpu))
  # m_gpu(x_gpu)
  # gradient(() -> sum(m_gpu(x_gpu)), Flux.params(m_gpu))
end

@testset "Two-streams Bilinear" begin
  x = zeros(Float32,10,9) |> gpu
  y = zeros(Float32,2,9) |> gpu
  b = Flux.Bilinear(10, 2, 3) |> gpu
  @test size(b(x,y)) == (3,9)
  @test sum(abs2, b(x,y)) ≈ 0f0
  gs_gpu = gradient(() -> sum(abs2.(b(x, y))), params(b))
  b_cpu, x_cpu, y_cpu = b |> cpu, x |> cpu, y |> cpu
  gs_cpu = gradient(() -> sum(abs2.(b_cpu(x_cpu, y_cpu))), params(b_cpu))
  for (pgpu, pcpu) in zip(params(b), params(b_cpu))
    @test gs_cpu[pcpu] ≈ Array(gs_gpu[pgpu])
  end
end

@testset "Two-streams Bilinear" begin
  x = zeros(Float32,10,9) |> gpu
  y = zeros(Float32,2,9) |> gpu
  b = Flux.Bilinear(10, 2, 3) |> gpu
  @test size(b(x,y)) == (3,9)
  @test sum(abs2, b(x,y)) ≈ 0f0
  gs_gpu = gradient(() -> sum(abs2.(b(x, y))), params(b))
  b_cpu, x_cpu, y_cpu = b |> cpu, x |> cpu, y |> cpu
  gs_cpu = gradient(() -> sum(abs2.(b_cpu(x_cpu, y_cpu))), params(b_cpu))
  for (pgpu, pcpu) in zip(params(b), params(b_cpu))
    @test gs_cpu[pcpu] ≈ Array(gs_gpu[pgpu])
  end
end

@testset "Parallel" begin
  @testset "zero sum" begin
    input = randn(10, 10, 10, 10) |> gpu
    layer_gpu = Parallel(+, zero, identity) |> gpu
    @test layer_gpu(input) == input
    @test layer_gpu(input) isa Flux.CUDA.CuArray
  end

  @testset "vararg input" begin
    inputs = (randn(10), randn(5), randn(4)) .|> gpu
    layer = Parallel(+, Dense(10, 2), Dense(5, 2), Dense(4, 2)) |> gpu
    @test size(layer(inputs)) == (2,)
  end

  @testset "gradient" begin
    input_cpu = randn(10, 10, 10, 10)
    input_gpu = input_cpu |> gpu
    layer_cpu = Parallel(+, x -> zero(x), identity)
    layer_gpu = layer_cpu |> gpu
    gs_cpu = gradient(() -> sum(abs2.(layer_cpu(input_cpu))), params(layer_cpu))
    gs_gpu = gradient(() -> sum(abs2.(layer_gpu(input_gpu))), params(layer_gpu))
    for (pgpu, pcpu) in zip(params(layer_cpu), params(layer_gpu))
      @test gs_cpu[pcpu] ≈ gs_gpu[pgpu]
    end
  end
end


text:
# setup.jl

@static if Sys.iswindows()

    # TODO better ways of finding the Stata binary
    function locate_executable()

        # 1.) Environment Variable
        if haskey(ENV,"STATA_BIN")
            return ENV["STATA_BIN"]
        end

        # if we're still here, we have not found the executable
        error("Could not find the Stata executable. Please set the
        \"STATA_BIN\" environment variable.")

    end

    stata_executable = locate_executable()

end

@static if Sys.isapple()

    function locate_executable()

        # try a few common directories
        if isfile("/usr/bin/stata")
            return "/usr/bin/stata"
       elseif isfile("/usr/local/bin/stata-se")
            return "/usr/local/bin/stata-se"
       elseif isfile("/usr/local/bin/stata-mp")
            return "/usr/local/bin/stata-mp"
        elseif isdir("/Applications/Stata/StataSE.app")
            appdir = "/Applications/Stata/StataSE.app"
        elseif isdir("/Applications/Stata/StataMP.app")
            appdir = "/Applications/Stata/StataMP.app"
        elseif isdir("/Applications/Stata/StataIC.app")
            appdir = "/Applications/Stata/StataIC.app"
        elseif isdir("/Applications/StataSE.app")
            appdir = "/Applications/StataSE.app"
        elseif isdir("/Applications/StataMP.app")
            appdir = "/Applications/StataMP.app"
        elseif isdir("/Applications/StataIC.app")
            appdir = "/Applications/StataIC.app"
        else
            appdir = ""
        end

        if (appdir != "") && (isfile(joinpath(appdir, "Contents", "MacOS", "StataMP")))
            return joinpath(appdir, "Contents", "MacOS", "StataMP")
        else

            # try environment variables
            if haskey(ENV,"STATA_BIN")
                return ENV["STATA_BIN"]
            end

            # if we're still here, we have not found the executable
            error("Could not find the Stata executable. Please set the
                \"STATA_BIN\" environment variable or create a symlink
                in /usr/bin/stata that points to the executable.")

        end

    end

    stata_executable = locate_executable()

end

@static if Sys.islinux()

    # TODO

end


text:
# ==============================================================================
# Wanda Feng
# last updated: March 2021
#
# Rabbit Food Nutritional Data
#
# ==============================================================================

# Packages =====================================================================
using DataFrames, DataFramesMeta, CSV
using Unicode
using Statistics

# Start timer ==================================================================
#@time begin

# Directories ==================================================================
dirdata = "/Users/wandafeng/Desktop/Rabbit_Food_data/"
USDA_dirdata = dirdata*"FoodData_Central_csv_2020-10-30/"

# Read files ===================================================================
# USDA all data ----------------------------------------------------------------
#       downloaded from: https://fdc.nal.usda.gov/download-datasets.html
#   food names & IDs
file1 = USDA_dirdata*"food.csv"
USDA_Food = DataFrame(CSV.File(file1)) # 533612x5
#   portion weights (g) & equivalents
file2 = USDA_dirdata*"food_portion.csv"
USDA_PortionData = DataFrame(CSV.File(file2)) # 56893x11
#   nutrients by ID
file3 = USDA_dirdata*"food_nutrient.csv"
USDA_NutrientData = DataFrame(CSV.File(file3)) # 7547844x11

# Isolate relevant data types ==================================================
#   FNDDS, Foundation, & SR Legacy contain relevant veg & fruit data
###### may want to add Branded Foods in future - e.g. "Spring mix"
USDA_Food_ed = @where(USDA_Food,
    in(["survey_fndds_food","foundation_food","sr_legacy_food"]).(:data_type)) # 15053x5

# Tabulate portion weights for rabbit safe veg & fruit =========================
#   combine food names, USDA IDs, + portion data -------------------------------
USDA_Food_Portion = innerjoin(USDA_Food_ed,USDA_PortionData, on = :fdc_id) # 46402x15
#CSV.write(dirdata*"USDA_Food_Portion.csv", USDA_Food_Portion)
USDA_Food_Portion_gr = groupby(USDA_Food_Portion, :data_type)
#   isolate foods where portion = 1 cup ----------------------------------------
FNDDS_FP = @where(USDA_Food_Portion_gr[2], in(["1 cup"]).(:portion_description)) # 3305x15
Foundation_FP = @where(USDA_Food_Portion_gr[1], in([1000]).(:measure_unit_id)) # 51x15
SRlegacy_FP = @where(USDA_Food_Portion_gr[3], in(["cup","cup, chopped",
    "cup shredded","cup, shredded","cup sliced","cup, sliced",
    "cup sprigs"]).(:modifier)) # 1880x15
#   search for food item by data type
function FindItemData(food_item, db)
    println("============================================")
    println("Finding portion values for "*food_item*" in "*db.data_type[1])
    #   make string comparison case insensitive
    food_item = Unicode.normalize(food_item*",", casefold=true)
    for j = 1:length(db.description)
        db_item = Unicode.normalize(db.description[j], casefold=true)
        if occursin(food_item,db_item) == true #&&
            #endswith(db.description[j],"raw") == true
            println("--- description: "*db.description[j])
        end
    end
end

#=
#   isolate FNDDS veg & fruit --------------------------------------------------
food_cat_id_begin = zeros(length(FNDDS_FP.food_category_id))
for i=1:length(FNDDS_FP.food_category_id)
    if startswith(string(FNDDS_FP.food_category_id[i]),"6") == true
        #println("index = "*string(i)*"--------------------------------------------")
        #println(string(FNDDS_FP.food_category_id[i]))
        food_cat_id_begin[i] = 1
    end
end
FNDDS_FP = hcat(FNDDS_FP,food_cat_id_begin)
FNDDS_FP = @where(FNDDS_FP, in([1.0]).(:x1)) #713x16
select!(FNDDS_FP,Not(:x1)) #713x15
=#
#   combine ed data ------------------------------------------------------------
USDA_Food_Portion = vcat(FNDDS_FP,Foundation_FP,SRlegacy_FP)

#   Rabbit safe veg & fruit lists ----------------------------------------------
#       Leafy Greens; Herbs & Flowers; Non Leafy; Fruit
LG_list = ["Arugula","Beet greens","Cabbage, Chinese","Carrot greens","Chicory greens",
    "Collards","Dandelion greens","Endive","Escarole","Kale",
    "Boston","Green leaf","Red leaf","Romaine","Mustard greens",
    "Radicchio","Radish greens","Spinach","Chard","Turnip greens","Watercress"]
    # Chinese Cabbage = Bok Choy
H_list = ["Basil","Caraway","Chamomile","Cilantro","Clover",
    "Dill weed","Echinacea","Goldenseal","Kava","Lavender",
    "Lemon balm","Lemongrass","Marigolds","Marjoram","Peppermint",
    "Spearmint","Nasturtium","Nettle","Oregano","Pansy",
    "Parsley","Purslane","Rose","Rosemary","Sage",
    "Shepherd's purse","Tarragon","Thyme","Wheatgrass"]
### FILTER FOR TSP & TBSP
NL_list = ["Beets","Pepper, sweet, green","Pepper, sweet, red", "Broccoli","Brussels",
    "Cabbage, green","Cabbage, red","Carrots","Cauliflower","Celery",
    "Fennel","Summer squash, green","Summer squash, yellow","zucchini"]
F_list = ["Apple","Apricot","Banana","Blackberries","Blueberries",
    "Cherries","Kiwi fruit","Mango","Nectarine","Papaya",
    "Peach","Pear","Pineapple","Plum","Raspberries",
    "Strawberries","Watermelon"]
Master_list = vcat(LG_list,H_list,NL_list,F_list) # 81
#   Find portion data for food list --------------------------------------------
function FindPortionData(food_list, verbose)
    #   empty dataframe - don't make df with 0s b/c multiple data vals for food items
    Food_PortionWts = DataFrame(food=[],data_type=[],fdc_id=[],description=[],gram_weight=[])
    for i = 1:length(food_list)
        if verbose == true
            println("============================================")
            println("Finding portion values for "*food_list[i]*"---------")
        end
        #   make string comparison case insensitive
        food_list_item = Unicode.normalize(food_list[i]*",", casefold=true)
        for j = 1:length(USDA_Food_Portion.description)
            USDA_Food_Portion_item = Unicode.normalize(USDA_Food_Portion.description[j], casefold=true)
            if occursin(food_list_item,USDA_Food_Portion_item) == true &&
                (endswith(USDA_Food_Portion.description[j],"raw") == true ||
                endswith(USDA_Food_Portion.description[j],"fresh") == true)
                # also ends with "fresh"
                if verbose == true
                    println("match in "*USDA_Food_Portion.data_type[j])
                end
                if USDA_Food_Portion.data_type[j] == "survey_fndds_food"
                    if verbose == true
                        println("--- description: "*USDA_Food_Portion.description[j])
                        println("--- g/cup = "*string(USDA_Food_Portion.gram_weight[j]))
                    end
                    push!(Food_PortionWts, [food_list[i], USDA_Food_Portion.data_type[j],
                        USDA_Food_Portion.fdc_id[j], USDA_Food_Portion.description[j],
                        USDA_Food_Portion.gram_weight[j]])
                elseif USDA_Food_Portion.data_type[j] == "foundation_food"
                    if verbose == true
                        println("--- description: "*USDA_Food_Portion.description[j])
                        println("--- "*string(USDA_Food_Portion.gram_weight[j])*"g/ "*
                            string(USDA_Food_Portion.amount[j])*" cup")
                    end
                    if USDA_Food_Portion.amount[j] == 1.0
                        push!(Food_PortionWts, [food_list[i], USDA_Food_Portion.data_type[j],
                            USDA_Food_Portion.fdc_id[j], USDA_Food_Portion.description[j],
                            USDA_Food_Portion.gram_weight[j]])
                    elseif USDA_Food_Portion.amount[j] == 0.5
                        push!(Food_PortionWts, [food_list[i], USDA_Food_Portion.data_type[j],
                            USDA_Food_Portion.fdc_id[j], USDA_Food_Portion.description[j],
                            2*USDA_Food_Portion.gram_weight[j]])
                    end
                elseif USDA_Food_Portion.data_type[j] == "sr_legacy_food"
                    if verbose == true
                        println("--- description: "*USDA_Food_Portion.description[j])
                        println("--- "*string(USDA_Food_Portion.gram_weight[j])*"g/ "*
                            string(USDA_Food_Portion.amount[j])*" cup")
                    end
                    if USDA_Food_Portion.amount[j] == 1.0
                        push!(Food_PortionWts, [food_list[i], USDA_Food_Portion.data_type[j],
                            USDA_Food_Portion.fdc_id[j], USDA_Food_Portion.description[j],
                            USDA_Food_Portion.gram_weight[j]])
                    elseif USDA_Food_Portion.amount[j] == 0.5
                        push!(Food_PortionWts, [food_list[i], USDA_Food_Portion.data_type[j],
                            USDA_Food_Portion.fdc_id[j], USDA_Food_Portion.description[j],
                            2*USDA_Food_Portion.gram_weight[j]])
                    end
                end
            end
        end
    end
    return Food_PortionWts
end

Rabbit_Food_PortionWts = FindPortionData(Master_list, true) # 84x5

#   remove outliers - e.g. Java plum result for Plum ---------------------------

filter(row -> !(row.description =="Java-plum, (jambolan), raw"), Rabbit_Food_PortionWts)

# Tabulate nutrient data [per 100g portion] for rabbit safe food ===============
#   combine portion weights + nutrient data =-----------------------------------
Rabbit_FP_Nutrient = innerjoin(Rabbit_Food_PortionWts,USDA_NutrientData, on = :fdc_id) # 5759x15
#   isolate nutrients of interest ----------------------------------------------
#       nutrient codes in "nutrient.csv"
#       Water [g] = 1051; Energy [kcal] = 1008; Protein [g] = 1003; Fat [g] = 1004;
#       Carb [g] = 1005; Fiber [g] = 1079; Sugars [g] = 2000; starch [g] = 1009;
#       Calcium [mg] = 1087
Rabbit_FP_Nutrient_ed = @where(Rabbit_FP_Nutrient,
    in([1051,1008,1003,1004,1005,1079,2000,1009,1087]).(:nutrient_id)) # 633x15
Rabbit_FP_Nutrient_gr = groupby(Rabbit_FP_Nutrient_ed,:fdc_id)
#   tabulate nutrients of interest by USDA ID ----------------------------------
Nutrient_Data = DataFrame(fdc_id=[],water=[],energy=[],protein=[],fat=[],
    carb=[],fiber=[],sugars=[],starch=[],Ca=[])
for i=1:length(Rabbit_FP_Nutrient_gr)
    println("============================================")
    println("Finding nutrient values for "*Rabbit_FP_Nutrient_gr[i].food[1]*"---------")
    #   find indices of nutrients
    water_i = first(indexin(1051,Rabbit_FP_Nutrient_gr[i].nutrient_id))
    energy_i = first(indexin(1008,Rabbit_FP_Nutrient_gr[i].nutrient_id))
    protein_i = first(indexin(1003,Rabbit_FP_Nutrient_gr[i].nutrient_id))
    fat_i = first(indexin(1004,Rabbit_FP_Nutrient_gr[i].nutrient_id))
    carb_i = first(indexin(1005,Rabbit_FP_Nutrient_gr[i].nutrient_id))
    fiber_i = first(indexin(1079,Rabbit_FP_Nutrient_gr[i].nutrient_id))
    sugar_i = first(indexin(2000,Rabbit_FP_Nutrient_gr[i].nutrient_id))
    starch_i = first(indexin(1009,Rabbit_FP_Nutrient_gr[i].nutrient_id))
    calcium_i = first(indexin(1087,Rabbit_FP_Nutrient_gr[i].nutrient_id))
    #=
    water_i = first(findall(x -> x == 1051,Rabbit_FP_Nutrient_gr[i].nutrient_id))
    energy_i = first(findall(x -> x == 1008,Rabbit_FP_Nutrient_gr[i].nutrient_id))
    protein_i = first(findall(x -> x == 1003,Rabbit_FP_Nutrient_gr[i].nutrient_id))
    fat_i = first(findall(x -> x == 1004,Rabbit_FP_Nutrient_gr[i].nutrient_id))
    carb_i = first(findall(x -> x == 1005,Rabbit_FP_Nutrient_gr[i].nutrient_id))
    fiber_i = first(findall(x -> x == 1079,Rabbit_FP_Nutrient_gr[i].nutrient_id))
    sugar_i = first(findall(x -> x == 2000,Rabbit_FP_Nutrient_gr[i].nutrient_id))
    calcium_i = first(findall(x -> x == 1087,Rabbit_FP_Nutrient_gr[i].nutrient_id))
    =#
    #   nutrient values
    water = Rabbit_FP_Nutrient_gr[i].amount[water_i]
    energy = Rabbit_FP_Nutrient_gr[i].amount[energy_i]
    protein = Rabbit_FP_Nutrient_gr[i].amount[protein_i]
    fat = Rabbit_FP_Nutrient_gr[i].amount[fat_i]
    carb = Rabbit_FP_Nutrient_gr[i].amount[carb_i]
    if fiber_i == nothing
        fiber = NaN
    elseif fiber_i > 0
        fiber = Rabbit_FP_Nutrient_gr[i].amount[fiber_i]
    end
    if sugar_i == nothing
        sugar = NaN
    elseif sugar_i > 0
        sugar = Rabbit_FP_Nutrient_gr[i].amount[sugar_i]
    end
    if starch_i == nothing
        starch = NaN
    elseif starch_i > 0
        starch = Rabbit_FP_Nutrient_gr[i].amount[starch_i]
    end
    calcium = Rabbit_FP_Nutrient_gr[i].amount[calcium_i]
    push!(Nutrient_Data, [Rabbit_FP_Nutrient_gr[i].fdc_id[1],
        water,energy,protein,fat,carb,fiber,sugar,starch,calcium])
end

# Combine portion & nutrient data ==============================================
Rabbit_Food_Data = innerjoin(Rabbit_Food_PortionWts, Nutrient_Data, on = :fdc_id)
#CSV.write(dirdata*"Rabbit Food Data.csv", Rabbit_Food_Data)
# remove ------------------------------

#   average values for duplicates ----------------------------------------------

# Convert nutrient data: per 100g to per 1 c. ==================================
water_bycup = Rabbit_Food_Data.water/100 .* Rabbit_Food_Data.gram_weight
energy_bycup = Rabbit_Food_Data.energy/100 .* Rabbit_Food_Data.gram_weight
protein_bycup = Rabbit_Food_Data.protein/100 .* Rabbit_Food_Data.gram_weight
fat_bycup = Rabbit_Food_Data.fat/100 .* Rabbit_Food_Data.gram_weight
carb_bycup = Rabbit_Food_Data.carb/100 .* Rabbit_Food_Data.gram_weight
fiber_bycup = Rabbit_Food_Data.fiber/100 .* Rabbit_Food_Data.gram_weight
sugars_bycup = Rabbit_Food_Data.sugars/100 .* Rabbit_Food_Data.gram_weight
starch_bycup = Rabbit_Food_Data.starch/100 .* Rabbit_Food_Data.gram_weight
ca_bycup = Rabbit_Food_Data.Ca/100 .* Rabbit_Food_Data.gram_weight

Rabbit_Food_Data_bycup = hcat(Rabbit_Food_Data.food, Rabbit_Food_Data.data_type,
    Rabbit_Food_Data.fdc_id,Rabbit_Food_Data.description,
    water_bycup,energy_bycup,protein_bycup,fat_bycup,carb_bycup,fiber_bycup,sugars_bycup,
    starch_bycup,ca_bycup)
Rabbit_Food_Data_bycup = DataFrame(Rabbit_Food_Data_bycup,
    [:"food",:"data_type",:"fdc_id",:"description",:"water",:"energy",:"protein",
    :"fat",:"carb",:"fiber",:"sugars",:"starch",:"Ca"])
#CSV.write(dirdata*"Rabbit Food Data by Cup.csv", Rabbit_Food_Data_bycup)

# end timer ====================================================================
#end

max_stars_repo_name: WandaVFeng/SomeBunnysDiet

max_issues_repo_name: WandaVFeng/SomeBunnysDiet

max_forks_repo_name: WandaVFeng/SomeBunnysDiet


text:
abstract type AbstractDynamicSplitter{F} end

struct DynamicBottleneckSplitter{F} <: AbstractDynamicSplitter{F}
    f::F
end

@inline _dynamic_splitter_combine(::DynamicBottleneckSplitter) = max

struct DynamicTotalSplitter{F} <: AbstractDynamicSplitter{F}
    f::F
end

@inline _dynamic_splitter_combine(::DynamicTotalSplitter) = +

function partition_stripe(A::SparseMatrixCSC{Tv, Ti}, K, method::AbstractDynamicSplitter, args...) where {Tv, Ti}
    #Reference Implementation
    @inbounds begin
        (m, n) = size(A)

        f = oracle_stripe(StepHint(), method.f, A, args...)
        g = _dynamic_splitter_combine(method)

        ptr = zeros(Ti, n + 1, K)
        cst = fill(typemax(cost_type(f)), n + 1, K)

        cst[1, 1] = f(1, 1, 1)
        ptr[1, 1] = 1
        for j′ = 2 : n + 1
            cst[j′, 1] = Step(f)(Same(1), Next(j′), Same(1))
            ptr[j′, 1] = 1
        end

        for k = 2:K
            j′₀ = k == K ? n + 1 : 1
            for j′ = j′₀ : n + 1
                cst[j′, k] = g(cst[1, k - 1], f(1, j′, k))
                ptr[j′, k] = 1
                for j = 2 : j′
                    c′ = g(cst[j, k - 1], Step(f)(Next(j), Same(j′), Same(k)))
                    if c′ <= cst[j′, k]
                        cst[j′, k] = c′
                        ptr[j′, k] = j
                    end
                end
            end
        end

        return unravel_splits(K, n, PermutedDimsArray(ptr, (2, 1)))
    end
end

function partition_stripe(A::SparseMatrixCSC{Tv, Ti}, K, method::AbstractDynamicChunker, args...) where {Tv, Ti}
    @inbounds begin
        (m, n) = size(A)

        f = oracle_stripe(StepHint(), method.f, A, args...)
        g = _dynamic_chunker_combine(method)

        ptr = zeros(Ti, K, n + 1)
        cst = fill(typemax(cost_type(f)), K, n + 1)

        for j′ = 1:n + 1
            K₁ = j′ == n + 1 ? K : K - 1
            Δc = f(1, j′)
            cst[1, j′] = Δc
            ptr[1, j′] = 1
            for k = 2:K₁
                c′ = g(cst[k - 1, 1], Δc)
                if c′ <= cst[k, j′]
                    cst[k, j′] = c′
                    ptr[k, j′] = 1
                end
            end
            for j = 2:j′
                Δc = Step(f)(Next(j), Same(j′))
                for k = 2:K₁
                    c′ = g(cst[k - 1, j], Δc)
                    if c′ <= cst[k, j′]
                        cst[k, j′] = c′
                        ptr[k, j′] = j
                    end
                end
            end
        end
        return unravel_splits(K, n, ptr)
    end
end

function unravel_splits(K, n, ptr)
    @inbounds begin
        spl = zeros(eltype(ptr), K + 1)
        spl[end] = n + 1
        for k = K:-1:1
            spl[k] = ptr[k, spl[k + 1]]
        end

        return SplitPartition(K, spl)
    end
end

struct WindowConstrainedMatrix{Tv, Ti} <: AbstractMatrix{Tv}
    z::Tv
    m::Int
    n::Int
    idx_lo::Vector{Ti}
    idx_hi::Vector{Ti}
    pos::Vector{Ti}
    val::Vector{Tv}
end

function WindowConstrainedMatrix{Tv, Ti}(z, m, n, idx_lo::Vector{Ti}, idx_hi::Vector{Ti}) where {Ti, Tv}
    @inbounds begin
        pos = undefs(Ti, n + 1)
        pos[1] = 1
        for j = 1:n
            pos[j + 1] = pos[j] + idx_hi[j] - idx_lo[j] + 1
        end
        WindowConstrainedMatrix{Tv, Ti}(z, m, n, idx_lo, idx_hi, pos)
    end
end

function WindowConstrainedMatrix{Tv, Ti}(z, m, n, idx_lo::Vector{Ti}, idx_hi::Vector{Ti}, pos::Vector{Ti}) where {Ti, Tv}
    return WindowConstrainedMatrix{Tv, Ti}(z, m, n, idx_lo, idx_hi, pos, undefs(Tv, pos[end]))
end

Base.size(A::WindowConstrainedMatrix) = (A.m, A.n)
@propagate_inbounds function Base.getindex(A::WindowConstrainedMatrix{Tv}, i, j) where {Tv}
    @boundscheck checkbounds(A, i, j)
    if A.idx_lo[j] <= i <= A.idx_hi[j]
        return A.val[A.pos[j] + i - A.idx_lo[j]]
    else
        return A.z
    end
end
@propagate_inbounds function Base.setindex!(A::WindowConstrainedMatrix{Tv}, v, i, j) where {Tv}
    @boundscheck checkbounds(A, i, j)
    if A.idx_lo[j] <= i <= A.idx_hi[j]
        return A.val[A.pos[j] + i - A.idx_lo[j]] = v
    else
        return A.z
    end
end

function column_constraints(A::SparseMatrixCSC{Tv, Ti}, K, w, w_max) where {Tv, Ti}
    @inbounds begin
        (m, n) = size(A)

        j′_lo = undefs(Ti, K)
        j′ = n + 1
        for k = K:-1:1
            j′_lo[k] = j′
            j = j′
            while j - 1 >= 1 && w(j - 1, j′, k) <= w_max
                j -= 1
            end
            j′ = j
        end

        j′_hi = undefs(Ti, K)
        j = 1
        for k = 1:K
            j′ = j
            while j′ + 1 <= n + 1 && w(j, j′ + 1, k) <= w_max
                j′ += 1
            end
            j′_hi[k] = j′
            j = j′
        end

        return (j′_lo, j′_hi)
    end
end

function part_constraints(A::SparseMatrixCSC{Tv, Ti}, K, w, w_max) where {Tv, Ti}
    @inbounds begin
        (m, n) = size(A)

        k_hi = zeros(Ti, n + 1)
        j′ = n + 1
        k_hi[end] = K
        for k = K:-1:1
            j = j′
            while j - 1 >= 1 && w(j - 1, j′, k) <= w_max
                j -= 1
                k_hi[j] = k - 1
            end
            j′ = j
        end

        k_lo = zeros(Ti, n + 1)
        j = 1
        k_lo[1] = 1
        for k = 1:K
            j′ = j
            while j′ + 1 <= n + 1 && w(j, j′ + 1, k) <= w_max
                j′ += 1
                k_lo[j′] = k
            end
            j = j′
        end

        return (k_lo, k_hi)
    end
end

function partition_stripe(A::SparseMatrixCSC{Tv, Ti}, K, method::AbstractDynamicSplitter{<:ConstrainedCost}, args...) where {Tv, Ti}
    @inbounds begin
        (m, n) = size(A)

        f = oracle_stripe(StepHint(), method.f.f, A, args...)
        w = oracle_stripe(StepHint(), method.f.w, A, args...)
        w_max = method.f.w_max
        g = _dynamic_splitter_combine(method)

        (j′_lo, j′_hi) = column_constraints(A, K, w, w_max)

        if j′_hi[K] < n + 1
            spl = ones(Ti, K + 1)
            spl[end] = n + 1
            #TODO throw(ArgumentError("infeasible"))
            return SplitPartition(K, spl)
        end

        ptr = WindowConstrainedMatrix{Ti, Ti}(zero(Ti), n + 1, K, j′_lo, j′_hi)
        cst = WindowConstrainedMatrix{cost_type(f), Ti}(typemax(cost_type(f)), n + 1, K, j′_lo, j′_hi, ptr.pos)

        for j′ = j′_lo[1] : j′_hi[1]
            cst[j′, 1] = f(1, j′, 1)
            ptr[j′, 1] = 1
        end

        for k = 2:K
            j₀ = j′_lo[k - 1]
            for j′ = j′_lo[k] : j′_hi[k]
                while w(j₀, j′, k) > w_max #TODO can this run over the end?
                    j₀ += 1
                end
                cst[j′, k] = g(cst[j₀, k - 1], f(j₀, j′, k))
                ptr[j′, k] = j₀
                for j = j₀ + 1 : min(j′, j′_hi[k - 1])
                    c′ = g(cst[j, k - 1], Step(f)(Next(j), Same(j′), Same(k)))
                    if c′ <= cst[j′, k]
                        cst[j′, k] = c′
                        ptr[j′, k] = j
                    end
                end
            end
        end

        return unravel_splits(K, n, PermutedDimsArray(ptr, (2, 1)))
    end
end

function partition_stripe(A::SparseMatrixCSC{Tv, Ti}, K, method::AbstractDynamicChunker{<:ConstrainedCost}, args...) where {Tv, Ti}
    @inbounds begin
        (m, n) = size(A)

        f = oracle_stripe(StepHint(), method.f.f, A, args...)
        w = oracle_stripe(StepHint(), method.f.w, A, args...)
        w_max = method.f.w_max
        g = _dynamic_chunker_combine(method)

        (k_lo, k_hi) = part_constraints(A, K, w, w_max)

        if k_lo[n + 1] == 0
            spl = ones(Ti, K + 1)
            spl[end] = n + 1
            #TODO throw(ArgumentError("infeasible"))
            return SplitPartition(K, spl)
        end

        ptr = WindowConstrainedMatrix{Ti, Ti}(zero(Ti), K, n + 1, k_lo, k_hi)
        cst = WindowConstrainedMatrix{cost_type(f), Ti}(typemax(cost_type(f)), K, n + 1, k_lo, k_hi, ptr.pos)

        # matrix notation...
        # i = 1:m rows, j = 1:n columns
        m, n = size(A)

        f = oracle_stripe(StepHint(), method.f.f, A)
        w = oracle_stripe(StepHint(), method.f.w, A, args...)
        w_max = method.f.w_max

        j₀ = 1
        for j′ = 1:n + 1
            while w(j₀, j′) > w_max
                j₀ += 1
            end
            @assert j₀ <= j′ #TODO infeasibility

            Δc = f(j₀, j′)
            if j₀ == 1
                cst[1, j′] = Δc
                ptr[1, j′] = 1
            end
            for k = max(k_lo[j₀] + 1, k_lo[j′]):min(k_hi[j₀] + 1, k_hi[j′])
                c′ = g(cst[k - 1, j₀], Δc)
                cst[k, j′] = c′
                ptr[k, j′] = j₀
            end
            for j = j₀ + 1:j′
                Δc = Step(f)(Next(j), Same(j′))
                for k = max(k_lo[j] + 1, k_lo[j′]):min(k_hi[j] + 1, k_hi[j′])
                    c′ = g(cst[k - 1, j], Δc)
                    if c′ <= cst[k, j′]
                        cst[k, j′] = c′
                        ptr[k, j′] = j
                    end
                end
            end
        end

        return unravel_splits(K, n, ptr)
    end
end


text:
"""
正方晶系相关
"""

"""
普通的正方格子
"""
function common_square_lattice(μ)
    disp(x, y) = -2*(cos(x) + cos(y)) + μ
    return QuadrateSystem{:COM}(
        [disp]
    )
end


text:
# This file is taken from BinaryBuilder.jl (license: MIT)
# https://github.com/JuliaPackaging/BinaryBuilder.jl

function create_or_update_pull_request(repo::GitHub.Repo, params; auth)
    try
        new_pr = GitHub.create_pull_request(repo; auth, params)
        @info("Created a new pull request: $(repo.full_name)#$(new_pr.number)")
        return new_pr
    catch ex
        # If it was already created, search for it so we can update it:
        if Registrator.CommentBot.is_pr_exists_exception(ex)
            prs, _ = GitHub.pull_requests(repo; auth=auth, params=Dict(
                "state" => "open",
                "base" => params["base"],
                "head" => string(split(repo.full_name, "/")[1], ":", params["head"]),
            ))
            existing_pr = GitHub.update_pull_request(repo, first(prs).number; auth, params)
            @info("Updated an existing pull request: $(repo.full_name)#$(existing_pr.number)")
            return existing_pr
        else
            rethrow(ex)
        end
    end
end


text:
using openQC
using Base.Test

include("utils/consts.jl")
include("utils/qusat.jl")
include("states/state.jl")
include("aqc/aqc.jl")

max_stars_repo_name: Roger-luo/openQC

max_issues_repo_name: Roger-luo/openQC

max_forks_repo_name: Roger-luo/openQC


text:
"""
    AbstractName

Abstract supertype for name wrappers.
"""
abstract type AbstractName end

Base.convert(T::Type{<:AbstractString}, name::AbstractName) = convert(T, string(name))

"""
    NoName <: AbstractName

    NoName()

NoName specifies an array is not named, and is the default `name`
value for all `AAbstractDimArray`s.
"""
struct NoName <: AbstractName end

Base.Symbol(::NoName) = Symbol("")
Base.string(::NoName) = ""

"""
    Name <: AbstractName

    Name(name::Union{Symbol,Name) => Name
    Name(name::NoName) => NoName

Name wrapper. This lets arrays keep symbol names when the array wrapper neeeds
to be `isbits, like for use on GPUs. It makes the name a property of the type.
It's not necessary to use in normal use, a symbol is probably easier.
"""
struct Name{X} <: AbstractName end
Name(name::Symbol) = Name{name}()
Name(name::NoName) = NoName()
Name(name::Name) = name

Base.Symbol(::Name{X}) where X = X
Base.string(::Name{X}) where X = string(X)

name(x::Name) = x
name(x) = name(typeof(x))
name(x::Type) = ""


text:
function LText(parent::Scene, text::String; kwargs...)
    LText(parent; text = text, kwargs...)
end

function LText(parent::Scene; bbox = nothing, kwargs...)
    attrs = merge!(Attributes(kwargs), default_attributes(LText))

    @extract attrs (text, textsize, font, color, visible, halign, valign,
        rotation, padding)

    sizeattrs = sizenode!(attrs.width, attrs.height)

    alignment = lift(tuple, halign, valign)

    autosizenode = Node{NTuple{2, Optional{Float32}}}((nothing, nothing))

    computedsize = computedsizenode!(sizeattrs, autosizenode)

    suggestedbbox = create_suggested_bboxnode(bbox)

    finalbbox = alignedbboxnode!(suggestedbbox, computedsize, alignment,
        sizeattrs, autosizenode)

    textpos = Node(Point3f0(0, 0, 0))

    t = text!(parent, text, position = textpos, textsize = textsize, font = font, color = color,
        visible = visible, align = (:center, :center), rotation = rotation, raw = true)[end]

    textbb = BBox(0, 1, 0, 1)

    onany(text, textsize, font, rotation, padding) do text, textsize, font, rotation, padding
        textbb = FRect2D(boundingbox(t))
        autowidth = width(textbb) + padding[1] + padding[2]
        autoheight = height(textbb) + padding[3] + padding[4]
        autosizenode[] = (autowidth, autoheight)
    end

    onany(finalbbox, padding) do bbox, padding

        tw = width(textbb)
        th = height(textbb)

        box = bbox.origin[1]
        boy = bbox.origin[2]

        tx = box + padding[1] + 0.5 * tw
        ty = boy + padding[3] + 0.5 * th

        textpos[] = Point3f0(tx, ty, 0)
    end

    # text has no protrusions
    protrusions = Node(RectSides(0f0, 0f0, 0f0, 0f0))

    layoutnodes = LayoutNodes{LText, GridLayout}(suggestedbbox, protrusions, computedsize, autosizenode, finalbbox, nothing)

    # trigger first update, otherwise bounds are wrong somehow
    text[] = text[]
    # trigger bbox
    suggestedbbox[] = suggestedbbox[]

    lt = LText(parent, layoutnodes, t, attrs)

    lt
end

defaultlayout(lt::LText) = ProtrusionLayout(lt)

function align_to_bbox!(lt::LText, bbox)
    lt.layoutnodes.suggestedbbox[] = bbox
end

computedsizenode(lt::LText) = lt.layoutnodes.computedsize
protrusionnode(lt::LText) = lt.layoutnodes.protrusions

function Base.getproperty(lt::LText, s::Symbol)
    if s in fieldnames(LText)
        getfield(lt, s)
    else
        lt.attributes[s]
    end
end

function Base.setproperty!(lt::LText, s::Symbol, value)
    if s in fieldnames(LText)
        setfield!(lt, s, value)
    else
        lt.attributes[s][] = value
    end
end

function Base.propertynames(lt::LText)
    [fieldnames(LText)..., keys(lt.attributes)...]
end

function Base.delete!(lt::LText)

    disconnect_layoutnodes!(lt.layoutnodes.gridcontent)
    remove_from_gridlayout!(lt.layoutnodes.gridcontent)
    empty!(lt.layoutnodes.suggestedbbox.listeners)
    empty!(lt.layoutnodes.computedbbox.listeners)
    empty!(lt.layoutnodes.computedsize.listeners)
    empty!(lt.layoutnodes.autosize.listeners)
    empty!(lt.layoutnodes.protrusions.listeners)

    # remove the plot object from the scene
    delete!(lt.parent, lt.text)
end


text:
using Statistics,QuadGK,FFTW,Distributions
using BenchmarkTools,Random

#FFT算放中提及到的特征函数变化
function M76_CharacFun(u::Float64,x0::Float64,T::Float64,r::Float64,σ::Float64,λ::Float64,μ::Float64,δ::Float64)
    ω=x0/T+r-0.5*(σ^2)-λ*(exp(μ+0.5*(δ^2))-1)
    tmp1=λ*(map(x->exp(x)-1,1im*u*μ .-0.5*(u.^2)*δ^2))
    value=map(x->exp(x*T),(1im*u*ω .-0.5*(u .^2)*σ^2 .+tmp1))

    return value
end

#M76模型下的FFT方法计算欧式看涨期权
function M76_value_call_FFT(S0::Float64,K::Float64,T::Float64,r::Float64,σ::Float64,λ::Float64,μ::Float64,δ::Float64)

    """

    # Arguments
    - `S0`: 初始标的值.
    - `K`: 行权值
    - `T`: 到期时间，年为单位
    - `r`: 无风险利率
    - `σ`: 波动率
    - `λ`: 泊松过程强度
    - `μ`: 跳跃过程中均值
    - `δ`: 跳跃过程中方差

    """

    k=log(K/S0)
    x0=log(1)
    g=2
    N=g*4096
    eps0=(g*150.0)^-1
    eta=2*pi/(N*eps0)
    b=0.5*N*eps0-k
    u=Vector(1:N)
    vo=eta*(u .-1)

    if S0>0.95*K
        α=1.5
        v=vo .-(α+1)*1im
        mod_char_fun=exp(-r*T)*M76_CharacFun(v,x0,T,r,σ,λ,μ,δ) ./(α^2+α .-(vo .^2)+1im*(2*α+1)*vo)
    else
        α=1.1
        v=(vo .-1im*α) .-1im
        mod_char_fun_1=exp(-r*T)*(1 ./(1 .+1im*(vo .-α*1im))-exp(r*T) ./(1im*(vo .-α*1im))-M76_CharacFun(v,x0,T,r,σ,λ,μ,δ) ./((vo .-α*1im) .^2-1im*(vo .-α*1im)))
        v=(vo .+α*1im) .-1im
        mod_char_fun_2=exp(-r*T)*(1 ./(1 .+1im*(vo .+α*1im))-exp(r*T) ./(1im*(vo .+α*1im))-M76_CharacFun(v,x0,T,r,σ,λ,μ,δ) ./((vo .+α*1im) .^2-1im*(vo .+α*1im)))
    end
    δ=zeros(N)
    δ[1]=1
    J=Vector(1:N)
    SimpsonW=((3 .+((-1) .^J))-δ)/3

    if S0>=0.95*K
        fft_func=eta*map(x->exp(x),1im*b*vo) .*mod_char_fun .*SimpsonW
        payoff=real(fft(fft_func))
        call_value_m=exp(-α*k)/pi*payoff
    else
        fft_func=0.5*eta*map(x->exp(x),1im*b*vo) .*(mod_char_fun_1-mod_char_fun_2) .*SimpsonW
        payoff=real(fft(fft_func))
        call_value_m=payoff/(sinh(α*k)*pi)
    end

    pos=round(Int64,(k+b)/eps0)
    call_value=call_value_m[pos]

    return call_value*S0
end

#M76模型下的蒙特卡洛方法计算欧式看涨期权
function M76_value_call_MCS(S0::Float64,K::Float64,T::Float64,r::Float64,σ::Float64,λ::Float64,μ::Float64,δ::Float64,
    M::Int,I::Int,dist::Int)

    """

    # Arguments
    - `S0`: 初始标的值.
    - `K`: 行权值
    - `T`: 到期时间，年为单位
    - `r`: 无风险利率
    - `σ`: 波动率
    - `λ`: 泊松过程强度
    - `μ`: 跳跃过程中均值
    - `δ`: 跳跃过程中方差
    - `I`: 模拟路径数
    - `M`: 时间步数
    """

    dt=T/M
    rj=λ*(exp(μ+0.5*δ^2)-1)

    S=zeros(Float64,M+1,I)
    S[1,:]=map(x->x=S0,S[1,:])

    Random.seed!(10000)

    rand1=rand(Normal(0,1),M+1,I)
    rand2=rand(Normal(0,1),M+1,I)
    rand3=rand(Poisson(λ*dt),M+1,I)

    @simd for t=2:M+1
        if dist==1
            p0=(1+(r-rj)*dt) .+σ*√(dt)*rand1[t,:]
            p1=map(x->exp(μ+δ*x)-1,rand2[t,:]) .*rand3[t,:]
            S[t,:]=S[t-1,:] .*(p0+p1)
        else
            p_0=map(x->exp((r-rj-0.5*σ^2)*dt+σ*√(dt)*x),rand1[t,:])
            p_1=map(x->exp(μ+δ*x)-1,rand2[t,:]) .*rand3[t,:]
            S[t,:]=S[t-1,:] .*(p_0+p_1)
        end
    end

    h=map(x->max(x-K,0),S[end,:])
    res=exp(-r*T)*sum(h)/I

    return res

end

S0=3225.93
K=3225.93
T=0.22
r=0.005
σ=0.113
λ=3.559
μ=-0.075
δ=0.041

M=50
I=200000

dist=2

# test=M76_value_call_FFT(S0,K,T,r,σ,λ,μ,δ)

test=M76_value_call_MCS(S0,K,T,r,σ,λ,μ,δ,M,I,dist)


text:
# All HTTP status codes, as a Dict of code => description
const STATUS_CODES = Dict([
    (100, "Continue"),
    (101, "Switching Protocols"),
    (102, "Processing"),                          # RFC 2518, obsoleted by RFC 4918
    (200, "OK"),
    (201, "Created"),
    (202, "Accepted"),
    (203, "Non-Authoritative Information"),
    (204, "No Content"),
    (205, "Reset Content"),
    (206, "Partial Content"),
    (207, "Multi-Status"),                        # RFC 4918
    (300, "Multiple Choices"),
    (301, "Moved Permanently"),
    (302, "Moved Temporarily"),
    (303, "See Other"),
    (304, "Not Modified"),
    (305, "Use Proxy"),
    (307, "Temporary Redirect"),
    (400, "Bad Request"),
    (401, "Unauthorized"),
    (402, "Payment Required"),
    (403, "Forbidden"),
    (404, "Not Found"),
    (405, "Method Not Allowed"),
    (406, "Not Acceptable"),
    (407, "Proxy Authentication Required"),
    (408, "Request Time-out"),
    (409, "Conflict"),
    (410, "Gone"),
    (411, "Length Required"),
    (412, "Precondition Failed"),
    (413, "Request Entity Too Large"),
    (414, "Request-URI Too Large"),
    (415, "Unsupported Media Type"),
    (416, "Requested Range Not Satisfiable"),
    (417, "Expectation Failed"),
    (418, "I'm a teapot"),                        # RFC 2324
    (422, "Unprocessable Entity"),                # RFC 4918
    (423, "Locked"),                              # RFC 4918
    (424, "Failed Dependency"),                   # RFC 4918
    (425, "Unordered Collection"),                # RFC 4918
    (426, "Upgrade Required"),                    # RFC 2817
    (428, "Precondition Required"),               # RFC 6585
    (429, "Too Many Requests"),                   # RFC 6585
    (431, "Request Header Fields Too Large"),     # RFC 6585
    (440, "Login Timeout"),
    (444, "nginx error: No Response"),
    (495, "nginx error: SSL Certificate Error"),
    (496, "nginx error: SSL Certificate Required"),
    (497, "nginx error: HTTP -> HTTPS"),
    (499, "nginx error or Antivirus intercepted request or ArcGIS error"),
    (500, "Internal Server Error"),
    (501, "Not Implemented"),
    (502, "Bad Gateway"),
    (503, "Service Unavailable"),
    (504, "Gateway Time-out"),
    (505, "HTTP Version Not Supported"),
    (506, "Variant Also Negotiates"),             # RFC 2295
    (507, "Insufficient Storage"),                # RFC 4918
    (509, "Bandwidth Limit Exceeded"),
    (510, "Not Extended"),                        # RFC 2774
    (511, "Network Authentication Required"),     # RFC 6585
    (520, "CloudFlare Server Error: Unknown"),
    (521, "CloudFlare Server Error: Connection Refused"),
    (522, "CloudFlare Server Error: Connection Timeout"),
    (523, "CloudFlare Server Error: Origin Server Unreachable"),
    (524, "CloudFlare Server Error: Connection Timeout"),
    (525, "CloudFlare Server Error: Connection Failed"),
    (526, "CloudFlare Server Error: Invalid SSL Ceritificate"),
    (527, "CloudFlare Server Error: Railgun Error"),
    (530, "Site Frozen")
])


text:
### Packages used
using Distributions

function writeIBW(filename,inputData)

    # Convert matrix input data to a vector with the two matrix columns after one another
    inputData = reshape(inputData, :, 1)

    A = open( filename, "w")
    write(A, ltoh(Int8(5)))
    close(A)

    A = open( filename, "w")
    ### Note that the following is just to copy the binary format of the IBW file.
    ### We are filling the bytes with garbage, the only point is that the garbage
    ### takes the exact amount of space as the header in the IBW file. This allows
    ### us to pull out the data at the end of the function.
    write(A, Int16(0))
    write(A,UInt16(0)) # Does not work

    for xLoop = 1:16
        write(A, Int32(0))
    end

    write(A, UInt32(0))
    write(A, UInt32(0))
    write(A, Int32(length(inputData)))  # npts, works but not yet calculated based on the matrix to write!

    write(A, Int16(2))   # type = single

    write(A, Int16(0))
    for aLoop = 1:6
        write(A, Char('x'))
    end
    write(A, Int16(0))

    for aLoop = 1:32
        write(A, Char('x'))
    end

    for aLoop = 1:6
        write(A, Int32(0))
    end

    write(A,Float64(0.0005)) # sFA[1]
    write(A,Float64(1.0))
    write(A,Float64(1.0))
    write(A,Float64(1.0))

    write(A,Float64(0.0))    # sFB[1]
    write(A,Float64(0.0))
    write(A,Float64(0.0))
    write(A,Float64(0.0))

    for aLoop = 1:4
        write(A, Char('x'))
    end

    for aLoop = 1:16
        write(A, Char('x'))
    end

    write(A, Int16(0))
    write(A, Int16(0))
    write(A, Float64(0.0))
    write(A, Float64(0.0))
    for aLoop = 1:26
        write(A, Int32(0))
    end
    for aLoop = 1:3
        write(A, Int16(0))
    end
    write(A, Char('x'))
    write(A, Char('x'))
    write(A, Int32(0))
    write(A, Int16(0))
    write(A, Int16(0))
    write(A, Int32(0))
    write(A, Int32(0))
    write(A, Int32(0))

    write(A, Float32.(inputData))
    # Here is where we write the real data.

    close(A)
    # Close file and end operation.
end

function generateSyntheticForceDisplacementSignal(matOfKeyPoints,sampleRate ; stdOfNoise = 0.0,sizeOfNoise =[0.0 0.0 0.0])
    # Should generate a well controlled force-displacement curve for testing.

    enrichedSignal = Array{Float64,2}(undef,0,3)
    for aLoop = 1:(size(matOfKeyPoints,1)-1)
        numelResults = length(matOfKeyPoints[aLoop,1]:1/sampleRate:matOfKeyPoints[aLoop+1,1])
        tmpMat = collect(matOfKeyPoints[aLoop,1]:1/sampleRate:matOfKeyPoints[aLoop+1,1])
        tmpMat2 = collect(range(matOfKeyPoints[aLoop,2], matOfKeyPoints[aLoop+1,2], length = numelResults))
        tmpMat3 = collect(range(matOfKeyPoints[aLoop,3], matOfKeyPoints[aLoop+1,3], length = numelResults))
        enrichedSignal = [enrichedSignal ; [tmpMat tmpMat2 tmpMat3]]

    end
    noise = sizeOfNoise.*rand(Normal(0,stdOfNoise),size(enrichedSignal))
    enrichedSignal = noise.+enrichedSignal
    return enrichedSignal[:, 3:-1:2]
end

max_stars_repo_name: abrandberg/afmIndentation

max_issues_repo_name: abrandberg/afmIndentation

max_forks_repo_name: abrandberg/afmIndentation


text:
using Documenter, MLKernels, LinearAlgebra

makedocs(
    format = :html,
    sitename="MLKernels.jl",
    authors = "Tim Thatcher",
    pages = Any[
        hide("Home" => "index.md"),
        "Interface" => "interface.md",
        "Kernels" => "kernels.md",
        "Kernel Theory" => "kernel-theory.md"
    ]
)

deploydocs(
    repo = "github.com/trthatcher/MLKernels.jl.git"
)


text:
# One step of the metropolis algorithm on a Ising's spin graph
function stepmetropolis!(grid::Array{Int, 2}; # Spin grid
                         h::Float64    = 0.0, # External field
                         temp::Float64 = 1.0) # Temperature

    # Randomly pick a position within the grid
    x   = rand(1:size(grid, 1))
    y   = rand(1:size(grid, 2))

    # Calculate the ΔE for switching the spin
    m      = nspins(grid, x, y) |> sum
    eplus  = -m - h
    ΔE     = -2eplus * grid[x,y]

    # Change spin accordingly
    if ΔE <= 0 || rand() < exp(-ΔE/temp)
        grid[x,y] *= -1
    end
end

# Several steps of the heat bath algorithm on a Ising's spin graph
function metropolis!(grid::Array{Int, 2};  # Spin grid
                     h::Float64=0.0,       # External field
                     temp::Float64=1.0,    # Temperature
                     iters::Integer=50000, # Number of iterations
                     plot::Bool=true,      # Plot flag
                     verbose::Bool=true)   # Verbose flag

    m = Float64[]
    for i in 1:iters
        stepmetropolis!(grid, h=h, temp=temp)
        push!(m, magnetization(grid))
        # Must find a better way to decide convergence
    end

    if verbose println("(T=$temp) Metropolis ended with magnetization $(m[end]) after $(length(m)) iterations") end
    if plot
        PyPlot.plot(1:length(m), m, "o", color="blue")
        PyPlot.plot(1:length(m), m, "-", color="blue")
        PyPlot.title("Metropolis for T=$temp for size $size ($(length(m)) iterations and H=$h)")
        PyPlot.xlabel("Number of Iterations")
        PyPlot.ylabel("Magnetization")
        PyPlot.ylim(0,1.1)
        PyPlot.savefig("metropolis_T=$temp.png")
        PyPlot.close()
    end

    return m
end


text:
BT = GenerativeBeliefMDP{FullyObservablePOMDP{AugScene,Int64},MergingUpdater,MergingBelief,Int64}
function POMDPs.initialstate(bmdp::BT, rng::AbstractRNG)
   s0 = initialstate(bmdp.pomdp, rng)
   b0 = initialize_belief(bmdp.updater, s0)
   return b0
end

function POMDPs.convert_s(t::Type{V}, b::MergingBelief, bmdp::BT) where V<:AbstractArray
    return convert_s(t, b, bmdp.pomdp.mdp)
end

function POMDPs.actionindex(bmdp::BT, a::Int64)
    return actionindex(bmdp.pomdp.mdp, a)
end

function Base.rand(rng::AbstractRNG, b::MergingBelief)
    return b.o
end

POMDPs.isterminal(bmdp::BT, b::MergingBelief) = isterminal(bmdp.pomdp.mdp, b.o)

POMDPs.n_actions(bmdp::BT) = n_actions(bmdp.pomdp.mdp)


text:
function binarize_proportion(g::Union{SimpleWeightedGraph,SimpleWeightedDiGraph}, prop)
    b = nonzeros(g.weights)
    N = round(Int, length(b) * prop)
    return binarize_absolute(g, sort(b)[N])
end

function compute_threshold_from_proportion(
    g::Union{SimpleWeightedGraph,SimpleWeightedDiGraph},
    prop::AbstractVector
)
    b = nonzeros(g.weights)
    N = length(b)
    sb = sort(b)
    return [sb[round(Int, N * p)] for p in prop]
end

function binarize_absolute(g, thresh)
    out = SimpleGraph(nv(g))
    for e in edges(g)
        if weight(e) > thresh
            add_edge!(out, src(e), dst(e))
        end
    end
    return out
end

each_threshold(op, g, thresh) = [op(binarize_absolute(g, t)) for t in thresh]
each_proportion(op, g, prop) = each_threshold(op, g, compute_threshold_from_proportion(g, prop))


text:
#=
We check the deposition step (i.e. computing \rho) works as expected.
=#
using Test

@testset "PIC deposition : uniform    " begin
    using LinearAlgebra
    using VlasovSolvers
    # using Plots

    include("PIC.jl")

    function samples_verif1(nsamples, kx::Float64)
        x0 = Array{Float64}(undef, nsamples)
        seq = SobolSeq(1)
        for i = 1:nsamples
            r = Sobol.next!(seq)
            x0[i] = 2π/kx * r[1]
        end
        v0 = rand(Uniform(-10, 10), nsamples)
        wei = 2π/kx/nsamples .* ones(Float64, nsamples)
        return x0, v0, wei
    end

    kx = 0.5
    L = 2π / kx
    np = 100000 # nb de particules
    nx = 32   # nb de mailles pour le calcul de l energie
    mesh1 = OneDGrid(CPU(), nx, 0, L)

    (x0, v0, w0) = samples_verif1(np, kx)
    p = Particles(x0, v0, w0, np)

    ρ, ρ_tot = compute_rho(p, mesh1)

    # arrx = mesh1.start:mesh1.step:(mesh1.stop - mesh1.step)
    # plot(arrx, ρ)

    @test maximum(abs.(ρ .- 1.)) < 10^-3
    @test ρ_tot ≈ 1.0
end

@testset "PIC deposition : 1 + 0.99cos" begin
    using LinearAlgebra, Roots
    using VlasovSolvers
    # using Plots

    include("PIC.jl")

    function samples_verif2(nsamples)
        x0 = Array{Float64}(undef, nsamples)
        seq = SobolSeq(1)
        for i = 1:nsamples
            r = Sobol.next!(seq)
            resx = find_zero(x-> x + 0.99 / kx * sin(kx*x) - r[1]*2π/kx, 0.0)
            x0[i] = resx
        end
        v0 = rand(Uniform(-10, 10), nsamples)
        wei = 2π/kx/nsamples .* ones(Float64, nsamples)
        return x0, v0, wei
    end

    kx = 0.5
    L = 2π / kx
    np = 100000 # nb de particules
    nx = 128   # nb de mailles pour le calcul de l energie
    mesh1 = OneDGrid(CPU(), nx, 0, L)

    (x0, v0, w0) = samples_verif2(np)
    p = Particles(x0, v0, w0, np)

    ρ, ρ_tot = compute_rho(p, mesh1)

    arrx = mesh1.start:mesh1.step:(mesh1.stop - mesh1.step)

    @test maximum(abs.(ρ - (1 .+ 0.99cos.(kx.*arrx)))) < 5.10^-3
    @test ρ_tot ≈ 1.0
end


text:
# ---
# title: 112. Path Sum
# id: problem112
# author: zhwang
# date: 2022-02-21
# difficulty: Easy
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/path-sum/description/>
# hidden: true
# ---
#
# Given a binary tree and a sum, determine if the tree has a root-to-leaf path
# such that adding up all the values along the path equals the given sum.
#
# **Note:**  A leaf is a node with no children.
#
# **Example:**
#
# Given the below binary tree and `sum = 22`,
#
#
#
#           **5**
#          **/** \
#         **4**   8
#        **/**   / \
#       **11**  13  4
#      /  **\**      \
#     7    **2**      1
#
#
# return true, as there exist a root-to-leaf path `5->4->11->2` which sum is 22.
#
#
## @lc code=start
using LeetCode

has_path_sum(::Nothing, ::Int) = false
function has_path_sum(root::TreeNode, target_sum::Int)::Bool
    root.val == target_sum && isnothing(root.left) && isnothing(root.right) && return true
    any(has_path_sum.([root.left, root.right], target_sum - root.val))
end

## @lc code=end
## @lc test=start
@testset "112.path-sum.jl" begin
    @test has_path_sum(TreeNode{Int}([5,4,8,11,nothing,13,4,7,2,nothing,nothing,nothing,1]), 22) == true
    @test has_path_sum(TreeNode{Int}([1,2]), 1) == false
    @test has_path_sum(nothing, 0) == false
end
## @lc test=end

## @lc info=start
# link: [solution to 112](https://leetcode-cn.com/problems/path-sum/solution/112-lu-jing-zong-he-by-rex-rs-qppi/)
## @lc info=end


text:
"""
    reorder(A::AbstractDimArray, order::Pair) => AbstractDimArray
    reorder(A::Dimension, order::Order) => AbstractDimArray

Reorder every dims index/array to `order`, or reorder index for
the the given dimension(s) to the `Order` they wrap.

`order` can be an [`Order`](@ref), or `Dimeension => Order` pairs.

If no axis reversal is required the same objects will be returned, without allocation.
"""
function reorder end

reorder(x, p::Pair, ps::Vararg{<:Pair}) = reorder(x, (p, ps...))
reorder(x, ps::Tuple{Vararg{<:Pair}}) = reorder(x, Dimensions.pairdims(ps...))
# Reorder specific dims.
reorder(x, dimwrappers::Tuple) = _reorder(x, dimwrappers)
# Reorder all dims.
reorder(x, ot::Order) = reorder(x, typeof(ot))
reorder(x, ot::Type{<:Order}) = _reorder(x, map(d -> rebuild(d, ot), dims(x)))
reorder(dim::Dimension, ot::Type{<:Order}) =
    ot <: basetypeof(order(dim)) ? dim : reverse(dim)

# Recursive reordering. x may be reversed here
_reorder(x, orderdims::DimTuple) = _reorder(reorder(x, orderdims[1]), tail(orderdims))
_reorder(x, orderdims::Tuple{}) = x

reorder(x, orderdim::Dimension) = _reorder(val(orderdim), x, dims(x, orderdim))

_reorder(neworder::Order, x, dim::DimOrDimType) = _reorder(basetypeof(neworder), x, dim)
# Reverse the dimension index
_reorder(::Type{O}, x, dim::DimOrDimType) where O<:Ordered =
    order(dim) isa O ? x : reverse(x; dims=dim)
_reorder(ot::Type{Unordered}, x, dim::DimOrDimType) = x

"""
    modify(f, A::AbstractDimArray) => AbstractDimArray
    modify(f, s::AbstractDimStack) => AbstractDimStack
    modify(f, dim::Dimension) => Dimension
    modify(f, x, lookupdim::Dimension) => typeof(x)

Modify the parent data, rebuilding the object wrapper without
change. `f` must return a `AbstractArray` of the same size as the original.

This method is mostly useful as a way of swapping the parent array type of
an object.

## Example

If we have a previously-defined `DimArray`, we can copy it to an Nvidia GPU with:

```julia
A = DimArray(rand(100, 100), (X, Y))
modify(CuArray, A)
```

This also works for all the data layers in a `DimStack`.
"""
function modify end
modify(f, s::AbstractDimStack) = map(a -> modify(f, a), s)
function modify(f, A::AbstractDimArray)
    newdata = f(parent(A))
    size(newdata) == size(A) || error("$f returns an array with a different size")
    rebuild(A, newdata)
end
modify(f, x, dim::DimOrDimType) = set(x, modify(f, dims(x, dim)))
modify(f, dim::Dimension) = rebuild(dim, modify(f, val(dim)))
function modify(f, lookup::LookupArray)
    newindex = modify(f, parent(lookup))
    rebuild(lookup; data=newindex)
end
function modify(f, index::AbstractArray)
    newindex = f(index)
    size(newindex) == size(index) || error("$f returns a vector with a different size")
    newindex
end

"""
    dimwise(f, A::AbstractDimArray{T,N}, B::AbstractDimArray{T2,M}) => AbstractDimArray{T3,N}

Dimension-wise application of function `f` to `A` and `B`.

## Arguments

- `a`: `AbstractDimArray` to broacast from, along dimensions not in `b`.
- `b`: `AbstractDimArray` to broadcast from all dimensions. Dimensions must be a subset of a.

This is like broadcasting over every slice of `A` if it is
sliced by the dimensions of `B`.
"""
function dimwise(f, A::AbstractDimArray, B::AbstractDimArray)
    dimwise!(f, similar(A, promote_type(eltype(A), eltype(B))), A, B)
end

"""
    dimwise!(f, dest::AbstractDimArray{T1,N}, A::AbstractDimArray{T2,N}, B::AbstractDimArray) => dest

Dimension-wise application of function `f`.

## Arguments

- `dest`: `AbstractDimArray` to update
- `A`: `AbstractDimArray` to broacast from, along dimensions not in `b`.
- `B`: `AbstractDimArray` to broadcast from all dimensions. Dimensions must be a subset of a.

This is like broadcasting over every slice of `A` if it is
sliced by the dimensions of `B`, and storing the value in `dest`.
"""
function dimwise!(
    f, dest::AbstractDimArray{T,N}, a::AbstractDimArray{TA,N}, b::AbstractDimArray{TB,NB}
) where {T,TA,TB,N,NB}
    N >= NB || error("B-array cannot have more dimensions than A array")
    comparedims(dest, a)
    common = commondims(a, dims(b))
    od = otherdims(a, common)
    # Lazily permute B dims to match the order in A, if required
    if !dimsmatch(common, dims(b))
        b = PermutedDimsArray(b, common)
    end
    # Broadcast over b for each combination of dimensional indices D
    if length(od) == 0
        dest .= f.(a, b)
    else
        map(DimIndices(od)) do D
            dest[D...] .= f.(a[D...], b)
        end
    end
    return dest
end

# Get a tuple of unique keys for DimArrays. If they have the same
# name we call them layerI.
function uniquekeys(das::Tuple{AbstractDimArray,Vararg{<:AbstractDimArray}})
    uniquekeys(Symbol.(map(name, das)))
end
function uniquekeys(keys::Tuple{Symbol,Vararg{<:Symbol}})
    ids = ntuple(x -> x, length(keys))
    map(keys, ids) do k, id
        count(k1 -> k == k1, keys) > 1 ? Symbol(:layer, id) : k
    end
end

_astuple(t::Tuple) = t
_astuple(x) = (x,)


text:
using Test
using SizedStrings

@testset "SizedStrings.jl" begin

    include("SizedString.jl")

    include("PaddedString.jl")

end


text:
"""generate a new sample from a probability density using slice sampling

Parameters
----------
init_x : array
logprob : callable, `lprob = logprob(x, *logprob_args)`
    A functions which returns the log probability at a given
    location

Returns
-------
new_x : float
    the sampled position
new_llh : float
    the log likelihood at the new position
Notes
-----
http://en.wikipedia.org/wiki/Slice_sampling
"""
function slicesample(init_x::Vector{Float64},
                     logprob::Function;
                     sigma = 1.0,
                     step_out=true,
                     max_steps_out=1000,
                     compwise=true,
                     numdir=2,
                     doubling_step = true,
                     verbose = true,
                     upper_bound = Inf,
                     lower_bound = -Inf)

    ## define a univariate directional slice sampler
    function direction_slice(direction, init_x)

        function dir_logprob(z)
            return logprob(direction*z + init_x)
        end

        function acceptable(z, llh_s, L, U)
            while (U-L) > 1.1*sigma
                middle = 0.5*(L+U)
                splits = ((middle > 0) & (z >= middle)) |
                         ((middle <= 0) & (z < middle))
                if z < middle
                    U = middle
                else
                    L = middle
                end
                # Probably these could be cached from the stepping out.
                if splits & (llh_s >= dir_logprob(U)) & (llh_s >= dir_logprob(L))
                    return false
                end
            end
            return true
        end

        # if we have upper/lower bounds on parameters, compute them in z space
        #if compwise:
        #    dir_upper_bound = upper_bound[direction==1.]
        #    dir_lower_bound = lower_bound[direction==1.]
        #else:
        #    dir_upper_bound = np.min(np.sign(direction)*(upper_bound - init_x)/direction)
        #    dir_lower_bound = np.max(np.sign(direction)*(lower_bound - init_x)/direction)

        # compute initial interval bounds
        upper = sigma * rand()
        lower = upper - sigma

        # sample uniformly under the probability at z = 0
        llh_s = log(rand()) + dir_logprob(0.0)

        # perform the stepping out or doubling procedure to compute interval I
        l_steps_out = 0
        u_steps_out = 0
        if step_out
            if doubling_step
                while ((dir_logprob(lower) > llh_s) | (dir_logprob(upper) > llh_s)) &
                      ((l_steps_out + u_steps_out) < max_steps_out)
                    if rand() < 0.5
                        l_steps_out += 1
                        lower       -= (upper-lower)
                    else
                        u_steps_out += 1
                        upper       += (upper-lower)
                    end
                end
            else
                while (dir_logprob(lower) > llh_s) & (l_steps_out < max_steps_out)
                    l_steps_out += 1
                    lower       -= sigma
                end
                while (dir_logprob(upper) > llh_s) & (u_steps_out < max_steps_out)
                    u_steps_out += 1
                    upper       += sigma
                end
            end
        end

        # uniformly sample - perform shrinkage (with accept check)
        start_upper = upper
        start_lower = lower
        steps_in = 0
        new_z    = 0.
        new_llh  = -Inf
        while true
            steps_in += 1
            if steps_in % 100 == 0
                println("shrinking, steps", steps_in)
            end

            new_z     = (upper - lower)*rand() + lower
            new_llh   = dir_logprob(new_z)
            if isnan(new_llh)
                println("new_z: ", new_z)
                println("new_th: ", direction*new_z + init_x)
                println("new_llh: ", new_llh)
                println("llh_s: ", llh_s)
                println("init_x", init_x)
                println("ll_init_x", logprob(init_x))
                throw("Slice sampler got a NaN")
            end

            # accept/or shrink
            if (new_llh > llh_s) & acceptable(new_z, llh_s, start_lower, start_upper)
                break
            elseif new_z < 0
                lower = new_z
            elseif new_z > 0
                upper = new_z
            else
                throw("Slice sampler shrank to zero!")
            end
        end
        #println(" new z, llh post: ", new_z, new_llh)

        #if verbose:
        #    print "Steps Out:", l_steps_out, u_steps_out, " Steps In:", steps_in
        return new_z*direction + init_x, new_llh
    end

    ## expand upper and lower bound to be same dimension as samples
    #if type(upper_bound) == float or isinstance(upper_bound, np.number):
    #    upper_bound = np.tile(upper_bound, init_x.shape)
    #if type(lower_bound) == float or isinstance(lower_bound, np.number):
    #    lower_bound = np.tile(lower_bound, init_x.shape)
    ## make sure we're starting in a legal spot
    #assert np.all(init_x < upper_bound), "init_x (%s) >= ub (%s)" % \
    #    (np.array_str(init_x), np.array_str(upper_bound))
    #assert np.all(init_x > lower_bound), "init_x (%s) <= lb (%s)" % \
    #    (np.array_str(init_x), np.array_str(lower_bound))

    # do either component-wise slice sampling, or random direction
    dims = length(init_x)
    new_x, new_llh = nothing, nothing
    if compwise
        ordering = shuffle(1:dims)
        new_x = copy(init_x)
        for d in ordering
            direction    = zeros(dims)
            direction[d] = 1.0
            new_x, new_llh = direction_slice(direction, new_x)
        end
    else
        new_x = copy(init_x)
        for d in 1:numdir
            direction    = randn(dims)
            direction    = direction / sqrt(sum(direction.^2))
            new_x, new_llh = direction_slice(direction, new_x)
        end
    end

    return new_x, new_llh
end

"""
Helper that draws N samples using slice sampling
"""
function slicesample_chain(lnpdf, th, N; print_skip=50)
    NUM_COLOR_COMPONENTS = length(th)
    samps = zeros(N, NUM_COLOR_COMPONENTS)
    lls   = zeros(N)
    tic()
    @printf "  iter : \t loglike \n"
    for n in 1:N
        th, ll = slicesample(th, lnpdf;
                             doubling_step=true,
                             compwise=true)
        samps[n,:] = th
        lls[n]     = ll
        if mod(n, print_skip) == 0
            @printf "   %d   : \t %2.4f \n" n ll
        end
    end
    elapsed = toc()
    @printf "%2.3f ms per sample (20k samples in %2.3f seconds) \n" 1000*(elapsed/N) elapsed
    return samps, lls
end


text:
import DocStringExtensions
import DocStringExtensions: SIGNATURES, TYPEDSIGNATURES, DOCSTRING,
    FIELDS, TYPEDFIELDS

DocStringExtensions.@template (FUNCTIONS, METHODS, MACROS) =
    """
    $(TYPEDSIGNATURES)

    $(DOCSTRING)
    """

DocStringExtensions.@template TYPES =
    """
    Fields:
    $(TYPEDFIELDS)

    ---

    $(DOCSTRING)
    """


text:
@testset "WeightedCovMatrix fit!" begin
    s, m, c = (map(x -> sum(x .* w), eachcol(x2)),
               map(x -> mean(x, weights(w)), eachcol(x2)),
               cov(x2, aweights(w), corrected = true))
    ma, ca = map(x -> mean(x, weights(w)), eachcol(x2)), cov(x2, aweights(w), corrected = true)
    mf, cf = map(x -> mean(x, weights(w)), eachcol(x2)), cov(x2, fweights(w), corrected = true)
    mp, cp = map(x -> mean(x, weights(w)), eachcol(x2)), cov(x2, pweights(w), corrected = true)

    o = WeightedCovMatrix()
    for i in 1:l
        fit!(o, x2[i,:], w[i])
    end
    o
    sfor, mfor, cfor = sum(o), mean(o), cov(o)

    o_32 = WeightedCovMatrix(Float32)
    for i in 1:l
        fit!(o_32, x2[i,:], w[i])
    end
    o_32
    sfor_32, mfor_32, cfor_32 = sum(o_32), mean(o_32), cov(o_32)

    sval, mval, cval = fit!(WeightedCovMatrix(), x2, w) |> x -> (sum(x), mean(x), cov(x))
    szip, mzip, czip = fit!(WeightedCovMatrix(), zip(eachrow(x2), w)) |> x -> (sum(x), mean(x), cov(x))

    mvala, cvala = fit!(WeightedCovMatrix(), x2, w) |>
        x -> (mean(x), cov(x, corrected = true, weight_type = :analytic))
    mvalf, cvalf = fit!(WeightedCovMatrix(), x2, w) |>
        x -> (mean(x), cov(x, corrected = true, weight_type = :frequency))
    @test_throws ArgumentError fit!(WeightedCovMatrix(), x2, w) |>
        x -> (mean(x), cov(x, corrected = true, weight_type = :something))

    @test c ≈ czip
    @test c ≈ cfor
    @test c ≈ cfor_32
    @test c ≈ cval
    @test m ≈ mzip
    @test m ≈ mfor
    @test m ≈ mfor_32
    @test m ≈ mval
    @test ca ≈ cvala
    @test cf ≈ cvalf
    @test ma ≈ mvala
    @test mf ≈ mvalf

    @test s ≈ sfor
    @test s ≈ sfor_32
    @test s ≈ sval
    @test s ≈ szip

    # After implementing :probability, these should pass/not throw any more:
    @test_throws ErrorException mvalp, vvalp = fit!(WeightedCovMatrix(), x2, w) |>
        x -> (mean(x), cov(x, corrected = true, weight_type = :probability))
    @test_broken vp ≈ vvalp
    @test_broken mp ≈ mvalp

    @test eltype(o) == Float64
    @test eltype(o_32) == Float32

    oold = copy(o)
    fit!(o, [missing, 1, 2, 3, 4], 1)
    @test o == oold
    fit!(o, [1, 2, 3, 4, 5], missing)
    @test o == oold

    # issue #29
    fit!(o, eachrow([missing 1 2 3 4;
                     missing 2 3 4 5]),
         [1.0, 1.0])
    @test o == oold
end

@testset "WeightedCovMatrix merge!" begin
    c = cov(x2, aweights(w), corrected = true)

    wc = fit!(WeightedCovMatrix(), x2, w)
    oc = map(eachrow(x2), w) do xi, wi
        fit!(WeightedCovMatrix(), xi, wi)
    end;
    rc = reduce(merge!, deepcopy(oc))
    rc2 = merge!(
        fit!(WeightedCovMatrix(), x2[1:end ÷ 2, :],       w[1:end ÷ 2]),
        fit!(WeightedCovMatrix(), x2[end ÷ 2 + 1:end, :], w[end ÷ 2 + 1:end]))

    rc_32 = reduce(merge!, deepcopy(oc), init = WeightedCovMatrix(Float32))
    rc2_32 = merge!(
        fit!(WeightedCovMatrix(Float32), x2[1:end ÷ 2, :],       w[1:end ÷ 2]),
        fit!(WeightedCovMatrix(),        x2[end ÷ 2 + 1:end, :], w[end ÷ 2 + 1:end]))

    @test rc.b ≈ wc.b
    @test rc.C ≈ wc.C
    @test rc.W ≈ wc.W
    @test rc.W2 ≈ wc.W2

    @test rc2.b ≈ wc.b
    @test rc2.C ≈ wc.C
    @test rc2.W ≈ wc.W
    @test rc2.W2 ≈ wc.W2

    @test cov(rc) ≈ c
    @test cov(rc2) ≈ c

    @test rc_32.b ≈ wc.b
    @test rc_32.C ≈ wc.C
    @test rc_32.W ≈ wc.W
    @test rc_32.W2 ≈ wc.W2

    @test rc2_32.b ≈ wc.b
    @test rc2_32.C ≈ wc.C
    @test rc2_32.W ≈ wc.W
    @test rc2_32.W2 ≈ wc.W2

    @test cov(rc_32) ≈ c
    @test cov(rc2_32) ≈ c

    @test eltype(rc) == Float64
    @test eltype(rc2) == Float64
    @test eltype(rc_32) == Float32
    @test eltype(rc2_32) == Float32
end

@testset "WeightedCovMatrix copy" begin
    o = fit!(WeightedCovMatrix(), x2, w)

    o2 = copy(o)

    @test Ref(o.C) != Ref(o2.C)
    @test Ref(o.A) != Ref(o2.A)
    @test Ref(o.b) != Ref(o2.b)
    @test Ref(o.W) != Ref(o2.W)
    @test Ref(o.W2) != Ref(o2.W2)
    @test Ref(o.n) != Ref(o2.n)

    @test o.C == o2.C
    @test o.A == o2.A
    @test o.b == o2.b
    @test o.W == o2.W
    @test o.W2 == o2.W2
    @test o.n == o2.n
end

@testset "WeightedCovMatrix constructor" begin
    @test WeightedCovMatrix{Float64}() == WeightedCovMatrix()
    @test WeightedCovMatrix{Float32}() == WeightedCovMatrix(Float32)
    @test WeightedCovMatrix() == WeightedCovMatrix(zeros(Float64, 0, 0),
                                                   zeros(Float64, 0, 0),
                                                   zeros(Float64, 0),
                                                   0.0, 0.0, 0)
    @test convert(WeightedCovMatrix{Float32}, WeightedCovMatrix()) ==
        WeightedCovMatrix{Float32}()
end