text:
import Base: convert, getindex, setindex!, similar, in
using NullableArrays: NullableArray

## Constructors and converters
## (special methods for AbstractArray{Nullable}, to avoid wrapping Nullable inside Nullable)

NullableCategoricalArray{T, N}(::Type{Nullable{T}}, dims::NTuple{N,Int}; ordered=false) =
    NullableCategoricalArray{T, N}(zeros(DefaultRefType, dims), CategoricalPool(ordered))
NullableCategoricalArray{T}(::Type{Nullable{T}}, dims::Int...; ordered=false) =
    NullableCategoricalArray(T, dims; ordered=ordered)

@compat (::Type{NullableCategoricalArray{Nullable{T}, N, R}}){T, N, R}(dims::NTuple{N,Int};
                                                                       ordered=false) =
    NullableCategoricalArray(zeros(R, dims), CategoricalPool{T, R}(ordered))
@compat (::Type{NullableCategoricalArray{Nullable{T}, N}}){T, N}(dims::NTuple{N,Int};
                                                                 ordered=false) =
    NullableCategoricalArray{T}(dims; ordered=ordered)
@compat (::Type{NullableCategoricalArray{Nullable{T}}}){T}(m::Int;
                                                           ordered=false) =
    NullableCategoricalArray{T}((m,); ordered=ordered)
@compat (::Type{NullableCategoricalArray{Nullable{T}}}){T}(m::Int, n::Int;
                                                           ordered=false) =
    NullableCategoricalArray{T}((m, n); ordered=ordered)
@compat (::Type{NullableCategoricalArray{Nullable{T}}}){T}(m::Int, n::Int, o::Int;
                                                           ordered=false) =
    NullableCategoricalArray{T}((m, n, o); ordered=ordered)

@compat (::Type{NullableCategoricalArray{Nullable{CategoricalValue{T, R}}, N, R}}){T, N, R}(dims::NTuple{N,Int};
                                                                                ordered=false) =
    NullableCategoricalArray{T, N, R}(dims; ordered=ordered)
@compat (::Type{NullableCategoricalArray{Nullable{CategoricalValue{T}}, N, R}}){T, N, R}(dims::NTuple{N,Int};
                                                                                ordered=false) =
    NullableCategoricalArray{T, N, R}(dims; ordered=ordered)
@compat (::Type{NullableCategoricalArray{Nullable{CategoricalValue{T, R}}, N}}){T, N, R}(dims::NTuple{N,Int};
                                                                                ordered=false) =
    NullableCategoricalArray{T, N, R}(dims; ordered=ordered)
@compat (::Type{NullableCategoricalArray{Nullable{CategoricalValue{T}}, N}}){T, N}(dims::NTuple{N,Int};
                                                                                ordered=false) =
    NullableCategoricalArray{T, N}(dims; ordered=ordered)
# @compat (::Type{NullableCategoricalArray{Nullable{CategoricalValue}, N}}){N}(dims::NTuple{N,Int};
#                                                                               ordered=false) =
#     NullableCategoricalArray{String, N}(dims; ordered=ordered)
# @compat (::Type{NullableCategoricalArray{Nullable{CategoricalValue}}}){N}(dims::NTuple{N,Int};
#                                                                            ordered=false) =
#     NullableCategoricalArray{String, N}(dims; ordered=ordered)

@compat (::Type{NullableCategoricalVector{Nullable{T}}}){T}(m::Int; ordered=false) =
    NullableCategoricalArray{T}((n,); ordered=ordered)
@compat (::Type{NullableCategoricalMatrix{Nullable{T}}}){T}(m::Int, n::Int; ordered=false) =
    NullableCategoricalArray{T}((m, n); ordered=ordered)

@compat (::Type{NullableCategoricalArray}){T<:Nullable}(A::AbstractArray{T};
                                                        ordered=_isordered(A)) =
    NullableCategoricalArray{eltype(T)}(A, ordered=ordered)
@compat (::Type{NullableCategoricalVector}){T<:Nullable}(A::AbstractVector{T};
                                                         ordered=_isordered(A)) =
    NullableCategoricalVector{eltype(T)}(A, ordered=ordered)
@compat (::Type{NullableCategoricalMatrix}){T<:Nullable}(A::AbstractMatrix{T};
                                                         ordered=_isordered(A)) =
    NullableCategoricalMatrix{eltype(T)}(A, ordered=ordered)

"""
    NullableCategoricalArray(A::AbstractArray, missing::AbstractArray{Bool};
                             ordered::Bool=false)

Similar to definition above, but marks as null entries for which the corresponding entry
in `missing` is `true`.
"""
function NullableCategoricalArray{T, N}(A::AbstractArray{T, N},
                                        missing::AbstractArray{Bool, N};
                                        ordered=false)
    res = NullableCategoricalArray{T, N}(size(A); ordered=ordered)
    @inbounds for (i, x, m) in zip(eachindex(res), A, missing)
        res[i] = ifelse(m, Nullable{T}(), x)
    end

    if method_exists(isless, (T, T))
        levels!(res, sort(levels(res)))
    end

    res
end

"""
    NullableCategoricalVector(A::AbstractVector, missing::AbstractVector{Bool};
                              ordered::Bool=false)

Similar to definition above, but marks as null entries for which the corresponding entry
in `missing` is `true`.
"""
NullableCategoricalVector{T}(A::AbstractVector{T},
                             missing::AbstractVector{Bool};
                             ordered=false) =
    NullableCategoricalArray(A, missing; ordered=ordered)

"""
    NullableCategoricalMatrix(A::AbstractMatrix, missing::AbstractMatrix{Bool};
                              ordered::Bool=false)

Similar to definition above, but marks as null entries for which the corresponding entry
in `missing` is `true`.
"""
NullableCategoricalMatrix{T}(A::AbstractMatrix{T},
                             missing::AbstractMatrix{Bool};
                             ordered=false) =
    NullableCategoricalArray(A, missing; ordered=ordered)

@inline function getindex(A::NullableCategoricalArray, I...)
    @boundscheck checkbounds(A, I...)
    # Let Array indexing code handle everything
    @inbounds r = A.refs[I...]

    if isa(r, Array)
        return ordered!(arraytype(A)(r, deepcopy(A.pool)),
                        isordered(A))
    else
        S = eltype(eltype(A))
        if r > 0
            @inbounds return Nullable{S}(A.pool[r])
        else
            return Nullable{S}()
        end
    end
end

@inline function setindex!(A::NullableCategoricalArray, v::Nullable, I::Real...)
    @boundscheck checkbounds(A, I...)
    if isnull(v)
        @inbounds A.refs[I...] = 0
    else
        @inbounds A[I...] = get(v)
    end
end

levels!(A::NullableCategoricalArray, newlevels::Vector; nullok=false) = _levels!(A, newlevels, nullok=nullok)

droplevels!(A::NullableCategoricalArray) = levels!(A, _unique(Array, A.refs, A.pool))

unique{T}(A::NullableCategoricalArray{T}) = _unique(NullableArray{T}, A.refs, A.pool)

function in{T, N, R}(x::Nullable, y::NullableCategoricalArray{T, N, R})
    ref = get(y.pool, get(x), zero(R))
    ref != 0 ? ref in y.refs : false
end

function in{S<:CategoricalValue, T, N, R}(x::Nullable{S}, y::NullableCategoricalArray{T, N, R})
    v = get(x)
    if v.pool === y.pool
        return v.level in y.refs
    else
        ref = get(y.pool, index(v.pool)[v.level], zero(R))
        return ref != 0 ? ref in y.refs : false
    end
end

in{T, N, R}(x::Any, y::NullableCategoricalArray{T, N, R}) = false
in{T, N, R}(x::CategoricalValue, y::NullableCategoricalArray{T, N, R}) = false


text:
module BenchSum

using BenchmarkTools
using BlockArrays: mortar
using FLoops
using FillArrays: Zeros

function sum_for(xs)
    acc = 0.0
    for x in xs
        acc += x
    end
    return acc
end

function sum_floop(xs)
    @floop begin
        acc = 0.0
        for x in xs
            acc += x
        end
    end
    return acc
end

const SUITE = BenchmarkGroup()

floats = randn(1000)
dataset = [
    "Vector" => floats,
    "filter" => Iterators.filter(!ismissing, ifelse.(floats .> 2, missing, floats)),
    "flatten" => Iterators.flatten((floats, (1, 2), 3:4, 5:0.2:6, Zeros(1000))),
    "BlockVector" => mortar([floats, floats]),
]

for (label, xs) in dataset
    @assert sum_for(xs) == sum_floop(xs)
    s1 = SUITE[:label=>label] = BenchmarkGroup()
    s1[:impl=>"for"] = @benchmarkable sum_for($xs)
    s1[:impl=>"floop"] = @benchmarkable sum_floop($xs)
end

end  # module
BenchSum.SUITE

ext: jl


text:
"""
A test for gage(RotationAngle, ...).

Arguments
----------
vector_1:Array
    input 2D vector 1

vecto_2:Array
    input 2D vector 2

solution:AbstractFloat
    answer key

msg="":AbstractString
    optional test message
"""
function test_gage(RotationAngle, vector_1::Array, vector_2::Array, solution::AbstractFloat, msg::AbstractString="")
    print_dashed_line(80)
    print_with_color(:blue, "Test gage(RotationAngle, ...)\n\n")
    if msg != ""
        print_with_color(:blue, "$(msg)\n\n")
    end

    @time answer = gage(RotationAngle,  vector_1, vector_2)
    @test_approx_eq answer solution

    print_with_color(:green, "VERIFIED!\n")
    print_dashed_line(80)
    return true
end

test_gage(RotationAngle, [1., 0.],[0., 1.],pi/2,"test 1")
test_gage(RotationAngle, [0., 0.],[0., 1.],0.,"test 2")
test_gage(RotationAngle, [1., 0.],[0., 0.],0.,"test 3")
test_gage(RotationAngle, [0., 0.],[0., 0.],0.,"test 4")
test_gage(RotationAngle, [-0.5005026560060004,-0.08853187256313988],[-5.908280911892572,-1.04509170227587],0.,"test 5")
test_gage(RotationAngle, -[-0.5005026560060004,-0.08853187256313988],[-5.908280911892572,-1.04509170227587],pi*1.,"test 6")
test_gage(RotationAngle, [1./sqrt(2),1./sqrt(2),1.],[-1./sqrt(2),-1./sqrt(2),1.],pi/2.,"test 7")


text:
using Catalyst
using Distributions: Normal, logpdf, pdf, MvNormal
using StatsBase
using PWS.Trajectories
using StaticArrays

corr_time_s = 100
mean_s = 100
corr_time_ratio = 5

λ = 1 / corr_time_s
κ = mean_s * λ
μ = corr_time_ratio / corr_time_s
ρ = μ
mean_x = mean_s

# see Tostevin, ten Wolde, eq. 27
sigma_squared_ss = mean_s
sigma_squared_sx = ρ * mean_s / (λ + μ)
sigma_squared_xx = mean_x * (1 + ρ / (λ + μ))

joint_stationary = MvNormal([mean_s, mean_x], [sigma_squared_ss sigma_squared_sx; sigma_squared_sx sigma_squared_xx])

signal_stationary = Normal(mean_s, sqrt(sigma_squared_ss))
response_stationary = Normal(mean_x, sqrt(sigma_squared_xx))

entropy(signal_stationary) + entropy(response_stationary) - entropy(joint_stationary)

sn = @reaction_network begin
    κ, ∅ --> S
    λ, S --> ∅
end κ λ

rn = @reaction_network begin
    ρ, S --> X + S
    μ, X --> ∅
end ρ μ

gen = Trajectories.configuration_generator(sn, rn, [κ, λ], [ρ, μ], mean_s, mean_x)
jump_problem = gen.joint_j_problem

using Plots

grid = vcat.(50.0:150.0, ( 50.0:150.0 )')
val = map(x -> pdf(gen.p0_dist, x), grid)
gr()
heatmap(50.0:150.0, 50.0:150.0, val, c=:bilbao)


text:
using GAIO
using StaticArrays
using Test

@testset "exported functionality" begin
    @testset "basics" begin
        center = SVector(0.0, 0.1)
        radius = SVector(10.0, 10.0)
        box = Box(center, radius)
        @test box.center == center
        @test box.radius == radius
    end
    @testset "types" begin
        center = SVector(0, 0, 1)
        radius = SVector(1.0, 0.1, 1.0)
        box = Box(center, radius)
        @test typeof(box.center) <: typeof(box.radius)
        @test typeof(box.radius) <: typeof(box.center)
        @test !(typeof(box.center) <: typeof(center))
    end
    @testset "containment" begin
        center = SVector(0.0, 0.0, 0.0)
        radius = SVector(1.0, 1.0, 1.0)
        box = Box(center, radius)
        inside = SVector(0.5, 0.5, 0.5)
        left = SVector(-1.0, -1.0, -1.0)
        right = SVector(1.0, 1.0, 1.0)
        on_boundary_left = SVector(0.0, 0.0, -1.0)
        on_boundary_right = SVector(0.0, 1.0, 0.0)
        outside_left = SVector(0.0, 0.0, -2.0)
        outside_right = SVector(0.0, 2.0, 0.0)
        @test inside ∈ box
        @test box.center ∈ box
        #boxes are half open to the right
        @test left ∈ box
        @test right ∉ box
        @test on_boundary_left ∈ box
        @test on_boundary_right ∉ box
        @test outside_left ∉ box
        @test outside_right ∉ box
    end
    @testset "non matching dimensions" begin
        center = SVector(0.0, 0.0, 0.0)
        radius = SVector(1.0, 1.0)
        @test_throws Exception Box(center, radius)
    end
    @testset "negative radii" begin
        center = SVector(0.0, 0.0)
        radius = SVector(1.0, -1.0)
        @test_throws Exception Box(center, radius)
    end
end
@testset "internal functionality" begin
    box = Box(SVector(0.0, 0.0), SVector(1.0, 1.0))
    @testset "integer point in box" begin
        point_int_outside = SVector(2, 2)
        point_int_inside = SVector(0, 0)
        @test point_int_inside ∈ box
        @test point_int_outside ∉ box
    end
    @test_throws DimensionMismatch SVector(0.0, 0.0, 0.0) ∈ box
end


text:
using Test, Pkg
using MPI

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    using CuArrays
    ArrayType = CuArray
else
    ArrayType = Array
end

MPI.Init()

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)
root = 0
isroot = rank == root

for T in Base.uniontypes(MPI.MPIDatatype)

    # Allocating version
    A = isroot ? ArrayType{T}(1:size) : ArrayType{T}(undef, 1)
    B = MPI.Scatter(A, 1, root, comm)
    @test B isa ArrayType{T}
    @test Array(B)[1] == T(rank+1)

    # Non Allocating version
    B = ArrayType{T}(undef, 1)
    MPI.Scatter!(A, B, 1, root, comm)
    @test Array(B)[1] == T(rank+1)

    # In place version
    B = isroot ? copy(A) : ArrayType{T}(undef, 1)
    if root == MPI.Comm_rank(comm)
        MPI.Scatter!(B, nothing, 1, root, comm)
    else
        MPI.Scatter!(nothing, B, 1, root, comm)
    end
    @test Array(B)[1] == T(rank+1)

    # Test throwing
    B = ArrayType{T}(undef, 0)
    @test_throws AssertionError MPI.Scatter!(A, B, 1, root, comm)
end

MPI.Finalize()
@test MPI.Finalized()


text:
using JLD
# using Pandas
using Base.Test
import Base.convert
using MAT
using ForwardDiff
using Polynomials
import ForwardDiff.DiffBase
using JSON
# using PyPlot

# local
# @everywhere include("/home/genis/cluster_home/model_fitting_julia/gradient_utils.jl")
# @everywhere include("/home/genis/cluster_home/model_fitting_julia/functions_julia.jl")
# path="/home/genis/cluster_home/model_fitting_julia/"

#portatil

path="/home/emma/fit_behaviour_julia/"
include(path*"functions/functions_julia.jl")
include(path*"functions/utils.jl")
include(path*"functions/analysis.jl")
include(path*"functions/parameter_functions.jl")
include(path*"functions/simulations.jl")

SUBJECTS=[24,25,36,37]

model_fit="StLogTimePi"

Tframe=0.5

model_f=Dict(:c2_f=>c2_const,:c4_f=>c4_const,:c6_f=>c6_zeros,:bias_d_f=>bias_d_const,
        :hbias_d_f=>history_bias_d_const,:st_f=>st_logistic_time,:x0_f=>initial_condition_bias_hbias,
        :ro_f=>read_out_perfect,:lambda=>0)

###### PARAM simulations ######

tau=1.0

#############################

subject=1

for subject in SUBJECTS
#subject=4
    path_load=path*"results/summary/model_fit"*model_fit*"/"
    filename_load=path_load*"rat"*string(subject)*"_05sec.jld"
    data_fit=load(filename_load)

    iminimum=find(x->x==minimum(data_fit["LL_training"]),data_fit["LL_training"])
    param_best=data_fit["param_all"][iminimum][1]
#    param_best = Dict("c2"=> [10, 0.0],"st"=> [1, 0.0, 0.0, 0.0],"sigma_a"=> [1.0],"x0"=> [0.0, 0.0],"bias_d"=> [0.0],"c4"=>[1],"hbias_d"=> [0.0])
	#println("param_ebst", param_best)
    PR_fit=data_fit["PR_training"][iminimum][1]
    path_data=path*"scripts_rats/data/"

    # data = Dict()
    open(path_data*"processed_data_rat"*string(subject)*"_dataset2_05.json", "r") do f

        global data
        dicttxt = readstring(f)  # file information to string
        data=JSON.parse(dicttxt)  # parse and transform data
    end
    n_stim_10 = data["n_stim_10"]
    n_stim_20 = data["n_stim_20"]

	#println("nstim_10! ",n_stim_10)
    #Nstim=100
    # stim=data["stim"][1:Nstim,:]
    # choices=data["choices"][1:Nstim]
    # i_sigma=data["i_sigma"][1:Nstim]

    stim=data["stim_10"]
    stim = permutedims(reshape(stim,(10,n_stim_10)))
    println("1stim ",stim[1,:])
    choices=data["choices_10"]
    i_sigmas = data["sigma_10"]
    Nframes=10

    tau=1.0
    Tframe=0.5
    T=Tframe*Nframes
    dt=0.04
    #NT=T/dt
    #println("T",T)
    #println(dt)
    #println(param_best)
    # stim=data["stim"]
    # choices=data["choices"]
    # i_sigmas=data["i_sigma"]
    Nstim=100
    Ntrials=500
    PR_sim=zeros(Nstim)
    for istim in 1:Nstim
        internal_noise_sigma=param_best["sigma_a"][1]*randn(Ntrials,50)
        stim_rep=repeat(transpose(stim[istim,:]),Ntrials)
        d,xfinal=simulation_Nframe_PI(param_best,stim_rep,internal_noise_sigma,T,tau,model_f)
        PR_sim[istim]=mean( (d+1)/2 )
    end

    save=Dict("PR_sim"=>PR_sim,"PR_fit"=>PR_fit)

    path_save=path*"figures/"
    filename_save=path_save*"check_PR_fit_PR_sim"*model_fit*"_subject"*string(subject)*"_05sec.jld"
    JLD.save(filename_save,save)
    json_string = json_string = JSON.json(save)
	open(path_save*"check_PR_fit_PR_sim"*model_fit*"_subject"*string(subject)*"_05sec.json","w") do f
		write(f, json_string)
	print("done")
	end

end

max_stars_repo_name: emmafdz/FinalDegreeProject

max_issues_repo_name: emmafdz/FinalDegreeProject

max_forks_repo_name: emmafdz/FinalDegreeProject


text:
# We have different sources for problem functions:
#  S1 = CEC 2013 competition on large-scale optimization
#  S2 = JADE paper http://150.214.190.154/EAMHCO/pdf/JADE.pdf
#  S3 = "Test Suite for the Special Issue of Soft Computing on Scalability of
#        Evolutionary Algorithms and other Metaheuristics for Large Scale
#        Continuous Optimization Problems", http://sci2s.ugr.es/eamhco/functions1-19.pdf
# Our primary focus is to implement all the problems from S1 since our
# focus is on large-scale optimization but these problems also can be used
# in lower dimensions.
include(joinpath(dirname(@__FILE__()), "single_objective_base_functions.jl"))
Shekel10 = minimization_problem(shekel10, "Shekel10", (0.0, 10.0), 4, -10.5364)
Shekel7 = minimization_problem(shekel7, "Shekel7", (0.0, 10.0), 4, -10.4029)
Shekel5 = minimization_problem(shekel5, "Shekel5", (0.0, 10.0), 4, -10.1532)
Hartman6 = minimization_problem(hartman6, "Hartman6", (0.0, 1.0), 6, -3.32237)
Hartman3 = minimization_problem(hartman3, "Hartman3", (0.0, 1.0), 3, -3.860038442)

"""
JADE collection of optimization problems.

We skip (for now) `f12` and `f13` in the JADE paper since they are penalized
functions which are quite nonstandard. We also skip `f8` since we are unsure
about its proper implementation.
"""
const JadeFunctionSet = Dict{Int,FunctionBasedProblemFamily}(
    1   => MinimizationProblemFamily(sphere, "Sphere", (-100.0, 100.0), 0.0),
    2   => MinimizationProblemFamily(schwefel2_22,  "Schwefel2.22",  ( -10.0,  10.0), 0.0),
    3   => MinimizationProblemFamily(schwefel1_2,   "Schwefel1.2",   (-100.0, 100.0), 0.0),
    4   => MinimizationProblemFamily(schwefel2_21,  "Schwefel2.21",  (-100.0, 100.0), 0.0),
    5   => MinimizationProblemFamily(rosenbrock,    "Rosenbrock",    ( -30.0,  30.0), 0.0),
    6   => MinimizationProblemFamily(s2_step,       "Step",          (-100.0, 100.0), 0.0),
    7   => MinimizationProblemFamily(noisy_quartic, "Noisy quartic", ( -30.0,  30.0)),
#    8   => MinimizationProblemFamily(schwefel2_26,  "Schwefel2.26",  (-500.0, 500.0)),
    9   => MinimizationProblemFamily(rastrigin,     "Rastrigin",     ( -5.12,  5.12), 0.0),
    10  => MinimizationProblemFamily(ackley,        "Ackley",        ( -32.0,  32.0), 0.0),
    11  => MinimizationProblemFamily(griewank,      "Griewank",      (-600.0, 600.0), 0.0)
)

#####################################################################
# S3 Base functions.
#####################################################################

#####################################################################
# S3 Transformations
#####################################################################

"""
`TransformedProblem{FS, P}` is an abstract class for optimization problems derived
from some original problem of type `P` by introducing a few changes.

The concrete derived types must implement:
* `objfunc()` method
* `fitness_scheme()` and `search_space()`
"""
abstract type TransformedProblem{FS, P <: OptimizationProblem} <: OptimizationProblem{FS} end

orig_problem(tp::TransformedProblem) = tp.orig_problem # default implementation
orig_problem_type(tp::Type{<:TransformedProblem{P}}) where P = P
name(tp::TransformedProblem) = "TransformedProblem("*name(sub_problem(tp))*")"

"""
A `TransformedProblem` subclass that shifts the objective argument and offsets the value:
``f(x - xshift) + fitshift`` → min/max, ``x`` ∈ ``X`` (``X`` stays intact).
"""
struct ShiftedAndBiasedProblem{FS, F, P <: OptimizationProblem} <: TransformedProblem{FS, P}
    orig_problem::P

    xshift::Vector{Float64}
    fitshift::F
    opt_value::F

    function ShiftedAndBiasedProblem(orig_problem::P;
        xshift = nothing,
        fitshift = nothing,
        opt_value = nothing
    ) where P <: OptimizationProblem{FS} where FS
        (xshift === nothing) && (xshift = fill(0.0, numdims(orig_problem)))
        F = fitness_type(FS)
        (fitshift === nothing) && (fitshift = zero(F)) # FIXME zerofitness() to support tuples
        # derive the default optimal value,
        # which would be not correct if `xshift` moves the optimum in/out of the search space
        (opt_value === nothing) && (opt_value = BlackBoxOptim.opt_value(orig_problem) + fitshift)
        new{FS, F, P}(orig_problem, copy(xshift), fitshift, opt_value)
    end
end

name(p::ShiftedAndBiasedProblem) = "ShiftedAndBiased("*name(orig_problem(p))*")"
fitness_scheme(tp::TransformedProblem) = fitness_scheme(orig_problem(tp))
search_space(tp::TransformedProblem) = search_space(orig_problem(tp)) # the space stays the same
opt_value(tp::TransformedProblem) = tp.opt_value

fitness(x, sp::ShiftedAndBiasedProblem) =
    fitness(x - sp.xshift, orig_problem(sp)) + sp.fitshift

#####################################################################
# S1 Base functions. Typically slightly transformed to break symmetry
#   and introduce irregularities.
#####################################################################
s1_sphere = sphere

function s1_elliptic(x)
    xt = t_irreg(x)
    elliptic(xt)
end

function s1_rastrigin(x)
    xt = t_diag(t_asy(t_irreg(x), 0.2), 10)
    rastrigin(xt)
end

function s1_ackley(x)
    xt = t_diag(t_asy(t_irreg(x), 0.2), 10)
    ackley(xt)
end

function s1_schwefel(x)
    xt = t_asy(t_irreg(x), 0.2)
    schwefel1_2(xt)
end

s1_rosenbrock = rosenbrock

#####################################################################
# S1 Transformations
#####################################################################

"""
Transform symmetric `f` into asymmetric objective function.
"""
function t_asy(f, beta)
    D = length(f)
    g = copy(f)
    temp = beta * range(0, stop=1, length=D)
    ind = collect(1:D)[f .> 0]
    t = f[ind] .^ (1 + temp[ind] .* sqrt(f[ind]))
    setindex!(g, t, ind)
    return g
end

"""
Transform `f` into objective function with ill-conditioning effect.
"""
function t_diag(f, alpha)
    D = length(f)
    scales = sqrt(alpha) .^ range(0, stop=1, length=D)
    return scales .* f
end

"""
Transform `f` into objective function with smooth local irregularities.
"""
function t_irreg(f)
    a = 0.1
    g = copy(f)
    indices = collect(1:length(f))

    idxp = indices[f .> 0]
    t = log(f[idxp])/a
    r = exp(t + 0.49*(sin(t) + sin(0.79*t))).^a
    setindex!(g, r, idxp)

    idxn = indices[f .< 0]
    t = log(-f[idxn])/a
    r = -exp(t + 0.49*(sin(0.55*t) + sin(0.31*t))).^a
    setindex!(g, r, idxn)

    return g
end

function xshifted(n, f)
    move = 10.0 * randn(n, 1)
    transformed_f(x) = f(x .- move)
end

function xrotatedandshifted(n, f, shiftAmplitude = 1.0, rotateAmplitude = 1.0)
    shift = shiftAmplitude * randn(n, 1)
    rotmatrix = rotateAmplitude * rand(n, n)
    transformed_f(x) = f(rotmatrix * (x .- shift))
end

const example_problems = Dict{String, Union{OptimizationProblem,FunctionBasedProblemFamily}}(
    "Sphere" => JadeFunctionSet[1],
    "Rosenbrock" => JadeFunctionSet[5],
    "Schwefel2.22" => JadeFunctionSet[2],
    "Schwefel1.2" => JadeFunctionSet[3],
    "Schwefel2.21" => JadeFunctionSet[4],
    "Step" => JadeFunctionSet[6],
    "Rastrigin" => JadeFunctionSet[9],
    "Ackley" => JadeFunctionSet[10],
    "Griewank" => JadeFunctionSet[11],
    "Ellipsoid" => MinimizationProblemFamily(ellipsoid, "Ellipsoid", (-65.536, 65.536), 0.0),
    "Cigar" => MinimizationProblemFamily(cigar, "Cigar", (-100.0, 100.0), 0.0),
    "DeceptiveCuccu2011_15_2" => MinimizationProblemFamily(deceptive_cuccu2011_15_2, "DeceptiveCuccu2011_15_2", (-100.0, 100.0), 0.0),
    "Shekel10" => Shekel10,
    "Shekel7" => Shekel7,
    "Shekel5" => Shekel5,
    "Hartman6" => Hartman6,
    "Hartman3" => Hartman3,
    "Tsallis1996" => MinimizationProblemFamily(energy_tsallis1996, "Tsallis1996", (-100.0, 100.0), 0.0),
)


text:
function Point2ID(X)
    Y = 0
    for ii = 1:length(X)
        Y += X[ii]*2^(ii-1);
    end
    Y
end

function GenNeighbor(X)
    N = length(X)
    Y = zeros(N)
    for ii = 1:N
        temp = copy(X)
        temp[ii] = 1-temp[ii]
        Y[ii] = Point2ID(temp) + 1
    end
    sort!(Y)
    Y
end

function GenCuts(X)
    N = length(X)
    Y = zeros(N)
    for ii = 1:N
        Y[ii] = ii + X[ii]*N
    end
    sort!(Y)
end

function Int2Bin(X,N)
    Y = zeros(N)
    for ii = N:-1:1
        Y[ii] = floor(X/2^(ii-1))
        X = X - Y[ii]*2^(ii-1)
    end
    Y
end
include("FindIntersectingPoint.jl")
include("Vert.jl")
include("PTVer.jl")


text:
export faces

function faces(p::Tetrahedron)
    v = vertices(p)
    return Vec(Triangle(v[1], v[2], v[3]),
               Triangle(v[1], v[2], v[4]),
               Triangle(v[2], v[3], v[4]),
               Triangle(v[3], v[1], v[4]))
end

function faces(p::Hexahedron)
    v = vertices(p)
    return Vec(Quadrilateral(v[1], v[2], v[3], v[4]),
               Quadrilateral(v[5], v[6], v[7], v[8]),
               Quadrilateral(v[1], v[2], v[6], v[5]),
               Quadrilateral(v[2], v[3], v[7], v[6]),
               Quadrilateral(v[3], v[4], v[8], v[7]),
               Quadrilateral(v[4], v[1], v[5], v[8]))
end

function faces(p::QuadraticTetrahedron)
    v = vertices(p)
    return Vec(QuadraticTriangle(v[1], v[2], v[3], v[5],  v[6],  v[7]),
               QuadraticTriangle(v[1], v[2], v[4], v[5],  v[9],  v[8]),
               QuadraticTriangle(v[2], v[3], v[4], v[6], v[10],  v[9]),
               QuadraticTriangle(v[3], v[1], v[4], v[7],  v[8], v[10]))
end

function faces(p::QuadraticHexahedron)
    v = vertices(p)
    return Vec(QuadraticQuadrilateral(v[1], v[2], v[3], v[4], v[ 9], v[10], v[11], v[12]),
               QuadraticQuadrilateral(v[5], v[6], v[7], v[8], v[13], v[14], v[15], v[16]),
               QuadraticQuadrilateral(v[1], v[2], v[6], v[5], v[ 9], v[18], v[13], v[17]),
               QuadraticQuadrilateral(v[2], v[3], v[7], v[6], v[10], v[19], v[14], v[18]),
               QuadraticQuadrilateral(v[3], v[4], v[8], v[7], v[11], v[20], v[15], v[19]),
               QuadraticQuadrilateral(v[4], v[1], v[5], v[8], v[12], v[17], v[16], v[20]))
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using BinaryBuilder

name = "s3gof3r"

version = v"0.5.0"
# Collection of sources required to build pprof
sources = [
    GitSource("https://github.com/rlmcpherson/s3gof3r", "31603a0dc94aefb822bfe2ceea75a6be6013b445"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/s3g/gof3r
export GO15VENDOREXPERIMENT=1
go get github.com/jessevdk/go-flags
go get github.com/rlmcpherson/s3gof3r/gof3r
mkdir -p ${bindir}
go build -o ${bindir}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("gof3r", :gof3r),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; compilers=[:c, :go])

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Autogenerated wrapper script for Luna_jll for i686-w64-mingw32-cxx03
export libluna, luna

using FFTW_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("Luna")
JLLWrappers.@declare_library_product(libluna, "libluna.dll")
JLLWrappers.@declare_executable_product(luna)
function __init__()
    JLLWrappers.@generate_init_header(FFTW_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libluna,
        "bin\\libluna.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        luna,
        "bin\\luna.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Statistics, DataFrames, CSV, HypothesisTests
	using Random, Distributions, Plots
	Random.seed!(0)
end;

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing6.01"

# ╔═╡ 1aa41898-3afa-11eb-2f49-570853ce6435
begin
	data = CSV.read(datadir("machine1.csv"), DataFrame; header=false)[:, 1]
	xBar, n = mean(data), length(data)
	sig = 1.2
	alpha = 0.1
	z = quantile(Normal(),1-alpha/2)

	Text("Calculating formula: $((xBar - z*sig/sqrt(n), xBar + z*sig/sqrt(n)))")
end

# ╔═╡ 07d8c2c0-3b2e-11eb-344d-c14ab2a82d0d
Text("Using confint() function: $(confint(OneSampleZTest(xBar,sig,n),alpha))")

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing6.01"

# ╔═╡ Cell order:
# ╠═ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═1aa41898-3afa-11eb-2f49-570853ce6435
# ╠═07d8c2c0-3b2e-11eb-344d-c14ab2a82d0d
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
for dimnames in [zip([:1, :2, :3], [:z, :y, :x]), zip([:1, :2], [:z, :x])]
    is = broadcast(x -> Symbol(string("i", x)), getindex.(collect(dimnames), 2))
    ist = Meta.parse(string("(", [string(s, ",") for s in is]..., ")"))
    N = Meta.parse(string(length(is)))

    for (idim, dim) in dimnames
        i = Symbol("i", string(dim))
        # velocity-free boundary conditions
        fname = Symbol("dirichlet", string(dim), "!")

        fdh=div(_fd.order,2)
        # mirror ghost cells
        ighost=vcat(
            [[replace(is, i => :($ifd)), replace(is, i => :($(_fd.order+1-ifd)))] for ifd = 1:fdh],
            [[replace(is, i => :(n + $(_fd.order - ifd))), replace(is, i => :(n + $(ifd-1)))] for ifd = 1:fdh]
        )

        isn = replace(is, i => :n)
        is1 = replace(is, i => :1)

        irest = Meta.parse(
            string("(", [string(s, ",") for s in filter(x -> x != i, is)]..., ")"),
        )
        v = Symbol("v", string(dim))
        vs = broadcast(x -> Symbol(string("v", x)), getindex.(collect(dimnames), 2))
        vrest = filter(x -> x != v, vs)

        @eval @parallel_indices($irest, function $fname($v::Data.Array{$N}, $(vrest...), n)
            # along other dimensions velocity grid matches tauii grid
            $((quote $vv[$(is1...)] = 0; $vv[$(isn...)] = 0 end for vv in vrest)...)
            # relative to the tauii grid, velocity at i=2 should be zero, so making use of ghost cells
            $((quote $v[$(ig[1]...)] = -$v[$(ig[2]...)] end for ig in ighost)...)
            return
        end)
    end
end


text:
module SimulateMotif

######################### Dependencies #########################
using Distributions
################################################################

################# Exported methods and types ###################
export single_part_motif,
       k_parts_motif,
       gapped_k_parts_motif,
       mixture_k_parts_motifs,
       mixture_gapped_k_parts_motifs,
       read_fasta_that_has_ground_truth,
       DNAdataset
################################################################

###################### Load files ##############################
include("types.jl")
include("helpers.jl")
include("sample.jl")
include("DNAdataset.jl")

end


text:
# Attributes #
# ========================== #
"""
    AbstractLShapedAttribute

Abstract supertype for attribute objects specific to the L-shaped algorithm.
"""
abstract type AbstractLShapedAttribute <: AbstractStructuredOptimizerAttribute end
"""
    FeasibilityCuts

An optimizer attribute for specifying if feasibility cuts should be used in the L-shaped algorithm. Options are `true` of `false`.
"""
struct FeasibilityCuts <: AbstractLShapedAttribute end
"""
    Regularizer

An optimizer attribute for specifying a regularization procedure to be used in the L-shaped algorithm. Options are:

- [`NoRegularization`](@ref):  L-shaped algorithm (default)
- [`RegularizedDecomposition`](@ref):  Regularized decomposition ?RegularizedDecomposition for parameter descriptions.
- [`TrustRegion`](@ref):  Trust-region ?TrustRegion for parameter descriptions.
- [`LevelSet`](@ref):  Level-set ?LevelSet for parameter descriptions.
"""
struct Regularizer <: AbstractLShapedAttribute end
"""
    Aggregator

An optimizer attribute for specifying an aggregation procedure to be used in the L-shaped algorithm. Options are:

- [`NoAggregation`](@ref):  Multi-cut L-shaped algorithm (default)
- [`PartialAggregation`](@ref):  ?PartialAggregation for parameter descriptions.
- [`FullAggregation`](@ref):  ?FullAggregation for parameter descriptions.
- [`DynamicAggregation`](@ref):  ?DynamicAggregation for parameter descriptions.
- [`ClusterAggregation`](@ref):  ?ClusterAggregation for parameter descriptions.
- [`HybridAggregation`](@ref):  ?HybridAggregation for parameter descriptions.
"""
struct Aggregator <: AbstractLShapedAttribute end
"""
    Consolidator

An optimizer attribute for specifying a consolidation procedure to be used in the L-shaped algorithm. Options are:

- [`NoConsolidation`](@ref)
- [`Consolidation`](@ref)
"""
struct Consolidator <: AbstractLShapedAttribute end
"""
    RegularizationParameter

Abstract supertype for regularization-specific attributes.
"""
abstract type RegularizationParameter <: AbstractLShapedAttribute end
"""
    AggregationParameter

Abstract supertype for aggregation-specific attributes.
"""
abstract type AggregationParameter <: AbstractLShapedAttribute end
"""
    ConsolidationParameter

Abstract supertype for consolidation-specific attributes.
"""
abstract type ConsolidationParameter <: AbstractLShapedAttribute end


text:
const Component = Union{PPCA, FactorAnalysis}

struct FinateMixtureModel{VS<:ValueSupport, D<:Distribution, C<:Component, CT<:Real} <: MultivariateMixture{VS,D}
    components::Vector{C}
    prior::Categorical{CT}

    function FinateMixtureModel{VS, D}(cs::Vector{C}, pri::Categorical{CT}) where {VS, C, D, CT}
        length(cs) == ncategories(pri) ||
            error("The number of components does not match the length of prior.")
        new{VS,D,C,CT}(cs, pri)
    end
end

FinateMixtureModel(components::Vector{C}, p::Vector{CT}) where {C<:Component, CT<:Real} =
    FinateMixtureModel{Continuous, MvNormal}(components, Categorical(p))

function show(io::IO, d::FinateMixtureModel)
    K = ncomponents(d)
    pr = probs(d)
    println(io, "Mixture of $(component_type(d))(K = $K)")
    Ks = min(K, 8)
    for i = 1:Ks
        print(io, "components[$i] (prior = $(pr[i])): ")
        println(io, component(d, i))
    end
    if Ks < K
        println(io, "The rest are omitted ...")
    end
end

# AbstractMixtureModel Interface

component_type(d::FinateMixtureModel) = typeof(first(d.components))
ncomponents(d::FinateMixtureModel) = length(d.components)
components(d::FinateMixtureModel) = d.components
component(d::FinateMixtureModel, k::Int) = d.components[k]
probs(d::FinateMixtureModel) = probs(d.prior)
length(d::FinateMixtureModel) = length(d.components[1])

function mean(d::FinateMixtureModel)
    K = ncomponents(d)
    π = probs(d)
    m = zeros(length(d))
    for i = 1:K
        πᵢ = π[i]
        if πᵢ > 0.0
            c = component(d, i)
            BLAS.axpy!(πᵢ, mean(c), m)
        end
    end
    return m
end


text:
#======================================================

  [NMAP XML DTD](https://nmap.org/book/nmap-dtd.html)

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

"""
    Timestamp
"""
Base.@kwdef mutable struct Timestamp <:Leaf
    time::Int       = 0
    date::DateTime  = unix2datetime(time)
end
Timestamp(unix::Int)    = Timestamp(time = unix)
Timestamp(unix::String) = Timestamp(parse(Int, unix))
function StructTypes.construct(::Type{Timestamp}, unix::T) where T<:Union{String, Int}
    return Timestamp(unix)
end
time(ts::Timestamp) = ts.time
date(ts::Timestamp) = ts.date

"""
    Task

Represents an nmap Task
"""
Base.@kwdef mutable struct Task <:Leaf
    task::String        = ""
    time::Timestamp     = Timestamp()
    extrainfo::String   = ""
end
time(task::Task) = task.time
date(task::Task) = time(task).date

include("xml/script.jl")
include("xml/os.jl")
include("xml/port.jl")
include("xml/ports.jl")
include("xml/host.jl")

"""
    Verbose

Verbosity level
Wraps the verbose xml element
"""
Base.@kwdef mutable struct Verbose <:Leaf
    level::Int = 0
end

"""
    Debugging

Wraps the debugging xml element
"""
Base.@kwdef mutable struct Debugging <:Leaf
    level::Int = 0
end
level(x::T) where T<:Union{Debugging, Verbose} = x.level

"""
    ScanInfo

Wraps the `scaninfo` xml element
"""
Base.@kwdef mutable struct ScanInfo <:Leaf
    numservices::Int    = 0
    protocol::String    = ""
    flags::String       = ""
    services::String    = ""
    type::String        = ""
end
StructTypes.names(::Type{ScanInfo}) = ((:flags, :scanflags),)

"""
    Finished

Wraps the `finished` xml element
"""
Base.@kwdef mutable struct Finished <:Leaf
    time::String        = ""
    timestr::String     = ""
    elapsed::String     = ""
    summary::String     = ""
    exit::String        = ""
    errormsg::String    = ""
end

"""
    Hosts

Wraps the `hosts` xml element
"""
Base.@kwdef mutable struct Hosts <:Leaf
    up::Int     = 0
    down::Int   = 0
    total::Int  = up + down
end

"""
    RunStats

Represents the runstats of a nmap scan
"""
Base.@kwdef mutable struct RunStats <:Leaf
    finished::Finished  = Finished()
    hosts::Hosts        = Hosts()
end

"""
    Target

How a target was specified when passed to nmap, its status and reasoning

Example xml
```xml
<target specification="domain.does.not.exist" status="skipped" reason="invalid"/>
```
"""
Base.@kwdef mutable struct Target <: Marsh
    specification::String   = ""
    status::String          = ""
    reason::String          = ""
end

"""
    Scan

Nmap executed scan
"""
Base.@kwdef mutable struct Scan <:Marsh
    #= File =#
    xml::String
    #= end File =#

    args::String             = ""
    scanner::String          = ""
    startstr::String         = ""
    version::String          = ""
    xmloutputversion::String = ""

    start::Timestamp         = Timestamp()
    verbose::Verbose         = Verbose()
    debugging::Debugging     = Debugging()
    runstats::RunStats       = RunStats()
    scaninfo::ScanInfo       = ScanInfo()
    hosts::Vector{Host}      = Vector{Host}()
    targets::Vector{Target}  = Vector{Target}()
    # Tasks
    taskbegin::Vector{Task}  = Vector{Task}()
    taskend::Vector{Task}    = Vector{Task}()
end
StructTypes.names(::Type{Scan}) = ((:hosts, :host), (:runstats, :runstats))

"""
    Scan(xmldict::XMLDict.XMLDictElement; rawxml=nothing) ::Scan

Build a `Scan` from an input `xml` document

For usage simplicity, this is implemented at constructor level

Examples
```julia
#= Import external scan results =#
scan = NMAP.Scan(read("myscan.xml"))
```
"""
function Scan(xmldict::XMLDict.XMLDictElement; xml::Union{Nothing, String}=nothing) ::Scan
    if isnothing(xml)
        # ~= initial xml string
        xml = XMLDict.dict_xml(XMLDict.xml_dict(xmldict))
    end
    this = Scan(xml = xml)

    outputversion = get(xmldict, :xmloutputversion, "0.0.0")
    if !(outputversion in OUTPUT_VERSIONS)
        @warn "Detected potentially unsupported XML output version $(outputversion), parsing may fail or produce inaccurate results"
    end
    for (k,v) in xmldict
        field = Marshalling.getfield(Scan, k)
        typemap = Dict([name=>type for (name, type) in zip(fieldnames(Scan), Scan.types)])
        if field in fieldnames(Scan)
            Base.setproperty!(this, field, StructTypes.construct(typemap[field], v))
        end
    end
    return this
end

"""
    Scan(xml::String) ::Scan

Build a `Scan` from an input `xml` string

Dispatcher method to allow `scan = NMAP.Scan("<xml>....</xml>")` and preserve the initial `xml` string
"""
function Scan(xml::String) ::Scan
    return Scan(XMLDict.parse_xml(xml), xml=xml)
end

"""
    Scan(data::Vector{UInt8}) ::Scan

Dispatcher method to allow `scan = NMAP.Scan(read("path/to/file.xml"))`
"""
Scan(data::Vector{UInt8}) = Scan(String(data))


text:
"""
    Publisher

A simple wrapper around a ZMQ publisher, but only publishes protobuf messages.

# Construction

    Publisher(context::ZMQ.Context, ipaddr, port; name)

To create a publisher, pass in a `ZMQ.Context`, which allows all related
publisher / subscribers to be collected in a "group." The publisher also
needs to be provided the IPv4 address (either as a string or as a `Sockets.IPv4`
object), and the port (either as an integer or a string).

A name can also be optionally provided via the `name` keyword, which can be used
to provide a helpful description about what the publisher is publishing. It defaults
to "publisher_#" where `#` is an increasing index.

If the port

# Usage
To publish a message, just use the `publish` method on a protobuf type:

    publish(pub::Publisher, proto_msg::ProtoBuf.ProtoType)
"""
struct ZmqPublisher <: Publisher
    socket::ZMQ.Socket
    port::Int64
    ipaddr::Sockets.IPv4
    buffer::IOBuffer
    name::String
    socket_lock::ReentrantLock
    flags::PublisherFlags

    function ZmqPublisher(
        ctx::ZMQ.Context,
        ipaddr::Sockets.IPv4,
        port::Integer;
        name = genpublishername(),
        buffersize = 255,
    )
        local socket
        @catchzmq(
            socket = ZMQ.Socket(ctx, ZMQ.PUB),
            "Could not create socket for publisher $name."
        )
        @catchzmq(
            ZMQ.bind(socket, "tcp://$ipaddr:$port"),
            "Could not bind publisher $name to $(tcpstring(ipaddr, port))"
        )

        @info "Publishing $name on: $(tcpstring(ipaddr, port)), isopen = $(isopen(socket))"
        new(
            socket,
            port,
            ipaddr,
            IOBuffer(zeros(UInt8, buffersize); read = true, write = true),
            name,
            ReentrantLock(),
            PublisherFlags(),
        )
    end
end

function ZmqPublisher(ctx::ZMQ.Context, ipaddr, port::Integer; kwargs...)
    if !(ipaddr isa Sockets.IPv4)
        ipaddr = Sockets.IPv4(ipaddr)
    end
    ZmqPublisher(ctx, ipaddr, port; kwargs...)
end

function ZmqPublisher(ctx::ZMQ.Context, ipaddr, port::AbstractString; kwargs...)
    ZmqPublisher(ctx, ipaddr, parse(Int, port); kwargs...)
end

getcomtype(::ZmqPublisher) = :zmq
function Base.isopen(pub::ZmqPublisher)
    return lock(() -> ZMQ.isopen(pub.socket), pub.socket_lock)
end

function Base.close(pub::ZmqPublisher)
    lock(pub.socket_lock) do
        if isopen(pub.socket)
            @debug "Closing ZmqPublisher: $(getname(pub))"
            ZMQ.close(pub.socket)
        end
    end
end

function forceclose(pub::ZmqPublisher)
    @warn "Force closing ZmqPublisher: $(getname(pub))"
    close(pub.socket)
end

function publish(pub::ZmqPublisher, msg::MercuryMessage)
    pub.flags.has_published = false

    bytes_sent = Int32(0)
    if isopen(pub)
        # Encode the message with protobuf
        msg_size = encode!(msg, pub.buffer)
        getflags(pub).bytespublished = msg_size

        # Create a new message to be sent and copy the encoded protobuf bytes
        zmsg = ZMQ.Message(msg_size)
        copyto!(zmsg, 1, pub.buffer.data, 1, msg_size)

        # Send over ZMQ
        # NOTE: ZMQ will de-allocate the message allocated above, so garbage
        # collection should not be an issue here
        ZMQ.send(pub.socket, zmsg)
        getflags(pub).has_published = true

        # Move to the beginning of the buffer
        seek(pub.buffer, 0)
    end
end

portstring(pub::ZmqPublisher) = tcpstring(pub.ipaddr, pub.port)


text:
using GraphRecipes # a graph vizualization package built ontop of Plots
using Plots
default(size=(1400,1000))

plot(AbstractFloat, method=:tree, fontsize=12, nodeshape=:ellipse)

savefig("types.png")

max_stars_repo_name: csimal/Julia-Unamur

max_issues_repo_name: csimal/Julia-Unamur

max_forks_repo_name: csimal/Julia-Unamur


text:
# This file contains some possible functions that are invariant under conjugation.

"""
    ord(g::SparseMatrixCSC{ℤ₂, Int64}, n::Int64)::Int64

Calculate the order of a given element `g`, with `g^order = identity`. Only works for elements with finite order ("order > 0").
"""
function ord(g::SparseMatrixCSC{ℤ₂, Int64}, n::Int64)::Int64
    i = 1
    I = id(n)
    while true
        if g^i == I
            return i
        end
        i += 1
    end
end

# trace in ℤ₂: LinearAlgebra.tr() works automatically
"""
    tr_of_power(g::SparseMatrixCSC{ℤ₂, Int64}, power::Int64)::Int64

Alias for the trace `tr(g^power)`.
"""
tr_of_power(g::SparseMatrixCSC{ℤ₂, Int64}, power::Int64)::Int64 = tr(g^power)

"""
    rank_of_power_plus_id(g::Matrix{ℤ₂}, power::Int64, n)::Int64

Alias for `rank(g^power + id(n))`.
"""
rank_of_power_plus_id(g::Matrix{ℤ₂}, power::Int64, n::Int64)::Int64 = rank(g^power + id(n))

# symmetry: LinearAlgebra.issymmetric() works automatically

# functions for the generalized Arf invariant:
## 1. Vectors in ℤ₂ⁿ
"""
    ℤ₂ⁿ(n::Int64)::Vector{Vector{ℤ₂}}

Calculate all `2^n` possible vectors in ``\\mathbb{Z}_{2}^n``.
"""
ℤ₂ⁿ(n::Int64)::Vector{Vector{ℤ₂}} = [[ℤ₂(parse(Bool,i)) for i in s] for s in SubString.(bitstring.(0:2^n-1), length(bitstring(1))-n+1)]

## 2. ker(g+I) = Vg
"""
    ker_g_plus_id(g::SparseMatrixCSC{ℤ₂, Int64}, n::Int64)::Vector{Vector{ℤ₂}}

We are interested in ``Vg = \\{\\vec{x} \\in \\mathbb{Z}_2^n \\, | \\, g \\cdot \\vec{x} = \\vec{x}\\} = \\{\\vec{x} \\in \\mathbb{Z}_2^n \\, | \\, (g \\pm \\mathbf{1}) \\cdot \\vec{x} = \\vec{0}\\}``.
This means, ``Vg`` is equivalent to the kernel, ``ker(g\\pm I)`` (``\\pm`` is equivalent in ℤ₂).
"""
function ker_g_plus_id(g::SparseMatrixCSC{ℤ₂, Int64}, n::Int64)::Vector{Vector{ℤ₂}}
    # define the null vector
    x0 = zeros(ℤ₂, n)
    Vg = Vector{ℤ₂}[]
    xs = ℤ₂ⁿ(n)
    I = id(n)
    for x in xs
        # test if vector is in kernel
        if (g + I) * x == x0
            # save vector x
            push!(Vg, x)
        end
    end
    return Vg
end

## 3. quadratic form of ℤ₂
"""
    quadratic_form(x::Vector{ℤ₂})::Int8

Quadratic form ``q: \\mathbb{Z}_2^n \\rightarrow \\mathbb{Z}_4,\\, \\vec{x} \\mapsto \\sum_i ( x_i \\mod{4})``.
"""
function quadratic_form(x::Vector{ℤ₂})::Int8
    # convert upwards ℤ₂ⁿ → ℤ₄ⁿ
    x′ = Vector{ℤ₄}(x)
    # sum all components in mod 4
    s = sum(x′)
    return s
end

## 4. GArf
"""
    GArf(f::Function, V::Vector; readable=true)

Implements the GArf-invariant by applying the function `f` to all vectors in `V` (List of Vectors), see:
``\\text{GArf}(q) = phase\\left\\{  \\sum_{\\vec x\\in V_g} i^{q(\\vec x)}  \\right\\}``
"""
function GArf(f::Function, V::Vector; readable=true)
    # define the imaginary unit
    𝑖 = 1im
    # calculate the phase of the sum of 𝑖^f(x) for all x ∈ V using broadcasting (".")
    φ = angle(sum(𝑖.^f.(V)))
    if readable
        # return as readable values
        return ("$(φ/π)π", "$(rad2deg(φ))°")
    else
        # return plain radian value
        return φ
    end
end

## 5. combine all steps to get the final invariant that is used
"""
    generalized_arf(g::SparseMatrixCSC{ℤ₂, Int64}, n::Int64; readable = true)

Combines the calculation of the kernel `Vg = ker_g_plus_id(g,n)` with the generalized Arf invariant `GArf(quadratic_form, Vg)`.
See [`ker_g_plus_id`](@ref), [`quadratic_form`](@ref) and [`GArf`](@ref).
"""
function generalized_arf(g::SparseMatrixCSC{ℤ₂, Int64}, n::Int64; readable = true)
    # calculate the kernel
    Vg = ker_g_plus_id(g, n)
    # calculate the GArf
    return GArf(quadratic_form, Vg; readable=readable)
end


text:
module VAWTAC

const modulepath = splitdir(@__FILE__)[1]

include("$modulepath/acmultiple.jl")
include("$modulepath/airfoilread.jl")

end

max_stars_repo_name: kevmoor/vawt-ac

max_issues_repo_name: kevmoor/vawt-ac

max_forks_repo_name: kevmoor/vawt-ac


text:
"""
    Represent one text item consisting of original text, translated texts
"""
mutable struct TranslationItem
    id::String
    plural::String
    context::String
    strings::Dict{Int,String}
    TranslationItem() = new("", "", "", Dict{Int,String}())
end

function init_ti!(ti::TranslationItem)
    ti.id = ""; ti.plural = ""; ti.context = ""
    empty!(ti.strings)
end

"""
    read_po_file'('f::AbstractString')'

Read a file, which contains text data according to the PO format of gettext
ref: //https://www.gnu.org/software/gettext/manual/html_node/PO-Files.html
Format is strictly line separated.
A line staring with '#' is a comment line and ignored.
A line starting with a keyword '(msg.*)' has to be followed by a string, enclosed
in '"' in the same line, optionally followed by string continuations in the
following lines. Those strings are concatenated. Characters '"' and '\\' have
to be escaped by a preceding '\\'. The usual C-escape sequences are honored.
Only keywords 'msgid', 'msgid_plural', 'msgstr', and 'msgstr[.*]' are supported.
The "msgid_plural" is used synonymous to the "msgid". If both entries are present,
two separate key-value pairs are generated, with identical '('typical array')' value.
"""
function read_po_file(file::Union{AbstractString,IO})

    in_sequence = false;
    in_keyword = false;
    dict = Vector{Pair{String,Union{String,Vector{String}}}}()
    buffer = IOBuffer()
    keyword = ""
    index = 0
    ti = TranslationItem()
    in_ti = 0

    function begin_sequence(str::AbstractString)
        truncate(buffer, 0)
        write(buffer, str)
        in_sequence = true
    end

    function continue_sequence(str::AbstractString)
        in_sequence || error("no continuation string expected")
        write(buffer, str)
    end

    function end_sequence()
        str = unescape_string(String(take!(buffer)))
        in_sequence = false
        in_keyword && end_keyword(str)
    end

    function begin_keyword(key::AbstractString, ix::Any, str::AbstractString)
        in_keyword = true
        keyword = key
        index = ix != nothing ? Meta.parse(ix) : -1
        begin_sequence(str)
    end

    function end_keyword(line::AbstractString)
        in_keyword = false
        process_keyword(keyword, index, line)
    end

    function process_keyword(key::AbstractString, index::Int, text::AbstractString)
        index == -1 || key == "msgstr" || error("not allowed $key[$index]")
        if key == "msgid"
            in_ti == 3 && process_translation_item(ti)
            in_ti <= 1 || error("unexpected $key '$text'")
            in_ti = 2
            isempty(ti.id) || error("several msgids in line ('$ti.id', '$text')")
            ti.id = text
        elseif key == "msgstr"
            in_ti != 0 || error("unexpected $key '$text'")
            in_ti = 3
            ti.strings[max(index,0)] = text
        elseif key == "msgid_plural"
            in_ti == 2 || error("unexpected $key '$text'")
            isempty(ti.plural) || error("several msgid_plural in line ('$ti.plural', '$text')")
            ti.plural = text
        elseif key == "msgctxt"
            in_ti == 3 && process_translation_item(ti)
            in_ti == 0 || error("unexpected $key '$text'")
            in_ti = 1
            isempty(ti.context) || error("several msgctx in line ('$ti.context', '$text')")
            ti.context = text
        end
    end

    function process_translation_item(ti::TranslationItem)
        add_translation_item!(dict, ti)
        init_ti!(ti)
        in_ti = 0
    end

    for line in eachline(file)
        if ( m = match(REG_KEYWORD, line) ) != nothing
            in_sequence && end_sequence()
            begin_keyword(m.captures[1], m.captures[3], m.captures[4])
        elseif ( m = match(REG_STRING, line) ) != nothing
            continue_sequence(m.captures[1])
        end
    end
    in_sequence && end_sequence()
    in_ti == 3 && process_translation_item(ti)

    dict
end

"""
    read_mo_file(f::AbstractString)

Read a file, which contains text data according to the MO format of gettext.
two separate key-value pairs are generated, with identical '('typical array')' value.
"""
function read_mo_file(f::AbstractString)
    MAGIC = 0x950412de # little-endian form
    NUL = '\x00' # separates different elements of msgid and msgstr (plural forms)
    EOT = '\x04' # separates msgctxt from msgid
    dict = Vector{Pair{String,Union{String,Vector{String}}}}()
    ti = TranslationItem()

    data = open(f, "r") do fp
        read(fp)
    end
    datal = sizeof(data)
    datal >= 28 || error("file too short - no MO file format")
    d = reinterpret(UInt32, data[1:28])
    le_machine = ENDIAN_BOM == 0x04030201
    le_file = d[1] == MAGIC
    le_file || ntoh(d[1]) == MAGIC || error("wrong magic number - no MO-file format")
    conv = le_file ? ltoh : ntoh
    le_machine != le_file && ( d = conv.(d) )
    magic, rev, n, origp, tranp, hsize, hashp = d[1:7]
    revma, revmi = rev >> 16, rev & 0xffff
    origp, tranp, hashp = origp÷4, tranp÷4, hashp÷4
    revma == 0 || error("revision id ($revma,$revmi) in MO-file - only supported: (0, x)")
    datal >= (tranp+2n)*4 || error("file too short - no MO file format")
    d = reinterpret(Int32, data[5:(tranp+2n)*4])
    le_machine != le_file && ( d = conv.(d) )
    for i = 0:2:2n-1
        leno = d[origp+i]
        ptro = d[origp+i+1]
        lent = d[tranp+i]
        ptrt = d[tranp+i+1]
        stro = String(data[ptro+1:ptro+leno])
        strt = String(data[ptrt+1:ptrt+lent])
        ix = firstnn(findfirst(isequal(EOT), stro), 0)
        if ix > 0
            ti.context = stro[1:prevind(stro, ix)]
            stro = stro[nextind(stro, ix):end]
        end
        ix = firstnn(findfirst(isequal(NUL), stro), 0)
        if ix > 0
            ti.plural = stro[nextind(stro, ix):end]
            stro = stro[1:prevind(stro, ix)]
        end
        ti.id = stro
        strtlist = string.(split(strt, NUL))
        ti.strings = Dict(enumerate(strtlist))
        add_translation_item!(dict, ti)
        init_ti!(ti)
    end
    dict
end

firstnn(a::Any) = a
firstnn(a::Nothing, b::Any...) = firstnn(b...)
firstnn(a::Any, b::Any...) = a

# add translation item to output vector
function add_translation_item!(dict::Vector, ti::TranslationItem)
    val = map(p->p.second, sort(collect(ti.strings)))
    if length(val) == 1 && isempty(ti.plural)
        val = val[1]
    end
    !isempty(ti.id) && push!(dict, skey(ti) => val)
    !isempty(ti.plural) && ti.plural != ti.id && push!(dict, pkey(ti) => val)
    isempty(ti.id) && isempty(ti.plural) && push!(dict, "" => val)
    ti
end

# Format strings with and without context information
key(ctx, id) = isempty(id) || isempty(ctx) ? id : string(SCONTEXT, ctx, SCONTEXT, id)
skey(ti::TranslationItem) = key(ti.context, ti.id)
pkey(ti::TranslationItem) = key(ti.context, ti.plural)

"""
    read_header(string)

Extract plural data (nplurals and function plural(n)) from string.
"""
function read_header(str::AbstractString)
    io = IOBuffer(str)
    for line in eachline(io)
        if ( m = match(REG_PLURAL, line) ) != nothing
            return translate_plural_data(m.captures[1])
        end
    end
    translate_plural_data("")
end

module Sandbox end
# clip output of eval(ex(n)) to interval [0, m)
create_plural(ex::Expr, m) = n -> max(min(Base.invokelatest(Sandbox.eval(ex), n), m-1), 0)

# avoid the following error when calling f by invokelatest:
# "MethodError: no method matching (::getfield(ResourceBundles, Symbol("...")))(::Int64)
# The applicable method may be too new: running in world age ~3, while current world is ~4."
# This is maybe an undesirable hack - looking for more elegant work around.

"""
    translate_plural_data(string)

Parses a string of form `"nplurals = expr1; plural = expr(n)"`.
`expr1` must be a constant integer expression.
`expr(n) must be an integer arthmetic expression with a variable `n`.

Outputs the evaluation of `Int(expr1)` and of `n::Int -> expr(n)`.
"""
function translate_plural_data(str::AbstractString)
    nplurals = 2
    plural = n -> n != 0
    str = replace(str, ':' => " : ") # surround : by blanks
    str = replace(str, '?' => " ? ") # surround ? by blanks
    str = replace(str, "/" => "÷") # use Julia integer division for /
    top = Meta.parse(str)
    isa(top, Expr) || return nplurals, plural
    for a in top.args
        if isa(a, Expr) && a.head == :(=)
            var, ex = a.args
            if var == :nplurals
                nplurals = Int(Sandbox.eval(ex)) # evaluate right hand side
            elseif var == :plural
                ex2 = Expr(:(->), :n, ex)
                plural = create_plural(ex2, nplurals)
                @assert plural(1) == 0 string(ex, " not 0 for n == 1")
            end
        end
    end
    nplurals, plural
end

const REG_KEYWORD = r"""^\s*([a-zA-Z_]+)(\[(\d+)\])?\s+"(.*)"\s*$"""
const REG_STRING = r"""^\s*"(.*)"\s*$"""
const REG_COMMENT = r"""^\s*#"""

const REG_PLURAL = r"""^\s*Plural-Forms:(.*)$"""


text:
#!/usr/bin/env julia
# del2z <delta.z@aliyun.com>

include("wiki.jl")

run(`pwd`)
parsewiki("/Users/del2z/Workspace/Corpora/wiki/demo",
          "temp.txt")
#parsewiki("../../../../Resources/Corpus/zhwiki-latest-pages-articles6.xml-p6231444p6382070",
#          "temp.txt")


text:
# Counts of discrete values

#################################################
#
#  counts on given levels
#
#################################################

const IntUnitRange{T<:Integer} = UnitRange{T}

if isdefined(Base, :ht_keyindex2)
    const ht_keyindex2! = Base.ht_keyindex2
else
    using Base: ht_keyindex2!
end

#### functions for counting a single list of integers (1D)
"""
    addcounts!(r, x, levels::UnitRange{<:Int}, [wv::AbstractWeights])

Add the number of occurrences in `x` of each value in `levels` to an existing
array `r`. If a weighting vector `wv` is specified, the sum of weights is used
rather than the raw counts.
"""
function addcounts!(r::AbstractArray, x::IntegerArray, levels::IntUnitRange)
    # add counts of integers from x to r

    k = length(levels)
    length(r) == k || throw(DimensionMismatch())

    m0 = levels[1]
    m1 = levels[end]
    b = m0 - 1

    @inbounds for i in 1 : length(x)
        xi = x[i]
        if m0 <= xi <= m1
            r[xi - b] += 1
        end
    end
    return r
end

function addcounts!(r::AbstractArray, x::IntegerArray, levels::IntUnitRange, wv::AbstractWeights)
    k = length(levels)
    length(r) == k || throw(DimensionMismatch())

    m0 = levels[1]
    m1 = levels[end]
    b = m0 - 1

    @inbounds for i in 1 : length(x)
        xi = x[i]
        if m0 <= xi <= m1
            r[xi - b] += wv[i]
        end
    end
    return r
end

"""
    counts(x, [wv::AbstractWeights])
    counts(x, levels::UnitRange{<:Integer}, [wv::AbstractWeights])
    counts(x, k::Integer, [wv::AbstractWeights])

Count the number of times each value in `x` occurs. If `levels` is provided, only values
falling in that range will be considered (the others will be ignored without
raising an error or a warning). If an integer `k` is provided, only values in the
range `1:k` will be considered.

If a weighting vector `wv` is specified, the sum of the weights is used rather than the
raw counts.

The output is a vector of length `length(levels)`.
"""
function counts end

counts(x::IntegerArray, levels::IntUnitRange) =
    addcounts!(zeros(Int, length(levels)), x, levels)
counts(x::IntegerArray, levels::IntUnitRange, wv::AbstractWeights) =
    addcounts!(zeros(eltype(wv), length(levels)), x, levels, wv)
counts(x::IntegerArray, k::Integer) = counts(x, 1:k)
counts(x::IntegerArray, k::Integer, wv::AbstractWeights) = counts(x, 1:k, wv)
counts(x::IntegerArray) = counts(x, span(x))
counts(x::IntegerArray, wv::AbstractWeights) = counts(x, span(x), wv)

"""
    proportions(x, levels=span(x), [wv::AbstractWeights])

Return the proportion of values in the range `levels` that occur in `x`.
Equivalent to `counts(x, levels) / length(x)`. If a weighting vector `wv`
is specified, the sum of the weights is used rather than the raw counts.
"""
proportions(x::IntegerArray, levels::IntUnitRange) = counts(x, levels) .* inv(length(x))
proportions(x::IntegerArray, levels::IntUnitRange, wv::AbstractWeights) =
    counts(x, levels, wv) .* inv(sum(wv))

"""
    proportions(x, k::Integer, [wv::AbstractWeights])

Return the proportion of integers in 1 to `k` that occur in `x`.
"""
proportions(x::IntegerArray, k::Integer) = proportions(x, 1:k)
proportions(x::IntegerArray, k::Integer, wv::AbstractWeights) = proportions(x, 1:k, wv)
proportions(x::IntegerArray) = proportions(x, span(x))
proportions(x::IntegerArray, wv::AbstractWeights) = proportions(x, span(x), wv)

#### functions for counting a single list of integers (2D)

function addcounts!(r::AbstractArray, x::IntegerArray, y::IntegerArray, levels::NTuple{2,IntUnitRange})
    # add counts of integers from x to r

    n = length(x)
    length(y) == n || throw(DimensionMismatch())

    xlevels, ylevels = levels

    kx = length(xlevels)
    ky = length(ylevels)
    size(r) == (kx, ky) || throw(DimensionMismatch())

    mx0 = xlevels[1]
    mx1 = xlevels[end]
    my0 = ylevels[1]
    my1 = ylevels[end]

    bx = mx0 - 1
    by = my0 - 1

    for i = 1:n
        xi = x[i]
        yi = y[i]
        if (mx0 <= xi <= mx1) && (my0 <= yi <= my1)
            r[xi - bx, yi - by] += 1
        end
    end
    return r
end

function addcounts!(r::AbstractArray, x::IntegerArray, y::IntegerArray,
                    levels::NTuple{2,IntUnitRange}, wv::AbstractWeights)
    # add counts of integers from x to r

    n = length(x)
    length(y) == length(wv) == n || throw(DimensionMismatch())

    xlevels, ylevels = levels

    kx = length(xlevels)
    ky = length(ylevels)
    size(r) == (kx, ky) || throw(DimensionMismatch())

    mx0 = xlevels[1]
    mx1 = xlevels[end]
    my0 = ylevels[1]
    my1 = ylevels[end]

    bx = mx0 - 1
    by = my0 - 1

    for i = 1:n
        xi = x[i]
        yi = y[i]
        if (mx0 <= xi <= mx1) && (my0 <= yi <= my1)
            r[xi - bx, yi - by] += wv[i]
        end
    end
    return r
end

# facet functions

function counts(x::IntegerArray, y::IntegerArray, levels::NTuple{2,IntUnitRange})
    addcounts!(zeros(Int, length(levels[1]), length(levels[2])), x, y, levels)
end

function counts(x::IntegerArray, y::IntegerArray, levels::NTuple{2,IntUnitRange}, wv::AbstractWeights)
    addcounts!(zeros(eltype(wv), length(levels[1]), length(levels[2])), x, y, levels, wv)
end

counts(x::IntegerArray, y::IntegerArray, levels::IntUnitRange) =
    counts(x, y, (levels, levels))
counts(x::IntegerArray, y::IntegerArray, levels::IntUnitRange, wv::AbstractWeights) =
    counts(x, y, (levels, levels), wv)

counts(x::IntegerArray, y::IntegerArray, ks::NTuple{2,Integer}) =
    counts(x, y, (1:ks[1], 1:ks[2]))
counts(x::IntegerArray, y::IntegerArray, ks::NTuple{2,Integer}, wv::AbstractWeights) =
    counts(x, y, (1:ks[1], 1:ks[2]), wv)
counts(x::IntegerArray, y::IntegerArray, k::Integer) = counts(x, y, (1:k, 1:k))
counts(x::IntegerArray, y::IntegerArray, k::Integer, wv::AbstractWeights) =
    counts(x, y, (1:k, 1:k), wv)
counts(x::IntegerArray, y::IntegerArray) = counts(x, y, (span(x), span(y)))
counts(x::IntegerArray, y::IntegerArray, wv::AbstractWeights) = counts(x, y, (span(x), span(y)), wv)

proportions(x::IntegerArray, y::IntegerArray, levels::NTuple{2,IntUnitRange}) =
    counts(x, y, levels) .* inv(length(x))
proportions(x::IntegerArray, y::IntegerArray, levels::NTuple{2,IntUnitRange}, wv::AbstractWeights) =
    counts(x, y, levels, wv) .* inv(sum(wv))

proportions(x::IntegerArray, y::IntegerArray, ks::NTuple{2,Integer}) =
    proportions(x, y, (1:ks[1], 1:ks[2]))
proportions(x::IntegerArray, y::IntegerArray, ks::NTuple{2,Integer}, wv::AbstractWeights) =
    proportions(x, y, (1:ks[1], 1:ks[2]), wv)
proportions(x::IntegerArray, y::IntegerArray, k::Integer) = proportions(x, y, (1:k, 1:k))
proportions(x::IntegerArray, y::IntegerArray, k::Integer, wv::AbstractWeights) =
    proportions(x, y, (1:k, 1:k), wv)
proportions(x::IntegerArray, y::IntegerArray) = proportions(x, y, (span(x), span(y)))
proportions(x::IntegerArray, y::IntegerArray, wv::AbstractWeights) =
    proportions(x, y, (span(x), span(y)), wv)

#################################################
#
#  countmap on unknown levels
#
#  These methods are based on dictionaries, and
#  can be used on any kind of hashable values.
#
#################################################

## auxiliary functions

function _normalize_countmap(cm::Dict{T}, s::Real) where T
    r = Dict{T,Float64}()
    for (k, c) in cm
        r[k] = c / s
    end
    return r
end

## 1D

"""
    addcounts!(dict, x[, wv]; alg = :auto)

Add counts based on `x` to a count map. New entries will be added if new values come up.
If a weighting vector `wv` is specified, the sum of the weights is used rather than the
raw counts.

`alg` can be one of:
- `:auto` (default): if `StatsBase.radixsort_safe(eltype(x)) == true` then use
                     `:radixsort`, otherwise use `:dict`.

- `:radixsort`:      if `radixsort_safe(eltype(x)) == true` then use the
                     [radix sort](https://en.wikipedia.org/wiki/Radix_sort)
                     algorithm to sort the input vector which will generally lead to
                     shorter running time. However the radix sort algorithm creates a
                     copy of the input vector and hence uses more RAM. Choose `:dict`
                     if the amount of available RAM is a limitation.

- `:dict`:           use `Dict`-based method which is generally slower but uses less
                     RAM and is safe for any data type.
"""
addcounts!(cm::Dict, x; alg = :auto) = _addcounts!(eltype(x), cm, x, alg = alg)

function _addcounts!(::Type{T}, cm::Dict, x; alg = :auto) where T
    # if it's safe to be sorted using radixsort then it should be faster
    # albeit using more RAM
    if radixsort_safe(T) && (alg == :auto || alg == :radixsort)
        addcounts_radixsort!(cm, x)
    elseif alg == :radixsort
        throw(ArgumentError("`alg = :radixsort` is chosen but type `radixsort_safe($T)` did not return `true`; use `alg = :auto` or `alg = :dict` instead"))
    else
        addcounts_dict!(cm,x)
    end
    return cm
end

"""Dict-based addcounts method"""
function addcounts_dict!(cm::Dict{T}, x) where T
    for v in x
        index = ht_keyindex2!(cm, v)
        if index > 0
            @inbounds cm.vals[index] += 1
        else
            @inbounds Base._setindex!(cm, 1, v, -index)
        end
    end
    return cm
end

# If the bits type is of small size i.e. it can have up to 65536 distinct values
# then it is always better to apply a counting-sort like reduce algorithm for
# faster results and less memory usage. However we still wish to enable others
# to write generic algorithms, therefore the methods below still accept the
# `alg` argument but it is ignored.
function _addcounts!(::Type{Bool}, cm::Dict{Bool}, x::AbstractArray{Bool}; alg = :ignored)
    sumx = sum(x)
    cm[true] = get(cm, true, 0) + sumx
    cm[false] = get(cm, false, 0) + length(x) - sumx
    cm
end

# specialized for `Bool` iterator
function _addcounts!(::Type{Bool}, cm::Dict{Bool}, x; alg = :ignored)
    sumx = 0
    len = 0
    for i in x
        sumx += i
        len += 1
    end
    cm[true] = get(cm, true, 0) + sumx
    cm[false] = get(cm, false, 0) + len - sumx
    cm
end

function _addcounts!(::Type{T}, cm::Dict{T}, x; alg = :ignored) where T <: Union{UInt8, UInt16, Int8, Int16}
    counts = zeros(Int, 2^(8sizeof(T)))

    @inbounds for xi in x
        counts[Int(xi) - typemin(T) + 1] += 1
    end

    for (i, c) in zip(typemin(T):typemax(T), counts)
        if c != 0
            index = ht_keyindex2!(cm, i)
            if index > 0
                @inbounds cm.vals[index] += c
            else
                @inbounds Base._setindex!(cm, c, i, -index)
            end
        end
    end
    cm
end

const BaseRadixSortSafeTypes = Union{Int8, Int16, Int32, Int64, Int128,
                                     UInt8, UInt16, UInt32, UInt64, UInt128,
                                     Float32, Float64}

"Can the type be safely sorted by radixsort"
radixsort_safe(::Type{T}) where T = T<:BaseRadixSortSafeTypes

function _addcounts_radix_sort_loop!(cm::Dict{T}, sx::AbstractArray{T}) where T
    isempty(sx) && return cm
    last_sx = sx[1]
    tmpcount = get(cm, last_sx, 0) + 1

    # now the data is sorted: can just run through and accumulate values before
    # adding into the Dict
    @inbounds for i in 2:length(sx)
        sxi = sx[i]
        if last_sx == sxi
            tmpcount += 1
        else
            cm[last_sx] = tmpcount
            last_sx = sxi
            tmpcount = get(cm, last_sx, 0) + 1
        end
    end

    cm[sx[end]] = tmpcount

    return cm
end

function addcounts_radixsort!(cm::Dict{T}, x::AbstractArray{T}) where T
    # sort the x using radixsort
    sx = sort(x, alg = RadixSort)

    # Delegate the loop to a separate function since sort might not
    # be inferred in Julia 0.6 after SortingAlgorithms is loaded.
    # It seems that sort is inferred in Julia 0.7.
    return _addcounts_radix_sort_loop!(cm, sx)
end

# fall-back for `x` an iterator
function addcounts_radixsort!(cm::Dict{T}, x) where T
    sx = sort!(collect(x), alg = RadixSort)
    return _addcounts_radix_sort_loop!(cm, sx)
end

function addcounts!(cm::Dict{T}, x::AbstractArray{T}, wv::AbstractVector{W}) where {T,W<:Real}
    n = length(x)
    length(wv) == n || throw(DimensionMismatch())
    z = zero(W)

    for i = 1 : n
        @inbounds xi = x[i]
        @inbounds wi = wv[i]
        cm[xi] = get(cm, xi, z) + wi
    end
    return cm
end

"""
    countmap(x; alg = :auto)
    countmap(x::AbstractVector, w::AbstractVector{<:Real}; alg = :auto)

Return a dictionary mapping each unique value in `x` to its number
of occurrences. A vector of weights `w` can be provided when `x` is a vector.

- `:auto` (default): if `StatsBase.radixsort_safe(eltype(x)) == true` then use
                     `:radixsort`, otherwise use `:dict`.

- `:radixsort`:      if `radixsort_safe(eltype(x)) == true` then use the
                     [radix sort](https://en.wikipedia.org/wiki/Radix_sort)
                     algorithm to sort the input vector which will generally lead to
                     shorter running time. However the radix sort algorithm creates a
                     copy of the input vector and hence uses more RAM. Choose `:dict`
                     if the amount of available RAM is a limitation.

- `:dict`:           use `Dict`-based method which is generally slower but uses less
                     RAM and is safe for any data type.
"""
countmap(x; alg = :auto) = addcounts!(Dict{eltype(x),Int}(), x; alg = alg)
countmap(x::AbstractArray{T}, wv::AbstractVector{W}) where {T,W<:Real} = addcounts!(Dict{T,W}(), x, wv)

"""
    proportionmap(x)

Return a dictionary mapping each unique value in `x` to its
proportion in `x`.
"""
proportionmap(x::AbstractArray) = _normalize_countmap(countmap(x), length(x))
proportionmap(x::AbstractArray, wv::AbstractWeights) = _normalize_countmap(countmap(x, wv), sum(wv))


text:
@testset "surface_mesh" begin
    dist = SurfaceDistribution(1e-3, 1e-4, .025, .05)
    af = NACA"0012"s
    x0, y0 = gen_airfoil(af)
    points = surface_mesh(dist, x0, y0)
    xs = getindex.(points, 1)
    ys = getindex.(points, 2)
    im = argmin(xs)
    x1 = xs[2:im-1]
    y1 = ys[2:im-1]
    dy = abs.(y_lower(af, x1)-y1)
    @test all(dy .< 1e-4)
    @test ys[2] < ys[end-1]
    @test points[1] === points[end]
    af = NACA"4412"
    x0, y0 = gen_airfoil(af)
    points = surface_mesh(dist, x0, y0)
    @test points[2][2] < points[end-1][2]
    @test points[1] === points[end]
end


text:
#######################################################################
# Coverage.jl
# Take Julia test coverage results and bundle them up in JSONs
# https://github.com/IainNZ/Coverage.jl
#######################################################################
module Coverage

    # process_cov
    # Given a .cov file, return the counts for each line, where the
    # lines that can't be counted are denoted with a -1
    export process_cov
    function process_cov(filename)
        fp = open(filename, "r")
        lines = readlines(fp)
        num_lines = length(lines)
        coverage = Array(Union(Nothing,Int), num_lines)
        for i = 1:num_lines
            cov_segment = lines[i][1:9]
            coverage[i] = cov_segment[9] == '-' ? nothing : int(cov_segment)
        end
        close(fp)
        return coverage
    end

    export Coveralls
    module Coveralls
        using Requests
        using Coverage
        using JSON

        # coveralls_process_file
        # Given a .jl file, return the Coveralls.io dictionary for this
        # file by reading in the file and its matching .cov. Don't convert
        # to JSON yet, just return dictionary.
        # https://coveralls.io/docs/api
        # {
        #   "name" : "$filename"
        #   "source": "...\n....\n...."
        #   "coverage": [null, 1, null]
        # }
        export process_file
        function process_file(filename)
            return ["name" => filename,
                    "source" => readall(filename),
                    "coverage" => process_cov(filename*".cov")]
        end

        # coveralls_process_src
        # Recursively walk through a Julia package's src/ folder
        # and collect coverage statistics
        export process_folder
        function process_folder(folder="src",source_files={})
            filelist = readdir(folder)
            for file in filelist
                fullfile = joinpath(folder,file)
                println(fullfile)
                if isfile(fullfile)
                    try
                        new_sf = process_file(fullfile)
                        push!(source_files, new_sf)
                    catch e
                        if !isa(e,SystemError)
                            rethrow(e)
                        end
                        # Skip
                        println("Skipped $fullfile")
                    end
                else isdir(fullfile)
                    process_folder(fullfile,source_files)
                end
            end
            if folder == "src"
                return source_files
            end
        end

        # submit
        # Submit coverage to Coveralls.io
        # https://coveralls.io/docs/api
        # {
        #   "service_job_id": "1234567890",
        #   "service_name": "travis-ci",
        #   "source_files": [
        #     {
        #       "name": "example.rb",
        #       "source": "def four\n  4\nend",
        #       "coverage": [null, 1, null]
        #     },
        #     {
        #       "name": "lib/two.rb",
        #       "source": "def seven\n  eight\n  nine\nend",
        #       "coverage": [null, 1, 0, null]
        #     }
        #   ]
        # }
        export submit, submit_token
        function submit(source_files)
            data = ["service_job_id" => ENV["TRAVIS_JOB_ID"],
                    "service_name" => "travis-ci",
                    "source_files" => source_files]
            r = Requests.post(URI("https://coveralls.io/api/v1/jobs"), files =
                [FileParam(JSON.json(data),"application/json","json_file","coverage.json")])
            dump(r.data)
        end

        function submit_token(source_files)
            data = ["repo_token" => ENV["REPO_TOKEN"],
                    "source_files" => source_files]
            r = post(URI("https://coveralls.io/api/v1/jobs"), files =
                [FileParam(JSON.json(data),"application/json","json_file","coverage.json")])
            dump(r.data)
        end
    end  # module Coveralls
end


text:
function coupledistal!(system::CVSystem,n::Int64,terms::Vector{Int64},hemoflag="no")
    # W1 at next time step
    for i = 1:length(terms)
        ret1 = Float64[0];
        CVModule.endinvariants!(ret1,system.branches.Q[terms[i]][system.solverparams.JL,n+1],
            system.branches.A[terms[i]][system.solverparams.JL,n+1],system.branches.Q[terms[i]][system.solverparams.JL-1,n+1],
            system.branches.A[terms[i]][system.solverparams.JL-1,n+1],system.branches.c0[terms[i]][end],
            system.branches.beta[terms[i]][end],system.solverparams.rho,system.solverparams.mu,
            system.solverparams.diffusioncoeff,system.branches.k[terms[i]],system.solverparams.h);
        system.branches.W1end[terms[i]] = ret1[1];
    end
    if hemoflag == "no"
        for i = 1:length(terms)
            yout = zeros(4);
            Vs = system.branches.term[terms[i]].V[n+1,1];
            vs = system.branches.c0[terms[i]][end];
            ts = system.branches.k[terms[i]]/vs;
            iters = Int64[0];
            CVModule.newtondist!(yout,iters,Vs,vs,ts,system.branches.term[terms[i]].V[n+1,1],
                system.branches.term[terms[i]].C[1],system.branches.term[terms[i]].Q[n+1,1],
                system.branches.term[terms[i]].V0[1],system.branches.beta[terms[i]][end],
                system.branches.c0[terms[i]][end],system.branches.A0[terms[i]][end],
                system.branches.W1end[terms[i]],system.solverparams.rho,CVModule.fdist,
                CVModule.Jdist,system.solverparams.maxiter,system.solverparams.epsJ,
                system.solverparams.epsN,system.solverparams.maxval,system.solverparams.h);
            system.solverparams.totaliter += iters[1];
            system.branches.term[terms[i]].V[n+2,1] = yout[1];
            system.branches.term[terms[i]].P[n+2,1] = yout[2];
            system.branches.A[terms[i]][system.solverparams.JL,n+2] = yout[3];
            system.branches.Q[terms[i]][system.solverparams.JL,n+2] = yout[4];
        end
    elseif hemoflag == "yes"
        for i = 1:length(terms)
            if system.hemo.injured[terms[i]] == false
                yout = zeros(4);
                Vs = system.branches.term[terms[i]].V[n+1,1];
                vs = system.branches.c0[terms[i]][end];
                ts = system.branches.k[terms[i]]/vs;
                iters = Int64[0];
                CVModule.newtondist!(yout,iters,Vs,vs,ts,system.branches.term[terms[i]].V[n+1,1],
                    system.branches.term[terms[i]].C[1],system.branches.term[terms[i]].Q[n+1,1],
                    system.branches.term[terms[i]].V0[1],system.branches.beta[terms[i]][end],
                    system.branches.c0[terms[i]][end],system.branches.A0[terms[i]][end],
                    system.branches.W1end[terms[i]],system.solverparams.rho,CVModule.fdist,
                    CVModule.Jdist,system.solverparams.maxiter,system.solverparams.epsJ,
                    system.solverparams.epsN,system.solverparams.maxval,system.solverparams.h);
                system.solverparams.totaliter += iters[1];
                system.branches.term[terms[i]].V[n+2,1] = yout[1];
                system.branches.term[terms[i]].P[n+2,1] = yout[2];
                system.branches.A[terms[i]][system.solverparams.JL,n+2] = yout[3];
                system.branches.Q[terms[i]][system.solverparams.JL,n+2] = yout[4];
            elseif system.hemo.injured[terms[i]] == true
                children = system.branches.children[terms[i]];
                if ~isempty(children) # bifurncation downstream
                    beta = zeros(length(children[terms[i]])+1);
                    A0 = zeros(length(children[terms[i]])+1);
                    c0 = zeros(length(children[terms[i]])+1);
                    W2 = zeros(length(children[terms[i]]));
                    yout = zeros(4+3*length(children));
                    beta[1] = system.branches.beta[terms[i]][end];
                    A0[1] = system.branches.A0[terms[i]][end];
                    c0[1] = system.branches.c0[terms[i]][end];
                    for i = 1:length(children)
                        beta[i+1] = system.branches.beta[children[i]][end];
                        A0[i+1] = system.branches.A0[children[i]][end];
                        c0[i+1] = system.branches.c0[children[i]][end];
                        W2[i] = system.branches.W2[children[i]];
                    end
                    CVModule.model1dhemo!(yout,children,system.hemo.Ph,beta,
                        system.solverparams.rho,c0,A0,system.solverparams.h,
                        system.branches.W1[terms[i]],W2);
                    system.branches.A[terms[i]][system.solverparams.JL,n+2] = yout[1];
                    system.branches.W2[terms[i]] = yout[2];
                    system.branches.Q[terms[i]][system.solverparams.JL,n+2] = yout[3];
                    system.hemo.Vloss += yout[4];
                    for i = 1:length(children)
                        system.branches.A[children[i]][1,n+2] = yout[5+3*(i-1)];
                        system.branches.W1[children[i]] = yout[6+3*(i-1)];
                        system.branches.Q[children[i]][1,n+2] = yout[7+3*(i-1)];
                    end
                else # terminal downstream
                    yout = zeros(5);
                    beta = [system.branches.beta[terms[i]][end]];
                    A0 = [system.branches.A0[terms[i]][end]];
                    c0 = [system.branches.c0[terms[i]][end]];
                    CVModule.model1dhemo!(yout,children,system.hemo.Ph,beta,
                        system.solverparams.rho,c0,A0,system.solverparams.h,
                        system.branches.W1end[terms[i]],0.,system.branches.term[terms[i]].V[n+1,1],
                        system.branches.term[terms[i]].Q[n+1,1]);
                    system.branches.A[terms[i]][system.solverparams.JL,n+2] = yout[1];
                    system.branches.Q[terms[i]][system.solverparams.JL,n+2] = yout[2];
                    system.hemo.Vloss += yout[3];
                    system.branches.term[terms[i]].V[n+2,1] = yout[4];
                    system.branches.term[terms[i]].P[n+2,1] = yout[5];
                end
            end
        end
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
using Flux: σ
using Zygote
using DNC: contentaddress

function generate_Sn(dims...)
    arr = rand(Float32, dims...)
    softmax(arr; dims=1)
end

function generate_Δn(dims...)
    fuzzysoftmax(arr; dims=1) =
        exp.(arr)./(sum(exp.(arr); dims=dims) .+ rand(Float32, size(arr))/2)
    arr = rand(Float32, dims...)
    fuzzysoftmax(arr; dims=1)
end

in_Sn(vec) = sum(vec) == 1.0
in_Δn(vec) = sum(vec) >= 0.0 && sum(vec) <= 1.0

N, W, R, B = 4, 8, 2, 4
keyr = rand(Float32, W, R, B)
keyw = rand(Float32, W, 1, B)
mem = rand(Float32, N, W, B)
βr = DNC.oneplus.(rand(Float32, R, B))
βw = DNC.oneplus.(rand(Float32, 1, B))
inputs = (wr=rand(Float32, N, R, B),
          f=σ.(rand(Float32, R, B)),
          ww=rand(Float32, N, 1, B),
          readmode=rand(Float32, 3, R, B))
state = DNC.State(
    zeros(Float32, N, N, B),
    zeros(Float32, N, B),
    zeros(Float32, N, B),
    zeros(Float32, N, 1, B),
    zeros(Float32, N, R, B))

@testset "Batch-training" begin
    cr= DNC.contentaddress(keyr, mem, βr)
    @test size(cr) == (N, R, B)
    @test eltype(cr) == Float32
    cw = DNC.contentaddress(keyw, mem, βw)
    @test size(cw) == (N, 1, B)
    memret = DNC.memoryretention(inputs.wr, inputs.f)
    @test eltype(memret) == Float32
    @test size(memret) == (N, B)
    u = DNC.usage(state.u, state.ww, state.wr, inputs.f)
    @test eltype(u) == Float32
    a = DNC.allocationweighting(u)
    @test eltype(a) == Float32
    gw, ga = rand(Float32, 1, B), rand(Float32, 1, B)
    ww = DNC.writeweight(cw, a, gw, ga)
    @test eltype(ww) == Float32
    @test size(ww) == (N, 1, B)
    forw = DNC.forwardweight(state.L, state.wr)
    backw = DNC.backwardweight(state.L, state.wr)
    wr = DNC.readweight(backw, cr, forw, inputs.readmode)
    @test eltype(forw) == Float32
    @test eltype(backw) == Float32
    @test eltype(wr) == Float32
    @test size(wr) == (N, R, B)
    p = DNC.precedenceweight(state.p, ww)
    @test size(p) == (N, B)
    L = zeros(Float32, N, N, B) # Don't interfere with later tests
    DNC.updatelinkmatrix!(L, p, ww)
    @test eltype(L) == Float32

end

@testset "Content-based addressing" begin
    M =reshape(Float32[0.1 0.5 1.5;
                       -1.2 0.8 0.0], 2, 3, 1)
    key = rand(Float32, 3, 1, 1)
    key[:, 1, 1] = M[1, : , 1]
    key[1, 1, 1] = 0 # Avoid exact match
    β = reshape([100.0f0, 100.0f0], 1, 2)
    # With a high sharpener β, the match is found
    res = contentaddress(key, M, β)
    @test res[1] ≈ 1
    @test eltype(res) == Float32
    # Test division by zero
    M[2, :] = zeros(Float32, 3)
    contentzero = contentaddress(key, M, β)
    @test !any(isnan.(contentzero))
    gs = gradient((k, m, b)->sum(contentaddress(k, m, b)), key, M, β)
    @test !any(isnan.(gs[1]))
    @test !any(isnan.(gs[2]))
    @test !any(isnan.(gs[3]))
end

@testset "Memory allocation" begin
    N, W, R, B = 3, 4, 1, 1
    wr1 = rand(Float32, N, R, B)
    wr1[:, 1, 1] = [0.5f0, 0.25f0, 0.25f0]
    f1 = ones(Float32, R, B)
    ww1 = rand(Float32, N, 1, B)
    ww1[:, 1, 1] = [0.25f0, 0.5f0, 0.25f0]
    case1 = (wr=wr1, f=f1, ww=ww1)
    R = 2
    wr2 = rand(Float32, N, R, B)
    wr2[:, :, 1] = Float32.([0.6 0.0;
                             0.3 0.5;
                             0.1 0.5])
    f2 = rand(Float32, R, B)
    f2[:, 1] = [1.f0, 0.f0]
    ww2 = rand(Float32, N, 1, B)
    ww2[:, 1, 1] = [0.25f0, 0.5f0, 0.25f0]
    case2 = (wr=wr2, f=f2, ww=ww2)
    u_prev = rand(Float32, N, B)
    u_prev[:, 1] = [1.0f0, 0.0f0, 0.0f0]
    state2 = DNC.State(N, R, B)
    state2.ww = case2.ww
    state2.wr = case2.wr
    state2.u = u_prev

    @testset "Memory retention 𝜓" begin
        memret = DNC.memoryretention(case1.wr, case1.f)
        @test memret[:, 1] == [0.5, 0.75, 0.75]
        @test eltype(memret) == Float32
        # Two read heads
        memret = DNC.memoryretention(case2.wr, case2.f)
        @test isapprox(memret[:, 1], [0.4, 0.7, 0.9], atol=1e-5)
        @test eltype(memret) == Float32
    end

    @testset "Usage u⃗" begin
        wr, f, ww = case1
        u = DNC.usage(u_prev, ww, wr, f)
        @test u[:, 1] == [1//2, 3//8, 3//16]
        @test eltype(u) == Float32
        # Two read heads
        wr, f, ww = case2
        u = DNC.usage(u_prev, ww, wr, f)
        @test isapprox(u[:, 1], [0.4, 0.35, 0.225], atol=1e-5)
        @test eltype(u) == Float32
    end

    @testset "Allocation a⃗" begin
        # Using approximation due to DNC's use of _EPSILON to avoid num. instability
        wr, f, ww = case2
        u = DNC.usage(u_prev,ww, wr, f)
        alloc = DNC.allocationweighting(u)
        @test eltype(alloc) == Float32
        @test isapprox(alloc[:, 1], [0.04725, 0.14625, 0.775]; atol=DNC._EPSILON*10)
        # Allocation is zero if all usages are 1
        allused = ones(5, B)
        @test isapprox(DNC.allocationweighting(allused), (zeros(5, B)); atol=DNC._EPSILON*10)
        # First location is allocated if all is available
        noneused = zeros(5, B)
        @test isapprox(DNC.allocationweighting(noneused)[:, 1], [1.0, 0, 0, 0, 0]; atol=DNC._EPSILON*10)
    end

    @testset "Write weighting ww" begin
        ga = zeros(Float32, 1, B)
        gw = ones(Float32, 1, B)
        a = generate_Δn(N, B)
        cw = generate_Sn(N, 1, B)
        ww = DNC.writeweight(cw, a, gw, ga)
        @test ww == cw
        @test eltype(ww) == Float32
    end
end

@testset "Temporal linkage" begin
    @testset "Precedence weights p⃗" begin
        N, B = 4, 1
        ww = rand(Float32, N, 1, B)
        ww[:, :, 1] = [0.1f0, 0.2f0, 0.3f0, 0.4f0]
        p0 = zeros(Float32, N, B)
        p = DNC.precedenceweight(p0, ww)
        @test p == dropdims(ww; dims=2)
        @test eltype(p) == Float32
    end

    # Write weights. Writes to location 1, 2, and 3 iteratively.
    N, R, B = 3, 1, 1
    ww = rand(Float32, N, 1, B)
    wwlist = [[1.0f0, 0.0f0, 0.0f0],
            [0.0f0, 1.0f0, 0.0f0],
            [0.0f0, 0.0f0, 1.0f0]]
    # Temporal link matrix
    L = zeros(Float32, N, N, B)
    # precedence weights
    p = zeros(Float32, N, B)
    # Expected evolution of L
    expected = [zeros(Float32, N, N, B) for i in 1:3]
    expected[2][2, 1, 1] = 1.0f0
    expected[3][2, 1, 1] = 1.0f0
    expected[3][3, 2, 1] = 1.0f0

    @testset "Link matrix L" begin
        for i in 1:length(wwlist)
            ww[:, :, 1] = wwlist[i]
            DNC.updatelinkmatrix!(L, p, ww)
            @test L == expected[i]
            @test eltype(L) == Float32
            p = DNC.precedenceweight(p, ww)
        end
        # Last read was location 2. The forward weight should point to 3, backward to 1.
        wr = rand(Float32, N, R, B)
        wr[:, :, 1] = [0.f0, 1, 0]
        forward = DNC.forwardweight(L, wr)
        backward = DNC.backwardweight(L, wr)
        @test  forward[:, 1] == [0, 0, 1]
        @test backward[:, 1] == [1, 0, 0]
        @test eltype(forward) == Float32
        @test eltype(backward) == Float32
    end

    @testset "Read weighting w_t^r" begin
        b = reshape([0.f0, 0, 1], 3, 1, 1)
        f = reshape([1.f0, 0, 0], 3, 1, 1)
        c = reshape([0.f0, 1, 0], 3, 1, 1)
        readmode = reshape([0.f0, 1, 0], 3, 1, 1)
        readw = DNC.readweight(b, c, f, readmode)
        @test readw == c
        @test in_Δn(readw)
        @test eltype(readw) == Float32
        readmode2 = reshape([0.5f0, 0.5f0, 0.0f0], 3, 1, 1)
        readw = DNC.readweight(b, c, f, readmode2)
        @test readw == 0.5f0b .+ 0.5f0c
    end

end


text:
########################
#
# These are models that I tried during the development of the kalman filter
# In the end I settled on the basic model included in kalman.jl
#
########################

# function EASSR_model(f::Number, fs::Number; x0::Vector=[0.0, 0, 0, 150, 0, 0, 0, 0, 0, 0, 1], W::Number=1.0)
#
#     Δt = 1 / fs
#     ω = 2 * pi * f
#
#     ## Process transition and noise covariance
#
#     a = zeros(11,11)
#     a[1:2, 1:2] = eye(2)
#     a[4:11, 3:10] = eye(8)
#     a[3, 11] = 1
#
#     F = a
#
#     V = diagm([1e-10, 1e-10, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001])  # We expect little variation in the state matrix
#
#     ## Observation and noise covariance
#
#     function G1(k);  cos(ω * Δt * k); end
#     function G2(k); -sin(ω * Δt * k); end
#     G = reshape([G1, G2, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 1, 11)
#     W = diagm([W])   # Noise Covariance
#
#     ## Inital guesses at state and error covariance
#
#     P0 = diagm([0.2, 0.2, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
#
#     ASSR_Model = StateSpaceModel(F, V, G, W, x0, P0)
#
# end

# function click_model(f::Number, fa::Number, fs::Number; x0::Vector=vcat(vec(0.1 * ones(6, 1)), vec(zeros(4, 1)), vec(20* ones(10, 1))),
#     W::Number=0.0001, V = diagm(vcat(vec(1e-10* ones(10, 1)), vec(0.001* ones(10, 1)))))
#
#     Δt = 1 / fs
#     ω = 2 * pi * f
#     ωa = 2 * pi * fa
#
#     ## Process transition and noise covariance
#
#     F = eye(20)
#     # V = diagm(vcat(vec(1e-10* ones(10, 1)), vec(0.001* ones(10, 1))))
#
#     ## Observation and noise covariance
#
#     names = (:G1, :G3, :G5, :G7, :G9)
#     freqs = (1,   2,   3,   4,   5)
#     for (name, freq) in zip(names, freqs)
#         @eval function ($name)(k)
#             cos($freq * $ωa * $Δt * k)
#         end
#     end
#
#     names = (:G2, :G4, :G6, :G8, :G10)
#     freqs = (1,   2,   3,   4,   5)
#     for (name, freq) in zip(names, freqs)
#         @eval function ($name)(k)
#             -sin($freq * $ωa * $Δt * k)
#         end
#     end
#
#     names = (:M1, :M3, :M5, :M7, :M9)
#     freqs = (1,   2,   3,   4,   5)
#     for (name, freq) in zip(names, freqs)
#         @eval function ($name)(k)
#             cos($freq * $ω * $Δt * k)
#         end
#     end
#
#     names = (:M2, :M4, :M6, :M8, :M10)
#     freqs = (1,   2,   3,   4,   5)
#     for (name, freq) in zip(names, freqs)
#         @eval function ($name)(k)
#             -sin($freq * $ω * $Δt * k)
#         end
#     end
#
#     G = reshape([M1, M2, M3, M4, M5, M6, M7, M8, M9, M10, G1, G2, G3, G4, G5, G6, G7, G8, G9, G10], 1, 20)
#     W = diagm([W])   # Noise Covariance
#
#     ## Inital guesses at state and error covariance
#
#     P0 = 20 * eye(20)
#
#     ASSR_Model = StateSpaceModel(F, V, G, W, x0, P0)
#
# end

# function ftest_model(f::Number, num_noise::Int, noise_freq::Number, fs::Number, W::Number;
#     x0::Vector = zeros(2 + (2 * num_noise), 1),  V = diagm(1e-10 * ones(2 + (2 * num_noise), 1)))
#
#     Δt = 1 / fs
#     ω = 2 * pi * f
#     ωa = 2 * pi * fa
#
#     ## Process transition and noise covariance
#
#     F = eye(20)
#     # V = diagm(vcat(vec(1e-10* ones(10, 1)), vec(0.001* ones(10, 1))))
#
#     ## Observation and noise covariance
#
#     names = (:G1, :G3, :G5, :G7, :G9)
#     freqs = (1,   2,   3,   4,   5)
#     for (name, freq) in zip(names, freqs)
#         @eval function ($name)(k)
#             cos($freq * $ωa * $Δt * k)
#         end
#     end
#
#     names = (:G2, :G4, :G6, :G8, :G10)
#     freqs = (1,   2,   3,   4,   5)
#     for (name, freq) in zip(names, freqs)
#         @eval function ($name)(k)
#             -sin($freq * $ωa * $Δt * k)
#         end
#     end
#
#     names = (:M1, :M3, :M5, :M7, :M9)
#     freqs = (1,   2,   3,   4,   5)
#     for (name, freq) in zip(names, freqs)
#         @eval function ($name)(k)
#             cos($freq * $ω * $Δt * k)
#         end
#     end
#
#     names = (:M2, :M4, :M6, :M8, :M10)
#     freqs = (1,   2,   3,   4,   5)
#     for (name, freq) in zip(names, freqs)
#         @eval function ($name)(k)
#             -sin($freq * $ω * $Δt * k)
#         end
#     end
#
#     G = reshape([M1, M2, M3, M4, M5, M6, M7, M8, M9, M10, G1, G2, G3, G4, G5, G6, G7, G8, G9, G10], 1, 20)
#     W = diagm([W])   # Noise Covariance
#
#     ## Inital guesses at state and error covariance
#
#     P0 = 20 * eye(20)
#
#     ASSR_Model = StateSpaceModel(F, V, G, W, x0, P0)
#
# end


text:
using Publish
using Pkg.Artifacts

using Metalhead

# override default theme
Publish.Themes.default() = artifact"flux-theme"

p = Publish.Project(Metalhead)

# serve documentation
serve(Metalhead)