text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineScaleSetVMProfile <: SwaggerModel
    osProfile::Any # spec type: Union{ Nothing, VirtualMachineScaleSetOSProfile } # spec name: osProfile
    storageProfile::Any # spec type: Union{ Nothing, VirtualMachineScaleSetStorageProfile } # spec name: storageProfile
    networkProfile::Any # spec type: Union{ Nothing, VirtualMachineScaleSetNetworkProfile } # spec name: networkProfile
    diagnosticsProfile::Any # spec type: Union{ Nothing, DiagnosticsProfile } # spec name: diagnosticsProfile
    extensionProfile::Any # spec type: Union{ Nothing, VirtualMachineScaleSetExtensionProfile } # spec name: extensionProfile
    licenseType::Any # spec type: Union{ Nothing, String } # spec name: licenseType

    function VirtualMachineScaleSetVMProfile(;osProfile=nothing, storageProfile=nothing, networkProfile=nothing, diagnosticsProfile=nothing, extensionProfile=nothing, licenseType=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetVMProfile, Symbol("osProfile"), osProfile)
        setfield!(o, Symbol("osProfile"), osProfile)
        validate_property(VirtualMachineScaleSetVMProfile, Symbol("storageProfile"), storageProfile)
        setfield!(o, Symbol("storageProfile"), storageProfile)
        validate_property(VirtualMachineScaleSetVMProfile, Symbol("networkProfile"), networkProfile)
        setfield!(o, Symbol("networkProfile"), networkProfile)
        validate_property(VirtualMachineScaleSetVMProfile, Symbol("diagnosticsProfile"), diagnosticsProfile)
        setfield!(o, Symbol("diagnosticsProfile"), diagnosticsProfile)
        validate_property(VirtualMachineScaleSetVMProfile, Symbol("extensionProfile"), extensionProfile)
        setfield!(o, Symbol("extensionProfile"), extensionProfile)
        validate_property(VirtualMachineScaleSetVMProfile, Symbol("licenseType"), licenseType)
        setfield!(o, Symbol("licenseType"), licenseType)
        o
    end
end # type VirtualMachineScaleSetVMProfile

const _property_map_VirtualMachineScaleSetVMProfile = Dict{Symbol,Symbol}(Symbol("osProfile")=>Symbol("osProfile"), Symbol("storageProfile")=>Symbol("storageProfile"), Symbol("networkProfile")=>Symbol("networkProfile"), Symbol("diagnosticsProfile")=>Symbol("diagnosticsProfile"), Symbol("extensionProfile")=>Symbol("extensionProfile"), Symbol("licenseType")=>Symbol("licenseType"))
const _property_types_VirtualMachineScaleSetVMProfile = Dict{Symbol,String}(Symbol("osProfile")=>"VirtualMachineScaleSetOSProfile", Symbol("storageProfile")=>"VirtualMachineScaleSetStorageProfile", Symbol("networkProfile")=>"VirtualMachineScaleSetNetworkProfile", Symbol("diagnosticsProfile")=>"DiagnosticsProfile", Symbol("extensionProfile")=>"VirtualMachineScaleSetExtensionProfile", Symbol("licenseType")=>"String")
Base.propertynames(::Type{ VirtualMachineScaleSetVMProfile }) = collect(keys(_property_map_VirtualMachineScaleSetVMProfile))
Swagger.property_type(::Type{ VirtualMachineScaleSetVMProfile }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineScaleSetVMProfile[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetVMProfile }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetVMProfile[property_name]

function check_required(o::VirtualMachineScaleSetVMProfile)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetVMProfile }, name::Symbol, val)
end


text:
if false
    Threads.nthreads()
end

function blockranges(nblocks, total_len)
    rem = total_len % nblocks
    main_len = div(total_len, nblocks)

    starts=Int[1]
    ends=Int[]
    for ii in 1:nblocks
        len = main_len
        if rem>0
            len+=1
            rem-=1
        end
        push!(ends, starts[end]+len-1)
        push!(starts, ends[end] + 1)
    end
    @assert ends[end] == total_len
    starts[1:end-1], ends
end

# import Hwloc
# topology = Hwloc.topology_load()
# l1cache = first(filter(t->t.type_==:Cache && t.attr.depth==1, topology)).attr

# l1cache = filter(t->t.type_==:Cache && t.attr.depth==1, topology)
# println("L1 cache information: $l2cache")
using BenchmarkTools, FastGroupBy, Base.Threads
import FastGroupBy.load_bits
function paradissort2!(svec::Vector{String}, rev = false, counter_type::Type{CT} = UInt) where CT <: Unsigned
    bits = load_bits.(UInt, svec)
    hist = zeros(CT, 1 << 16, 4, nthreads())
    @threads for j = 1:length(bits)
        for i = 0:3
            @inbounds hist[1+Int((bits[j] >> (i << 4)) & 0xffff), i+1, threadid()] += 1
        end
    end
    @threads for j in 1:4
        for i = 2:nthreads()
           @inbounds hist[:, j, 1] .+= hist[:, j, i]
        end
    end
    hist[:,:,1]
end

# this is provided for timing control
# DO NOT DELETE
function paradissort2a!(svec::Vector{String}, rev = false)
    bits = load_bits.(UInt, svec)
    hist = zeros(UInt, 65536, 4)
    for x1 in bits
        for i in 0:3
            @inbounds hist[1+Int((x1 >> (i << 4)) & 0xffff), i+1] += 1
        end
    end
    hist
end

function paradissort2a8!(svec::Vector{String}, rev = false)
    bits = load_bits.(UInt, svec)
    hist = zeros(UInt, 1<<16, 4)
    hist_hot = zeros(UInt8, 1<<16)
    @inbounds for i in 0:3
        for x1 in bits
            idx = 1+ Int((x1 >> (i << 4)) & 0xffff)
            hist_hot[idx] += UInt8(1)
            if hist_hot[idx] == 0
                hist[idx, i+1] += (1<<8)
            end
        end
        for idx in 1:(1<<16)
            hist[idx,i+1] += hist_hot[idx]
        end
        fill!(hist_hot, zero(UInt8))
    end
    hist
end

function paradissort2a8t!(svec::Vector{String}, rev = false)
    bits = load_bits.(UInt, svec)
    hist = zeros(UInt64, 1<<16, 4)
    hist_hot = zeros(UInt8, 1<<16, 4)
    @threads for i in 0:3
        @inbounds for x1 in bits
            idx = 1+ Int((x1 >> (i << 4)) & 0xffff)
            hist_hot[idx,i+1] += UInt8(1)
            if hist_hot[idx,i+1] == 0
                hist[idx, i+1] += (1<<8)
            end
        end
        @inbounds for idx in 1:(1<<16)
            hist[idx,i+1] += hist_hot[idx, i+1]
        end
    end
    hist
end

# svec = rand("id".*dec.(1:1_000_000,10), 100_000_000);
svec = rand([randstring(8) for i =1:1_000_000], 100_000_000);
using BenchmarkTools
t1 = @belapsed a = paradissort2!($svec, false, UInt16);
t2 = @belapsed b = paradissort2a!($svec);
t3 = @belapsed c = paradissort2a8!($svec);
t4 = @belapsed d = paradissort2a8t!($svec);
(a .== b) |> all
(a .== c) |> all
(a .== d) |> all
using Plots
bar(["naive threaded", "naive", "count8", "count8t"], [t1,t2,t3, t4])

@time radixsort!(svec)

ss = 1
ee = 25_000_000
x = svec[ss:ee];
function abc(x)
    @time FastGroupBy.radixsort!(x)
end
@time abc(x);

function def(x)
    @time FastGroupBy.radixsort!(@view(svec[ss:ee]))
end
@time def(x);

@time xx = paradissort2!(svec);

ext: jl


text:
module Parsers

export GML

include("gml/GML.jl")
include("dot/DOT.jl")

end


text:
module PalmerPenguins

using CSV: CSV
using DataDeps: DataDeps

const DEPNAME = "PalmerPenguins"
const LINK = "https://cdn.jsdelivr.net/gh/allisonhorst/palmerpenguins@433439c8b013eff3d36c847bb7a27fa0d7e353d8/inst/extdata/"
const DATA = ["penguins.csv", "penguins_raw.csv"]
const INFO = """
             Dataset: The Palmer penguins dataset
             Authors: Allison Horst, Alison Hill, Kristen Gorman
             Website: https://allisonhorst.github.io/palmerpenguins/index.html

             The Palmer penguins dataset is a dataset for data exploration & visualization, as an
             alternative to the Iris dataset.

             The dataset contains data for 344 penguins. There are 3 different species of penguins in
             this dataset, collected from 3 islands in the Palmer Archipelago, Antarctica.

             Data were collected and made available by
             [Dr. Kristen Gorman](https://www.uaf.edu/cfos/people/faculty/detail/kristen-gorman.php)
             and the [Palmer Station, Antarctica LTER](https://pal.lternet.edu/), a member of the
             [Long Term Ecological Research Network](https://lternet.edu/).

             Data are available by
             [CC-0](https://creativecommons.org/share-your-work/public-domain/cc0/) license in
             accordance with the [Palmer Station LTER Data Policy](http://pal.lternet.edu/data/policies)
             and the
             [LTER Data Access Policy for Type I data](https://lternet.edu/data-access-policy/).

             Please include this citation if you plan to use this database:
             [Gorman et al., 2014]
                 Gorman KB, Williams TD, Fraser WR (2014) Ecological Sexual Dimorphism and Environmental
                 Variability within a Community of Antarctic Penguins (Genus Pygoscelis). PLoS ONE 9(3):
                 e90081. https://doi.org/10.1371/journal.pone.0090081
             """

"""
    load([; raw = false])

Load the Palmer penguins dataset.

If `raw` is `true`, then the raw data is returned. Otherwise the simplified version of
the data is loaded.

# References

[Gorman et al., 2014]
    Gorman KB, Williams TD, Fraser WR (2014) Ecological Sexual Dimorphism and Environmental
    Variability within a Community of Antarctic Penguins (Genus Pygoscelis). PLoS ONE 9(3):
    e90081. https://doi.org/10.1371/journal.pone.0090081
"""
function load(; raw::Bool=false)
    file = DataDeps.@datadep_str DEPNAME * "/" * DATA[1 + Int(raw)]
    return CSV.File(
        file;
        missingstring="NA",
        truestrings=["Yes"],
        falsestrings=["No"],
        dateformat="y-m-d",
    )
end

function __init__()
    return DataDeps.register(
        DataDeps.DataDep(
            DEPNAME,
            INFO,
            LINK .* DATA,
            "839b058be09b164f7dfa0d030a959cb3c5426dea281e8a34d34c557e398bd01f",
        ),
    )
end

end


text:
"""
    RectangularCluster

Fills the landscape with rectangles containing a random value. The size of each
rectangle/patch is between `minimum` and `maximum` (the two can be equal for a
fixed size rectangle).
"""
struct RectangularCluster <: NeutralLandscapeMaker
    minimum::Integer
    maximum::Integer
    function RectangularCluster(x::T, y::T) where {T <: Integer}
        @assert 0 < x <= y
        new(x, y)
    end
end

RectangularCluster() = RectangularCluster(2,4)

function _landscape!(mat, alg::RectangularCluster)
    mat .= -1.0
    while minimum(mat) == -1.0
        width, height = rand(alg.minimum:alg.maximum, 2)
        row = rand(1:(size(mat,1)-(width-1)))
        col = rand(1:(size(mat,2)-(height-1)))
        mat[row:(row+(width-1)) , col:(col+(height-1))] .= rand()
    end
    return mat
end


text:
using Plots, FreqTables, CategoricalArrays

function sunflowerplot(x::AbstractVector{T}, y::AbstractVector{T};
        initialize=true, r=0, color=:black, rotate=false, title="", xlabel="", ylabel="") where T<:Real
    function petals(x, y, n, r, color, rotate)
        δ = rotate ? 180rand(1)[1] : 0
        for θ in 90+δ:(360/n):450+δ
            plot!([x, x + r*cosd(θ)], [y, y + r*sind(θ)], color=color, linewidth=2)
        end
    end
    gr(label="", grid=false, tick_direction=:out,
        markerstrokewidth=0, aspect_ratio=1)
    r == 0 && (r = min(diff(vcat(extrema(x)...)), diff(vcat(extrema(y)...)))[1]/40)
    tbl = freqtable(x, y)
    xval, yval = names(tbl)
    nx, ny = size(tbl)
    initialize && (p = plot())
    p = plot!(title=title, xlabel=xlabel, ylabel=ylabel)
    for i in 1:nx
        for j in 1:ny
            n = tbl[i, j]
            if n == 1
                scatter!([xval[i]], [yval[j]], color=color, markersize=3)
            elseif n != 0
                scatter!([xval[i]], [yval[j]], color=color, markersize=2)
                 n > 1 && petals(xval[i], yval[j], n, r, color, rotate)
            end
        end
    end
    r
    p
end

max_stars_repo_name: r-de-r/stats

max_issues_repo_name: r-de-r/stats

max_forks_repo_name: r-de-r/stats


text:
L = 20
labels = rand(Bool, L);
scores = add_noise.(labels, 0.6)

curve = DET(labels, scores)
B = BinaryMetrics(labels, scores)
curve = DET(B)
println(curve)

OP = op(B)
@test op(curve) == OP
@test fpr(curve) == fpr(B,OP)
@test tpr(curve) == tpr(B,OP)


text:
"""
    get_PFFRG_timers() :: Nothing

Test current performance of FRG solver.
"""
function get_PFFRG_timers() :: Nothing

    println()
    println("Testing performance of FRG solver ...")
    println()

    # time lattice
    println("Testing lattice performance ...")
    println()
    get_lattice_timers()
    println()
    println()

    # time frequencies
    println("Testing frequency performance ...")
    println()
    get_frequency_timers()
    println()
    println()

    # time action
    println("Testing action performance ...")
    println()
    get_action_timers()
    println()
    println()

    # time flow equations
    println("Testing flow performance ...")
    get_flow_timers()
    println()
    println()

    # time observables
    println("Testing observable performance ...")
    println()
    get_observable_timers()
    println()
    println()

    println("Done.")

    return nothing
end


text:
"""
    vgg_block(ifilters, ofilters, depth, batchnorm)

A VGG block of convolution layers
([reference](https://arxiv.org/abs/1409.1556v6)).

# Arguments

  - `ifilters`: number of input feature maps
  - `ofilters`: number of output feature maps
  - `depth`: number of convolution/convolution + batch norm layers
  - `batchnorm`: set to `true` to include batch normalization after each convolution
"""
function vgg_block(ifilters, ofilters, depth, batchnorm)
    k = (3, 3)
    p = (1, 1)
    layers = []
    for _ in 1:depth
        if batchnorm
            append!(layers, conv_bn(k, ifilters, ofilters; pad = p, bias = false))
        else
            push!(layers, Conv(k, ifilters => ofilters, relu; pad = p))
        end
        ifilters = ofilters
    end
    ifilters = ofilters
    return layers
end

"""
    vgg_convolutional_layers(config, batchnorm, inchannels)

Create VGG convolution layers
([reference](https://arxiv.org/abs/1409.1556v6)).

# Arguments

  - `config`: vector of tuples `(output_channels, num_convolutions)`
    for each block (see [`Metalhead.vgg_block`](#))
  - `batchnorm`: set to `true` to include batch normalization after each convolution
  - `inchannels`: number of input channels
"""
function vgg_convolutional_layers(config, batchnorm, inchannels)
    layers = []
    ifilters = inchannels
    for c in config
        append!(layers, vgg_block(ifilters, c..., batchnorm))
        push!(layers, MaxPool((2, 2); stride = 2))
        ifilters, _ = c
    end
    return layers
end

"""
    vgg_classifier_layers(imsize, nclasses, fcsize, dropout)

Create VGG classifier (fully connected) layers
([reference](https://arxiv.org/abs/1409.1556v6)).

# Arguments

  - `imsize`: tuple `(width, height, channels)` indicating the size after
    the convolution layers (see [`Metalhead.vgg_convolutional_layers`](#))
  - `nclasses`: number of output classes
  - `fcsize`: input and output size of the intermediate fully connected layer
  - `dropout`: the dropout level between each fully connected layer
"""
function vgg_classifier_layers(imsize, nclasses, fcsize, dropout)
    return Chain(MLUtils.flatten,
                 Dense(Int(prod(imsize)), fcsize, relu),
                 Dropout(dropout),
                 Dense(fcsize, fcsize, relu),
                 Dropout(dropout),
                 Dense(fcsize, nclasses))
end

"""
    vgg(imsize; config, inchannels, batchnorm = false, nclasses, fcsize, dropout)

Create a VGG model
([reference](https://arxiv.org/abs/1409.1556v6)).

# Arguments

  - `imsize`: input image width and height as a tuple
  - `config`: the configuration for the convolution layers
    (see [`Metalhead.vgg_convolutional_layers`](#))
  - `inchannels`: number of input channels
  - `batchnorm`: set to `true` to use batch normalization after each convolution
  - `nclasses`: number of output classes
  - `fcsize`: intermediate fully connected layer size
    (see [`Metalhead.vgg_classifier_layers`](#))
  - `dropout`: dropout level between fully connected layers
"""
function vgg(imsize; config, inchannels, batchnorm = false, nclasses, fcsize, dropout)
    conv = vgg_convolutional_layers(config, batchnorm, inchannels)
    imsize = outputsize(conv, (imsize..., inchannels); padbatch = true)[1:3]
    class = vgg_classifier_layers(imsize, nclasses, fcsize, dropout)
    return Chain(Chain(conv), class)
end

const vgg_conv_config = Dict(:A => [(64, 1), (128, 1), (256, 2), (512, 2), (512, 2)],
                             :B => [(64, 2), (128, 2), (256, 2), (512, 2), (512, 2)],
                             :D => [(64, 2), (128, 2), (256, 3), (512, 3), (512, 3)],
                             :E => [(64, 2), (128, 2), (256, 4), (512, 4), (512, 4)])

const vgg_config = Dict(11 => :A,
                        13 => :B,
                        16 => :D,
                        19 => :E)

struct VGG
    layers::Any
end

"""
    VGG(imsize::Dims{2}; config, inchannels, batchnorm = false, nclasses, fcsize, dropout)

Construct a VGG model with the specified input image size. Typically, the image size is `(224, 224)`.

## Keyword Arguments:

  - `config` : VGG convolutional block configuration. It is defined as a vector of tuples `(output_channels, num_convolutions)` for each block
  - `inchannels`::Integer : number of input channels
  - `batchnorm`::Bool : set to `true` to use batch normalization after each convolution
  - `nclasses`::Integer : number of output classes
  - `fcsize`: intermediate fully connected layer size
    (see [`Metalhead.vgg_classifier_layers`](#))
  - `dropout`: dropout level between fully connected layers
"""
function VGG(imsize::Dims{2};
             config, inchannels, batchnorm = false, nclasses, fcsize, dropout)
    layers = vgg(imsize; config = config,
                 inchannels = inchannels,
                 batchnorm = batchnorm,
                 nclasses = nclasses,
                 fcsize = fcsize,
                 dropout = dropout)

    return VGG(layers)
end

@functor VGG

(m::VGG)(x) = m.layers(x)

backbone(m::VGG) = m.layers[1]
classifier(m::VGG) = m.layers[2]

"""
    VGG(depth::Integer = 16; pretrain = false, batchnorm = false)

Create a VGG style model with specified `depth`. Available values include (11, 13, 16, 19).
([reference](https://arxiv.org/abs/1409.1556v6)).
See also [`VGG`](#).

!!! warning

    `VGG` does not currently support pretrained weights.

# Arguments

  - `pretrain`: set to `true` to load pre-trained model weights for ImageNet
"""
function VGG(depth::Integer = 16; pretrain = false, batchnorm = false, nclasses = 1000)
    @assert depth in keys(vgg_config) "depth must be from one in $(sort(collect(keys(vgg_config))))"
    model = VGG((224, 224); config = vgg_conv_config[vgg_config[depth]],
                inchannels = 3,
                batchnorm = batchnorm,
                nclasses = nclasses,
                fcsize = 4096,
                dropout = 0.5)
    if pretrain && !batchnorm
        loadpretrain!(model, string("VGG", depth))
    elseif pretrain
        loadpretrain!(model, "VGG$(depth)-BN)")
    end
    return model
end

# deprecations
@deprecate VGG11(; kw...) VGG(11; kw...)
@deprecate VGG13(; kw...) VGG(13; kw...)
@deprecate VGG16(; kw...) VGG(16; kw...)
@deprecate VGG19(; kw...) VGG(19; kw...)


text:
const water_fraction_to_M = 1.0m^3*m^-3 * 1kg*L^-1 / 18.0g*mol^-1

using Base: tail

abstract type AbstractLayerInterpolator end

struct LinearLayerInterpolator <: AbstractLayerInterpolator
    layer::Int
    frac::Float64
end

struct LinearLayerInterpolators{H} <: AbstractLayerInterpolator
    range::LinearLayerInterpolator
    increment::LinearLayerInterpolator
    height::H
end
LinearLayerInterpolators(env, height) =
    LinearLayerInterpolators(range_interpolator(env, height), increment_interpolator(env, height), height)

# Use recursion so the compiler does this with tuples. Huge performance improvement.
@inline layer_sizes(z, I::Tuple, NL) = (lsize(z, I[1], I[2]), layer_sizes(I, NL)...)
@inline layer_sizes(I::Tuple, NL) = (lsize(I[1], I[2], I[3]), layer_sizes(tail(I), NL)...)
@inline layer_sizes(I::Tuple{X,X}, NL) where X = (lsize(I[1], I[2], NL),)

@inline lsize(a, b, c) = (c + b)  / 2 - (b + a) / 2

@inline layer_bounds(inc::Tuple{X,Vararg}, next) where X = (bound(inc[1], inc[2]), layer_bounds(tail(inc), next)...)
@inline layer_bounds(inc::Tuple{X}, next) where X = (bound(inc[1], next),)

@inline bound(a::Number, b::Number) = (a + b) / 2

@inline max_height(env) = sum(layer_sizes(env))
@inline layer_props(env) = layer_sizes(env) ./ max_height(env)

" Calculate current interpolation layers and fraction from NicheMapR data"
@inline increment_interpolator(env, height) = begin
    increments = get_increments(env)
    for (i, height_upper) in enumerate(increments)
        if height_upper > height
            height_lower = increments[i - 1]
            frac = (height - height_lower) / (height_upper - height_lower)
            return LinearLayerInterpolator(i - 1, frac)
        end
    end
    # Otherwise it's taller/deeper than we have data, so use the largest we have.
    LinearLayerInterpolator(lastindex(increments) - 1, 1.0)
end

@inline range_interpolator(env, height) = begin
    lower, upper = get_range(env)
    h = min(max(height, lower), upper)
    frac = (h - lower)/(upper - lower)
    LinearLayerInterpolator(1, frac)
end

" Interpolate between two layers of environmental data. "
@inline interp_layer(layers::Matrix, t, interp) =
    lin_interp(layers, t, interp.layer) * (oneunit(interp.frac) - interp.frac) +
    lin_interp(layers, t, interp.layer + 1) * interp.frac

" Linear interpolation "
@inline lin_interp(vector::Vector, t) = begin
    t_int = floor(Int64, t)
    frac = t - t_int
    vector[t_int] * (oneunit(frac) - frac) + vector[t_int + 1] * frac
end
@inline lin_interp(matrix::Matrix, t, l) = begin
    t_int = floor(Int64, t)
    frac = t - t_int
    matrix[t_int, l] * (oneunit(frac) - frac) + matrix[t_int + 1, l] * frac
end
@inline lin_interp(vector::Vector, t::Int) = vector[t]
@inline lin_interp(matrix::Matrix, t::Int, l) = matrix[t, l]

@inline weightedmean(env, layers, t, inpterp::LinearLayerInterpolators) =
    weightedmean(env, layers, t, inpterp.height)
weightedmean(env, layers, t, height::Number) = begin
    wmean = zero(layers[1][1])*m
    lbounds = layer_bounds(env)
    lsizes = layer_sizes(env)
    h = min(height, max_height(env))
    # sum all layers < height, after size adjustment
    for l = 1:length(lbounds)
        if h > lbounds[l]
            # add the whole layer
            wmean += lin_interp(layers, t, l) * lsizes[l]
        else
            # add the last fractional layer and break the loop
            frac = 1 - ((lbounds[l] - h) / lsizes[l])
            wmean += lin_interp(layers, t, l) * lsizes[l] * frac
            break
        end
    end
    wmean / h
end

@inline layermax(layers, t, interp) = begin
    t_int = floor(Int, t)
    val = layers[t_int, 1]
    for l = 2:interp.layer
        val = max(val, layers[t_int, l])
    end
    val
end


text:
RecipesBase.@recipe function f(var::ScalarVariable{N}) where N
    marker_z --> unwrapdata(var)
    markersize --> 2.5
    markerstrokewidth --> 0
    aspect_ratio --> 1

    RecipesBase.@series begin
        seriestype := scatter
        x = getdomain(var)[1]
        y = getdomain(var)[2]
        if N > 2
            z = getdomain(var)[3]

            x, y, z
        else
            x, y
        end
    end
end

RecipesBase.@recipe function f(::Type{T}, f::T) where T <: ScalarField
    unwrapdata(f)
end


text:
export DiffCorrectAxisymmetric

#-------#
# TYPES #
#-------#

abstract type Abstract_AxisymmetricCorrector <: Abstract_DifferentialCorrector end

struct DiffCorrectAxisymmetric{D<:DiffEqBase.DEAlgorithm,F<:Abstract_ReferenceFrame} <: Abstract_AxisymmetricCorrector
    u0_free :: SVector{2,Integer}  # Indices of the free variables in the initial condition
    u1_free :: SVector{3,Integer}  # Indices of the free variables in the final condition
    alg :: D
    frame :: F

    # Search perturbations for axisymmetric case [Russell 2006]
    # NOTE: Only valid if the reference sol is already near-axisymmetric
    DiffCorrectAxisymmetric(u0_free, u1_free, alg, frame) = new{typeof(alg),typeof(frame)}(u0_free, u1_free, alg, frame)
    DiffCorrectAxisymmetric() = DiffCorrectAxisymmetric((5,6), (2,3,4), DEFAULT_ALG, SynodicFrame(true))
end

#-----------#
# FUNCTIONS #
#-----------#

abstract type AxisymmetricException <: Exception end
struct DidNotCross <: AxisymmetricException end
struct Crashed <: AxisymmetricException end
struct BadTspan <: AxisymmetricException end

function orbit_xcrossing(u, t, integrator)
    state = get(integrator.opts.userdata, :real_state, nothing)
    isnothing(state) && error("orbit_xcrossing expects a `real_state` to be provided in userdata.")
    new_state = remake(state; u0=u, tspan=(t, t))
    synodic_state = convert_to_frame(new_state, default_synodic_reference_frame(state.model))

    # XXX: This cannot return the Dual, or the downstream Duals (i.e. sol.u) get very broken when this triggers!
    return @inbounds ForwardDiff.value.(synodic_state.prob.u0[2])
end

function terminate_after_N_crossings!(integrator)
    t = integrator.t
    tspan = integrator.sol.prob.tspan
    if (t - tspan[begin]) / (tspan[end] - tspan[begin]) >= 0.01
        N = get!(integrator.opts.userdata, :crossings, 1)
        integrator.opts.userdata[:crossings] = N - 1
        if integrator.opts.userdata[:crossings] == 0
            terminate!(integrator, :NCrossings)
        end
    end
end

crossed(sol::Trajectory) = sol.retcode == :NCrossings

function corrector_callback(::Abstract_AxisymmetricCorrector, system::EphemerisNBP)
    # TODO: Play with interp_points, interp_points=0 halves runtime/memory, but might cause issues due to oscillation around y-axis
    # TODO: idxs=[2] makes things very slow, even though it seems it should speed it up.
    return VectorContinuousCallback(
        (integrator, idx) -> begin
            if idx == 3 # Crossing!
                terminate_after_N_crossings!(integrator)
            else
                terminate!(integrator, :Crashed)
            end
        end, 3;
        interp_points=10, save_positions=(false, false)) do du, u, t, integrator
            diam1 = maximum(bodvrd(String(primary_body(system)), "RADII"))  # km
            diam2 = maximum(bodvrd(String(secondary_body(system)), "RADII"))  # km
            du[1] = check_distance(u, t, system, primary_body(system), diam1)
            du[2] = check_distance(u, t, system, secondary_body(system), diam2)
            du[3] = orbit_xcrossing(u, t, integrator)

            # XXX: I think that some of the underlying SPICE functions are slightly discontinuous,
            # maybe there are differences in number of bits being used between here and in the
            # SPICE library. Basically, this causes issues in event callback detection, because it
            # gets stuck in a value close, but not quite equal to zero (e.g. on a value like -2.0e-14),
            # despite incrementing the timestep -- see https://github.com/SciML/DiffEqBase.jl/issues/646.
            # This fix is very hacky, but I have no idea what might be a better option for now.
            @. du[isapprox(du, 0; atol=1e-4)] = 0.
        end
end

function corrector_callback(::Abstract_AxisymmetricCorrector, ::Abstract_DynamicalModel)
    # TODO: Play with interp_points, interp_points=0 halves runtime/memory, but might cause issues due to oscillation around y-axis
    # TODO: idxs=[2] makes things very slow, even though it seems it should speed it up.
    return ContinuousCallback(orbit_xcrossing, terminate_after_N_crossings!; interp_points=10, save_positions=(false, false))
end

function corrector_solve(corrector::Abstract_AxisymmetricCorrector, state::State; strict=false, verbose=false, crossings=1, kwargs...)
    callback = strict ? corrector_callback(corrector, state.model) : nothing
    userdata = Dict{Symbol,Any}(:crossings => crossings)  # XXX: Dict type specified for easy merging

    sol = sensitivity_trace(AD, state, corrector.frame, corrector.alg; trace_time=true, abstol=1e-12, reltol=1e-12, callback, userdata, kwargs...)
    sol_tspan = ForwardDiff.value.((sol.t[begin], sol.t[end]))

    # Error-checking
    if strict
        try
            crashed(sol) && throw(Crashed())
            !crossed(sol) && throw(DidNotCross())
        catch ex
            if verbose && isa(ex, AxisymmetricException)
                display(plot(solve(remake(state, tspan=sol_tspan)), corrector.frame; title="$(ex)"))
            end
            rethrow(ex)
        end
    end

    if verbose
        display(plot(sol, corrector.frame))
    end

    return convert_to_frame(sol, corrector.frame), sol_tspan
end

function (corrector::Abstract_AxisymmetricCorrector)(state::State, F, J, x; kwargs...)
    tspan = (state.tspan[begin], state.tspan[end] + x[end])
    if (tspan[end] - tspan[begin]) < 0 || (tspan[end] - tspan[begin]) > 1e8
        # Don't like this tspan, it's either negative time direction, or too long
        throw(BadTspan())
    end

    # Recreate the state with the parameters provided in x
    correction_state = convert_to_frame(state, corrector.frame)
    u0 = deepcopy(correction_state.u0)
    u0[corrector.u0_free] .= x[1:length(corrector.u0_free)]
    correction_state = remake(correction_state; u0, tspan)

    # Compute sensitivities!
    sol, sol_tspan = corrector_solve(corrector, correction_state; strict=false, kwargs...)

    if !isnothing(F)
        # Calculate how far away we are from the desirable final state (i.e. the residual).
        end_u0 = ForwardDiff.value.(sol.u[end])
        F .= [end_u0[corrector.u1_free]...]
    end

    if !isnothing(J)
        # As per [Russell.2006, Eq.(10) to (13)], we need to add the time variation term,
        # which requires the values of [ẏ, ż, ẍ]. Here, these are provided in the STM
        # on the "end" column (since we're calling corrector_solve with trace_time=true).
        STM = extract_STMs([sol.u[end]])[1]
        J .= STM[[corrector.u1_free...], [corrector.u0_free..., end]]
    end
end


text:
# This file calculates the third application in Y.-S. Chan, A. C. Fannjiang, and G. H. Paulino,
# Integral equations with hypersingular kernels -- theory and applications to fracture mechanics,
# Int. J. Eng. Sci., 41:683--720, 2003.

using ApproxFun, SingularIntegralEquations

x = Fun(identity)
w = 1/sqrt(1-x^2)
d = domain(x)
d2 = d^2
B = [dirichlet(d);neumann(d)]
H2 = Hilbert(d,2)
H4 = Hilbert(d,4)
ϵ = 0.2
L = -6ϵ^2*H4[w] + H2[w]
f = -Fun(one)

uSIE = [B;L]\[zeros(4);f]
@time uSIE = [B;L]\[zeros(4);f]
println("The length of uSIE is: ",length(uSIE))
println("The extrema of uSIE are: ",extrema(uSIE))

test0 = uSIE/(1-x^2)^2
test = Fun(x->test0(x),Ultraspherical(1,d))
temp = 0.0
[temp+= i*test.coefficients[i] for i=1:length(test)];
println("The normalized generalized Stress Intensity Factors are: ",3ϵ*temp)


text:
# ManagePic.jl

using FileIO, Images, ImageShow

"""
Convert a file to color matrices.

# Arguments
- `filename`::String`: name of file to load.
"""
function generateMatricesRGB(filename)
    img = load(filename)
    r = (Float64.(red.(img)))
    g = (Float64.(green.(img)))
    b = (Float64.(blue.(img)))
    return [r,g,b]
end # RGB

"""
Convert a loaded image to color matrices.

# Arguments
- `img`::Array{RGB{Normed{UInt8,8}},2}`: raw image.
"""
function generateMatricesRGB(img)
    r = (Float64.(red.(img)))
    g = (Float64.(green.(img)))
    b = (Float64.(blue.(img)))
    return [r,g,b]
end # RGB

"""
Convert color matrices o RGB matrice.

# Arguments
- `r`::Array{Float, 2}` full-size matrix of red color.
- `g`::Array{Float, 2}` full-size matrix of green color.
- `b`::Array{Float, 2}` full-size matrix of blue color.
"""
function matriceRGB(r, g, b)
    return ([RGB(r[i,j], g[i, j], b[i,j]) for i in 1:size(r)[1], j in 1:size(r)[2] ])
end # matriceRGB

"""
Save a r g b matrices as a picture.

# Arguments
- `newFilename`::String`: name of file.
- `r`::Array{Float, 2}` full-size matrix of red color.
- `g`::Array{Float, 2}` full-size matrix of green color.
- `b`::Array{Float, 2}` full-size matrix of blue color.
"""
function savePicture(newFilename, r, g, b)
    aa = matriceRGB(r, g, b)
    save(newFilename, aa)
end # save

"""
Save RGB matrice as a picture.

# Arguments
- `newFilename`::String`: name of file.
- `RGB`::Array{Array{Float64,2},1}`: RGB matrices list.
"""
function savePictureRGB(newFilename, RGB)
    savePicture(newFilename, RGB[1], RGB[2], RGB[3])
end # save2

"""
Convert RGB to HSL.
Return image as HSL matrices list.

# Arguments
- `rgb`::Array{Array{Float64,2},1}`: r, g, b matrices list.
"""
function rgb2hsl(rgb)
    r, g, b = rgb

    arraySize = size(r)
    h = zeros(Float64, arraySize[1], arraySize[2])
    s = zeros(Float64, arraySize[1], arraySize[2])
    l = zeros(Float64, arraySize[1], arraySize[2])

    # L:Lightness - value of white light
    l = lightness(rgb)

    Cmax = max.(r, g, b)
    Cmin = min.(r, g, b)
    Δ = Cmax .- Cmin

    for i in 1:arraySize[1], j in 1:arraySize[2]

        # H:Hue - value of color
        if Δ[i,j] == 0
            h[i,j] = 0
        elseif Cmax[i,j] == r[i,j]
            h[i,j] = 60 * (((g[i,j] - b[i,j]) ./ Δ[i,j]) .% 6)
        elseif Cmax[i,j] == g[i,j]
            h[i,j] = 60 * (((b[i,j] - r[i,j]) ./ Δ[i,j]) .+ 2)
        else
            h[i,j] = 60 * (((r[i,j] - g[i,j]) ./ Δ[i,j]) .+ 4)
        end

        # S:Saturation -  color saturation
        if Δ[i,j] == 0
            s[i,j] = 0
        else
            s[i,j] = Δ[i,j] ./ (1 - abs(2 * l[i,j] - 1))
        end
    end
    return [h,s,l]
end

"""
Convert HSL to RGB.
Return image as RGB matrices list.

#Arguments
- `hsl`::Array{Array{Float64,2},1}`: HSL matrices list.
"""
function hsl2rgb(hsl)
    h, s, l = hsl
    arraySize = size(h)
    r = zeros(Float64, arraySize[1], arraySize[2])
    g = zeros(Float64, arraySize[1], arraySize[2])
    b = zeros(Float64, arraySize[1], arraySize[2])

    for i in 1:arraySize[1], j in 1:arraySize[2]
        C = (1 - abs(2 * l[i,j] - 1)) .* s[i,j]
        X = C .* (1 - abs((h[i,j] ./ 60) % 2 - 1))
        m = l[i,j] - C ./ 2

        if h[i,j] >= 0 && h[i,j] < 60
            r[i,j] = C + m
            g[i,j] = X + m
            b[i,j] = m
        elseif h[i,j] >= 60 && h[i,j] < 120
            r[i,j] = X + m
            g[i,j] = C + m
            b[i,j] = m
        elseif h[i,j] >= 120 && h[i,j] < 180
            r[i,j] = m
            g[i,j] = m + C
            b[i,j] = m + X
        elseif h[i,j] >= 180 && h[i,j] < 240
            r[i,j] = m
            g[i,j] = X + m
            b[i,j] = C + m
        elseif h[i,j] >= 240 && h[i,j] < 300
            r[i,j] = X + m
            g[i,j] = m
            b[i,j] = C + m
        else
            r[i,j] = C + m
            g[i,j] = m
            b[i,j] = X + m
        end
    end
    return [r,g,b]
end

"""
Count lightness.

Arduments:
- `rgb`::Array{Array{Float64,2},1}`: r, g, b matrices list.
"""
function lightness(rgb)
    r, g, b = rgb
    l = (max.(r, g, b) .+ min.(r, g, b)) ./ 2
    return l
end

max_stars_repo_name: o-Mateo-o/ImageEditorApp

max_issues_repo_name: o-Mateo-o/ImageEditorApp

max_forks_repo_name: o-Mateo-o/ImageEditorApp


text:
## Reactive signal update performance: TimeArray slices
using TimeSeries
using Reactive

nma = 50
ohlc_BA = TimeSeries.readtimearray("../../data/OHLC_BA.csv")

s_ohlc = Reactive.Input(ohlc_BA[1:nma])
s_close = lift(s -> values(s["Close"])[end], s_ohlc, typ=Float64)

vsig = [s_close.value]
function vsigupd!(vsig::Vector{Float64}, nma::Int64, nmax::Int64)
    for i in (nma + 1):nmax
        push!(s_ohlc, ohlc_BA[i-nma:i])
        push!(vsig, s_close.value)
    end
end

println("OHLC updates TimeArray slices version:")
println("reactive signal updates timing")
expr_run = parse("@time vsigupd!(vsig, nma, length(ohlc_BA))")
for i = 1:5
    eval(expr_run)
end

### results
# elapsed time: 1.062034687 seconds (367636432 bytes allocated, 33.78% gc time)

Profile.clear()  # in case we have any previous profiling data
@profile vsigupd!(vsig, nma, length(ohlc_BA))
Profile.print()


text:
module PumasDocs

end


text:
using Test
using NeXLCore
using NeXLDatabase
using SQLite
using Dates
using Statistics
using Pkg.Artifacts

function tifffilename(i)
    res = "$i.tif"
    return repeat('0', max(0, 9 - length(res))) * res
end

# Download the necessary data using the Artifact mechanism from Google Drive
zip = artifact"shooter0"
dbname = ":memory:"
# dbname = "C:\\Users\\nritchie\\Desktop\\test.db" # tempname()
(dbname!=":memory:") && (rm(dbname, force=true)==nothing)
db = openNeXLDatabase(dbname)
ld1 = write(db, DBPerson, "Butterbean Q. Grenfeld, III", "bqg@gmail.com")
ld2 = write(db, DBPerson, "Jeanne I. Bottle", "jib@gmail.com")
ld3 = write(db, DBPerson, "Harvey Sun Beetle", "hsb@gmail.com")

l1 = write(db, DBLaboratory, "OpenLab", read(db, DBPerson, ld1))
l2 = write(db, DBLaboratory, "ClosedLab", read(db, DBPerson, ld2))

write(db, DBMember, l1, 1)
write(db, DBMember, l2, 2)
write(db, DBMember, l2, ld3)

i1 = write(db, DBInstrument, l1, "JEOL", "JXA-7001", "H103")
i2 = write(db, DBInstrument, l1, "TESCAN", "MIRA3", "H119")
i3 = write(db, DBInstrument, l2, "Hitachi", "SU-2000", "Gr19")
i4 = write(db, DBInstrument, l2, "Cameca", "SX-5", "Gr22")

d1 = write(db, DBDetector, i1, "Bruker", "Esprit 6|30", "30 mm² SDD", 125.0, 106, -495.0, 5.0, 4, 21, 55, 94)
d2 = write(db, DBDetector, i2, "Oxford", "XMAX 100", "100 mm² SDD", 130.0, 20, 0.0, 5.0, 4, 18, 50, 90)
d3 = write(db, DBDetector, i3, "EDAX", "Octane", "30 mm² SDD", 128.0, 20, 0.0, 5.0, 4, 21, 55, 94)
d4 = write(db, DBDetector, i4, "Thermo", "UltraDry", "100 mm² SDD", 126.0, 20, 0.0, 10.0, 4, 21, 55, 94)
d5 = write(db, DBDetector, i2, "Pulsetor", "Torrent", "4 × 30 mm² SDD", 132.0, 10, 0.0, 10.0, 4, 21, 55, 94)

mats = NeXLCore.compositionlibrary()
for (name, mat) in mats
    write(db, mat)
end

@testset "NeXLDatabase" begin
    # db = openNeXLDatabase(dbname)
    det3 = convert(BasicEDS, read(db, DBDetector, d3), 4096)
    @test isapprox(resolution(energy(n"Mn K-L3"), det3), 128.0, atol = 0.001)
    @test energy(1, det3) == 0.0
    @test length(NeXLSpectrum.isvisible(characteristic(n"Ca", ltransitions), det3)) == 0
    @test read(db,DBPerson,"dale.newbury@nist.gov").name == "Dale Newbury"

    k240 = read(db, Material, "K240")
    lab = read(db, DBLaboratory, l1)
    s1 = write(db, DBSample, l1, "SPI REP", "Rare Earth Phosphates")
    hsb = read(db,DBPerson,"hsb@gmail.com")
    testproj = read(db, DBProject, write(db, DBProject, "Tests", "Test projects", hsb))
    @test testproj.name == "Tests"
    @test testproj.description == "Test projects"
    @test testproj.createdBy.name == hsb.name
    proj = read(db, DBProject, write(db, DBProject, "REP", "Rare Earth Phosphates", hsb, testproj))
    @test proj.name == "REP"
    @test proj.description == "Rare Earth Phosphates"
    @test proj.parent.name == testproj.name
    for mat in ("CeP5O14", "LaP5O14", "NdP5O14", "PrP5O14", "SmP5O14", "YP5O14")
        midx = write(db, parse(Material, mat))
        sidx = write(
            db,
            DBSpectrum,
            d5,
            20.0e3,
            midx,
            ld2,
            s1,
            now(),
            mat * " standard",
            joinpath(joinpath(@__DIR__, "spectra"), mat * " std.msa"),
            "EMSA",
        )
        write(db, NeXLDatabase.DBProjectSpectrum, proj.pkey, sidx)
    end
    @testset "Base DB" begin
        @test isapprox(value(k240[n"Ba"]), 0.2687, atol = 0.000001)
        @test find(db, DBPerson, "jib@gmail.com") == ld2
        jib = read(db, DBPerson, find(db, DBPerson, "jib@gmail.com"))
        @test jib.name == "Jeanne I. Bottle"
        @test jib.email == "jib@gmail.com"
        @test jib.pkey == ld2
        @test find(db, DBPerson, "unknown@gmail.com") == -1
        @test find(db, DBPerson, "bqg@gmail.com") == ld1
        @test find(db, DBPerson, "nicholas.ritchie@nist.gov") == 1

        lab = read(db, DBLaboratory, l1)
        @test lab.pkey == l1
        @test lab.name == "OpenLab"
        @test lab.contact.pkey == find(db, DBPerson, "bqg@gmail.com")
    end

    @testset "Amy's GSR test" begin
        SQLite.transaction(db) do
            @test find(db, DBProject, 0, "Tests") == testproj.pkey
            write(db, DBPerson, "Amy", "amy@thelab.com")

            proj = write(db, DBProject, "GSR", "From Amy's GSR", read(db, DBPerson, "amy@thelab.com"), testproj)
            s2 = write(db, DBSample, l1, "Shooter #1 - Zero time", "A sample collected by the Boston Police")
            sidx = write(
                db,
                DBSpectrum,
                d5,
                20.0e3,
                1,
                ld1,
                s2,
                now(),
                "03272.tif",
                joinpath(joinpath(@__DIR__, "spectra"), "03272.tif"),
                "ASPEX",
            )
            write(db, NeXLDatabase.DBProjectSpectrum, proj, sidx)

            specs = read(db, DBProject, DBSpectrum, proj)
            @test length(specs) == 1
            sps = asa.(Spectrum, specs)
            @test sps[1][480] == 110.0

            for i in 1:100
                fn = tifffilename(i)
                sidx = write(
                    db,
                    DBSpectrum,
                    d5,
                    20.0e3,
                    1,
                    ld1,
                    s2,
                    now(),
                    fn,
                    joinpath(zip,fn),
                    "ASPEX",
                )
                write(db, NeXLDatabase.DBProjectSpectrum, proj, sidx)
            end
        end
        spec = asa(Spectrum, read(db, DBSpectrum, 95))
        @test spec[369] == 22.0
        @test max(spec[1:500]...) == 733.0
    end

    @testset "ADM-6005a" begin
        SQLite.transaction(db) do
            path = joinpath(@__DIR__, "ADM-6005a")
            project = write(db, DBProject, "ADM-6005a Test", "Description of ADM-6005a Test", read(db, DBPerson, "hsb@gmail.com"), testproj)
            det, person, e0, comp = 1, 2, 20.0e3, find(db, Material, "ADM6005a")
            campaign =
                write(db, NeXLDatabase.DBCampaign, project, det, [n"O", n"Al", n"Si", n"Ca", n"Ti", n"Zn", n"Ge"])
            sample = write(db, DBSample, 1, "ADM-6005a block", "ADM-6005a prepared by Eric Windsor")
            dt = DateTime(Date(2019, 7, 12), Time(9, 30, 0))
            for i in 1:15
                fn = "$path\\ADM-6005a_$(i).msa"
                spec = write(db, DBSpectrum, det, e0, comp, person, sample, dt, "ADM-6005a[$i]", fn, "EMSA")
                write(db, NeXLDatabase.DBFitSpectrum, campaign, spec)
            end
            block1 = write(db, DBSample, 1, "Standard Block C", "NIST Standard Block C")
            block2 = write(db, DBSample, 1, "High TC Block", "NIST High Temperature Superconductor Block")
            block3 = write(db, DBSample, 1, "Copper QC", "Copper QC Standard")
            for ref in ("Al", "Fe", "Ge", "Si", "SiO2", "Ti", "Zn")
                fn = "$path\\$(ref) std.msa"
                comp = find(db, Material, ref)
                if comp < 0
                    comp = write(db, parse(Material, ref))
                end
                sample = write(db, DBSample, block1, 1, ref, "$ref standard")
                spec = write(db, DBSpectrum, det, e0, comp, person, sample, dt, "$ref std", fn, "EMSA")
                ref = write(db, NeXLDatabase.DBReference, campaign, spec, [keys(parse(Material, ref))...])
            end
            for ref in ("CaF2",)
                fn = "$path\\$(ref) std.msa"
                comp = find(db, Material, ref)
                if comp < 0
                    comp = write(db, parse(Material, ref))
                end
                sample = write(db, DBSample, block2, 1, ref, "$ref standard")
                spec = write(db, DBSpectrum, det, e0, comp, person, sample, dt, "$ref std", fn, "EMSA")
                ref = write(db, NeXLDatabase.DBReference, campaign, spec, [keys(parse(Material, ref))...])
            end
        end

        fs = read(db, NeXLDatabase.DBCampaign, 1)
        @test fs.detector.vendor == "Bruker"
        @test fs.project.name == "ADM-6005a Test"
        @test !(n"Fe" in fs.elements)
        @test n"Ge" in fs.elements
        @test length(fs.fitspectrum) == 15
        @test length(fs.refspectrum) == 8
    end
end
@testset "Fit ADM from DB"  begin
    person=read(db,DBPerson,"dale.newbury@nist.gov")
    krp = read(db, DBProject, 0, "K-ratio Project")
    project = read(db, DBProject, write(db, DBProject, "Quant Test", "Description of Quant Test", person, krp))
    lab = findall(db,DBLaboratory, person)[1]
    sample = read(db, DBSample, write(db, DBSample, lab, "ADM-6005a block", "ADM-6005a prepared by Eric Windsor"))
    detector = findall(db,DBDetector,findall(db,DBInstrument, lab)[1])[1]
    unkComp = read(db,Material,"ADM6005a")

    blkC = read(db, DBSample, write(db, DBSample, lab, "Standard Block C", "NIST Standard Block C"))
    hiTC = read(db, DBSample, write(db, DBSample, lab, "High TC Block", "NIST High Temperature Superconductor Block"))
    e0 = 20.0e3

    path = joinpath(@__DIR__, "ADM-6005a")

    cfs = write(db, DBCampaign, project, sample, unkComp, detector, person, e0,
        [ joinpath(path,"ADM-6005a_$(i).msa") for i in 1:15 ], [
        # DBSample, Union{Material, Missing}, String, Float64, Vector{Element}}
        ( blkC, parse(Material,"Al"), joinpath(path,"Al std.msa"), e0, [ n"Al" ]),
        ( blkC, parse(Material,"Fe"), joinpath(path,"Fe std.msa"), e0, [ n"Fe" ]),
        ( blkC, parse(Material,"Ge"), joinpath(path,"Ge std.msa"), e0, [ n"Ge" ]),
        ( blkC, parse(Material,"Si"), joinpath(path,"Si std.msa"), e0, [ n"Si" ]),
        ( blkC, parse(Material,"SiO2"), joinpath(path,"SiO2 std.msa"), e0, [ n"Si", n"O", n"C" ]),
        ( blkC, parse(Material,"Ti"), joinpath(path,"Ti std.msa"), e0, [ n"Ti" ]),
        ( blkC, parse(Material,"Zn"), joinpath(path,"Zn std.msa"), e0, [ n"Zn" ]),
        ( hiTC, parse(Material,"CaF2"), joinpath(path,"CaF2 std.msa"), e0, [ n"Ca", n"F", n"C" ])
        ], [ n"C" ])
    unkcomp = read(db, Material, find(db,Material,"ADM6005a"))
    ffrs = fit_spectrum(db, DBCampaign, cfs, unkcomp)

    @test isapprox(mean(value.(kratio(n"O K-L3", ffr) for ffr in ffrs)), 0.4896, atol = 0.002)
    @test isapprox(mean(value.(kratio(n"Si K-L3", ffr) for ffr in ffrs)), 0.0214, atol = 0.013)
    @test isapprox(mean(value.(kratio(n"Al K-L3", ffr) for ffr in ffrs)), 0.0281, atol = 0.002)
    @test isapprox(mean(value.(kratio(n"Ca K-L3", ffr) for ffr in ffrs)), 0.1211, rtol = 0.003)
    @test isapprox(mean(value.(kratio(n"Zn L3-M5", ffr) for ffr in ffrs)), 0.0700, rtol = 0.05)
    @test isapprox(mean(value.(kratio(n"Zn K-L3", ffr) for ffr in ffrs)), 0.1115, atol = 0.0005)
    @test isapprox(mean(value.(kratio(n"Zn K-M3", ffr) for ffr in ffrs)), 0.1197, atol = 0.002)
    @test isapprox(mean(value.(kratio(n"Ti L3-M3", ffr) for ffr in ffrs)), 0.0541, atol = 0.22)
    @test isapprox(mean(value.(kratio(n"Ti K-L3", ffr) for ffr in ffrs)), 0.064, atol = 0.001)
    @test isapprox(mean(value.(kratio(n"Ti K-M3", ffr) for ffr in ffrs)), 0.064, rtol = 0.06)
    @test isapprox(mean(value.(kratio(n"Ge L3-M5", ffr) for ffr in ffrs)), 0.1789, rtol = 0.01)
    @test isapprox(mean(value.(kratio(n"Ge K-L3", ffr) for ffr in ffrs)), 0.2628, atol = 0.001)
    @test isapprox(mean(value.(kratio(n"Ge K-M3", ffr) for ffr in ffrs)), 0.279, atol = 0.011)
end

db=nothing
GC.gc()


text:
include("utils.jl")
#just a holder for the z partitions.
#to allow split_model to work correctly

abstract type SAFTgammaMieModel <: GCSAFTModel end

struct γMieZ <: EoSModel
    components::Vector{String}
    z::SingleParam{Vector{Float64}}
end

struct SAFTgammaMieParam <: EoSParam
    segment::SingleParam{Int}
    mixedsegment::SingleParam{Vector{Float64}}
    shapefactor::SingleParam{Float64}
    lambda_a::PairParam{Float64}
    lambda_r::PairParam{Float64}
    sigma::PairParam{Float64}
    epsilon::PairParam{Float64}
    epsilon_assoc::AssocParam{Float64}
    bondvol::AssocParam{Float64}
end

@registermodel γMieZ

struct SAFTgammaMie{I,VR} <: SAFTgammaMieModel
    components::Vector{String}
    groups::GroupParam
    sites::SiteParam
    params::SAFTgammaMieParam
    idealmodel::I
    vrmodel::VR
    mie_zfractions::γMieZ
    assoc_options::AssocOptions
    references::Array{String,1}
end

function split_model(model::SAFTgammaMie,subset = nothing)
    if !(subset === nothing)
        error("SAFTgammaMie does not support custom subsets")
    end
    pures = auto_split_model(model)
    for (i,pure) in pairs(pures)
        filter!(!iszero,pure.mie_zfractions.z.values[1])
        mixm = pure.params.mixedsegment.values
        for mi in mixm
            xi = mi[i]
            resize!(mi, 1)
            mi[1] = xi
        end
    end
    return pures
end

function SAFTgammaMie(components;
    idealmodel=BasicIdeal,
    userlocations=String[],
    ideal_userlocations=String[],
    verbose=false,
    assoc_options = AssocOptions())

    groups = GroupParam(components, ["SAFT/SAFTgammaMie/SAFTgammaMie_groups.csv"]; verbose=verbose)
    params,sites = getparams(groups, ["SAFT/SAFTgammaMie","properties/molarmass_groups.csv"]; userlocations=userlocations, verbose=verbose)
    _sites = sites.i_sites
    components = groups.components
    gc = groups.i_flattenedgroups
    comps = 1:length(components)
    gc_segment = params["vst"]
    shapefactor = params["S"]
    S = shapefactor.values
    vst = gc_segment.values
    v  = groups.n_flattenedgroups
    comp_segment = zeros(Float64,length(comps))
    for i ∈ comps
        res_i = 0.0
        vi = v[i]
        groups_i = groups.i_groups[i]
        for idx ∈ 1:length(groups_i)
            k = groups_i[idx]
            res_i += vi[k]*S[k]*vst[k]
        end
        comp_segment[i] = res_i
    end
    segment = SingleParam("segment",components,comp_segment)

    #used in x_S:
    #x_S(group i) = dot(z,mixsegment[i])/dot(z,m_vr)
    mixsegment =  [[v[i][k]*vst[k]*S[k] for i ∈ comps] for k ∈ gc]
    gc_mixedsegment = SingleParam("mixed segment",groups.flattenedgroups,mixsegment,[false for i ∈ gc],String[],String[])
    gc_sigma = sigma_LorentzBerthelot(params["sigma"])
    gc_sigma.values .*= 1E-10
    gc_epsilon = epsilon_HudsenMcCoubrey(params["epsilon"], gc_sigma)
    ϵ = gc_epsilon.values
    σ = gc_sigma.values
     #helper functions
     function ẑ(i, k)
        return v[i][k]*vst[k]*S[k] / ∑(v[i][l]*vst[l]*S[l] for l ∈ gc)
    end

    zz = [[0.0 for k ∈ 1:length(v[i])] for i ∈ comps]
    zzparam = SingleParam("z fraction",components,zz,[false for i ∈ 1:length(components)],String[],String[])
    for i ∈ 1:length(zz)
        zzi = zz[i]
        for k ∈ 1:length(zzi)
            zzi[k] = ẑ(i, k)
        end
    end

    function σ̄(i)
        return cbrt(∑(∑(ẑ(i,k)*ẑ(i,l)*σ[k,l]^3 for l ∈ gc) for k ∈ gc))
    end

    function σ̄(i, j)
        return (σ̄(i) + σ̄(j))/2
    end

    function ϵ̄(i)
        return ∑(∑(ẑ(i,k)*ẑ(i,l)*ϵ[k,l] for l ∈ gc) for k ∈ gc)
    end

    function ϵ̄(i, j)
        if i == j
            return ϵ̄(i)
        else
            return sqrt(σ̄(i)*σ̄(j))/σ̄(i,j) * sqrt(ϵ̄(i)*ϵ̄(j))
        end
    end

    comp_ϵ = [ϵ̄(i, j) for (i,j) ∈ Iterators.product(comps,comps)]
    epsilon = PairParam("epsilon",components,comp_ϵ)

    comp_σ = [σ̄(i, j) for (i,j) ∈ Iterators.product(comps,comps)]
    sigma = PairParam("sigma",components,comp_σ)

    gc_lambda_a = lambda_LorentzBerthelot(params["lambda_a"])
    gc_lambda_r = lambda_LorentzBerthelot(params["lambda_r"])
    #@show gc_lambda_a
    gc_λa = gc_lambda_a.values
    gc_λr = gc_lambda_r.values

    function λi(ll,i)
        return ∑(∑(ẑ(i,k)*ẑ(i,l)*ll[k,l] for l ∈ gc) for k ∈ gc)
    end
    comp_lambda_a = [λi(gc_λa,i) for i ∈ comps]
    comp_lambda_r = [λi(gc_λr,i) for i ∈ comps]
    lambda_a = lambda_LorentzBerthelot(SingleParam("lambda_a",components,comp_lambda_a))
    lambda_r = lambda_LorentzBerthelot(SingleParam("lambda_r",components,comp_lambda_r))

    #assoc
    gc_epsilon_assoc = params["epsilon_assoc"]
    gc_bondvol = params["bondvol"]
    comp_sites,idx_dict = gc_to_comp_sites(sites,groups)
    assoc_idx = gc_to_comp_assoc_idx(gc_bondvol,comp_sites,idx_dict)
    assoc_idxs,outer,inner,outer_size,inner_size = assoc_idx.values,assoc_idx.outer_indices,assoc_idx.inner_indices,assoc_idx.outer_size,assoc_idx.inner_size
    _comp_bondvol = [gc_bondvol.values.values[i] for i ∈ assoc_idxs]
    _comp_epsilon_assoc = [gc_epsilon_assoc.values.values[i] for i ∈ assoc_idxs]
    compval_bondvol = Compressed4DMatrix(_comp_bondvol,outer,inner,outer_size,inner_size)
    compval_epsilon_assoc = Compressed4DMatrix(_comp_epsilon_assoc,outer,inner,outer_size,inner_size)
    comp_bondvol = AssocParam{Float64}("epsilon assoc",components,compval_bondvol,comp_sites.sites,String[],String[])
    comp_epsilon_assoc = AssocParam{Float64}("bondvol",components,compval_epsilon_assoc,comp_sites.sites,String[],String[])
    #mixing of Mw
    comp_mw = zeros(Float64,length(components))
    gc_mw = params["Mw"].values
    for i ∈ 1:length(components)
        vi = v[i]
        gi = groups.i_groups[i]
        mwi = 0.0
        for idx ∈ 1:length(gi)
            j = gi[idx]
            mwi += gc_mw[j]*vi[j]
        end
        comp_mw[i] =mwi
    end
    _mw = SingleParam("molecular_weight",components,comp_mw)
    gcparams = SAFTgammaMieParam(gc_segment,gc_mixedsegment, shapefactor,gc_lambda_a,gc_lambda_r,gc_sigma,gc_epsilon,gc_epsilon_assoc,gc_bondvol)
    vrparams = SAFTVRMieParam(segment,sigma,lambda_a,lambda_r,epsilon,comp_epsilon_assoc,comp_bondvol,_mw)
    idmodel = init_model(idealmodel,components,ideal_userlocations,verbose)
    vr = SAFTVRMie(vrparams, comp_sites, idmodel; ideal_userlocations, verbose, assoc_options)
    mie_z = γMieZ(components,zzparam)
    γmierefs = ["10.1063/1.4851455", "10.1021/je500248h"]
    gmie = SAFTgammaMie(components,groups,sites,gcparams,idmodel,vr,mie_z,assoc_options,γmierefs)
    return gmie
end
@registermodel SAFTgammaMie

const SAFTγMie = SAFTgammaMie
export SAFTgammaMie,SAFTγMie

include("equations.jl")


text:
"""
    mutable struct Structure_Ca  <: Str

This contains information for model constructing of C and its compounds.

# Supported list:
- structure_caco3()
"""
mutable struct Structure_Ca <: Str
    atom_vec::Matrix{Float64}
    cell_vec::Matrix{Float64}
    atom_type
    atom_name
    atom_charge
    para_mass
    num_atoms
    num_atom_types
	vec_type_id
	res_name
end

"""
    structure_caco3()

Do this will generate a Structure_Ca type which contains all information needed to build a CaCO3 moedel.
"""
function structure_caco3()
    atom_vec = [
		0.772116  0.360801   0.0909051
		0.32908   0.721614   0.454547
		0.329085  0.0        0.818179
		0.329085  0.0        0.272726
		0.772116  0.360801   0.636358
		0.32908   0.721614   1.0
		0.329085  0.0        0.0
		0.772116  0.360801   0.363631
		0.32908   0.721614   0.727274
		0.32908   0.721614   0.181821
		0.329085  0.0        0.545453
		0.772116  0.360801   0.909084
		1.0       0.360801   0.0909051
		0.556973  0.721614   0.454547
		0.556969  0.0        0.818179
		0.658165  0.0824151  0.0909051
		0.215143  0.443217   0.454547
		0.658169  0.804029   0.818179
		0.658169  0.639199   0.0909051
		0.215138  1.0        0.454547
		0.215143  0.278386   0.818179
		0.987254  0.0        0.272726
		0.544223  0.360801   0.636358
		0.101196  0.721614   1.0
		0.443027  0.278386   0.272726
		0.886054  0.639199   0.636358
		0.443031  1.0        1.0
		0.0       0.804029   0.272726
		0.886058  0.0824151  0.636358
		0.443027  0.443217   1.0
    ]
    cell_vec = [
        5.69943  0.0       0.0
		0.0      4.04044   0.0
		0.0      0.0      15.8088
    ]
    atom_type = vcat([1 for n=1:6], [2 for n=1:6], [3 for n=1:18])
    atom_charge = 0 .* atom_type
    para_mass = [40.078, 12.011, 15.999]
    atom_name = split("Ca C O")
    num_atoms = length(atom_type)
	num_atom_types = length(atom_name)
	res_name = "CaCO"
    Structure_Ca(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms, num_atom_types, [1], res_name)
end


text:
using NumericalMethodsforEngineers

println("\n'Fortran' version:\n")
a = [16.; 4.; 5.; 8.; -4.; 22.]
kdiag = [1; 3; 6]

amat = fromSkyline(a, kdiag)
amat |> display
println()

b = [4., 2., 5.]

nfix = 1
no = [2;]
val = [5.0;]
penalty = 1.0e20

a[kdiag[no]] .+= penalty
b[no] = a[kdiag[no]] .* val

sparin!(a, kdiag)
spabac!(a, b, kdiag)
d = copy(b)
println("Solution Vector: \n $d")

#------------------------------------

println("\n'Julia' version:\n")
a = amat
b = [4., 2., 5.]

a[2, 2] += 1.0e20
b[2] = a[2, 2] * 5.0

a = sparse(a)

F = cholfact(a)
upper = F[:U]
upper |> display

lower = F[:L]
lower |> display

y = lower \ b

c = upper \ y
println("Solution Vector: \n $c")

@test c == a \ b
@test c == d


text:
#
#
#  UNTIL THE IniFile  PKG GETS UPDATED METADATA PASTE IT HERE
#
#

import Base.get,
       Base.haskey,
       Base.read,
       Base.write,
       Base.show

typealias INIVAL Union(String,Number,Bool,Nothing)
typealias HTSS Dict{String,INIVAL}

type Inifile
    sections::Dict{String,HTSS}
    defaults::HTSS
end

Inifile() = Inifile((String=>HTSS)[], HTSS())

defaults(inifile::Inifile) = inifile.defaults

sections(inifile::Inifile) = inifile.sections

function read(inifile::Inifile, stream::IO)
    current_section = inifile.defaults
    for line in EachLine(stream)
        s = strip(line)
        # comments start with # or ;
        if length(s) < 3 || s[1] == '#' || s[1] == ';'
            continue
        elseif s[1] == '[' && s[end] == ']'
            section = s[2:end-1]
            if !haskey(inifile.sections, section)
                inifile.sections[section] = HTSS()
            end
            current_section = inifile.sections[section]
        else
            i = search(s, '=')
            j = search(s, ':')
            if i == 0 && j == 0
                # TODO: allow multiline values
                println("skipping malformed line: $s")
            else
                idx = min(i, j)
                if idx == 0
                    idx = max(i, j)
                end
                key = rstrip(s[1:idx-1])
                val = lstrip(s[idx+1:end])
                current_section[key] = val
            end
        end
    end
    inifile
end

function read(inifile::Inifile, filename::String)
    open(filename) do f
        read(inifile, f)
    end
    inifile
end

show(io::IO, inifile::Inifile) = write(io, inifile)
function write(io::IO, inifile::Inifile)
    for (key, value) in defaults(inifile)
        println(io, "$key=$value")
    end
    for (name, htss) in sections(inifile)
        println(io, "[$name]")
        for (key, value) in htss
            println(io, "$key=$value")
        end
    end
end

get(inifile::Inifile, section::String, key::String) = get(inifile, section, key, :notfound)

function get(inifile::Inifile, section::String, key::String, notfound)
    if haskey(inifile.sections, section) && haskey(inifile.sections[section], key)
        return inifile.sections[section][key]
    elseif haskey(inifile.defaults, key)
        return inifile.defaults[key]
    end
    notfound
end

   # This needs to be fixed to set the "defaults" when necessary
   # If there is nothing in the section just set the default
function set(inifile::Inifile, section::String, key::String, val::INIVAL)
    if !haskey(inifile.sections, section)
        (val == nothing) && return val
        inifile.sections[section] = HTSS()
    end

    sec = inifile.sections[section]
    if val == nothing
        if haskey(sec, key)
            delete!(sec, key)
        end
        return val
    end
    sec[key] = val
    val
end

#
#
#  ===============================================================================
#
#

function LOCAL_STORE_PATH_STATUS( )
    ini = Inifile()
    conf_path = Pkg.dir("MethylUtils", "conf", "conf.ini")
    read(ini,conf_path)
    local_store_path = get(ini,"","LOCAL_STORE_PATH")
    if(local_store_path == "nothing")
         return "unset"
    else
         return "set"
    end
end

function LOCAL_STORE_PATH()
    ini = Inifile()
    conf_path = Pkg.dir("MethylUtils", "conf", "conf.ini")
    read(ini,conf_path)
    #println("$conf_path ",ini)
    local_store_path = get(ini,"","LOCAL_STORE_PATH")
    #println("GOT PATH $local_store_path")
    if( local_store_path == "nothing" )
       println("going to set local path")
       local_store_path = SET_LOCAL_STORE_PATH()
    end

    return local_store_path
end

function SET_LOCAL_STORE_PATH(path="")
     default_path = joinpath(homedir(), ".methylUtil", "annotation")

     if( path == "" )
         path = default_path
     end

     conf_path = Pkg.dir("MethylUtils", "conf", "conf.ini")
     ini = Inifile()
     read(ini,conf_path)
     println("Setting $path")
     ini.defaults["LOCAL_STORE_PATH"] = path
     io = open(conf_path,"w")
     write(io, ini)
     close(io)

     return path
end

function current_annotation_chr_names()
    CURRENT_ANNOTATION_ORG_NAMES = load_current_annotation_organisms()
    CURRENT_ANNOTATION_CHR_NAMES = load_current_annotation_chr_names(CURRENT_ANNOTATION_ORG_NAMES)
end

function current_annotation_chr_index()
    CURRENT_ANNOTATION_ORG_NAMES = load_current_annotation_organisms()
    CURRENT_ANNOTATION_CHR_NAMES = load_current_annotation_chr_index(CURRENT_ANNOTATION_ORG_NAMES)
end

function load_current_annotation_organisms()
     # reads a file with a list of genomes on system
    path = joinpath( LOCAL_STORE_PATH(), "organism.names")
    if(!isfile(path))
       error("Expected an organism.names file to exist - path trying to access here: $path")
    end

    iostream = open(path)
    orgs = Array(ASCIIString, 0)
    for line in eachline(iostream)
        push!(orgs, chomp(line) )
    end
    close(iostream)
    return orgs
end

function load_current_annotation_chr_names(orgs)
         # We want to load
         # CURRENT_ANNOTATION_CHR_NAMES
         # CURRENT_ANNOTATION_CHR_INDEX
         # for each folder - look for chr.names and load
         # joinpath(LOCAL_STORE_PATH(),
    org_chr_names = Dict()
    for org in orgs
        path = joinpath(LOCAL_STORE_PATH(), org,"chr.names")
        iostream = open(path)
        if(!isfile(path))
            error("Expected a chr.names file to exist for $org - path trying to access here: $path")
        end
        chrs = Array(ASCIIString,0)
        for line in eachline(iostream)
            #trim whitespace
            push!(chrs, chomp( line ) )
        end
        org_chr_names[org] = chrs
    end
    return org_chr_names
end

function load_current_annotation_chr_index(orgs)

    CURRENT_ANNOTATION_CHR_NAMES = load_current_annotation_chr_names(orgs)
    CURRENT_ANNOTATION_CHR_INDEX = Dict()
    for organism_id in keys(CURRENT_ANNOTATION_CHR_NAMES)
        chromnames = CURRENT_ANNOTATION_CHR_NAMES[organism_id]
        h = Dict()
        idx = 1
        for chr in chromnames
            h[chr] = idx
            idx +=1
        end
        CURRENT_ANNOTATION_CHR_INDEX[organism_id] = h
        idx=1
    end
    return CURRENT_ANNOTATION_CHR_INDEX
end

#LOCAL_CURRENT_ANNOTATIONS = joinpath(LOCAL_STORE_PATH(),"current_annotation.jl")
#include(LOCAL_CURRENT_ANNOTATIONS)


text:
# ===============================
# Written by AAE
# TU Delft, Spring 2014
# simulkade.com
# Last edited: 30 December, 2014
# ===============================

# =====================================================================
# 2014-12-30: - cell variables can be created including bounddary cells
# =====================================================================

# %% ===================== Cell Variables ==================================
"""
Creates a cell variable on the mesh and assigns value `phi0` to it:

```julia
var = createCellVariable(mesh, phi0)
```

Inputs:

    + mesh: is a mesh structure created by one of `create*Mesh` functions
    + phi0 (real): is an initial value of the variable on each cell

Output:

    + var: cell variable

"""
function createCellVariable(m::MeshStructure, phi0::Real)
    CellValue(m, phi0*ones(tuple(m.dims.+2...)))
end

"""
Creates a cell variable on the mesh and assigns value `phi0` to it.
Cell variable is forced to satisfy provided boundary conditions

```julia
var = createCellVariable(mesh, phi0, bc)
```

Inputs:

    + mesh: is a mesh structure created by one of `create*Mesh` functions
    + phi0 (real): is an initial value of the variable on each cell
    + bc: boundary conditions created on the mesh

Output:

    + var: cell variable

"""
function createCellVariable(m::MeshStructure, phi0::Real, BC::BoundaryCondition)
    phi = CellValue(m, phi0*ones(tuple(m.dims.+2...)))
    cellBoundary!(phi, BC)
end

"""
Creates a cell variable on the mesh from the existing values in the
array that follows the mesh structure shape:

```julia
var = createCellVariable(mesh, array)
```

Inputs:
    + mesh: is a mesh structure created by one of `create*Mesh` functions
    + array (array of reals): array of initial values

Output:

    + var: cell variable
"""
function createCellVariable{T<:Real}(m::MeshStructure, phi0::Array{T})
    if prod(m.dims+2)==length(phi0)
      CellValue(m, phi0)
    elseif prod(m.dims)==length(phi0)
      d=m.dimension
      phival = zeros(tuple(m.dims.+2...))
      BC = createBC(m) # Neumann boundaries
      if d==1 || d==1.5
        phival[2:end-1] = phi0
      elseif d==2 d==2 || d==2.5 || d==2.8
        phival[2:end-1, 2:end-1] = phi0
      elseif d==3 || d==3.2
        phival[2:end-1,2:end-1,2:end-1] = phi0
      end
      phi = CellValue(m, phival)
      cellBoundary!(phi, BC)
    else
      error("JFVMM: Matrix must be the same size as the domain.")
    end
end

"""
Creates a cell variable on the mesh from the existing values in the
array that follows the mesh structure shape. The cell variable is forced to
satisfy provided boundary conditions:

```julia
var = createCellVariable(mesh, array, bc)
```

Inputs:
    + mesh: is a mesh structure created by one of `create*Mesh` functions
    + array (array of reals): array of initial values
    + bc: boundary conditions on the mesh

Output:

    + var: cell variable
"""
function createCellVariable{T<:Real}(m::MeshStructure, phi0::Array{T}, BC::BoundaryCondition)
    if prod(m.dims+2)==length(phi0)
      error("jFVT: Matrix must be the same size as the domain.")
    elseif prod(m.dims)==length(phi0)
      d=m.dimension
      phival = zeros(tuple(m.dims.+2...))
      if d==1 || d==1.5
        phival[2:end-1] = phi0
      elseif d==2 || d==2.5 || d==2.8
        phival[2:end-1, 2:end-1] = phi0
      elseif d==3 || d==3.2
        phival[2:end-1,2:end-1,2:end-1] = phi0
      end
      phi = CellValue(m, phival)
      cellBoundary!(phi, BC)
    else
      error("JFVMM: Matrix must be the same size as the domain.")
    end
end

# %% ============== Face Variavle ======================================
"""
Creates a face variable (vector) on the mesh and assigns value `phi0` to it:

```julia
var = createFaceVariable(mesh, phi0)
```

Inputs:

    + mesh: is a mesh structure created by one of `create*Mesh` functions
    + phi0 (vector): is an initial value of the variable on each face;
        the dimension of phi0 must be the same as the mesh (1,2 or 3)

Output:

    + var: face variable

"""
function createFaceVariable{T<:Real}(m::MeshStructure, phi0::Array{T,1})
d=m.dimension
  if d==1 || d==1.5
    FaceValue(m,
	      ones(m.dims[1]+1)*phi0[1],
	      [1.0],
	      [1.0])
  elseif d==2 || d==2.5 || d==2.8
    FaceValue(m,
	      ones(m.dims[1]+1, m.dims[2])*phi0[1],
	      ones(m.dims[1], m.dims[2]+1)*phi0[2],
	      [1.0])
  elseif d==3 || d==3.2
    FaceValue(m,
	      ones(m.dims[1]+1, m.dims[2], m.dims[3])*phi0[1],
	      ones(m.dims[1], m.dims[2]+1, m.dims[3])*phi0[2],
	      ones(m.dims[1], m.dims[2], m.dims[3]+1)*phi0[3])
  end
end

"""
Creates a face variable (vector) on the mesh and assigns value `phi0` to it:

```julia
var = createFaceVariable(mesh, phi0)
```

Inputs:

    + mesh: is a mesh structure created by one of `create*Mesh` functions
    + phi0 (real): is an initial value of the variable on each face;
        each active direction in the face value will get this value

Output:

    + var: face variable

"""
function createFaceVariable(m::MeshStructure, phi0::Real)
# creates a face variable based on the mesh structure
d=m.dimension
  if d==1 || d==1.5
    FaceValue(m,
	      ones(m.dims[1]+1)*phi0,
	      [1.0],
	      [1.0])
  elseif d==2 || d==2.5 || d==2.8
    FaceValue(m,
	      ones(m.dims[1]+1, m.dims[2])*phi0,
	      ones(m.dims[1], m.dims[2]+1)*phi0,
	      [1.0])
  elseif d==3 || d==3.2
    FaceValue(m,
	      ones(m.dims[1]+1, m.dims[2], m.dims[3])*phi0,
	      ones(m.dims[1], m.dims[2]+1, m.dims[3])*phi0,
	      ones(m.dims[1], m.dims[2], m.dims[3]+1)*phi0)
  end
end

# %% ============== Cell Vector =====================================
# TODO: docstring for createCellVector
function createCellVector{T<:Real}(m::MeshStructure, phi0::Array{T,1})
    d=m.dimension
    if d==1 || d==1.5
        CellVector(m,
    	      ones(m.dims[1])*phi0[1],
    	      [1.0],
    	      [1.0])
    elseif d==2 || d==2.5 || d==2.8
        CellVector(m,
    	      ones(m.dims[1], m.dims[2])*phi0[1],
    	      ones(m.dims[1], m.dims[2])*phi0[2],
    	      [1.0])
    elseif d==3 || d==3.2
        CellVector(m,
    	      ones(m.dims[1], m.dims[2], m.dims[3])*phi0[1],
    	      ones(m.dims[1], m.dims[2], m.dims[3])*phi0[2],
    	      ones(m.dims[1], m.dims[2], m.dims[3])*phi0[3])
    end
end

function createCellVector(m::MeshStructure, phi0::Real)
    d=m.dimension
    if d==1 || d==1.5
        CellVector(m,
    	      ones(m.dims[1])*phi0,
    	      [1.0],
    	      [1.0])
    elseif d==2 || d==2.5 || d==2.8
        CellVector(m,
    	      ones(m.dims[1], m.dims[2])*phi0,
    	      ones(m.dims[1], m.dims[2])*phi0,
    	      [1.0])
    elseif d==3 || d==3.2
        CellVector(m,
    	      ones(m.dims[1], m.dims[2], m.dims[3])*phi0,
    	      ones(m.dims[1], m.dims[2], m.dims[3])*phi0,
    	      ones(m.dims[1], m.dims[2], m.dims[3])*phi0)
    end
end

# %% Copy cell
"""
Copy cell value
"""
copyCell(phi::CellValue) = CellValue(phi.domain, Base.copy(phi.value))


text:
immutable BoundingBox{P1 <: Vec, P2 <: Vec}
    x0::P1
    a::P2
end

BoundingBox{P <: Vec2, T1 <: Measure, T2 <: Measure}(x0::P, width::T1, height::T2) =
    BoundingBox{P, Tuple{T1, T2}}(x0, (width, height))

BoundingBox{P <: Vec3, T1 <: Measure, T2 <: Measure, T3 <: Measure}(x0::P, width::T1, height::T2, depth::T3) =
    BoundingBox{P, Tuple{T1, T2, T3}}(x0, (width, height, depth))

#BoundingBox{X}(x0::Vec{2, X}, width::Measure, height::Measure) =
    #BoundingBox{2, X, Measure}(x0, (width, height))

BoundingBox(x0::Measure, y0::Measure, width::Measure, height::Measure) =
    BoundingBox((x0, y0), width, height)

BoundingBox(x0::Measure, y0::Measure, z0::Measure, width::Measure, height::Measure, depth::Measure) =
    BoundingBox((x0, y0, z0), width, height, depth)

BoundingBox() = BoundingBox(0mm, 0mm, 1w, 1h)
BoundingBox(width, height) = BoundingBox(0mm, 0mm, width, height)
BoundingBox(width, height, depth) = BoundingBox(0mm, 0mm, 0mm, width, height, depth)

isabsolute{P1, P2}(b::BoundingBox{P1, P2}) = isabsolute(b.x0) && isabsolute(b.a)

typealias AbsoluteBox{N}   BoundingBox{NTuple{N, Length{:mm, Float64}},
                                       NTuple{N, Length{:mm, Float64}}}
typealias Absolute2DBox    BoundingBox{Tuple{AbsoluteLength, AbsoluteLength},
                                       Tuple{AbsoluteLength, AbsoluteLength}}
typealias Absolute3DBox    BoundingBox{
                                Tuple{
                                    AbsoluteLength,
                                    AbsoluteLength,
                                    AbsoluteLength
                                },
                                Tuple{
                                    AbsoluteLength,
                                    AbsoluteLength,
                                    AbsoluteLength
                                }
                            }

width(x::BoundingBox)  = x.a[1]
height(x::BoundingBox) = x.a[2]
depth(x::BoundingBox)  = x.a[3]


text:
__precompile__()

module BufferedStreams

export BufferedInputStream,
       BufferedOutputStream,
       peek,
       peekbytes!,
       fillbuffer!,
       isanchored,
       anchor!,
       upanchor!,
       takeanchored!

using Compat
import Compat: view

# default buffer size is 128 KiB
const default_buffer_size = 128 * 2^10

# pretty printer of data size
function _datasize(nbytes)
    if nbytes < 1024
        return string(nbytes, " B")
    end
    k = 1
    for suffix in [" KiB", " MiB", " GiB", " TiB"]
        if nbytes < 1024^(k+1)
            return string(floor(nbytes / 1024^k, 1), suffix)
        end
        k += 1
    end
    return string(floor(nbytes / 1024^k, 1), " PiB")
end

include("bufferedinputstream.jl")
include("bufferedoutputstream.jl")
include("emptystream.jl")
include("io.jl")
include("deprecated.jl")

end # module BufferedStreams


text: # TODO: implement an ARMA-based statistical model.


text:
# Glove model
struct Model{T<:AbstractFloat}
    W_main::Matrix{T}
    W_ctx::Matrix{T}
    b_main::Vector{T}
    b_ctx::Vector{T}
    W_main_grad::Matrix{T}
    W_ctx_grad::Matrix{T}
    b_main_grad::Vector{T}
    b_ctx_grad::Vector{T}
end

"""
Each vocab word in associated with a word vector and a context vector.
The paper initializes the weights to values [-0.5, 0.5] / vecsize+1 and
the gradients to 1.0.

The +1 term is for the bias.
"""
function Model(::Type{T}, vocabsize::Int, vecsize::Int) where T<:AbstractFloat
    shift = T(0.5)
    Model(
        (rand(T, vecsize, vocabsize) .- shift) ./ T(vecsize + 1),
        (rand(T, vecsize, vocabsize) .- shift) ./ T(vecsize + 1),
        (rand(T, vocabsize) .- shift) ./ T(vecsize + 1),
        (rand(T, vocabsize) .- shift) ./ T(vecsize + 1),
        ones(T, vecsize, vocabsize),
        ones(T, vecsize, vocabsize),
        ones(T, vocabsize),
        ones(T, vocabsize),
    )
end

function adagrad!(
  ::Type{T},
  m::Model{T},
  covec::CooccurenceVector{T},
  epochs::Int;
  lrate=T(1e-2),
  xmax::Int=100,
  alpha=T(0.75)) where T<:AbstractFloat

    # store the cost per iteration
    costs = zeros(T, epochs)
    vecsize = size(m.W_main, 1)

    for n = 1:epochs
        for i = 1:length(covec)
            @inbounds pair = covec[i]
            midx = pair[1] # main index
            cidx = pair[2] # context index
            val = pair[3] # value

            diff = m.b_main[midx] + m.b_ctx[cidx] - log(val)
            @inbounds @simd for j = 1:vecsize
              diff += m.W_main[j, midx] * m.W_ctx[j, cidx]
            end

            fdiff = ifelse(val < xmax, (val / xmax) ^ alpha, T(1)) * diff
            costs[n] = T(0.5) * fdiff * diff

            # Adaptive learning gradient updates
            fdiff *= lrate
            @inbounds @simd for j = 1:vecsize
                tm = fdiff * m.W_ctx[j, cidx]
                tc = fdiff * m.W_main[j, midx]
                m.W_main[j, midx] -= tm / sqrt(m.W_main_grad[j, midx])
                m.W_ctx[j, cidx] -= tc / sqrt(m.W_ctx_grad[j, cidx])
                m.W_main_grad[j, midx] += tm * tm
                m.W_ctx_grad[j, cidx] += tc * tc
            end

            # bias updates
            m.b_main[midx] -= fdiff / sqrt(m.b_main_grad[midx])
            m.b_ctx[cidx] -= fdiff / sqrt(m.b_ctx_grad[cidx])
            fdiff *= fdiff
            m.b_main_grad[midx] += fdiff
            m.b_ctx_grad[cidx] += fdiff
        end
    end
    costs
end


text:
using Documenter, MiniBee

makedocs(;
    modules=[MiniBee],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/serenity4/MiniBee.jl/blob/{commit}{path}#L{line}",
    sitename="MiniBee.jl",
    authors="serenity4 <cedric.bel@hotmail.fr>",
    assets=String[],
)


text:
# ---
# title: 1514. Path with Maximum Probability
# id: problem1514
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Graph
# link: <https://leetcode.com/problems/path-with-maximum-probability/description/>
# hidden: true
# ---
#
# You are given an undirected weighted graph of `n` nodes (0-indexed),
# represented by an edge list where `edges[i] = [a, b]` is an undirected edge
# connecting the nodes `a` and `b` with a probability of success of traversing
# that edge `succProb[i]`.
#
# Given two nodes `start` and `end`, find the path with the maximum probability
# of success to go from `start` to `end` and return its success probability.
#
# If there is no path from `start` to `end`, **return  0**. Your answer will be
# accepted if it differs from the correct answer by at most **1e-5**.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2019/09/20/1558_ex1.png)**
#
#
#
#     Input: n = 3, edges = [[0,1],[1,2],[0,2]], succProb = [0.5,0.5,0.2], start = 0, end = 2
#     Output: 0.25000
#     Explanation:  There are two paths from start to end, one having a probability of success = 0.2 and the other has 0.5 * 0.5 = 0.25.
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2019/09/20/1558_ex2.png)**
#
#
#
#     Input: n = 3, edges = [[0,1],[1,2],[0,2]], succProb = [0.5,0.5,0.3], start = 0, end = 2
#     Output: 0.30000
#
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2019/09/20/1558_ex3.png)**
#
#
#
#     Input: n = 3, edges = [[0,1]], succProb = [0.5], start = 0, end = 2
#     Output: 0.00000
#     Explanation:  There is no path between 0 and 2.
#
#
#
#
# **Constraints:**
#
#   * `2 <= n <= 10^4`
#   * `0 <= start, end < n`
#   * `start != end`
#   * `0 <= a, b < n`
#   * `a != b`
#   * `0 <= succProb.length == edges.length <= 2*10^4`
#   * `0 <= succProb[i] <= 1`
#   * There is at most one edge between every two nodes.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
export scalar_solvers

scalar_solvers = []

#ARC methods
push!(scalar_solvers, ARC_Nwt)
push!(scalar_solvers, ARC_Cub)
push!(scalar_solvers, ARC_Sec)
push!(scalar_solvers, ARC_SecA)

#TR methods
push!(scalar_solvers, TR_Nwt)
push!(scalar_solvers, TR_Cub)
push!(scalar_solvers, TR_Sec)
push!(scalar_solvers, TR_SecA)

#zoom methods
push!(scalar_solvers, zoom_Nwt)
push!(scalar_solvers, zoom_Cub)
push!(scalar_solvers, zoom_Sec)
push!(scalar_solvers, zoom_SecA)

#find interval and bissection methods
push!(scalar_solvers, bissect)
push!(scalar_solvers, bissect_Cub)
push!(scalar_solvers, bissect_nwt)
push!(scalar_solvers, bissect_sec)
push!(scalar_solvers, bissect_secA)


text:
"""
A generic, standalone implementation of Monte Carlo Tree Search.
It can be used on any game that implements `GameInterface`
and with any external oracle.

Both a synchronous and an asynchronous version are implemented, which
share most of their code. When browsing the sources for the first time,
we recommend that you study the sychronous version first.
"""
module MCTS

using Distributions: Categorical, Dirichlet

using ..Util: @printing_errors, apply_temperature
import ..GI, ..GameType

#####
##### Interface for External Oracles
#####

"""
    MCTS.Oracle{Game}

Abstract base type for an oracle. Oracles must implement
[`MCTS.evaluate`](@ref) and [`MCTS.evaluate_batch`](@ref).
"""
abstract type Oracle{Game} end

"""
    MCTS.evaluate(oracle::Oracle, state)

Evaluate a single state from the current player's perspective.

Return a pair `(P, V)` where:

  - `P` is a probability vector on `GI.available_actions(Game(state))`
  - `V` is a scalar estimating the value or win probability for white.
"""
function evaluate end

"""
    MCTS.evaluate_batch(oracle::Oracle, states)

Evaluate a batch of states.

Expect a vector of states and return a vector of `(P, V)` pairs.

A default implementation is provided that calls [`MCTS.evaluate`](@ref)
sequentially on each position.
"""
function evaluate_batch(oracle::Oracle, states)
  return [evaluate(oracle, b) for b in states]
end

GameType(::Oracle{Game}) where Game = Game

#####
##### Some Example Oracles
#####

"""
    MCTS.RolloutOracle{Game}(γ=1.) <: MCTS.Oracle{Game}

This oracle estimates the value of a position by simulating a random game
from it (a rollout). Moreover, it puts a uniform prior on available actions.
Therefore, it can be used to implement the "vanilla" MCTS algorithm.
"""
struct RolloutOracle{Game} <: Oracle{Game}
  gamma :: Float64
  RolloutOracle{G}(γ=1.) where G = new{G}(γ)
end

function rollout!(game, γ=1.)
  r = 0.
  while !GI.game_terminated(game)
    action = rand(GI.available_actions(game))
    GI.play!(game, action)
    r = γ * r + GI.white_reward(game)
  end
  return r
end

function evaluate(r::RolloutOracle{Game}, state) where Game
  game = Game(state)
  wp = GI.white_playing(game)
  n = length(GI.available_actions(game))
  P = ones(n) ./ n
  wr = rollout!(game, r.gamma)
  V = wp ? wr : -wr
  return P, V
end

struct RandomOracle{Game} <: Oracle{Game} end

function evaluate(::RandomOracle{Game}, state) where Game
  s = Game(state)
  n = length(GI.available_actions(s))
  P = ones(n) ./ n
  V = 0.
  return P, V
end

#####
##### Workers and Queries
#####

# Workers can send three types of messages to the inference server:
# - Done: there is no more work available for the worker
# - Query: the worker requests the server to evaluate a position
# - NoQuery: the worker finished an MCTS simulation without sending
#   an evaluation query (see `worker_yield!` for why this matters)

abstract type Message{S} end
struct Done{S} <: Message{S} end
struct NoQuery{S} <: Message{S} end
struct Query{S} <: Message{S} state :: S end
const AnyMessage{S} = Union{Done{S}, NoQuery{S}, Query{S}}
const EvaluationResult{R} = Tuple{Vector{R}, R}

mutable struct Worker{State}
  id :: Int # useful for debugging purposes
  send :: Channel{AnyMessage{State}}
  recv :: Channel{EvaluationResult{Float32}}
  queried :: Bool # the worker queried the server during the current simulation

  function Worker{S}(id) where S
    send = Channel{AnyMessage{S}}(1)
    recv = Channel{EvaluationResult{Float32}}(1)
    new{S}(id, send, recv, false)
  end
end

#####
##### State Statistics
#####

struct ActionStats
  P :: Float32 # Prior probability as given by the oracle
  W :: Float64 # Cumulated Q-value for the action (Q = W/N)
  N :: Int # Number of times the action has been visited
  nworkers :: UInt16 # Number of workers currently exploring this branch
end

struct StateInfo
  stats :: Vector{ActionStats}
  Vest  :: Float32 # Value estimate given by the oracle
end

Ntot(b::StateInfo) = sum(s.N for s in b.stats)

#####
##### MCTS Environment
#####

"""
    MCTS.Env{Game}(oracle; <keyword args>) where Game

Create and initialize an MCTS environment with a given `oracle`.

## Keyword Arguments

  - `nworkers=1`: numbers of asynchronous workers (see below)
  - `fill_batches=false`: if true, a constant batch size is enforced for
     evaluation requests, by completing batches with dummy entries if necessary
  - `gamma=1.`: the reward discount factor
  - `cpuct=1.`: exploration constant in the UCT formula
  - `noise_ϵ=0., noise_α=1.`: parameters for the dirichlet exploration noise
     (see below)
  - `prior_temperature=1.`: temperature to apply to the oracle's output
     to get the prior probability vector used by MCTS.

## Asynchronous MCTS

  - If `nworkers == 1`, MCTS is run in a synchronous fashion and the oracle is
    invoked through [`MCTS.evaluate`](@ref).

  - If `nworkers > 1`, `nworkers` asynchronous workers are spawned,
    along with an additional task to serve state evaluation requests.
    Such requests are processed by batches of
    size `nworkers` using [`MCTS.evaluate_batch`](@ref).

## Dirichlet Noise

A naive way to ensure exploration during training is to adopt an ϵ-greedy
policy, playing a random move at every turn instead of using the policy
prescribed by [`MCTS.policy`](@ref) with probability ϵ.
The problem with this naive strategy is that it may lead the player to make
terrible moves at critical moments, thereby biasing the policy evaluation
mechanism.

A superior alternative is to add a random bias to the neural prior for the root
node during MCTS exploration: instead of considering the policy ``p`` output
by the neural network in the UCT formula, one uses ``(1-ϵ)p + ϵη`` where ``η``
is drawn once per call to [`MCTS.explore!`](@ref) from a Dirichlet distribution
of parameter ``α``.
"""
mutable struct Env{Game, State, Oracle}
  # Store (nonterminal) state statistics assuming player one is to play
  tree :: Dict{State, StateInfo}
  # External oracle to evaluate positions
  oracle :: Oracle
  # Workers
  workers :: Vector{Worker{State}}
  global_lock :: ReentrantLock # This impl. is asynchronous but sequential
  remaining :: Int # counts the number of remaining simulations to do
  # Parameters
  fill_batches :: Bool
  gamma :: Float64 # Discount factor
  cpuct :: Float64
  noise_ϵ :: Float64
  noise_α :: Float64
  prior_temperature :: Float64
  # Performance statistics
  total_time :: Float64
  inference_time :: Float64
  total_simulations :: Int64
  total_nodes_traversed :: Int64

  function Env{G}(oracle;
      nworkers=1, fill_batches=false, gamma=1.,
      cpuct=1., noise_ϵ=0., noise_α=1., prior_temperature=1.) where G
    S = GI.State(G)
    tree = Dict{S, StateInfo}()
    total_time = 0.
    inference_time = 0.
    total_simulations = 0
    total_nodes_traversed = 0
    lock = ReentrantLock()
    remaining = 0
    workers = [Worker{S}(i) for i in 1:nworkers]
    new{G, S, typeof(oracle)}(
      tree, oracle, workers, lock, remaining, fill_batches, gamma,
      cpuct, noise_ϵ, noise_α, prior_temperature,
      total_time, inference_time, total_simulations, total_nodes_traversed)
  end
end

asynchronous(env::Env) = length(env.workers) > 1

GameType(::Env{Game}) where Game = Game

Done(::Env{G,S}) where {G,S} = Done{S}()
NoQuery(::Env{G,S}) where {G,S} = NoQuery{S}()

#####
##### Access and initialize state information
#####

function init_state_info(P, V, prior_temperature)
  P = apply_temperature(P, prior_temperature)
  stats = [ActionStats(p, 0, 0, 0) for p in P]
  return StateInfo(stats, V)
end

# Returns statistics for the current player, along with a boolean indicating
# whether or not a new node has been created.
# Synchronous version
function state_info_sync(env, worker, state)
  if haskey(env.tree, state)
    return (env.tree[state], false)
  else
    (P, V), time = @timed evaluate(env.oracle, state)
    env.inference_time += time
    info = init_state_info(P, V, env.prior_temperature)
    env.tree[state] = info
    return (info, true)
  end
end

# Equivalent of `state_info_sync` for asynchronous MCTS
function state_info_async(env, worker, state)
  if haskey(env.tree, state)
    return (env.tree[state], false)
  else
    # Send a request to the inference server
    put!(worker.send, Query(state))
    unlock(env.global_lock)
    P, V = take!(worker.recv)
    lock(env.global_lock)
    worker.queried = true
    # Another worker may have sent the same request and initialized
    # the node before. Therefore, we have to test membership again.
    if !haskey(env.tree, state)
      info = init_state_info(P, V, env.prior_temperature)
      env.tree[state] = info
      return (info, true)
    else
      # The inference result is ignored and we proceed as if
      # the node was already in the tree.
      return (env.tree[state], false)
    end
  end
end

function state_info(env, worker, state)
  asynchronous(env) ?
    state_info_async(env, worker, state) :
    state_info_sync(env, worker, state)
end

#####
##### Main algorithm
#####

function uct_scores(info::StateInfo, cpuct, ϵ, η)
  @assert iszero(ϵ) || length(η) == length(info.stats)
  sqrtNtot = sqrt(Ntot(info))
  return map(enumerate(info.stats)) do (i, a)
    Q = (a.W - a.nworkers) / max(a.N, 1)
    P = iszero(ϵ) ? a.P : (1-ϵ) * a.P + ϵ * η[i]
    Q + cpuct * P * sqrtNtot / (a.N + 1)
  end
end

function increment_visit_counter!(env, state, action_id)
  stats = env.tree[state].stats
  astats = stats[action_id]
  stats[action_id] = ActionStats(
    astats.P, astats.W, astats.N + 1, astats.nworkers + 1)
end

# Also decreases the visit count
function update_state_info!(env, state, action_id, q)
  stats = env.tree[state].stats
  astats = stats[action_id]
  stats[action_id] = ActionStats(
    astats.P, astats.W + q, astats.N, astats.nworkers - 1)
end

# It is important to guarantee that a worker sends one request to the
# inference server per simulation. Otherwise, a worker could block all
# the others if it repeatedly doesn't need to make evaluation queries.
function worker_yield!(env::Env, worker::Worker)
  if asynchronous(env) && !worker.queried
    put!(worker.send, NoQuery(env))
    unlock(env.global_lock)
    take!(worker.recv)
    lock(env.global_lock)
  end
end

# Run a single MCTS simulation, updating the statistics of all traversed states.
# Return the estimated Q-value for the current player.
# Leave the game unchanged (a copy is made while visiting the root).
function worker_explore!(env::Env, worker::Worker, game, η, root=true)
  if root
    worker.queried = false
    env.total_simulations += 1
    game = copy(game)
  end
  if GI.game_terminated(game)
    worker_yield!(env, worker)
    return 0.
  else
    state = GI.current_state(game)
    actions = GI.available_actions(game)
    info, new_node = state_info(env, worker, state)
    if new_node
      worker_yield!(env, worker)
      return info.Vest
    else
      ϵ = root ? env.noise_ϵ : 0.
      scores = uct_scores(info, env.cpuct, ϵ, η)
      action_id = argmax(scores)
      action = actions[action_id]
      wp = GI.white_playing(game)
      GI.play!(game, action)
      wr = GI.white_reward(game)
      r = wp ? wr : -wr
      pswitch = wp != GI.white_playing(game)
      increment_visit_counter!(env, state, action_id)
      qnext = worker_explore!(env, worker, game, η, false)
      qnext = pswitch ? -qnext : qnext
      q = r + env.gamma * qnext
      update_state_info!(env, state, action_id, q)
      env.total_nodes_traversed += 1
      return q
    end
  end
end

function inference_server(env::Env{G, S, A}) where {G, S, A}
  to_watch = env.workers
  while true
    # Update the list of workers to watch (the ones that are not done)
    requests = [take!(w.send) for w in to_watch]
    done = [isa(r, Done) for r in requests]
    active = .~ done
    any(active) || break
    to_watch = to_watch[active]
    requests = requests[active]
    # Gather a batch of queries and evaluate them
    batch = [q.state for q in requests if isa(q, Query)]
    if isempty(batch)
      answers, time = EvaluationResult{Float32}[], 0.
    else
      if env.fill_batches
        nmissing = length(env.workers) - length(batch)
        if nmissing > 0
          append!(batch, [batch[1] for i in 1:nmissing])
        end
        @assert length(batch) == length(env.workers)
      end
      answers, time = @timed evaluate_batch(env.oracle, batch)
    end
    env.inference_time += time
    # Send responses to all workers. Workers that sent `NoQuery` receive
    # a dummy answer, which they ignore.
    dummy_answer = (Float32[], 0f0)
    for (i, q) in enumerate(requests)
      if isa(q, Query)
        put!(to_watch[i].recv, popfirst!(answers))
      else
        put!(to_watch[i].recv, dummy_answer)
      end
    end
  end
end

function dirichlet_noise(game, α)
  actions = GI.available_actions(game)
  n = length(actions)
  return rand(Dirichlet(n, α))
end

function explore_sync!(env::Env, game, nsims)
  η = dirichlet_noise(game, env.noise_α)
  elapsed = @elapsed for i in 1:nsims
    worker_explore!(env, env.workers[1], game, η)
  end
  env.total_time += elapsed
end

# Spawn an inference server along with `env.nworkers` workers.
# Each worker runs simulations in a loop until `env.remaining` reaches zero.
function explore_async!(env::Env, game, nsims)
  env.remaining = nsims
  η = dirichlet_noise(game, env.noise_α)
  elapsed = @elapsed begin
    @sync begin
      @async @printing_errors inference_server(env)
      for w in env.workers
        @async @printing_errors begin
          lock(env.global_lock)
          while env.remaining > 0
            env.remaining -= 1
            worker_explore!(env, w, game, η)
          end
          put!(w.send, Done(env))
          unlock(env.global_lock)
        end
      end
    end
  end
  env.total_time += elapsed
  return
end

"""
    MCTS.explore!(env, game, nsims)

Run `nsims` MCTS simulations from the current state.
"""
function explore!(env::Env, game, nsims)
  asynchronous(env) ?
    explore_async!(env, game, nsims) :
    explore_sync!(env, game, nsims)
end

"""
    MCTS.policy(env, game)

Return the recommended stochastic policy on the current state.

A call to this function must always be preceded by
a call to [`MCTS.explore!`](@ref).
"""
function policy(env::Env, game)
  actions = GI.available_actions(game)
  state = GI.current_state(game)
  info =
    try env.tree[state]
    catch e
      if isa(e, KeyError)
        error("MCTS.explore! must be called before MCTS.policy")
      else
        rethrow(e)
      end
    end
  Ntot = sum(a.N for a in info.stats)
  π = [a.N / Ntot for a in info.stats]
  π ./= sum(π)
  return actions, π
end

"""
    MCTS.reset!(env)

Empty the MCTS tree.
"""
function reset!(env)
  empty!(env.tree)
  GC.gc(true)
end

#####
##### Profiling Utilities
#####

"""
    MCTS.inference_time_ratio(env)

Return the ratio of time spent by [`MCTS.explore!`](@ref)
on position evaluation (through functions [`MCTS.evaluate`](@ref) or
[`MCTS.evaluate_batch`](@ref)) since the environment's creation.
"""
function inference_time_ratio(env)
  T = env.total_time
  iszero(T) ? 0. : env.inference_time / T
end

"""
    MCTS.average_exploration_depth(env)

Return the average number of nodes that are traversed during an
MCTS simulation, not counting the root.
"""
function average_exploration_depth(env)
  return env.total_nodes_traversed / env.total_simulations
end

"""
    MCTS.memory_footprint_per_node(env)

Return an estimate of the memory footprint of a single node
of the MCTS tree (in bytes).
"""
function memory_footprint_per_node(env::Env{G}) where G
  # The hashtable is at most twice the number of stored elements
  # For every element, a state and a pointer are stored
  size_key = 2 * (GI.state_memsize(G) + sizeof(Int))
  dummy_stats = StateInfo([
    ActionStats(0, 0, 0, 0) for i in 1:GI.num_actions(G)], 0)
  size_stats = Base.summarysize(dummy_stats)
  return size_key + size_stats
end

"""
    MCTS.approximate_memory_footprint(env)

Return an estimate of the memory footprint of the MCTS tree (in bytes).
"""
function approximate_memory_footprint(env::Env)
  return memory_footprint_per_node(env) * length(env.tree)
end

# Possibly very slow for large trees
memory_footprint(env::Env) = Base.summarysize(env.tree)

end


text:
########## DIRECTORIES ##########

"""
gives the .mrbl dir
"""
function mrbldir()
    if "MARBLE_HOME" ∈ keys(ENV) # Would rather use `get_key`, but it doesnt work with ENV
        return "$(ENV["MARBLE_HOME"])"
    else
        return "$(homedir())/.mrbl"
    end
end
mrbldir(path...) = joinpath(mrbldir(), path...)

"""
Runs code snippets within the specified directory.
Due to the changeing directories, beware of relatvie paths
"""
function runindir(f::Function, path)
    curdir = pwd()
    cd(path)
    output = nothing
    try
        output = f()
    finally
        cd(curdir)
    end
    return output
end

########## PATH MANIPULATION  ##########

"""
Removes the file extension from a file or directory
"""
trimext(path) = splitext(path)[1]
trimext(doc::MarbleDoc) = trimext(doc.docname)

"""
Takes a given filename and returns its extension
Trims the leading `.`
"""
ext(path) = splitext(path)[2][2:end]

"""
Create a basename from a given path (filename without extension)
"""
get_basename(path) = splitext(splitdir(path)[2])[1]
get_basename(e::MarbleDoc) = get_basename(e.docname)

# Also use `basename`
# Maybe `dirname`
# splitext, splitdir

"""
Returns whether a document is a Markdown document
Should the default extension be .md or .mrbl
"""
ismarkdown(path, extension=["md"]) = ext(path) ∈ extension ? true : false
ismarkdown(path, s::SettingsBundle) = ismarkdown(path, extension=s["md_extension"])


text:
abstract type AbstractSequenceReplay <: AbstractReplay end

mutable struct EpisodicSequenceReplay{CB} <: AbstractSequenceReplay
    buffer::CB
    place::Int64
    terminal_locs::Vector{Int}
    terminal_symbol::Symbol
end

function EpisodicSequenceReplay(size, types, shapes, column_names; terminal_symbol = :t)
    cb = CircularBuffer(size, types, shapes, column_names)
    EpisodicSequenceReplay(cb, 1, Int[], terminal_symbol)
end

Base.length(er::EpisodicSequenceReplay) = length(er.buffer)
Base.getindex(er::EpisodicSequenceReplay, idx) =
    er.buffer[idx]

Base.view(er::EpisodicSequenceReplay, idx) =
    @view er.buffer[idx]

function Base.push!(er::EpisodicSequenceReplay, experience)
    if er.buffer._full
        er.place = (er.place % capacity(er.buffer)) + 1
    end
    push!(er.buffer, experience)
end

function get_episode_ends(er::EpisodicSequenceReplay)
    # TODO: n-computations. Maybe store in a cache?
    findall((exp)->exp::Bool, er.buffer._stg_tuple[er.terminal_symbol])
end

function get_valid_starting_range(s, e, seq_length)
    if e - seq_length <= s
        s:s
    else
        (s:e-seq_length)
    end
end

function get_valid_indicies(er::EpisodicSequenceReplay, seq_length)
    # episode_ends = get_episode_ends(er)

    1:(length(er) + 1 - seq_length)

end

function get_sequence(er::EpisodicSequenceReplay, start_ind, max_seq_length)
    # ret = [view(er, start_ind)]
    ret = [er[start_ind]]
    er_size = length(er)
    if ((start_ind + 1 - 1) % er_size) + 1 == er.place || ret[end][er.terminal_symbol][]::Bool
        return ret
    end
    for i ∈ 1:(max_seq_length-1)
        push!(ret, er[(((start_ind + i - 1) % er_size) + 1)])
        if ret[end][er.terminal_symbol][]::Bool || ((start_ind + i + 1 - 1) % er_size) + 1 == er.place
            break
        end
    end
    ret
end

sample(er::EpisodicSequenceReplay,
       batch_size,
       max_seq_length) =
           sample(Random.GLOBAL_RNG, er, batch_size, max_seq_length)

function sample(rng::Random.AbstractRNG,
                er::EpisodicSequenceReplay,
                batch_size,
                max_seq_length)
    # get valid starting indicies
    valid_inx = get_valid_indicies(er, 1)
    start_inx = rand(rng, valid_inx, batch_size)
    exp = [get_sequence(er, si, max_seq_length) for si ∈ start_inx]
    start_inx, exp
    # padding and batching handled by agent.
end