text:
## Discretize given percentiles (unweighted)
function discretize_given_percentiles(inV :: Vector{Float64}, pctV :: Vector{Float64}, dbg :: Bool = false)
    edgeV = bin_edges_from_percentiles(inV, pctV, dbg);
    return discretize(inV, edgeV, dbg);
end

## Discretize given percentiles (weighted)
# function discretize_given_percentiles(inV :: Vector{Float64}, wtV :: AbstractWeights,
#    pctV :: Vector{Float64}, dbg :: Bool = false)

#    edgeV = bin_edges_from_percentiles(inV, wtV, pctV, dbg);
#    return discretize(inV, edgeV, dbg);
# end

## Discretize given bounds
# Equivalent to Discretizers LinearDiscretizer
function discretize(inV :: Vector{Float64}, edgeV :: Vector{Float64}, dbg :: Bool = false)
    n = length(edgeV);
    outV = zeros(Int64, size(inV));
    for i1 = n : -1 : 2
        outV[inV .<= edgeV[i1]] .= i1 - 1;
    end
    outV[inV .<= edgeV[1]] .= 0;
    return outV
end

"""
Bin edges from percentiles (unweighted)
Lowest bin includes minimum of inV
"""
function bin_edges_from_percentiles(inV :: Vector{Float64}, pctV :: Vector{Float64}, dbg :: Bool = false)
    edgeV = quantile(inV, [0.0; pctV]);
    # Ensure that lowest point is inside edges
    edgeV[1] -= 1e-8;
    if dbg
        @assert (minimum(inV) >= edgeV[1])
    end
    return edgeV
end

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


text:
abstract type FastText{LANG} <: EmbeddingSystem{LANG} end
abstract type FastText_Text{LANG} <: FastText{LANG} end
abstract type FastText_Bin{LANG} <: FastText{LANG} end

function _load_embeddings(::Type{<:FastText_Bin}, embedding_file, max_vocab_size, keep_words)
    error("FastText Binary Format not supported. If anyone knows how to parse it please feel encouraged to make a PR.")
end

function _load_embeddings(::Type{<:FastText_Text}, embedding_file, max_vocab_size, keep_words)
   #If there are any words in keep_words, then only those are kept, otherwise all are kept
    local LL, indexed_words, index

    if length(keep_words) > 0
        max_stored_vocab_size = length(keep_words)
    end

    open(embedding_file,"r") do fh
        vocab_size, vector_size = parse.(Int64, split(readline(fh)))
        max_stored_vocab_size = min(max_vocab_size, vocab_size)

        indexed_words = Vector{String}(undef, max_stored_vocab_size)
        LL = Array{Float32}(undef, vector_size, max_stored_vocab_size)
        index = 1
        @inbounds for _ in 1:vocab_size
            line = readline(fh)
            toks = split(line)
            word = first(toks)
            if length(keep_words)==0 || word in keep_words
                indexed_words[index]=word
                LL[:,index] .= parse.(Float32, @view toks[2:end])

                index+=1
                if index>max_stored_vocab_size
                    break
                end
            end
        end
    end

    LL = LL[:,1:index-1] #throw away unused columns
    indexed_words = indexed_words[1:index-1] #throw away unused columns
    LL, indexed_words
end

function init(::Type{FastText})

    #########
    # English
    for (source, name, hashstring) in [
            ("Common Crawl", "crawl-300d-2M.vec", "5bfffffbabdab299d4c9165c47275e8f982807a6eaca37ee1f71d3a79ddb544d"),
            ("Wiki News", "wiki-news-300d-1M.vec","bdeb85f44892c505953e3654183e9cb0d792ee51be0992460593e27198d746f8")
        ]

        push!(language_files(FastText_Text{:en}), "FastText $(source)/$(name)")
        register(DataDep("FastText $(source)",
            """
            Dataset: FastText Word Embeddings for English (original release)
            Author: Bojanowski et. al. (Facebook)
            License: CC-SA 3.0
            Website: https://fasttext.cc/docs/en/english-vectors.html

            1 million 300 dimentional  word vectors trained on Wikipedia 2017, UMBC webbase corpus and statmt.org news dataset (16B tokens)
            Citation: P. Bojanowski*, E. Grave*, A. Joulin, T. Mikolov, Enriching Word Vectors with Subword Information

            Notice: this file is ~ 1GB
            """,
            "https://s3-us-west-1.amazonaws.com/fasttext-vectors/$(name).zip",
            hashstring,
            post_fetch_method=DataDeps.unpack
        ));
    end

    #########################
    # Common Crawl
    for (lang, text_hashstring, bin_hashstring) in fast_commoncrawl_languages_and_hashes
        push!(language_files(FastText_Bin{lang}), "FastText $lang CommonCrawl Binary/cc.$(lang).300.bin")
        push!(language_files(FastText_Text{lang}), "FastText $lang CommonCrawl Text/cc.$(lang).300.vec")

        for (mode, hashstring, ext) in [("Text", text_hashstring, "vec"), ("Binary", bin_hashstring, "bin")]
            register(DataDep("FastText $lang CommonCrawl $mode",
                """
                Dataset: 300 dimentional FastText Word Embeddings, for $lang trained on Wikipedia and the CommonCrawl
                Website: https://fasttext.cc/docs/en/crawl-vectors.html
                Author:  Grave et. al. (Facebook)
                License: CC-SA 3.0
                Citation: E. Grave*, P. Bojanowski*, P. Gupta, A. Joulin, T. Mikolov, Learning Word Vectors for 157 Languages
                """,
                "https://dl.fbaipublicfiles.com/fasttext/vectors-crawl/cc.$(lang).300.$(ext).gz",
                hashstring;
                post_fetch_method=DataDeps.unpack
            ));
        end
    end

    for (lang, hashstring) in fast_text_wiki_languages_and_hashes
        # TODO Add Binary files as well
        push!(language_files(FastText_Text{lang}), "FastText $lang Wiki Text/wiki.$(lang).vec")
        register(DataDep("FastText $lang Wiki Text",
            """
            Dataset: 300 dimentional FastText Word Embeddings for $lang, trained on Wikipedia
            Website: https://fasttext.cc/docs/en/pretrained-vectors.html
            Author:  Bojanowski et. al. (Facebook)
            License: CC-SA 3.0
            Citation: P. Bojanowski*, E. Grave*, A. Joulin, T. Mikolov, Enriching Word Vectors with Subword Information
            """,
            "https://dl.fbaipublicfiles.com/fasttext/vectors-wiki/wiki.$(lang).vec",
            hashstring;
        ));
    end
end

# Lang, TextHash, BinHash
const fast_commoncrawl_languages_and_hashes = [
	(:af, nothing, nothing),
	(:als, nothing, nothing),
	(:am, nothing, nothing),
	(:an, nothing, nothing),
	(:ar, nothing, nothing),
	(:arz, nothing, nothing),
	(:as, nothing, nothing),
	(:ast, nothing, nothing),
	(:az, nothing, nothing),
	(:azb, nothing, nothing),
	(:ba, nothing, nothing),
	(:bar, nothing, nothing),
	(:bcl, nothing, nothing),
	(:be, nothing, nothing),
	(:bg, nothing, nothing),
	(:bh, nothing, nothing),
	(:bn, nothing, nothing),
	(:bo, nothing, nothing),
	(:bpy, nothing, nothing),
	(:br, nothing, nothing),
	(:bs, nothing, nothing),
	(:ca, nothing, nothing),
	(:ce, nothing, nothing),
	(:ceb, nothing, nothing),
	(:ckb, nothing, nothing),
	(:co, nothing, nothing),
	(:cs, nothing, nothing),
	(:cv, nothing, nothing),
	(:cy, nothing, nothing),
	(:da, nothing, nothing),
	(:de, nothing, nothing),
	(:diq, nothing, nothing),
	(:dv, nothing, nothing),
	(:el, nothing, nothing),
	(:eml, nothing, nothing),
	(:eo, nothing, nothing),
	(:es, nothing, nothing),
	(:et, nothing, nothing),
	(:eu, nothing, nothing),
	(:fa, nothing, nothing),
	(:fi, nothing, nothing),
	(:fr, "ab0dca4ef2a8a38d97ca119b491c701c57e17178cfc2f032f3de973e86fe87aa", nothing),
	(:frr, nothing, nothing),
	(:fy, nothing, nothing),
	(:ga, nothing, nothing),
	(:gd, nothing, nothing),
	(:gl, nothing, nothing),
	(:gom, nothing, nothing),
	(:gu, nothing, nothing),
	(:gv, nothing, nothing),
	(:he, nothing, nothing),
	(:hi, nothing, nothing),
	(:hif, nothing, nothing),
	(:hr, nothing, nothing),
	(:hsb, nothing, nothing),
	(:ht, nothing, nothing),
	(:hu, nothing, nothing),
	(:hy, nothing, nothing),
	(:ia, nothing, nothing),
	(:id, nothing, nothing),
	(:ilo, nothing, nothing),
	(:io, nothing, nothing),
	(:is, nothing, nothing),
	(:it, nothing, nothing),
	(:ja, nothing, nothing),
	(:jv, nothing, nothing),
	(:ka, nothing, nothing),
	(:kk, nothing, nothing),
	(:km, nothing, nothing),
	(:kn, nothing, nothing),
	(:ko, nothing, nothing),
	(:ku, nothing, nothing),
	(:ky, nothing, nothing),
	(:la, nothing, nothing),
	(:lb, nothing, nothing),
	(:li, nothing, nothing),
	(:lmo, nothing, nothing),
	(:lt, nothing, nothing),
	(:lv, nothing, nothing),
	(:mai, nothing, nothing),
	(:mg, nothing, nothing),
	(:mhr, nothing, nothing),
	(:min, nothing, nothing),
	(:mk, nothing, nothing),
	(:ml, nothing, nothing),
	(:mn, nothing, nothing),
	(:mr, nothing, nothing),
	(:mrj, nothing, nothing),
	(:ms, nothing, nothing),
	(:mt, nothing, nothing),
	(:mwl, nothing, nothing),
	(:my, nothing, nothing),
	(:myv, nothing, nothing),
	(:mzn, nothing, nothing),
	(:nah, nothing, nothing),
	(:nap, nothing, nothing),
	(:nds, nothing, nothing),
	(:ne, nothing, nothing),
	(:new, nothing, nothing),
	(:nl, nothing, nothing),
	(:nn, nothing, nothing),
	(:no, nothing, nothing),
	(:nso, nothing, nothing),
	(:oc, nothing, nothing),
	(:or, nothing, nothing),
	(:os, nothing, nothing),
	(:pa, nothing, nothing),
	(:pam, nothing, nothing),
	(:pfl, nothing, nothing),
	(:pl, nothing, nothing),
	(:pms, nothing, nothing),
	(:pnb, nothing, nothing),
	(:ps, nothing, nothing),
	(:pt, nothing, nothing),
	(:qu, nothing, nothing),
	(:rm, nothing, nothing),
	(:ro, nothing, nothing),
	(:ru, nothing, nothing),
	(:sa, nothing, nothing),
	(:sah, nothing, nothing),
	(:sc, nothing, nothing),
	(:scn, nothing, nothing),
	(:sco, nothing, nothing),
	(:sd, nothing, nothing),
	(:sh, nothing, nothing),
	(:si, nothing, nothing),
	(:sk, nothing, nothing),
	(:sl, nothing, nothing),
	(:so, nothing, nothing),
	(:sq, nothing, nothing),
	(:sr, nothing, nothing),
	(:su, nothing, nothing),
	(:sv, nothing, nothing),
	(:sw, nothing, nothing),
	(:ta, nothing, nothing),
	(:te, nothing, nothing),
	(:tg, nothing, nothing),
	(:th, nothing, nothing),
	(:tk, nothing, nothing),
	(:tl, nothing, nothing),
	(:tr, nothing, nothing),
	(:tt, nothing, nothing),
	(:ug, nothing, nothing),
	(:uk, nothing, nothing),
	(:ur, nothing, nothing),
	(:uz, nothing, nothing),
	(:vec, nothing, nothing),
	(:vi, nothing, nothing),
	(:vls, nothing, nothing),
	(:vo, nothing, nothing),
	(:wa, nothing, nothing),
	(:war, nothing, nothing),
	(:xmf, nothing, nothing),
	(:yi, nothing, nothing),
	(:yo, nothing, nothing),
	(:zea, nothing, nothing),
	(:zh, nothing, nothing),
]

const fast_text_wiki_languages_and_hashes = [
    (:en, "f4d87723baad28804f89c0ecf74fd0f52ac2ae194c270cb4c89b0a84f0bcf53b"),
    # PR Welcome to add Hashes for the ones below
    (:ab, nothing),
    (:ace, nothing),
    (:ady, nothing),
    (:aa, nothing),
    (:af, nothing),
    (:ak, nothing),
    (:sq, nothing),
    (:als, nothing),
    (:am, nothing),
    (:ang, nothing),
    (:ar, nothing),
    (:an, nothing),
    (:arc, nothing),
    (:hy, nothing),
    (:roa_rup, nothing),
    (:as, nothing),
    (:ast, nothing),
    (:av, nothing),
    (:ay, nothing),
    (:az, nothing),
    (:bm, nothing),
    (:bjn, nothing),
    (:map_bms, nothing),
    (:ba, nothing),
    (:eu, nothing),
    (:bar, nothing),
    (:be, nothing),
    (:bn, nothing),
    (:bh, nothing),
    (:bpy, nothing),
    (:bi, nothing),
    (:bs, nothing),
    (:br, nothing),
    (:bug, nothing),
    (:bg, nothing),
    (:my, nothing),
    (:bxr, nothing),
    (:zh_yue, nothing),
    (:ca, nothing),
    (:ceb, nothing),
    (:bcl, nothing),
    (:ch, nothing),
    (:cbk_zam, nothing),
    (:ce, nothing),
    (:chr, nothing),
    (:chy, nothing),
    (:ny, nothing),
    (:zh, nothing),
    (:cho, nothing),
    (:cv, nothing),
    (:zh_classical, nothing),
    (:kw, nothing),
    (:co, nothing),
    (:cr, nothing),
    (:crh, nothing),
    (:hr, nothing),
    (:cs, nothing),
    (:da, nothing),
    (:dv, nothing),
    (:nl, nothing),
    (:nds_nl, nothing),
    (:dz, nothing),
    (:pa, nothing),
    (:arz, nothing),
    (:eml, nothing),
    (:myv, nothing),
    (:eo, nothing),
    (:et, nothing),
    (:ee, nothing),
    (:ext, nothing),
    (:fo, nothing),
    (:hif, nothing),
    (:fj, nothing),
    (:fi, nothing),
    (:frp, nothing),
    (:fr, nothing),
    (:fur, nothing),
    (:ff, nothing),
    (:gag, nothing),
    (:gl, nothing),
    (:gan, nothing),
    (:ka, nothing),
    (:de, nothing),
    (:glk, nothing),
    (:gom, nothing),
    (:got, nothing),
    (:el, nothing),
    (:kl, nothing),
    (:gn, nothing),
    (:gu, nothing),
    (:ht, nothing),
    (:hak, nothing),
    (:ha, nothing),
    (:haw, nothing),
    (:he, nothing),
    (:hz, nothing),
    (:mrj, nothing),
    (:hi, nothing),
    (:ho, nothing),
    (:hu, nothing),
    (:is, nothing),
    (:io, nothing),
    (:ig, nothing),
    (:ilo, nothing),
    (:id, nothing),
    (:ia, nothing),
    (:ie, nothing),
    (:iu, nothing),
    (:ik, nothing),
    (:ga, nothing),
    (:it, nothing),
    (:jam, nothing),
    (:ja, nothing),
    (:jv, nothing),
    (:kbd, nothing),
    (:kab, nothing),
    (:xal, nothing),
    (:kn, nothing),
    (:kr, nothing),
    (:pam, nothing),
    (:krc, nothing),
    (:kaa, nothing),
    (:ks, nothing),
    (:csb, nothing),
    (:kk, nothing),
    (:km, nothing),
    (:ki, nothing),
    (:rw, nothing),
    (:ky, nothing),
    (:rn, nothing),
    (:kv, nothing),
    (:koi, nothing),
    (:kg, nothing),
    (:ko, nothing),
    (:kj, nothing),
    (:ku, nothing),
    (:ckb, nothing),
    (:lad, nothing),
    (:lbe, nothing),
    (:lo, nothing),
    (:ltg, nothing),
    (:la, nothing),
    (:lv, nothing),
    (:lez, nothing),
    (:lij, nothing),
    (:li, nothing),
    (:ln, nothing),
    (:lt, nothing),
    (:olo, nothing),
    (:jbo, nothing),
    (:lmo, nothing),
    (:nds, nothing),
    (:dsb, nothing),
    (:lg, nothing),
    (:lb, nothing),
    (:mk, nothing),
    (:mai, nothing),
    (:mg, nothing),
    (:ms, nothing),
    (:ml, nothing),
    (:mt, nothing),
    (:gv, nothing),
    (:mi, nothing),
    (:mr, nothing),
    (:mh, nothing),
    (:mzn, nothing),
    (:mhr, nothing),
    (:cdo, nothing),
    (:zh_min_nan, nothing),
    (:min, nothing),
    (:xmf, nothing),
    (:mwl, nothing),
    (:mdf, nothing),
    (:mo, nothing),
    (:mn, nothing),
    (:mus, nothing),
    (:nah, nothing),
    (:na, nothing),
    (:nv, nothing),
    (:ng, nothing),
    (:nap, nothing),
    (:ne, nothing),
    (:new, nothing),
    (:pih, nothing),
    (:nrm, nothing),
    (:frr, nothing),
    (:lrc, nothing),
    (:se, nothing),
    (:nso, nothing),
    (:no, nothing),
    (:nn, nothing),
    (:nov, nothing),
    (:ii, nothing),
    (:oc, nothing),
    (:cu, nothing),
    (:or, nothing),
    (:om, nothing),
    (:os, nothing),
    (:pfl, nothing),
    (:pi, nothing),
    (:pag, nothing),
    (:pap, nothing),
    (:ps, nothing),
    (:pdc, nothing),
    (:fa, nothing),
    (:pcd, nothing),
    (:pms, nothing),
    (:pl, nothing),
    (:pnt, nothing),
    (:pt, nothing),
    (:qu, nothing),
    (:ksh, nothing),
    (:rmy, nothing),
    (:ro, nothing),
    (:rm, nothing),
    (:ru, nothing),
    (:rue, nothing),
    (:sah, nothing),
    (:sm, nothing),
    (:bat_smg, nothing),
    (:sg, nothing),
    (:sa, nothing),
    (:sc, nothing),
    (:stq, nothing),
    (:sco, nothing),
    (:gd, nothing),
    (:sr, nothing),
    (:sh, nothing),
    (:st, nothing),
    (:sn, nothing),
    (:scn, nothing),
    (:szl, nothing),
    (:simple, nothing),
    (:sd, nothing),
    (:si, nothing),
    (:sk, nothing),
    (:sl, nothing),
    (:so, nothing),
    (:azb, nothing),
    (:es, nothing),
    (:srn, nothing),
    (:su, nothing),
    (:sw, nothing),
    (:ss, nothing),
    (:sv, nothing),
    (:tl, nothing),
    (:ty, nothing),
    (:tg, nothing),
    (:ta, nothing),
    (:roa_tara, nothing),
    (:tt, nothing),
    (:te, nothing),
    (:tet, nothing),
    (:th, nothing),
    (:bo, nothing),
    (:ti, nothing),
    (:tpi, nothing),
    (:to, nothing),
    (:ts, nothing),
    (:tn, nothing),
    (:tcy, nothing),
    (:tum, nothing),
    (:tr, nothing),
    (:tk, nothing),
    (:tyv, nothing),
    (:tw, nothing),
    (:udm, nothing),
    (:uk, nothing),
    (:hsb, nothing),
    (:ur, nothing),
    (:ug, nothing),
    (:uz, nothing),
    (:ve, nothing),
    (:vec, nothing),
    (:vep, nothing),
    (:vi, nothing),
    (:vo, nothing),
    (:fiu_vro, nothing),
    (:wa, nothing),
    (:war, nothing),
    (:cy, nothing),
    (:vls, nothing),
    (:fy, nothing),
    (:pnb, nothing),
    (:wo, nothing),
    (:wuu, nothing),
    (:xh, nothing),
    (:yi, nothing),
    (:yo, nothing),
    (:diq, nothing),
    (:zea, nothing),
    (:za, nothing),
    (:zu, nothing),
]

ext: jl


text:
mutable struct CheckpointPopulation
    # Manages all the species
    population::Vector{MLPChromosome}
    popsize::Int64
    species::Vector{Species}
    species_history::Vector{Vector{Tuple{Int, Int}}}#Vector{<:Union{Array{Int64,1}, Tuple{Int, Int}}}
    generation::Int64
    avg_fitness::Vector{Float64}
    best_fitness::Vector{MLPChromosome}
    number_species::Int
    simulation_time::Float64
    #setup::Evaluation
    #evaluate::Function # Evaluates population. Override this method in your experiments
    function CheckpointPopulation(p::Population)
        checkpoint_population =
            new(p.population,
                p.popsize,
                p.species,
                p.species_history,
                p.generation,
                p.avg_fitness,
                p.best_fitness,
                p.number_species,
                p.simulation_time
            )
        return checkpoint_population
    end
end

function generate_population(p::CheckpointPopulation, g::Global)
    Population(g,
               p.population,
               p.popsize,
               p.species,
               p.species_history,
               p.generation,
               p.avg_fitness,
               p.best_fitness,
               p.number_species,
               p.simulation_time
    )
end

function generate_population(p::Population, g::Global)
    p
end

function get_population(filename::String)
    data = jldopen(filename, "r")["simulation"]
    generate_population(data["Population"], data["Global"])
end

function get_simulation_state(filename::String)
    data = jldopen(filename, "r")["simulation"]
    SimulationState(data["Global"], data["Population"])
end

function resume_checkpoint(checkpoint_filename::String)
    # Resumes the simulation from a previous saved point.
    checkpoint = string(checkpoint_filename, ".jld2")
    if isfile(checkpoint)
        #return loadFromFile(checkpoint, "Population")
        return get_population(checkpoint)
    else
        error("There is no checkpoint file '$checkpoint'")
    end
end

#Dates.format(now(), "HH:MM:SS-E:U:yyyy")
using Dates

function create_checkpoint(g::Global, p::Population, report::Bool)
        filename = string(g.cf.checkpoint_filename, ".jld2")
        # Saves global variables.
        addToFile(filename, "Global", g)
        # Saves the current simulation state.
        #population = CheckpointPopulation(p)
        population = p
        addToFile(filename, "Population", population)
        # Saves current time and date
        current_date = Dates.now()
        time = Dates.format(current_date, "HH:MM:SS")
        date = Dates.format(current_date, "dd/mm/yyyy")
        addToFile(filename, "Time", time)
        addToFile(filename, "Date", date)

        if report println("Creating checkpoint file at generation: ", p.generation)  end
end

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
using Documenter, VSB

makedocs(
    format = Documenter.HTML(),
    sitename="VSB.jl",
    pages = ["index.md", "guide.md", "docstrings.md"],
    modules = [VSB]
)

# deploydocs(
#     repo = "github.com/camperD/VSB.jl.git",
# )


text:
using Test
using LoopVectorization
using LinearAlgebra
using BQCESubroutine
#using CheapThreads
using YaoLocations: YaoLocations, plain, Locations, CtrlLocations, merge_locations
using BQCESubroutine: ctrl_offset, bcomspace, bsubspace, broutine!, log2dim
using ThreadingUtilities
using BenchmarkTools

function naive_broutine!(st::AbstractVector, U, locs)
    n = log2dim(st)
    subspace = bsubspace(n, sort(locs))
    comspace = bcomspace(n, sort(locs))
    indices = [idx + 1 for idx in comspace]
    for k in subspace
        idx = indices .+ k
        st[idx] = U * st[idx]
    end
    return st
end

function naive_broutine!(st::AbstractMatrix, U, locs)
    n = log2dim(st)
    subspace = bsubspace(n, sort(locs))
    comspace = bcomspace(n, sort(locs))
    indices = [idx + 1 for idx in comspace]
    for b in 1:size(st, 1)
        for k in subspace
            idx = indices .+ k
            st[b, idx] = U * st[b, idx]
        end
    end
    return st
end

function broutine2x2!(S::AbstractMatrix{Complex{T}}, U::AbstractMatrix, locs::Locations) where T
    plain_locs = plain(locs)::Int
    step_1 = 1 << (plain(locs) - 1)
    step_2 = 1 << plain(locs)

    U11 = U[1, 1]; U12 = U[1, 2];
    U21 = U[2, 1]; U22 = U[2, 2];

    Hmax = size(S, 2)
    nthreads = Threads.nthreads()
    space = BQCESubroutine.CartesianSpace(0:step_1-1, 1:size(S, 1), 0:(Hmax>>>plain_locs)-1)

    if step_1 == 1
        @inbounds @batch for _j in 0:(Hmax>>>plain_locs)-1, b in axes(S, 1)
            j = _j << plain_locs

            ST1 = U11 * S[b, j+1] + U12 * S[b, j+step_1+1]
            ST2 = U21 * S[b, j+1] + U22 * S[b, j+step_1+1]

            S[b, j+1] = ST1
            S[b, j+step_1+1] = ST2
        end
        return S
    end

    if step_1 == 2
        @inbounds @batch for _j in 0:(Hmax>>>plain_locs)-1, b in axes(S, 1)
            j = _j << plain_locs
            Base.Cartesian.@nexprs 2 k->begin
                ST1 = U11 * S[b, j+k] + U12 * S[b, j+k+step_1]
                ST2 = U21 * S[b, j+k] + U22 * S[b, j+k+step_1]

                S[b, j+k] = ST1
                S[b, j+k+step_1] = ST2
            end
        end
        return S
    end
    @inbounds @batch for tid in 1:nthreads
        range = BQCESubroutine.schedule_task(space, tid, nthreads)
        broutine2x2_task!(S, (U11, U12, U21, U22), range, step_1, plain_locs)
    end

    # @inbounds @batch for j in 0:step_2:size(S, 2)-step_1, b in 1:size(S, 1), i in j:j+step_1-1
    #     ST1 = U11 * S[b, i+1] + U12 * S[b, i+step_1+1]
    #     ST2 = U21 * S[b, i+1] + U22 * S[b, i+step_1+1]

    #     S[b, i+1] = ST1
    #     S[b, i+step_1+1] = ST2
    # end
    return S
end

function broutine2x2_task!(S, (U11, U12, U21, U22), range, step_1, plain_locs)
    for (k, b, _j) in range
        j = _j << plain_locs
        i = j+k
        ST1 = U11 * S[b, i+1] + U12 * S[b, i+step_1+1]
        ST2 = U21 * S[b, i+1] + U22 * S[b, i+step_1+1]

        S[b, i+1] = ST1
        S[b, i+step_1+1] = ST2
    end
end

N = 20
S = rand(ComplexF64, 1, 1 << N)
U = rand(ComplexF64, 2, 2)
locs = Locations(19)
S1 = broutine2x2!(copy(S), U, locs);
S2 = naive_broutine!(copy(S), U, locs)
S3 = broutine!(copy(S), U, locs)
S1 ≈ S2

@benchmark broutine2x2!(S, U, locs) setup=(S=$(copy(S)))
@benchmark broutine!(S, U, locs) setup=(S=$(copy(S)))

step_1 = 1 << 0
step_2 = 1 << 1

l1 = 0:2:size(S, 2)-step_1
length(l1)
size(S, 2) >>> 1

l2 = 0:(size(S, 2)>>>1)-1

l2[1] << 1
l2[2] << 1
l2[end] << 1

# S1 ≈ S2
# S1 ≈ S3

# 1 << 1

# @benchmark broutine2x2!($(copy(S)), U, locs)
# @benchmark broutine!($(copy(S)), U, locs)


text:
include("helper.jl")

# Locate and/or download Yeppp!
if isfile("$(pkgdir)deps\\downloads\\yeppp-1.0.0.tar.bz2") || (Libdl.find_library(["libyeppp"]) != "")
else
    info("====== Installing Yeppp! into local directory ======")
    trycmd(`cmd /C mkdir $(pkgdir)deps\\downloads`, err="Unable to create Vectorize/deps/downloads")
    trycmd(`cmd /C mkdir $(pkgdir)deps\\src`, err="Unable to create Vectorize/src/")
    trycmd(`wget64 -P downloads http://bitbucket.org/MDukhan/yeppp/downloads/yeppp-1.0.0.tar.bz2`, err="Unable to download Yeppp!")
    trycmd(`7z x $(pkgdir)deps\\downloads\\yeppp-1.0.0.tar.bz2 -o$(pkgdir)deps\\downloads\\`, err="Unable to decompress tar archive")
    trycmd(`7z x $(pkgdir)deps\\downloads\\yeppp-1.0.0.tar -o$(pkgdir)deps\\src\\`, err="Unable to expand tar archive")
    trycmd(`cmd /C move $(pkgdir)deps\\src\\yeppp-1.0.0 $(pkgdir)deps\\src\\yeppp`, err="Unable to rename Yeppp folder")
    info("====== Successfully installed Yeppp! ======")
end


text:
immutable FEMDiffEqPoisson <: AbstractPoissonFEMAlgorithm end

immutable FEMDiffEqHeatEuler <: AbstractHeatFEMAlgorithm end
immutable FEMDiffEqHeatImplicitEuler <: AbstractHeatFEMAlgorithm end
immutable FEMDiffEqHeatCrankNicholson <: AbstractHeatFEMAlgorithm end
immutable FEMDiffEqHeatSemiImplicitEuler <: AbstractHeatFEMAlgorithm end
immutable FEMDiffEqHeatSemiImplicitCrankNicholson <: AbstractHeatFEMAlgorithm end


text:
module RemoteGtkREPL

    using Base64, Sockets, Serialization, Distributed, Pkg, REPL
    using DataFrames #if I don't use DataFrames here autocompletion doesn't work (not sure why)

    import Distributed: remotecall_fetch, remotecall
    import REPL.REPLCompletions.completions

    export eval_command_remotely, isdone, interrupt_task, run_task, eval_symbol

    """ Reference to the socket to GtkREPL and the console index in GtkREPL's
    ConsoleManager """
    struct GtkREPL_Server
        socket::Sockets.TCPSocket
        console_idx::Int
        remote_module::String
    end

    function estalbish_connection(ARGS)
        gtkrepl_port = parse(Int, ARGS[1])
        idx = parse(Int, ARGS[2]) # console/worker id
        remote_module = ARGS[3]   # the module calling us as a String

        estalbish_connection(gtkrepl_port, idx, remote_module)
    end

    function estalbish_connection(gtkrepl_port, idx, remote_module)

        port, server = start_server()
        socket = connect(gtkrepl_port)
        global gtkrepl = GtkREPL_Server(socket, idx, remote_module)

        remotecall_fetch(include_string, socket, Main,
            "$(remote_module).add_remote_console_cb($(idx), $(port))"
        )

        @async begin
            read_stdout, wr = redirect_stdout()
            watch_stdio_task = @async watch_stream(read_stdout, socket, idx, remote_module)
        end
    end

    function __init__()
        global _run_task = @async begin end
        global serialize_lock = ReentrantLock()
    end

    include("eval.jl")
    include("server.jl")
    include("doc.jl")
    include("stdio.jl")

end # module


text:
# This file contains functions for the Profit Efficiency Reverse DDF DEA model
"""
    ProfitReverseDDFDEAModel
An data structure representing a profit DEA model.
"""
struct ProfitReverseDDFDEAModel <: AbstractProfitDEAModel
    n::Int64
    m::Int64
    s::Int64
    measure::Symbol
    monetary::Bool
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    lambda::SparseMatrixCSC{Float64, Int64}
    techeff::Vector
    alloceff::Vector
    normalization::Vector
    Xtarget::Matrix
    Ytarget::Matrix
    Gxrddf::Union{Vector,Matrix}
    Gyrddf::Union{Vector,Matrix}
end

"""
    deaprofitrddf(X, Y, W, P, measure)
Compute profit efficiency using data envelopment analysis Reverse DDF model for
inputs `X`, outputs `Y`, price of inputs `W`, price of outputs `P`, and efficiency `measure`.

# Measure specification:
- `:ERG`: Enhanced Russell Graph Slack Based Measure.
- `:MDDF`: Modified Directional Distance Function.

# Direction specification:
For the Modified Directional Distance Function, the directions `Gx` and `Gy` can be one of the following symbols.
- `:Observed`: use observed values.

Alternatively, a vector or matrix with the desired directions can be supplied.

# Optional Arguments
- `atol=1e-6`: tolerance for DMU to be considered efficient.
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [2; 4; 8; 12; 6; 14; 14; 9.412];

julia> Y = [1; 5; 8; 9; 3; 7; 9; 2.353];

julia> W = [1; 1; 1; 1; 1; 1; 1; 1];

julia> P = [2; 2; 2; 2; 2; 2; 2; 2];

julia> deaprofitrddf(X, Y, W, P, :ERG)
Profit Reverse DDF DEA Model
DMUs = 8; Inputs = 1; Outputs = 1
Returns to Scale = VRS
Associated efficiency measure = ERG
──────────────────────────────────
     Profit  Technical  Allocative
──────────────────────────────────
1  4.0        0.0         4.0
2  0.5        0.0         0.5
3  0.0        0.0         0.0
4  0.166667   0.0         0.166667
5  0.8        0.6         0.2
6  0.571429   0.52381     0.047619
7  0.285714   0.142857    0.142857
8  0.949449   0.8         0.149449
──────────────────────────────────
```
"""
function deaprofitrddf(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector},
    W::Union{Matrix,Vector}, P::Union{Matrix,Vector},
    measure::Symbol;
    Gx::Union{Symbol,Matrix,Vector,Nothing} = nothing, Gy::Union{Symbol,Matrix,Vector,Nothing} = nothing,
    monetary::Bool = false, atol::Float64 = 1e-6,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::ProfitReverseDDFDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    nw, mw = size(W, 1), size(W, 2)
    np, sp = size(P, 1), size(P, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nw != nx
        throw(DimensionMismatch("number of rows in W and X ($nw, $nx) are not equal"));
    end
    if np != ny
        throw(DimensionMismatch("number of rows in P and Y ($np, $ny) are not equal"));
    end
    if mw != m
        throw(DimensionMismatch("number of columns in W and X ($mw, $m) are not equal"));
    end
    if sp != s
        throw(DimensionMismatch("number of columns in P and Y ($sp, $s) are not equal"));
    end

    n = nx

    # Get directions based on efficiency measure
    Gxrddf, Gyrddf = DataEnvelopmentAnalysis.dearddfdirections(X, Y, measure, Gx = Gx, Gy = Gy, orient = :Graph, rts = :VRS, optimizer = optimizer)

    # Get normalization of efficiency measure model
    if measure == :ERG
        measureprofitmodel = deaprofiterg(X, Y, W, P, optimizer = optimizer)
    elseif measure == :MDDF
        measureprofitmodel = deaprofitmddf(X, Y, W, P, Gx = Gx, Gy = Gy, optimizer = optimizer)
    end

    normmeasureprof = normfactor(measureprofitmodel)
    pefficiencymeasure = efficiency(measureprofitmodel)

    # Profit model using calculated directions
    rddfprofitmodel = deaprofit(X, Y, W, P, Gx = Gxrddf, Gy = Gyrddf, optimizer = optimizer)

    pefficiency = efficiency(rddfprofitmodel)
    plambdaeff = peersmatrix(rddfprofitmodel)
    techefficiency = efficiency(rddfprofitmodel, :Technical)
    allocefficiency = efficiency(rddfprofitmodel, :Allocative)
    normalization = normfactor(rddfprofitmodel)
    Xtarget = targets(rddfprofitmodel, :X)
    Ytarget = targets(rddfprofitmodel, :Y)

    # Change normalization factor of efficient DMUs
    for i=1:n
        if abs(techefficiency[i] ) <= atol
            pefficiency[i] = pefficiencymeasure[i]
            normalization[i] = normmeasureprof[i]
        end
    end

    if monetary
        pefficiency = pefficiency .* normalization
        techefficiency = techefficiency .* normalization
    end

    allocefficiency = pefficiency - techefficiency

    return ProfitReverseDDFDEAModel(n, m, s, measure, monetary, names, pefficiency, plambdaeff, techefficiency, allocefficiency, normalization, Xtarget, Ytarget, Gxrddf, Gyrddf)
end

function Base.show(io::IO, x::ProfitReverseDDFDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    dmunames = names(x)

    eff = efficiency(x)
    techeff = efficiency(x, :Technical)
    alloceff = efficiency(x, :Allocative)

    if !compact
        print(io, "Profit Reverse DDF DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Returns to Scale = VRS")
        print(io, "\n")
        print(io, "Associated efficiency measure = ", string(x.measure))
        print(io, "\n")
        show(io, CoefTable(hcat(eff, techeff, alloceff), ["Profit", "Technical", "Allocative"], dmunames))
    end

end


text:
using RecursiveArrayTools

@inline function lazy_glue(A::AbstractArray{IT,N}, code::Tuple) where {IT,N}
    gluecodecheck(A, code)
    # if code == (:,*)
    #     VectorOfArray(A)
    # elseif
    if code == (*,:)
        transpose(VectorOfArray(A))
    # elseif count(isequal(*), code) == 1 && code[end] == (*)
    #     VectorOfArray(A)
    elseif iscodesorted(code)
        flat = VectorOfArray(vec(A))
        finalsize = (size(first(A))..., size(A)...)
        reshape(flat, finalsize)
        # always reshape so that it prints like an array,
        # also makes copy() produce an array
    else
        error("can't glue code = $code with VectorOfArray")
    end
end


text:
# Tests that init_score works for multiclass

# Create a 3-class problem
N1 = 1000
N2 = 1000
N3 = 1000
N = N1 + N2 + N3

numClasses = 3
numFeats = 5

X = rand((N, numFeats))
y = vcat(zeros(N1), ones(N2), 2*ones(N3))

estimator = LightGBM.LGBMMulticlass(num_iterations = 50,
                                    learning_rate = .5,
                                    feature_fraction = 0.5,
                                    bagging_fraction = 1.0,
                                    bagging_freq = 1,
                                    num_leaves = 5,
                                    metric = ["multi_logloss"],
                                    num_class = numClasses,
                                    min_data_in_leaf=1,
                                    min_sum_hessian_in_leaf=1,
                                    early_stopping_round = 1);

# using the incorrect size of init_score should throw
try
    LightGBM.fit(estimator, X, y; init_score = [1.2, 3.4])
    @test false  # LightGBM.fit did not throw with incorrect init_score size
catch
end

# Create init score as if it were already predicting very well the data
init_score = -10000.0 * ones(numClasses, N)
for i in 1:N
    init_score[Int(y[i] + 1), i] *= -1
end
init_score = reshape(init_score, length(init_score))

LightGBM.fit(estimator, X, y; init_score = init_score)

# then it should have early stopped without training any tree
@test LightGBM.LGBM_BoosterGetCurrentIteration(estimator.booster) == 0


text:
using Rsvg: handle_new_from_data
using Cairo: CairoPDFSurface
import Cairo, Rsvg

function _savefig_cairo(
        p::Plot, raw_svg::AbstractString, fn::AbstractString,
        suf::AbstractString
    )
    if suf == "pdf"
        r = Rsvg.handle_new_from_data(raw_svg)
        cs = Cairo.CairoPDFSurface(fn, size(p)...)
        ctx = Cairo.CairoContext(cs)
        Rsvg.handle_render_cairo(ctx, r)
        Cairo.show_page(ctx)
        Cairo.finish(cs)
    elseif suf == "png"
        r = Rsvg.handle_new_from_data(raw_svg)
        cs = Cairo.CairoImageSurface(size(p)..., Cairo.FORMAT_ARGB32)
        ctx = Cairo.CairoContext(cs)
        Rsvg.handle_render_cairo(ctx, r)
        Cairo.write_to_png(cs, fn)
    elseif suf == "eps"
        r = Rsvg.handle_new_from_data(raw_svg)
        cs = Cairo.CairoEPSSurface(fn, size(p)...)
        ctx = Cairo.CairoContext(cs)
        Rsvg.handle_render_cairo(ctx, r)
        Cairo.show_page(ctx)
        Cairo.finish(cs)
    else
        error("Unknown file-suffix $suf")
    end
end


text:
module Spglib

using spglib_jll: libsymspg
using UnPack: @unpack

# All public methods
export get_symmetry,
    get_symmetry!,
    get_symmetry_with_collinear_spin,
    get_symmetry_with_collinear_spin!,
    get_hall_number_from_symmetry,
    get_dataset,
    get_spacegroup_number,
    get_spacegroup_type,
    get_international,
    get_schoenflies,
    standardize_cell,
    find_primitive,
    refine_cell,
    niggli_reduce,
    delaunay_reduce,
    get_multiplicity,
    get_ir_reciprocal_mesh,
    get_stabilized_reciprocal_mesh,
    get_version

include("model.jl")
include("symmetry.jl")
include("standardize.jl")
include("reduce.jl")
include("reciprocal.jl")

end


text:
# Core defines the core WebSocket types, such as a frame and opcodes.

# Description of a WebSocket frame from https://tools.ietf.org/html/rfc6455, chapter 5.2.
#
#      0                   1                   2                   3
#      0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1
#     +-+-+-+-+-------+-+-------------+-------------------------------+
#     |F|R|R|R| opcode|M| Payload len |    Extended payload length    |
#     |I|S|S|S|  (4)  |A|     (7)     |             (16/64)           |
#     |N|V|V|V|       |S|             |   (if payload len==126/127)   |
#     | |1|2|3|       |K|             |                               |
#     +-+-+-+-+-------+-+-------------+ - - - - - - - - - - - - - - - +
#     |     Extended payload length continued, if payload len == 127  |
#     + - - - - - - - - - - - - - - - +-------------------------------+
#     |                               |Masking-key, if MASK set to 1  |
#     +-------------------------------+-------------------------------+
#     | Masking-key (continued)       |          Payload Data         |
#     +-------------------------------- - - - - - - - - - - - - - - - +
#     :                     Payload Data continued ...                :
#     + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +
#     |                     Payload Data continued ...                |
#     +---------------------------------------------------------------+

# TODO: Rethink what we're exporting here. Should we export anything? None of this
#       should be part of the interface to the user.
export Frame,
       Opcode,
       OPCODE_CONTINUATION,
       OPCODE_TEXT,
       OPCODE_BINARY,
       OPCODE_CLOSE,
       OPCODE_PING,
       OPCODE_PONG

import Base.==

# TODO: Documentation.

struct Opcode
    op::UInt8
end

# Requirement
# @5_2-3 Opcodes
const OPCODE_CONTINUATION = Opcode(0)
const OPCODE_TEXT         = Opcode(1)
const OPCODE_BINARY       = Opcode(2)
const OPCODE_CLOSE        = Opcode(8)
const OPCODE_PING         = Opcode(9)
const OPCODE_PONG         = Opcode(10)

==(a::Opcode, b::Opcode) = a.op == b.op

struct Frame
    fin::Bool
    rsv1::Bool
    rsv2::Bool
    rsv3::Bool
    opcode::Opcode
    ismasked::Bool
    len::UInt8 # Is actually 7 bits.
    extended_len::UInt64
    mask::AbstractArray{UInt8, 1}
    payload::AbstractArray{UInt8, 1}
end

Frame(fin::Bool, opcode::Opcode, ismasked::Bool, len::Int,
      extended_len::Int, mask::AbstractArray{UInt8, 1}, payload::AbstractArray{UInt8, 1}) =
  Frame(fin, false, false, false, opcode, ismasked, len, extended_len, mask, payload)

==(a::Frame, b::Frame) = a.fin == b.fin && a.rsv1 == b.rsv1 && a.rsv2 == b.rsv2 &&
  a.opcode == b.opcode && a.ismasked == b.ismasked && a.len == b.len && a.extended_len == b.extended_len &&
  a.mask == b.mask && a.payload == b.payload

function Base.read(s::IO, ::Type{Frame})
  x    = read(s, UInt8)
  fin  = x & 0b1000_0000 != 0
  rsv1 = x & 0b0100_0000 != 0
  rsv2 = x & 0b0010_0000 != 0
  rsv3 = x & 0b0001_0000 != 0
  op   = x & 0b0000_1111

  # Requirement
  # @5_2-5 The mask bit is 1 if set, 0 if not.
  y        = read(s, UInt8)
  ismasked = y & 0b1000_0000 != 0
  len      = y & 0b0111_1111

  extended_len::UInt64 = 0
  if len == 126
    extended_len = ntoh(read(s, UInt16))
  elseif len == 127
    extended_len = ntoh(read(s, UInt64))
  end

  mask = Array{UInt8,1}()
  if ismasked
    mask = read(s, 4)
  end

  payload_length = extended_len != 0 ? extended_len : len
  payload = read(s, payload_length)

  Frame(fin, rsv1, rsv2, rsv3, Opcode(op), ismasked, len, extended_len, mask, payload)
end

function Base.write(s::IO, frame::Frame)
  x1 =
    UInt8(frame.fin)  << 7 |
    UInt8(frame.rsv1) << 6 |
    UInt8(frame.rsv2) << 5 |
    UInt8(frame.rsv3) << 4 |
    frame.opcode.op & 0b0000_1111

  # Requirement
  # @5_2-5 The mask bit is 1 if set, 0 if not.
  x2 = UInt8(frame.ismasked) << 7 |
       frame.len & 0b0111_1111

  write(s, x1)
  write(s, x2)

  if frame.len == 126
    write(s, hton(UInt16(frame.extended_len)))
  elseif frame.len == 127
    write(s, hton(frame.extended_len))
  end

  if frame.ismasked
    write(s, frame.mask)
  end

  write(s, frame.payload)
end


text:
# # Classification of the Wisconsin breast cancer dataset
#
# In this case study we will deal with the Wisconsin breast cancer dataset which can be
# browsed freely on the [UCI website](https://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Wisconsin+(Diagnostic)).
#
# In particular, this dataset contains *10 features* and 699 instances. In the work we
# will do here, however, we will skip some instances due to some missing values.
#
# The dataset contains only two classes, and the purpose is to use all ten features to
# answer a simple question:
# > Does the subject have a benign or malign tumor?
# To answer this question, we will train a Least Squares Support Vector Machine as
# implemented in `LeastSquaresSVM`.
#
# First, we need to import all the necessary packages.

using MLJ, MLJBase
using DataFrames, CSV
using CategoricalArrays
using Random, Statistics
using LeastSquaresSVM

# We then need to specify a seed to enable reproducibility of the results.
rng = MersenneTwister(801239);

# Here we are creating a list with all the headers.
headers = [
	"id", "Clump Thickness",
	"Uniformity of Cell Size", "Uniformity of Cell Shape",
	"Marginal Adhesion", "Single Epithelial Cell Size",
	"Bare Nuclei", "Bland Chromatin",
	"Normal Nucleoli", "Mitoses", "class"
];

# We define the path were the dataset is located
path = joinpath("src", "examples", "wbc.csv");

# We load the csv file and convert it to a `DataFrame`. Note that we are specifying
# to the file reader to replace the string `?` to a `missing` value. This dataset contains
# the the string `?` when there is a value missing.
data = CSV.File(path; header=headers, missingstring="?") |> DataFrame;

# We can display the first 10 rows from the dataset
first(data, 10)

# We can see that all the features have been added correctly, we can see that we have
# an unncessary feature called `id`, so we will remove it.

select!(data, Not(:id));

# We also need to remove all the missing data from the `DataFrame`
data = dropmissing(data);

# The `class` column should be of type `categorical`, following the `MLJ` API, so we
# encode it here.
transform!(data, :class => categorical, renamecols=false);

# Check statistics per column.
describe(data)

# Split the dataset into training and testing.
y, X = unpack(data, ==(:class), colname -> true);

# We will use only 2/3 for training.
train, test = partition(eachindex(y), 2 / 3, shuffle=true, rng=rng);

# Always remove mean and set the standard deviation to 1.0 when dealing with SVMs.
stand1 = Standardizer(count=true);
X = MLJBase.transform(fit!(machine(stand1, X)), X);

# Check statistics per column again to ensure standardization, but remember to do it now
# with the `X` matrix.
describe(X)

# Good, now every column has a mean very close to zero, so the standardization was
# done correctly.

# We now create our model with `LeastSquaresSVM`
model = LeastSquaresSVM.LSSVClassifier();

# These are the values for the hyperparameter grid search. We need to find the best subset
# from this set of parameters.
# Although I will not do this here, the best approach is to find a set of good hyperparameters
# and then refine the search space around that set. That way we can ensure we will always get
# the best results.
sigma_values = [0.5, 5.0, 10.0, 15.0, 25.0, 50.0, 100.0, 250.0, 500.0];
r1 = MLJBase.range(model, :σ, values=sigma_values);
gamma_values = [0.01, 0.05, 0.1, 0.5, 1.0, 5.0, 10.0, 50.0, 100.0, 500.0, 1000.0];
r2 = MLJBase.range(model, :γ, values=gamma_values);

# We now create a `TunedModel` that will use a 10-folds stratified cross validation scheme
# in order to find the best set of hyperparameters. The stratification is needed because
# the classes are somewhat imbalanced:
# - Benign: 458 (65.5%)
# - Malignant: 241 (34.5%)

self_tuning_model = TunedModel(
    model=model,
    tuning=Grid(rng=rng),
    resampling=StratifiedCV(nfolds=10),
    range=[r1, r2],
    measure=accuracy,
    acceleration=CPUThreads(), # We use this to enable multithreading
);

# Once the best model is found, we create a `machine` with it, and fit it
mach = machine(self_tuning_model, X, y);
fit!(mach, rows=train, verbosity=0);

# We can now show the best hyperparameters found.
fitted_params(mach).best_model

# And we test the trained model. We expect somewhere around 94%-96% accuracy.
results = predict(mach, rows=test);
acc = accuracy(results, y[test]);

# Show the accuracy for the testing set
println(acc * 100.0)

# As you can see, it is fairly easy to use `LeastSquaresSVM` together with MLJ. We got a good
# accuracy result and this proves that the implementation is actually correct. This
# dataset is commonly used as a benchmark dataset to test new algorithms.

max_stars_repo_name: gmagannaDevelop/Elysivm

max_issues_repo_name: gmagannaDevelop/Elysivm

max_forks_repo_name: gmagannaDevelop/Elysivm


text:
# xaxis(args...; kw...) = Axis(:x, args...; kw...)
# yaxis(args...; kw...) = Axis(:y, args...; kw...)
# zaxis(args...; kw...) = Axis(:z, args...; kw...)

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

function Axis(sp::Subplot, letter::Symbol, args...; kw...)
    # init with values from _plot_defaults
    d = KW(
        :letter => letter,
        # :extrema => (Inf, -Inf),
        :extrema => Extrema(),
        :discrete_map => Dict(),   # map discrete values to discrete indices
        :continuous_values => zeros(0),
        :use_minor => false,
        :show => true,  # show or hide the axis? (useful for linked subplots)
    )

    # get defaults from letter version, unless match
    for (k,v) in _axis_defaults
        lk = Symbol(letter, k)
        lv = _axis_defaults_byletter[lk]
        d[k] = (lv == :match ? v : lv)
    end

    # merge!(d, _axis_defaults)
    d[:discrete_values] = []

    # update the defaults
    update!(Axis(sp, d), args...; kw...)
end

function get_axis(sp::Subplot, letter::Symbol)
    axissym = Symbol(letter, :axis)
    if haskey(sp.attr, axissym)
        sp.attr[axissym]
    else
        sp.attr[axissym] = Axis(sp, letter)
    end
end

function process_axis_arg!(d::KW, arg, letter = "")
    T = typeof(arg)
    arg = get(_scaleAliases, arg, arg)

    if typeof(arg) <: Font
        d[Symbol(letter,:tickfont)] = arg
        d[Symbol(letter,:guidefont)] = arg

    elseif arg in _allScales
        d[Symbol(letter,:scale)] = arg

    elseif arg in (:flip, :invert, :inverted)
        d[Symbol(letter,:flip)] = true

    elseif T <: AbstractString
        d[Symbol(letter,:guide)] = arg

    # xlims/ylims
    elseif (T <: Tuple || T <: AVec) && length(arg) == 2
        sym = typeof(arg[1]) <: Number ? :lims : :ticks
        d[Symbol(letter,sym)] = arg

    # xticks/yticks
    elseif T <: AVec
        d[Symbol(letter,:ticks)] = arg

    elseif arg == nothing
        d[Symbol(letter,:ticks)] = []

    elseif typeof(arg) <: Number
        d[Symbol(letter,:rotation)] = arg

    else
        warn("Skipped $(letter)axis arg $arg")

    end
end

# update an Axis object with magic args and keywords
function update!(axis::Axis, args...; kw...)
    # first process args
    d = axis.d
    for arg in args
        process_axis_arg!(d, arg)
    end

    # then override for any keywords... only those keywords that already exists in d
    for (k,v) in kw
        if haskey(d, k)
            if k == :discrete_values
                # add these discrete values to the axis
                for vi in v
                    discrete_value!(axis, vi)
                end
            else
                d[k] = v
            end
        end
    end

    # replace scale aliases
    if haskey(_scaleAliases, d[:scale])
        d[:scale] = _scaleAliases[d[:scale]]
    end

    axis
end

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

Base.show(io::IO, axis::Axis) = dumpdict(axis.d, "Axis", true)
# Base.getindex(axis::Axis, k::Symbol) = getindex(axis.d, k)
Base.setindex!(axis::Axis, v, ks::Symbol...) = setindex!(axis.d, v, ks...)
Base.haskey(axis::Axis, k::Symbol) = haskey(axis.d, k)
Base.extrema(axis::Axis) = (ex = axis[:extrema]; (ex.emin, ex.emax))

# get discrete ticks, or not
function get_ticks(axis::Axis)
    ticks = axis[:ticks]
    dvals = axis[:discrete_values]
    if !isempty(dvals) && ticks == :auto
        axis[:continuous_values], dvals
    else
        ticks
    end
end

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

function expand_extrema!(ex::Extrema, v::Number)
    ex.emin = min(v, ex.emin)
    ex.emax = max(v, ex.emax)
    ex
end

function expand_extrema!(axis::Axis, v::Number)
    expand_extrema!(axis[:extrema], v)
end
function expand_extrema!{MIN<:Number,MAX<:Number}(axis::Axis, v::Tuple{MIN,MAX})
    ex = axis[:extrema]
    ex.emin = min(v[1], ex.emin)
    ex.emax = max(v[2], ex.emax)
    ex
end
function expand_extrema!{N<:Number}(axis::Axis, v::AVec{N})
    ex = axis[:extrema]
    for vi in v
        expand_extrema!(ex, vi)
    end
    ex
end

function expand_extrema!(sp::Subplot, d::KW)
    # first expand for the data
    for letter in (:x, :y, :z)
        data = d[letter]
        axis = sp.attr[Symbol(letter, "axis")]
        if eltype(data) <: Number
            expand_extrema!(axis, data)
        elseif isa(data, Surface) && eltype(data.surf) <: Number
            expand_extrema!(axis, data)
        elseif data != nothing
            # TODO: need more here... gotta track the discrete reference value
            #       as well as any coord offset (think of boxplot shape coords... they all
            #       correspond to the same x-value)
            # @show letter,eltype(data),typeof(data)
            d[letter], d[Symbol(letter,"_discrete_indices")] = discrete_value!(axis, data)
        end
    end

    # # expand for fillrange/bar_width
    # fillaxis, baraxis = sp.attr[:yaxis], sp.attr[:xaxis]
    # if isvertical(d)
    #     fillaxis, baraxis = baraxis, fillaxis
    # end

    # expand for fillrange
    vert = isvertical(d)
    fr = d[:fillrange]
    if fr == nothing && d[:seriestype] == :bar
        fr = 0.0
    end
    if fr != nothing
        axis = sp.attr[vert ? :yaxis : :xaxis]
        if typeof(fr) <: Tuple
            for fri in fr
                expand_extrema!(axis, fri)
            end
        else
            expand_extrema!(axis, fr)
        end
    end

    # expand for bar_width
    if d[:seriestype] == :bar
        dsym = vert ? :x : :y
        data = d[dsym]

        bw = d[:bar_width]
        if bw == nothing
            bw = d[:bar_width] = mean(diff(data))
        end
        # @show data bw

        axis = sp.attr[Symbol(dsym, :axis)]
        expand_extrema!(axis, maximum(data) + 0.5maximum(bw))
        expand_extrema!(axis, minimum(data) - 0.5minimum(bw))
    end

end

function expand_extrema!(sp::Subplot, xmin, xmax, ymin, ymax)
    expand_extrema!(sp[:xaxis], (xmin, xmax))
    expand_extrema!(sp[:yaxis], (ymin, ymax))
end

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

# push the limits out slightly
function widen(lmin, lmax)
    span = lmax - lmin
    # eps = max(1e-16, min(1e-2span, 1e-10))
    eps = max(1e-16, 0.03span)
    lmin-eps, lmax+eps
end

# figure out if widening is a good idea.  if there's a scale set it's too tricky,
# so lazy out and don't widen
function default_should_widen(axis::Axis)
    should_widen = false
    if axis[:scale] == :identity
        for series in series_list(axis.sp)
            if series.d[:seriestype] in (:scatter,) || series.d[:markershape] != :none
                should_widen = true
            end
        end
    end
    should_widen
end

# using the axis extrema and limit overrides, return the min/max value for this axis
function axis_limits(axis::Axis, should_widen::Bool = default_should_widen(axis))
    ex = axis[:extrema]
    amin, amax = ex.emin, ex.emax
    lims = axis[:lims]
    if (isa(lims, Tuple) || isa(lims, AVec)) && length(lims) == 2
        if isfinite(lims[1])
            amin = lims[1]
        end
        if isfinite(lims[2])
            amax = lims[2]
        end
    end
    if amax <= amin && isfinite(amin)
        amax = amin + 1.0
    end
    if should_widen
        widen(amin, amax)
    else
        amin, amax
    end
end

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

# these methods track the discrete (categorical) values which correspond to axis continuous values (cv)
# whenever we have discrete values, we automatically set the ticks to match.
# we return (continuous_value, discrete_index)
function discrete_value!(axis::Axis, dv)
    cv_idx = get(axis[:discrete_map], dv, -1)
    # @show axis[:discrete_map], axis[:discrete_values], dv
    if cv_idx == -1
        ex = axis[:extrema]
        cv = max(0.5, ex.emax + 1.0)
        expand_extrema!(axis, cv)
        push!(axis[:discrete_values], dv)
        push!(axis[:continuous_values], cv)
        cv_idx = length(axis[:discrete_values])
        axis[:discrete_map][dv] = cv_idx
        cv, cv_idx
    else
        cv = axis[:continuous_values][cv_idx]
        cv, cv_idx
    end
end

# continuous value... just pass back with axis negative index
function discrete_value!(axis::Axis, cv::Number)
    cv, -1
end

# add the discrete value for each item.  return the continuous values and the indices
function discrete_value!(axis::Axis, v::AVec)
    n = length(v)
    cvec = zeros(n)
    discrete_indices = zeros(Int, n)
    for i=1:n
        cvec[i], discrete_indices[i] = discrete_value!(axis, v[i])
    end
    cvec, discrete_indices
end

# add the discrete value for each item.  return the continuous values and the indices
function discrete_value!(axis::Axis, v::AMat)
    n,m = size(v)
    cmat = zeros(n,m)
    discrete_indices = zeros(Int, n, m)
    for i=1:n, j=1:m
        cmat[i,j], discrete_indices[i,j] = discrete_value!(axis, v[i,j])
    end
    cmat, discrete_indices
end

function discrete_value!(axis::Axis, v::Surface)
    map(Surface, discrete_value!(axis, v.surf))
end

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

function pie_labels(sp::Subplot, series::Series)
    d = series.d
    if haskey(d,:x_discrete_indices)
        dvals = sp.attr[:xaxis].d[:discrete_values]
        [dvals[idx] for idx in d[:x_discrete_indices]]
    else
        d[:x]
    end
end


text:
## ENV["toolkit"] = "Tk"
using JGUI
using Base.Test

tstore = treestore(Int, Float64, String)

t1 = (1, 1.0, "one")
t11 = (11, 11.0, "one-one")
t12 = (12, 12.0, "one-two")
t2  = (2, 2.0, "two")

node = insert!(tstore, nothing, 1, "label1", t1)
insert!(tstore, node, 1, "label11", t11)
insert!(tstore, node, 2, "label12", t12)
node = insert!(tstore, nothing, 2, "label2", t2)

w = window()
view = treeview(w, tstore)      # have to add tpl if tstore is empty
view[:names] = ["Int", "Float", "String"]
view[:widths] =  [50,50,-1]

tv[:keyname] = "key"
tv[:keywidth] = 100

push!(view); raise(w)

## value is path of selected item
view[:value] = [1,2]
@test view[:value] == [1,2]

node = path_to_node(tstore, [1])
expand_node(view, node)

## remove node
node = path_to_node(tstore, [1,1])
pop!(tstore, node)

##


text:
using ModelVisualisations, InvPendulum, RLExploration

dt = 1/25
t = 0:dt:5
X = zeros(4, length(t))
U = [NaN for i = 1:length(t)]
temp = GaussianProcess()

X[:, 1] .= [0.0, 0.0, 0.0, -pi]
for i = 1:length(t)-1
    U[i] = sample(temp)
    X[:, i+1] .= invpendulum_rk4(X[:, i], U[i], dt)
end

options = Dict(:fps => 25, :size => (1920, 1080), :U => U)
visualise(:InvPendulum, X, t, "invpendulum.mp4"; options = options)

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
int main() {
  int[] a;
  a = new int[50];
  a[0] = true;
  return 0;
}

max_stars_repo_name: fredefox/javalette

max_issues_repo_name: fredefox/javalette

max_forks_repo_name: fredefox/javalette


text:
using GPUArrays, CuArrays, FillArrays
using CuArrays:threadIdx, blockIdx, blockDim
import Base:minimum, maximum, findmax, findmin

import Flux.Optimise: apply!, Descent, InvDecay

#####
# patch for maximum/minimum/findmax/findmin with dims==1
# considering that discrete actions in RL are usually small,
# the implementation here may be not that slow?
# watch https://github.com/JuliaGPU/CuArrays.jl/issues/304
#####
function findminmaxcol!(f, values::GPUArray{Float32, 2}, vals::GPUArray{Float32, 2}, args::GPUArray{CartesianIndex{2}, 1})
    function kernel(state, values, vals, args)
        i = linear_index(state)
        if i <= length(args)
            for j in axes(values, 1)
                @inbounds if f(values[j, i], vals[i])
                    vals[i] = values[j, i]
                    args[i] = CartesianIndex(j, i)
                end
            end
        end
        return
    end
    gpu_call(kernel, args, (values, vals, args))
    vals, args
end

function Base._findmax(values::CuArray{Float32, 2}, dims::Int)
    if dims == 1
        findminmaxcol!(Base.isgreater, values, values[1:1, :], cu(CartesianIndex.(1, axes(values, 2))))
    else
        error("unsupported yet for dims=$dims")
    end
end

function Base._findmin(values::CuArray{Float32, 2}, dims::Int)
    if dims == 1
        findminmaxcol!(Base.isless, values, values[1:1, :], cu(CartesianIndex.(1, axes(values, 2))))
    else
        error("unsupported yet for dims=$dims")
    end
end

function Base._maximum(x::CuArray{Float32, 2}, dims::Int)
    if dims == 1
        _findmax(x, dims)[1]
    else
        error("unsupported yet for dims=$dims")
    end
end

function Base._minimum(x::CuArray{Float32, 2}, dims::Int)
    if dims == 1
        _findmin(x, dims)[1]
    else
        error("unsupported yet for dims=$dims")
    end
end

#####
# Cartesian indexing of CuArray
#####

function Base.getindex(xs::CuArray{T, N}, indices::CuArray{CartesianIndex{N}, 1}) where {T, N}
  n = length(indices)
  ys = CuArray{T}(undef, n)

  if n > 0
    num_threads = min(n, 256)
    num_blocks = ceil(Int, n / num_threads)

    function kernel(ys::CuArrays.CuDeviceArray{T}, xs::CuArrays.CuDeviceArray{T}, indices)
        i = threadIdx().x + (blockIdx().x - 1) * blockDim().x

        if i <= length(ys)
            ind = indices[i]
            ys[i] = xs[ind]
        end

        return
    end

    CuArrays.@cuda blocks=num_blocks threads=num_threads kernel(ys, xs, indices)
  end

  return ys
end

function Base.setindex!(xs::CuArray{T, N}, v::CuArray{T}, indices::CuArray{CartesianIndex{N}, 1}) where {T, N}
    # @assert length(indices) == length(v) "$xs, $(size(xs)), $v, $(size(v)), $indices, $(size(indices))"
    n = length(indices)

    if n > 0
        num_threads = min(n, 256)
        num_blocks = ceil(Int, n / num_threads)

        function kernel(xs::CuArrays.CuDeviceArray{T}, indices, v)
            i = threadIdx().x + (blockIdx().x - 1) * blockDim().x

            if i <= length(indices)
                ind = indices[i]
                xs[ind] = v[i]
            end

            return
        end

        CuArrays.@cuda blocks=num_blocks threads=num_threads kernel(xs, indices, v)
    end
    return v
end

function Base.setindex!(xs::CuArray{T, N}, v::T, indices::CuArray{CartesianIndex{N}, 1}) where {T, N}
  n = length(indices)

  if n > 0
    num_threads = min(n, 256)
    num_blocks = ceil(Int, n / num_threads)

    function kernel(xs::CuArrays.CuDeviceArray{T}, indices, v)
        i = threadIdx().x + (blockIdx().x - 1) * blockDim().x

        if i <= length(indices)
            ind = indices[i]
            xs[ind] = v
        end

        return
    end

    CuArrays.@cuda blocks=num_blocks threads=num_threads kernel(xs, indices, v)
  end
  return v
end

Base.setindex!(xs::CuArray{T, N}, v::Fill{T}, indices::CuArray{CartesianIndex{N}, 1}) where {T, N} = setindex!(xs, v.value, indices)


text:
using MAT
using MultipleViewGeometry, Test
using MultipleViewGeometry.ModuleCostFunction
using MultipleViewGeometry.ModuleTypes
using BenchmarkTools
using StaticArrays, LinearAlgebra, SparseArrays

file = matopen("debug/mPts1.mat")
mpts1 = read(file,"mPts1")
close(file)

file = matopen("debug/mPts2.mat")
mpts2 = read(file,"mPts2")
close(file)

file = matopen("debug/covs.mat")
covs = read(file,"covs")
close(file)

file = matopen("debug/Ct.mat")
Ct = read(file,"Ct")
close(file)

npts, dim  = size(mpts1)

s²  = 1e-7
s = sqrt(s² )
Λ₁ =  [SMatrix{3,3}(s² * Matrix(Diagonal([1.0,1.0,0.0]))) for i = 1:npts]
Λ₂ =  [SMatrix{3,3}(s² *Matrix(Diagonal([1.0,1.0,0.0]))) for i = 1:npts]

𝒪 = [Point2DH(vcat(mpts1[i,:],[1])) for i =1:npts]
𝒪ʹ= [Point2DH(vcat(mpts2[i,:],[1])) for i =1:npts]

𝐅₀ = estimate(FundamentalMatrix(),DirectLinearTransform(), (𝒪, 𝒪ʹ))
𝐅 = estimate(FundamentalMatrix(),
                           FundamentalNumericalScheme(reshape(𝐅₀,9,1), 5, 1e-10),
                                                             (Λ₁,Λ₂), (𝒪, 𝒪ʹ))
𝐟 = reshape(𝐅,9,1)
𝐟 = 𝐟 / norm(𝐟)

# # Validate the covariance matrix of an estimate based on the AML cost function.
C1 = covariance_matrix(AML(), HessianApproximation(), FundamentalMatrix(), reshape(𝐅,9,1), (Λ₁,Λ₂), (𝒪 , 𝒪ʹ))
C2 = covariance_matrix(AML(), CanonicalApproximation(), FundamentalMatrix(),  𝐟, (Λ₁,Λ₂), (𝒪 , 𝒪ʹ))

@test norm(C1-Ct) / norm(Ct) * 100 < 0.6
@test norm(C2-Ct) / norm(Ct) * 100 < 0.6


text:
function define_indices_d(::Type{Gridded{Constant}}, d, pad)
    symix, symx = Symbol("ix_",d), Symbol("x_",d)
    symk, symkix = Symbol("k_",d), Symbol("kix_",d)
    # Choose the closer of k[ix] and k[ix+1]
    quote
        $symix = clamp($symix, 1, size(itp, $d)-1)
        $symkix = $symk[$symix]
        cmp = $symk[$symix+1]
        $symix = abs($symx-$symkix) < abs($symx-cmp) ? $symix : $symix+1
    end
end

function coefficients(::Type{Gridded{Constant}}, N, d)
    sym, symx = Symbol("c_",d), Symbol("x_",d)
    :($sym = 1)
end

function gradient_coefficients(::Type{Gridded{Constant}}, N, d)
    sym, symx = Symbol("c_",d), Symbol("x_",d)
    :($sym = 0)
end

function index_gen{IT<:DimSpec{Gridded}}(::Type{Gridded{Constant}}, ::Type{IT}, N::Integer, offsets...)
    if (length(offsets) < N)
        d = length(offsets)+1
        sym = Symbol("c_", d)
        return :($sym * $(index_gen(IT, N, offsets..., 0)))
    else
        indices = [offsetsym(offsets[d], d) for d = 1:N]
        return :(itp.coefs[$(indices...)])
    end
end


text:
# Test model-fitting
function build(theta)
    F = diagm(theta[1])
    V = diagm(exp(theta[2]))
    G = reshape(theta[3:5], 3, 1)
    W = diagm(exp(theta[6:8]))
    x0 = [theta[9]]
    P0 = diagm(1e7)
    StateSpaceModel(F, V, G, W, x0, P0)
end

@testset "Parameter Estimation" begin

    #=
    context("Nonlinear solver") do
        mod1 = build_model()
        x, y = simulate(mod1, 100)
        fit(y, build, zeros(9))
    end
    =#

    #=
    context("Unspecified constraints default behaviour") do

        context("Without inputs") do
            mod1 = build_model()
            x, y = simulate(mod1, 100)
            #p, fitmod = fit(y, build(randn(9)))
        end

        context("With inputs") do
            #test
        end

    end
    =#

    @testset "Parameter estimation with constraints" begin

        m, b, dt = 5, 2, .1
        t = 0:dt:200
        y_true = m*t + b

        @testset "Without exogenous inputs" begin

            s = 2
            y_noisy = y_true + s*randn(length(t))

            lm = ParametrizedSSM(
                    parametrize_none([1 0.1; 0 1]), #A
                    parametrize_none(ones(1,1)), #Q
                    parametrize_none([1. 0]), #C
                    parametrize_diag([1.]), #R
                    parametrize_none([1. 1.]'), #x1
                    parametrize_none(100*eye(2)), #P1
                    G = _->zeros(2,1) #G
                )
            lm_params = SSMParameters(1., R=rand(1))
            fitm_params, fitm = fit(y_noisy, lm, lm_params)
            @test sqrt(fitm_params.R[1]) ≈ s atol=0.1

            A = [.5 .1 .4; .25 .8 .5; .25 .1 .1]
            model = StateSpaceModel(A, diagm([.01,.01,.01]),
                                        eye(3), zeros(3,3), ones(3)/3, 0eye(3))
            _, y = simulate(model, 100)
            lm = ParametrizedSSM(
                    ParametrizedMatrix([0, .25, .75, .1, 0, .9, 0, .5, .5], [
                        1. 0 0; 0 0 0; -1 0 0; 0 0 0; 0 1 0; 0 -1 0; 0 0 1; 0 0 0; 0 0 -1
                    ], (3,3)), #A
                    parametrize_diag(ones(3)), #Q
                    parametrize_none(eye(3)), #C
                    parametrize_none(ones(1,1)), #R
                    parametrize_none(ones(3,1)/3), #x1
                    parametrize_none(100*eye(3)), #P1
                    H = _->zeros(3,1) #H
                )

            lm_params = SSMParameters(1., A=rand(3), Q=rand(3))
            fitm_params, fitm = fit(y, lm, lm_params)
            @test fitm_params.A ≈ [.5, .8, .4] atol=0.2
            @test fitm_params.Q ≈ [.01, .01, .01] atol=0.05

            coeffs = randn(3)
            X = randn(200,3)
            Y = X * coeffs + .4*randn(200)
            lm = ParametrizedSSM(
                parametrize_none(eye(3)), #A
                parametrize_none(eye(1)), #Q
                parametrize_none(eye(3)), #C2
                parametrize_diag(ones(1)), #R
                parametrize_full(zeros(3,1)), #x1
                parametrize_none(eye(3)), #S
                G=_->zeros(3,1), C1=t->X[t:t, :]
            )

            lm_params = SSMParameters(1., R=ones(1), x1=100*randn(3))
            fitm_params, fitm = fit(Y, lm, lm_params)
            @test fitm_params.x1   ≈ coeffs atol=0.1
            @test fitm_params.R[1] ≈ .16    atol=0.05

            lm = ParametrizedSSM(
                parametrize_none(eye(3)), #A
                parametrize_none(eye(1)), #Q
                parametrize_none(eye(3)), #C2
                parametrize_none(.16*ones(1,1)), #R
                parametrize_none(coeffs+[.2,.3,.5].*randn(3,1)), #x1
                parametrize_diag(ones(3)), #S
                G=_->zeros(3,1), C1=t->X[t:t, :]
            )
            lm_params = SSMParameters(1., S=100*ones(3))
            fitm_params, fitm = fit(Y, lm, lm_params)
            @test fitm_params.S ≈ [.2,.3,.5].^2 atol=0.15

        end

        @testset "With exogenous inputs" begin

            s1, s2, s3 = 1., 2., 3.
            input = 100.*[sin.(t./2 .+ 0.1) sin.(t./4 .+ .1) cos.(t./2 .+ .1) cos.(t./4 .+ .1)] .+ 10
            y_noisy = [y_true zeros(length(t)) -y_true] +
                        [input[:,1] + input[:,2] input[:,1] + input[:,3] input[:,3] + input[:,4]] +
                        [s1 s2 s3] .* randn(length(t), 3)

            lm = ParametrizedSSM(
                  parametrize_none([1 0.1; 0 1]), #A
                  parametrize_none(eye(1)), #Q
                  parametrize_none([1. 0; 0 0; -1 0]), #C
                  parametrize_diag(ones(3)), #R
                  parametrize_none([2. 5.]'), #x1
                  parametrize_none(0.001*eye(2)), #P1
                  B2=parametrize_none(zeros(2,4)), #B
                  G=_->zeros(2,1), #G
                  D2=parametrize_full(randn(3,4)) #D
            )
            lm_params = SSMParameters(1., R=rand(3), D=randn(12))
            fitm_params, fitm = fit(y_noisy, lm, lm_params, u=input)
            @test fitm_params.D       ≈ vec([1. 1 0 0; 1 0 1 0; 0 0 1 1]) atol=.1
            @test sqrt.(fitm_params.R) ≈ [s1, s2, s3] atol=.15

            y_noisy = [0 0 0;
                  [input[:,1]+input[:,2] input[:,1]+input[:,3] input[:,3]+input[:,4]][1:end-1, :]] +
                  [s1 s2 s3] .* randn(length(t), 3)

            lm = ParametrizedSSM(
                  parametrize_none(zeros(3,3)), #A
                  parametrize_diag(ones(3)), #Q
                  parametrize_none(eye(3)), #C
                  parametrize_none(eye(1)), #R
                  parametrize_none(zeros(3,1)), #x1
                  parametrize_none(0.001*eye(3)), #P1
                  B2=parametrize_none([1. 1 0 0; 1 0 1 0; 0 0 1 1]), #B
                  H=_->zeros(3,1), #H
                  D2=parametrize_none(zeros(3,4)), #D
            )
            lm_params = SSMParameters(1., Q=rand(3))
            fitm_params, fitm = fit(y_noisy, lm, lm_params, u=input)
            @test sqrt.(fitm_params.Q) ≈ [s1, s2, s3] atol=.1

            lm = ParametrizedSSM(
                  parametrize_none(zeros(3,3)), #A
                  parametrize_none(diagm([s1, s2, s3])), #Q
                  parametrize_none(eye(3)), #C
                  parametrize_none(eye(1)), #R
                  parametrize_none(zeros(3,1)), #x1
                  parametrize_none(0.001*eye(3)), #P1
                  B2=parametrize_full(randn(3,4)), #B
                  H=_->zeros(3,1), #H
                  D2=parametrize_none(zeros(3,4)), #D
            )
            lm_params = SSMParameters(1., B=randn(12))
            fitm_params, fitm = fit(y_noisy, lm, lm_params, u=input)
            @test fitm_params.B ≈ vec([1. 1 0 0; 1 0 1 0; 0 0 1 1]) atol=.11

            y_noisy = [0; .5input[1:end-1,1] - .8input[1:end-1,2] - .3input[1:end-1,3] + .7input[1:end-1,4]] + s2*randn(length(t))

            lm = ParametrizedSSM(
                  parametrize_none(zeros(4,4)), #A
                  parametrize_none(eye(1)), #Q
                  parametrize_full(ones(1,4)), #C
                  parametrize_full(diagm(s2)), #R
                  parametrize_none(zeros(4,1)), #x1
                  parametrize_none(0.001*eye(4)), #P1
                  B2=parametrize_none(eye(4)), #B
                  G=_->zeros(4,1), #H
                  D2=parametrize_none(zeros(1,4)), #D
            )
            lm_params = SSMParameters(1., C=randn(4), R=rand(1))
            fitm_params, fitm = fit(y_noisy, lm, lm_params, u=input)
            @test fitm_params.C ≈ [.5, -.8, -.3, .7] atol=.1
            @test sqrt(fitm_params.R[1]) ≈ s2 atol=.1

        end

    end

end


text:
using FeynmanDiagram
using Lehmann
using Test, LinearAlgebra, Random, StaticArrays, Printf, Parameters, Documenter
using AbstractTrees

# @testset "doctest" begin
#     DocMeta.setdocmeta!(FeynmanDiagram, :DocTestSetup, :(using FeynmanDiagram); recursive = true)
#     doctest(FeynmanDiagram)
# end

include("common.jl")
include("diagram_tree.jl")
include("expression_tree.jl")
include("parquet_builder.jl")


text:
export Positions, GridPositions, RegularGridPositions, ChebyshevGridPositions,
       MeanderingGridPositions, UniformRandomPositions, ArbitraryPositions,
       SphericalTDesign, BreakpointGridPositions
export SpatialDomain, AxisAlignedBox, Ball
export loadTDesign, getPermutation
export fieldOfView, fieldOfViewCenter, shape
export idxToPos, posToIdx, posToLinIdx, spacing, isSubgrid, deriveSubgrid, toDict

abstract type Positions end
abstract type GridPositions<:Positions end

function Positions(file::HDF5.File)
  if haskey(file, "/positionsBreakpoint")
    return BreakpointGridPositions(file)
  end

  typ = read(file, "/positionsType")
  if typ == "RegularGridPositions"
    positions = RegularGridPositions(file)
  elseif typ == "ChebyshevGridPositions"
    positions = ChebyshevGridPositions(file)
  elseif typ == "SphericalTDesign"
    positions = SphericalTDesign(file)
  elseif typ == "UniformRandomPositions"
    positions = UniformRandomPositions(file)
  elseif typ == "ArbitraryPositions"
    positions = ArbitraryPositions(file)
  else
    throw(ErrorException("No grid found to load from $file"))
  end

  if haskey(file, "/positionsMeandering") && typ in ["RegularGridPositions","ChebyshevGridPositions"] && read(file, "/positionsMeandering") == Int8(1)
    positions = MeanderingGridPositions(positions)
  end

  return positions
end

function Positions(params::Dict)
  if haskey(params, "positionsBreakpoint")
    return BreakpointGridPositions(params)
  end

  typ = params["positionsType"]
  if typ == "RegularGridPositions"
    positions = RegularGridPositions(params)
  elseif typ == "ChebyshevGridPositions"
    positions = ChebyshevGridPositions(params)
  elseif typ == "SphericalTDesign"
    positions = SphericalTDesign(params)
  elseif typ == "UniformRandomPositions"
    positions = UniformRandomPositions(params)
  elseif typ == "ArbitraryPositions"
    positions = ArbitraryPositions(params)
  else
    throw(ErrorException("No grid found to load from dict $params"))
  end

  if haskey(params, "positionsMeandering") && typ in ["RegularGridPositions","ChebyshevGridPositions"] && params["positionsMeandering"]
    positions = MeanderingGridPositions(positions)
  end

  return positions
end

# Cartesian grid
mutable struct RegularGridPositions{T} <: GridPositions where {T<:Unitful.Length}
  shape::Vector{Int}
  fov::Vector{T}
  center::Vector{T}
  sign::Vector{Int}
end

function range(grid::RegularGridPositions, dim::Int)
  if grid.shape[dim] > 1
    sp = spacing(grid)
    return range(grid.center[dim] - grid.fov[dim]/2 + sp[dim]/2,
                 step=sp[dim], length=grid.shape[dim])
  else
    return 1:1
  end
end

RegularGridPositions(shape, fov, center) = RegularGridPositions(shape, fov, center, ones(Int,length(shape)))

function RegularGridPositions(file::HDF5.File)
  shape = read(file, "/positionsShape")
  fov = read(file, "/positionsFov")*Unitful.m
  center = read(file, "/positionsCenter")*Unitful.m

  return RegularGridPositions(shape, fov, center)
end

function RegularGridPositions(params::Dict)
  shape = params["positionsShape"]

  fov = params["positionsFov"]
  if !(eltype(fov) <: Quantity)
    fov = fov.*Unitful.m
  end

  center = params["positionsCenter"]
  if !(eltype(center) <: Quantity)
    center = center.*Unitful.m
  end

  return RegularGridPositions(shape, fov, center)
end

# Find a joint grid
function RegularGridPositions(positions::Vector{T}) where T<:RegularGridPositions
  posMin = positions[1].center .- 0.5*positions[1].fov
  posMax = positions[1].center .+ 0.5*positions[1].fov
  minSpacing = spacing(positions[1])
  for position in positions
    sp = spacing(position)
    for d=1:length(posMin)
      posMin[d] = min(posMin[d], position.center[d] - 0.5*position.fov[d])
      posMax[d] = max(posMax[d], position.center[d] + 0.5*position.fov[d])
      minSpacing[d] = min(minSpacing[d],sp[d])
    end
  end
  center = (posMin .+ posMax)/2
  fov = posMax .- posMin
  shape = round.(Int64,fov./minSpacing)
  fov = shape .* minSpacing
  return RegularGridPositions(shape, fov, center)
end

function isSubgrid(grid::RegularGridPositions, subgrid::RegularGridPositions)
  if any(fieldOfView(grid) .- fieldOfView(subgrid) .< 0) ||
     any(spacing(grid) .!= spacing(subgrid))
    return false
  else
    centerPosIdx = posToIdxFloat(grid,subgrid[ones(Int,length(subgrid.shape))])
    return all(isapprox.(centerPosIdx, round.(Int,centerPosIdx) ,rtol=1e-5))
  end
end

function deriveSubgrid(grid::RegularGridPositions, subgrid::RegularGridPositions)
  minI = ones(Int,length(subgrid.shape))
  maxI = copy(subgrid.shape)
  for d=1:length(minI)
    if subgrid.sign[d] == -1
      minI[d] = subgrid.shape[d]-minI[d]+1
      maxI[d] = subgrid.shape[d]-maxI[d]+1
    end
  end
  minPos = subgrid[ minI ]
  maxPos = subgrid[ maxI ]

  minIdx = posToIdx(grid,minPos)
  maxIdx = posToIdx(grid,maxPos)
  #shp = maxIdx-minIdx+ones(Int,length(subgrid.shape))
  shp = shape(subgrid)
  #center = (grid[minIdx].+grid[maxIdx])/2
  # TODO round properly
  center = (grid[minIdx].+grid[minIdx.+shp.-1])/2
  fov = shp.*spacing(grid)
  return RegularGridPositions(shp,fov,center,subgrid.sign)
end

function write(file::HDF5.File, positions::RegularGridPositions)
  write(file,"/positionsType", "RegularGridPositions")
  write(file, "/positionsShape", positions.shape)
  write(file, "/positionsFov", Float64.(ustrip.(uconvert.(Unitful.m, positions.fov))) )
  write(file, "/positionsCenter", Float64.(ustrip.(uconvert.(Unitful.m, positions.center))) )
end

function toDict(positions::RegularGridPositions)
  params = Dict{String,Any}()
  params["positionsType"] = "RegularGridPositions"
  params["positionsShape"] = positions.shape
  params["positionsFov"] = Float64.(ustrip.(uconvert.(Unitful.m, positions.fov)))
  params["positionsCenter"] = Float64.(ustrip.(uconvert.(Unitful.m, positions.center)))
  return params
end

function getindex(grid::RegularGridPositions, i::Integer)
  if i>length(grid) || i<1
     throw(BoundsError(grid,i))
  end

  #idx = collect(ind2sub(tuple(shape(grid)...), i))
  if length(grid.shape) == 1 #Very ugly but improves compile time
    idx = [i]
  elseif length(grid.shape) == 2
    idx = collect(Tuple((CartesianIndices(tuple(grid.shape[1],grid.shape[2])))[i]))
  elseif length(grid.shape) == 3
    idx = collect(Tuple((CartesianIndices(tuple(grid.shape[1],grid.shape[2],grid.shape[3])))[i]))
  else
    idx = collect(Tuple((CartesianIndices(tuple(grid.shape...)))[i]))
  end

  for d=1:length(idx)
    if grid.sign[d] == -1
      idx[d] = grid.shape[d]-idx[d]+1
    end
  end
  return ((-shape(grid).+(2 .*idx.-1))./shape(grid)).*fieldOfView(grid)./2 + fieldOfViewCenter(grid)
end

function getindex(grid::RegularGridPositions, idx::Vector{T}) where T<:Number
  for d=1:length(idx)
    if grid.sign[d] == -1
      idx[d] = grid.shape[d]-idx[d]+1
    end
  end
  return 0.5.*fieldOfView(grid) .* (-1 .+ (2 .* idx .- 1) ./ shape(grid)) .+ fieldOfViewCenter(grid)
end

function posToIdxFloat(grid::RegularGridPositions,pos::Vector)
  idx = 0.5 .* (shape(grid) .* ((pos .- fieldOfViewCenter(grid)) ./
              ( 0.5 .* fieldOfView(grid) ) .+ 1) .+ 1)
  return idx
end

function posToIdx(grid::RegularGridPositions,pos::Vector)
  idx = round.(Int64, posToIdxFloat(grid,pos))
  for d=1:length(idx)
    if grid.sign[d] == -1
      idx[d] = grid.shape[d]-idx[d]+1
    end
  end
  return idx
end

function posToLinIdx(grid::RegularGridPositions,pos::Vector)
  return (LinearIndices(tuple(shape(grid)...)))[posToIdx(grid,pos)...]
end

# Chebyshev Grid
mutable struct ChebyshevGridPositions{S,T} <: GridPositions where {S,T<:Unitful.Length}
  shape::Vector{Int}
  fov::Vector{S}
  center::Vector{T}
end

function write(file::HDF5.File, positions::ChebyshevGridPositions)
  write(file,"/positionsType", "ChebyshevGridPositions")
  write(file, "/positionsShape", positions.shape)
  write(file, "/positionsFov", Float64.(ustrip.(uconvert.(Unitful.m, positions.fov))) )
  write(file, "/positionsCenter", Float64.(ustrip.(uconvert.(Unitful.m, positions.center))) )
end

function toDict(positions::ChebyshevGridPositions)
  params = Dict{String,Any}()
  params["positionsType"] = "ChebyshevGridPositions"
  params["positionsShape"] = positions.shape
  params["positionsFov"] = Float64.(ustrip.(uconvert.(Unitful.m, positions.fov)))
  params["positionsCenter"] = Float64.(ustrip.(uconvert.(Unitful.m, positions.center)))
  return params
end

function ChebyshevGridPositions(file::HDF5.File)
  shape = read(file, "/positionsShape")
  fov = read(file, "/positionsFov")*Unitful.m
  center = read(file, "/positionsCenter")*Unitful.m
  return ChebyshevGridPositions(shape,fov,center)
end

function ChebyshevGridPositions(params::Dict)
  shape = params["positionsShape"]

  fov = params["positionsFov"]
  if !(eltype(fov) <: Quantity)
    fov = fov.*Unitful.m
  end

  center = params["positionsCenter"]
  if !(eltype(center) <: Quantity)
    center = center.*Unitful.m
  end

  return ChebyshevGridPositions(shape,fov,center)
end

function getindex(grid::ChebyshevGridPositions, i::Integer)
  if i>length(grid) || i<1
    throw(BoundsError(grid,i))
  else
    idx = collect(Tuple(CartesianIndices(tuple(shape(grid)...))[i]))
    return -cos.((idx .- 0.5) .* pi ./ shape(grid)) .* fieldOfView(grid) ./ 2 .+ fieldOfViewCenter(grid)
  end
end

# Meander regular grid positions
mutable struct MeanderingGridPositions{T} <: GridPositions where {T<:GridPositions}
  grid::T
end

function MeanderingGridPositions(file::HDF5.File)
  typ = read(file, "/positionsType")
  if typ == "RegularGridPositions"
    grid = RegularGridPositions(file)
    return MeanderingGridPositions(grid)
  elseif typ == "ChebyshevGridPositions"
    grid = ChebyshevGridPositions(file)
    return MeanderingGridPositions(grid)
  end
end

function MeanderingGridPositions(params::Dict)
  typ = params["positionsType"]
  if typ == "RegularGridPositions"
    grid = RegularGridPositions(params)
    return MeanderingGridPositions(grid)
  elseif typ == "ChebyshevGridPositions"
    grid = ChebyshevGridPositions(params)
    return MeanderingGridPositions(grid)
  end
end

function write(file::HDF5.File, positions::MeanderingGridPositions)
  write(file,"/positionsMeandering", Int8(1))
  write(file, positions.grid)
end

function toDict(positions::MeanderingGridPositions)
  params = toDict(positions.grid)
  params["positionsMeandering"] = true
  return params
end

function indexPermutation(grid::MeanderingGridPositions, i::Integer)
  dims = tuple(shape(grid)...)
  ndims = length(dims)
  idx = collect(Tuple(CartesianIndices(dims)[i]))
    for d=1:ndims-1
      if isodd(sum(idx[d+1:end])-(ndims-d))
      idx[d] = -idx[d] + shape(grid)[d] + 1
    end
  end
  linidx = (LinearIndices(dims))[idx...]
end

function getindex(grid::MeanderingGridPositions, i::Integer)
  iperm = indexPermutation(grid,i)
  return grid.grid[iperm]
end

function getPermutation(grid::MeanderingGridPositions)
  N = length(grid)
  perm = Array{Int}(undef,N)

  for i in eachindex(perm)
    perm[i] = indexPermutation(grid,i)
  end
  return vec(perm)
end

mutable struct BreakpointGridPositions{T,S} <: GridPositions where {T<:GridPositions}
  grid::T
  breakpointIndices::Vector{Int64}
  breakpointPosition::Vector{S}
end

function BreakpointGridPositions(file::HDF5.File)
  typ = read(file, "/positionsType")
  breakpointPosition = read(file, "/positionsBreakpoint") * Unitful.m
  breakpointIndices = read(file, "/indicesBreakpoint")

  if typ == "MeanderingGridPositions"
    grid = MeanderingGridPositions(file)
    return BreakpointGridPositions(grid,breakpointIndices, breakpointPosition)
  elseif typ == "RegularGridPositions"
    grid = RegularGridPositions(file)
    return BreakpointGridPositions(grid, breakpointIndices, breakpointPosition)
  elseif typ == "ChebyshevGridPositions"
    grid = ChebyshevGridPositions(file)
    return BreakpointGridPositions(grid, breakpointIndices, breakpointPosition)
  end
end

function BreakpointGridPositions(params::Dict)
  typ = params["positionsType"]

  breakpointPosition = params["positionsBreakpoint"]
  if !(eltype(breakpointPosition) <: Quantity)
    breakpointPosition = breakpointPosition.*Unitful.m
  end

  breakpointIndices = params["indicesBreakpoint"]

  if typ == "MeanderingGridPositions"
    grid = MeanderingGridPositions(params)
    return BreakpointGridPositions(grid, breakpointIndices, breakpointPosition)
  elseif typ == "RegularGridPositions"
    grid = RegularGridPositions(params)
    return BreakpointGridPositions(grid, breakpointIndices, breakpointPosition)
  elseif typ == "ChebyshevGridPositions"
    grid = ChebyshevGridPositions(params)
    return BreakpointGridPositions(grid, breakpointIndices, breakpointPosition)
  end
end

function write(file::HDF5.File, positions::BreakpointGridPositions)
  write(file,"/positionsBreakpoint", Float64.(ustrip.(uconvert.(Unitful.m, positions.breakpointPosition))))
  write(file,"/indicesBreakpoint", positions.breakpointIndices)
  write(file, positions.grid)
end

function toDict(positions::BreakpointGridPositions)
  params = toDict(positions.grid)
  params["positionsBreakpoint"] = Float64.(ustrip.(uconvert.(Unitful.m, positions.breakpointPosition)))
  params["indicesBreakpoint"] = positions.breakpointIndices
  return params
end

function getmask(grid::BreakpointGridPositions)
  bgind=grid.breakpointIndices
  mask = zeros(Bool, length(grid.grid)+length(bgind))
  mask[bgind] .= true
  return mask
end

function getindex(grid::BreakpointGridPositions, i::Integer)

  bgind=grid.breakpointIndices

  if i>(length(grid.grid)+length(bgind)) || i<1
    return throw(BoundsError(grid,i))
  elseif any(i .== bgind)
    return grid.breakpointPosition
  else
    pastBgind = sum(i .> bgind)
    return grid.grid[i-pastBgind]
  end
end

# Uniform random distributed positions
abstract type SpatialDomain end

struct AxisAlignedBox <: SpatialDomain
  fov::Vector{S} where {S<:Unitful.Length}
  center::Vector{T} where {T<:Unitful.Length}
end

function write(file::HDF5.File, domain::AxisAlignedBox)
  write(file, "/positionsDomain", "AxisAlignedBox")
  write(file, "/positionsDomainFieldOfView", Float64.(ustrip.(uconvert.(Unitful.m, domain.fov))) )
  write(file, "/positionsDomainCenter", Float64.(ustrip.(uconvert.(Unitful.m, domain.center))) )
end

function AxisAlignedBox(file::HDF5.File)
  fov = read(file, "/positionsDomainFieldOfView")*Unitful.m
  center = read(file, "/positionsDomainCenter")*Unitful.m
  return AxisAlignedBox(fov,center)
end

struct Ball <: SpatialDomain
  radius::S where {S<:Unitful.Length}
  center::Vector{T} where {T<:Unitful.Length}
end

function write(file::HDF5.File, domain::Ball)
  write(file, "/positionsDomain", "Ball")
  write(file, "/positionsDomainRadius", Float64.(ustrip.(uconvert.(Unitful.m, domain.radius))) )
  write(file, "/positionsDomainCenter", Float64.(ustrip.(uconvert.(Unitful.m, domain.center))) )
end

function Ball(file::HDF5.File)
  radius = read(file, "/positionsDomainRadius")*Unitful.m
  center = read(file, "/positionsDomainCenter")*Unitful.m
  return Ball(radius,center)
end

mutable struct UniformRandomPositions{T} <: Positions where {T<:SpatialDomain}
  N::UInt
  seed::UInt32
  domain::T
end

radius(rpos::UniformRandomPositions{Ball}) = rpos.domain.radius
seed(rpos::UniformRandomPositions) = rpos.seed

function getindex(rpos::UniformRandomPositions{AxisAlignedBox}, i::Integer)
  if i>length(rpos) || i<1
    throw(BoundsError(rpos,i))
  else
    # make sure Positions are randomly generated from given seed
    mersenneTwister = MersenneTwister(seed(rpos))
    rP = rand(mersenneTwister, 3, i)[:,i]
    return (rP.-0.5).*fieldOfView(rpos)+fieldOfViewCenter(rpos)
  end
end

function getindex(rpos::UniformRandomPositions{Ball}, i::Integer)
  if i>length(rpos) || i<1
    throw(BoundsError(rpos,i))
  else
    # make sure Positions are randomly generated from given seed
    mersenneTwister = MersenneTwister(seed(rpos))
    D = rand(mersenneTwister, i)[i]
    P = randn(mersenneTwister, 3, i)[:,i]
    return radius(rpos)*D^(1/3)*normalize(P)+fieldOfViewCenter(rpos)
  end
end

function write(file::HDF5.File, positions::UniformRandomPositions{T}) where {T<:SpatialDomain}
  write(file, "/positionsType", "UniformRandomPositions")
  write(file, "/positionsN", positions.N)
  write(file, "/positionsSeed", positions.seed)
  write(file, positions.domain)
end

function UniformRandomPositions(file::HDF5.File)
  N = read(file, "/positionsN")
  seed = read(file, "/positionsSeed")
  dom = read(file,"/positionsDomain")
  if dom=="Ball"
    domain = Ball(file)
    return UniformRandomPositions(N,seed,domain)
  elseif dom=="AxisAlignedBox"
    domain = AxisAlignedBox(file)
    return UniformRandomPositions(N,seed,domain)
  else
    throw(ErrorException("No method to read domain $domain"))
  end
end

# TODO fix conversion methods
#=
function convert(::Type{UniformRandomPositions}, N::Integer,seed::UInt32,fov::Vector{S},center::Vector{T}) where {S,T<:Unitful.Length}
  if N<1
    throw(DomainError())
  else
    uN = convert(UInt,N)
    return UniformRandomPositions(uN,seed,fov,center)
  end
end

function convert(::Type{UniformRandomPositions}, N::Integer,fov::Vector,center::Vector)
  return UniformRandomPositions(N,rand(UInt32),fov,center)
end
=#

# General functions for handling grids
fieldOfView(grid::GridPositions) = grid.fov
fieldOfView(grid::UniformRandomPositions{AxisAlignedBox}) = grid.domain.fov
fieldOfView(mgrid::MeanderingGridPositions) = fieldOfView(mgrid.grid)
fieldOfView(bgrid::BreakpointGridPositions) = fieldOfView(bgrid.grid)
shape(grid::GridPositions) = grid.shape
shape(mgrid::MeanderingGridPositions) = shape(mgrid.grid)
shape(bgrid::BreakpointGridPositions) = shape(bgrid.grid)
fieldOfViewCenter(grid::GridPositions) = grid.center
fieldOfViewCenter(grid::UniformRandomPositions) = grid.domain.center
fieldOfViewCenter(mgrid::MeanderingGridPositions) = fieldOfViewCenter(mgrid.grid)
fieldOfViewCenter(bgrid::BreakpointGridPositions) = fieldOfViewCenter(bgrid.grid)

spacing(grid::GridPositions) = grid.fov ./ grid.shape

mutable struct SphericalTDesign{S,V} <: Positions where {S,V<:Unitful.Length}
  T::Unsigned
  radius::S
  positions::Matrix
  center::Vector{V}
end

function SphericalTDesign(file::HDF5.File)
  T = read(file, "/positionsTDesignT")
  N = read(file, "/positionsTDesignN")
  radius = read(file, "/positionsTDesignRadius")*Unitful.m
  center = read(file, "/positionsCenter")*Unitful.m
  return loadTDesign(Int64(T),N,radius,center)
end

function write(file::HDF5.File, positions::SphericalTDesign)
  write(file,"/positionsType", "SphericalTDesign")
  write(file, "/positionsTDesignT", positions.T)
  write(file, "/positionsTDesignN", size(positions.positions,2))
  write(file, "/positionsTDesignRadius", Float64.(ustrip.(uconvert.(Unitful.m, positions.radius))) )
  write(file, "/positionsCenter", Float64.(ustrip.(uconvert.(Unitful.m, positions.center))) )
end

function toDict(positions::SphericalTDesign)
  params = Dict{String,Any}()
  params["positionsType"] = "SphericalTDesign"
  params["positionsTDesignT"] = positions.T
  params["positionsTDesignN"] = size(positions.positions,2)
  params["positionsTDesignRadius"] = Float64.(ustrip.(uconvert.(Unitful.m, positions.radius)))
  params["positionsCenter"] = Float64.(ustrip.(uconvert.(Unitful.m, positions.center)))
  return params
end

getindex(tdes::SphericalTDesign, i::Integer) = tdes.radius.*tdes.positions[:,i] + tdes.center

"""
    loadTDesign(t::Int64, N::Int64, radius::S=10Unitful.mm, center::Vector{V}=[0.0,0.0,0.0]Unitful.mm, filename::String=joinpath(@__DIR__, "TDesigns.hd5")) where {S,V<:Unitful.Length}
*Description:* Returns the t-design array for chosen degree t and number of points N\\
\\
*Input:*
- `t` - degree
- `N` - number of points
- `radius` - radius of the sphere (default: 10mm)
- `center` - center of the sphere (default: [0.0,0.0,0.0]mm)
- `filename` - name of the file containing the t-designs (default loads TDesign.hd5)

*Output:*
- t-design of type SphericalTDesign in Cartesian coordinates containing t, radius, center and positions (which are located on the unit sphere unless `getindex(tdes,i)` is used)
"""
function loadTDesign(t::Int64, N::Int64, radius::S=10Unitful.mm, center::Vector{V}=[0.0,0.0,0.0]Unitful.mm, filename::String=joinpath(@__DIR__, "TDesigns.hd5")) where {S,V<:Unitful.Length}
  h5file = h5open(filename, "r")
  address = "/$t-Design/$N"

  if haskey(h5file, address)
    positions = copy(transpose(read(h5file, address)))
    return SphericalTDesign(UInt(t),radius,positions, center)
  else
    if haskey(h5file, "/$t-Design/")
      Ns = Int[]
      for N in keys(read(h5file, string("/$t-Design")))
	push!(Ns,parse(Int,N))
      end
      sort!(Ns)
      @info "No spherical $t-Design with $N points availible!\nThere are spherical $t-Designs with following N:" Ns
      throw(DomainError(1))
    else
      ts = Int[]
      for d in keys(read(h5file))
	m = match(r"(\d{1,})-(Design)",d)
	if m != nothing
	  push!(ts,parse(Int,m[1]))
        end
      end
      sort!(ts)
      @info "No spherical $t-Design availible!\n Choose another t."
      throw(DomainError(1))
    end
  end
end

# Unstructured collection of positions
mutable struct ArbitraryPositions{T} <: Positions where {T<:Unitful.Length}
  positions::Matrix{T}
end

getindex(apos::ArbitraryPositions, i::Integer) = apos.positions[:,i]

function ArbitraryPositions(grid::GridPositions)
  T = eltype(grid.fov)
  positions = zeros(T,3,length(grid))
  for i=1:length(grid)
    positions[:,i] = grid[i]
  end
  return ArbitraryPositions(positions)
end

function write(file::HDF5.File, apos::ArbitraryPositions,)
  write(file,"/positionsType", "ArbitraryPositions")
  write(file, "/positionsPositions", Float64.(ustrip.(uconvert.(Unitful.m, apos.positions))) )
end

function ArbitraryPositions(file::HDF5.File)
  pos = read(file, "/positionsPositions")*Unitful.m
  return ArbitraryPositions(pos)
end

# TODO: Specialize to make it faster
function Base.:(==)(val1::Positions, val2::Positions)
  for (x, y) in zip(val1, val2)
    if !(upreferred.(x) == upreferred.(y))
      return false
    end
  end

  return true
end

# fuction related to looping
length(tdes::SphericalTDesign) = size(tdes.positions,2)
length(apos::ArbitraryPositions) = size(apos.positions,2)
length(grid::GridPositions) = prod(grid.shape)
length(rpos::UniformRandomPositions) = rpos.N
length(mgrid::MeanderingGridPositions) = length(mgrid.grid)
length(bgrid::BreakpointGridPositions) = length(bgrid.grid)+length(bgrid.breakpointIndices)

start_(grid::Positions) = 1
next_(grid::Positions,state) = (grid[state],state+1)
done_(grid::Positions,state) = state > length(grid)
iterate(grid::Positions, s=start_(grid)) = done_(grid, s) ? nothing : next_(grid, s)

include("Interpolation.jl")


text:
using LinearAlgebra
using DataStructures

struct P{T<:Real}
    x::T
    y::T
end

function Base.:+(A::P,B::P)
 return P(A.x+B.x,A.y+B.y)
end
function Base.:-(A::P,B::P)
 return P(A.x-B.x,A.y-B.y)
end

function LinearAlgebra.norm(P::P)
    sqrt(P.x^2+P.y^2)
end

struct Polygon{T}
    name::Symbol
    edges::T
end

function Base.show(io::IO,  P::Polygon)
    print(io, P.name)
end
universe = [Polygon(:box,(P(0.0,0.0),P(0.0,5.0),P(10.0,5.0),P(10.0,0.0))),
            Polygon(:square,(P(0.0,0.0),P(0.0,10.0),P(10.0,10.0),P(10.0,0.0))),
            Polygon(:isosceles,(P(0.0,0.0),P(5.0,5.0),P(10.0,0.0))),
            Polygon(:right_triangle,(P(0.0,0.0),P(4.0,3.0),P(4.0,0.0))),
            Polygon(:trapezoid,(P(0.0,0.0),P(5.0,5.0),P(20.0,5.0),P(25.0,0.0))),
            Polygon(:parallelogram,(P(0.0, 0.0), P(5.0, 5.0), P(15.0, 5.0), P(10.0, 0.0))),
            Polygon(:rhombus,(P(0.0, 0.0), P(4., 3.), P(9., 3.), P(5., 0.))),
            Polygon(:multi,(P(0., 0.), P(0., 10.), P(4., 15.), P(10., 15.), P(15., 10.), P(15., 4.), P(10., 0.) )),
            Polygon(:line,(P(0.0,0.0),P(10.0,0.0))),
            Polygon(:dot,(P(0.0,0.0),))]

"""
    many_sides(P::Polygon)

Returns T if P has more than 6 sides.
"""
function many_sides(P::Polygon)
    return length(P.edges) > 6;
end

"""
    nonzero_area(P::Polygon)

Returns T if P has nonzero area.
"""
function nonzero_area(P::Polygon)
    return area(P)!=0
end # function

"""
    area(P::Polygon)

    Computes the area enclosed by polygon P.
    Taken from https://rosettacode.org/wiki/Shoelace_formula_for_polygonal_area#Kotlin
"""
function area(P::Polygon)
    n = length(P.edges)
    a = 0.0
    for i in range(1,stop=n-1)
        a+= P.edges[i].x * P.edges[i+1].y - P.edges[i+1].x * P.edges[i].y
    end
    return abs(a + P.edges[n].x * P.edges[1].y - P.edges[1].x * P.edges[n].y) / 2.0
end # function

"""
    equal_sides(P::Polygon)

EQUAL-SIDES is true if all sides of P have equal length
"""
function equal_sides(P::Polygon)
    n = length(P.edges)
    len = norm(P.edges[1]-P.edges[n])
    for i=1:n-1
        if norm(P.edges[i]-P.edges[i+1]) != len
            return false
        end
    end
    return true
end # function

predicates = [many_sides,nonzero_area,equal_sides]

reporting = true

agenda = PriorityQueue{Any,Float64}(Base.Order.Reverse);

Singleton = World(universe,predicates,reporting,agenda)

Object = Concept(Singleton, "Abstract Object", (x)->true, 50, [], [], missing, [])

Singleton.agenda[(find_examples_of, Object)] = 50.0;
Singleton.agenda[(make_specialization, Object)] = 25.0;


text:
include("nanop.jl")

using NCDatasets
using Formatting
using Statistics

Re  = 6371000.0 #  m
c_p  = 1004.0    #  J / kg K
g   = 9.81      #  m / s^2
Lw  = 2.26e6    #  J / kg

casenames = [
    "lowres_SSM_SOM",
    "lowres_SSM_NK",
    "lowres_SSM_SOM_noQflux",
]

linestyle = ["r-", "g-", "b-"]

Dataset("domain.lnd.fv4x5_gx3v7.091218.nc", "r") do ds
    global lat, lon, lev, ilev, rec

    lat = replace(ds["yc"][1,:], missing=>NaN)
    lon = replace(ds["xc"][:,1], missing=>NaN)

end

using PyPlot

fig, ax = plt[:subplots](1, 1, figsize=(12,8))
ax[:plot]([-90, 90], [0, 0], linewidth=2, color="#cccccc")

for i = 1:length(casenames)

    casename = casenames[i]

    Dataset("$casename.ocn.h.ma.fv4x5.0001-0020.nc", "r") do ds

        t_rng = 240-5*12+1:240

        println(ds.dim["time"])

        if match(r"SSM_NK", casename) != nothing
            T_ML = ds["T"][:, :, 1, t_rng]
        elseif match(r"SSM_SOM", casename) != nothing
            T_ML = ds["T_ML"][:, :, t_rng]
        else
            throw(ErrorException("ERROR"))
        end

        T_ML = replace(T_ML, missing=>NaN)
        T_ML = nanmean(T_ML, dims=(1,3,))[1, :, 1] .- 273.15

        ax[:plot](lat, T_ML, linestyle[i], linewidth=2,  label=casename)

    end

end

ax[:legend]()
ax[:set_xlabel]("Latitude [deg]")
ax[:set_ylabel]("Sea Surface Temperature [degC]")
ax[:set_xlim](-90, 90)
ax[:set_ylim](-5, 30)

fig[:savefig]("compare_models_SST.png", dpi=200)

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
module TOMLIO

using Pkg.TOML

function stack(io::IO,msg::Vector{UInt8})
    frontbytes = reinterpret(UInt8,Int16[length(msg)])
    item = UInt8[frontbytes...,msg...]
    write(io,item)
end

function unstack(io::IO)
    sizebytes = [read(io,UInt8),read(io,UInt8)]
    size = reinterpret(Int16,sizebytes)[1]

    msg = UInt8[]
    for i in 1:size
        push!(msg,read(io,UInt8))
    end
    return msg
end

import Base: GenericIOBuffer

function serialize(io::GenericIOBuffer,x::Any) #:Certificate{MemberID})
    dict = Dict(x)
    TOML.print(io, dict)
end

function deserialize(io::GenericIOBuffer,::Type{T}) where T
    str = String(take!(io))
    dict = TOML.parse(str)
    #return Certificate{MemberID}(dict)
    return T(dict)
end

function serialize(io::IO,x::Any)
    msg = IOBuffer()
    serialize(msg,x)
    stack(io,take!(msg))
end

function deserialize(io::IO,x::Type{T}) where T
    msg = unstack(io)
    return deserialize(IOBuffer(msg),x)
end

export serialize, deserialize

end


text:
"""
    evaluate

Evaluate whether the mission is success or failure.

# Notes
jld2: a loaded dictionary from `JLD2.load`.
"""
function evaluate(jld2::Dict, t1::Real, threshold::Real)
    @assert threshold > 0
    @unpack df, traj_des, t0, tf = jld2
    @assert t0 < t1 && t1 < tf
    df_filtered = filter(:time => t -> t >= t1, df)
    poss = df_filtered.sol |> Map(datum -> datum.plant.state.p) |> collect
    poss_des = df_filtered.time |> Map(traj_des) |> collect
    error_max = poss - poss_des |> Map(norm) |> collect |> maximum
    is_success = error_max <= threshold
end