text:
c1 = [-1; -4]
c2 = [-2; -3]

dim_x = length(c1)
dim_u = length(c2)

b = [-2; -3]

A1 = [1 -3;
     -1 -3]
A2 = [1 -2;
     -1 -1]

M = 1000;

# Loading the necessary packages
#-------------------------------
using JuMP
using GLPK
using LinearAlgebra
using Test

# Master Problem Description
# --------------------------

master_problem_model = Model(GLPK.Optimizer)

# Variable Definition
# ----------------------------------------------------------------
@variable(master_problem_model, 0 <= x[1:dim_x] <= 1e6, Int)
@variable(master_problem_model, t <= 1e6)

# Objective Setting
# -----------------
@objective(master_problem_model, Max, t)
global iter_num = 1

print(master_problem_model)

iter_num = 1

while(true)
    println("\n-----------------------")
    println("Iteration number = ", iter_num)
    println("-----------------------\n")
    println("The current master problem is")
    print(master_problem_model)

    optimize!(master_problem_model)

    t_status = termination_status(master_problem_model)
    p_status = primal_status(master_problem_model)

    if p_status == MOI.INFEASIBLE_POINT
        println("The problem is infeasible :-(")
        break
    end

    (fm_current, x_current) = if t_status == MOI.INFEASIBLE_OR_UNBOUNDED
        (M, M * ones(dim_x))
    elseif p_status == MOI.FEASIBLE_POINT
        (value(t), value.(x))
    else
        error("Unexpected status: $((t_status, p_status))")
    end

    println("Status of the master problem is ", t_status,
            "\nwith fm_current = ", fm_current,
            "\nx_current = ", x_current)

    sub_problem_model = Model(GLPK.Optimizer)

    c_sub = b - A1 * x_current

    local u = @variable(sub_problem_model, u[1:dim_u] >= 0)

    @constraint(sub_problem_model, constr_ref_subproblem[j = 1:size(A2, 2)], dot(A2[:, j], u) >= c2[j])
    # The second argument of @constraint macro,
    # constr_ref_subproblem[j=1:size(A2,2)] means that the j-th constraint is
    # referenced by constr_ref_subproblem[j].

    @objective(sub_problem_model, Min, dot(c1, x_current) + dot(c_sub, u))

    print("\nThe current subproblem model is \n", sub_problem_model)

    optimize!(sub_problem_model)

    t_status_sub = termination_status(sub_problem_model)
    p_status_sub = primal_status(sub_problem_model)

    fs_x_current = objective_value(sub_problem_model)

    u_current = value.(u)

    γ = dot(b, u_current)

    println("Status of the subproblem is ", t_status_sub,
        "\nwith fs_x_current = ", fs_x_current,
        "\nand fm_current = ", fm_current)

    if p_status_sub == MOI.FEASIBLE_POINT && fs_x_current == fm_current # we are done
        println("\n################################################")
        println("Optimal solution of the original problem found")
        println("The optimal objective value t is ", fm_current)
        println("The optimal x is ", x_current)
                println("The optimal v is ", dual.(constr_ref_subproblem))
        println("################################################\n")
        break
    end

    if p_status_sub == MOI.FEASIBLE_POINT && fs_x_current < fm_current
        println("\nThere is a suboptimal vertex, add the corresponding constraint")
        cv = A1' * u_current - c1
        @constraint(master_problem_model, t + dot(cv, x) <= γ)
        println("t + ", cv, "ᵀ x <= ", γ)
    end

    if t_status_sub == MOI.INFEASIBLE_OR_UNBOUNDED
        println("\nThere is an  extreme ray, adding the corresponding constraint")
        ce = A1'* u_current
        @constraint(master_problem_model, dot(ce, x) <= γ)
        println(ce, "ᵀ x <= ", γ)
    end

    global iter_num += 1
end

@test value(t) ≈ -4


text:
""" Computing in the category of finite sets and functions, and its skeleton.
"""
module FinSets
export FinOrd, FinOrdFunction, FinOrdFunc, FinOrdVector, force,
  terminal, product, equalizer, pullback, limit,
  initial, coproduct, coequalizer, pushout, colimit

using AutoHashEquals
using DataStructures: IntDisjointSets, union!, find_root

using ...GAT
using ...Theories: Category
import ...Theories: dom, codom, id, compose, ⋅, ∘,
  terminal, product, equalizer, initial, coproduct, coequalizer
using ..ShapeDiagrams

# Category of finite ordinals
#############################

""" Finite ordinal (natural number).

An object in the category of finite ordinals, which is the skeleton of the
category of finite sets.
"""
@auto_hash_equals struct FinOrd
  n::Int
end

""" Function between sets in the form of finite ordinals.

A morphism in the category of finite ordinals, which is the skeleton of the
category of finite sets. The function can be defined implicitly by an arbitrary
Julia function, in which case it is evaluated lazily, or explictly by a vector
of integers. In the latter case, the function (1↦1, 2↦3, 3↦2, 4↦3), for example,
is represented by the vector [1,3,2,3].
"""
abstract type FinOrdFunction end

FinOrdFunction(f, dom::FinOrd, codom::FinOrd) =
  FinOrdFunction(f, dom.n, codom.n)

""" Function in FinOrd defined by an arbitrary Julia function.

To be evaluated lazily unless forced.
"""
@auto_hash_equals struct FinOrdFunc <: FinOrdFunction
  func::Function
  dom::Int
  codom::Int
end
FinOrdFunction(f::Function, dom::Int, codom::Int) =
  FinOrdFunc(f, dom, codom)

(f::FinOrdFunc)(x) = f.func(x)

""" Function in FinOrd represented explicitly by a vector.
"""
@auto_hash_equals struct FinOrdVector{T<:AbstractVector{Int}} <: FinOrdFunction
  func::T
  codom::Int
end
FinOrdFunction(f::AbstractVector) = FinOrdVector(f, maximum(f))
FinOrdFunction(f::AbstractVector, codom::Int) = FinOrdVector(f, codom)

function FinOrdFunction(f::AbstractVector, dom::Int, codom::Int)
  @assert length(f) == dom
  FinOrdVector(f, codom)
end

(f::FinOrdVector)(x) = f.func[x]

""" Force evaluation of lazy function or relation.
"""
force(f::FinOrdFunction) = FinOrdVector(map(f, 1:dom(f).n), codom(f).n)
force(f::FinOrdVector) = f

Base.collect(f::FinOrdFunction) = force(f).func

""" Category of finite ordinals and functions.
"""
@instance Category(FinOrd, FinOrdFunction) begin
  dom(f::FinOrdFunction) = FinOrd(f.dom)
  codom(f::FinOrdFunction) = FinOrd(f.codom)

  id(A::FinOrd) = FinOrdFunction(identity, A, A)

  function compose(f::FinOrdFunction, g::FinOrdFunction)
    @assert codom(f) == dom(g)
    FinOrdFunction(compose_impl(f,g), dom(f), codom(g))
  end
end

dom(f::FinOrdVector) = FinOrd(length(f.func))
compose_impl(f::FinOrdFunction, g::FinOrdFunction) = g ∘ f
compose_impl(f::FinOrdVector, g::FinOrdVector) = g.func[f.func]

# Limits
########

terminal(::Type{FinOrd}) = FinOrd(1)

function product(A::FinOrd, B::FinOrd)
  m, n = A.n, B.n
  indices = CartesianIndices((m, n))
  π1 = FinOrdFunction(i -> indices[i][1], m*n, m)
  π2 = FinOrdFunction(i -> indices[i][2], m*n, n)
  Span(π1, π2)
end

function product(Xs::Vector{<:FinOrd})
  ns = Int[X.n for X in Xs]
  indices = CartesianIndices(tuple(ns...))
  apex = prod(ns)
  πs = [FinOrdFunction(i -> indices[i][j],apex,ns[j]) for j in 1:length(ns)]
  Cone(FinOrd(apex),πs)
end

function equalizer(f::FinOrdFunction, g::FinOrdFunction)
  @assert dom(f) == dom(g) && codom(f) == codom(g)
  m = dom(f).n
  FinOrdFunction(filter(i -> f(i) == g(i), 1:m), m)
end

function equalizer(fs::Vector{<:FinOrdFunction})
  @assert length(fs) >= 1
  f1 = fs[1]
  frest = fs[2:end]
  @assert all(dom(f) == dom(f1) && codom(f) == codom(f1) for f in frest)
  m = dom(f1).n
  FinOrdFunction(filter(i -> all(f1(i) == f(i) for f in frest), 1:m),m)
end

""" Pullback of cospan of functions between finite ordinals.

TODO: This logic is completely generic. Make it independent of FinOrd.
"""
function pullback(cospan::Cospan{<:FinOrdFunction,<:FinOrdFunction})
  f, g = left(cospan), right(cospan)
  prod = product(dom(f), dom(g))
  π1, π2 = left(prod), right(prod)
  eq = equalizer(π1⋅f, π2⋅g)
  Span(eq⋅π1, eq⋅π2)
end

function limit(d::Diagram{<:FinOrd, <:FinOrdFunction})
  p = product(d.obs)
  n = apex(p).n
  satisfy((s,t,g),x) = g(leg(p,s)(x)) == leg(p,t)(x)
  f = FinOrdFunction(filter(i -> all(satisfy(h,i) for h in d.homs), 1:n), n)
  Cone(dom(f),[compose(f,leg(p,i)) for i in 1:length(d.obs)])
end

# Colimits
##########

initial(::Type{FinOrd}) = FinOrd(0)

function coproduct(A::FinOrd, B::FinOrd)
  m, n = A.n, B.n
  ι1 = FinOrdFunction(1:m, m, m+n)
  ι2 = FinOrdFunction(m+1:m+n, n, m+n)
  Cospan(ι1, ι2)
end

function coproduct(Xs::Vector{<:FinOrd})
  ns = Int[X.n for X in Xs]
  base = sum(ns)
  offsets = [0,cumsum(ns)...]
  πs = [FinOrdFunction((1:ns[j]) .+ offsets[j],ns[j],base) for j in 1:length(ns)]
  Cocone(FinOrd(base),πs)
end

function coequalizer(f::FinOrdFunction, g::FinOrdFunction)
  @assert dom(f) == dom(g) && codom(f) == codom(g)
  m, n = dom(f).n, codom(f).n
  sets = IntDisjointSets(n)
  for i in 1:m
    union!(sets, f(i), g(i))
  end
  h = [ find_root(sets, i) for i in 1:n ]
  roots = unique!(sort(h))
  FinOrdFunction([ searchsortedfirst(roots, r) for r in h], length(roots))
end

function coequalizer(fs::Vector{<:FinOrdFunction})
  @assert length(fs) >= 1
  f1 = fs[1]
  frest = fs[2:end]
  @assert all(dom(f) == dom(f1) && codom(f) == codom(f1) for f in frest)
  m,n = dom(f1).n, codom(f1).n
  sets = IntDisjointSets(n)
  for i in 1:m
    for f in frest
      union!(sets, f1(i), f(i))
    end
  end
  h = [ find_root(sets, i) for i in 1:n ]
  roots = unique!(sort(h))
  FinOrdFunction([searchsortedfirst(roots, r) for r in h], length(roots))
end

""" Pushout of span of functions between finite ordinals.

TODO: This logic is completely generic. Make it independent of FinOrd.
"""
function pushout(span::Span{<:FinOrdFunction,<:FinOrdFunction})
  f, g = left(span), right(span)
  coprod = coproduct(codom(f), codom(g))
  ι1, ι2 = left(coprod), right(coprod)
  coeq = coequalizer(f⋅ι1, g⋅ι2)
  Cospan(ι1⋅coeq, ι2⋅coeq)
end

function colimit(d::Diagram{<:FinOrd, <:FinOrdFunction})
  cp = coproduct(d.obs)
  n = base(cp).n
  sets = IntDisjointSets(n)
  for (s,t,h) in d.homs
    for i in 1:d.obs[s].n
      union!(sets,leg(cp,s)(i),leg(cp,t)(h(i)))
    end
  end
  h = [ find_root(sets, i) for i in 1:n ]
  roots = unique!(sort(h))
  m = length(roots)
  f = FinOrdFunction([searchsortedfirst(roots, r) for r in h], m)
  Cocone(FinOrd(m),[compose(leg(cp,i),f) for i in 1:length(d.obs)])
end

end

ext: jl


text:
export associate!

function update_async!(net::HopfieldNetwork, activatefunc::Function)
	for i in randperm(length(net.state))
		net.state[i] = activatefunc(dot(net.weight[:, i], net.state))
	end
	return
end

function update_sync!(net::HopfieldNetwork, activatefunc::Function)
 	state = copy(net.state)
	for i=1:length(state)
		net.state[i] = activatefunc(dot(net.weight[:, i], state))
	end
	return
end

function associate!{T <: Number}(net::HopfieldNetwork{T}, pattern::AbstractVector{T}, activatefunc::Function, iterations::Integer = 1_000, sync::Bool = false)
	copy!(net.state, pattern)
	associate!(net,activatefunc,iterations,sync)
end

function associate!(net::HopfieldNetwork, activatefunc::Function, iterations::Integer = 1_000, sync::Bool = false)
	ite = iterations
	state = copy(net.state)

	if sync
		state2 = copy(net.state)
		for i in 1:iterations
			update_sync!(net,activatefunc)
			if state2 == net.state
				ite = i
				break
			end
			copy!(state2,state)
			copy!(state,net.state)
		end
	else
		for i in 1:iterations
			update_async!(net,activatefunc)
			if state == net.state
				ite = i
				break
			end
			copy!(state,net.state)
		end
	end

	return ite
end


text:
#####################################################################
#####                 STACK COEFFS PARSING                       ####
#####################################################################

function Qstack_coefs(mach)
    models = getfield(mach.model.Q̅, :models)
    coeffs = [p[2] for p in fitted_params(mach).Q̅.metalearner.coefs]
    return [m => c for (m,c) in zip(models, coeffs)]
end

repr_Qstack_coefs(mach; digits=2) = join([m => round(c, digits=digits) for (m, c) in Qstack_coefs(mach)], ", ")

#####################################################################
#####                     QUERY PARSING                          ####
#####################################################################

function parse_queries(queryfile::String)
    config = TOML.parsefile(queryfile)
    queries = Pair[]
    for (queryname, querydict) in config
        if lowercase(queryname) ∉ ("threshold", "snps")
            rsids = collect(keys(querydict))
            vals = [split(filter(x->!isspace(x), querydict[rsid]), "->") for rsid in rsids]
            rsids_symbols = Tuple(Symbol(x) for x in rsids)
            push!(queries, queryname => NamedTuple{rsids_symbols}(vals))
        end
    end
    return queries
end

function querystring(query)
    querystring_ = ""
    for (key, val) in pairs(query)
        querystring_ *= string(key, ": ", join(val, " -> "),  " & ")
    end
    return querystring_[1:length(querystring_) - 3]
end

#####################################################################
#####                 PHENOTYPE PARSING                          ####
#####################################################################

phenotypesnames(phenotypefile::String) =
    keys(only(CSV.File(phenotypefile, limit=1, drop=["eid"])))

phenotypes_list(phenotype_listfile::Nothing, done_phenotypes, allnames) =
    filter(x -> x ∉ done_phenotypes, allnames)

function phenotypes_list(phenotype_listfile::String, done_phenotypes, allnames)
    phenotypes_list = CSV.File(phenotype_listfile,
                               header=[:PHENOTYPES],
                               type=Symbol)
    filter(x -> (x ∈ Set(phenotypes_list.PHENOTYPES)) & (x ∉ done_phenotypes), allnames)
end

function init_or_retrieve_results(outfile, run_fn::typeof(PhenotypeTMLEEpistasis))
    if isfile(outfile)
        df = CSV.File(outfile, select=[:PHENOTYPE], type=Symbol) |> DataFrame
    else
        df = DataFrame(
            PHENOTYPE=Symbol[],
            QUERYNAME=String[],
            QUERYSTRING=String[],
            ESTIMATE=Float64[],
            PVALUE=Float64[],
            LOWER_BOUND=Float64[],
            UPPER_BOUND=Float64[],
            STD_ERROR=Float64[],
            QSTACK_COEFS=String[],
            NROWS=Int[],
            TCOUNTS=String[]
            )
        CSV.write(outfile, df)
    end
    return Set(df.PHENOTYPE)
end

function init_or_retrieve_results(outfile, run_fn::typeof(PhenotypeCrossValidation))
    if isfile(outfile)
        df = CSV.File(outfile, select=[:PHENOTYPE], type=Symbol) |> DataFrame
    else
        df = DataFrame(
            PHENOTYPE=Symbol[],
            VARIANTS=String[],
            Q_RESULTSTRING=String[],
            G_RESULTSTRING=String[],
            )
        CSV.write(outfile, df)
    end
    return Set(df.PHENOTYPE)
end

#####################################################################
#####                TMLE REPORT PARSING                         ####
#####################################################################

get_query_report(reports::NamedTuple, i) = reports
get_query_report(reports::Vector, i) = reports[i]

#####################################################################
#####                TREATMENT ENCODING                          ####
#####################################################################

function encode(T)
    n, p = size(T)
    if p == 1
        return T[:, 1]
    end

    joint_levels_it = Iterators.product((levels(T[!, name]) for name in names(T))...)
    encoding = Dict(Tuple(jl) => i for (i, jl) in enumerate(joint_levels_it))

    t_target = Vector{Int}(undef, n)
    for i in 1:n
        t_target[i] = encoding[values(T[i, :])]
    end

    return categorical(t_target; levels=collect(values(encoding)))

end

#####################################################################
#####                 TREATMENT COUNTS                           ####
#####################################################################

function TreatmentCountsRepr(T)
    countrepr = join(names(T), " & ") * ": "
    for row in eachrow(combine(groupby(T, names(T)), nrow))
        countrepr *= join(row, " ") * " | "
    end
    return countrepr[1:end-3]
end


text:
## Simple diet demo for real.

using DataFrames
using CSV
using JuMP
using Clp

## These are very simple codes
function from_csv(filename)
    #=
    Given the name of a csv file whose first line are headers,
    return a list of dictionaries, one for each row of the file,
    whose keys are the header for that column.
    =#
    content = CSV.read(filename)
    headers = names(content)

    m = size(content)[1]
    table = []
    for i = 1 : m
        line = content[i,:]
        push!(table, Dict(zip(headers, line)))
    end
    return table
end

function clean_up_dict(dict)
    #=
    Given dictionary, clean up the dict (set Missing data to be 0.0)
    =#
    for pair in dict
        if typeof(pair[2]) == Missing
            dict[pair[1]] = 0.0
        end
    end
    return dict
end
function clean_up_table(table)
    #=
    Given table, clean up the table (set Missing data to be 0.0)
    =#
    m = size(table)[1]
    for i = 1 : m
        clean_up_dict(table[i])
    end
end

function create_variable_dict(table)
    #=
    mapping the food_name with variables
    =#
    m = size(table)[1]
    variable_dict = Dict()
    @variable(myModel, 0 <= x[1:m] <= 10)
    for i = 1 : m
        variable_dict[table[i][:description]] = x[i]
    end
    return variable_dict
end

function foods_for_nutrient(nutrient_name)
    #=
    Get the nutrient_amount for variables
    =#
    nutrients_list = []
    var_list = []
    for row in nutrients_table
        var = variable_dict[row[:description]]
        nutrient_amount = row[Symbol(nutrient_name)]
        if nutrient_amount > 0
            push!(nutrients_list, nutrient_amount)
            push!(var_list, var)
        end
    end
    return nutrients_list, var_list
end

function create_percent_constraint(nutrient_name, lower, upper, calories_per_gram)
    #=
    Compute the constraint that says the total consumed nutrient
    must be between 'lower' and 'upper' percent of the total calories.
    =#
    calories_lower_list, var_calories_lower = foods_for_nutrient(calories_name)
    calories_upper_list, var_calories_upper = foods_for_nutrient(calories_name)
    nutrient_list, var_nutrient_list = foods_for_nutrient(nutrient_name)
    m = size(calories_lower_list)[1]
    n = size(calories_upper_list)[1]
    t = size(nutrient_list)[1]

    @constraint(myModel, sum(calories_lower_list[j] * var_calories_lower[j] * lower/100.0 for j = 1 : m ) <= sum(nutrient_list[i] * var_nutrient_list[i] * calories_per_gram for i = 1 : t))
    @constraint(myModel, sum(nutrient_list[i] * var_nutrient_list[i] * calories_per_gram for i = 1 : t) <= sum(calories_upper_list[j] * var_calories_upper[j] * upper/100.0 for j = 1 : n ))
end

function create_constraint(nutrient_name, lower, upper)
    #=
     Each constraint is a lower and upper bound on the
     sum of all food variables, scaled by how much of the
     relevant nutrient is in that food.
    =#
    if lower == 0.0
        return
    end
    if haskey(percent_constraints, nutrient_name)
        calories_per_gram = percent_constraints[nutrient_name]["calories_per_gram"]
        create_percent_constraint(nutrient_name, lower, upper, calories_per_gram)
        return
    end

    foods_list, var_foods = foods_for_nutrient(nutrient_name)
    m = size(foods_list)[1]
    @constraint(myModel, lower <= sum(foods_list[i] * var_foods[i] for i = 1:m))

    if upper == 0.0
        return
    end
    @constraint(myModel, sum(foods_list[i] * var_foods[i] for i = 1:m) <= upper)

end

function create_constraints(table)
    #=
    Create constraints
    =#
    constraint_bounds = Dict()
    for row in table
        nutrient = row[:nutrient]
        lower_bound = row[:lower_bound]
        upper_bound = row[:upper_bound]
        constraint_bounds[nutrient] = (lower_bound, upper_bound)
        create_constraint(nutrient, lower_bound, upper_bound)
    end
    return constraint_bounds
end

#MODEL CONSTRUCTION
#------------------------------------------------------------------------------------------

myModel = Model(with_optimizer(Clp.Optimizer))
# solver=GLPKsolverLP() means that to solve the optimization problem we will use GLPK solver.

calories_name = "energy (kcal)"

nutrients_table = from_csv("nutrients-kun.csv")
constraints_table = from_csv("diet-kun.csv")

clean_up_table(nutrients_table)
clean_up_table(constraints_table)

# next step is to create_variable_dict()
variable_dict = create_variable_dict(nutrients_table)

#next is to create constraints
percent_constraints = Dict("total fat (g)"=> Dict("calories_per_gram"=> 9))
constraint_bounds = create_constraints(constraints_table)

# now create objective
var_list = []
calories_list = []
for row in nutrients_table
    name = row[:description]
    var = variable_dict[name]
    push!(var_list, var)
    calories_in_food = row[Symbol(calories_name)]
    push!(calories_list, calories_in_food)
end
num = size(var_list)[1]
@objective(myModel, Min, sum(var_list[i] * calories_list[i] for i = 1:num))

# Solve it!
optimize!(myModel)

##
for item in keys(variable_dict)
    var = variable_dict[item]
    println(item, ";", value(var))
end


text:
# using OpenAIGym
# import OpenAIGym
using Flux
using JuliaRL
using Random
import LearnBase.IntervalSet

# rand(rng, int_set::)
Base.rand(rng::AbstractRNG, s::IntervalSet{T}) where T <: AbstractVector = Float64[rand(rng) * (s.hi[i] - s.lo[i]) + s.lo[i] for i=1:length(s)]

function test_gym(num_episodes, num_steps, to_render::Bool=false)

    rng = MersenneTwister(1)
    env = GymEnv(:BipedalWalker, :v2; seed=1)
    println(typeof(env))
    for i ∈ 1:num_episodes
        T = 0
        R = 0.0
        action_set = get_actions(env)
        st = start!(env)
        for i in 1:num_steps
            action = rand(rng, action_set)
            step!(env, action)
            if to_render
                render(env)
            end
            R += get_reward(env)
            T += 1
        end

        @info("Episode $i finished after $T steps. Total reward: $R")
    end

    # You can explicitly close the connection for the python environment here. But will be taken care of by the garbage collector eventually.
    close(env)

end

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
"turns top level arrays into dicts"
function arrays_to_dicts!(data::Dict{String,<:Any})
    # update lookup structure
    for (k,v) in data
        if isa(v, Array) && length(v) > 0 && isa(v[1], Dict)
            #println("updating $(k)")
            dict = Dict{String,Any}()
            for (i,item) in enumerate(v)
                if haskey(item, "index")
                    key = string(item["index"])
                else
                    key = string(i)
                end

                if !(haskey(dict, key))
                    dict[key] = item
                else
                    @warn "skipping component $(item["index"]) from the $(k) table because a component with the same id already exists"
                end
            end
            data[k] = dict
        end
    end
end

"takes a row from a matrix and assigns the values names and types"
function row_to_typed_dict(row_data, columns)
    dict_data = Dict{String,Any}()
    for (i,v) in enumerate(row_data)
        if i <= length(columns)
            name, typ = columns[i]
            dict_data[name] = check_type(typ, v)
        else
            dict_data["col_$(i)"] = v
        end
    end
    return dict_data
end

"takes a row from a matrix and assigns the values names"
function row_to_dict(row_data, columns)
    dict_data = Dict{String,Any}()
    for (i,v) in enumerate(row_data)
        if i <= length(columns)
            dict_data[columns[i]] = v
        else
            dict_data["col_$(i)"] = v
        end
    end
    return dict_data
end

row_to_dict(row_data) = row_to_dict(row_data, [])


text:
using Measurements

@testset "ccm89" begin

    # NOTE: Test is only to precision of 0.015 because there is a discrepancy
    # of 0.014 for the B band wavelength of unknown origin (and up to 0.002 in
    # other bands).
    #
    # Note that a and b can be obtained with:
    # b = ccm89(wave, 0.)
    # a = ccm89(wave, 1.) - b
    #
    # These differ from the values tablulated in the original paper.
    # Could be due to floating point errors in the original paper?

    # U, B, V, R, I, J, H, K band effective wavelengths from CCM '89 table 3
    x_inv_microns = [10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.6, 4.0, 3.0, 2.78, 2.27, 1.82, 1.43, 1.11, 0.8, 0.63, 0.46]
    wave = 1e4 ./ x_inv_microns

    # A(lambda)/A(V) for different R_V from Table 3 of CCM '89
    ref_values = Dict(3.1 => [5.23835484, 4.13406452, 3.33685933, 2.77962453, 2.52195399, 2.84252644, 3.18598916, 2.31531711, 1.64254927, 1.56880904, 1.32257836, 1.0, 0.75125994, 0.4780346, 0.28206957, 0.19200814, 0.11572348],
        2.0 => [9.407, 7.3065, 5.76223881, 4.60825807, 4.01559036, 4.43845534, 4.93952892, 3.39275574, 2.068771, 1.9075018, 1.49999733, 1.0, 0.68650255, 0.36750326, 0.21678862, 0.14757062, 0.08894094],
        3.0 => [5.491, 4.32633333, 3.48385202, 2.8904508, 2.6124774, 2.9392494, 3.2922643, 2.38061642, 1.66838089, 1.58933588, 1.33333103, 1.0, 0.74733525, 0.47133573, 0.27811315, 0.18931496, 0.11410029],
        4.0 => [3.533, 2.83625, 2.34465863, 2.03154717, 1.91092092, 2.18964643, 2.46863199, 1.87454675, 1.46818583, 1.43025292, 1.24999788, 1.0, 0.7777516, 0.52325196, 0.30877542, 0.21018713, 0.12667997],
        5.0 => [2.3582, 1.9422, 1.66114259, 1.51620499, 1.48998704, 1.73988465, 1.97445261, 1.57090496, 1.3480688, 1.33480314, 1.19999799, 1.0, 0.79600141, 0.5544017, 0.32717278, 0.22271044, 0.13422778],
        6.0 => [1.575, 1.34616667, 1.20546523, 1.17264354, 1.20936444, 1.44004346, 1.64499968, 1.36847709, 1.26799077, 1.27116996, 1.16666472, 1.0, 0.80816794, 0.5751682, 0.33943769, 0.23105931, 0.13925965])

    # Test defaults
    @test @inferred(broadcast(ccm89, wave)) ≈ ref_values[3.1] rtol = 0.016
    @test_deprecated ccm89(wave)

    for rv in [2.0, 3.0, 3.1, 4.0, 5.0, 6.0]
        output = @inferred broadcast(ccm89, wave, rv)
        @test output ≈ ref_values[rv] rtol = 0.016

        # Test deprecated array syntax
        @test_deprecated ccm89(wave, rv)

        bad_waves = [100, 4e4]
        @test @inferred(broadcast(ccm89, bad_waves, rv)) == zeros(length(bad_waves))

        # uncertainties
        noise = randn(length(wave)) .* 0.01
        wave_unc = wave .± noise
        reddening = @inferred broadcast(ccm89, wave_unc, rv)
        @test Measurements.value.(reddening) ≈ ref_values[rv] rtol = 0.016

        # Unitful
        wave_u = wave * u"angstrom"
        reddening = @inferred broadcast(ccm89, wave_u, rv)
        @test eltype(reddening) <: Gain
        @test ustrip.(reddening) ≈ ref_values[rv] rtol = 0.016
    end
end

@testset "od94" begin
    # NOTE: The tabulated values go to 0.001, but the test is only for matching
    # at the 0.005 level, because there is currently a discrepancy up to 0.0047
    # of unknown origin.

    # Tests od94() at Rv = 3.1 against the widely used values tabulated in
    # Schlegel, Finkbeiner and Davis (1998)
    # http://adsabs.harvard.edu/abs/1998ApJ...500..525S

    # This is tested by evaluating the extinction curve at a (given)
    # effective wavelength, since these effective wavelengths:
    # "... represent(s) that wavelength on the extinction curve
    # with the same extinction as the full passband."

    # The test does not include UKIRT L' (which, at 3.8 microns) is
    # beyond the range of wavelengths allowed by the function
    # or the APM b_J filter which is defined in a non-standard way.

    # The SFD98 tabulated values go to 1e-3, so we should be able to match at
    # that level.

    wave = [3372., 4404., 5428., 6509., 8090.,
            3683., 4393., 5519., 6602., 8046.,
            12660., 16732., 22152.,
            5244., 6707., 7985., 9055.,
            6993.,
            3502., 4676., 4127.,
            4861., 5479.,
            3546., 4925., 6335., 7799., 9294.,
            3047., 4711., 5498.,
            6042., 7068., 8066.,
            4814., 6571., 8183.]

    ref_values = [1.664, 1.321, 1.015, 0.819, 0.594,
                1.521, 1.324, 0.992, 0.807, 0.601,
                0.276, 0.176, 0.112,
                1.065, 0.793, 0.610, 0.472,
                0.755,
                1.602, 1.240, 1.394,
                1.182, 1.004,
                1.579, 1.161, 0.843, 0.639, 0.453,
                1.791, 1.229, 0.996,
                0.885, 0.746, 0.597,
                1.197, 0.811, 0.580]

    reddening = @inferred broadcast(od94, wave, 3.1)
    @test reddening ≈ ref_values rtol = 0.016

    # Test deprecated array syntax
    @test_deprecated od94(wave)

    # Test out of bounds
    bad_waves = [100, 4e4]
    @test @inferred(broadcast(od94, bad_waves, 3.1)) == zeros(length(bad_waves))

    @testset "uncertainties" begin
        noise = randn(length(wave)) .* 10
        wave_unc = wave .± noise
        reddening = @inferred broadcast(od94, wave_unc, 3.1)
        @test Measurements.value.(reddening) ≈ ref_values rtol = 0.016
    end
    @testset "unitful" begin
        wave_u = wave * u"angstrom"
        reddening = @inferred broadcast(od94, wave_u, 3.1)
        @test eltype(reddening) <: Gain
        @test ustrip.(reddening) ≈ ref_values rtol = 0.016
    end
end

@testset "cal00" begin

    refwave = [3090.90909091,  4561.61616162,  6872.72727273,  9604.04040404,
    14646.46464646, 14646.46464646, 15486.86868687, 18218.18181818,
    20529.29292929, 20949.49494949]

    refmag = Dict(3.1 => [ 1.88010678,  1.27137591,  0.70513192,  0.33600273,  0.01622915,
    0.01622915, -0.01682162, -0.10317769, -0.15830055, -0.16701622],
        2.4 => [ 2.13680458,  1.35052722,  0.61912873,  0.14233687, -0.27070401,
    -0.27070401, -0.3133946 , -0.42493784, -0.49613821, -0.50739595])

    # test defaults
    @test cal00.(refwave) ≈ refmag[3.1]

    for rv in [2.4, 3.1]
        reddening = @inferred broadcast(cal00, refwave, rv)
        @test reddening ≈ refmag[rv]

        bad_waves = [
            1e2,
            3e4,
            4e4
        ]
        @test @inferred(broadcast(cal00, bad_waves, rv)) == zeros(length(bad_waves))

        # Uncertainties
        noise = randn(length(refwave)) .* 10
        wave_unc = refwave .± noise
        reddening = @inferred broadcast(cal00, refwave, rv)
        @test Measurements.value.(reddening) ≈ refmag[rv]

        # Unitful
        wave_u = refwave * u"angstrom"
        reddening = @inferred broadcast(cal00, wave_u, rv)
        @test eltype(reddening) <: Gain
        @test ustrip.(reddening) ≈ refmag[rv]
    end
end


text:
struct TwoLoopVars{T, TR, M}
    d::T # search direction
    S::TR # change in x
    Y::TR # change in successive gradients
    α::M
    ρ::M
end
function TwoLoopVars(x, memory)
    d = copy(x)

    S = [copy(x) for i = 1:memory]
    Y = [copy(x) for i = 1:memory]
    α = similar(x, memory)
    ρ = similar(x, memory)
    TwoLoopVars(d, S, Y, α, ρ)
end
lbfgs_vars(method, x) = TwoLoopVars(x, method.memory)
function solve(problem::OptimizationProblem, x0, approach::LineSearch{<:LBFGS{<:Inverse, <:TwoLoop}, <:Any, <:QNScaling}, options::OptimizationOptions, cache=preallocate_qn_caches(mstyle(problem), x0))
    solve(problem, (x0, nothing), approach, options, cache)
end
function solve(problem::OptimizationProblem, s0::Tuple, approach::LineSearch{<:LBFGS{<:Inverse, <:TwoLoop}, <:Any, <:QNScaling}, options::OptimizationOptions, cache=preallocate_qn_caches(mstyle(problem), first(s0)))
    _solve(InPlace(), problem, s0, approach, options, cache)
end

function _solve(mstyle, problem::OptimizationProblem, s0::Tuple, approach::LineSearch{<:LBFGS{<:Inverse, <:TwoLoop}, <:Any, <:QNScaling}, options::OptimizationOptions, cache)
    t0 = time()

    #==============
         Setup
    ==============#
    x0, B0 = s0
    T = eltype(x0)
    # Remove B from all of this
    objvars = prepare_variables(problem, approach, x0, copy(x0), B0)

    qnvars = lbfgs_vars(modelscheme(approach), objvars.x)
    ∇f0 = norm(objvars.∇fz, Inf)
    f0 = objvars.fz
    P = initial_preconditioner(approach, x0)

    #========================
         First iteration
    ========================#
    objvars, qnvars, P = iterate(mstyle, 1, qnvars, objvars, P, approach, problem, problem, options)
    iter = 1
    # Check for gradient convergence
    is_converged = converged(approach, objvars, ∇f0, options)
    while iter <= options.maxiter && !any(is_converged)
        iter += 1
        # take a step and update approximation
        objvars, qnvars, P = iterate(mstyle, iter, qnvars, objvars, P, approach, problem, problem, options, false)
        # Check for gradient convergence
        is_converged = converged(approach, objvars, ∇f0, options)
    end
    x, fx, ∇fx, z, fz, ∇fz, B, Pg = objvars
    return ConvergenceInfo(approach, (P=P,B=B, ρs=norm(x.-z), ρx=norm(x), minimizer=z, fx=fx, minimum=fz, ∇fz=∇fz, f0=f0, ∇f0=∇f0, iter=iter, time=time()-t0), options)
end
function iterate(mstyle::InPlace, iter::Integer, qnvars::TwoLoopVars, objvars, P, approach::LineSearch{<:LBFGS{<:Inverse, <:TwoLoop}, <:Any, <:QNScaling}, prob::OptimizationProblem, problem::OptimizationProblem, options::OptimizationOptions, is_first=nothing)
    # split up the approach into the hessian approximation scheme and line search
    scheme, linesearch, K = modelscheme(approach), algorithm(approach), approach.scaling # make grabber and get a better name than K
    current_memory = min(iter-1, scheme.memory)
    x, fx, ∇fx, z, fz, ∇fz, B, Pg = objvars
    Tf = typeof(fz)

    # Move nexts into currs
    fx = fz
    copyto!(x, z)
    copyto!(∇fx, ∇fz)

    # Update preconditioner
    P = update_preconditioner(scheme, x, P) # returns nothing????
    # Update current gradient and calculate the search direction
    d = find_direction!(scheme, copy(∇fz), qnvars, current_memory, K, P) # solve Bd = -∇fx
    φ = _lineobjective(mstyle, problem, ∇fz, z, x, d, fx, real(dot(∇fx, d))) # real is needed to convert complex dots to actually being real

    # Perform line search along d
    α, f_α, ls_success = find_steplength(mstyle, linesearch, φ, Tf(1))

    s = current_memory == length(qnvars.S) ? qnvars.S[1] : qnvars.S[1+current_memory]
    @. s = α * d
    z = retract(problem, z, x, s)

    # Update approximation
    fz, ∇fz, qnvars = update_obj!(problem, qnvars, α, x, ∇fx, z, ∇fz, current_memory, scheme)
    return (x=x, fx=fx, ∇fx=∇fx, z=z, fz=fz, ∇fz=∇fz, B=nothing, Pg=nothing), qnvars, P
end

function iterate(mstyle::OutOfPlace, cache, objvars, P, approach::LineSearch{<:LBFGS{<:Inverse, <:TwoLoop}, <:Any, <:QNScaling}, prob::OptimizationProblem, problem::OptimizationProblem, options::OptimizationOptions, is_first=nothing)
    # split up the approach into the hessian approximation scheme and line search
    scheme, linesearch = modelscheme(approach), algorithm(approach)
    # Move nexts into currs
    x, fx, ∇fx, z, fz, ∇fz, B, Pg = objvars
    Tf = typeof(fz)
    fx = fz
    x = copy(z)
    ∇fx = copy(∇fz)
    # Update preconditioner
    P = update_preconditioner(scheme, x, P)
    # Update current gradient and calculate the search direction
    d = find_direction(B, P, ∇fx, scheme) # solve Bd = -∇fx
    φ = _lineobjective(mstyle, problem, problem, ∇fz, z, x, d, fx, dot(∇fx, d))

    # # Perform line search along d
    α, f_α, ls_success = find_steplength(linesearch, φ, Tf(1))

    # # Calculate final step vector and update the state
    s = @. α * d
    z = retract(problem, z, x, s)

    # Update approximation
    fz, ∇fz, B = update_obj(problem, s, ∇fx, z, ∇fz, B, scheme, is_first)

    return x, fx, ∇fx, z, fz, ∇fz, B, P
end


text:
using PyTensorFlow
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

if Sys.islinux()
py"""
import tensorflow as tf
libPoissonOp = tf.load_op_library('../Ops/Poisson/build/libPoissonOp.so')
@tf.custom_gradient
def poisson_op(coef,g,h,rhograv,index):
    p = libPoissonOp.poisson_op(coef,g,h,rhograv,index)
    def grad(dy):
        return libPoissonOp.poisson_op_grad(dy, p, coef, g, h, rhograv, index)
    return p, grad
"""
elseif Sys.isapple()
py"""
import tensorflow as tf
libPoissonOp = tf.load_op_library('../Ops/Poisson/build/libPoissonOp.dylib')
@tf.custom_gradient
def poisson_op(coef,g,h,rhograv,index):
    p = libPoissonOp.poisson_op(coef,g,h,rhograv,index)
    def grad(dy):
        return libPoissonOp.poisson_op_grad(dy, p, coef, g, h, rhograv, index)
    return p, grad
"""
elseif Sys.iswindows()
py"""
import tensorflow as tf
libPoissonOp = tf.load_op_library('../Ops/Poisson/build/libPoissonOp.dll')
@tf.custom_gradient
def poisson_op(coef,g,h,rhograv,index):
    p = libPoissonOp.poisson_op(coef,g,h,rhograv,index)
    def grad(dy):
        return libPoissonOp.poisson_op_grad(dy, p, coef, g, h, rhograv, index)
    return p, grad
"""
end

poisson_op = py"poisson_op"

if Sys.islinux()
py"""
import tensorflow as tf
libUpwpsOp = tf.load_op_library('../Ops/Upwps/build/libUpwpsOp.so')
@tf.custom_gradient
def upwps_op(permi,mobi,src,funcref,h,rhograv,index):
    pres = libUpwpsOp.upwps_op(permi,mobi,src,funcref,h,rhograv,index)
    def grad(dy):
        return libUpwpsOp.upwps_op_grad(dy, pres, permi,mobi,src,funcref,h,rhograv,index)
    return pres, grad
"""
elseif Sys.isapple()
py"""
import tensorflow as tf
libUpwpsOp = tf.load_op_library('../Ops/Upwps/build/libUpwpsOp.dylib')
@tf.custom_gradient
def upwps_op(permi,mobi,src,funcref,h,rhograv,index):
    pres = libUpwpsOp.upwps_op(permi,mobi,src,funcref,h,rhograv,index)
    def grad(dy):
        return libUpwpsOp.upwps_op_grad(dy, pres, permi,mobi,src,funcref,h,rhograv,index)
    return pres, grad
"""
elseif Sys.iswindows()
py"""
import tensorflow as tf
libUpwpsOp = tf.load_op_library('../Ops/Upwps/build/libUpwpsOp.dll')
@tf.custom_gradient
def upwps_op(permi,mobi,src,funcref,h,rhograv,index):
    pres = libUpwpsOp.upwps_op(permi,mobi,src,funcref,h,rhograv,index)
    def grad(dy):
        return libUpwpsOp.upwps_op_grad(dy, pres, permi,mobi,src,funcref,h,rhograv,index)
    return pres, grad
"""
end

upwps_op = py"upwps_op"


text:
# Step 1: Home page prompts login.
function index(::HTTP.Request)
    links = map(collect(PROVIDERS)) do p
        link = p.first
        name = p.second.name
        """
Log in to $name
"""
    end
    return html(join(links, "
"))
end


text:
module event

const _NC_EVENT_TIMEOUT_MSEC = 1
const _NC_EVENT_FILE = 2
const _NC_EVENT_FILE_READABLE = 2
#  Not supported yet... *
const _NC_EVENT_FILE_WRITABLE = 4
const _NC_EVENT_FILE_EXCEPTION = 8

immutable _nc_event
  _type::Cint
  type data
    timeout_msec::Clong
    type fev
        flags::Cuint
        fd::Cint
        result::Cuint
    end
  end
end

immutable _nc_eventlist
  count::Cint
  result_flags::Cint
  events::Ptr{Array{_nc_event, 1}}
end

end


text:
"""
    GA_reproduction(pa::AbstractVector{T},
                    pb::AbstractVector{T},
                    bounds;
                    η_cr = 20,
                    η_m  = 15,
                    p_cr = 0.9,
                    p_m  = 0.1)

Crate two solutions by applying SBX to parents pa and pb and polynomial mutation to
offspring. Return two vectors.
"""
function GA_reproduction(pa::AbstractVector{T},
        pb::AbstractVector{T},
        bounds;
        η_cr = 20,
        η_m  = 15,
        p_cr = 0.9,
        p_m  = 0.1
    ) where T <: AbstractFloat

    # crossover
    c1, c2 = SBX_crossover(pa, pb, bounds, η_cr, p_cr)

    # mutation
    polynomial_mutation!(c1, bounds,η_m, p_m)
    polynomial_mutation!(c2, bounds,η_m, p_m)

    # rapair solutions if necesary
    reset_to_violated_bounds!(c1, bounds)
    reset_to_violated_bounds!(c2, bounds)

    return c1, c2
end

"""
    GA_reproduction_half(pa::AbstractVector{T},
                    pb::AbstractVector{T},
                    bounds;
                    η_cr = 20,
                    η_m  = 15,
                    p_cr = 0.9,
                    p_m  = 0.1)

Same that `GA_reproduction` but only returns one offspring.
"""
function GA_reproduction_half(pa::AbstractVector{T},
        pb::AbstractVector{T},
        bounds;
        η_cr = 20,
        η_m  = 15,
        p_cr = 0.9,
        p_m  = 0.1
    ) where T <: AbstractFloat

    # crossover
    _, c2 = SBX_crossover(pa, pb, bounds, η_cr, p_cr)

    # mutation
    polynomial_mutation!(c, bounds,η_m, p_m)

    # rapair solution if necesary
    reset_to_violated_bounds!(c, bounds)

    return c
end


text:
using Turing, Random, Test
using Turing.Inference: split_var_str

dir = splitdir(splitdir(pathof(Turing))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")

@testset "io.jl" begin
    @testset "threaded sampling" begin
        # Only test threading if 1.3+.
        if VERSION > v"1.2"
            # Test that chains with the same seed will sample identically.
            chain1 = psample(Random.seed!(5), gdemo_default, HMC(0.1, 7), 1000, 4)
            chain2 = psample(Random.seed!(5), gdemo_default, HMC(0.1, 7), 1000, 4)
            @test all(chain1.value .== chain2.value)
            check_gdemo(chain1)

            # Smoke test for default psample call.
            chain = psample(gdemo_default, HMC(0.1, 7), 1000, 4)
            check_gdemo(chain)
        end
    end
    @testset "chain save/resume" begin
        Random.seed!(1234)

        alg1 = HMCDA(1000, 0.65, 0.15)
        alg2 = PG(20)
        alg3 = Gibbs(PG(30, :s), HMCDA(500, 0.65, 0.05, :m))

        chn1 = sample(gdemo_default, alg1, 3000; save_state=true)
        check_gdemo(chn1)

        chn1_resumed = Turing.Inference.resume(chn1, 1000)
        check_gdemo(chn1_resumed)

        chn1_contd = sample(gdemo_default, alg1, 1000; resume_from=chn1)
        check_gdemo(chn1_contd)

        chn1_contd2 = sample(gdemo_default, alg1, 1000; resume_from=chn1, reuse_spl_n=1000)
        check_gdemo(chn1_contd2)

        chn2 = sample(gdemo_default, alg2, 500; save_state=true)
        check_gdemo(chn2)

        chn2_contd = sample(gdemo_default, alg2, 500; resume_from=chn2)
        check_gdemo(chn2_contd)

        chn3 = sample(gdemo_default, alg3, 500; save_state=true)
        check_gdemo(chn3)

        chn3_contd = sample(gdemo_default, alg3, 500; resume_from=chn3)
        check_gdemo(chn3_contd)
    end
    @testset "split var string" begin
        var_str = "x"
        sym, inds = split_var_str(var_str)
        @test sym == "x"
        @test inds == Vector{String}[]

        var_str = "x[1,1][2,3]"
        sym, inds = split_var_str(var_str)
        @test sym == "x"
        @test inds[1] == ["1", "1"]
        @test inds[2] == ["2", "3"]

        var_str = "x[Colon(),1][2,Colon()]"
        sym, inds = split_var_str(var_str)
        @test sym == "x"
        @test inds[1] == ["Colon()", "1"]
        @test inds[2] == ["2", "Colon()"]

        var_str = "x[2:3,1][2,1:2]"
        sym, inds = split_var_str(var_str)
        @test sym == "x"
        @test inds[1] == ["2:3", "1"]
        @test inds[2] == ["2", "1:2"]

        var_str = "x[2:3,2:3][[1,2],[1,2]]"
        sym, inds = split_var_str(var_str)
        @test sym == "x"
        @test inds[1] == ["2:3", "2:3"]
        @test inds[2] == ["[1,2]", "[1,2]"]
    end
end


text:
################################
###    Summary statistics    ###
################################
"""
	poissonFixation(observedValues,λds, λdn)

Divergence sampling from Poisson distribution. The expected neutral and selected fixations are subset through their relative expected rates ([`fixNeut`](@ref), [`fixNegB`](@ref), [`fixPosSim`](@ref)). Empirical values are used are used to simulate the locus *L* along a branch of time *T* from which the expected *Ds* and *Dn* raw count estimated given the mutation rate (``\\mu``). Random number generation is used to subset samples arbitrarily given the success rate ``\\lambda`` in the distribution.

```math
\\mathbb{E}[D_N] = X \\in Poisson\\left(\\lambda = D \\times \\left[\\frac{\\mathbb{E}[D_+] + \\mathbb{E}[D_-]}{\\mathbb{E}[D_+] + \\mathbb{E}[D_-] + \\mathbb{E}[D_0]}\\right]\\right)
```
```math
\\mathbb{E}[D_S] = X \\in Poisson\\left(\\lambda = D \\times \\left[\\frac{\\mathbb{E}[D_0]}{\\mathbb{E}[D_+] + \\mathbb{E}[D_-] + \\mathbb{E}[D_0]}\\right]\\right)
```
# Arguments
 - `observedValues::Array`: Array containing the total observed divergence.
 - ` λds::Float64`: expected neutral fixations rate.
 - ` λdn::Float64`: expected selected fixations rate.
# Returns
 - `Array{Int64,1}` containing the expected count of neutral and selected fixations.

"""
function poissonFixation(;observedValues::Array, λds::Array, λdn::Array,λweak::Array,λstrong::Array)

	ds = @. λds / (λds + λdn) * observedValues
	dn = @. λdn / (λds + λdn) * observedValues
	dweak = @. λweak / (λds + λdn) * observedValues
	dstrong = @. λstrong / (λds + λdn) * observedValues

	sampledDs     = pois_rand.(ds)
	sampledDn     = pois_rand.(dn)
	sampledWeak   = pois_rand.(dweak)
	sampledStrong = pois_rand.(dstrong)

	alphas = @. [sampledWeak/sampledDn sampledStrong/sampledDn (sampledWeak+sampledStrong)/sampledDn]

	out = alphas,sampledDn, sampledDs
	return out
end

"""
	poissonPolymorphism(observedValues,λps,λpn)

Polymorphism sampling from Poisson distributions. The total expected neutral and selected polimorphism are subset through the relative expected rates at the frequency spectrum ([`fixNeut`](@ref), [`DiscSFSNeutDown`](@ref),). Empirical sfs are used to simulate the locus *L* along a branch of time *T* from which the expected *Ps* and *Pn* raw count are estimated given the mutation rate (``\\mu``). Random number generation is used to subset samples arbitrarily from the whole sfs given each frequency success rate ``\\lambda`` in the distribution.

The success rate managing the Poisson distribution by the observed count each frequency.  We considered both sampling variance and process variance is affecting the number of variable alleles we sample from SFS. This variance arises from the random mutation-fixation process along the branch. To incorporate this variance we do one sample per frequency-bin and use the total sampled variation and the SFS for the summary statistics.

```math
\\mathbb{E}[P_N] = \\sum_{x=0}^{x=1} X \\in Poisson\\left(\\lambda = SFS_{(x)} \\times \\left[\\frac{\\mathbb{E}[P_{+(x)}] + \\mathbb{E}[P_{-(x)}]}{\\mathbb{E}[P_{+(x)}] + \\mathbb{E}[P_{-(x)}] + \\mathbb{E}[P_{0(x)}]}\\right]\\right)
```

```math
\\mathbb{E}[P_S] = \\sum_{x=0}^{x=1} X \\in Poisson\\left(\\lambda = SFS_{(x)} \\times \\left[\\frac{\\mathbb{E}[P_{0(x)}]}{\\mathbb{E}[P_{+(x)}] + \\mathbb{E}[P_{-(x)}] + \\mathbb{E}[P_{0(x)}]}\\right]\\right)
```

# Arguments
 - `observedValues::Array{Int64,1}`: Array containing the total observed divergence.
 - ` λps::Array{Float64,1} `: expected neutral site frequency spectrum rate.
 - ` λpn::Array{Float64,1} `: expected selected site frequency spectrum rate.
# Returns
 - `Array{Int64,2}` containing the expected total count of neutral and selected polymorphism.

"""
function poissonPolymorphism(;observedValues::Array, λps::Array, λpn::Array)

	# Neutral λ;
	λ1 = @. λps / (λps + λpn) * observedValues
	# Selected λ;
	λ2 = @. λpn / (λps + λpn) * observedValues

	# Replace negative and 0 rates values. Strange behaviour depending on θ and Γ values
	# Relative rates output NaN values due to 0 divisons.

	replace!(λ1,NaN=>1)
	replace!(λ2,NaN=>1)

	sampledPs = pois_rand.(λ1)
	sampledPn = pois_rand.(λ2)

	return (sampledPn, sampledPs)
end

"""
	sampledAlpha(observedValues,λds, λdn)

Ouput the expected values from the Poisson sampling process. Please check [`poissonFixation`](@ref) and [`poissonPolymorphism`](@ref) to understand the samplingn process. α(x) is estimated through the expected values of Dn, Ds, Pn and Ps.

# Arguments
 - `param::parameters`: Array containing the total observed divergence.
 - `d::Array`: observed divergence.
 - `afs::Array`: observed polymorphism.
 - ` λdiv::Array{Float64,2}`: expected fixations rate.
 - ` λdiv::Array{Float64,2}`: expected site frequency spectrum rates.
# Returns
αS,expDn,expDs,expPn,expPs,ssAlpha
 - `Array{Int64,2}` containing α(x) values.
 - `Array{Int64,1}` expected non-synonymous divergence.
 - `Array{Int64,1}` expected synonymous divergence.
 - `Array{Int64,1}` expected non-synonymous polymorphism.
 - `Array{Int64,1}` expected synonymous polymorphism.
 - `Array{Int64,1}` expected synonymous polymorphism.
 - `Array{Int64,1}` expected synonymous polymorphism.
 - `Array{Int64,2}` containing α(x) binned values.

"""
function sampledAlpha(;d::Array,afs::Array,observed::Matrix,λdiv::Array,λpol::Array,dac::Array)

	## Outputs
	alphas, expDn, expDs = poissonFixation(observedValues=d,λds=λdiv[1],λdn=λdiv[2],λweak=λdiv[3],λstrong=λdiv[4])
	expPn, expPs         = poissonPolymorphism(observedValues=afs,λps=λpol[1],λpn=λpol[2])

	θ_pn = @. observed[1,1] / (expPn[1,:]);
	θ_ps = @. observed[1,2] / (expPs[1,:]);

	r_pn = ones(size(expPn));
	r_ps = ones(size(expPs));
	for j = 2:size(dac,1)
		r_pn[j,:] .= @. (observed[j,1]) / (expPn[j,:]);
		r_pn[j,:] = @. r_pn[j,:] / θ_pn;

		r_ps[j,:] .= @. (observed[j,2]) / (expPs[j,:]);
		r_ps[j,:] .= @. r_ps[j,:] / θ_ps;

	end

	expPn_rj = zeros(size(expPn))
	expPs_rj = zeros(size(expPs))

	expPn_rj[1,:] = @. expPn[1,:] * θ_pn;
	expPs_rj[1,:] = @. expPs[1,:] * θ_ps;
	for j = 2:size(dac,1)
		expPn_rj[j,:] = @. expPn[j,:] * θ_pn * r_pn[j,:];
		expPs_rj[j,:] = @. expPs[j,:] * θ_ps * r_ps[j,:];
	end

	## Alpha from expected values. Used as summary statistics
	αS = @. round(1 - ((expDs/expDn) * (expPn/expPs)'),digits=5)
	αS_v2 = @. round(1 - ((expDs/expDn) * (expPn_rj/expPs_rj)'),digits=5)

	return αS,alphas,expDn,expDs,expPn,expPs
end

"""
	samplingFromRates(gammaL,gammaH,pposL,pposH,observedData,nopos)
"""
function samplingFromRates(m::Array,s::Array,d::Array,nt::Array,sl::Array,x::Array,observed::Matrix,dac::Array)

	ds             = x[:,1]
	dn             = x[:,2]
	dweak          = x[:,3]
	dstrong        = x[:,4]
	gn             = abs.(m[:,4])
	sh             = round.(m[:,end-1],digits=5)

	alxSummStat, alphasDiv, expectedDn, expectedDs, expectedPn, expectedPs = sampledAlpha(d=d,afs=s,observed=observed,λdiv=[ds,dn,dweak,dstrong],λpol=[permutedims(nt),permutedims(sl)],dac=dac)
	expectedValues = hcat(round.(alphasDiv,digits=5),gn,sh,alxSummStat)

	return(expectedValues)
end

"""
	summaryStatsFromRates(param::parameters,rates::JLD2.JLDFile,analysisFolder::String,summstatSize::Int64,replicas::Int64,bootstrap::Bool)

Estimate summary statistics using observed data and analytical rates. *analysisFolder* will check for the SFS and divergence file and will be used to output summary statistics

# Arguments
 - `param::parameters` : Mutable structure containing the models
 - `rates::JLD2.JLDFile` : HDF5 containing solved models, fixation and polymorphic rates
 - `analysisFolder::String` : Folder containing the SFS and divergence files. It will be used to output the observed data and summary estatistics.
 - `summstatSize::Int64` : Number of summary statistics
 - `replicas::Int64` : Number of bootstrap replicas
 - `bootstrap::Bool` : Boolean to perform or not bootstrapping
# Output
 - Obserded data and summary statistics to ABC inference
"""
function summaryStatsFromRates(;param::parameters,h5file::JLD2.JLDFile,analysisFolder::String,summstatSize::Int64,replicas::Int64,bootstrap::Bool)

	#Opening files
	sFile   = filter(x -> occursin("sfs",x), readdir(analysisFolder,join=true));
	dFile   = filter(x -> occursin("div",x), readdir(analysisFolder,join=true));

	sfs,d,α,observed = openSfsDiv(sFile,dFile,param.dac,replicas,bootstrap);

	#Open rates
	tmp    = h5file[string(param.N) * "/" * string(param.n)]

	#Subset index
	idx    = sample(1:size(tmp["models"],1),summstatSize,replace=false)
	models = Array(view(tmp["models"],idx,:));
	dsdn   = Array(view(tmp["dsdn"],idx,:));

	# Filtering polymorphic rate by dac
	n    = hcat(map(x -> view(tmp["neut"][x],:),param.dac)...);
	s    = hcat(map(x -> view(tmp["sel"][x],:),param.dac)...);
	neut = Array(view(n,idx,:));
	sel  = Array(view(s,idx,:));

	#Making summaries
	f(x,y,o,m=models,nt=neut,sl=sel,d=dsdn,dac=param.dac) = samplingFromRates(m,x,y,nt,sl,d,o,dac);
	expectedValues = f.(sfs,d,observed);

	w(x,name) = CSV.write(name,DataFrame(x,:auto),delim='\t',header=false);

	# Controling outlier cases
	fltInf(e) = replace!(e, -Inf=>NaN)
	expectedValues = fltInf.(expectedValues)
	fltNan(e) = e[vec(.!any(isnan.(e),dims=2)),:]
	expectedValues = fltInf.(expectedValues)
	expectedValues = fltNan.(expectedValues)
	expectedValues = map(x -> x[(x[:,3] .> 0 ) .& (x[:,3] .<1 ),:],expectedValues)

	# Writting ABCreg input
	#=w(vcat(α...), analysisFolder * "/alphas.txt");
	w(expectedValues, analysisFolder * "/summstat.txt");=#
	progress_map(w,α, analysisFolder * "/alphas_" .* string.(1:size(sFile,1)) .* ".txt";progress= Progress(size(sFile,1),desc="Writting α "));
	progress_pmap(w,expectedValues,  analysisFolder * "/summstat_" .* string.(1:size(sFile,1)) .* ".txt";progress= Progress(size(sFile,1),desc="Writting summaries "));

	return(expectedValues)
end


text:
## Data structure for the simulations

using Parameters

@with_kw struct Data

  # particle data
  eps :: Float64  = 5.0
  sig :: Float64  = 0.5
  mass :: Float64 = 1.0

  # box side and cutoff information
  cutoff :: Float64 = 3.
  side :: Float64   = 150.

  # number of particles
  N :: Int64 = 400

  # Data for computation of ljpair
  eps4 :: Float64  = 4*eps
  sig6 :: Float64  = sig^6
  sig12 :: Float64 = sig^12

end


text:
using Glob: glob
using Random: shuffle
using Distributed: procs, addprocs
using Distributed: @everywhere

const learning_rate = .01
const hm_epochs     = 50
const batch_size    = 2

const in_size       = 52
const layers        = [52]
const out_size      = 52

if length(procs()) <= 2
    hm_procs = batch_size > Sys.CPU_THREADS ?
        Sys.CPU_THREADS-1 : batch_size-1
    addprocs(hm_procs) ;println("$(length(procs())) cores running.")
end
@everywhere include("gru_api(char_to_char).jl")
@everywhere include("utils.jl")

create_model_definitions(layers)

const model = Model(mk_layers(in_size, layers, out_size)...)

const data = import_data(glob("class*.txt")[1])

for i in 1:hm_epochs
    print("epoch: $i ")
    train!(model, shuffle(data), learning_rate, batch_size)
    end
end

save_model(model)

print("Hit enter to continue..")
chomp(readline())


text:
function __init__cifar100()
    DEPNAME = "CIFAR100"

    register(DataDep(
        DEPNAME,
        """
        Dataset: The CIFAR-100 dataset
        Authors: Alex Krizhevsky, Vinod Nair, Geoffrey Hinton
        Website: https://www.cs.toronto.edu/~kriz/cifar.html
        Reference: https://www.cs.toronto.edu/~kriz/learning-features-2009-TR.pdf

        [Krizhevsky, 2009]
            Alex Krizhevsky.
            "Learning Multiple Layers of Features from Tiny Images",
            Tech Report, 2009.

        The CIFAR-100 dataset is a labeled subsets of the 80
        million tiny images dataset. It consists of 60000
        32x32 colour images in 100 classes. Specifically, it
        has 100 classes containing 600 images each. There are
        500 training images and 100 testing images per class.
        The 100 classes in the CIFAR-100 are grouped into 20
        superclasses. Each image comes with a "fine" label
        (the class to which it belongs) and a "coarse" label
        (the superclass to which it belongs).

        The compressed archive file that contains the
        complete dataset is available for download at the
        offical website linked above; specifically the binary
        version for C programs. Note that using the data
        responsibly and respecting copyright remains your
        responsibility. The authors of CIFAR-10 aren't really
        explicit about any terms of use, so please read the
        website to make sure you want to download the
        dataset.
        """,
        "https://www.cs.toronto.edu/~kriz/cifar-100-binary.tar.gz",
        "58a81ae192c23a4be8b1804d68e518ed807d710a4eb253b1f2a199162a40d8ec",
        post_fetch_method = file -> (run(BinDeps.unpack_cmd(file, dirname(file), ".gz", ".tar")); rm(file))
    ))
end

"""
    CIFAR100(; Tx=Float32, split=:train, dir=nothing)
    CIFAR100([Tx, split])

The CIFAR100 dataset is a labeled subsets of the 80
million tiny images dataset. It consists of 60000
32x32 colour images in 10 classes, with 6000 images
per class.

Return the CIFAR-100 **trainset** labels (coarse and fine)
corresponding to the given `indices` as a tuple of two `Int` or
two `Vector{Int}`. The variables returned are the coarse label(s)
(`Yc`) and the fine label(s) (`Yf`) respectively.

# Arguments

$ARGUMENTS_SUPERVISED_ARRAY
- `split`: selects the data partition. Can take the values `:train:` or `:test`.

# Fields

$FIELDS_SUPERVISED_ARRAY
- `split`.

# Methods

$METHODS_SUPERVISED_ARRAY
- [`convert2image`](@ref) converts features to `RGB` images.

# Examples

```julia-repl
julia> dataset = CIFAR100()
CIFAR100:
  metadata    =>    Dict{String, Any} with 3 entries
  split       =>    :train
  features    =>    32×32×3×50000 Array{Float32, 4}
  targets     =>    (coarse = "50000-element Vector{Int64}", fine = "50000-element Vector{Int64}")

julia> dataset[1:5].targets
(coarse = [11, 15, 4, 14, 1], fine = [19, 29, 0, 11, 1])

julia> X, y = dataset[];

julia> dataset.metadata
Dict{String, Any} with 3 entries:
  "n_observations"     => 50000
  "class_names_coarse" => ["aquatic_mammals", "fish", "flowers", "food_containers", "fruit_and_vegetables", "household_electrical_devices", "household_furniture", "insects", "large_carnivores", "large_man-made_…
  "class_names_fine"   => ["apple", "aquarium_fish", "baby", "bear", "beaver", "bed", "bee", "beetle", "bicycle", "bottle"  …  "train", "trout", "tulip", "turtle", "wardrobe", "whale", "willow_tree", "wolf", "w…
```
"""
struct CIFAR100 <: SupervisedDataset
    metadata::Dict{String, Any}
    split::Symbol
    features::Array{<:Any, 4}
    targets::NamedTuple{(:coarse, :fine), Tuple{Vector{Int}, Vector{Int}}}
end

CIFAR100(; split=:train, Tx=Float32, dir=nothing) = CIFAR100(Tx, split; dir)

function CIFAR100(Tx::Type, split::Symbol=:train; dir=nothing)
    DEPNAME = "CIFAR100"
    TRAINSET_FILENAME = joinpath("cifar-100-binary", "train.bin")
    TESTSET_FILENAME  = joinpath("cifar-100-binary", "test.bin")
    COARSE_FILENAME = joinpath("cifar-100-binary", "coarse_label_names.txt")
    FINE_FILENAME = joinpath("cifar-100-binary", "fine_label_names.txt")
    TRAINSET_SIZE = 50_000
    TESTSET_SIZE  = 10_000

    @assert split ∈ (:train, :test)

    if split == :train
        file_path = datafile(DEPNAME, TRAINSET_FILENAME, dir)
        images, labels_c, labels_f = CIFAR100Reader.readdata(file_path, TRAINSET_SIZE)
    else
        file_path = datafile(DEPNAME, TESTSET_FILENAME, dir)
        images, labels_c, labels_f = CIFAR100Reader.readdata(file_path, TESTSET_SIZE)
    end
    features = bytes_to_type(Tx, images)
    targets = (coarse = labels_c, fine = labels_f)

    metadata = Dict{String, Any}()
    metadata["class_names_coarse"] = readlines(datafile(DEPNAME, COARSE_FILENAME, dir))
    metadata["class_names_fine"] = readlines(datafile(DEPNAME, FINE_FILENAME, dir))
    metadata["n_observations"] = size(features)[end]

    return CIFAR100(metadata, split, features, targets)
end

convert2image(::Type{<:CIFAR100}, x) = convert2image(CIFAR10, x)

# DEPRECATED INTERFACE, REMOVE IN v0.7 (or 0.6.x)
function Base.getproperty(::Type{CIFAR100}, s::Symbol)
    if s == :traintensor
        @warn "CIFAR100.traintensor() is deprecated, use `CIFAR100(split=:train).features` instead." maxlog=2
        traintensor(T::Type=N0f8; kws...) = traintensor(T, :; kws...)
        traintensor(i; kws...) = traintensor(N0f8, i; kws...)
        function traintensor(T::Type, i; dir=nothing)
            CIFAR100(; split=:train, Tx=T, dir)[i][1]
        end
        return traintensor
    elseif s == :testtensor
        @warn "CIFAR100.testtensor() is deprecated, use `CIFAR100(split=:test).features` instead."  maxlog=2
        testtensor(T::Type=N0f8; kws...) = testtensor(T, :; kws...)
        testtensor(i; kws...) = testtensor(N0f8, i; kws...)
        function testtensor(T::Type, i; dir=nothing)
            CIFAR100(; split=:test, Tx=T, dir)[i][1]
        end
        return testtensor
    elseif s == :trainlabels
        @warn "CIFAR100.trainlabels() is deprecated, use `CIFAR100(split=:train).targets` instead."  maxlog=2
        trainlabels(; kws...) = trainlabels(:; kws...)
        function trainlabels(i; dir=nothing)
            yc, yf = CIFAR100(; split=:train, dir)[i][2]
            yc, yf
        end
        return trainlabels
    elseif s == :testlabels
        @warn "CIFAR100.testlabels() is deprecated, use `CIFAR100(split=:test).targets` instead." maxlog=2
        testlabels(; kws...) = testlabels(:; kws...)
        function testlabels(i; dir=nothing)
            yc, yf = CIFAR100(; split=:test, dir)[i][2]
            yc, yf
        end
        return testlabels
    elseif s == :traindata
        @warn "CIFAR100.traindata() is deprecated, use `CIFAR100(split=:train)[]` instead." maxlog=2
        traindata(T::Type=N0f8; kws...) = traindata(T, :; kws...)
        traindata(i; kws...) = traindata(N0f8, i; kws...)
        function traindata(T::Type, i; dir=nothing)
            x, (yc, yf) = CIFAR100(; split=:train, Tx=T, dir)[i]
            x, yc, yf
        end
        return traindata
    elseif s == :testdata
        @warn "CIFAR100.testdata() is deprecated, use `CIFAR100(split=:test)[]` instead."  maxlog=2
        testdata(T::Type=N0f8; kws...) = testdata(T, :; kws...)
        testdata(i; kws...) = testdata(N0f8, i; kws...)
        function testdata(T::Type, i; dir=nothing)
            x, (yc, yf) = CIFAR100(; split=:test, Tx=T, dir)[i]
            x, yc, yf
        end
        return testdata
    elseif s == :convert2image
        @warn "CIFAR100.convert2image(x) is deprecated, use `convert2image(CIFAR100, x)` instead"
        return x -> convert2image(CIFAR100, x)
    elseif s == :classnames_fine
        @warn "CIFAR100.classnames_fine() is deprecated, use `CIFAR100().metadata[\"class_names_fine\"]` instead"
        return () -> CIFAR100().metadata["class_names_fine"]
    elseif s == :classnames_coarse
        @warn "CIFAR100.classnames_coarse() is deprecated, use `CIFAR100().metadata[\"class_names_coarse\"]` instead"
        return () -> CIFAR100().metadata["class_names_coarse"]
    else
        return getfield(CIFAR100, s)
    end
end


text:
# Pairwise Scalar & Vector Operation  ======================================================

@inline base_initiate(::BaseFunction, ::Type{T}) where {T} = zero(T)
@inline base_return(::BaseFunction, s::T) where {T} = s

function base_evaluate(f::BaseFunction, x::T, y::T) where {T<:AbstractFloat}
    base_return(f, base_aggregate(f, base_initiate(f,T), x, y))
end

# Note: no checks, assumes length(x) == length(y) >= 1
function unsafe_base_evaluate(
        f::BaseFunction,
        x::AbstractArray{T},
        y::AbstractArray{T}
    ) where {T<:AbstractFloat}
    s = base_initiate(f, T)
    @simd for I in eachindex(x, y)
        @inbounds xi = x[I]
        @inbounds yi = y[I]
        s = base_aggregate(f, s, xi, yi)
    end
    base_return(f, s)
end

function base_evaluate(
        f::BaseFunction,
        x::AbstractArray{T},
        y::AbstractArray{T}
    ) where {T<:AbstractFloat}
    if (n = length(x)) != length(y)
        throw(DimensionMismatch("Arrays x and y must have the same length."))
    elseif n == 0
        throw(DimensionMismatch("Arrays x and y must be at least of length 1."))
    end
    unsafe_base_evaluate(f, x, y)
end

# Pairwise Matrix Calculation  =============================================================

for orientation in (:row, :col)

    row_oriented = orientation == :row
    dim_obs      = row_oriented ? 1 : 2
    dim_param    = row_oriented ? 2 : 1

    @eval begin

        @inline function subvector(
                ::Val{$(Meta.quot(orientation))},
                X::AbstractMatrix,
                i::Integer
            )
            $(row_oriented ? :(view(X, i, :)) : :(view(X, :, i)))
        end

        @inline function allocate_basematrix(
                ::Val{$(Meta.quot(orientation))},
                X::AbstractMatrix{T}
            ) where {T<:AbstractFloat}
            Array{T}(undef, size(X,$dim_obs), size(X,$dim_obs))
        end

        @inline function allocate_basematrix(
                ::Val{$(Meta.quot(orientation))},
                X::AbstractMatrix{T},
                Y::AbstractMatrix{T}
            ) where {T<:AbstractFloat}
            Array{T}(undef, size(X,$dim_obs), size(Y,$dim_obs))
        end

        function checkdimensions(
                ::Val{$(Meta.quot(orientation))},
                P::Matrix,
                X::AbstractMatrix
            )
            n = size(P,1)
            if size(P,2) != n
                throw(DimensionMismatch("Pairwise matrix P must be square"))
            elseif size(X, $dim_obs) != n
                errorstring = string("Dimensions of P must match dimension ", $dim_obs, " of X")
                throw(DimensionMismatch(errorstring))
            end
            return n
        end

        function checkdimensions(
                ::Val{$(Meta.quot(orientation))},
                P::Matrix,
                X::AbstractMatrix,
                Y::AbstractMatrix
            )
            n = size(X, $dim_obs)
            m = size(Y, $dim_obs)
            if n != size(P,1)
                errorstring = string("Dimension 1 of P must match dimension ", $dim_obs, "of X")
                throw(DimensionMismatch(errorstring))
            elseif m != size(P,2)
                errorstring = string("Dimension 2 of P must match dimension ", $dim_obs, "of Y")
                throw(DimensionMismatch(errorstring))
            end
            if size(X, $dim_param) != size(Y, $dim_param)
                errorstring = string("Dimension $($dim_param) of X and Y must match")
                throw(DimensionMismatch(errorstring))
            end
            return (n, m)
        end
    end
end

function basematrix!(
        σ::Orientation,
        P::Matrix{T},
        f::BaseFunction,
        X::AbstractMatrix{T},
        symmetrize::Bool
    ) where {T<:AbstractFloat}
    n = checkdimensions(σ, P, X)
    for j = 1:n
        xj = subvector(σ, X, j)
        for i = 1:j
            xi = subvector(σ, X, i)
            @inbounds P[i,j] = unsafe_base_evaluate(f, xi, xj)
        end
    end
    symmetrize ? LinearAlgebra.copytri!(P, 'U', false) : P
end

function basematrix!(
        σ::Orientation,
        P::Matrix{T},
        f::BaseFunction,
        X::AbstractMatrix{T},
        Y::AbstractMatrix{T},
    ) where {T<:AbstractFloat}
    n, m = checkdimensions(σ, P, X, Y)
    for j = 1:m
        yj = subvector(σ, Y, j)
        for i = 1:n
            xi = subvector(σ, X, i)
            @inbounds P[i,j] = unsafe_base_evaluate(f, xi, yj)
        end
    end
    P
end

# ScalarProduct using BLAS/Built-In methods ================================================

@inline function basematrix!(
        σ::Orientation,
        P::Matrix{T},
        f::ScalarProduct,
        X::Matrix{T},
        symmetrize::Bool
    ) where {T<:LinearAlgebra.BLAS.BlasReal}
    gramian!(σ, P, X, symmetrize)
end

@inline function basematrix!(
        σ::Orientation,
        P::Matrix{T},
        f::ScalarProduct,
        X::Matrix{T},
        Y::Matrix{T},
    ) where {T<:LinearAlgebra.BLAS.BlasReal}
    gramian!(σ, P, X, Y)
end

# SquaredDistance using BLAS/Built-In methods ==============================================

function squared_distance!(
        G::Matrix{T},
        xᵀx::Vector{T},
        symmetrize::Bool
    ) where {T<:AbstractFloat}
    if !((n = size(G,1)) == size(G,2))
        throw(DimensionMismatch("Gramian matrix must be square."))
    end
    if length(xᵀx) != n
        throw(DimensionMismatch("Length of xᵀx must match order of G"))
    end
    @inbounds for j = 1:n
        xᵀx_j = xᵀx[j]
        for i = 1:(j-1)
            G[i,j] = (xᵀx[i] + xᵀx_j) - 2G[i,j]
        end
        G[j,j] = zero(T)
    end
    symmetrize ? LinearAlgebra.copytri!(G, 'U') : G
end

function squared_distance!(
        G::Matrix{T},
        xᵀx::Vector{T},
        yᵀy::Vector{T}
    ) where {T<:AbstractFloat}
    n, m = size(G)
    if length(xᵀx) != n
        throw(DimensionMismatch("Length of xᵀx must match rows of G"))
    elseif length(yᵀy) != m
        throw(DimensionMismatch("Length of yᵀy must match columns of G"))
    end
    @inbounds for j = 1:m
        yᵀy_j = yᵀy[j]
        for i = 1:n
            G[i,j] = (xᵀx[i] + yᵀy_j) - 2G[i,j]
        end
    end
    G
end

function fix_negatives!(
        σ::Orientation,
        D::Matrix{T},
        X::Matrix{T},
        symmetrize::Bool,
        ϵ::T=zero(T)
    ) where {T<:AbstractFloat}
    if !((n = size(D,1)) == size(D,2))
        throw(DimensionMismatch("Distance matrix must be square."))
    end
    for j = 1:n
        xj = subvector(σ, X, j)
        for i = 1:(j-1)
            if D[i,j] < ϵ
                xi = subvector(σ, X, i)
                D[i,j] = unsafe_base_evaluate(SquaredEuclidean(), xi, xj)
            end
        end
    end
    symmetrize ? LinearAlgebra.copytri!(D, 'U') : D
end

function fix_negatives!(
        σ::Orientation,
        D::Matrix{T},
        X::Matrix{T},
        Y::Matrix{T},
        ϵ::T=zero(T)
    ) where {T<:AbstractFloat}
    n, m = size(D)
    for j = 1:m
        yj = subvector(σ, Y, j)
        for i = 1:n
            if D[i,j] < ϵ
                xi = subvector(σ, X, i)
                D[i,j] = unsafe_base_evaluate(SquaredEuclidean(), xi, yj)
            end
        end
    end
    D
end

function basematrix!(
        σ::Orientation,
        P::Matrix{T},
        f::SquaredEuclidean,
        X::Matrix{T},
        symmetrize::Bool
    ) where {T<:LinearAlgebra.BLAS.BlasReal}
    gramian!(σ, P, X, false)
    xᵀx = dotvectors(σ, X)
    squared_distance!(P, xᵀx, false)
    fix_negatives!(σ, P, X, true)
    symmetrize ? LinearAlgebra.copytri!(P, 'U') : P
end

function basematrix!(
        σ::Orientation,
        P::Matrix{T},
        f::SquaredEuclidean,
        X::Matrix{T},
        Y::Matrix{T},
    ) where {T<:LinearAlgebra.BLAS.BlasReal}
    gramian!(σ, P, X, Y)
    xᵀx = dotvectors(σ, X)
    yᵀy = dotvectors(σ, Y)
    squared_distance!(P, xᵀx, yᵀy)
    fix_negatives!(σ, P, X, Y)
end


text:
using JLD2 # save, load

# "Compressed sparse column format document-term matrix"
# struct DocumentTermMatrix
#     "ordered term vector"
#     terms::Vector{String}
#     "term to column index mapping"
#     termidx::Dict{String,Int}
#     "document-term matrix"
#     matrix::SparseMatrixCSC{Int,Int}
# end

abstract type AbstractDocument end

"Generic document metadata container"
struct GenericDocument <: AbstractDocument
    filename::String
    metadata::Dict{String,Any}
end

"Simple document container identified with filename"
struct SimpleDocument <: AbstractDocument
    filename::String
end

"Extralinguistic document metadata for complexity project"
struct ComplexityDocument <: AbstractDocument
    filename::String
    id::String
    lang::String
    corpus::String
    author::String
    title::String
    year::Int8
    genre::String
    brow::String
    narrative::String
end

tokens(d::AbstractDocument) = String[
    token
    for line in eachline(d.filename)
    for token in split(line, " ")
    if token != "" && token != "<EOS>" && token != "<PGB>" && token != "<CHB>" && !all(ispunct, token)
]

function sentences(d::AbstractDocument)
    sens = Vector{Vector{String}}()
    tokens = Vector{String}()
    for line in eachline(d.filename)
        if !isempty(tokens) && line == "<EOS>" || line == "<PGB>" || line == "<CHB>"
            push!(sens, tokens)
            tokens = Vector{String}()
        else
            push!(tokens, line)
        end
    end
    sens
end

struct CorpusMetadata
    metadata::Dict{String,Any}
end

mutable struct Corpus{T <: AbstractDocument}
    metadata::CorpusMetadata
    documents::Vector{T}
    ntokens::UInt
    index::Dict{String,UInt}
    inverseindex::Dict{String,Vector{UInt}}
    dtm::SparseMatrixCSC{UInt,UInt}
end

Base.start(c::Corpus) = c.documents[1]
Base.next(c::Corpus, i) = (c.documents[i], i+1)
function Base.getindex(c::Corpus, i::Int)
    1 <= i <= c.ntokens || throw(BoundsError(c, i))
    c.documents[i]
end
Base.getindex(c::Corpus, i::Number) = c[convert(Int, i)]
Base.getindex(c::Corpus, I) = [c[i] for i in I]
Base.endof(c::Corpus) = c.ntokens
Base.getindex(c::Corpus, docname::AbstractString) = get(c.document_names, docname)

"""
    Corpus

"""
function Corpus(cmeta::CorpusMetadata, documents::Vector{T}) where T <: AbstractDocument
    ntokens = 0
    index = Dict{String,UInt}
    inverseindex = Dict{String,Vector{UInt}}

    for doc in documents
        ts = tokens(doc.filename)
        ntokens += length(ts)

    end

    Corpus(
        metadata,
        documents,
        ntokens,
        Dict{String,Int}(),
        Dict{String,Vector{Int}}()
    )
end

metadata(d::AbstractDocument) = d.metadata


text:
#=
rdb-26_002.jl
based on
https://gitlab.com/fMRI/toppe/blob/master/+toppe/+utils/read_rdb_hdr.m
=#

"""
`s = read_rdb_hdr_26_002(fid)`

Read GE raw (RDB) header for MRI scan file

in
- `fid::IOStream` from open(file, "r")

out
- `s::NamedTuple` with header values, accessed by s.key

Matlab version is:
Copyright (c) 2012 by General Electric Company. All rights reserved.

2019-05-22 Julia version by Jeff Fessler
"""
function read_rdb_hdr_26_002(fid::IOStream)
	seek(fid, 0)

	read1 = (T::DataType) -> read(fid, T)

	fread = (n::Integer, T::DataType) ->
		begin
			data = Array{T}(undef, n)
			read!(fid, data)
			return data
		end

	freadc = (n::Integer) -> rstrip(String(fread(n, UInt8)))

	return ( # NamedTuple
		rdbm_rev = read1(Float32),
		off_data = read1(Int32),
		off_per_pass = read1(Int32),
		off_unlock_raw = read1(Int32),
		off_data_acq_tab = read1(Int32),
		off_nex_tab = read1(Int32),
		off_nex_abort_tab = read1(Int32),
		off_tool = read1(Int32),
		off_exam = read1(Int32),
		off_series = read1(Int32),
		off_image = read1(Int32),
		off_ps = read1(Int32),
		off_grad_data = read1(Int32),
		off_CTT_data = read1(Int32),
		off_spare1 = read1(Int32),
		off_spare2 = read1(Int32),
		off_spare3 = read1(Int32),
		off_spare4 = read1(Int32),
		off_spare5 = read1(Int32),
		off_spare6 = read1(Int32),
		run_int = read1(Int32),
		scan_seq = read1(Int16),
		run_char = freadc(6),
		scan_date = freadc(10),
		scan_time = freadc(8),
		logo = freadc(10),
		file_contents = read1(Int16),
		lock_mode = read1(Int16),
		dacq_ctrl = read1(Int16),
		recon_ctrl = read1(UInt16),
		exec_ctrl = read1(UInt16),
		scan_type = read1(Int16),
		data_collect_type = read1(Int16),
		data_format = read1(Int16),
		recon = read1(Int16),
		datacq = read1(Int16),
		npasses = read1(Int16),
		npomp = read1(Int16),
		nslices = read1(UInt16),
		nechoes = read1(Int16),
		navs = read1(Int16),
		nframes = read1(Int16),
		baseline_views = read1(Int16),
		hnover = read1(Int16),
		frame_size = read1(UInt16),
		point_size = read1(Int16),
		vquant = read1(Int16),
		cheart = read1(Int16),
		ctr = read1(Float32),
		ctrr = read1(Float32),
		initpass = read1(Int16),
		incrpass = read1(Int16),
		method_ctrl = read1(Int16),
		da_xres = read1(UInt16),
		da_yres = read1(Int16),
		rc_xres = read1(Int16),
		rc_yres = read1(Int16),
		im_size = read1(Int16),
		rc_zres = read1(Int32),
		fermi_radius = read1(Float32),
		fermi_width = read1(Float32),
		fermi_ecc = read1(Float32),
		clip_min = read1(Float32),
		clip_max = read1(Float32),
		default_offset = read1(Float32),
		xoff = read1(Float32),
		yoff = read1(Float32),
		nwin = read1(Float32),
		ntran = read1(Float32),
		scalei = read1(Float32),
		scaleq = read1(Float32),
		rotation = read1(Int16),
		transpose = read1(Int16),
		kissoff_views = read1(Int16),
		slblank = read1(Int16),
		gradcoil = read1(Int16),
		ddaover = read1(Int16),
		sarr = read1(Int16),
		fd_tr = read1(Int16),
		fd_te = read1(Int16),
		fd_ctrl = read1(Int16),
		algor_num = read1(Int16),
		fd_df_dec = read1(Int16),
		dab = fread(8, Int16),
		user0 = read1(Float32),
		user1 = read1(Float32),
		user2 = read1(Float32),
		user3 = read1(Float32),
		user4 = read1(Float32),
		user5 = read1(Float32),
		user6 = read1(Float32),
		user7 = read1(Float32),
		user8 = read1(Float32),
		user9 = read1(Float32),
		user10 = read1(Float32),
		user11 = read1(Float32),
		user12 = read1(Float32),
		user13 = read1(Float32),
		user14 = read1(Float32),
		user15 = read1(Float32),
		user16 = read1(Float32),
		user17 = read1(Float32),
		user18 = read1(Float32),
		user19 = read1(Float32),
		v_type = read1(Int32),
		v_coefxa = read1(Float32),
		v_coefxb = read1(Float32),
		v_coefxc = read1(Float32),
		v_coefxd = read1(Float32),
		v_coefya = read1(Float32),
		v_coefyb = read1(Float32),
		v_coefyc = read1(Float32),
		v_coefyd = read1(Float32),
		v_coefza = read1(Float32),
		v_coefzb = read1(Float32),
		v_coefzc = read1(Float32),
		v_coefzd = read1(Float32),
		vm_coef1 = read1(Float32),
		vm_coef2 = read1(Float32),
		vm_coef3 = read1(Float32),
		vm_coef4 = read1(Float32),
		v_venc = read1(Float32),
		spectral_width = read1(Float32),
		csi_dims = read1(Int16),
		xcsi = read1(Int16),
		ycsi = read1(Int16),
		zcsi = read1(Int16),
		roilenx = read1(Float32),
		roileny = read1(Float32),
		roilenz = read1(Float32),
		roilocx = read1(Float32),
		roilocy = read1(Float32),
		roilocz = read1(Float32),
		numdwell = read1(Float32),
		ps_command = read1(Int32),
		ps_mps_r1 = read1(Int32),
		ps_mps_r2 = read1(Int32),
		ps_mps_tg = read1(Int32),
		ps_mps_freq = read1(UInt32),
		ps_aps_r1 = read1(Int32),
		ps_aps_r2 = read1(Int32),
		ps_aps_tg = read1(Int32),
		ps_aps_freq = read1(UInt32),
		ps_scalei = read1(Float32),
		ps_scaleq = read1(Float32),
		ps_snr_warning = read1(Int32),
		ps_aps_or_mps = read1(Int32),
		ps_mps_bitmap = read1(Int32),
		ps_powerspec = freadc(256),
		ps_filler1 = read1(Int32),
		ps_filler2 = read1(Int32),
		halfecho = read1(Int16),
		im_size_y = read1(Int16),
		data_collect_type1 = read1(Int32),
		freq_scale = read1(Float32),
		phase_scale = read1(Float32),
		ovl = read1(Int16),
		pclin = read1(Int16),
		pclinnpts = read1(Int16),
		pclinorder = read1(Int16),
		pclinavg = read1(Int16),
		pccon = read1(Int16),
		pcconnpts = read1(Int16),
		pcconorder = read1(Int16),
		pcextcorr = read1(Int16),
		pcgraph = read1(Int16),
		pcileave = read1(Int16),
		hdbestky = read1(Int16),
		pcctrl = read1(Int16),
		pcthrespts = read1(Int16),
		pcdiscbeg = read1(Int16),
		pcdiscmid = read1(Int16),
		pcdiscend = read1(Int16),
		pcthrespct = read1(Int16),
		pcspacial = read1(Int16),
		pctemporal = read1(Int16),
		pcspare = read1(Int16),
		ileaves = read1(Int16),
		kydir = read1(Int16),
		alt = read1(Int16),
		reps = read1(Int16),
		ref = read1(Int16),
		pcconnorm = read1(Float32),
		pcconfitwt = read1(Float32),
		pclinnorm = read1(Float32),
		pclinfitwt = read1(Float32),
		pcbestky = read1(Float32),
		vrgf = read1(Int32),
		vrgfxres = read1(Int32),
		bp_corr = read1(Int32),
		recv_freq_s = read1(Float32),
		recv_freq_e = read1(Float32),
		hniter = read1(Int32),
		fast_rec = read1(Int32),
		refframes = read1(Int32),
		refframep = read1(Int32),
		scnframe = read1(Int32),
		pasframe = read1(Int32),
		user_usage_tag = read1(UInt32),
		user_fill_mapMSW = read1(UInt32),
		user_fill_mapLSW = read1(UInt32),
		user20 = read1(Float32),
		user21 = read1(Float32),
		user22 = read1(Float32),
		user23 = read1(Float32),
		user24 = read1(Float32),
		user25 = read1(Float32),
		user26 = read1(Float32),
		user27 = read1(Float32),
		user28 = read1(Float32),
		user29 = read1(Float32),
		user30 = read1(Float32),
		user31 = read1(Float32),
		user32 = read1(Float32),
		user33 = read1(Float32),
		user34 = read1(Float32),
		user35 = read1(Float32),
		user36 = read1(Float32),
		user37 = read1(Float32),
		user38 = read1(Float32),
		user39 = read1(Float32),
		user40 = read1(Float32),
		user41 = read1(Float32),
		user42 = read1(Float32),
		user43 = read1(Float32),
		user44 = read1(Float32),
		user45 = read1(Float32),
		user46 = read1(Float32),
		user47 = read1(Float32),
		user48 = read1(Float32),
		pcfitorig = read1(Int16),
		pcshotfirst = read1(Int16),
		pcshotlast = read1(Int16),
		pcmultegrp = read1(Int16),
		pclinfix = read1(Int16),
		pcconfix = read1(Int16),
		pclinslope = read1(Float32),
		pcconslope = read1(Float32),
		pccoil = read1(Int16),
		vvsmode = read1(Int16),
		vvsaimgs = read1(Int16),
		vvstr = read1(Int16),
		vvsgender = read1(Int16),
		zip_factor = read1(Int16),
		maxcoef1a = read1(Float32),
		maxcoef1b = read1(Float32),
		maxcoef1c = read1(Float32),
		maxcoef1d = read1(Float32),
		maxcoef2a = read1(Float32),
		maxcoef2b = read1(Float32),
		maxcoef2c = read1(Float32),
		maxcoef2d = read1(Float32),
		maxcoef3a = read1(Float32),
		maxcoef3b = read1(Float32),
		maxcoef3c = read1(Float32),
		maxcoef3d = read1(Float32),
		ut_ctrl = read1(Int32),
		dp_type = read1(Int16),
		arw = read1(Int16),
		vps = read1(Int16),
		mcReconEnable = read1(Int16),
		fov = read1(Float32),
		te = read1(Int32),
		te2 = read1(Int32),
		dfmrbw = read1(Float32),
		dfmctrl = read1(Int32),
		raw_nex = read1(Int32),
		navs_per_pass = read1(Int32),
		dfmxres = read1(Int32),
		dfmptsize = read1(Int32),
		navs_per_view = read1(Int32),
		dfmdebug = read1(Int32),
		dfmthreshold = read1(Float32),
		grid_control = read1(Int16),
		b0map = read1(Int16),
		grid_tediff = read1(Int16),
		grid_motion_comp = read1(Int16),
		grid_radius_a = read1(Float32),
		grid_radius_b = read1(Float32),
		grid_max_gradient = read1(Float32),
		grid_max_slew = read1(Float32),
		grid_scan_fov = read1(Float32),
		grid_a2d_time = read1(Float32),
		grid_density_factor = read1(Float32),
		grid_display_fov = read1(Float32),
		fatwater = read1(Int16),
		fiestamlf = read1(Int16),
		app = read1(Int16),
		rhncoilsel = read1(Int16),
		rhncoillimit = read1(Int16),
		app_option = read1(Int16),
		grad_mode = read1(Int16),
		pfile_passes = read1(Int16),
		asset = read1(Int32),
		asset_calthresh = read1(Int32),
		asset_R = read1(Float32),
		coilConfigUID = read1(UInt32),
		asset_phases = read1(Int32),
		scancent = read1(Float32),
		position = read1(Int32),
		entry = read1(Int32),
		lmhor = read1(Float32),
		last_slice_num = read1(Int32),
		asset_slice_R = read1(Float32),
		asset_slabwrap = read1(Float32),
		dwnav_coeff = read1(Float32),
		dwnav_cor = read1(Int16),
		dwnav_view = read1(Int16),
		dwnav_corecho = read1(Int16),
		dwnav_sview = read1(Int16),
		dwnav_eview = read1(Int16),
		dwnav_sshot = read1(Int16),
		dwnav_eshot = read1(Int16),
		a3dwin_type = read1(Int16),
		a3dwin_apod = read1(Float32),
		a3dwin_q = read1(Float32),
		ime_scic_enable = read1(Int16),
		clariview_type = read1(Int16),
		ime_scic_edge = read1(Float32),
		ime_scic_smooth = read1(Float32),
		ime_scic_focus = read1(Float32),
		clariview_edge = read1(Float32),
		clariview_smooth = read1(Float32),
		clariview_focus = read1(Float32),
		scic_reduction = read1(Float32),
		scic_gauss = read1(Float32),
		scic_threshold = read1(Float32),
		ectricks_no_regions = read1(Int32),
		ectricks_input_regions = read1(Int32),
		psc_reuse = read1(Int16),
		left_blank = read1(Int16),
		right_blank = read1(Int16),
		acquire_type = read1(Int16),
		retro_control = read1(Int16),
		etl = read1(Int16),
		pcref_start = read1(Int16),
		pcref_stop = read1(Int16),
		ref_skip = read1(Int16),
		extra_frames_top = read1(Int16),
		extra_frames_bot = read1(Int16),
		multiphase_type = read1(Int16),
		nphases = read1(Int16),
		pure = read1(Int16),
		pure_scale = read1(Float32),
		new_wnd_level_flag = read1(Int32),
		wnd_image_hist_area = read1(Int32),
		wnd_high_hist = read1(Float32),
		wnd_lower_hist = read1(Float32),
		pure_filter = read1(Int16),
		cfg_pure_filter = read1(Int16),
		cfg_pure_fit_order = read1(Int16),
		cfg_pure_kernelsize_z = read1(Int16),
		cfg_pure_kernelsize_xy = read1(Int16),
		cfg_pure_weight_radius = read1(Int16),
		cfg_pure_intensity_scale = read1(Int16),
		cfg_pure_noise_threshold = read1(Int16),
		wienera = read1(Float32),
		wienerb = read1(Float32),
		wienert2 = read1(Float32),
		wieneresp = read1(Float32),
		wiener = read1(Int16),
		flipfilter = read1(Int16),
		dbgrecon = read1(Int16),
		ech2skip = read1(Int16),
		tricks_type = read1(Int32),
		lcfiesta_phase = read1(Float32),
		lcfiesta = read1(Int16),
		herawflt = read1(Int16),
		herawflt_befnwin = read1(Int16),
		herawflt_befntran = read1(Int16),
		herawflt_befamp = read1(Float32),
		herawflt_hpfamp = read1(Float32),
		heover = read1(Int16),
		pure_correction_threshold = read1(Int16),
		swiftenable = read1(Int32),
		numslabs = read1(Int16),
		numCoilConfigs = read1(UInt16),
		ps_autoshim_status = read1(Int32),
		dynaplan_numphases = read1(Int32),
		medal_cfg = read1(Int16),
		medal_nstack = read1(Int16),
		medal_echo_order = read1(Int16),
		medal_kernel_up = read1(Int16),
		medal_kernel_down = read1(Int16),
		medal_kernel_smooth = read1(Int16),
		medal_start = read1(Int16),
		medal_end = read1(Int16),
		rcideal = read1(UInt32),
		rcdixproc = read1(UInt32),
		df = read1(Float32),
		bw = read1(Float32),
		feextra = read1(Int32),
		raw_pass_size = read1(UInt64),
		sspsave = read1(UInt64),
		udasave = read1(UInt64),
		vibrant = read1(Int16),
		asset_torso = read1(Int16),
		asset_alt_cal = read1(Int32),
		kacq_uid = read1(Int32),
		psc_ta = read1(Int32),
		disk_acq_ctrl = read1(Int32),
		asset_localTx = read1(Int32),
		a3dscale = read1(Float32),
		broad_band_select = read1(Int32),
		scanner_mode = read1(Int16),
		numbvals = read1(Int16),
		numdifdirs = read1(Int16),
		difnext2 = read1(Int16),
		difnextab = fread(100, Int16),
		channel_combine_method = read1(Int16),
		nexForUnacquiredEncodes = read1(Int16),
		a2dscale = read1(Float32),
		dd_mode = read1(Int16),
		dd_q_ta_offset = read1(Int16),
		dd_q_phase_delay = read1(Float32),
		dacq_ctrl_chksum = read1(UInt32),
		patient_checksum = read1(UInt32),
		rcidealctrl = read1(UInt32),
		echotimes = fread(64, Float32),
		asl_perf_weighted_scale = read1(Int16),
		echo_pc_extra_frames_bot = read1(UInt16),
		echo_pc_ctrl = read1(UInt32),
		echo_pc_ylines = read1(UInt16),
		echo_pc_primary_yfirst = read1(UInt16),
		echo_pc_reverse_yfirst = read1(UInt16),
		echo_pc_zlines = read1(UInt16),
		echo_pc_yxfitorder = read1(UInt16),
		channel_combine_filter_type = read1(Int16),
		mavric_control = read1(UInt32),
		mavric_ImageType = read1(UInt32),
		mavric_bin_separation = read1(Int32),
		mavric_b0_offset = fread(40, Float32),
		channel_combine_filter_width = read1(Float32),
		channel_combine_filter_beta = read1(Float32),
		low_pass_nex_filter_width = read1(Float32),
		aps_tg_status = read1(UInt32),
		cal_pass_set_vector = read1(Int32),
		cal_nex_vector = read1(Int32),
		cal_weight_vector = read1(Int32),
		pure_filtering_mode = read1(Int32),
		pure_lambda = read1(Float32),
		pure_tuning_factor_surface = read1(Float32),
		pure_tuning_factor_body = read1(Float32),
		pure_derived_cal_fraction = read1(Float32),
		pure_derived_cal_reapodization = read1(Float32),
		noncart_grid_factor = read1(Float32),
		noncart_dual_traj = read1(Int32),
		noncart_traj_kmax_ratio = read1(Int32),
		noncart_traj_merge_start = read1(Int32),
		noncart_traj_merge_end = read1(Int32),
		oversamplingfactor = read1(Float32),
		nspokes_highres = read1(Int32),
		nspokes_lowres = read1(Int32),
		nrefslices = read1(Int32),
		psmde_pixel_offset = read1(Int32),
		hoecc = read1(Int32),
		hoec_fit_order = read1(Int32),
		esp = read1(Int32),
		viewSharing3D = read1(UInt32),
		gradwarp_interp_type = read1(Int32),
		pure_blur_enable = read1(Int32),
		pure_blur = read1(Float32),
		slice_info_tab_checksum = read1(UInt32),
		mb_factor = read1(UInt32),
		mb_slice_fov_shift_factor = read1(UInt32),
		mb_slice_order_sign = read1(Int32),
		cal_options = read1(UInt32),
		cs_factor = read1(Float32),
		cs_flag = read1(Int32),
		cs_maxiter = read1(Int32),
		cs_consistency = read1(Float32),
		cs_ph_stride = read1(Int32),
		cs_sl_stride = read1(Int32),
		noncart_traj_mode = read1(Int32),
		pure_alpha = read1(Float32),
		pure_otsu_class = read1(Int32),
		pure_exp_wt = read1(Float32),
		pure_erode_dist = read1(Int32),
		pure_dilate_dist = read1(Int32),
		pure_aniso_blur = read1(Int32),
		pure_aniso_erode_dist = read1(Int32),
		pure_aniso_dilate_dist = read1(Int32),
		calmode = read1(Int32),
		cine_ctrl = read1(UInt16),
		cine_ncycles = read1(UInt16),
		cine_recon_phases = read1(UInt16),
		cine_acq_phases = read1(UInt16),
		cine_virtual_coils = read1(UInt16),
		cine_nrejects = read1(UInt16),
		cine_low_rr = read1(UInt16),
		cine_high_rr = read1(UInt16),
		cine_avg_rr = read1(UInt16),
		cine_bpm = read1(UInt16),
		kat_arc_cal_size_ky = read1(UInt16),
		kat_arc_cal_size_kz = read1(UInt16),
		kat_arc_accel_outer = read1(UInt16),
		kat_arc_accel_center = read1(UInt16),
		cardt1map_ctrl = read1(UInt16),
		cardt1map_ti = read1(UInt16),
		propellerCtrl = read1(Int32),
		daviewsPerBlade = read1(Int32),
		rotationThreshold = read1(Float32),
		shiftThreshold = read1(Float32),
		correlationThreshold = read1(Float32),
		npwfactor = read1(Float32),
		viewshd_num = read1(Int32),
		bvalstab = fread(100, Float32),
		synbvalstab = fread(10, Float32),
		numsynbvals = read1(UInt16),
		position_detection = read1(UInt16),
		medal_param = read1(Float32),
		clariview_filter_x_nr_val = read1(UInt16),
		clariview_filter_x_sh_val = read1(UInt16),
		clariview_filter_x_is_factory = read1(UInt16),
		clariview_filter_x_option = read1(UInt16),
		grad_intensity_thres = read1(Float32),
		anne_mask_dilation_length = read1(Int32),
		anne_intensity_thres = read1(Float32),
		anne_channel_percentage = read1(Float32),
		ec3_iterations = read1(Int32),
		combined_coil_map_thres = read1(Float32),
		coil_map_smooth_factor = read1(Float32),
		coil_map_2_filter_width = read1(Float32),
		ec3_iteration_method = read1(Int32),
		img_ctrl = read1(Int32),
		moco_ctrl = read1(Int32),
		distcorr_ctrl = read1(Int32),
		distcorr_B0_filter_size = read1(Int32),
		distcorr_B0_filter_std_dev = read1(Float32),
		tensor_file_number = read1(Int32),
		muse = read1(UInt16),
		muse_nav = read1(UInt16),
		muse_nav_accel = read1(UInt16),
		muse_nav_hnover = read1(UInt16),
		muse_nav_yres = read1(UInt16),
		optimal_snr_noise_type = read1(UInt16),
		virtual_chan_policy = read1(UInt16),
		virtual_chan_sens_type = read1(UInt16),
		virtual_chan_apply_type = read1(UInt16),
		virtual_chan_uniformity_type = read1(UInt16),
		num_virtual_channels = read1(UInt16),
		excess = fread(507, Int16),
	)
end


text:
module RobotZoo

using StaticArrays
using RobotDynamics
using Rotations
using LinearAlgebra
using ForwardDiff, FiniteDiff
const RD = RobotDynamics

using RobotDynamics: ContinuousDynamics, RigidBody, LieGroupModel, @autodiff
import RobotDynamics: dynamics!, dynamics, jacobian!, forces, moments, wrenches, inertia, inertia_inv, orientation
import RobotDynamics: state_dim, control_dim

include("acrobot.jl")
include("car.jl")
include("cartpole.jl")
include("double_integrator.jl")
include("pendulum.jl")
include("quadrotor.jl")
include("planar_quad.jl")
include("yak_plane.jl")
include("satellite.jl")
include("freebody.jl")
include("rocket.jl")

include("random_linear.jl")
include("linearmodel.jl")

include("LinearModels.jl")

end # module


text:
import Base: +, -, *, /, abs, sqrt, ≈, conj, real, imag, exp, ==

struct ComplexNumber <: Number
    real::Real
    imag::Real
end

const jm = ComplexNumber(0, 1)

abs(x::ComplexNumber) = sqrt(x.real^2 + x.imag^2)

+(x::ComplexNumber, y::ComplexNumber) = ComplexNumber(x.real + y.real, x.imag + y.imag)

+(x::Real, y::ComplexNumber) = ComplexNumber(x + y.real, y.imag)

-(x::ComplexNumber, y::ComplexNumber) = ComplexNumber(x.real - y.real, x.imag - y.imag)

function *(x::ComplexNumber, y::ComplexNumber)
    real = x.real * y.real - x.imag * y.imag
    imag = x.imag * y.real + x.real * y.imag
    return ComplexNumber(real, imag)
end

*(x::Real, y::ComplexNumber) = ComplexNumber(x * y.real, x * y.imag)

function /(x::ComplexNumber, y::ComplexNumber)
    real = (x.real * y.real + x.imag * y.imag) / (y.real^2 + y.imag^2)
    imag = (x.imag * y.real - x.real * y.imag) / (y.real^2 + y.imag^2)
    return ComplexNumber(real, imag)
end

exp(x::ComplexNumber) = exp(x.real) * ComplexNumber(cos(x.imag), sin(x.imag))

==(x::ComplexNumber, y::ComplexNumber) = (x.real == y.real) && (x.imag == y.imag)

≈(x::ComplexNumber, y::ComplexNumber) = isapprox(x.real, y.real, atol=1e-6) && isapprox(x.imag, y.imag, atol=1e-6)

conj(x::ComplexNumber) = ComplexNumber(x.real, -x.imag)

real(x::ComplexNumber) = x.real

imag(x::ComplexNumber) = x.imag

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
using DynamicGrids, Test
import DynamicGrids: neighbors, SimData

init = [0 0 0 1 1 1;
        1 0 1 1 0 1;
        0 1 1 1 1 1;
        0 1 0 0 1 0;
        0 0 0 0 1 1;
        0 1 0 1 1 0]

moore = RadialNeighborhood{1}()
vonneumann = VonNeumannNeighborhood()
t = 1

buf = [0 0 0
       0 1 0
       0 0 0]
state = buf[2, 2]
@test neighbors(moore, nothing, buf, state) == 0
@test neighbors(vonneumann, nothing, buf, state) == 0

buf = [1 1 1
       1 0 1
       1 1 1]
state = buf[2, 2]
@test neighbors(moore, nothing, buf, state) == 8
@test neighbors(vonneumann, nothing, buf, state) == 4

buf = [1 1 1
       0 0 1
       0 0 1]
state = buf[2, 2]
@test neighbors(moore, nothing, buf, state) == 5
@test neighbors(vonneumann, nothing, buf, state) == 2

buf = [0 1 0 0 1
       0 0 1 0 0
       0 0 0 1 1
       0 0 1 0 1
       1 0 1 0 1]
state = buf[3, 3]
custom1 = CustomNeighborhood(((-1,-1), (2,-2), (2,2), (-1,2), (0,0)))
custom2 = CustomNeighborhood(((-1,-1), (0,-1), (1,-1), (2,-1), (0,0)))
layered = LayeredCustomNeighborhood((CustomNeighborhood((-1,1), (-2,2)), CustomNeighborhood((1,2), (2,2))))

@test neighbors(custom1, nothing, buf, state) == 2
@test neighbors(custom2, nothing, buf, state) == 0
@test neighbors(layered, nothing, buf, state) == (1, 2)


text:
using HTTP: HTTP, request, URI
using JSON3: JSON3
using DataFrames: AbstractDataFrame, DataFrames, DataFrame!, DataFrame, rename!, groupby, combine, nrow
using Dates: Dates, Date, format, unix2datetime
using CSV: CSV, File

donativos = DataFrame([String, Float16, String, String, String, String, Date],
                      [:nombre, :cantidad, :afiliación, :comité, :ciudad, :método, :fecha],
                      0)
periods = Date("2019-01-01"):Day(1):(today() - Day(1))
"""
    busqueda_de_donativos!(datos::AbstractDataFrame, día::Date)
"""
function busqueda_de_donativos!(datos::AbstractDataFrame, día::Date)
    día = format(día, "mm/dd/yyyy")
    respuesta = request("POST",
                        URI(scheme = "https",
                            host = "serviciosenlinea.oce.pr.gov",
                            path = "/PublishedDonor/FillPublishedDonorResultsGrid",
                            query = "dateFrom=$día&dateTo=$día&loadResult=true"))
    json = JSON3.read(respuesta.body)
    json.Total ≥ 1_000 && throw(ArgumentError("$día tiene $(json.Total)"))
    isempty(json.Data) && return datos
    data = DataFrame!(json.Data)
    data[!,:Date] .= Date.(unix2datetime.(parse.(Int, SubString.(data[!,:Date], 7, 16))))
    rename!(data, [:nombre, :cantidad, :afiliación, :comité, :ciudad, :método, :fecha])
    append!(datos, data)
end
for period in periods
    sleep(0.1)
    busqueda_de_donativos!(donativos, period)
end
sort!(donativos, :fecha)
CSV.write(joinpath("data", "original", "donativos.tsv"), donativos, delim = '\t')


text:
module Camera
  using CG;
  export ICamera, OrthoCamera, PinholeCamera;

  abstract ICamera;

  # OrthoCamera for ex2:
  type OrthoCamera <: ICamera
    camToWorld::Transformation;
    worldToCam::Transformation;
  end

  OrthoCamera(rc::Vec4f, rv::Vec4f, ru::Vec4f) =
  begin
    a = unitize(xproduct(rv, ru));
    b = unitize(ru); # pointing up
    c = unitize(rv); # pointing backwards
    camToWorld = Transformation(
      Vec4f(a[1], b[1], c[1], rc[1]),
      Vec4f(a[2], b[2], c[2], rc[2]),
      Vec4f(a[3], b[3], c[3], rc[3]),
      Vec4f(0, 0, 0, 1.0)
    );

    worldToCam = inv(camToWorld);
    return OrthoCamera(camToWorld, worldToCam);
  end

  # PinholeCamera for ex3:
  type PinholeCamera <: ICamera
    camToWorld::Transformation;
    worldToCam::Transformation;
    # screen resolution
    nx::Int;
    ny::Int;
    # screen dimensions
    w::Float32;
    h::Float32;
    # distance eye to screen
    d::Float32;

    PinholeCamera(camToWorld::Transformation, worldToCam::Transformation) =
      new(camToWorld, worldToCam, 800, 800, 2.0f0, 2.0f0, 2.0f0);
  end

  function PinholeCamera(rc::Vec4f, rv::Vec4f, ru::Vec4f)
    oc = OrthoCamera(rc, rv, ru);
    return PinholeCamera(oc.camToWorld, oc.worldToCam);
  end
end

max_stars_repo_name: fpoppinga/cg

max_issues_repo_name: fpoppinga/cg

max_forks_repo_name: fpoppinga/cg


text:
using Taste
using FactCheck

facts("Have a beer") do
	todays_selection = [
	    Lambic = Bitter(rand(11:23))
	    Blonde = Bitter(rand(15:30))
	    Kölsch = Bitter(rand(18:25))
	    Märzen = Bitter(rand(18:25))
	    English = Bitter(rand(20:35))
	    Porter = Bitter(rand(20:40))
	    Brown = Bitter(rand(15:25))
	    NorthAmericanBrown = Bitter(rand(25:45))
	    Pilsener = Bitter(rand(30:45))
	    IPA = Bitter(rand(40:140))
	    IrishStout = Bitter(rand(25:60))
	]

	randbeer = rand(todays_selection)
    info(randbeer)

    istasty(beer::Bitter) = beer.IBU < 200
	@fact istasty(randbeer) => true
end

facts("Test Mac Tastebud Sensor") do
    if OS_NAME == :Darwin
        @fact measureumami() >= 0 => true
    else
        @fact_throws ErrorException measureumami()
    end
end


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

mutable struct RedirectTypeEnum <: SwaggerModel

    function RedirectTypeEnum(;)
        o = new()
        o
    end
end # type RedirectTypeEnum

const _property_map_RedirectTypeEnum = Dict{Symbol,Symbol}()
const _property_types_RedirectTypeEnum = Dict{Symbol,String}()
Base.propertynames(::Type{ RedirectTypeEnum }) = collect(keys(_property_map_RedirectTypeEnum))
Swagger.property_type(::Type{ RedirectTypeEnum }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_RedirectTypeEnum[name]))}
Swagger.field_name(::Type{ RedirectTypeEnum }, property_name::Symbol) =  _property_map_RedirectTypeEnum[property_name]

function check_required(o::RedirectTypeEnum)
    true
end

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