text:
# SOCP Lasso Testproblem

using COSMO, Test, LinearAlgebra, SparseArrays, Random

# This test is precision agnostic and will be run with T precision
if !@isdefined(UnitTestFloats)
    UnitTestFloats = [Float64] #if not run in full test setup, just do it for one float type
end

@testset "SOCP - Lasso" begin
  for T in UnitTestFloats
      if T != BigFloat
      # generate problem data
      rng = Random.MersenneTwister(12345)
      n = 8
      m = 50 * n
      F = rand(rng, T, m, n)

      vtrue = sprand(rng, T, n, 1, 0.1)
      noise = T(0.1) * rand(rng, T, m, 1)
      b = F * vtrue + noise
      μMax = norm(F' * b, Inf)
      μ = T(0.1) * μMax

      # define lasso problem as SOCP

      A1 = -sparse([1 zeros(T, 1, 2 * n + 1) 1 zeros(T, 1, m);
        -1 zeros(T, 1, 2 * n) 1 zeros(T, 1, m + 1);
        zeros(T, m, 1) -2 * F zeros(T, m, n + 2) diagm( 0 => ones(T, m))])

      A2 = -sparse([zeros(T, n, 1) diagm(0 => ones(T, n)) -diagm(0 => ones(T, n)) zeros(T, n, m + 2);
        zeros(T, n, 1) -diagm(0 => ones(T, n)) -diagm(0 => ones(T, n)) zeros(T, n, m + 2)])
      A3 = -sparse([zeros(T, 1, 2 * n + 1) -one(T) zeros(T, 1, m + 1);
       zeros(T, 1, 2 * n + 2) -one(T) zeros(T, 1, m);
       zeros(T, m, 2 * n + 3) -diagm( 0 => ones(T, m))])
      b1 = T[1; 1; -2 * b]
      b2 = zeros(T, 2 * n)
      b3 = zeros(T, m + 2)

      q = [one(T); zeros(T, n); μ * ones(T, n); zeros(T, m + 2)]
      P = spzeros(T, length(q), length(q))

      cs1 = COSMO.Constraint(A1 ,b1, COSMO.ZeroSet)
      cs2 = COSMO.Constraint(A2 ,b2, COSMO.Nonnegatives)
      cs3 = COSMO.Constraint(A3 ,b3, COSMO.SecondOrderCone)

      # Solve with OSSDP
      model = COSMO.Model{T}()
      assemble!(model, P, q, [cs1; cs2; cs3])
      res = COSMO.optimize!(model);

      @testset "SOCP - Lasso (T = $(T))" begin
      @test res.status == :Solved
      @test isapprox(res.obj_val, 0.4422849814458825,atol=1e-2)
      end
    end
  end
end
nothing


text:
using Cairo
include("shape_functions.jl")

function test_all(;tmax = 2.0, save_flag = false)
    print("cairo version: ",Cairo.libcairo_version,"\n");
    size_surface = [512];#,512,1024]; #three sizes of a surface
    paint_width = [0.5,1.0,3.0,5.0];
    shapes = [ddots1, ddots2, ddots3, ddots4, ddots5,
              rdots1, rdots2, rdots3, rdots4, rdots5, rdots6,
              lines0, lines1, lines2, lines3, lines4];
    n_elements = [100,300,1000,3000,10000,30000,100000];

    n_tests = length(size_surface) * length(paint_width) * length(shapes);
    test_data = zeros(n_tests+1,3 + length(n_elements));
    test_data[1,3+(1:length(n_elements))] = n_elements;
    i = 1;

    if save_flag && !isdir("test_results")
        mkdir("test_results")
    end

    for s_size in size_surface
        s = Cairo.CairoARGBSurface(s_size,s_size);
        c = Cairo.CairoContext(s);
        print("Surface Size: ",s_size,"\n")
        for w in paint_width
            print("Paint Width: ",w,"\n")
            for m_index in 1:length(shapes)
                test_data[i,1] = s_size;
                test_data[i,2] = w;
                test_data[i,3] = m_index;
                m = shapes[m_index];
                for n_index in 1:length(n_elements)
                    n = n_elements[n_index];
                    t = @elapsed(m(c,s_size,s_size,w,n));
                    print(m," ",@sprintf("%6d",n),"  elapsed time: ",@sprintf("%f seconds",t),"\n");
                    test_data[i+1,3+n_index] = t;
                    if save_flag
                        write_to_png(s,@sprintf("test_results/test_speed_%2.1f_%02d_%02d.png",w,m_index,n_index));
                    end
                    if t > tmax && n_index < length(n_elements)
                        println("Remaining tests in this sequence skipped")
                        break
                    end
                end
                print("\n")
                i += 1;
            end
        end
    end
end

test_all()

ext: jl


text:
@testset "bayesquad" begin
    k = SqExponentialKernel()

    σ = 4.0
    @test BQ.scale(BayesQuad(σ * k)) == σ

    test_convergence(BayesQuad(SqExponentialKernel()), PriorSampling())
end


text:
include("limbdark/src/transit_poly_struct.jl")

# produces an array of impact parameters for a transit assuming constant velocity
function b_array(t, t0, b0, rp, duration)
    b = zeros(Real, length(t))
    chordlength = 2*sqrt((1+rp)^2-b0^2)
    v = chordlength/duration

    for i in 1:length(t)
        if abs(t[i]-t0) > chordlength/(2*v)
            b[i] = 1+rp
        else
            x = v*(t[i]-t0)
            b[i] = sqrt(b0^2 + x^2)
        end
    end
    return b
end

# transit light curve at times t (center=t0, maximum impact parameter = b0, radius=rp, duration=d, limb darkening parameters = u)
function transit(t, t0, b0, rp, d, u)
    b = b_array(t, t0, b0, rp, d)
    trans = transit_init(rp, b0, u, true)
    flux = zeros(Real, length(b))
    gradflux = zeros(Real, length(b), 2 + length(u))
    for i in 1:length(b)
        trans.b = b[i]
        flux[i] = transit_poly!(trans)
        gradflux[i,:] = trans.dfdrbu
    end
    return flux, gradflux
end


text:
@testset "Cylinder" begin

C = GRAYTR.Cylinder(scaling(1, 1, 0.5))
@test C.id == 1
@test GRAYTR.can_intersect(C)
@test GRAYTR.area(C) ≈ 4π

@test GRAYTR.obj_bounds(C) == GRAYTR.BoundingBox(Point3(-1, -1, -1), Point3(1, 1, 1))
@test GRAYTR.world_bounds(C) == GRAYTR.BoundingBox(Point3(-1, -1, -0.5), Point3(1, 1, 0.5))

T = translation(0, 0, 0.5)
C2 = T(C)
@test GRAYTR.world_bounds(C2) == GRAYTR.BoundingBox(Point3(-1, -1, 0), Point3(1, 1, 1))

R = GRAYTR.Ray(Point3(5, 0, 0), Vector3(-1, 0, 0))
@test GRAYTR.intersectP(R, C)

R = GRAYTR.Ray(Point3(5, 4, 0), Vector3(-1, 0, 0))
@test !GRAYTR.intersectP(R, C)

R = GRAYTR.Ray(Point3(0, 0, 3), Vector3(0, 0, -1))
@test !GRAYTR.intersectP(R, C)

R = GRAYTR.Ray(Point3(5, 0, 0), Vector3(-1, 0, 0), 15.0, Inf, 1)
@test !GRAYTR.intersectP(R, C)

R = GRAYTR.Ray(Point3(5, 0, 0), Vector3(-1, 0, 0))
dg, t, e = GRAYTR.shape_intersect(R, C)
@test dg.p ≈ Point3(1, 0, 0)
@test t ≈ 4

R = GRAYTR.Ray(Point3(5, 4, 0), Vector3(-1, 0, 0))
dg, t, e = GRAYTR.shape_intersect(R, C)
@test dg == nothing
@test isnan(t)
@test isnan(e)

R = GRAYTR.Ray(Point3(0, 0, 3), Vector3(0, 0, -1))
dg, t, e = GRAYTR.shape_intersect(R, C)
@test dg == nothing
@test isnan(t)
@test isnan(e)

R = GRAYTR.Ray(Point3(5, 0, 0), Vector3(-1, 0, 0), 15.0, Inf, 1)
dg, t, e = GRAYTR.shape_intersect(R, C)
@test dg == nothing
@test isnan(t)
@test isnan(e)

end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
module test_julia_cfunction

using Test

f() = 2
g = @cfunction f Int ()
@test ccall(g, Int, ()) == 2

record = Any[]
function cb_func(delta::Cdouble, msg::Ptr{Cuchar}, cb_info::Ptr{Cvoid})::Cvoid
    push!(record, (delta, msg, cb_info))
    nothing
end
h = @cfunction cb_func Cvoid (Cdouble, Ptr{Cuchar}, Ptr{Cvoid})

msg = pointer("abc")
cb_info = C_NULL
ccall(h, Cvoid, (Cdouble, Ptr{Cuchar}, Ptr{Cvoid}), Cdouble(1), msg, cb_info)

@test record == [(1.0, msg, cb_info)]

end # module test_julia_cfunction


text:
# ---
# title: 136. Single Number
# id: problem136
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Hash Table, Bit Manipulation
# link: <https://leetcode.com/problems/single-number/description/>
# hidden: true
# ---
#
# Given a **non-empty**  array of integers `nums`, every element appears _twice_
# except for one. Find that single one.
#
# **Follow up:**  Could you implement a solution with a linear runtime
# complexity and without using extra memory?
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [2,2,1]
#     Output: 1
#
#
# **Example 2:**
#
#
#
#     Input: nums = [4,1,2,1,2]
#     Output: 4
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 3 * 104`
#   * `-3 * 104 <= nums[i] <= 3 * 104`
#   * Each element in the array appears twice except for one element which appears only once.
#
#
## @lc code=start
using LeetCode

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


text:
module SBML

using SBML_jll, Libdl
using SparseArrays
using Symbolics
using IfElse
using Unitful

include("types.jl")
include("structs.jl")
include("version.jl")

include("converters.jl")
include("math.jl")
include("readsbml.jl")
include("symbolics.jl")
include("utils.jl")

sbml(sym::Symbol) = dlsym(SBML_jll.libsbml_handle, sym)

export readSBML, getS, getLBs, getUBs, getOCs
export set_level_and_version, libsbml_convert, convert_simplify_math

end # module


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -0.27492067550218746 <= q <= 2.084201684402027)

                     add_NL_constraint(m, :(swish(-0.0416569430659357 + 0.13675627334278007*swish(0.1752604319892206 + -0.8557289632032381*swish(-0.028885055816546323 + 0.15091441453709065*swish(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + -0.8760277243880852*swish(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))) + 0.21149394336428795*swish(-0.29077418744555983 + -0.11119024401630329*swish(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + 0.4817613487326513*swish(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3])))) + 0.703085236566185*swish(-0.4914467550189556 + 0.5078124686195609*swish(-0.028885055816546323 + 0.15091441453709065*swish(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + -0.8760277243880852*swish(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))) + 0.7523099505151869*swish(-0.29077418744555983 + -0.11119024401630329*swish(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + 0.4817613487326513*swish(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))) + swish(0.6991595157542867 + 0.5132010992296534*swish(0.1752604319892206 + -0.8557289632032381*swish(-0.028885055816546323 + 0.15091441453709065*swish(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + -0.8760277243880852*swish(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))) + 0.21149394336428795*swish(-0.29077418744555983 + -0.11119024401630329*swish(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + 0.4817613487326513*swish(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3])))) + 0.6682334830395109*swish(-0.4914467550189556 + 0.5078124686195609*swish(-0.028885055816546323 + 0.15091441453709065*swish(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + -0.8760277243880852*swish(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))) + 0.7523099505151869*swish(-0.29077418744555983 + -0.11119024401630329*swish(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + 0.4817613487326513*swish(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
######################################################
# Jacobian of `f(::AbstractArray...)::AbstractArray` #
######################################################

"""
    ReverseDiff.jacobian(f, input, cfg::JacobianConfig = JacobianConfig(input))

If `input` is an `AbstractArray`, assume `f` has the form
`f(::AbstractArray{<:Real})::AbstractArray{<:Real}` and return `J(f)(input)`.

If `input` is a tuple of `AbstractArray`s, assume `f` has the form
`f(::AbstractArray{<:Real}...)::AbstractArray{<:Real}` (such that it can be called as
`f(input...)`) and return a `Tuple` where the `i`th element is the  Jacobian of `f` w.r.t.
`input[i].`

Note that `cfg` can be preallocated and reused for subsequent calls.

If possible, it is highly recommended to use `ReverseDiff.JacobianTape` to prerecord `f`.
Otherwise, this method will have to re-record `f`'s execution trace for every subsequent
call.
"""
function jacobian(f, input, cfg::JacobianConfig = JacobianConfig(input))
    tape = JacobianTape(f, input, cfg)
    isa(input, TrackedArray) && empty!(input.tape)
    result = jacobian!(tape, input)
    empty!(tape.tape)
    return result
end

"""
    ReverseDiff.jacobian!(result, f, input, cfg::JacobianConfig = JacobianConfig(input))

Returns `result`. This method is exactly like `ReverseDiff.jacobian(f, input, cfg)`, except
it stores the resulting Jacobian(s) in `result` rather than allocating new memory.

`result` can be an `AbstractArray` or a `Tuple` of `AbstractArray`s. The `result` (or any
of its elements, if `isa(result, Tuple)`), can also be a `DiffResults.DiffResult`, in which
case the primal value `f(input)` (or `f(input...)`, if `isa(input, Tuple)`) will be stored
in it as well.
"""
function jacobian!(result, f, input, cfg::JacobianConfig = JacobianConfig(input))
    tape = JacobianTape(f, input, cfg)
    isa(input, TrackedArray) && empty!(input.tape)
    jacobian!(result, tape, input)
    empty!(tape.tape)
    return result
end

#########################################################
# Jacobian of `f!(::AbstractArray, ::AbstractArray...)` #
#########################################################

"""
    ReverseDiff.jacobian(f!, output, input, cfg::JacobianConfig = JacobianConfig(output, input))

Exactly like `ReverseDiff.jacobian(f, input, cfg)`, except the target function has the
form `f!(output::AbstractArray{<:Real}, input::AbstractArray{<:Real}...)`.
"""
function jacobian(f!, output, input, cfg::JacobianConfig = JacobianConfig(output, input))
    tape = JacobianTape(f!, output, input, cfg)
    isa(input, TrackedArray) && empty!(input.tape)
    result = jacobian!(tape, input)
    extract_result_value!(output, output_hook(tape))
    empty!(tape.tape)
    return result
end

"""
    ReverseDiff.jacobian!(result, f!, output, input, cfg::JacobianConfig = JacobianConfig(output, input))

Exactly like `ReverseDiff.jacobian!(result, f, input, cfg)`, except the target function has the
form `f!(output::AbstractArray{<:Real}, input::AbstractArray{<:Real}...)`.
"""
function jacobian!(result, f!, output, input, cfg::JacobianConfig = JacobianConfig(output, input))
    tape = JacobianTape(f!, output, input, cfg)
    isa(input, TrackedArray) && empty!(input.tape)
    jacobian!(result, tape, input)
    extract_result_value!(output, output_hook(tape))
    empty!(tape.tape)
    return result
end

###########################
# Executing JacobianTapes #
###########################

"""
    ReverseDiff.jacobian!(tape::Union{JacobianTape,CompiledJacobian}, input)

If `input` is an `AbstractArray`, assume `tape` represents a function of the form
`f(::AbstractArray{<:Real})::AbstractArray{<:Real}` or `f!(::AbstractArray{<:Real},
::AbstractArray{<:Real})` and return `tape`'s Jacobian w.r.t. `input`.

If `input` is a tuple of `AbstractArray`s, assume `tape` represents a function of the form
`f(::AbstractArray{<:Real}...)::AbstractArray{<:Real}` or `f!(::AbstractArray{<:Real},
::AbstractArray{<:Real}...)` and return a `Tuple` where the `i`th element is `tape`'s Jacobian
w.r.t. `input[i].`

Note that if `tape` represents a function of the form `f!(output, input...)`, you can only
execute `tape` with new `input` values. There is no way to re-run `tape`'s tape with new
`output` values; since `f!` can mutate `output`, there exists no stable "hook" for loading
new `output` values into the tape.
"""
function jacobian!(tape::Union{JacobianTape,CompiledJacobian}, input)
    result = construct_result(output_hook(tape), input_hook(tape))
    jacobian!(result, tape, input)
    return result
end

"""
    ReverseDiff.jacobian!(result, tape::Union{JacobianTape,CompiledJacobian}, input)

Returns `result`. This method is exactly like `ReverseDiff.jacobian!(tape, input)`, except it
stores the resulting Jacobian(s) in `result` rather than allocating new memory.

`result` can be an `AbstractArray` or a `Tuple` of `AbstractArray`s. The `result` (or any
of its elements, if `isa(result, Tuple)`), can also be a `DiffResults.DiffResult`, in which
case the primal value of the target function will be stored in it as well.
"""
function jacobian!(result, tape::Union{JacobianTape,CompiledJacobian}, input)
    seeded_forward_pass!(tape, input)
    seeded_reverse_pass!(result, tape)
    return result
end

##################################################
# unused (but faster) versions of the above code #
##################################################

#=
These commented-out versions of `jacobian` are faster than the ones we're
actually using above, because they avoid a redundant forward pass. This extra
forward pass should be unneccesary - since no input values are changing,
the record pass should be sufficient on its own. However, for some unknown
reason, getting rid of the superfluous forward pass breaks nested
differentation.
=#

# function jacobian(f, input, cfg::JacobianConfig = JacobianConfig(input))
#     tape = JacobianTape(f, input, cfg)
#     result = construct_result(output_hook(tape), input_hook(tape))
#     seeded_reverse_pass!(result, output_hook(tape), input_hook(tape), tape.tape)
#     empty!(tape.tape)
#     return result
# end
#
# function jacobian!(result, f, input, cfg::JacobianConfig = JacobianConfig(input))
#     tape = JacobianTape(f, input, cfg)
#     seeded_reverse_pass!(result, output_hook(tape), input_hook(tape), tape.tape)
#     empty!(tape.tape)
#     return result
# end
#
# function jacobian(f!, output, input, cfg::JacobianConfig = JacobianConfig(output, input))
#     tape = JacobianTape(f!, output, input, cfg)
#     result = construct_result(output_hook(tape), input_hook(tape))
#     seeded_reverse_pass!(result, output_hook(tape), input_hook(tape), tape.tape)
#     extract_result_value!(output, output_hook(tape))
#     empty!(tape.tape)
#     return result
# end
#
# function jacobian!(result, f!, output, input, cfg::JacobianConfig = JacobianConfig(output, input))
#     tape = JacobianTape(f!, output, input, cfg)
#     seeded_reverse_pass!(result, output_hook(tape), input_hook(tape), tape.tape)
#     extract_result_value!(output, output_hook(tape))
#     empty!(tape.tape)
#     return result
# end


text:
export EpsilonGreedySelector, get_prob

using StatsBase: sample
using .Utils: findallmax

"""
    EpsilonGreedySelector{T}(;kwargs...)

> Epsilon-greedy strategy: The best lever is selected for a proportion `1 - epsilon` of the trials, and a lever is selected at random (with uniform probability) for a proportion epsilon . [Multi-armed_bandit](https://en.wikipedia.org/wiki/Multi-armed_bandit)

Two kinds of epsilon-decreasing strategy are implmented here (`linear` and `exp`).

> Epsilon-decreasing strategy: Similar to the epsilon-greedy strategy, except that the value of epsilon decreases as the experiment progresses, resulting in highly explorative behaviour at the start and highly exploitative behaviour at the finish.  - [Multi-armed_bandit](https://en.wikipedia.org/wiki/Multi-armed_bandit)

# Keywords

- `T::Symbol`: defines how to calculate the epsilon in the warmup steps. Supported values are `linear` and `exp`.
- `step::Int = 1`: record the current step.
- `ϵ_init::Float64 = 1.0`: initial epsilon.
- `warmup_steps::Int=0`: the number of steps to use `ϵ_init`.
- `decay_steps::Int=0`: the number of steps for epsilon to decay from `ϵ_init` to `ϵ_stable`.
- `ϵ_stable::Float64`: the epsilon after `warmup_steps + decay_steps`.

# Example

```julia
s = EpsilonGreedySelector{:linear}(ϵ_init=0.9, ϵ_stable=0.1, warmup_steps=100, decay_steps=100)
plot([RL.get_ϵ(s, i) for i in 1:500], label="linear epsilon")
```
![](../assets/img/linear_epsilon_greedy_selector.png)

```julia
s = EpsilonGreedySelector{:exp}(ϵ_init=0.9, ϵ_stable=0.1, warmup_steps=100, decay_steps=100)
plot([RL.get_ϵ(s, i) for i in 1:500], label="exp epsilon")
```
![](../assets/img/exp_epsilon_greedy_selector.png)
"""
Base.@kwdef mutable struct EpsilonGreedySelector{T} <: AbstractDiscreteActionSelector
    ϵ_stable::Float64
    ϵ_init::Float64 = 1.0
    warmup_steps::Int = 0
    decay_steps::Int = 0
    step::Int = 1
end

"""
    EpsilonGreedySelector(ϵ) -> EpsilonGreedySelector{:linear}(; ϵ_stable = ϵ)
"""
EpsilonGreedySelector(ϵ) = EpsilonGreedySelector{:linear}(; ϵ_stable = ϵ)

function get_ϵ(s::EpsilonGreedySelector{:linear}, step)
    if step <= s.warmup_steps
        s.ϵ_init
    elseif step >= (s.warmup_steps + s.decay_steps)
        s.ϵ_stable
    else
        steps_left = s.warmup_steps + s.decay_steps - step
        s.ϵ_stable + steps_left / s.decay_steps * (s.ϵ_init - s.ϵ_stable)
    end
end

function get_ϵ(s::EpsilonGreedySelector{:exp}, step)
    if step <= s.warmup_steps
        s.ϵ_init
    else
        n = step - s.warmup_steps
        s.ϵ_stable + (s.ϵ_init - s.ϵ_stable) * exp(-1.0 * n / s.decay_steps)
    end
end

get_ϵ(s::EpsilonGreedySelector) = get_ϵ(s, s.step)

"""
    (s::EpsilonGreedySelector)(values; step) where T

!!! note
    If multiple values with the same maximum value are found.
    Then a random one will be returned!

    `NaN` will be filtered unless all the values are `NaN`.
    In that case, a random one will be returned.
"""
function (s::EpsilonGreedySelector)(values)
    ϵ = get_ϵ(s)
    s.step += 1
    rand() > ϵ ? sample(findallmax(values)[2]) : rand(1:length(values))
end

"""
    get_prob(s::EpsilonGreedySelector, values)

Return the probability of selecting each action given the estimated `values` of each action.
"""
function get_prob(s::EpsilonGreedySelector, values)
    ϵ, n = get_ϵ(s), length(values)
    probs = fill(ϵ / n, n)
    max_val_inds = findallmax(values)[2]
    for ind in max_val_inds
        probs[ind] += (1 - ϵ) / length(max_val_inds)
    end
    probs
end


text:
module microcoverage

## Utility function to construct an Expr from head and args.

function mkex(head::Symbol, args...)
    v = Expr(head)
    for arg in args
        push!(v.args, arg)
    end
    v
end

## We annotate the source file with a sequence of ints and strings.

typealias AnnotationIndexType Array{Union(Int,ASCIIString), 1}

## Constructor for an annotation index.

initannotationindex() = (Union(Int,ASCIIString))[]

# This dictionary maps source file names to their
# annotation indices.

const annotationindexdict = Dict{ASCIIString, AnnotationIndexType}()

# Trackarray holds the counters that track the line
# numbers and calls.

const trackarray = (Int)[]

## incrtrackarray increments the appropriate entry of the
## track array

function incrtrackarray(subscr)
    global trackarray::Array{Int,1}
    trackarray[subscr] += 1
    nothing
end

## makecallinst: Generate an expression that corresponds to the
## statement
##    Main.microcoverage.incrtrackarray(trknum)

makecallinst(trknum::Int) =
mkex(:call,
     mkex(:.,
          mkex(:.,
                    :Main,
               QuoteNode(:microcoverage)),
          QuoteNode(:incrtrackarray)),
     trknum)

## coverage_rewrite_recursive
## Take an expression inprogram (could be a whole program),
## the name of the tracking array, and an initial position
## in the tracking array.
##
## Produce an outprogram, which is the rewritten inprogram
## with calls to the coverage-statement incrementers.
## Also return the updated trknum (position in the tracking array).
## As a side effect, generate annotationindex, which tracks what should
## be printed in the output after coverage checking is complete.
##
## There are three versions of routine; the correct one is
## selected by multiple dispatch on the first argument.

function coverage_rewrite_recursive(inprogram::LineNumberNode,
                                    starttrknum::Int,
                                    linenumberoffset::Int,
                                    annotationindex::AnnotationIndexType)

    # If the argument is of type LineNumberNode, i.e.,
    # a new line number in the source code,
    # replace it with a block that consists of the
    # the line number and a tracking array incrementer
    trknum = starttrknum
    newline = inprogram.line + linenumberoffset
    push!(annotationindex, -newline, trknum)
    outprogram = mkex(:block,
                      deepcopy(inprogram),
                      makecallinst(trknum))
    trknum += 1
    return outprogram, trknum
end

function coverage_rewrite_recursive(inprogram::Any,
                                    starttrknum::Int,
                                    linenumberoffset::Int,
                                    ::AnnotationIndexType)

    ## This is the default version of coverage_rewrite_recursive
    ## that doesn't do anything

    return inprogram, starttrknum
end

function coverage_rewrite_recursive(inprogram::Expr,
                                    starttrknum::Int,
                                    linenumberoffset::Int,
                                    annotationindex::AnnotationIndexType)
    ## This is the primary version of coverage_rewrite
    ## recursive.  It takes an expression and inserts
    ## tracking statements for line numbers and internal branches
    ## in expressions.
    trknum = starttrknum
    if inprogram.head == :line
        # If the expression is an expression of type :line, i.e.,
        # a new line number in the source code,
        # replace it with a block that consists of the
        # the line number and a tracking array incrementer
        newline = inprogram.args[1] + linenumberoffset
        push!(annotationindex, -newline, trknum)
        outprogram = mkex(:block,
                          deepcopy(inprogram),
                          makecallinst(trknum))
        trknum += 1
    elseif inprogram.head == :if && (!(typeof(inprogram.args[2]) <: Expr) ||
        inprogram.args[2].head != :block)
        # If the expression is of the form a? b : c
        # then generate tracking statements for b and c
        outprogram = Expr(:if)
        outprogram1, trknum = coverage_rewrite_recursive(inprogram.args[1],
                                                         trknum,
                                                         linenumberoffset,
                                                         annotationindex)
        push!(outprogram.args, outprogram1)
        push!(annotationindex, " ? ")
        @assert(length(inprogram.args) == 3)
        for k = 2 : 3
            if k > 2
                push!(annotationindex, " : ")
            end
            a2 = inprogram.args[k]
            push!(annotationindex, "(", trknum)
            callinst = makecallinst(trknum)
            trknum += 1
            outprogram1, trknum = coverage_rewrite_recursive(a2,
                                                             trknum,
                                                             linenumberoffset,
                                                             annotationindex)
            push!(outprogram.args, mkex(:block,
                                        callinst,
                                        outprogram1))
            push!(annotationindex, ")")
        end
    elseif inprogram.head == :(=) && typeof(inprogram.args[1]) <: Expr &&
        inprogram.args[1].head == :call && inprogram.args[1].args[1] != :eval
        # If the line is a statement-function definition other than the
        # definition of "eval", then insert
        # a tracking statement into the function body.
        @assert length(inprogram.args) == 2
        outprogram1, trknum = coverage_rewrite_recursive(inprogram.args[1],
                                                         trknum,
                                                         linenumberoffset,
                                                         annotationindex)
        savetrknum = trknum
        callinst = makecallinst(trknum)
        trknum += 1
        outprogram2, trknum = coverage_rewrite_recursive(inprogram.args[2],
                                                         trknum,
                                                         linenumberoffset,
                                                         annotationindex)
        push!(annotationindex, "(called", savetrknum, "time(s))")
        outprogram = mkex(:(=), outprogram1, mkex(:block,
                                                  callinst,
                                                  outprogram2))
    elseif inprogram.head == :|| || inprogram.head == :&&
        ## If the expression is the || or && operator, generate
        ## a tracking statement for each branch.
        @assert length(inprogram.args) == 2
        outprogram = Expr(inprogram.head)
        for k = 1 : 2
            callinst = makecallinst(trknum)
            push!(annotationindex, "(", trknum)
            trknum += 1
            outprogram1, trknum = coverage_rewrite_recursive(inprogram.args[k],
                                                             trknum,
                                                             linenumberoffset,
                                                             annotationindex)
            push!(outprogram.args, mkex(:block,
                                        callinst,
                                        outprogram1))
            push!(annotationindex, ")")
            if k == 1
                if inprogram.head == :||
                    push!(annotationindex, " || ")
                else
                    push!(annotationindex, " && ")
                end
            end
        end

    elseif inprogram.head == :global || inprogram.head == :import ||
        inprogram.head == :importall || inprogram.head == :export ||
        inprogram.head == :typealias || inprogram.head == :abstract ||
        inprogram.head == :using

        outprogram = inprogram
    elseif inprogram.head == :immutable || inprogram.head == :type
        outprogram = Expr(inprogram.head)
        for expr1 in inprogram.args
            if typeof(expr1) <: Expr &&
                (expr1.head == :(=) || expr1.head == :function)
                outprogram1, trknum = coverage_rewrite_recursive(expr1,
                                                                 trknum,
                                                                 linenumberoffset,
                                                                 annotationindex)
                push!(outprogram.args, outprogram1)
            else
                push!(outprogram.args, expr1)
            end
        end
    else
        ## For all other expression types, just make the output same as
        ## the input (with recursive calls)
        outprogram = Expr(inprogram.head)
        for expr1 in inprogram.args
            outprogram1, trknum = coverage_rewrite_recursive(expr1,
                                                             trknum,
                                                             linenumberoffset,
                                                             annotationindex)
            push!(outprogram.args, outprogram1)
        end
    end
    outprogram, trknum
end

function linecount(string)
    count = 0
    pos = 1
    while true
        a = search(string, '\n', pos)
        a == 0 && break
        count += 1
        pos = a + 1
    end
    count
end

filepreamble = "# Automatically generated by microcoverage.jl-- will be automatically deleted upon completion"

## This function takes a sourcefile name.  It renames
## it to <oldname>.orig.  It parses the original file
## and inserts tracking statements.
## Then it generates a new
## sourcefile with the same name as the old.  The new
## file eval's the parsed version of the old file with
## tracking statements.

function begintrack(sourcefilename::ASCIIString)
    println("reading $sourcefilename")
    src = ""
    open(sourcefilename, "r") do h
        src = convert(ASCIIString, readbytes(h))
    end
    annotationindex = initannotationindex()
    global trackarray::Array{Int,1}
    lasttrknum = length(trackarray)
    initsize = lasttrknum
    srcpos = 1
    src_parse_rewrite = (Expr)[]
    println("parsing")
    linenumberoffset = 0
    while srcpos <= length(src)
        if isspace(src[srcpos])
            if src[srcpos] == '\n'
                linenumberoffset += 1
            end
            srcpos += 1
        elseif src[srcpos] == '#'
            eolpos = search(src, '\n', srcpos)
            srcpos = eolpos + 1
            linenumberoffset += 1
        else
            src_parse, srcposfinal = parse(src, srcpos)
            rewrite1,lasttrknum = coverage_rewrite_recursive(src_parse,
                                                             lasttrknum + 1,
                                                             linenumberoffset,
                                                             annotationindex)
            linenumberoffset += linecount(src[srcpos : srcposfinal - 1])
            srcpos = srcposfinal
            push!(src_parse_rewrite, rewrite1)
        end
    end
    resize!(trackarray, lasttrknum)
    for j = initsize + 1 : lasttrknum
        trackarray[j] = 0
    end
    global annotationindexdict::Dict{ASCIIString,AnnotationIndexType}
    annotationindexdict[sourcefilename] = annotationindex
    renamed = sourcefilename * ".orig"
    if stat(renamed).size > 0
        error("Cannot rename original -- file already exists with the name $renamed")
    end
    println("renaming $sourcefilename to $renamed")
    mv(sourcefilename, renamed)
    println("saving machine-generated code in $sourcefilename")
    open(sourcefilename,"w") do h2
        global filepreamble::ASCIIString
        println(h2, filepreamble)
        for rewrite in src_parse_rewrite
            ss = IOBuffer()
            serialize(ss, rewrite)
            ser_rewrite = takebuf_array(ss)
            numbyte = length(ser_rewrite)
            println(h2, "eval(deserialize(IOBuffer((UInt8)[")
            for count = 1 : numbyte
                byte = ser_rewrite[count]
                show(h2, byte)
                count < numbyte && print(h2, ", ")
                count % 8 == 0 && println(h2)
            end
            println(h2, "])))")
        end
    end
end

# Small routine to make an ASCIIString consisting of i spaces
spaces(i::Int) = convert(ASCIIString, 32*ones(UInt8,i))

## The next four functions are handlers for items in annotationindex.
## This first one handles an integer-- dispatches on whether it is positive
## or negative

function printmcov(item::Int,
                   lastprint::Int,
                   curcol::Int,
                   newline::Int,
                   horig::IO,
                   hcov::IO)
    if item < 0
        printmcovli(-item, lastprint,curcol,newline,horig,hcov)
    else
        printmcovtn(item, lastprint,curcol,newline,horig,hcov)
    end
end

## Prints a line number and source line and advances the file
## to that line number, printing more source lines as necessary.

function printmcovli(lineno::Int,
                     lastprint::Int,
                     curcol::Int,
                     newline::Int,
                     horig::IO,
                     hcov::IO)
    oldnewline = newline
    newline = lineno
    if newline < lastprint
        error("Line numbers out of order in tracking info")
    end
    for count = lastprint + 1 : newline - 1
        s = chomp(readline(horig))
        if curcol > 16
            println(hcov)
            curcol = 0
        end
        println(hcov, spaces(16-curcol), "* ", s)
        curcol = 0
    end
    lastprint = newline - 1
    if newline != oldnewline
        nls = "$newline"
        print(hcov, "L", nls, spaces(8 - length(nls)))
        curcol += 9
    end
    lastprint, curcol, newline
end

## Print a tracking number

function printmcovtn(tracknum::Int,
                     lastprint::Int,
                     curcol::Int,
                     newline::Int,
                     ::IO,
                     hcov::IO)

    global trackarray::Array{Int,1}
    cov = trackarray[tracknum]
    cs = " $cov "
    print(hcov, cs)
    curcol += length(cs)
    lastprint, curcol, newline
end

## Print a string

function printmcov(outstring::ASCIIString,
                   lastprint::Int,
                   curcol::Int,
                   newline::Int,
                   ::IO,
                   hcov::IO)
    print(hcov, outstring)
    curcol += length(outstring)
    lastprint, curcol, newline
end

## endtrack is called when the trackign is finished
## It produces a coverage report in a file called <origfilename>.mcov.
## Then it renames the files back to how they were before
## begintrack was called.

function endtrack(sourcefilename::ASCIIString)
    renamed = sourcefilename * ".orig"
    covfilename = sourcefilename * ".mcov"
    open(renamed, "r") do horig
        open(covfilename, "w") do hcov
            println("Writing coverage information to $covfilename")
            global annotationindexdict::Dict{ASCIIString,AnnotationIndexType}
            lastprint = 0
            curcol = 0
            newline = -1
            for item in annotationindexdict[sourcefilename]
                lastprint, curcol, newline = printmcov(item,
                                                       lastprint,
                                                       curcol,
                                                       newline,
                                                       horig,
                                                       hcov)
            end
            if curcol > 16
                curcol = 0
                println(hcov)
            end
            while !eof(horig)
                s = chomp(readline(horig))
                println(hcov, spaces(16-curcol), "* ", s)
                curcol = 0
            end
        end
    end
    restore(sourcefilename)
end

function restore(sourcefilename::ASCIIString)
    renamed = sourcefilename * ".orig"
    if stat(renamed).size == 0
        error("No file called $renamed")
    end
    s = ""
    open(sourcefilename, "r") do hrewr
        s = chomp(readline(hrewr))
    end

    global filepreamble::ASCIIString
    if s != filepreamble
        error("Cannot overwrite $sourcefilename; missing preamble statement")
    end
    println("renaming $renamed to $sourcefilename; machine-generated $sourcefilename overwritten")
    mv(renamed, sourcefilename, remove_destination=true)
end

export restore
export begintrack
export endtrack

end

max_stars_repo_name: StephenVavasis/microcoverage

max_issues_repo_name: StephenVavasis/microcoverage

max_forks_repo_name: StephenVavasis/microcoverage


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

using Libdl

global snaplib

function loadsnap(mdppath::String)
    libexist = 0
    if Sys.isapple()
        libpath = mdppath * "src/Julia/Potential/cpuSnap.dylib";
    elseif Sys.isunix()
        libpath = mdppath * "src/Julia/Potential/cpuSnap.so";
    elseif Sys.windows()
        libpath = mdppath * "src/Julia/Potential/cpuSnap.so";
    end
    if isfile(libpath)
        libexist = 1;
    end
    if libexist==0
        cdir = pwd();
        cd(mdppath * "src/Julia/Potential")
        cstr = `g++ -std=c++11 -Wall -Wextra -pedantic -c -fPIC cpuSnap.cpp -o cpusnap.o`
        run(cstr)
        if Sys.isapple()
            cstr = `g++ -shared cpusnap.o -o cpuSnap.dylib`
        else
            cstr = `g++ -shared cpusnap.o -o cpuSnap.so`
        end
        run(cstr)
        cd(cdir)
    end
    global snaplib = Libdl.dlopen(libpath)
end

function closesnap()
    Libdl.dlclose(snaplib)
end

function bnum(twojmax::Int32)
    count = ccall(Libdl.dlsym(snaplib, :idxbcount), Cint, (Cint,), twojmax)
    return count
end

function unum(twojmax::Int32)
    count = ccall(Libdl.dlsym(snaplib, :idxucount), Cint, (Cint,), twojmax)
    return count
end

function znum(twojmax::Int32)
    count = ccall(Libdl.dlsym(snaplib, :idxzcount), Cint, (Cint,), twojmax)
    return count
end

function cgnum(twojmax::Int32)
    count = ccall(Libdl.dlsym(snaplib, :idxcgcount), Cint, (Cint,), twojmax)
    return count
end

function factable()

    facto = [
        1,
        1,
        2,
        6,
        24,
        120,
        720,
        5040,
        40320,
        362880,
        3628800,
        39916800,
        479001600,
        6227020800,
        87178291200,
        1307674368000,
        20922789888000,
        355687428096000,
        6.402373705728e+15,
        1.21645100408832e+17,
        2.43290200817664e+18,
        5.10909421717094e+19,
        1.12400072777761e+21,
        2.5852016738885e+22,
        6.20448401733239e+23,
        1.5511210043331e+25,
        4.03291461126606e+26,
        1.08888694504184e+28,
        3.04888344611714e+29,
        8.8417619937397e+30,
        2.65252859812191e+32,
        8.22283865417792e+33,
        2.63130836933694e+35,
        8.68331761881189e+36,
        2.95232799039604e+38,
        1.03331479663861e+40,
        3.71993326789901e+41,
        1.37637530912263e+43,
        5.23022617466601e+44,
        2.03978820811974e+46,
        8.15915283247898e+47,
        3.34525266131638e+49,
        1.40500611775288e+51,
        6.04152630633738e+52,
        2.65827157478845e+54,
        1.1962222086548e+56,
        5.50262215981209e+57,
        2.58623241511168e+59,
        1.24139155925361e+61,
        6.08281864034268e+62,
        3.04140932017134e+64,
        1.55111875328738e+66,
        8.06581751709439e+67,
        4.27488328406003e+69,
        2.30843697339241e+71,
        1.26964033536583e+73,
        7.10998587804863e+74,
        4.05269195048772e+76,
        2.35056133128288e+78,
        1.3868311854569e+80,
        8.32098711274139e+81,
        5.07580213877225e+83,
        3.14699732603879e+85,
        1.98260831540444e+87,
        1.26886932185884e+89,
        8.24765059208247e+90,
        5.44344939077443e+92,
        3.64711109181887e+94,
        2.48003554243683e+96,
        1.71122452428141e+98,
        1.19785716699699e+100,
        8.50478588567862e+101,
        6.12344583768861e+103,
        4.47011546151268e+105,
        3.30788544151939e+107,
        2.48091408113954e+109,
        1.88549470166605e+111,
        1.45183092028286e+113,
        1.13242811782063e+115,
        8.94618213078297e+116,
        7.15694570462638e+118,
        5.79712602074737e+120,
        4.75364333701284e+122,
        3.94552396972066e+124,
        3.31424013456535e+126,
        2.81710411438055e+128,
        2.42270953836727e+130,
        2.10775729837953e+132,
        1.85482642257398e+134,
        1.65079551609085e+136,
        1.48571596448176e+138,
        1.3520015276784e+140,
        1.24384140546413e+142,
        1.15677250708164e+144,
        1.08736615665674e+146,
        1.03299784882391e+148,
        9.91677934870949e+149,
        9.61927596824821e+151,
        9.42689044888324e+153,
        9.33262154439441e+155,
        9.33262154439441e+157,
        9.42594775983835e+159,
        9.61446671503512e+161,
        9.90290071648618e+163,
        1.02990167451456e+166,
        1.08139675824029e+168,
        1.14628056373471e+170,
        1.22652020319614e+172,
        1.32464181945183e+174,
        1.44385958320249e+176,
        1.58824554152274e+178,
        1.76295255109024e+180,
        1.97450685722107e+182,
        2.23119274865981e+184,
        2.54355973347219e+186,
        2.92509369349301e+188,
        3.3931086844519e+190,
        3.96993716080872e+192,
        4.68452584975429e+194,
        5.5745857612076e+196,
        6.68950291344912e+198,
        8.09429852527344e+200,
        9.8750442008336e+202,
        1.21463043670253e+205,
        1.50614174151114e+207,
        1.88267717688893e+209,
        2.37217324288005e+211,
        3.01266001845766e+213,
        3.8562048236258e+215,
        4.97450422247729e+217,
        6.46685548922047e+219,
        8.47158069087882e+221,
        1.118248651196e+224,
        1.48727070609069e+226,
        1.99294274616152e+228,
        2.69047270731805e+230,
        3.65904288195255e+232,
        5.01288874827499e+234,
        6.91778647261949e+236,
        9.61572319694109e+238,
        1.34620124757175e+241,
        1.89814375907617e+243,
        2.69536413788816e+245,
        3.85437071718007e+247,
        5.5502938327393e+249,
        8.04792605747199e+251,
        1.17499720439091e+254,
        1.72724589045464e+256,
        2.55632391787286e+258,
        3.80892263763057e+260,
        5.71338395644585e+262,
        8.62720977423323e+264,
        1.31133588568345e+267,
        2.00634390509568e+269,
        3.08976961384735e+271,
        4.78914290146339e+273,
        7.47106292628289e+275,
        1.17295687942641e+278,
        1.85327186949373e+280,
        2.94670227249504e+282,
        4.71472363599206e+284,
        7.59070505394721e+286,
        1.22969421873945e+289,
        2.0044015765453e+291,
        3.28721858553429e+293,
        5.42391066613159e+295,
        9.00369170577843e+297,
        1.503616514865e+300,
    ]

    facto = Float64.(Array(facto))

    return facto

end

function initsna(param, elemradius, elemweight, snapcoeff=nothing)

    ntypes = Int32(param[1]);
    twojmax = Int32(param[2]);
    rcutfac = param[3];
    rfac0 = param[4];
    rmin0 = param[5];
    bzeroflag = Int32(param[6]);
    switchflag = Int32(param[7]);
    quadraticflag = Int32(param[8]);
    chemflag = Int32(param[9]);
    bnormflag = Int32(param[10]);
    wselfallflag = Int32(param[11]);

    sna = initializesna(ntypes, twojmax, chemflag)

    facto = factable()
    ccall(Libdl.dlsym(snaplib, :cpuInitSna), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
            Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint),
            sna.rootpqarray, sna.cglist, facto, sna.idx_max, sna.idxz, sna.idxz_block,
            sna.idxb, sna.idxb_block, sna.idxu_block, sna.idxcg_block, sna.twojmax)

    sna.radelem[2:end] = elemradius
    sna.wjelem[2:end] = elemweight
    sna.map[2:end] = Int32.(Array(1:ntypes) .- 1)

    if snapcoeff === nothing
        sna.coeffelem = zeros(Float64, sna.ncoeffall)
    elseif length(snapcoeff) == sna.ntypes*sna.ncoeffall
        sna.coeffelem = snapcoeff
    else
        display([length(snapcoeff) sna.ntypes*sna.ncoeffall])
        error("number of snap coefficients is incompatible with the snap paramaters")
    end

    rcutmax = 0.0
    for ielem = 1:ntypes
        rcutmax = max(2.0*elemradius[ielem]*rcutfac,rcutmax)
    end

    for i = 1:ntypes
        for j = 1:ntypes
            cut = (elemradius[i] + elemradius[j])*rcutfac
            sna.rcutsq[j + (i-1)*(ntypes)] = cut*cut
        end
    end

    wself = 1.0;
    www = wself*wself*wself;
    for j = 0:twojmax
        if (bnormflag == 1)
            sna.bzero[j+1] = www;
        else
            sna.bzero[j+1] = www*(j+1);
        end
    end

    sna.nperdim = sna.ncoeffall-1;
    sna.bnormflag = bnormflag;
    sna.chemflag = chemflag;
    sna.quadraticflag = quadraticflag;
    sna.switchflag = switchflag;
    sna.bzeroflag = bzeroflag;
    sna.wselfallflag = wselfallflag;
    sna.wself = wself;
    sna.rmin0 = rmin0;
    sna.rfac0 = rfac0;
    sna.rcutfac = rcutfac;
    sna.rcutmax = rcutmax;

    return sna
end

function snappotential(x, t, a, b, c, pbc, param, elemradius, elemweight, snapcoeff)

    sna = initsna(param, elemradius, elemweight, snapcoeff)
    eatom, fatom, vatom = snappotential(x, t, a, b, c, pbc, sna)
    return eatom, fatom, vatom
end

function snappotential(x, t, a, b, c, pbc, sna)

    dim, N = size(x)

    #idxcg_max = sna.idxcg_max;
    idxu_max = sna.idxu_max;
    idxb_max = sna.idxb_max;
    idxz_max = sna.idxz_max;
    twojmax = sna.twojmax;
    #ncoeff = sna.ncoeff;
    ncoeffall = sna.ncoeffall;
    ntypes = sna.ntypes;
    nelem = sna.nelements;
    # ndoubles = sna.ndoubles;
    # ntriples = sna.ntriples;
    # nperdim = sna.nperdim;
    bnormflag = sna.bnormflag;
    chemflag = sna.chemflag;
    quadraticflag = sna.quadraticflag;
    switchflag = sna.switchflag;
    bzeroflag = sna.bzeroflag;
    wselfallflag = sna.wselfallflag;

    map = sna.map;
    idxz = sna.idxz;
    #idxz_block = sna.idxz_block;
    idxb = sna.idxb;
    idxb_block = sna.idxb_block;
    idxu_block = sna.idxu_block;
    idxcg_block = sna.idxcg_block;

    rcutmax = sna.rcutmax
    wself = sna.wself;
    rmin0 = sna.rmin0;
    rfac0 = sna.rfac0;
    rcutfac = sna.rcutfac;
    rcutmax = sna.rcutmax;
    bzero = sna.bzero;
    rootpqarray = sna.rootpqarray;
    cglist = sna.cglist;
    rcutsq = sna.rcutsq;
    radelem = sna.radelem;
    wjelem = sna.wjelem;
    coeffelem = sna.coeffelem;

    y, alist, neighlist, neighnum = fullneighborlist(x, a, b, c, pbc, rcutmax);
    ilist = Int32.(Array(1:N));
    atomtype = Int32.(t[:])
    alist = Int32.(alist)
    neighlist = Int32.(neighlist)
    neighnum = Int32.(neighnum)

    #pairlist, pairnum = neighpairlist(y, ilist, neighlist, neighnum, rcutmax*rcutmax);
    pairlist, pairnum = neighpairlist(y, ilist, alist, atomtype, neighlist, neighnum, rcutsq, ntypes)
    rij, ai, aj, ti, tj = neighpairs(y, pairlist, pairnum, atomtype, ilist, alist);
    ijnum = pairnum[end]

    # offset 1 for C++ code
    ilist = ilist .- Int32(1)
    alist = alist .- Int32(1)
    ai = ai .- Int32(1)
    aj = aj .- Int32(1)
    N = Int32(N)

    ulisttot_r = zeros(N*idxu_max*nelem)
    ulisttot_i = zeros(N*idxu_max*nelem)
    ccall(Libdl.dlsym(snaplib, :cpuSnapComputeUi), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
        Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, Cdouble, Ptr{Cint}, Ptr{Cint},
        Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint, Cint),
        ulisttot_r, ulisttot_i, rootpqarray, rij, wjelem, radelem, rmin0, rfac0, rcutfac,
        map, ai, ti, tj, twojmax, idxu_max, N, ijnum, switchflag, chemflag)

    ccall(Libdl.dlsym(snaplib, :cpuAddWself2Ui), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Cdouble,
        Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint, Cint),
        ulisttot_r, ulisttot_i, wself, idxu_block, atomtype, map, wselfallflag, chemflag,
        idxu_max, nelem, twojmax, N)

    eatom = zeros(N)
    ccall(Libdl.dlsym(snaplib, :cpuSnapComputeEi), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
        Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
        Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint),
        eatom, ulisttot_r, ulisttot_i, cglist, bzero, coeffelem, ilist, map, atomtype,
        idxb, idxcg_block, idxu_block, twojmax, idxb_max, idxu_max, nelem,
        ncoeffall, bnormflag, bzeroflag, wselfallflag, chemflag, N)

    ylist_r = zeros(N*idxu_max*nelem)
    ylist_i = zeros(N*idxu_max*nelem)
    ccall(Libdl.dlsym(snaplib, :cpuSnapComputeYi), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
        Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
        Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint),
        ylist_r, ylist_i, ulisttot_r, ulisttot_i, cglist, coeffelem, map, atomtype, idxz, idxb_block,
        idxu_block, idxcg_block, twojmax, idxb_max, idxu_max, idxz_max, nelem, ncoeffall, bnormflag, chemflag, N)

    fatom = zeros(3,N)
    vatom = zeros(6,N)
    ccall(Libdl.dlsym(snaplib, :cpuSnapComputeFi), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
        Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, Cdouble,
        Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint,
        Cint, Cint, Cint, Cint),
        fatom, vatom, ylist_r, ylist_i, rootpqarray, rij, wjelem, radelem, rmin0, rfac0, rcutfac,
        map, ai, ai, aj, ti, tj, twojmax, idxu_max, N, N, ijnum, switchflag, chemflag)

    return eatom, fatom, vatom
end

function snappotential(x, t, alist, neighlist, neighnum, sna)

    dim = size(x,1);
    N = length(t)

    #idxcg_max = sna.idxcg_max;
    idxu_max = sna.idxu_max;
    idxb_max = sna.idxb_max;
    idxz_max = sna.idxz_max;
    twojmax = sna.twojmax;
    #ncoeff = sna.ncoeff;
    ncoeffall = sna.ncoeffall;
    ntypes = sna.ntypes;
    nelem = sna.nelements;
    # ndoubles = sna.ndoubles;
    # ntriples = sna.ntriples;
    # nperdim = sna.nperdim;
    bnormflag = sna.bnormflag;
    chemflag = sna.chemflag;
    quadraticflag = sna.quadraticflag;
    switchflag = sna.switchflag;
    bzeroflag = sna.bzeroflag;
    wselfallflag = sna.wselfallflag;

    map = sna.map;
    idxz = sna.idxz;
    #idxz_block = sna.idxz_block;
    idxb = sna.idxb;
    idxb_block = sna.idxb_block;
    idxu_block = sna.idxu_block;
    idxcg_block = sna.idxcg_block;

    rcutmax = sna.rcutmax
    wself = sna.wself;
    rmin0 = sna.rmin0;
    rfac0 = sna.rfac0;
    rcutfac = sna.rcutfac;
    rcutmax = sna.rcutmax;
    bzero = sna.bzero;
    rootpqarray = sna.rootpqarray;
    cglist = sna.cglist;
    rcutsq = sna.rcutsq;
    radelem = sna.radelem;
    wjelem = sna.wjelem;
    coeffelem = sna.coeffelem;

    #y, alist, neighlist, neighnum = fullneighborlist(x, a, b, c, pbc, rcutmax);
    ilist = Int32.(Array(1:N));
    atomtype = Int32.(t[:])
    alist = Int32.(alist)
    neighlist = Int32.(neighlist)
    neighnum = Int32.(neighnum)

    pairlist, pairnum = neighpairlist(x, ilist, alist, atomtype, neighlist, neighnum, rcutsq, ntypes)
    rij, ai, aj, ti, tj = neighpairs(x, pairlist, pairnum, atomtype, ilist, alist);
    ijnum = pairnum[end]

    # offset 1 for C++ code
    ilist = ilist .- Int32(1)
    alist = alist .- Int32(1)
    ai = ai .- Int32(1)
    aj = aj .- Int32(1)

    ulisttot_r = zeros(N*idxu_max*nelem)
    ulisttot_i = zeros(N*idxu_max*nelem)
    ccall(Libdl.dlsym(snaplib, :cpuSnapComputeUi), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
        Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, Cdouble, Ptr{Cint}, Ptr{Cint},
        Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint, Cint),
        ulisttot_r, ulisttot_i, rootpqarray, rij, wjelem, radelem, rmin0, rfac0, rcutfac,
        map, ai, ti, tj, twojmax, idxu_max, N, ijnum, switchflag, chemflag)

    ccall(Libdl.dlsym(snaplib, :cpuAddWself2Ui), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Cdouble,
        Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint, Cint),
        ulisttot_r, ulisttot_i, wself, idxu_block, atomtype, map, wselfallflag, chemflag,
        idxu_max, nelem, twojmax, N)

    eatom = zeros(N)
    ccall(Libdl.dlsym(snaplib, :cpuSnapComputeEi), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
        Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
        Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint),
        eatom, ulisttot_r, ulisttot_i, cglist, bzero, coeffelem, ilist, map, atomtype,
        idxb, idxcg_block, idxu_block, twojmax, idxb_max, idxu_max, nelem,
        ncoeffall, bnormflag, bzeroflag, wselfallflag, chemflag, N)

    ylist_r = zeros(N*idxu_max*nelem)
    ylist_i = zeros(N*idxu_max*nelem)
    ccall(Libdl.dlsym(snaplib, :cpuSnapComputeYi), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
        Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
        Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint),
        ylist_r, ylist_i, ulisttot_r, ulisttot_i, cglist, coeffelem, map, atomtype, idxz, idxb_block,
        idxu_block, idxcg_block, twojmax, idxb_max, idxu_max, idxz_max, nelem, ncoeffall, bnormflag, chemflag, N)

    fatom = zeros(3,N)
    vatom = zeros(6,N)
    ccall(Libdl.dlsym(snaplib, :cpuSnapComputeFi), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
        Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, Cdouble,
        Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint,
        Cint, Cint, Cint, Cint),
        fatom, vatom, ylist_r, ylist_i, rootpqarray, rij, wjelem, radelem, rmin0, rfac0, rcutfac,
        map, ai, ai, aj, ti, tj, twojmax, idxu_max, N, N, ijnum, switchflag, chemflag)

    return eatom, fatom, vatom
end

function snapdescriptors(x, t, a, b, c, pbc, param, elemradius, elemweight)

    sna = initsna(param, elemradius, elemweight)
    bi, bd, bv = snapdescriptors(x, t, a, b, c, pbc, sna)
    return bi, bd, bv
end

function snapdescriptors(x, t, a, b, c, pbc, sna)

    dim, N = size(x)

    #idxcg_max = sna.idxcg_max;
    idxu_max = sna.idxu_max;
    idxb_max = sna.idxb_max;
    idxz_max = sna.idxz_max;
    twojmax = sna.twojmax;
    ncoeff = sna.ncoeff;
    #ncoeffall = sna.ncoeffall;
    ntypes = sna.ntypes;
    nelem = sna.nelements;
    ndoubles = sna.ndoubles;
    ntriples = sna.ntriples;
    #nperdim = sna.nperdim;
    bnormflag = sna.bnormflag;
    chemflag = sna.chemflag;
    quadraticflag = sna.quadraticflag;
    switchflag = sna.switchflag;
    bzeroflag = sna.bzeroflag;
    wselfallflag = sna.wselfallflag;

    map = sna.map;
    idxz = sna.idxz;
    idxz_block = sna.idxz_block;
    idxb = sna.idxb;
    #idxb_block = sna.idxb_block;
    idxu_block = sna.idxu_block;
    idxcg_block = sna.idxcg_block;

    rcutmax = sna.rcutmax
    wself = sna.wself;
    rmin0 = sna.rmin0;
    rfac0 = sna.rfac0;
    rcutfac = sna.rcutfac;
    rcutmax = sna.rcutmax;
    bzero = sna.bzero;
    rootpqarray = sna.rootpqarray;
    cglist = sna.cglist;
    rcutsq = sna.rcutsq;
    radelem = sna.radelem;
    wjelem = sna.wjelem;
    #coeffelem = sna.coeffelem;

    y, alist, neighlist, neighnum = fullneighborlist(x, a, b, c, pbc, rcutmax);
    ilist = Int32.(Array(1:N));
    atomtype = Int32.(t[:])
    alist = Int32.(alist)
    neighlist = Int32.(neighlist)
    neighnum = Int32.(neighnum)

    pairlist, pairnum = neighpairlist(y, ilist, alist, atomtype, neighlist, neighnum, rcutsq, ntypes)
    rij, ai, aj, ti, tj = neighpairs(y, pairlist, pairnum, atomtype, ilist, alist);
    ijnum = pairnum[end]

    # offset 1 for C++ code
    ilist = ilist .- Int32(1)
    alist = alist .- Int32(1)
    ai = ai .- Int32(1)
    aj = aj .- Int32(1)

    n = idxu_max*ijnum
    ulist_r = zeros(n)
    ulist_i = zeros(n)
    dulist_r = zeros(3*n)
    dulist_i = zeros(3*n)
    ccall(Libdl.dlsym(snaplib, :cpuComputeUij), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
        Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble,
        Cdouble, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint, Cint),
        ulist_r, ulist_i, dulist_r, dulist_i, rootpqarray, rij, wjelem, radelem, rmin0,
        rfac0, rcutfac, idxu_block, ti, tj, twojmax, idxu_max, ijnum, switchflag)

    ulisttot_r = zeros(N*idxu_max*nelem)
    ulisttot_i = zeros(N*idxu_max*nelem)
    ccall(Libdl.dlsym(snaplib, :cpuZeroUarraytot2), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Cdouble,
        Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},  Cint, Cint, Cint, Cint, Cint, Cint),
        ulisttot_r, ulisttot_i, wself, idxu_block, atomtype, map, ai, wselfallflag, chemflag,
        idxu_max, nelem, twojmax, N)

    ccall(Libdl.dlsym(snaplib, :cpuAddUarraytot), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
        Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint, Cint),
        ulisttot_r, ulisttot_i, ulist_r, ulist_i, map, ai, tj, idxu_max, N, ijnum, chemflag)

    zlist_r = zeros(idxz_max*ndoubles*N)
    zlist_i = zeros(idxz_max*ndoubles*N)
    ccall(Libdl.dlsym(snaplib, :cpuComputeZi2), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
        Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint, Cint),
        zlist_r, zlist_i, ulisttot_r, ulisttot_i, cglist, idxz, idxu_block, idxcg_block, twojmax,
        idxu_max, idxz_max, nelem, bnormflag, N)

    blist = zeros(idxb_max*ntriples*N)
    ccall(Libdl.dlsym(snaplib, :cpuComputeBi2), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
        Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
        Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint),
        blist, zlist_r, zlist_i, ulisttot_r, ulisttot_i, bzero, ilist, atomtype,
        map, idxb, idxu_block, idxz_block, twojmax, idxb_max, idxu_max, idxz_max,
        nelem, bzeroflag, wselfallflag, chemflag, N)

    dblist = zeros(idxb_max*ntriples*dim*ijnum)
    ccall(Libdl.dlsym(snaplib, :cpuComputeDbidrj), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
        Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
        Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint),
        dblist, zlist_r, zlist_i, dulist_r, dulist_i, idxb, idxu_block, idxz_block, map, ai, tj,
        twojmax, idxb_max, idxu_max, idxz_max, nelem, bnormflag, chemflag, N, ijnum)

    bi = zeros(ntypes*ncoeff)
    ccall(Libdl.dlsym(snaplib, :cpuSnapTallyBispectrum), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble},
        Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint),
        bi, blist, ilist, atomtype, N, ncoeff, ncoeff, ntypes, quadraticflag)

    bd = zeros(ntypes*ncoeff*3*N)
    ccall(Libdl.dlsym(snaplib, :cpuSnapTallyBispectrumDeriv), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble},
        Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint, Cint),
        bd, blist, dblist, ai, ai, aj, ti, N, ijnum, ncoeff, ncoeff, ntypes, quadraticflag)

    bv = zeros(6*ntypes*ncoeff)
    ccall(Libdl.dlsym(snaplib, :cpuSnapTallyBispectrumVirial), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
        Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint, Cint),
        bv, blist, dblist, rij, ai, ti, N, ijnum, ncoeff, ncoeff, ntypes, quadraticflag)

    bd = reshape(bd, (3,N,ntypes*ncoeff))
    bv = reshape(bv, (6,ntypes*ncoeff))
    bv = bv[[1; 2; 3; 6; 5; 4],:]

    return bi, bd, bv
end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using Test
using EzXML
using BacenDRM

function document_example()
    # create document with basic information
    doc = BacenDRM.Documento(
        "2060",              # id_docto::String,
        "v1",                # id_docto_versao::String,
        "2020-06",           # data_base::String,
        123456,              # id_inst_financ::Int64,
        :I,                  # tipo_arq::Symbol,
        "Fulano",            # nome_contato::String,
        "555-1234",          # fone_contato::String
    )

    # create ItemCarteira, which is an index to document data
    item_carteira = BacenDRM.ItemCarteira(
        :A20,
        nothing,
        :JM1,
        :offshore,
        :banking
    )

    # repeat for each vertice
    codigo_vertice = 1
    vertice_1 = BacenDRM.Vertice(
        100_000.0, # valor_alocado::Float64
          0.0      # valor_mam::Float64
    )
    push!(doc, item_carteira, codigo_vertice, vertice_1) # add Vertice
    # end repeat

    # two vertices
    ic = BacenDRM.ItemCarteira(:A30, nothing, :ME1, :offshore, :banking)
    push!(doc, ic,  3, BacenDRM.Vertice(100_000.0,  0.0))
    push!(doc, ic, 12, BacenDRM.Vertice(200_000.0, 10_000.0))

    # you can also build a vertice adding incrementally its values
    ic = BacenDRM.ItemCarteira(:P30, nothing, :JM1, :onshore_sem_clearing, :trading)
    push!(doc, ic, 1, BacenDRM.Vertice(70_000.0, 0.0))
    push!(doc, ic, 1, BacenDRM.Vertice(230_000.0, 0.0))
    # will result in: Vertice(300_000.0, 0.0)

    # one-liner
    push!(doc, BacenDRM.ItemCarteira(:D41, :C, :JM1, :onshore_clearing, :banking), 1, BacenDRM.Vertice(400_000.0, 0.0))

    # hedge fund assets: data goes into ativo_fundo section
    ic = BacenDRM.ItemCarteira(
        :A90,
        nothing,
        :JM1,
        :offshore,
        :banking,
        true      # is_ativo_fundo::Bool=false
    )
    push!(doc, ic, 1, BacenDRM.Vertice(500_000.0, 0.0))
    push!(doc, BacenDRM.ItemCarteira(:P10, nothing, :JI1, :offshore, :trading, true), 1, BacenDRM.Vertice(700_000.0, 0.0))
    push!(doc, BacenDRM.ItemCarteira(:D20, nothing, :ME3, :offshore, :trading, true), 1, BacenDRM.Vertice(800_000.0, 0.0))

    # atividade_financeira
    push!(doc, BacenDRM.ItemCarteira(:AFC, :V, :JM1, nothing, :banking), 1, BacenDRM.Vertice(600_000.0, 0.0))

    return doc
end

@testset "Vertice" begin

    v1 = BacenDRM.Vertice(100_000.0, 0.0)
    @test v1.valor_alocado == 100_000.0
    @test v1.valor_mam == 0.0

    v12 = BacenDRM.Vertice(200_000.0, 20_000.0)
    @test v12.valor_alocado == 200_000.0
    @test v12.valor_mam == 20_000.0

end

@testset "Fluxo" begin

    fluxo = BacenDRM.Fluxos()
    fluxo[1] = BacenDRM.Vertice(10_000.0, 0.0)

    @test fluxo[1].valor_alocado == 10_000.0
    @test fluxo[Symbol("01")].valor_mam == 0.0

    fluxo[Symbol("12")] = BacenDRM.Vertice(10.0, 20.0)

    @test_throws AssertionError fluxo[2]  = BacenDRM.Vertice(0.0, 1_000.0) # invalid valor_mam
    @test_throws AssertionError fluxo[Symbol("3")]  = BacenDRM.Vertice(0.0, 1_000.0) # invalid valor_mam
    @test_throws AssertionError fluxo[13] = BacenDRM.Vertice(0.0, 1_000.0) # invalid codigo_vertice
    @test_throws AssertionError fluxo[Symbol("14")] = BacenDRM.Vertice(0.0, 1_000.0) # invalid codigo_vertice

end

@testset "Item Carteira" begin

    ic = BacenDRM.ItemCarteira(
        :A20,              # item::Symbol
        nothing,           # id_posicao::SymbolOrNothing
        :JM1,              # fator_risco::Symbol
        :onshore_clearing, # local_registro::SymbolOrNothing
        :trading          # carteira_negoc::Symbol
    )

    @test ic.item == :A20
    @test BacenDRM.is_ativo(ic)
    @test !BacenDRM.is_passivo(ic)
    @test !BacenDRM.is_derivativo(ic)
    @test !BacenDRM.is_ativo_fundo(ic)
    @test !BacenDRM.is_passivo_fundo(ic)
    @test !BacenDRM.is_derivativo_fundo(ic)
    @test !BacenDRM.is_atividade_financeira(ic)

    ic2 = BacenDRM.ItemCarteira(
        :D41,              # item::Symbol
        :C,                # id_posicao::SymbolOrNothing
        :JM1,              # fator_risco::Symbol
        :onshore_clearing, # local_registro::SymbolOrNothing
        :trading,          # carteira_negoc::Symbol
        true               # is_abertura_fundo::Bool=false
    )

    @test ic2.fator_risco == :JM1
    @test !BacenDRM.is_ativo_fundo(ic2)
    @test !BacenDRM.is_passivo_fundo(ic2)
    @test BacenDRM.is_derivativo_fundo(ic2)
    @test !BacenDRM.is_ativo(ic2)
    @test !BacenDRM.is_passivo(ic2)
    @test !BacenDRM.is_derivativo(ic2)
    @test !BacenDRM.is_atividade_financeira(ic2)

end

@testset "Documento" begin

    doc = BacenDRM.Documento(
        "2060",              # id_docto::String,
        "v1",                # id_docto_versao::String,
        "2020-06",           # data_base::String,
        123456,              # id_inst_financ::Int64,
        :I,                  # tipo_arq::Symbol,
        "Fulano",            # nome_contato::String,
        "555-1234"           # fone_contato::String,
    )

    # ativo :A20
    ic = BacenDRM.ItemCarteira(:A20, nothing, :JM1, :offshore, :banking)
    push!(doc, ic, 1, BacenDRM.Vertice(100_000.0, 0.0))
    # ativo :A30
    ic = BacenDRM.ItemCarteira(:A30, nothing, :ME1, :offshore, :banking)
    push!(doc, ic,            3, BacenDRM.Vertice(100_000.0, 0.0))  # integer
    push!(doc, ic, Symbol("12"), BacenDRM.Vertice(100_000.0, 10_000.0)) # or symbol

    # passivo
    push!(doc, BacenDRM.ItemCarteira(:P30, nothing, :JM1, :onshore_sem_clearing, :trading), 1, BacenDRM.Vertice(100_000.0, 0.0))
    # derivativo
    push!(doc, BacenDRM.ItemCarteira(:D41, :C, :JM1, :onshore_clearing, :banking), 1, BacenDRM.Vertice(100_000.0, 0.0))
    # ativo_fundo
    ic = BacenDRM.ItemCarteira(:A90, nothing, :JM1, :offshore, :banking, true) # is_abertura_fundo=true
    push!(doc, ic, 1, BacenDRM.Vertice(100_000.0, 0.0))
    # atividade_financeira
    push!(doc, BacenDRM.ItemCarteira(:AFC, :V, :JM1, nothing, :banking), 1, BacenDRM.Vertice(100_000.0, 0.0))

    @test doc.id_docto == "2060"
    @test length(doc.ativo) == 2

end

@testset "Write XML" begin

    doc = document_example()

    function _is_equal(doc::BacenDRM.Documento, str_xml::String)
        b = IOBuffer()
        BacenDRM.write_xml(b, doc)
        xml1 = String(take!(b))

        b2 = IOBuffer()
        parsed = parsexml(str_xml)
        prettyprint(b2, parsed)
        xml2 = String(take!(b2))

        #println(xml1) # debug
        #println(xml2) # debug

        return xml1 == xml2

    end

    str_xml = """
<?xml version="1.0" encoding="UTF-8"?>
<DocDRM>
  <IdDocto>2060</IdDocto>
  <IdDoctoVersao>v1</IdDoctoVersao>
  <DataBase>2020-06</DataBase>
  <IdInstFinanc>123456</IdInstFinanc>
  <TipoArq>I</TipoArq>
  <NomeContato>Fulano</NomeContato>
  <FoneContato>555-1234</FoneContato>
  <Ativo>
    <ItemCarteira Item="A20" FatorRisco="JM1" LocalRegistro="03" CarteiraNegoc="02">
      <FluxoVertice CodVertice="01" ValorAlocado="100"/>
    </ItemCarteira>
    <ItemCarteira Item="A30" FatorRisco="ME1" LocalRegistro="03" CarteiraNegoc="02">
      <FluxoVertice CodVertice="03" ValorAlocado="100"/>
      <FluxoVertice CodVertice="12" ValorAlocado="200" ValorMaM="10"/>
    </ItemCarteira>
  </Ativo>
  <Passivo>
    <ItemCarteira Item="P30" FatorRisco="JM1" LocalRegistro="02" CarteiraNegoc="01">
      <FluxoVertice CodVertice="01" ValorAlocado="300"/>
    </ItemCarteira>
  </Passivo>
  <Derivativo>
    <ItemCarteira Item="D41" IdPosicao="C" FatorRisco="JM1" LocalRegistro="01" CarteiraNegoc="02">
      <FluxoVertice CodVertice="01" ValorAlocado="400"/>
    </ItemCarteira>
  </Derivativo>
  <AtivoFundo>
    <ItemCarteira Item="A90" FatorRisco="JM1" LocalRegistro="03" CarteiraNegoc="02">
      <FluxoVertice CodVertice="01" ValorAlocado="500"/>
    </ItemCarteira>
  </AtivoFundo>
  <PassivoFundo>
    <ItemCarteira Item="P10" FatorRisco="JI1" LocalRegistro="03" CarteiraNegoc="01">
      <FluxoVertice CodVertice="01" ValorAlocado="700"/>
    </ItemCarteira>
  </PassivoFundo>
  <DerivativoFundo>
    <ItemCarteira Item="D20" FatorRisco="ME3" LocalRegistro="03" CarteiraNegoc="01">
      <FluxoVertice CodVertice="01" ValorAlocado="800"/>
    </ItemCarteira>
  </DerivativoFundo>
  <AtividadeFinanceira>
    <ItemCarteira Item="AFC" IdPosicao="V" FatorRisco="JM1" CarteiraNegoc="02">
      <FluxoVertice CodVertice="01" ValorAlocado="600"/>
    </ItemCarteira>
  </AtividadeFinanceira>
</DocDRM>
    """

    @test _is_equal(doc, str_xml)

    ic = BacenDRM.ItemCarteira(:P30, nothing, :JM1, :onshore_sem_clearing, :trading)
    @test doc.passivo[ic][Symbol("01")].valor_alocado == 300_000.0

    #BacenDRM.write_xml(Base.stdout, doc) # debug

    file = tempname()
    BacenDRM.write_xml(file, doc)
    rm(file)

end

@testset "Write CSV" begin

    doc = document_example()

    #BacenDRM.write_csv(Base.stdout, doc) # debug
    #BacenDRM.write_csv("drm.csv", doc)   # debug

    file = tempname()
    BacenDRM.write_csv(file, doc)
    @test isfile(file)
    rm(file)

end

@testset "ItemCarteira Helpers" begin

    # item::Symbol
    # id_posicao::SymbolOrNothing
    # fator_risco::Symbol
    # local_registro::SymbolOrNothing
    # carteira_negoc::Symbol
    # BacenDRM.Fluxos::Vector{FluxoVertice}

    # equal
    @test BacenDRM.ItemCarteira(:A10, nothing, :JJ1, :offshore, :trading) ==
          BacenDRM.ItemCarteira(:A10, nothing, :JJ1, :offshore, :trading)

    @test BacenDRM.ItemCarteira(:A10, nothing, :JJ1, nothing, :trading) ==
          BacenDRM.ItemCarteira(:A10, nothing, :JJ1, nothing, :trading)

    # different
    @test BacenDRM.ItemCarteira(:A10, nothing, :JJ1, :offshore, :trading) !=
          BacenDRM.ItemCarteira(:P10, nothing, :JJ1, :offshore, :trading)

    @test BacenDRM.ItemCarteira(:A10, nothing, :JJ1, :offshore, :trading) !=
          BacenDRM.ItemCarteira(:A10, :C, :JJ1, :offshore, :trading)

    @test BacenDRM.ItemCarteira(:A10, nothing, :JJ1, :offshore, :trading) !=
          BacenDRM.ItemCarteira(:A10, nothing, :JP1, :offshore, :trading)

    @test BacenDRM.ItemCarteira(:A10, nothing, :JJ1, :offshore, :trading) !=
          BacenDRM.ItemCarteira(:A10, nothing, :JJ1, :onshore_clearing, :trading)

    @test BacenDRM.ItemCarteira(:A10, nothing, :JJ1, :offshore, :trading) !=
          BacenDRM.ItemCarteira(:A10, nothing, :JJ1, :offshore, :banking)

    @test BacenDRM.ItemCarteira(:A10, nothing, :JJ1, :offshore, :trading) !=
          BacenDRM.ItemCarteira(:A10, nothing, :JJ1, nothing, :trading)

    @test BacenDRM.ItemCarteira(:A10, :C, :JJ1, :offshore, :trading) !=
          BacenDRM.ItemCarteira(:A10, :V, :JJ1, :offshore, :trading)

end


text:
function get_model(re::Chain, p)
  m = re(p)
  reset_state!(m, p)
  (x, p) -> m(x)
end
function get_model(m::FastChain, p)
  reset_state!(m, p)
  m
end

function get_model(re, p)
  m = re(p)
  reset_state!(m, p)
  (x, p) -> m(x)
end

function loss_seq(p, m, x, y)
  m = get_model(m,p)

  ŷb = Flux.Zygote.Buffer([y[1]], size(y,1))
  @inbounds for i in 1:size(x,1)
    xi = x[i]
    ŷi = m(xi, p)
    (Inf32 ∈ ŷi || NaN ∈ ŷi) && return Inf32, y, y # TODO: what if a layer after MTKRecur can handle Infs?
    ŷb[i] = ŷi
  end
  ŷ = copy(ŷb)

  # ŷ = m.(x,[p])
  # Inf32 ∈ ŷ && return Inf32, ŷ, y

  return mean(Flux.Losses.mse.(ŷ,y, agg=mean)), ŷ, y
end

# function loss_3d_seq(p, m, x, y)
#   m = get_model(m,p)
#   ŷ = m(x, p)
#   Inf ∈ ŷ && return Inf32, ŷ, y
#   l = mean(Flux.Losses.mse.(ŷ,y, agg=mean))
#   return l, ŷ, y
# end


text:
"""
    MLEstimator

Simple maximum likelihood estimator for response rate `p`.
"""
struct MLEstimator{TD<:Design} <: Estimator

  design::TD

  MLEstimator{TD}(design::TD) where {TD<:Design} = new(design)

end # MLEStimator

"""
    MLEstimator(design::TD) where {TD<:Design}

Create maximum likelihood estimator for response rate under given design.
"""
MLEstimator(design::TD) where {TD<:Design} = MLEstimator{TD}(design)

Base.show(io::IO, estimator::MLEstimator) = print("MLEstimator")

function estimate(estimator::MLEstimator, x1::T, x2::T) where {T<:Integer}

    checkx1x2(x1, x2, design(estimator))
    return (x1 + x2)/samplesize(design(estimator), x1)

end

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
# Autogenerated wrapper script for libspecfemUtils_jll for x86_64-apple-darwin-libgfortran4
export libspecfemUtils

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("libspecfemUtils")
JLLWrappers.@declare_library_product(libspecfemUtils, "@rpath/libspecfemUtils.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libspecfemUtils,
        "lib/libspecfemUtils.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
abstract type AbstractTerminator end

struct Terminator <: AbstractTerminator
    proc::Base.Process
end

struct NullTerminator <: AbstractTerminator end

function check_supported()
    # "OS" APIs used in ./server.jl
    isdir("/proc/$(getpid())") || return false
    try
        read(`kill -l`)
    catch
        return false
    end
    return true
end

const IS_SUPPORTED = check_supported()

function Terminators.open(f)
    p = Terminators.open()
    try
        f()
    finally
        close(p)
    end
end

function Terminators.open()
    IS_SUPPORTED || return NullTerminator()
    julia = Base.julia_cmd()
    script = joinpath(@__DIR__, "server.jl")
    cmd = `$julia --startup-file=no $script -- $(getpid())`
    proc = open(cmd; read = true, write = true)
    return Terminator(proc)
end

function Base.close(p::Terminator)
    close(p.proc)
    wait(p.proc)
end

Base.close(::NullTerminator) = nothing

const LOCK = ReentrantLock()
const REF = Base.Ref{Union{Nothing,Terminator}}(nothing)

function global_terminator()
    IS_SUPPORTED || return NullTerminator()
    lock(LOCK) do
        p = REF[]
        p isa Terminator && return p
        p = REF[] = Terminators.open()
        atexit() do
            close(p)
        end
        return p
    end
end

const TIMER_ID = Ref{UInt}(0)

function wait_ack(p::Terminator, id::UInt)
    ln = readline(p.proc)
    @assert ln == "ack $id"
end

function Terminators.start(
    timeout::Number,
    p::AbstractTerminator = global_terminator();
    label = "",
)
    p isa NullTerminator && return UInt(0)
    timeout = Float64(timeout)
    id = TIMER_ID[] += 1
    write(p.proc, "start $id $timeout $label\n")
    flush(p.proc)
    wait_ack(p, id)
    return id
end

function Terminators.stop(id::UInt, p::AbstractTerminator = global_terminator())
    p isa NullTerminator && return UInt(0)
    write(p.proc, "stop $id 0 \n")
    flush(p.proc)
    wait_ack(p, id)
    return id
end

function Terminators.withtimeout(
    f,
    timeout::Number,
    p::AbstractTerminator = global_terminator();
    kwargs...,
)
    id = Terminators.start(timeout, p; kwargs...)
    try
        return f()
    finally
        Terminators.stop(id, p)
    end
end


text:
# OffsetArrays uses the Delegation pattern to reuse existing
# functionalities from AbstractArray.

# A regular array is indexed from 1 to length of the array

x = rand(3)
x[1:3]
#=
julia> x = rand(3)
3-element Array{Float64,1}:
 0.8261642367634516
 0.8348724419409639
 0.9849442213235109

julia> x[1:3]
3-element Array{Float64,1}:
 0.8261642367634516
 0.8348724419409639
 0.9849442213235109
=#

# An offset array can be indexed with different values
# Here's an example of a 0-based index

using OffsetArrays
y = OffsetArray(rand(3), 0:2)
y[0:2]
#=
julia> using OffsetArrays

julia> y = OffsetArray(rand(3), 0:2)
OffsetArray(::Array{Float64,1}, 0:2) with eltype Float64 with indices 0:2:
 0.5101137754975911
 0.18269142048901643
 0.9359844942154083

julia> y[0:2]
3-element Array{Float64,1}:
 0.5101137754975911
 0.18269142048901643
 0.9359844942154083
=#

# code sniplet from https://github.com/JuliaArrays/OffsetArrays.jl
#=

struct OffsetArray{T,N,AA<:AbstractArray} <: AbstractArray{T,N}
    parent::AA
    offsets::NTuple{N,Int}
end

Base.parent(A::OffsetArray) = A.parent

Base.size(A::OffsetArray) = size(parent(A))
Base.size(A::OffsetArray, d) = size(parent(A), d)

Base.eachindex(::IndexCartesian, A::OffsetArray) = CartesianIndices(axes(A))
Base.eachindex(::IndexLinear, A::OffsetVector)   = axes(A, 1)

=#


text:
using OpenCL, CLBLAS, CLBlast, BenchmarkTools
@static if VERSION < v"0.7-"
    LA = Base.LinAlg
else
    using LinearAlgebra, Random, Printf
    LA = LinearAlgebra
end

CLBLAS.setup()

function run_gemv(m::Integer, n::Integer, T=Float32)
    srand(12345)
    A = rand(T, m, n)
    y = rand(T, m)
    x = rand(T, n)
    α = rand(T)
    β = rand(T)

    @printf("\nm = %d, n = %d, eltype = %s\n", m, n, T)

    println("BLAS:")
    y_true = copy(y)
    LinAlg.BLAS.gemv!('N', α, A, x, β, y_true)
    yy = copy(y)
    display(@benchmark LinAlg.BLAS.gemv!($('N'), $α, $A, $x, $β, $yy))
    println()

    for device in cl.devices()
        println("-"^70)
        @printf("Platform name   : %s\n", device[:platform][:name])
        @printf("Platform version: %s\n", device[:platform][:version])
        @printf("Device name     : %s\n", device[:name])
        @printf("Device type     : %s\n", device[:device_type])
        println()

        ctx = cl.Context(device)
        queue = cl.CmdQueue(ctx)

        println("CLBLAS:")
        A_cl = cl.CLArray(queue, A)
        y_cl = cl.CLArray(queue, y)
        x_cl = cl.CLArray(queue, x)
        CLBLAS.gemv!('N', α, A_cl, x_cl, β, y_cl)
        @assert cl.to_host(y_cl) ≈ y_true
        display(@benchmark CLBLAS.gemv!($('N'), $α, $A_cl, $x_cl, $β, $y_cl))
        println()

        println("CLBlast:")
        A_cl = cl.CLArray(queue, A)
        y_cl = cl.CLArray(queue, y)
        x_cl = cl.CLArray(queue, x)
        CLBlast.gemv!('N', α, A_cl, x_cl, β, y_cl)
        @assert cl.to_host(y_cl) ≈ y_true
        display(@benchmark CLBlast.gemv!($('N'), $α, $A_cl, $x_cl, $β, $y_cl))
        println()
    end
end

run_gemv(2^10, 2^10, Float32)


text:
@block AnthonyWang [documentation] begin
    @cell "pong" [animated, scatter, updating] begin
        # init speed and velocity vector
        xyvec = rand(Point2f0, (2)) .* 5 .+ 1
        velvec = rand(Point2f0, (2)) .* 10
        # define some other parameters
        t = 0
        ts = 0.03
        balldiameter = 1
        origin = Point2f0(0, 0)
        xybounds = Point2f0(10, 10)
        N = 200
        scene = scatter(
            xyvec,
            markersize = balldiameter,
            color = rand(RGBf0, 2),
            limits = FRect(0, 0, xybounds)
        )
        s = scene[end] # last plot in scene

        record(scene, @outputfile(mp4), 1:N) do i
            # calculate new ball position
            global t = t + ts
            global xyvec = xyvec .+ velvec .* ts
            global velvec = map(xyvec, xybounds, origin, velvec) do p, b, o, vel
                boolvec = ((p .+ balldiameter/2) .> b) .| ((p .- balldiameter/2) .< o)
                velvec = map(boolvec, vel) do b, v
                    b ? -v : v
                end
            end
            # plot
            s[1] = xyvec
        end
    end

    @cell "pulsing marker" [animated, scatter, markersize, updating] begin
        N = 100
        scene = scatter([0], [0], marker = '❤', markersize = 0.5, color = :red, raw = true)
        s = scene[end] # last plot in scene
        record(scene, @outputfile(mp4), range(0, stop = 10pi, length = N)) do i
            s[:markersize] = (cos(i) + 1) / 4 + 0.2
        end
    end

    @cell "Travelling wave" [animated, lines, updating, interaction] begin
        scene = Scene()
        time = Node(0.0)
        f(v, t) = sin(v + t)
        scene = lines!(
            scene,
            lift(t -> f.(range(0, stop = 2pi, length = 50), t), time),
            color = :blue
        )
        p1 = scene[end];
        N = 100
        record(scene, @outputfile(mp4), range(0, stop = 4pi, length = N)) do i
            time[] = i
        end
    end

    @cell "Viridis scatter" ["2d", scatter, color, viridis, colormap] begin
        N = 30
        scatter(1:N, 1:N, markersize = 2, color = to_colormap(:viridis, N))
    end

    @cell "Viridis meshscatter" ["3d", scatter, color, viridis, colormap] begin
        N = 30
        R = 2
        theta = 4pi
        h = 5
        x = [R .* (t/3) .* cos(t) for t = range(0, stop = theta, length = N)]
        y = [R .* (t/3) .* sin(t) for t = range(0, stop = theta, length = N)]
        z = range(0, stop = h, length = N)
        meshscatter(x, y, z, markersize = 0.5, color = to_colormap(:viridis, N))
    end

    @cell "Marker sizes + Marker colors" ["2d", scatter, markersize, color] begin
        scatter(
            rand(20), rand(20),
            markersize = rand(20) ./20 .+ 0.02,
            color = rand(RGBf0, 20)
        )
    end

    @cell "Marker offset" [scatter, marker_offset] begin
        scene = Scene(@resolution)
        points = Point2f0[(0,0), (1,1), (2,2)]
        offset = rand(Point2f0, 3)./5
        scatter!(scene, points)
        scatter!(scene, points, marker_offset = offset, color = :red)
    end

    @cell "colormaps" [image, translate, colormap, colorbrewer, meta] begin
        h = 0.0
        offset = 0.1
        scene = Scene()
        cam2d!(scene)
        plot = map(AbstractPlotting.colorbrewer_names) do cmap
            global h
            c = to_colormap(cmap)
            cbar = image!(
                scene,
                range(0, stop = 10, length = length(c)),
                range(0, stop = 1, length = length(c)),
                reshape(c, (1, length(c))),
                show_axis = false
            )[end]
            text!(
                scene,
                string(cmap, ":"),
                position = Point2f0(-0.1, 0.5 + h),
                align = (:right, :center),
                show_axis = false,
                textsize = 0.4
            )
            translate!(cbar, 0, h, 0)
            h -= (1 + offset)
        end
        scene
    end

    @cell "Available markers" [annotations, markers, meta] begin
        using GeometryTypes
        scene = Scene()
        cam2d!(scene)
        marker = collect(AbstractPlotting._marker_map)
        positions = Point2f0.(0, 1:length(marker))
        scatter!(
            scene,
            positions,
            marker = last.(marker),
            markersize = 0.8,
            raw = true,
            marker_offset = Vec2f0(0.5, -0.4)
        )
        annotations!(
            scene,
            string.(":", first.(marker)),
            positions,
            align = (:right, :center),
            textsize = 0.4,
            raw = true
        )
    end

    #TODO: this doesn't work starting on step 9
    # @cell "Theming" [theme, scatter, surface, set_theme, stepper] begin
    #     new_theme = Theme(
    #         resolution = (500, 500),
    #         linewidth = 3,
    #         colormap = :RdYlGn,
    #         color = :red,
    #         scatter = Theme(
    #             marker = '⊝',
    #             markersize = 0.03,
    #             strokecolor = :black,
    #             strokewidth = 0.1,
    #         ),
    #     )
    #     AbstractPlotting.set_theme!(new_theme)
    #     scene = scatter(rand(100), rand(100))
    #     st = Stepper(scene, @outputfile)
    #     step!(st)
    #     new_theme[:color] = :blue
    #     step!(st)
    #     new_theme[:scatter, :marker] = '◍'
    #     step!(st)
    #     new_theme[:scatter, :markersize] = 0.05
    #     step!(st)
    #     new_theme[:scatter, :strokewidth] = 0.1
    #     step!(st)
    #     new_theme[:scatter, :strokecolor] = :green
    #     step!(st)
    #     empty!(scene)
    #     scene = scatter!(rand(100), rand(100))
    #     step!(st)
    #     scene[end][:marker] = 'π'
    #     step!(st)
    #
    #     r = range(-0.5pi, stop = pi + pi/4, length = 100)
    #     AbstractPlotting.set_theme!(new_theme)
    #     empty!(scene)
    #     scene = surface!(r, r, (x, y)-> sin(2x) + cos(2y))
    #     step!(st)
    #     scene[end][:colormap] = :PuOr
    #     step!(st)
    #     surface!(r + 2pi - pi/4, r, (x, y)-> sin(2x) + cos(2y))
    #     step!(st)
    #     AbstractPlotting.set_theme!(resolution = (500, 500))
    #     empty!(scene)
    #     surface!(r + 2pi - pi/4, r, (x, y)-> sin(2x) + cos(2y))
    #     step!(st)
    #     st
    # end

    @cell "Axis theming" [stepper, axis, lines, stepper] begin
        using GeometryTypes
        scene = Scene()
        points = decompose(Point2f0, Circle(Point2f0(10), 10f0), 9)
        lines!(
            scene,
            points,
            linewidth = 8,
            color = :black
        )

        axis = scene[Axis] # get axis
        scene

        st = Stepper(scene, @outputfile)
        step!(st)
        axis[:frame][:linewidth] = 5
        step!(st)
        axis[:grid][:linewidth] = (1, 5)
        step!(st)
        axis[:grid][:linecolor] = ((:red, 0.3), (:blue, 0.5))
        step!(st)
        axis[:names][:axisnames] = ("x", "y   ")
        step!(st)
        axis[:ticks][:title_gap] = 1
        step!(st)
        axis[:names][:rotation] = (0.0, -3/8*pi)
        step!(st)
        axis[:names][:textcolor] = ((:red, 1.0), (:blue, 1.0))
        step!(st)
        axis[:ticks][:font] = ("Dejavu Sans", "Helvetica")
        step!(st)
        axis[:ticks][:rotation] = (0.0, -pi/2)
        step!(st)
        axis[:ticks][:textsize] = (3, 7)
        step!(st)
        axis[:ticks][:gap] = 5
        step!(st)
    end

    @cell "Stepper demo" [stepper, text, annotation] begin

        function stepper_demo()
            scene = Scene()
            pos = (50, 50)
            steps = ["Step 1", "Step 2", "Step 3"]
            colors = AbstractPlotting.to_colormap(:Set1, length(steps))
            lines!(scene, Rect(0,0,500,500), linewidth = 0.0001)
            # initialize the stepper and give it an output destination
            st = Stepper(scene, @outputfile)

            for i = 1:length(steps)
                text!(
                    scene,
                    steps[i],
                    position = pos,
                    align = (:left, :bottom),
                    textsize = 100,
                    font = "Blackchancery",
                    color = colors[i],
                    scale_plot = false
                )
                pos = pos .+ 100
                step!(st) # saves the step and increments the step by one
            end
            st
        end
        stepper_demo()
    end
end
#
# using Makie
#
# scene = scatter(rand(100), rand(100))
# empty!(scene)
# scene = scatter!(rand(100), rand(100))
#
# r = range(-0.5pi, stop = pi + pi/4, length = 100)
# empty!(scene)
# scene = surface!(r, r, (x, y)-> sin(2x) + cos(2y))
# scene.camera_controls[]
#
# scene[end][:colormap] = :PuOr
# surface!(r + 2pi - pi/4, r, (x, y)-> sin(2x) + cos(2y))
# AbstractPlotting.set_theme!(resolution = (500, 500))
# empty!(scene)
# surface!(r + 2pi - pi/4, r, (x, y)-> sin(2x) + cos(2y))


text:
# Autogenerated wrapper script for tree_sitter_highlight_jll for powerpc64le-linux-gnu
export libtree_sitter_highlight

JLLWrappers.@generate_wrapper_header("tree_sitter_highlight")
JLLWrappers.@declare_library_product(libtree_sitter_highlight, "libtree_sitter_highlight.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libtree_sitter_highlight,
        "lib/libtree_sitter_highlight.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

using Latte
using FactCheck

#=
facts("Testing GRU layer") do
    net = Net(8; time_steps=5)
    data, data_value   = MemoryDataLayer(net, :data, (2,))
    inputs = rand(Float32, size(data_value)..., 8, 5)
    gru1 = GRULayer(:gru1, net, data, 2)
    init(net)
    for t = 1:5
        data_value[:,:] = inputs[:,:,t]
        forward(net; t=t)
    end
    for t in 2:5
        W_z = get_buffer(net, :gru1W_z, 1)
        U_z = get_buffer(net, :gru1U_z, 1)
        b_z = get_buffer(net, :gru1b_z, 1)

        W_r = get_buffer(net, :gru1W_r, 1)
        U_r = get_buffer(net, :gru1U_r, 1)
        b_r = get_buffer(net, :gru1b_r, 1)

        W_h = get_buffer(net, :gru1W_h, 1)
        U_h = get_buffer(net, :gru1U_h, 1)
        b_h = get_buffer(net, :gru1b_h, 1)

        x_z = W_z' * inputs[:,:,t] .+ b_z[:]
        x_r = W_r' * inputs[:,:,t] .+ b_r[:]
        x_h = W_h' * inputs[:,:,t] .+ b_h[:]

        h_tm1 = get_buffer(net, :gru1value, t-1)
        @pending h_tm1 --> get_buffer(net, :gru1inputs2, t)
        sigmoid(x) = 1.0f0 ./ (1.0f0 .+ exp(-x))
        z = sigmoid(x_z .+ (U_z' * h_tm1))
        r = sigmoid(x_r .+ (U_r' * h_tm1))
        hh = tanh(x_h .+ (U_h' * (r .* h_tm1)))
        h = z .* h_tm1 + (1-z) .* hh
        eps = 1e-2
        @pending all(-eps .< get_buffer(net, :gru1value, t) - h .< eps) --> true
    end
    backward(net)
end

FactCheck.exitstatus()
=#


text:
struct EnsembleTestSolution{T,N,S} <: AbstractEnsembleSolution{T,N}
  u::S
  errors::Dict{Symbol,Vector{T}}
  weak_errors::Dict{Symbol,T}
  error_means::Dict{Symbol,T}
  error_medians::Dict{Symbol,T}
  elapsedTime::Float64
  converged::Bool
end
function EnsembleTestSolution(sim::AbstractEnsembleSolution{T,N},errors,weak_errors,error_means,error_medians,elapsedTime,converged) where {T,N}
  EnsembleTestSolution{T,N,typeof(sim.u)}(sim.u,errors,weak_errors,error_means,error_medians,sim.elapsedTime,sim.converged)
end
function EnsembleTestSolution(u,errors,weak_errors,error_means,error_medians,elapsedTime,converged)
  EnsembleTestSolution(EnsembleSolution(u,elapsedTime,converged),errors,weak_errors,error_means,error_medians,elapsedTime,converged)
end

struct EnsembleSolution{T,N,S} <: AbstractEnsembleSolution{T,N}
  u::S
  elapsedTime::Float64
  converged::Bool
end
EnsembleSolution(sim, dims::NTuple{N},elapsedTime,converged) where {N} =
                  EnsembleSolution{eltype(eltype(sim)), N, typeof(sim)}(sim,elapsedTime,converged)
EnsembleSolution(sim,elapsedTime,converged) =
             EnsembleSolution(sim, (length(sim),),elapsedTime,converged) # Vector of some type which is not an array
EnsembleSolution(sim::T,elapsedTime,converged) where T <: AbstractVector{T2} where T2 <: AbstractArray =
             EnsembleSolution(sim, (size(sim[1])..., length(sim)),elapsedTime,converged) # Requires `size` defined on `sim`

struct EnsembleSummary{T,N,Tt,S,S2,S3,S4} <: AbstractEnsembleSolution{T,N}
  t::Tt
  u::S
  v::S2
  qlow::S3
  qhigh::S4
  num_monte::Int
  elapsedTime::Float64
  converged::Bool
end

function calculate_ensemble_errors(sim::AbstractEnsembleSolution;kwargs...)
  calculate_ensemble_errors(sim.u;elapsedTime=sim.elapsedTime,converged=sim.converged,kwargs...)
end

function calculate_ensemble_errors(u;elapsedTime=0.0,converged=false,
                                weak_timeseries_errors=false,weak_dense_errors=false)
  errors = Dict{Symbol,Vector{eltype(u[1].u[1])}}() #Should add type information
  error_means  = Dict{Symbol,eltype(u[1].u[1])}()
  error_medians= Dict{Symbol,eltype(u[1].u[1])}()
  for k in keys(u[1].errors)
    errors[k] = [sol.errors[k] for sol in u]
    error_means[k] = mean(errors[k])
    error_medians[k]=median(errors[k])
  end
  # Now Calculate Weak Errors
  weak_errors =  Dict{Symbol,eltype(u[1].u[1])}()
  # Final
  m_final = mean([s[end] for s in u])
  m_final_analytic = mean([s.u_analytic[end] for s in u])
  res = norm(m_final - m_final_analytic)
  weak_errors[:weak_final] = res
  if weak_timeseries_errors
    ts_weak_errors = [mean([u[j][i] - u[j].u_analytic[i] for j in 1:length(u)]) for i in 1:length(u[1])]
    ts_l2_errors = [sqrt.(sum(abs2,err)/length(err)) for err in ts_weak_errors]
    l2_tmp = sqrt(sum(abs2,ts_l2_errors)/length(ts_l2_errors))
    max_tmp = maximum([maximum(abs.(err)) for err in ts_weak_errors])
    weak_errors[:weak_l2] = l2_tmp
    weak_errors[:weak_l∞] = max_tmp
  end
  if weak_dense_errors
    densetimes = collect(range(u[1].t[1], stop=u[1].t[end], length=100))
    u_analytic = [[sol.prob.f(Val{:analytic},sol.prob.u0,sol.prob.p,densetimes[i],sol.W(densetimes[i])[1]) for i in eachindex(densetimes)] for sol in u]
    dense_weak_errors = [mean([u[j](densetimes)[i] - u_analytic[j][i] for j in 1:length(u)]) for i in eachindex(densetimes)]
    dense_L2_errors = [sqrt.(sum(abs2,err)/length(err)) for err in dense_weak_errors]
    L2_tmp = sqrt(sum(abs2,dense_L2_errors)/length(dense_L2_errors))
    max_tmp = maximum([maximum(abs.(err)) for err in dense_weak_errors])
    weak_errors[:weak_L2] = L2_tmp
    weak_errors[:weak_L∞] = max_tmp
  end
  return EnsembleTestSolution(u,errors,weak_errors,error_means,error_medians,elapsedTime,converged)
end

### Displays

Base.summary(A::AbstractEnsembleSolution) = string("EnsembleSolution Solution of length ",length(A.u)," with uType:\n",eltype(A.u))
function Base.show(io::IO, A::AbstractEnsembleSolution)
  print(io,summary(A))
end
function Base.show(io::IO, m::MIME"text/plain", A::AbstractEnsembleSolution)
  print(io,summary(A))
end

### Plot Recipes

@recipe function f(sim::AbstractEnsembleSolution;
                   zcolors = typeof(sim.u)<:AbstractArray ? fill(nothing, length(sim.u)) : nothing,
                   idxs = typeof(sim.u)<:AbstractArray ? eachindex(sim.u) : 1)
  for i in idxs
    size(sim[i].u, 1) == 0 && continue
    @series begin
      legend := false
      xlims --> (-Inf,Inf)
      ylims --> (-Inf,Inf)
      zlims --> (-Inf,Inf)
      zcolor --> zcolors[i]
      sim[i]
    end
  end
end

@recipe function f(sim::EnsembleSummary;
                   idxs= typeof(sim.u[1])<:AbstractArray ? eachindex(sim.u[1]) : 1,
                   error_style=:ribbon,ci_type=:quantile)
  if typeof(sim.u[1])<:AbstractArray
    u = vecarr_to_vectors(sim.u)
  else
    u = [sim.u.u]
  end
  if ci_type == :SEM
    if typeof(sim.u[1])<:AbstractArray
      ci_low = vecarr_to_vectors(VectorOfArray([sqrt.(sim.v[i]/sim.num_monte).*1.96 for i in 1:length(sim.v)]))
      ci_high = ci_low
    else
      ci_low = [[sqrt(sim.v[i]/length(sim.num_monte)).*1.96 for i in 1:length(sim.v)]]
      ci_high = ci_low
    end
  elseif ci_type == :quantile
    if typeof(sim.u[1])<:AbstractArray
      ci_low = u - vecarr_to_vectors(sim.qlow)
      ci_high = vecarr_to_vectors(sim.qhigh) - u
    else
      ci_low = [u[1] - sim.qlow.u]
      ci_high = [sim.qhigh.u - u[1]]
    end
  else
    error("ci_type choice not valid. Must be :variance or :quantile")
  end
  for i in idxs
    @series begin
      legend --> false
      lw --> 3
      fillalpha --> 0.2
      if error_style == :ribbon
        ribbon --> (ci_low[i],ci_high[i])
      elseif error_style == :bars
        yerror --> (ci_low[i],ci_high[i])
      elseif error_style == :none
        nothing
      else
        error("error_style not recognized")
      end
      x := sim.t
      y := u[i]
    end
  end
end