text:
## Functions to complete the model structure
function petlion(cathode::Function;kwargs...)
    θ = Dict{Symbol,Float64}()
    funcs = _funcs_numerical()

    anode, system = cathode(θ,funcs)
    if anode isa Function && system isa Function
        anode(θ,funcs)
        θ, bounds, opts, N, numerics = system(θ, funcs, cathode, anode; kwargs...)

        return initialize_param(θ, bounds, opts, N, numerics)
    else
        error("Cathode function must return both anode and system functions.")
    end
end
function petlion(;cathode=cathode,anode=anode,system=system, # Input chemistry - can be modified
    kwargs... # keyword arguments for system
    )
    θ = Dict{Symbol,Float64}()
    funcs = _funcs_numerical()

    cathode(θ, funcs)
    anode(θ, funcs)
    θ, bounds, opts, N, numerics = system(θ, funcs; kwargs...)

    p = initialize_param(θ, bounds, opts, N, numerics)

    return p
end

function initialize_param(θ, bounds, opts, _N, numerics)

    ind, N_diff, N_alg, N_tot = state_indices(_N, numerics)
    N = discretizations_per_section(_N.p, _N.s, _N.n, _N.a, _N.z, _N.r_p, _N.r_n, N_diff, N_alg, N_tot)
    cache = build_cache(θ, ind, N, numerics, opts)

    check_errors_initial(θ, numerics, N)

    # while θ must contain Float64s, we need it to possibly
    # contain any type during the function generation process
    θ_any = Dict{Symbol,Any}()
    @inbounds for key in keys(θ)
        θ_any[key] = θ[key]
    end
    θ[:I1C] = calc_I1C(θ)

    ## temporary params with no functions
    _p = model_skeleton(θ_any, numerics, N, ind, opts, bounds, cache)

    funcs = load_functions(_p)

    ## Real params with functions
    p = model(
        θ,
        numerics,
        N,
        ind,
        opts,
        bounds,
        cache,
        funcs,
    )
    return p
end

function build_cache(θ, ind, N, numerics, opts)
    """
    Creates the cache struct for _petlion
    """
    outputs_tot = Symbol[]
    @inbounds for (name,_type) in zip(fieldnames(typeof(ind)), fieldtypes(typeof(ind)))
        if !(_type <: Tuple) push!(outputs_tot, name) end
    end
    outputs_tot = (outputs_tot...,)

    function variables_in_indices()
        """
        Defines all the possible outputs from simulate
        """

        var = Symbol[]
        tot = Int64[]

        for field in outputs_tot
            ind_var = getproperty(ind, field)
            if ind_var.var_type ∈ (:differential, :algebraic)
                push!(var, field)
                push!(tot, ind_var.start)
            end
        end
        (var[sortperm(tot)]...,)
    end

    function variable_labels()
        """
        Returns a vector of labels for all the states. Same length as Y and YP
        """
        labels = Symbol[]

        function add_label!(var)
            x = getproperty(ind, var)
            sections = x.sections
            if !isempty(x.sections)
                for section in sections
                    for i in 1:length(getproperty(x, section))
                        push!(labels, Symbol(var, :_, section, i))
                    end
                end
            else
                push!(labels, var)
            end
        end

        for var in variables_in_indices()
            x = getproperty(ind, var)
            x.var_type ∈ (:differential,:algebraic) ? add_label!(var) : nothing
        end

        @assert length(labels) == N.tot

        return labels
    end

    θ_tot =  Float64[]
    θ_keys =  Symbol[]

    vars = variables_in_indices()

    opts.var_keep = model_states_logic(opts.outputs)[1]

    Y0 = zeros(Float64, N.tot)
    YP0 = zeros(Float64, N.tot)
    res = zeros(Float64, N.alg)
    Y_alg = zeros(Float64, N.alg)

    id = [
        ones(Int64, N.diff)
        zeros(Int64,N.alg)
        ]

    save_start_dict = Dict{save_start_info,Vector{Float64}}()

    cache = cache_run(
        θ_tot,
        θ_keys,
        variable_labels(),
        vars,
        save_start_dict,
        Y0,
        YP0,
        res,
        Y_alg,
        id,
        )

    return cache
end

Base.@kwdef mutable struct state_sections{T} <: AbstractArray{T,1}
    tot::AbstractArray{T,1} = nothing
    a = nothing
    p = nothing
    s = nothing
    n = nothing
    z = nothing
    sections::Tuple = ()
    var_type::Symbol = :NA
end
Base.size(state::state_sections) = size(state.tot)
Base.IndexStyle(::Type{<:state_sections}) = IndexLinear()
Base.getindex(state::state_sections, i::Int64)= state.tot[i]
Base.setindex!(state::state_sections, v, i::Int64)= (state.tot[i] = v)

function retrieve_states(Y::AbstractArray, p::AbstractModel)
    """
    Creates a dictionary of variables based on an input array and the indices in p.ind
    """
    if length(Y) < p.N.tot error("Input vector must be ≥ the number of indices.") end

    states = Dict{Symbol,Any}()
    ind = p.ind
    vars_in_use = p.cache.vars

    vars = fieldnames(solution)[findall(fieldtypes(solution) .<: AbstractArray{<:Number})]

    sections = Symbol[]
    @inbounds for (field,_type) in zip(fieldnames(index_state), fieldtypes(index_state))
        if _type <: AbstractUnitRange && !(field ∈ (:start,:stop))
            push!(sections, field)
        end
    end
    sections = (sections...,)

    @views @inbounds function add_state!(var)
        ind_var = getproperty(ind, var)

        section_values = Any[]

        for section in sections
            section_name = Symbol(var, :_, section)
            ind_section = getproperty(ind_var, section)

            push!(section_values, section ∈ ind_var.sections ? Y[ind_section] : nothing)
        end

        state = state_sections(
            var ∈ vars_in_use ? Y[ind_var] : [],
            section_values...,
            ind_var.sections,
            ind_var.var_type,
        )

        states[var] = state

        return nothing
    end

    @inbounds for var in vars
        add_state!(var)
    end

    return states
end

function state_new(x, sections::Tuple, p;
        var_type=:NA,
        section_length = [getproperty(p.N, section) for section in sections],
    )

    x_vec = Any[]
    start = 0
    @inbounds for (i,section) in enumerate(sections)
        N = section_length[i]
        ind = (1:N) .+ start
        push!(x_vec, @views @inbounds x[ind])

        start += N
    end

    state = state_sections(;
        tot = x,
        NamedTuple{sections}((x_vec...,))...,
        sections = sections,
        var_type = var_type,
    )

    return state
end

function state_indices(N, numerics)
    """
    Creates the struct which has the indices of all the states.
    There can be more items in the struct than the simulation actually uses –
    this is because states are needed despite not being a function in the DAE,
    e.g., temperature is still used in isothermal simulations
    """

    N_diff = 0
    N_alg = 0

    state_vars = Symbol[]
    function add(var::Symbol, tot, vars::Tuple, var_type::Symbol=:NA;
        radial::Bool = false, replace = 0:0)
        """
        Adds a state to the indices dictionary
        """
        if tot === nothing
            return index_state()
        elseif tot isa Int
            tot = tot:tot
        end

        if var_type ∈ (:differential, :algebraic)
            tot = tot .+ (N_diff+N_alg)
            if     var_type == :differential
                N_diff += length(tot)
            elseif var_type == :algebraic
                N_alg  += length(tot)
            end
        end

        ind_a = 1:N.a
        ind_p = 1:(radial ? N.p*N.r_p : N.p)
        ind_s = 1:N.s
        ind_n = 1:(radial ? N.n*N.r_n : N.n)
        ind_z = 1:N.z

        sections = Symbol[]
        ind_start = 0
        :a ∈ vars ? (a = tot[ind_a .+ ind_start]; push!(sections, :a); ind_start += length(ind_a)) : a = replace
        :p ∈ vars ? (p = tot[ind_p .+ ind_start]; push!(sections, :p); ind_start += length(ind_p)) : p = replace
        :s ∈ vars ? (s = tot[ind_s .+ ind_start]; push!(sections, :s); ind_start += length(ind_s)) : s = replace
        :n ∈ vars ? (n = tot[ind_n .+ ind_start]; push!(sections, :n); ind_start += length(ind_n)) : n = replace
        :z ∈ vars ? (z = tot[ind_z .+ ind_start]; push!(sections, :z); ind_start += length(ind_z)) : z = replace

        start = tot[1]
        stop = tot[end]

        push!(state_vars, var)

        return index_state(start, stop, a, p, s, n, z, (sections...,), var_type)
    end

    c_e_tot     = 1:(N.p+N.s+N.n)
    c_s_avg_tot = numerics.solid_diffusion == :Fickian ? (1:N.p*N.r_p + N.n*N.r_n) : (1:(N.p+N.n))
    T_tot       = numerics.temperature ? (1:(N.p+N.s+N.n) + (N.a+N.z)) : nothing
    film_tot    = numerics.aging == :SEI ? (1:N.n) : nothing
    Q_tot       = numerics.solid_diffusion == :polynomial ? (1:(N.p+N.n)) : nothing
    j_tot       = 1:(N.p+N.n)
    j_s_tot     = numerics.aging ∈ (:SEI, :R_aging) ? (1:N.n) : nothing
    SOH_tot     = numerics.aging ∈ (:SEI, :R_aging) ? 1 : nothing
    Φ_e_tot     = 1:(N.p+N.s+N.n)
    Φ_s_tot     = 1:(N.p+N.n)
    I_tot       = 1

    c_e     = add(:c_e,     c_e_tot,     (:p, :s, :n),         :differential)
    c_s_avg = add(:c_s_avg, c_s_avg_tot, (:p, :n),             :differential; radial = numerics.solid_diffusion == :Fickian)
    T       = add(:T,       T_tot,       (:a, :p, :s, :n, :z), :differential)
    film    = add(:film,    film_tot,    (:n,),                :differential)
    Q       = add(:Q,       Q_tot,       (:p, :n),             :differential)
    SOH     = add(:SOH,     SOH_tot,     (),                   :differential)
    j       = add(:j,       j_tot,       (:p, :n),             :algebraic)
    j_s     = add(:j_s,     j_s_tot,     (:n,),                :algebraic)
    Φ_e     = add(:Φ_e,     Φ_e_tot,     (:p, :s, :n),         :algebraic)
    Φ_s     = add(:Φ_s,     Φ_s_tot,     (:p, :n),             :algebraic)
    I       = add(:I,       I_tot,       (),                   :algebraic)

    N_tot = N_diff + N_alg

    # These are the rest of the fields in the model_states struct that, while must be input, are unused
    Y = YP = t = V = P = SOC = index_state()
    state_vars = (state_vars...,)

    ind = indices_states(Y, YP, c_e, c_s_avg, T, film, Q, j, j_s, Φ_e, Φ_s, I, t, V, P, SOC, SOH, state_vars)

    return ind, N_diff, N_alg, N_tot
end

@inline function guess_init(p::AbstractModel, X_applied=0.0)
    """
    Get the initial guess in the DAE initialization.
    This function is made symbolic by Symbolics and saved as
    `initial_guess.jl`
    """

    Y0  = zeros(eltype(p.θ[:c_e₀]), p.N.tot)
    YP0 = zeros(eltype(p.θ[:c_e₀]), p.N.tot)

    states = retrieve_states(Y0, p)

    build_T!(states, p)

    states[:c_s_avg].p .= p.θ[:c_max_p] * (p.opts.SOC*(p.θ[:θ_max_p] - p.θ[:θ_min_p]) + p.θ[:θ_min_p])
    states[:c_s_avg].n .= p.θ[:c_max_n] * (p.opts.SOC*(p.θ[:θ_max_n] - p.θ[:θ_min_n]) + p.θ[:θ_min_n])

    build_c_s_star!(states, p)

    build_OCV!(states, p)

    # differential
    states[:c_e] = repeat([p.θ[:c_e₀]], (p.N.p+p.N.s+p.N.n))

    states[:T] = repeat([p.θ[:T₀]], (p.N.p+p.N.s+p.N.n)+p.N.a+p.N.z)

    states[:film] = zeros(p.N.n)

    states[:Q] = zeros(p.N.p+p.N.n)

    if p.numerics.aging ∈ (:SEI, :R_aging) states[:SOH] = 1.0 end

    # algebraic
    states[:j] = 0.0

    states[:Φ_e] = 0.0

    states[:Φ_s] = states[:U]

    states[:I] = X_applied

    if !isempty(states[:j_s]) states[:j_s] .= 0.0 end

    build_residuals!(Y0, states, p)

    return Y0, YP0
end

model_info(p::AbstractModel) = model_info(p.N, p.numerics)
function model_info(N::T1,numerics::T2) where {T1<:discretizations_per_section,T2<:options_numerical}
    version = "PETLION version: v"*join(Symbol.(PETLION_VERSION),".")

    numerical = ["$field: $(getproperty(numerics,field))" for field in fieldnames(T2)]

    discretization = [
        "N.p: $(N.p)";
        "N.s: $(N.s)";
        "N.n: $(N.n)";
        numerics.temperature                  ? "N.a: $(N.a)\nN.z: $(N.z)" : "";
        numerics.solid_diffusion == :Fickian ? "N.r_p: $(N.r_p)\nN.r_n: $(N.r_n)" : "";
        ]
    filter!(!isempty,discretization)

    str  = version * "\n"^2
    str *= "options_numerical\n" * join(numerical, "\n") * "\n"^2
    str *= "discretizations_per_section\n" * join(discretization, "\n")

    return str
end

function strings_directory_func(N::discretizations_per_section, numerics::T; create_dir=false) where T<:options_numerical

    dir_saved_models = "saved_models"
    # If the file directory is not specified, use the current working directory
    if !isnothing(options[:FILE_DIRECTORY])
        dir_saved_models = joinpath(options[:FILE_DIRECTORY], dir_saved_models)
    end

    if create_dir && !isdir(dir_saved_models)
        mkdir(dir_saved_models)
    end

    dir_cell = "$dir_saved_models/$(Symbol(numerics.cathode))_$(Symbol(numerics.anode))"

    str = join(
        [
            ["$(getproperty(numerics,field))" for field in filter(x->!(x ∈ (:cathode,:anode)), fieldnames(T))];
            "Np$(N.p)";
            "Ns$(N.s)";
            "Nn$(N.n)";
            numerics.temperature                  ? "Na$(N.a)_Nz$(N.z)" : "";
            numerics.solid_diffusion == :Fickian ? "Nr_p$(N.r_p)_Nr_n$(N.r_n)" : "";
        ],
        "_"
    )

    str = Base.bytes2hex(sha1(str))

    dir = "$dir_cell/$str"

    if create_dir
        for x in (dir_saved_models, dir_cell, dir)
            if !isdir(x)
                mkdir(x)
            end
        end
    end

    return dir
end

function strings_directory_func(p::AbstractModel; create_dir=false)
    strings_directory_func(p.N, p.numerics; create_dir=create_dir)
end

function strings_directory_func(p::AbstractModel, x; kw...)
    strings_directory = string("$(strings_directory_func(p; kw...))/$x.jl")

    return strings_directory
end

@inline function trapz(x::T1,y::T2) where {T1<:AbstractVector{<:Number},T2<:AbstractVector{<:Number}}
    """
    Trapezoidal rule with SIMD vectorization
    """
    @assert length(x) == length(y)
    out = 0.0
    @inbounds @simd for i in 2:length(x)
        out += 0.5*(x[i] - x[i-1])*(y[i] + y[i-1])
    end
    return out
end
@inline function cumtrapz(x::T1,y::T2) where {T1<:AbstractVector{<:Number},T2<:AbstractVector{<:Number}}
    # Check matching vector length
    @assert length(x) == length(y)

    # Initialize Output
    out = similar(x)
    out[1] = 0

    # Iterate over arrays
    @inbounds for i in 2:length(x)
        out[i] = out[i-1] + 0.5*(x[i] - x[i-1])*(y[i] + y[i-1])
    end
    # Return output
    return out
end

function extrap_x_0(x::AbstractVector,y::AbstractVector)
    @inbounds y[1] - ((y[3] - y[1] - (((x[2] - x[1])^-1)*(x[3] - x[1])*(y[2] - y[1])))*((x[3]^2 - (x[1]^2) - (((x[2] - x[1])^-1)*(x[2]^2 - (x[1]^2))*(x[3] - x[1])))^-1)*(x[1]^2)) - ((y[2] - y[1] - (((x[3]^2 - (x[1]^2) - (((x[2] - x[1])^-1)*(x[2]^2 - (x[1]^2))*(x[3] - x[1])))^-1)*(x[2]^2 - (x[1]^2))*(y[3] - y[1] - (((x[2] - x[1])^-1)*(x[3] - x[1])*(y[2] - y[1])))))*((x[2] - x[1])^-1)*x[1])
end
function extrapolate_section(y::AbstractVector,p::AbstractModel,section::Symbol)
    """
    Extrapolate to the edges of the FVM sections with a second-order polynomial
    """

    N = getproperty(p.N, section)

    x_range = [
        0
        collect(range(1/2N,1-1/2N,length=N))
        1
    ]

    x_interp = @views @inbounds x_range[2:4]

    y_range = [
        extrap_x_0(x_interp,y)
        y
        extrap_x_0(x_interp,(@inbounds @views y[end:-1:end-2]))
    ]

    x_range *= p.θ[Symbol(:l_,section)]

    return x_range, y_range
end


text:
using GeoDataFrames; const GDF = GeoDataFrames
using Test
using DataFrames
using GeoFormatTypes; const GFT = GeoFormatTypes
using Pkg.PlatformEngines
using Dates

# Use ArchGDAL datasets to test with
probe_platform_engines!()  # for download

const testdatadir = @__DIR__

REPO_URL = "https://github.com/yeesian/ArchGDALDatasets/blob/master/"

remotefiles = [
    ("ospy/data1/sites.dbf", "7df95edea06c46418287ae3430887f44f9116b29715783f7d1a11b2b931d6e7d"),
    ("ospy/data1/sites.prj", "81fb1a246728609a446b25b0df9ede41c3e7b6a133ce78f10edbd2647fc38ce1"),
    ("ospy/data1/sites.sbn", "198d9d695f3e7a0a0ac0ebfd6afbe044b78db3e685fffd241a32396e8b341ed3"),
    ("ospy/data1/sites.sbx", "49bbe1942b899d52cf1d1b01ea10bd481ec40bdc4c94ff866aece5e81f2261f6"),
    ("ospy/data1/sites.shp", "69af5a6184053f0b71f266dc54c944f1ec02013fb66dbb33412d8b1976d5ea2b"),
    ("ospy/data1/sites.shx", "1f3da459ccb151958743171e41e6a01810b2a007305d55666e01d680da7bbf08"),
]
@info "Downloading test files..."
for (f, sha) in remotefiles
    localfn = joinpath(testdatadir, basename(f))
    url = REPO_URL * f * "?raw=true"
    PlatformEngines.download_verify(url, sha, localfn; force=true)
end

@testset "GeoDataFrames.jl" begin
    fn = joinpath(testdatadir, "sites.shp")
    coords = zip(rand(10), rand(10))

    @testset "Read shapefile" begin
        t = GDF.read(fn)
        @test nrow(t) == 42
        @test "ID" in names(t)
    end

    @testset "Read shapefile with layer id" begin
        t = GDF.read(fn, 0)
        @test nrow(t) == 42
        @test "ID" in names(t)
    end

    @testset "Read shapefile with layer name" begin
        t = GDF.read(fn, "sites")
        @test nrow(t) == 42
        @test "ID" in names(t)
    end

    @testset "Read shapefile with non-existing layer name" begin
        @test_throws ArgumentError GDF.read(fn, "foo")
    end

    # @testset "Read shapefile with NULLs" begin
    #     fnn = joinpath(testdatadir, "null.gpkg")
    #     t = GDF.read(fnn)
    #     @test nrow(t) == 2
    #     @test "ID" in names(t)
    #     @test t.t == [1,2]
    # end

    @testset "Read self written file" begin
        # Save table with a few random points
        table = DataFrame(geom=createpoint.(coords), name="test")
        GDF.write("test_points.shp", table)
        GDF.write("test_points.gpkg", table, layer_name="test_points")
        GDF.write("test_points.geojson", table, layer_name="test_points")

        ntable = GDF.read("test_points.shp")
        @test nrow(ntable) == 10
        ntable = GDF.read("test_points.gpkg")
        @test nrow(ntable) == 10
        ntable = GDF.read("test_points.geojson")
        @test nrow(ntable) == 10
    end

    @testset "Write shapefile" begin

        t = GDF.read(fn)

        # Save table from reading
        GDF.write("test_read.shp", t, layer_name="test_coastline")
        GDF.write("test_read.gpkg", t, layer_name="test_coastline")
        GDF.write("test_read.geojson", t, layer_name="test_coastline")

    end

    @testset "Write shapefile with non-GDAL types" begin
        coords = zip(rand(Float32, 2), rand(Float32, 2))
        t = DataFrame(geom=createpoint.(coords), name=["test", "test2"], check=[false, true], z=[Float32(8), Float32(-1)], odd=[1, missing], date=[now(), now()])

        GDF.write("test_exotic.shp", t)
        GDF.write("test_exotic.gpkg", t)
        GDF.write("test_exotic.geojson", t)
    end

    @testset "Read shapefile with non-GDAL types" begin
        t = GDF.read("test_exotic.shp")
        GDF.read("test_exotic.gpkg")
        GDF.read("test_exotic.geojson")
    end

    @testset "Spatial operations" begin
        table = DataFrame(geom=createpoint.(coords), name="test")

        # Buffer to also write polygons
        table.geom = buffer(table.geom, 10)
        GDF.write("test_polygons.shp", table)
        GDF.write("test_polygons.gpkg", table)
        GDF.write("test_polygons.geojson", table)
    end

    @testset "Reproject" begin
        table = DataFrame(geom=createpoint.([[0,0,0]]), name="test")
        reproject(table.geom, GFT.EPSG(4326), GFT.EPSG(28992))
        @test GDF.AG.getpoint(table.geom[1], 0)[1] ≈ -587791.596556932
        GDF.write("test_reprojection.gpkg", table, crs=GFT.EPSG(28992))
    end
end

ext: jl


text:
import LogDensityProblems
import ProgressMeter
import TransformVariables

const default_learning_rate = 1e-4
const default_max_iterations = 1_000_000
const default_show_progress_meter = true
const default_throw_convergence_exception = true
const default_tolerance = 1e-10

"""
    fit(transformed_gradient_problem, theta_hat_initial; kwargs...)

Find the maximum likelihood estimatator for the parameters `theta`.

# Arguments
- `transformed_gradient_problem`
- `theta_hat_initial`

# Keyword Arguments
- `learning_rate`. Default value: $(default_learning_rate)
- `max_iterations`. Default value: $(default_max_iterations)
- `show_progress_meter`. Default value: $(default_show_progress_meter)
- `throw_convergence_exception`. Default value: $(default_throw_convergence_exception)
- `tolerance`. Default value: $(default_tolerance)

See the documentation for fully worked-out examples.
"""
function fit(transformed_gradient_problem,
             theta_hat_initial;
             learning_rate = default_learning_rate,
             max_iterations = default_max_iterations,
             show_progress_meter = default_show_progress_meter,
             throw_convergence_exception = default_throw_convergence_exception,
             tolerance = default_tolerance)
    transformed_log_density = parent(transformed_gradient_problem)
    transformation = transformed_log_density.transformation
    theta_hat_initial_deepcopy = deepcopy(theta_hat_initial)
    theta_hat_initial_deepcopy_inversetransformed = TransformVariables.inverse(transformation,
                                                                               theta_hat_initial_deepcopy)
    theta_hat_initial_deepcopy_inversetransformed_deepcopy = deepcopy(theta_hat_initial_deepcopy_inversetransformed)
    theta_hat_new_inversetransformed = _fit(transformed_gradient_problem,
                                            theta_hat_initial_deepcopy_inversetransformed_deepcopy;
                                            learning_rate = learning_rate,
                                            max_iterations = max_iterations,
                                            show_progress_meter = show_progress_meter,
                                            throw_convergence_exception = throw_convergence_exception,
                                            tolerance = tolerance)
    theta_hat_new = TransformVariables.transform(transformation,
                                                 theta_hat_new_inversetransformed)
    return theta_hat_new
end

function _fit(transformed_gradient_problem,
              theta_hat_initial;
              learning_rate,
              max_iterations,
              show_progress_meter,
              throw_convergence_exception,
              tolerance)
    theta_hat_new = deepcopy(theta_hat_initial)
    if show_progress_meter
        progress = ProgressMeter.ProgressThresh(tolerance, "Maximizing:")
    end
    for iter = 1:max_iterations
        theta_hat_old = theta_hat_new
        log_likelihood_value, gradient = LogDensityProblems.logdensity_and_gradient(transformed_gradient_problem,
                                                                                theta_hat_old)
        update = learning_rate * gradient
        theta_hat_new = theta_hat_old + update
        update_norm = sum(abs, update)
        if show_progress_meter
            showvalues = [(:iteration, iter),
                          (:loglikelihood, log_likelihood_value)]
            ProgressMeter.update!(progress,
                                  update_norm;
                                  showvalues = showvalues)
        end
        if update_norm < tolerance
            @info("Successfully converged after $(iter) iterations.")
            return theta_hat_new
        end
    end
    if show_progress_meter
        ProgressMeter.finish!(progress)
    end
    if throw_convergence_exception
        @error("Failed to converge after $(max_iterations) iterations.")
        throw(ConvergenceException())
    else
        @warn("Failed to converge after $(max_iterations) iterations.")
    end
    return theta_hat_new
end


text:
using TM1, JSON, HTTP, MbedTLS
using Dates, Test, Base64
using TM1: name

include("utils/auth_tests.jl")

include("dimensions/dimensions_tests.jl")
include("dimensions/hierarchies_tests.jl")
include("dimensions/elements_tests.jl")
include("dimensions/subsets_tests.jl")

include("cubes/cubes_tests.jl")

include("processes/processes_tests.jl")
include("processes/chores_tests.jl")

include("cellsets/cellsets_tests.jl")


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.

function random_aff_expr(N, vars::Vector)
    ex = Expr(:call, :+)
    for _ in 1:N
        vl, vr = Any[], Any[]
        for i in randperm(length(vars))
            v = vars[i]
            if rand(Bool)
                push!(vl, :($(rand()-0.5)))
            end
            if rand(Bool)
                push!(vr, :($(rand()-0.5)))
            end
            if rand(Bool)
                push!(vl, :($(rand()-0.5) * $v))
            end
            if rand(Bool)
                push!(vr, :($(rand()-0.5) * $v))
            end
            if rand(Bool)
                push!(vl, :($v * $(rand()-0.5)))
            end
            if rand(Bool)
                push!(vl, :($v * $(rand()-0.5)))
            end
        end
        if rand(Bool)
            push!(vl, :($(rand()-0.5)))
        end
        if rand(Bool)
            push!(vr, :($(rand()-0.5)))
        end
        if !isempty(vl) || !isempty(vr)
            if isempty(vl)
                tmp = Expr(:call, :+, vr...)
            elseif isempty(vr)
                tmp = Expr(:call, :+, vl...)
            else
                tmp = Expr(:call, :*,
                        Expr(:call, :+, vl...),
                        Expr(:call, :+, vr...))
            end
            push!(ex.args, tmp)
        end
    end
    return ex
end

m = Model()
@variable(m, x)
@variable(m, y)
@variable(m, z)
@variable(m, w)
@variable(m, v)

N = 5
vars = [:x, :y, :z, :w, :v, :(identity(x)), :(identity(y)), :(identity(z))]

const ε = 10eps()
nvars = length(vars)

function test_approx_equal_exprs(ex1, ex2)
    res = true
    # test constant term
    abs(ex1.aff.constant - ex2.aff.constant) < ε || (res = false)

    # test aff terms
    vals = zeros(nvars)
    for i in 1:length(ex1.aff.vars)
        vals[ex1.aff.vars[i].col] += ex1.aff.coeffs[i]
    end
    for i in 1:length(ex2.aff.vars)
        vals[ex2.aff.vars[i].col] -= ex2.aff.coeffs[i]
    end
    for v in vals
        abs(v) < ε || (res = false)
    end

    # test quad terms
    qvals = zeros(nvars,nvars)
    for i in 1:length(ex1.qcoeffs)
        j,k = sort([ex1.qvars1[i].col, ex1.qvars2[i].col])
        qvals[j,k] += ex1.qcoeffs[i]
    end
    for i in 1:length(ex2.qcoeffs)
        j,k = sort([ex2.qvars1[i].col, ex2.qvars2[i].col])
        qvals[j,k] -= ex2.qcoeffs[i]
    end
    for v in vals
        abs(v) < ε || (res = false)
    end
    if !res
        str = "The following expression did not pass the fuzzer:\n    ex1 = QuadExpr($(ex1.qvars1),$(ex1.qvars2),$(ex1.qcoeffs),AffExpr($(ex1.aff.vars),$(ex1.aff.coeffs),$(ex1.aff.constant)))\n    ex2 = QuadExpr($(ex2.qvars1),$(ex2.qvars2),$(ex2.qcoeffs),AffExpr($(ex2.aff.vars),$(ex2.aff.coeffs),$(ex2.aff.constant)))"
        warn(str)
    end
    return res
end

println("[fuzzer] Check macros for expression construction")

for _ in 1:100
    raff = random_aff_expr(N, vars)
    ex = @eval @JuMP.Expression($raff)
    @fact test_approx_equal_exprs(ex, eval(raff)) --> true
end


text:
using MLBase
using Base.Test

import StatsBase
import StatsBase: harmmean

## correctrate & errorrate

a = [1, 1, 1, 2, 2, 2, 3, 3]
b = [1, 1, 2, 2, 2, 3, 3, 3]

@test correctrate(a, b) == 0.75
@test errorrate(a, b) == 0.25

## confusmat

@test confusmat(3, a, b) == [2 1 0; 0 2 1; 0 0 2]

## counthits & hitrates

gt = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2]
rs = [1 2 2 1 3 2 1 1 3 3;
      2 3 1 3 1 1 2 2 1 2;
      4 4 3 2 2 3 3 3 4 1;
      3 1 4 4 4 4 4 4 2 4]

@test [counthits(gt, rs, k) for k=1:5] == [3, 8, 8, 10, 10]

@test counthits(gt, rs, 1:3) == [3, 8, 8]
@test counthits(gt, rs, [2, 4]) == [8, 10]
@test counthits(gt, rs, 1:2:5) == [3, 8, 10]

@test_approx_eq [hitrate(gt, rs, k) for k=1:5] [0.3, 0.8, 0.8, 1.0, 1.0]
@test_approx_eq hitrates(gt, rs, 1:3) [0.3, 0.8, 0.8]
@test_approx_eq hitrates(gt, rs, [2, 4]) [0.8, 1.0]
@test_approx_eq hitrates(gt, rs, 1:2:5) [0.3, 0.8, 1.0]

## ROCNums

r = ROCNums{Int}(
    100, # p == tp + fn
    200, # n == tn + fp
    80,  # tp
    150, # tn
    50,  # fp
    20)  # fn

@test true_positive(r) == 80
@test true_negative(r) == 150
@test false_positive(r) == 50
@test false_negative(r) == 20

@test true_positive_rate(r) == 0.80
@test true_negative_rate(r) == 0.75
@test false_positive_rate(r) == 0.25
@test false_negative_rate(r) == 0.20

@test recall(r) == 0.80
@test sensitivity(r) == 0.80
@test specificity(r) == 0.75
@test precision(r) == (8/13)
@test positive_predictive_value(r) == (8/13)
@test negative_predictive_value(r) == (15/17)
@test false_discovery_rate(r) == (5/13)
@test accuracy(r) == (23/30)
@test_approx_eq f1score(r) harmmean([recall(r), precision(r)])

## auxiliary: find_thresbin & lin_threshold

ts = [2, 4, 6, 8]
ft = [MLBase.find_thresbin(i, ts) for i = 1:9]
@test isa(ft, Vector{Int})
@test ft == [1, 2, 2, 3, 3, 4, 4, 5, 5]

@test MLBase.lin_thresholds([1, 5], 5, Forward) == 1.0:1.0:5.0
@test MLBase.lin_thresholds([1, 5], 5, Reverse) == 5.0:-1.0:1.0

## roc

gt = [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]
ss = [1., 2., 3., 4., 5., 4., 5., 6., 7., 8., 9.]

r2 = ROCNums{Int}(6, 5, 6, 1, 4, 0)
r4 = ROCNums{Int}(6, 5, 6, 3, 2, 0)
r6 = ROCNums{Int}(6, 5, 4, 5, 0, 2)

@test roc(gt, ss, 2.0) == r2
@test roc(gt, ss, 4.0) == r4
@test roc(gt, ss, 6.0) == r6
@test roc(gt, ss, 2.0, Forward) == r2
@test roc(gt, ss, 4.0, Forward) == r4
@test roc(gt, ss, 6.0, Forward) == r6

@test roc(gt, ss, [2.0, 4.0, 6.0], Forward) == [r2, r4, r6]
@test roc(gt, ss, [2.0, 4.0, 6.0]) == [r2, r4, r6]

@test roc(gt, ss, 5) == roc(gt, ss, [1., 3., 5., 7., 9.])
@test roc(gt, ss) == roc(gt, ss, 100)

gt = [1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0]
ss = [1., 2., 3., 4., 5., 4., 5., 6., 7., 8., 9.]

r2 = ROCNums{Int}(5, 6, 2, 6, 0, 3)
r4 = ROCNums{Int}(5, 6, 4, 5, 1, 1)
r6 = ROCNums{Int}(5, 6, 5, 3, 3, 0)

@test roc(gt, ss, 2.0, Reverse) == r2
@test roc(gt, ss, 4.0, Reverse) == r4
@test roc(gt, ss, 6.0, Reverse) == r6

@test roc(gt, ss, [6.0, 4.0, 2.0], Reverse) == [r6, r4, r2]

gt = [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2]
pr = [1, 1, 1, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1]
ss = [0.2, 0.2, 0.2, 0.3, 0.3, 0.6, 0.6, 0.6, 0.6, 0.6, 0.8, 0.8, 0.8, 0.8, 0.8]

r00 = roc(gt, (pr, ss), 0.00)
r25 = roc(gt, (pr, ss), 0.25)
r50 = roc(gt, (pr, ss), 0.50)
r75 = roc(gt, (pr, ss), 0.75)
r100 = roc(gt, (pr, ss), 1.00)

@test r00 == ROCNums{Int}(10, 5, 8, 0, 5, 0)
@test r25 == ROCNums{Int}(10, 5, 8, 3, 2, 0)
@test r50 == ROCNums{Int}(10, 5, 8, 5, 0, 0)
@test r75 == ROCNums{Int}(10, 5, 4, 5, 0, 5)
@test r100 == ROCNums{Int}(10, 5, 0, 5, 0, 10)

@test roc(gt, (pr, ss), 0.0:0.25:1.0) == [r00, r25, r50, r75, r100]
# @test roc(gt, (pr, ss), 7) == roc(gt, (pr, ss), 0.2:0.1:0.8, Forward)
@test roc(gt, (pr, ss)) == roc(gt, (pr, ss), MLBase.lin_thresholds([0.2, 0.8], 100, Forward))


text:
function printpath(d::AbstractDict, s::Int, v::Int; pred = 2)::Nothing
    if s == v
        println(s)
    elseif d[v][pred] == -1
        error("no path from $s to $v exists.")
    else
        printpath(d, s, d[v][pred])
        println(v)
    end
end


text:
export fields

add_verbose_scope(:Fields)
add_assert_scope(:Fields)

add_verbose_scope(:FieldsNonFancy)

################################################################################
#
#  Types
#
################################################################################

mutable struct cocycle_ctx
  projection::GAP.GapObj
  inclusion::GAP.GapObj
  cocycle::Function
  values_cyclic::Function
  gen_kernel::GAP.GapObj
  inclusion_of_pSylow::GAP.GapObj

  function cocycle_ctx(proj::GAP.GapObj, incl::GAP.GapObj, cocycle::Function)
    z = new()
    z.projection = proj
    z.inclusion = incl
    z.cocycle = cocycle
    return z
  end
end

mutable struct FieldsTower
  field::AnticNumberField
  generators_of_automorphisms::Vector{NfToNfMor}
  subfields::Vector{NfToNfMor}
  ramified_primes::Vector{fmpz}
  isabelian::Bool
  #Solvable embedding problems for the extension
  #They are here to improve the conductor computation
  isomorphism::Dict{NfToNfMor, GAP.GapObj}
  admissible_cocycles::Vector{cocycle_ctx}
  projections_for_conductors::Vector{GAP.GapObj}

  function FieldsTower(K::AnticNumberField, auts::Vector{NfToNfMor}, subfields::Vector{NfToNfMor})
    z = new()
    z.field = K
    z.generators_of_automorphisms = auts
    z.subfields = subfields
    z.isabelian = false
    return z
  end

end

function Base.show(io::IO, F::FieldsTower)
  print(io, "Field context for the number field defined by $(F.field.pol)")
  return nothing
end

include("./merge.jl")
include("./abelian_layer.jl")
include("./read_write.jl")
include("./conductors.jl")
include("./brauer.jl")
include("./chain.jl")
include("./maximal_abelian_subextension.jl")
include("./non_normal.jl")

Generic.degree(F::FieldsTower) = degree(F.field)
Hecke.maximal_order(F::FieldsTower) = maximal_order(F.field)
number_field(F::FieldsTower) = F.field

function ramified_primes(F::FieldsTower)
  if !isdefined(F, :ramified_primes)
    f = factor(discriminant(maximal_order(F.field)))
    F.ramified_primes = collect(keys(f.fac))
  end
  return F.ramified_primes
end

################################################################################
#
#  Creation of the context
#
################################################################################

function field_context(K::AnticNumberField)
  layers = Vector{NfToNfMor}[]
  autsK = automorphisms(K, copy = false)
  permGC = _from_autos_to_perm(autsK)
  G = _perm_to_gap_grp(permGC)
  D2 = Vector{Tuple{GAP.GapObj, NfToNfMor}}(undef, length(autsK))
  for i = 1:length(autsK)
    p =  _perm_to_gap_perm(permGC[i])
    D2[i] = (p, autsK[i])
  end
  @assert GAP.Globals.IsSolvable(G)
  L = GAP.Globals.DerivedSeries(G)
  embs = Vector{NfToNfMor}(undef, length(L)-1)
  F = K
  for i = length(L)-1:-1:2
    H = L[i]
    gensGAP = GAP.Globals.GeneratorsOfGroup(H)
    ggs = NfToNfMor[ x[2] for x in D2 if GAP.Globals.IN(x[1], gensGAP)]
    push!(layers, closure(ggs))
    F, mF = fixed_field(K, ggs)
    F, mS = simplify(F)
    mF = mS * mF
    embs[i] = mF
  end
  H = L[1]
  gensGAP = GAP.Globals.GeneratorsOfGroup(H)
  ggs = NfToNfMor[ x[2] for x in D2 if GAP.Globals.IN(x[1], gensGAP)]
  push!(layers, closure(ggs))
  auts = small_generating_set(layers[1])
  for i = 2:length(layers)
    auts_layers = NfToNfMor[x for x in layers[i] if !(x in layers[i-1])]
    append!(auts, small_generating_set(auts_layers))
  end
  KQ = rationals_as_number_field()[1]
  embs[1] = hom(KQ, F, one(F))
  return FieldsTower(K, reverse(auts), embs)
end

################################################################################
#
#  Assure has automorphisms
#
################################################################################

function assure_automorphisms(T::FieldsTower)
  assure_automorphisms(T.field, T.generators_of_automorphisms)
end

function assure_automorphisms(K::AnticNumberField, gens::Vector{NfToNfMor})
  if isautomorphisms_known(K)
    return nothing
  end
  auts = closure(gens, degree(K))
  _set_automorphisms_nf(K, auts)
  return nothing
end

###############################################################################
#
#  From automorphisms to permutation groups
#
###############################################################################

function Base.push!(G::AbstractAlgebra.Generic.geobucket{T}, p::T) where {T <: AbstractAlgebra.MPolyElem}
   R = parent(p)
   i = max(1, ndigits(length(p), base=4))
   l = length(G.buckets)
   if length(G.buckets) < i
     resize!(G.buckets, i)
     for j in (l + 1):i
       G.buckets[j] = zero(R)
     end
   end
   G.buckets[i] = addeq!(G.buckets[i], p)
   while i <= G.len
      if length(G.buckets[i]) >= 4^i
         G.buckets[i + 1] = addeq!(G.buckets[i + 1], G.buckets[i])
         G.buckets[i] = R()
         i += 1
      end
      break
   end
   if i == G.len + 1
      Base.push!(G.buckets, R())
      G.len += 1
   end
end

function permutation_group(G::Vector{Hecke.NfRelNSToNfRelNSMor_nf_elem})
  permutations = permutation_group1(G)
  return _perm_to_gap_grp(permutations)
end

function permutations(G::Array{Hecke.NfToNfMor, 1})
  K = domain(G[1])
  n = length(G)
  dK = degree(K)
  d = numerator(discriminant(K.pol))
  p = 11
  R = GF(p, cached = false)
  Rx, x = PolynomialRing(R, "x", cached = false)
  fmod = Rx(K.pol)
  while iszero(discriminant(fmod))
    p = next_prime(p)
    R = GF(p, cached = false)
    Rx, x = PolynomialRing(R, "x", cached = false)
    fmod = Rx(K.pol)
  end
  pols = gfp_poly[x]
  gpol = Rx(image_primitive_element(G[1]))
  if gpol != x
    push!(pols, gpol)
    gpol = compose_mod(gpol, pols[2], fmod)
    while gpol != x
      push!(pols, gpol)
      gpol = compose_mod(gpol, pols[2], fmod)
    end
  end
  order = length(pols)

  for i in 2:n
    pi = Rx(image_primitive_element(G[i]))
    if !(pi in pols)
      previous_order = order
      order = order + 1
      push!(pols, pi)
      for j in 2:previous_order
        order = order + 1
        push!(pols, compose_mod(pols[j], pi, fmod))
      end
      if order == dK
        break
      end
      rep_pos = previous_order + 1
      while rep_pos <= order
        for k in 1:i
          po = Rx(image_primitive_element(G[k]))
          att = compose_mod(pols[rep_pos], po, fmod)
          if !(att in pols)
            order = order + 1
            push!(pols, att)
            for j in 2:previous_order
              order = order + 1
              push!(pols, compose_mod(pols[j], att, fmod))
            end
            if order == dK
              break
            end
          end
          if order == dK
            break
          end
        end
        if order == dK
          break
        end
        rep_pos = rep_pos + previous_order
      end
    end
  end
  #Now, I have the images mod p
  Dcreation = Vector{Tuple{gfp_poly, Int}}(undef, length(pols))
  for i = 1:length(pols)
    Dcreation[i] = (pols[i], i)
  end
  perms = Array{Array{Int, 1},1}(undef, n)
  for i = 1:n
    perms[i] = Vector{Int}(undef, dK)
  end
  gen_pols = gfp_poly[Rx(image_primitive_element(s)) for s in G]
  D = Dict{gfp_poly, Int}(Dcreation)
  for s = 1:n
    for i = 1:length(pols)
      perms[s][i] = D[Hecke.compose_mod(gen_pols[s], pols[i], fmod)]
    end
  end
  return perms
end

function permutation_group(G::Array{Hecke.NfToNfMor, 1})
  return _perm_to_gap_grp(permutations(G))
end

function _from_autos_to_perm(G::Array{Hecke.NfToNfMor,1})

  K = domain(G[1])
  @assert degree(K) == length(G)
  n = length(G)
  #First, find a good prime
  p = 3
  R = GF(p, cached = false)
  Rx, x = PolynomialRing(R, "x", cached = false)
  fmod = Rx(K.pol)
  while iszero(discriminant(fmod))
    p = next_prime(p)
    R = GF(p, cached = false)
    Rx, x = PolynomialRing(R, "x", cached = false)
    fmod = Rx(K.pol)
  end
  pols = Vector{Tuple{gfp_poly, Int}}(undef, n)
  for i = 1:n
    pols[i] = (Rx(image_primitive_element(G[i])), i)
  end
  D = Dict{gfp_poly, Int}(pols)
  permutations = Array{Array{Int, 1},1}(undef, n)
  for s = 1:n
    perm = Array{Int, 1}(undef, n)
    for i = 1:n
      perm[i] = D[Hecke.compose_mod(pols[i][1], pols[s][1], fmod)]
    end
    permutations[s] = perm
  end
  return permutations

end

function _perm_to_gap_grp(perm::Array{Array{Int, 1},1})
  g = GAP.GapObj[]
  for x in perm
    z = _perm_to_gap_perm(x)
    push!(g, z)
  end
  g1 = GAP.julia_to_gap(g)
  return GAP.Globals.Group(g1)
end

function _perm_to_gap_perm(x::Array{Int, 1})
  x1 = GAP.julia_to_gap(x)
  z = GAP.Globals.PermList(x1)
  return z
end

function IdGroup(autos::Array{NfToNfMor, 1})
  G = permutation_group(autos)
  return GAP.Globals.IdGroup(G)
end

###############################################################################
#
#  Split Extensions
#
###############################################################################

function _split_extension(G::Array{Hecke.NfToNfMor, 1}, mats::Array{Hecke.GrpAbFinGenMap, 1})

  gtype = map(Int, domain(mats[1]).snf)
  G1 = permutation_group(G)
  gensG1 = GAP.Globals.GeneratorsOfGroup(G1)
  A = GAP.Globals.AbelianGroup(GAP.julia_to_gap(gtype))
  gens = GAP.Globals.GeneratorsOfGroup(A)
  auts = Array{GAP.GapObj, 1}(undef, length(mats))
  for i = 1:length(mats)
    images = Array{GAP.GapObj, 1}(undef, length(gtype))
    for j = 1:length(gtype)
      g = GAP.Globals.Identity(A)
      for k = 1:length(gtype)
        if !iszero(mats[i].map[j,k])
          g *= gens[k]^Int(mats[i].map[j,k])
        end
      end
      images[j] = g
    end
    auts[i] = GAP.Globals.GroupHomomorphismByImages(A, A, gens, GAP.julia_to_gap(images))
  end
  AutGrp = GAP.Globals.Group(GAP.julia_to_gap(auts))
  mp = GAP.Globals.GroupHomomorphismByImages(G1, AutGrp, gensG1, GAP.julia_to_gap(auts))
  return GAP.Globals.SplitExtension(G1, mp, A)

end

###############################################################################
#
#  Check group extension
#
###############################################################################

function check_group_extension(TargetGroup::GAP.GapObj, autos::Array{NfToNfMor, 1}, res_act::Array{GrpAbFinGenMap, 1})

  GS = domain(res_act[1])
  @assert issnf(GS)
  expo = Int(GS.snf[end])
  K = domain(autos[1])
  d = degree(K)
  com, uncom = ppio(expo, d)

  if com == 1
    # I only need to check the split extension, since the second cohomology group is
    # trivial, regardless of the action
    if length(res_act) == 1 && isprime(order(GS)) == 1 && isprime(degree(K)) && iscoprime(d, order(GS))
      #Just need to check if the action is non trivial
      return !isone(mod(res_act[1].map[1, 1], GS.snf[1]))
    end
    H = _split_extension(autos, res_act)
    if GAP.Globals.IdGroup(H) == TargetGroup
      return true
    else
      return false
    end
  end

  if uncom == 1
    #Need a cohomological check. Only useful in the prime power case.
    return true
  end

  # I check the split extension related to only uncom
  #Now, I have to check if the split extension is isomorphic to IdH
  Qn, mQn = quo(GS, uncom, false)
  S1, mS1 = snf(Qn)
  new_res_act = Array{GrpAbFinGenMap, 1}(undef, length(res_act))
  for i = 1:length(res_act)
    Mat = mS1.map*mQn.imap*res_act[i].map*mQn.map*mS1.imap
    Hecke.reduce_mod_snf!(Mat, S1.snf)
    new_res_act[i] = hom(S1, S1, Mat)
  end
  H = _split_extension(autos, new_res_act)
  if GAP.Globals.IdGroup(H) == TargetGroup
    return true
  else
    return false
  end

end

###############################################################################
#
#  Interface to find Fields
#
###############################################################################

function field_extensions(list::Vector{FieldsTower}, bound::fmpz, IsoE1::GAP.GapObj, l::Array{Int, 1}, only_real::Bool)

  grp_to_be_checked = Dict{Int, GAP.GapObj}()
  d = degree(list[1])
  n = lcm(l)
  com, uncom = ppio(n, d)
  fac = factor(n)
  for (p, v) in fac
    grp_to_be_checked[Int(p)] = _construct_grp(IsoE1, Int(p)^v)
  end
  if uncom != 1
    IsoCheck = _construct_grp(IsoE1, uncom)
  else
    IsoCheck = IsoE1
  end
  final_list = FieldsTower[]
  for (j, x) in enumerate(list)
    @vprint :Fields 1 "Field $(j)/$(length(list)): $(x.field.pol)"
    @vprint :FieldsNonFancy 1 "Field $(j)/$(length(list)): $(x.field.pol)\n"
    append!(final_list, field_extensions(x, bound, IsoCheck, l, only_real, grp_to_be_checked, IsoE1))
  end
  return final_list

end

function field_extensions(x::FieldsTower, bound::fmpz, IsoE1::GAP.GapObj, l::Array{Int, 1}, only_real::Bool, grp_to_be_checked::Dict{Int, GAP.GapObj}, IsoG::GAP.GapObj)

  list_cfields = _abelian_normal_extensions(x, l, bound, IsoE1, only_real, IsoG)
  if isempty(list_cfields)
    @vprint :Fields 1 "\e[1F$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())Number of new fields found: 0\n\n"
    @vprint :FieldsNonFancy 1 "Number of new fields found: 0\n\n"
    return Vector{FieldsTower}()
  end
  list = from_class_fields_to_fields(list_cfields, x.generators_of_automorphisms, grp_to_be_checked, IsoG)
  @vprint :Fields 1 "\e[1F$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())Computing maximal orders"
  @vprint :FieldsNonFancy 1 "Computing maximal orders\n"
  final_list = Vector{FieldsTower}(undef, length(list))
  for j = 1:length(list)
    fld, autos, embed = _relative_to_absolute(list[j][1], list[j][2])
    previous_fields = Array{NfToNfMor, 1}(undef, length(x.subfields)+1)
    for s = 1:length(x.subfields)
      previous_fields[s] = x.subfields[s]
    end
    previous_fields[end] = embed
    final_list[j] = FieldsTower(fld, autos, previous_fields)
  end

  @vprint :Fields 1 "$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())"
  @vprint :Fields 1 "Number of new fields found: $(length(final_list))\n\n"
  @vprint :FieldsNonFancy 1 "Number of new fields found: $(length(final_list))\n\n"
  return final_list

end

###############################################################################
#
#  Interface
#
###############################################################################

function fields(a::Int, b::Int, list::Vector{FieldsTower}, absolute_bound::fmpz; only_real::Bool = false)
  G = GAP.Globals.SmallGroup(a, b)
  L = GAP.Globals.DerivedSeries(G)
  lvl = _real_level(L)
  first = true
  for i = 2:length(L)-1
    G1 = GAP.Globals.FactorGroup(L[1], L[i])
    if first && GAP.Globals.Size(G1) != degree(list[1].field)
      continue
    end
    first = false
    E1 = GAP.Globals.FactorGroup(L[1], L[i+1])
    H1 = GAP.Globals.FactorGroup(L[i], L[i+1])
    l = GAP.gap_to_julia(Vector{Int64}, GAP.Globals.AbelianInvariants(H1))
    @vprint :Fields 1 "contructing abelian extensions with invariants $l \n"
    @vprint :FieldsNonFancy 1 "contructing abelian extensions with invariants $l \n"
    o = divexact(GAP.Globals.Size(G), GAP.Globals.Size(E1))
    bound = root(absolute_bound, o)
    IsoE1 = GAP.Globals.IdGroup(E1)
    @vprint :Fields 1 "Number of fields at the $i -th step: $(length(list)) \n"
    @vprint :FieldsNonFancy 1 "Number of fields at the $i -th step: $(length(list)) \n"
    lG = snf(abelian_group(l))[1]
    invariants = map(Int, lG.snf)
    onlyreal = (lvl > i || only_real)
    #First, I search for obstruction.
    @vprint :Fields 1 "Computing obstructions\n"
    @vprint :FieldsNonFancy 1 "Computing obstructions\n"
    #@vtime :Fields 1
    list = check_obstruction(list, L, i, invariants)
    @vprint :Fields 1 "Fields to check: $(length(list))\n\n"
    @vprint :FieldsNonFancy 1 "Fields to check: $(length(list))\n\n"
    if isempty(list)
      return FieldsTower[]
    end
    list = field_extensions(list, bound, IsoE1, invariants, onlyreal)
    @vprint :Fields 1 "Step $i completed\n"
    @vprint :FieldsNonFancy 1 "Step $i completed\n"
    if isempty(list)
      return FieldsTower[]
    end
  end
  if first
    error("The fields given can not be extended!")
  end
  return list
end

function fields_direct_product(g1, g2, red, redfirst, absolute_bound; only_real = false)
  b1 = root(absolute_bound, g2[1])
  b2 = root(absolute_bound, g1[1])
  @vprint :Fields 1 "The group is the product of $(g1) and $(g2)\n"
  l2 = fields(g2[1], g2[2], b2, only_real = only_real)
  if isempty(l2)
    return FieldsTower[]
  end
  if g1 == g2
    return _merge(l2, l2, absolute_bound, red, redfirst, g1, g2)
  end
  l1 = fields(g1[1], g1[2], b1, only_real = only_real)
  if isempty(l1)
    return FieldsTower[]
  end
  return _merge(l1, l2, absolute_bound, red, redfirst, g1, g2)
end

function fields(a::Int, b::Int, absolute_bound::fmpz; using_direct_product::Bool = true, only_real::Bool = false)
  if a == 1
    @assert b == 1
    K = rationals_as_number_field()[1]
    g = hom(K, K, K(1))
    return FieldsTower[FieldsTower(K, NfToNfMor[g], Array{NfToNfMor, 1}())]
  end
  G = GAP.Globals.SmallGroup(a, b)
  if using_direct_product
    g1, g2, red, redfirst = direct_product_decomposition(G, (a, b))
    if g2 != (1, 1)
      @vprint :Fields 1 "computing extensions with Galois group ($a, $b) and bound ~10^$(clog(absolute_bound, 10))\n"
      return fields_direct_product(g1, g2, red, redfirst, absolute_bound; only_real = only_real)
    end
  end
  L = GAP.Globals.DerivedSeries(G)
  G1 = GAP.Globals.FactorGroup(L[1], L[end-1])
  invariants = GAP.gap_to_julia(Vector{Int}, GAP.Globals.AbelianInvariants(L[end-1]))
  lG = snf(abelian_group(invariants))[1]
  invariants = map(Int, lG.snf)
  if GAP.Globals.IsAbelian(G)
    @vprint :Fields 1 "computing abelian extension of Q with invariants $(invariants) and bound ~10^$(clog(absolute_bound, 10))\n"
    @vprint :FieldsNonFancy 1 "Doing Group ($a, $b) with bound $absolute_bound\n"
    return abelian_extensionsQQ(invariants, absolute_bound, only_real)
  end

  lvl = _real_level(L)
  IdGroupGAP = GAP.Globals.IdGroup(G1)
  IdGroup = GAP.gap_to_julia(Vector{Int}, IdGroupGAP)
  bound = root(absolute_bound, prod(invariants))
  list = fields(IdGroup[1], IdGroup[2], bound; using_direct_product = using_direct_product, only_real = (only_real || lvl == length(L)-1))
  if isempty(list)
    return list
  end
  @vprint :Fields 1 "computing extensions with Galois group ($a, $b) and bound ~10^$(clog(absolute_bound, 10))\n"
  @vprint :Fields 1 "Abelian invariants of the relative extension: $(invariants)\n"
  @vprint :Fields 1 "Number of fields at this step: $(length(list)) \n"
  @vprint :FieldsNonFancy 1 "Number of fields at this step: $(length(list)) \n"

  @vprint :Fields 1 "Computing obstructions\n"
  @vprint :FieldsNonFancy 1 "Computing obstructions\n"
  #@vtime :Fields 1
  list = check_obstruction(list, L, length(L)-1, invariants)
  @vprint :Fields 1 "Fields to check: $(length(list))\n\n"
  @vprint :FieldsNonFancy 1 "Fields to check: $(length(list))\n\n"
  if isempty(list)
    return FieldsTower[]
  end
  Id = GAP.Globals.IdGroup(G)
  return field_extensions(list, absolute_bound, Id, invariants, only_real)
end


text:
struct gvNode <: GraphvizPoperties
    id::Int64
    name::String
    attributes::Properties
end

gvNode(id::Int64) = gvNode(id, String(id), Properties())

const gvNodes = Vector{gvNode}

struct gvEdge <: GraphvizPoperties
    from::Int
    to::Int
    attributes::Properties
end

gvEdge(from::Int, to::Int) = gvEdge(from, to, Properties())

const gvEdges = Vector{gvEdge}

struct gvSubGraph
    type::String
    plot_options::Properties
    graph_options::Properties
    node_options::Properties
    edge_options::Properties
    nodes::gvNodes
    edges::gvEdges
end

# empty outer constructor:
function gvSubGraph(type::String)
    type = type
    plot_options = Properties()
    graph_options = Properties()
    node_options = Properties()
    edge_options = Properties()
    node = gvNodes()
    edges = gvEdges()

    return gvSubGraph(type, plot_options, graph_options, node_options, edge_options, node, edges)
end

const gvSubGraphs = Vector{gvSubGraph}

# the main struct!
struct GraphvizAttributes
    plot_options::Properties
    graph_options::Properties
    node_options::Properties
    edge_options::Properties
    subgraphs::gvSubGraphs
    nodes::gvNodes
    edges::gvEdges
end

# empty outer constructor:
function GraphvizAttributes()
    plot_options = Properties()
    graph_options = Properties()
    node_options = Properties()
    edge_options = Properties()
    subgraphs = gvSubGraphs()
    node = gvNodes()
    edges = gvEdges()

    return GraphvizAttributes(plot_options, graph_options, node_options, edge_options, subgraphs, node, edges)
end

# derived from SimpleWeightedGraph
function GraphvizAttributes(graph::AbstractSimpleWeightedGraph; node_label::Bool=true, edge_label::Bool=false)
    directed = Graphs.is_directed(graph)
    n = nv(graph)
    large_graph = 200

    plot_options = [Property("weights", (edge_label) ? "true" : "false"),
        Property("largenet", "200")]

    graph_options = [Property("center", "\"1,1\""),
        Property("overlap", "scale"),
        Property("concentrate", "true"),
        Property("layout", (directed) ? "dot" : "neato"),
        Property("size", (n < 20) ? "3.0" : ((n < 100) ? "7.0" : "10.0"))]  # scaling function missing!
    node_options = [Property("color", "Turquoise"),
        Property("fontsize", (node_label) ? ((n < 100) ? "7.0" : "5.0") : "1.0"),
        Property("width", (node_label) ? "0.25" : "0.20"),
        Property("height", (node_label) ? "0.25" : "0.20"),
        Property("fixedsize", "true"),
        Property("shape", (node_label) ? "circle" : "point")]
    edge_options = [Property("arrowsize", "0.5"),
        Property("arrowtype", "normal"),
        Property("fontsize", (edge_label) ? "8.0" : "1.0")]

    nodes = [gvNode(i, "$i", Properties()) for i = 1:n]

    edges = []
    for i = 1:n, j = 1:n
        if !(graph.weights[i, j] == 0)
            push!(edges, gvEdge(i, j, [Property("xlabel", graph.weights[i, j])]))
        end
    end

    gv_attr = GraphvizAttributes(plot_options, graph_options, node_options, edge_options, gvSubGraphs(), nodes, edges)

    # modifier by large networks:
    if n > large_graph
        gv_attr = mod_attr_large_network!(gv_attr)
    end

    return gv_attr
end

function mod_attr_large_network!(attrs::GraphvizAttributes)
    set!(attrs.graph_options, "fontsize", "1")
    set!(attrs.graph_options, "concetrate", "true")
    set!(attrs.graph_options, "layout", "sfdp")
    set!(attrs.plot_options, "weights", "false")
    set!(attrs.node_options, "shape", "point")
    set!(attrs.node_options, "color", "black")

    return attrs
end

#### SETTER:

"""
    set!(edges::gvEdges, from::Int, to::Int, attribute::Property)
    set!(edges::gvNodes, id::String, attribute::Property))
    set!(nodes::gvNodes, id::Int, attribute::Property)
    set!(nodes::gvNodes, name::String, attribute::Property)

    set (Graphviz) Property of nodes or edges

"""
function set!(edges::gvEdges, from::Int, to::Int, attribute::Property; override=true)
    if !isempty(edges)
        for e in edges
            if (e.from == from) && (e.to == to) && (override == true)
                set!(e.attributes, attribute.key, attribute.value)
            end
        end
    else
        push!(edges, gvEdge(from, to, [attribute]))
    end
end

function set!(nodes::gvNodes, id::Int, attribute::Property)
    if !isempty(nodes)
        for n in nodes
            if n.id == id
                set!(n.attributes, attribute.key, attribute.value)
            end
        end
    else
        push!(nodes, gvNode(id, String(id), [attribute]))
    end
end

function set!(nodes::gvNodes, name::String, attribute::Property)
    if !isempty(nodes)
        for n in nodes
            if n.name == name
                set!(n.attributes, attribute.key, attribute.value)
            end
        end
    else
        max = _max_id(nodes::gvNodes)
        push!(nodes, gvNode((max + 1), name, [attribute]))
    end
end

##### REMOVER:

function rm!(nodes::gvNodes, key::String, attribute::Property)
    if !isempty(nodes)
        for n in nodes
            if n.name == key
                rm!(n.attributes, attribute.key)
            end
        end
    end
end

function rm!(nodes::gvNodes, id::Int, attribute::Property)
    if !isempty(nodes)
        for n in nodes
            if n.id == id
                rm!(n.attributes, attribute.key)
            end
        end
    end
end

function rm!(edges::gvEdges, from::Int, to::Int, attribute::Property)
    if !isempty(edges)
        for e in edges
            if (e.from == from) && (e.to == to)
                rm!(e.attributes, attribute.key)
            end
        end
    end
end

##### GETTER:
function val(edges::gvEdges, from::Int, to::Int, key::String)
    if !isempty(edges)
        for e in edges
            if (e.from == from) && (e.to == to)
                return val(e.attributes, key)
            end
        end
    end
    return []
end

function val(nodes, id::Int64, key::String)
    if !isempty(nodes)
        for n in nodes
            if n.id == id
                return val(n.attributes, key)
            end
        end
    end
    return []
end

#### special functions:

function _get_label(nodes::gvNodes, id::Int64)
    if !isempty(nodes)
        for n in nodes
            if n.id == id
                label = val(n.attributes, "label")
                if !isempty(label)
                    return label
                end
            end
        end
    end
    return []
end

function _max_id(nodes::gvNodes)
    max = 0
    if !isempty(nodes)
        for n in nodes
            (n.id > max) ? max = n.id : nothing
        end
    end
    return max
end

function get_id(nodes, str)
    if !isempty(nodes)
        for n in nodes
            #  labl = val_node(nodes, n.id, "label")
            if n.name == str
                return n.id
            end
        end
    end
    return 0
end

function get_name(nodes, id)
    if !isempty(nodes)
        for n in nodes
            #  labl = val_node(nodes, n.id, "label")
            if n.id == id
                return n.name
            end
        end
    end
    return 0
end

function get_node(nodes::gvNodes, id::Int)
    for node in nodes
        if node.id == id
            return node
        end
    end
    return []
end

max_stars_repo_name: tragisch/PlotGraphviz

max_issues_repo_name: tragisch/PlotGraphviz


text:
"""
    ispangram(input)

Return `true` if `input` contains every alphabetic character (case insensitive).

"""
function clean_input(input)
    # Remove unnecessary characters
    return replace(lowercase(input), r"[0-9 _.\"]" => s"")
end

function ispangram(input)
    input_chars = Set(collect(clean_input(input)))
    all_alphabets = Set('a':'z')
    return all_alphabets == input_chars
end

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism


text:
using CSV
using Dates
using DataFrames
include("../src/TrueSkill.jl")
using .TrueSkill
global const ttt = TrueSkill
using DataFrames

base = Dates.value(Date("1970-01-01") - Date("1900-01-01"))

data = CSV.read("input/history.csv", DataFrame)
djokovic = "d643"
federer = "f324"
sampras = "s402"
lendl = "l018"
connors = "c044"
nadal = "n409"
john_mcenroe = "m047"
bjorn_borg = "b058"
aggasi = "a092"
hewitt = "h432"
edberg = "e004"
vilas = "v028"
nastase = "n008"

courier = "c243"
kuerten = "k293"
murray = "mc10"
wilander = "w023"
roddick = "r485"
#data.w1_id[occursin.("roddick",data.w1_name)]
#data.w1_name[occursin.("pi74",data.w1_id)]

days = Dates.value.(data[:,"time_start"] .- Date("1900-01-01"))
#dates = Date.(Dates.UTD.(days .+ Dates.value(Date("1900-01-01"))))

composition = [ r.double == "t" ? [[r.w1_id,r.w2_id],[r.l1_id,r.l2_id]] : [[r.w1_id],[r.l1_id]] for r in eachrow(data) ]
function fit()
    h = ttt.History(composition=composition, times = days, sigma = 1.6, gamma = 0.036)
    ttt.convergence(h,epsilon=0.01, iterations=10)
    return h
end
h = ttt.History(composition=composition, times = days, sigma = 1.6, gamma = 0.036)
lc = ttt.learning_curves(h)

dict = Dict{String,Vector}()
dict["federer_mu"] = [tp[2].mu for tp in lc["f324"]]
dict["federer_sigma"]= [tp[2].sigma for tp in lc["f324"]]
dict["federer_time"] = [tp[1] for tp in lc["f324"]]
dict["nadal_mu"] = [tp[2].mu for tp in lc["n409"]]
dict["nadal_sigma"] = [tp[2].sigma for tp in lc["n409"]]
dict["nadal_time"] = [tp[1] for tp in lc["n409"]]
dict["djokovic_mu"] = [tp[2].mu for tp in lc["d643"]]
dict["djokovic_sigma"] = [tp[2].sigma for tp in lc["d643"]]
dict["djokovic_time"] = [tp[1] for tp in lc["d643"]]
dict["sampras_mu"] = [tp[2].mu for tp in lc["s402"]]
dict["sampras_sigma"] = [tp[2].sigma for tp in lc["s402"]]
dict["sampras_time"] = [tp[1] for tp in lc["s402"]]
dict["aggasi_mu"] = [tp[2].mu for tp in lc["a092"]]
dict["aggasi_sigma"] = [tp[2].sigma for tp in lc["a092"]]
dict["aggasi_time"] = [tp[1] for tp in lc["a092"]]
dict["vilas_mu"] = [tp[2].mu for tp in lc["v028"]]
dict["vilas_sigma"] = [tp[2].sigma for tp in lc["v028"]]
dict["vilas_time"] = [tp[1] for tp in lc["v028"]]
dict["borg_mu"] = [tp[2].mu for tp in lc["b058"]]
dict["borg_sigma"] = [tp[2].sigma for tp in lc["b058"]]
dict["borg_time"] = [tp[1] for tp in lc["b058"]]
dict["mcenroe_mu"] = [tp[2].mu for tp in lc["m047"]]
dict["mcenroe_sigma"] = [tp[2].sigma for tp in lc["m047"]]
dict["mcenroe_time"] = [tp[1] for tp in lc["m047"]]

dict["lendl_mu"] = [tp[2].mu for tp in lc["l018"]]
dict["lendl_sigma"] = [tp[2].sigma for tp in lc["l018"]]
dict["lendl_time"] = [tp[1] for tp in lc["l018"]]

dict["connors_mu"] = [tp[2].mu for tp in lc["c044"]]
dict["connors_sigma"] = [tp[2].sigma for tp in lc["c044"]]
dict["connors_time"] = [tp[1] for tp in lc["c044"]]

dict["hewitt_mu"] = [tp[2].mu for tp in lc["h432"]]
dict["hewitt_sigma"] = [tp[2].sigma for tp in lc["h432"]]
dict["hewitt_time"] = [tp[1] for tp in lc["h432"]]

dict["edberg_mu"] = [tp[2].mu for tp in lc["e004"]]
dict["edberg_sigma"] = [tp[2].sigma for tp in lc["e004"]]
dict["edberg_time"] = [tp[1] for tp in lc["e004"]]

dict["nastase_mu"] = [tp[2].mu for tp in lc["n008"]]
dict["nastase_sigma"] = [tp[2].sigma for tp in lc["n008"]]
dict["nastase_time"] = [tp[1] for tp in lc["n008"]]

dict["courier_mu"] = [tp[2].mu for tp in lc["c243"]]
dict["courier_sigma"] = [tp[2].sigma for tp in lc["c243"]]
dict["courier_time"] = [tp[1] for tp in lc["c243"]]

dict["kuerten_mu"] = [tp[2].mu for tp in lc["k293"]]
dict["kuerten_sigma"] = [tp[2].sigma for tp in lc["k293"]]
dict["kuerten_time"] = [tp[1] for tp in lc["k293"]]

dict["murray_mu"] = [tp[2].mu for tp in lc["mc10"]]
dict["murray_sigma"] = [tp[2].sigma for tp in lc["mc10"]]
dict["murray_time"] = [tp[1] for tp in lc["mc10"]]

dict["wilander_mu"] = [tp[2].mu for tp in lc["w023"]]
dict["wilander_sigma"] = [tp[2].sigma for tp in lc["w023"]]
dict["wilander_time"] = [tp[1] for tp in lc["w023"]]

dict["roddick_mu"] = [tp[2].mu for tp in lc["r485"]]
dict["roddick_sigma"] = [tp[2].sigma for tp in lc["r485"]]
dict["roddick_time"] = [tp[1] for tp in lc["r485"]]

maxlen = maximum([length(value) for (key, value) in dict])

df = DataFrame(Dict(key => [value;repeat([missing],maxlen-length(value))] for (key, value) in dict))

CSV.write("output/atp_learning_curves_trueskill.csv", df; header=true)

ttt.convergence(h,epsilon=0.01, iterations=10)
#@time h = fit()
lc = ttt.learning_curves(h)

top = Set([key  for (key, value) in lc for s in value if ((s[2].mu > 6.0) & (s[1]>base))])

dict = Dict{String,Vector}()
dict["federer_mu"] = [tp[2].mu for tp in lc["f324"]]
dict["federer_sigma"]= [tp[2].sigma for tp in lc["f324"]]
dict["federer_time"] = [tp[1] for tp in lc["f324"]]
dict["nadal_mu"] = [tp[2].mu for tp in lc["n409"]]
dict["nadal_sigma"] = [tp[2].sigma for tp in lc["n409"]]
dict["nadal_time"] = [tp[1] for tp in lc["n409"]]
dict["djokovic_mu"] = [tp[2].mu for tp in lc["d643"]]
dict["djokovic_sigma"] = [tp[2].sigma for tp in lc["d643"]]
dict["djokovic_time"] = [tp[1] for tp in lc["d643"]]
dict["sampras_mu"] = [tp[2].mu for tp in lc["s402"]]
dict["sampras_sigma"] = [tp[2].sigma for tp in lc["s402"]]
dict["sampras_time"] = [tp[1] for tp in lc["s402"]]
dict["aggasi_mu"] = [tp[2].mu for tp in lc["a092"]]
dict["aggasi_sigma"] = [tp[2].sigma for tp in lc["a092"]]
dict["aggasi_time"] = [tp[1] for tp in lc["a092"]]
dict["vilas_mu"] = [tp[2].mu for tp in lc["v028"]]
dict["vilas_sigma"] = [tp[2].sigma for tp in lc["v028"]]
dict["vilas_time"] = [tp[1] for tp in lc["v028"]]
dict["borg_mu"] = [tp[2].mu for tp in lc["b058"]]
dict["borg_sigma"] = [tp[2].sigma for tp in lc["b058"]]
dict["borg_time"] = [tp[1] for tp in lc["b058"]]
dict["mcenroe_mu"] = [tp[2].mu for tp in lc["m047"]]
dict["mcenroe_sigma"] = [tp[2].sigma for tp in lc["m047"]]
dict["mcenroe_time"] = [tp[1] for tp in lc["m047"]]

dict["lendl_mu"] = [tp[2].mu for tp in lc["l018"]]
dict["lendl_sigma"] = [tp[2].sigma for tp in lc["l018"]]
dict["lendl_time"] = [tp[1] for tp in lc["l018"]]

dict["connors_mu"] = [tp[2].mu for tp in lc["c044"]]
dict["connors_sigma"] = [tp[2].sigma for tp in lc["c044"]]
dict["connors_time"] = [tp[1] for tp in lc["c044"]]

dict["hewitt_mu"] = [tp[2].mu for tp in lc["h432"]]
dict["hewitt_sigma"] = [tp[2].sigma for tp in lc["h432"]]
dict["hewitt_time"] = [tp[1] for tp in lc["h432"]]

dict["edberg_mu"] = [tp[2].mu for tp in lc["e004"]]
dict["edberg_sigma"] = [tp[2].sigma for tp in lc["e004"]]
dict["edberg_time"] = [tp[1] for tp in lc["e004"]]

dict["nastase_mu"] = [tp[2].mu for tp in lc["n008"]]
dict["nastase_sigma"] = [tp[2].sigma for tp in lc["n008"]]
dict["nastase_time"] = [tp[1] for tp in lc["n008"]]

dict["courier_mu"] = [tp[2].mu for tp in lc["c243"]]
dict["courier_sigma"] = [tp[2].sigma for tp in lc["c243"]]
dict["courier_time"] = [tp[1] for tp in lc["c243"]]

dict["kuerten_mu"] = [tp[2].mu for tp in lc["k293"]]
dict["kuerten_sigma"] = [tp[2].sigma for tp in lc["k293"]]
dict["kuerten_time"] = [tp[1] for tp in lc["k293"]]

dict["murray_mu"] = [tp[2].mu for tp in lc["mc10"]]
dict["murray_sigma"] = [tp[2].sigma for tp in lc["mc10"]]
dict["murray_time"] = [tp[1] for tp in lc["mc10"]]

dict["wilander_mu"] = [tp[2].mu for tp in lc["w023"]]
dict["wilander_sigma"] = [tp[2].sigma for tp in lc["w023"]]
dict["wilander_time"] = [tp[1] for tp in lc["w023"]]

dict["roddick_mu"] = [tp[2].mu for tp in lc["r485"]]
dict["roddick_sigma"] = [tp[2].sigma for tp in lc["r485"]]
dict["roddick_time"] = [tp[1] for tp in lc["r485"]]

maxlen = maximum([length(value) for (key, value) in dict])

using DataFrames
df = DataFrame(Dict(key => [value;repeat([missing],maxlen-length(value))] for (key, value) in dict))

CSV.write("output/atp_learning_curves.csv", df; header=true)


text:
# VERY SIMILAR TO 20190319_broadcastbenchmarks.jl
# EXTENDED FOR BEST CASE VS FUNCTION WRAPPERS

using BenchmarkTools

# see this section of docs https://docs.julialang.org/en/v1/base/arrays/#Broadcast-and-vectorization-1

function G_single(t::T, params::Vector{T}) where T<:Real

   G₀, G₁, η₁, G₂, η₂ = params

   G₀ + G₁*exp(-G₁/η₁*t) + G₂*exp(-G₂/η₂*t)

end

function G_vec(t::Vector{T}, params::Vector{T}) where T<:Real

   G₀, G₁, η₁, G₂, η₂ = params

   G₀ .+ G₁ .* exp.(-G₁./η₁.*t) .+ G₂.*exp.(-G₂./η₂.*t)

end

function G_loop(t::Vector{T}, params::Vector{T}) where T<:Real

    out = Array{T}(undef, length(t))

    @inbounds for (i,ti) in enumerate(t)
        out[i] = G_single(ti, params)
    end

    out
end

t = collect(0.0:0.1:10000)
params = [1.0, 1.0, 1.0, 1.0, 1.0]

@btime G_single.($t, $(params,)) 				# 1.971 ms (2 allocations: 781.39 KiB)
@btime broadcast($G_single, $t, $(params,)) 	# 1.869 ms (2 allocations: 781.39 KiB)
@btime [G_single(ti, $params) for ti in $t]		# 1.887 ms (2 allocations: 781.39 KiB)
@btime G_vec($t, $params) 						# 1.587 ms (2 allocations: 781.39 KiB)
@btime G_loop($t, $params) 						# 1.880 ms (2 allocations: 781.39 KiB)

nothing

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
function prime_factors(n)
    found = false

    v = 2

    i = 2

    while i < n - 1 && !found
        if mod(n, i) == 0
            found = true

            v = i
        end

        i += 1
    end

    if found
        factors = vcat([v], prime_factors(div(n, v)))
    else
        factors = [n]
    end

    factors
end

function remove_first(item, collection)
    index = findfirst(x -> x == item, collection)

    if index != nothing
        deleteat!(collection, index)
    else
        collection
    end
end

function overlap(a, b)
    a_copy = a

    for n in b
        a_copy = remove_first(n, a_copy)
    end

    sort(vcat(a_copy, b))
end

println(prod(reduce(overlap, map(i -> prime_factors(i), 1:20); init=[1])))

max_stars_repo_name: seaneshbaugh/rosetta-euler

max_issues_repo_name: seaneshbaugh/rosetta-euler

max_forks_repo_name: seaneshbaugh/rosetta-euler


text:
using Flux, MLDataPattern, Mill, JsonGrinder, JSON, Statistics, IterTools
using Serialization
using JsonGrinder: suggestextractor, ExtractDict
using Mill: reflectinmodel

###############################################################
# start by loading all samples
###############################################################
samples = map(readlines("/Users/tomas.pevny/Work/Presentations/JuliaMeetup/dataset/train.json")) do s
	JSON.parse(s)
end;
JSON.print(samples[3],2)

labelkey = "device_class"
minibatchsize = 100
iterations = 10_000
neurons = 20 		# neurons per layer

targets = map(i -> i[labelkey], samples)
foreach(i -> delete!(i, labelkey), samples)
foreach(i -> delete!(i, "id"), samples)

#####
#  Create the schema and extractor
#####
sch = JsonGrinder.schema(samples)
extractor = suggestextractor(sch)

#####
#  Convert samples to Mill structure and extract targets
#####
data = tmap(extractor, samples)
labelnames = unique(targets)

#####
#  Create the model
#####
model = reflectinmodel(sch, extractor,
	k -> Dense(k, neurons, relu),
	d -> SegmentedMeanMax(d),
	b = Dict("" => k -> Dense(k, length(labelnames))),
)

#####
#  Train the model
#####
function minibatch()
	idx = sample(1:length(data), minibatchsize, replace = false)
	reduce(catobs, data[idx]), Flux.onehotbatch(targets[idx], labelnames)
end

accuracy(x,y) = mean(map(xy -> labelnames[argmax(model(xy[1]).data[:])] == xy[2], zip(x, y)))

cb = () -> println("accuracy = ", accuracy(data, targets))
ps = Flux.params(model)
loss = (x,y) -> Flux.logitcrossentropy(model(x).data, y)
Flux.Optimise.train!(loss, ps, repeatedly(minibatch, iterations), ADAM(), cb = Flux.throttle(cb, 2))

###############################################################
#  Classify test data
###############################################################
test_samples = map(readlines("/Users/tomas.pevny/Work/Presentations/JuliaMeetup/dataset/test.json")) do s
	extractor(JSON.parse(s))
end
o = Flux.onecold(model(reduce(catobs, test_samples)).data);
ns = extract_target[:device_class].keyvalemap
ns = Dict([ v => k for (k,v) in ns]...)
o = [ns[i] for i in o]
# Id,Predicted
# crossentropy = 0.028738448 accuracy = 0.999067454149829. = F-score 0.86994
# writedlm("/Users/tomas.pevny/Work/Presentations/JuliaMeetup/dataset/kaggle.csv",hcat(id, o), ",")


text:
#=
export AbstractSymmetry

"""
    AbstractSymmetry{OperationType<:AbstractSpaceSymmetryOperation}

Abstract symmetry whose elements are of type `OperationType`.
"""
abstract type AbstractSymmetry{OperationType<:AbstractSpaceSymmetryOperation} end
=#


text:
module CheckedArithmetic

using CheckedArithmeticCore
import CheckedArithmeticCore: safearg_type, safearg, safeconvert, accumulatortype, acc

using Base.Meta: isexpr
using LinearAlgebra: Factorization, UniformScaling
using Random: AbstractRNG
using Dates

export @checked, @check
export accumulatortype, acc # re-export

const op_checked = Dict(
    Symbol("unary-") => :(Base.Checked.checked_neg),
    :abs => :(Base.Checked.checked_abs),
    :+   => :(Base.Checked.checked_add),
    :-   => :(Base.Checked.checked_sub),
    :*   => :(Base.Checked.checked_mul),
    :÷   => :(Base.Checked.checked_div),
    :div => :(Base.Checked.checked_div),
    :%   => :(Base.Checked.checked_rem),
    :rem => :(Base.Checked.checked_rem),
    :fld => :(Base.Checked.checked_fld),
    :mod => :(Base.Checked.checked_mod),
    :cld => :(Base.Checked.checked_cld),
    )

function replace_op!(expr::Expr, op_map::Dict)
    if expr.head == :call
        f, len = expr.args[1], length(expr.args)
        op = isexpr(f, :.) ? f.args[2].value : f # handle module-scoped functions
        if op === :+ && len == 2                 # unary +
            # no action required
        elseif op === :- && len == 2             # unary -
            op = get(op_map, Symbol("unary-"), op)
            if isexpr(f, :.)
                f.args[2].value = op
                expr.args[1] = f
            else
                expr.args[1] = op
            end
        else                                     # arbitrary call
            op = get(op_map, op, op)
            if isexpr(f, :.)
                f.args[2].value = op
                expr.args[1] = f
            else
                expr.args[1] = op
            end
        end
        for a in Iterators.drop(expr.args, 1)
            if isa(a, Expr)
                replace_op!(a, op_map)
            end
        end
    else
        for a in expr.args
            if isa(a, Expr)
                replace_op!(a, op_map)
            end
        end
    end
    return expr
end

"""
    @checked expr

Perform all the operations in `expr` using checked arithmetic.

# Examples

```jldoctest
julia> 0xff + 0x10    # operation that overflows
0x0f

julia> @checked 0xff + 0x10
ERROR: OverflowError: 255 + 16 overflowed for type UInt8
```

You can also wrap method definitions (or blocks of code) in `@checked`:

```jldoctest
julia> plus(x, y) = x + y; minus(x, y) = x - y
minus (generic function with 1 method)

julia> @show plus(0xff, 0x10) minus(0x10, 0x20);
plus(0xff, 0x10) = 0x0f
minus(0x10, 0x20) = 0xf0

julia> @checked (plus(x, y) = x + y; minus(x, y) = x - y)
minus (generic function with 1 method)

julia> plus(0xff, 0x10)
ERROR: OverflowError: 255 + 16 overflowed for type UInt8

julia> minus(0x10, 0x20)
ERROR: OverflowError: 16 - 32 overflowed for type UInt8
```
"""
macro checked(expr)
    isa(expr, Expr) || return expr
    expr = copy(expr)
    return esc(replace_op!(expr, op_checked))
end

macro check(expr, kws...)
    isexpr(expr, :call) || error("expected :call expression, got ",
                                 isa(expr, Expr) ? QuoteNode(expr.head) : typeof(expr))
    safeexpr = copy(expr)
    for i = 2:length(expr.args)
        safeexpr.args[i] = Expr(:call, :(CheckedArithmetic.safearg), expr.args[i])
    end
    cmpexpr = isempty(kws) ? :(val == valcmp || error(val, " is not equal to ", valcmp)) :
                             :(isapprox(val, valcmp; $(esc(kws...))) || error(val, " is not approximately equal to ", valcmp))
    return quote
        local val = $(esc(expr))
        local valcmp = CheckedArithmetic.safeconvert(typeof(val), $(esc(safeexpr)))
        if ismissing(val) && ismissing(valcmp)
            val
        else
            $cmpexpr
            val
        end
    end
end

# safearg_type
# ------------
safearg_type(::Type{BigInt})  = BigInt
safearg_type(::Type{Int128})  = Int128
safearg_type(::Type{Int64})   = Int128
safearg_type(::Type{Int32})   = Int128
safearg_type(::Type{Int16})   = Int128
safearg_type(::Type{Int8})    = Int128
safearg_type(::Type{UInt128}) = UInt128
safearg_type(::Type{UInt64})  = UInt128
safearg_type(::Type{UInt32})  = UInt128
safearg_type(::Type{UInt16})  = UInt128
safearg_type(::Type{UInt8})   = UInt128

safearg_type(::Type{Bool}) = Bool   # these have a special meaning that must be preserved

safearg_type(::Type{BigFloat}) = BigFloat
safearg_type(::Type{Float64})  = Float64
safearg_type(::Type{Float32})  = Float64
safearg_type(::Type{Float16})  = Float64
safearg_type(::Type{T}) where T<:Base.TwicePrecision = T

safearg_type(::Type{<:Rational}) = Float64

# safearg
#--------
# Containers
safearg(t::Tuple) = map(safearg, t)
safearg(t::NamedTuple) = map(safearg, t)
safearg(p::Pair) = safearg(p.first) => safearg(p.second)
## AbstractArrays and similar
safearg(A::AbstractArray{T}) where T<:Number = convert(AbstractArray{safearg_type(T)}, A)
safearg(A::BitArray) = A
safearg(A::AbstractArray{Bool}) = A
safearg(rng::LinRange) = LinRange(safearg(rng.start), safearg(rng.stop), rng.len)
safearg(rng::StepRangeLen) = StepRangeLen(safearg(rng.ref), safearg(rng.step), rng.len, rng.offset)
safearg(rng::StepRange) = StepRange(safearg(rng.start), safearg(rng.step), safearg(rng.stop))
safearg(rng::AbstractUnitRange{Int}) = rng  # AbstractUnitRange{Int} is usually used for indexing, preserve this
safearg(rng::AbstractUnitRange{T}) where T<:Integer = convert(AbstractUnitRange{safearg_type(T)}, rng)
safearg(I::UniformScaling) = UniformScaling(safearg(I.λ))
safearg(F::Factorization{Float64}) = F
safearg(ref::Ref) = Ref(safearg(ref[]))
## AbstractDicts
safearg(d::Dict) = Dict(safearg(p) for p in d)
safearg(d::Base.EnvDict) = d
safearg(d::Base.ImmutableDict) = Base.ImmutableDict(safearg(p) for p in d)
safearg(d::Iterators.Pairs) = Iterators.Pairs(safearg(d.data), d.itr)
safearg(d::IdDict) = IdDict(safearg(p) for p in d)
safearg(d::WeakKeyDict) = WeakKeyDict(k=>safearg(v) for (k, v) in d)  # do not convert keys
## AbstractSets
safearg(s::Set) = Set(safearg(x) for x in s)
safearg(s::Base.IdSet) = Base.IdSet(safearg(x) for x in s)
safearg(s::BitSet) = s

# Other common types
safearg(::Nothing) = nothing
safearg(m::Missing) = m
safearg(s::Some) = Some(safearg(s.value))
safearg(rng::AbstractRNG) = rng
safearg(mode::RoundingMode) = mode
safearg(str::AbstractString) = str
safearg(c::AbstractChar) = c
safearg(sym::Symbol) = sym
safearg(mime::MIME) = mime
safearg(rex::Regex) = rex
safearg(rex::RegexMatch) = rex
safearg(chan::AbstractChannel) = chan
safearg(i::Base.AbstractCartesianIndex) = i
safearg(i::Base.UUID) = i
safearg(cmd::Base.AbstractCmd) = cmd
if isdefined(Base, :AbstractLock)
    safearg(lock::Base.AbstractLock) = lock
end
safearg(f::Function) = f
safearg(io::IO) = io
safearg(m::Module) = m
safearg(v::Val) = v
safearg(t::Task) = t
safearg(ord::Base.Order.Ordering) = ord
safearg(t::Timer) = t
safearg(d::Dates.AbstractDateTime) = d
safearg(t::Dates.AbstractTime) = t
safearg(t::Dates.AbstractDateToken) = t

# safeconvert
# -----------
safeconvert(::Type{T}, x) where T<:Integer = round(T, x)
safeconvert(::Type{T}, x) where T<:AbstractFloat = T(x)
safeconvert(::Type{AA}, A::AbstractArray) where AA<:AbstractArray{T} where T<:Integer = round.(T, A)

# accumulatortype
# ---------------
const SignPreserving = Union{typeof(+), typeof(*)}
const ArithmeticOp = Union{SignPreserving, typeof(-)}

accumulatortype(::ArithmeticOp, ::Type{BigInt})    = BigInt
accumulatortype(::ArithmeticOp, ::Type{Int128})    = Int128
accumulatortype(::ArithmeticOp, ::Type{Int64})     = Int64
accumulatortype(::ArithmeticOp, ::Type{Int32})     = Int
accumulatortype(::ArithmeticOp, ::Type{Int16})     = Int
accumulatortype(::ArithmeticOp, ::Type{Int8})      = Int
accumulatortype(::ArithmeticOp, ::Type{Bool})      = Int
accumulatortype(::SignPreserving, ::Type{UInt128}) = UInt128
accumulatortype(::SignPreserving, ::Type{UInt64})  = UInt64
accumulatortype(::SignPreserving, ::Type{UInt32})  = UInt
accumulatortype(::SignPreserving, ::Type{UInt16})  = UInt
accumulatortype(::SignPreserving, ::Type{UInt8})   = UInt
accumulatortype(::typeof(-), ::Type{UInt128})      = Int128
accumulatortype(::typeof(-), ::Type{UInt64})       = Int64
accumulatortype(::typeof(-), ::Type{UInt32})       = Int
accumulatortype(::typeof(-), ::Type{UInt16})       = Int
accumulatortype(::typeof(-), ::Type{UInt8})        = Int

accumulatortype(::ArithmeticOp, ::Type{BigFloat})  = BigFloat
accumulatortype(::ArithmeticOp, ::Type{Float64})   = Float64
accumulatortype(::ArithmeticOp, ::Type{Float32})   = Float64
accumulatortype(::ArithmeticOp, ::Type{Float16})   = Float64

end # module


text:
function merge(a::Array, b::Array)
  len = length(a) + length(b)
  B = Array(Int64, len)

  for k = 1:len

    if first(a) <= first(b)
      B[k] = shift!(a)
     else
       B[k] = shift!(b)
    end

    if isempty(a)
      return transfer_tail(B, b, k)
    elseif isempty(b)
      return transfer_tail(B, a, k)
    end

  end
end

function transfer_tail(arr::Array, tail::Array, count::Int64)
  for k = (count + 1):(length(tail) + count)
    arr[k] = shift!(tail)
  end
  return arr
end

function merge_sort(arr::Array)
    if length(arr) <= 1
       return arr
    end

    q = div(length(arr),2)

    a = merge_sort(arr[1 : q])
    b = merge_sort(arr[q + 1 : end])
    merge(a,b)
end

A = [3,6,5,8,99,2,1,5,9]
println(A)

A = merge_sort(A)
println(A)

max_stars_repo_name: wchmielarz/NobodyDoesItLikeJuliet

max_issues_repo_name: wchmielarz/NobodyDoesItLikeJuliet

max_forks_repo_name: wchmielarz/NobodyDoesItLikeJuliet


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

const max_ccall_threads = parse(Int, get(ENV, "UV_THREADPOOL_SIZE", "4"))
const thread_notifiers = Union{Condition, Nothing}[nothing for i in 1:max_ccall_threads]
const threadcall_restrictor = Semaphore(max_ccall_threads)

"""
    @threadcall((cfunc, clib), rettype, (argtypes...), argvals...)

The `@threadcall` macro is called in the same way as [`ccall`](@ref) but does the work
in a different thread. This is useful when you want to call a blocking C
function without causing the main `julia` thread to become blocked. Concurrency
is limited by size of the libuv thread pool, which defaults to 4 threads but
can be increased by setting the `UV_THREADPOOL_SIZE` environment variable and
restarting the `julia` process.

Note that the called function should never call back into Julia.
"""
macro threadcall(f, rettype, argtypes, argvals...)
    # check for usage errors
    isa(argtypes,Expr) && argtypes.head == :tuple ||
        error("threadcall: argument types must be a tuple")
    length(argtypes.args) == length(argvals) ||
        error("threadcall: wrong number of arguments to C function")

    # hygiene escape arguments
    f = esc(f)
    rettype = esc(rettype)
    argtypes = map(esc, argtypes.args)
    argvals = map(esc, argvals)

    # construct non-allocating wrapper to call C function
    wrapper = :(function (args_ptr::Ptr{Cvoid}, retval_ptr::Ptr{Cvoid})
        p = args_ptr
        # the rest of the body is created below
    end)
    body = wrapper.args[2].args
    args = Symbol[]
    for (i, T) in enumerate(argtypes)
        arg = Symbol("arg", i)
        push!(body, :($arg = unsafe_load(convert(Ptr{$T}, p))))
        push!(body, :(p += Core.sizeof($T)))
        push!(args, arg)
    end
    push!(body, :(ret = ccall($f, $rettype, ($(argtypes...),), $(args...))))
    push!(body, :(unsafe_store!(convert(Ptr{$rettype}, retval_ptr), ret)))
    push!(body, :(return Int(Core.sizeof($rettype))))

    # return code to generate wrapper function and send work request thread queue
    wrapper = Expr(Symbol("hygienic-scope"), wrapper, @__MODULE__)
    return :(let fun_ptr = @cfunction($wrapper, Int, (Ptr{Cvoid}, Ptr{Cvoid}))
        do_threadcall(fun_ptr, $rettype, Any[$(argtypes...)], Any[$(argvals...)])
    end)
end

function do_threadcall(fun_ptr::Ptr{Cvoid}, rettype::Type, argtypes::Vector, argvals::Vector)
    # generate function pointer
    c_notify_fun = @cfunction(
        function notify_fun(idx)
            global thread_notifiers
            notify(thread_notifiers[idx])
            return
        end, Cvoid, (Cint,))

    # cconvert, root and unsafe_convert arguments
    roots = Any[]
    args_size = isempty(argtypes) ? 0 : sum(Core.sizeof, argtypes)
    args_arr = Vector{UInt8}(undef, args_size)
    ptr = pointer(args_arr)
    for (T, x) in zip(argtypes, argvals)
        isbits(T) || throw(ArgumentError("threadcall requires isbits argument types"))
        y = cconvert(T, x)
        push!(roots, y)
        unsafe_store!(convert(Ptr{T}, ptr), unsafe_convert(T, y)::T)
        ptr += Core.sizeof(T)
    end

    # create return buffer
    ret_arr = Vector{UInt8}(undef, Core.sizeof(rettype))

    # wait for a worker thread to be available
    acquire(threadcall_restrictor)
    idx = findfirst(isequal(nothing), thread_notifiers)::Int
    thread_notifiers[idx] = Condition()

    GC.@preserve args_arr ret_arr roots begin
        # queue up the work to be done
        ccall(:jl_queue_work, Cvoid,
            (Ptr{Cvoid}, Ptr{UInt8}, Ptr{UInt8}, Ptr{Cvoid}, Cint),
            fun_ptr, args_arr, ret_arr, c_notify_fun, idx)

        # wait for a result & return it
        wait(thread_notifiers[idx])
        thread_notifiers[idx] = nothing
        release(threadcall_restrictor)

        r = unsafe_load(convert(Ptr{rettype}, pointer(ret_arr)))
    end
    return r
end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
"    BufferedHDF5Dataset2D{T}(g::Union{HDF5File,HDF5Group}, name, nbases, chunksize)"
mutable struct BufferedHDF5Dataset2D{T}
  ds::HDF5Dataset
  v::Vector{Vector{T}} #naievley you want a 2d array here, but there is no
  # equivalent of push!, so it uses less copying to have a vector of vectors
  # see https://github.com/JuliaLang/julia/issues/10546
  lastrow::Int64 # last index in hdf5 dataset
end
function BufferedHDF5Dataset2D{T}(g::Union{HDF5File,HDF5Group}, name, nbases, chunksize) where T
  ds = d_create(g, name, T, ((nbases,1), (nbases,-1)), "chunk", (nbases,chunksize))
  BufferedHDF5Dataset2D{T}(ds, Vector{Vector{T}}(),0)
end
function write_to_hdf5(b::BufferedHDF5Dataset2D)
  r = b.lastrow .+ (1:length(b.v))
  if length(r)>0
    set_dims!(b.ds, (size(b.ds,1), last(r)) )
    b.ds[:,r] = hcat(b.v...)
    empty!(b.v)
    b.lastrow=last(r)
  end
  return
end
Base.write(b::BufferedHDF5Dataset2D{T},x::Vector{T}) where T = push!(b.v,x)
Base.write(b::BufferedHDF5Dataset2D{T},x::Vector{Vector{T}}) where T = append!(b.v,x)

# d_create(h5, "ds", Float32, ((6,1), (6,-1)), "chunk", (6,chunksize))
# a single pulse is stored in a column, and accessed a[:,1]

struct BasisDataProduct <: DataProduct
    reduced::Vector{Float32} # reduced pulse, aka projection of pulse into subspace defined by basis
    residual_std::Float32
    frame1index::Int
    timestamp_usec::Int
    first_rising_sample::Int32
    nsamples::Int32
end
function Base.write(io::IO, d::BasisDataProduct)
  write(io, d.reduced)
  write(io, d.residual_std)
  write(io, d.frame1index)
  write(io, d.timestamp_usec)
  write(io, d.first_rising_sample)
  write(io, d.nsamples)
end

"""
    BasisAnalyzer(projectors,basis)
    `projectors` has size `(nbases,nsamples)`, basis has size '(nsamples, nabses)'
    comonly called as `BasisAnalyzer(projectors, projectors')` for testing or white noise case
An implementaion of AbstractBasisAnalyzer must have a field `projectors` that is an array with size ``(nbases, nsamples)` such that
the reduced pulse is calculated as `projectors*data` where `data` is a `Vector` of with length `nsamples`. It must also overload calling itself
to call `record2dataproduct`. See Julia issue #14919 for why this is neccesary.
"""
abstract type AbstractBasisAnalyzer end
struct BasisAnalyzer <: AbstractBasisAnalyzer
    projectors::Array{Float32,2} # size (nbases,nsamples)
    basis::Array{Float32,2} # size (nsamples,nbases)
    function BasisAnalyzer(projectors,basis)
      size(projectors) == reverse(size(basis)) || error("size(projectors) must == size(basis')")
      new(projectors,basis)
    end
end
nbases(a::AbstractBasisAnalyzer) = size(a.projectors,1)
nsamples(a::AbstractBasisAnalyzer) = size(a.projectors,2)
check_compatability(a::AbstractBasisAnalyzer, ljh) = nothing
(a::BasisAnalyzer)(record) = record2dataproduct(a,record)
modelreduce(a::AbstractBasisAnalyzer,data::AbstractVector)=a.projectors*data
modelpulse(a::AbstractBasisAnalyzer,reduced)=a.basis*reduced

function record2dataproduct(a::AbstractBasisAnalyzer,record::Union{LJH.LJHRecord, LJH.LJH3Record})
  reduced = modelreduce(a,LJH.data(record))
  modeled = modelpulse(a,reduced)
  residual_std = std(modeled-LJH.data(record))
  BasisDataProduct(reduced, residual_std, LJH.frame1index(record),
    LJH.timestamp_usec(record), LJH.first_rising_sample(record), length(record))
end

mutable struct BasisBufferedWriter <: BufferedWriter
    endchannel           ::Channel{Bool}
    timeout_s            ::Float64
    task                 ::Task
    reduced              ::BufferedHDF5Dataset2D{Float32}
    residual_std         ::BufferedHDF5Dataset{Float32}
    frame1index          ::BufferedHDF5Dataset{Int}
    timestamp_usec       ::BufferedHDF5Dataset{Int}
    first_rising_sample  ::BufferedHDF5Dataset{Int32}
    nsamples             ::BufferedHDF5Dataset{Int32}
end
function BasisBufferedWriter(h5::HDF5File, channel_number, nbases, chunksize, timeout_s;start=true)
  BasisBufferedWriter(g_create(h5,"$channel_number"),nbases,chunksize,timeout_s,start=start)
end
function BasisBufferedWriter(g::HDF5Group, nbases, chunksize, timeout_s;start=true)
  b=BasisBufferedWriter(Channel{Bool}(1), timeout_s, Task(nothing),
  BufferedHDF5Dataset2D{Float32}(g,"reduced", nbases, chunksize),
  BufferedHDF5Dataset{Float32}(g,"residual_std", chunksize),
  BufferedHDF5Dataset{Int}(g,"frame1index", chunksize),
  BufferedHDF5Dataset{Int}(g,"timestamp_usec", chunksize),
  BufferedHDF5Dataset{Int32}(g,"first_rising_sample", chunksize),
  BufferedHDF5Dataset{Int32}(g,"nsamples", chunksize), )
  start && schedule(b)
  b
end
"    hdf5file(b::BasisBufferedWriter)
Return the filename of the hdf5 file associated with `b`."
hdf5file(b::BasisBufferedWriter) = file(b.residual_std.ds)
function write_to_hdf5(b::BasisBufferedWriter)
  write_to_hdf5(b.reduced);write_to_hdf5(b.residual_std);write_to_hdf5(b.frame1index)
  write_to_hdf5(b.timestamp_usec);write_to_hdf5(b.first_rising_sample);write_to_hdf5(b.nsamples)
end
function Base.write(b::BasisBufferedWriter,x::BasisDataProduct)
  write(b.reduced, x.reduced)
  write(b.residual_std, x.residual_std)
  write(b.frame1index, x.frame1index)
  write(b.timestamp_usec, x.timestamp_usec)
  write(b.first_rising_sample, x.first_rising_sample)
  write(b.nsamples, x.nsamples)
  return nothing
end

function write_header(d::BasisBufferedWriter,r)
  channelgroup = parent(d.residual_std.ds)
  channelgroup["header"]="header"
end
function write_header_end(d::BasisBufferedWriter,ljh,analyzer::AbstractBasisAnalyzer)
  # dont add datasets, groups or attributes after SWMR writing is started
  # channelgroup = parent(d.filt_value.ds)
  # channelattrs = attrs(channelgroup)
end
function write_header_allchannel(d::BasisBufferedWriter, r::LJHReaderFeb2017)
  # should be called exact once when `schedule` is called on a Readers
  h5 = hdf5file(d)
  h5["BasisBufferedWriterHeader"]="BasisBufferedWriterHeader"
  h5["jointheader"]="jointheader"
  # execute only once for the whole HDF5 file
  flush(h5)
  try
    HDF5.start_swmr_write(h5)
  catch
    println("SKIPPING START_SWMR_WRITE, HDF5 VERSION TOO LOW")
  end
end

function make_buffered_hdf5_writer(h5, channel_number, analyzer::AbstractBasisAnalyzer, chunksize=1000, timeout_s=1.0)
  g = g_require(h5,"$channel_number")
  BasisBufferedWriter(g,nbases(analyzer), chunksize, timeout_s,start=true)
end


text:
# ------------------------------------------------------------------------------
# Jacobi preconditioner
# ------------------------------------------------------------------------------

"""
```julia
Jacobi(operator)
```

Jacobi diagonal preconditioner for finite element operators

# Arguments:
- `operator`:  finite element operator to precondition

# Returns:
- Jacobi preconditioner object

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
mass = GalleryOperator("mass", 4, 4, mesh);

# preconditioner
jacobi = Jacobi(mass);

# verify
println(jacobi)
println(jacobi.operator)

# output
jacobi preconditioner
finite element operator:
2d mesh:
    dx: 1.0
    dy: 1.0

2 inputs:
operator field:
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 4
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    interpolation
operator field:
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 4
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    quadratureweights

1 output:
operator field:
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 4
    numbercomponents: 1
    dimension: 2
  evaluation mode:
    interpolation
```
"""
mutable struct Jacobi <: AbstractPreconditioner
    # data never changed
    operator::Operator

    # data empty until assembled
    operatordiagonalinverse::AbstractArray{Float64}

    # inner constructor
    Jacobi(operator::Operator) = new(operator)
end

# printing
# COV_EXCL_START
Base.show(io::IO, preconditioner::Jacobi) = print(io, "jacobi preconditioner")
# COV_EXCL_STOP

# ------------------------------------------------------------------------------
# data for computing symbols
# ------------------------------------------------------------------------------

"""
```julia
getoperatordiagonalinverse(preconditioner)
```

Compute or retrieve the inverse of the symbol matrix diagonal for a Jacobi
    preconditioner

# Arguments:
- `preconditioner`:  preconditioner to compute diagonal inverse

# Returns:
- Symbol matrix diagonal inverse for the operator

# Example:
```jldoctest
# setup
mesh = Mesh1D(1.0);
diffusion = GalleryOperator("diffusion", 3, 3, mesh);

# preconditioner
jacobi = Jacobi(diffusion)

# note: either syntax works
diagonalinverse = LFAToolkit.getoperatordiagonalinverse(jacobi);
diagonalinverse = jacobi.operatordiagonalinverse;

# verify
@assert diagonalinverse ≈ [6/7 0; 0 3/4]

# output

```
"""
function getoperatordiagonalinverse(preconditioner::Jacobi)
    # assemble if needed
    if !isdefined(preconditioner, :operatordiagonalinverse)
        # retrieve diagonal and invert
        diagonalinverse = preconditioner.operator.diagonal^-1

        # store
        preconditioner.operatordiagonalinverse = diagonalinverse
    end

    # return
    return getfield(preconditioner, :operatordiagonalinverse)
end

# ------------------------------------------------------------------------------
# get/set property
# ------------------------------------------------------------------------------

function Base.getproperty(preconditioner::Jacobi, f::Symbol)
    if f == :operatordiagonalinverse
        return getoperatordiagonalinverse(preconditioner)
    else
        return getfield(preconditioner, f)
    end
end

function Base.setproperty!(preconditioner::Jacobi, f::Symbol, value)
    if f == :operator
        throw(ReadOnlyMemoryError()) # COV_EXCL_LINE
    else
        return setfield!(preconditioner, f, value)
    end
end

# ------------------------------------------------------------------------------
# compute symbols
# ------------------------------------------------------------------------------

"""
```julia
computesymbols(preconditioner, ω, θ)
```

Compute or retrieve the symbol matrix for a Jacobi preconditioned operator

# Arguments:
- `preconditioner`:  Jacobi preconditioner to compute symbol matrix for
- `ω`:               smoothing weighting factor array
- `θ`:               Fourier mode frequency array (one frequency per dimension)

# Returns:
- Symbol matrix for the Jacobi preconditioned operator

# Example:
```jldoctest
using LinearAlgebra

for dimension in 1:3
    # setup
    mesh = []
    if dimension == 1
        mesh = Mesh1D(1.0);
    elseif dimension == 2
        mesh = Mesh2D(1.0, 1.0);
    elseif dimension == 3
        mesh = Mesh3D(1.0, 1.0, 1.0);
    end
    diffusion = GalleryOperator("diffusion", 3, 3, mesh);

    # preconditioner
    jacobi = Jacobi(diffusion);

    # compute symbols
    A = computesymbols(jacobi, [1.0], π*ones(dimension));

    # verify
    using LinearAlgebra;
    eigenvalues = real(eigvals(A));
    if dimension == 1
        @assert maximum(eigenvalues) ≈ 1/7
    elseif dimension == 2
        @assert minimum(eigenvalues) ≈ -1/14
    elseif dimension == 3
        @assert minimum(eigenvalues) ≈ -0.33928571428571486
    end
end

# output

```
"""
function computesymbols(preconditioner::Jacobi, ω::Array, θ::Array)
    # validate number of parameters
    if length(ω) != 1
        Throw(error("exactly one parameter required for Jacobi smoothing")) # COV_EXCL_LINE
    end

    # return
    return I -
           ω[1] *
           preconditioner.operatordiagonalinverse *
           computesymbols(preconditioner.operator, θ)
end

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


text:
using .Mongoc

"""
The client and collection to use for this local cache.
"""
mutable struct LocalCache
  client::Mongoc.Client
  collection #::MongoCollection
end

global __localCache = nothing

global __forceOnlyLocalCache = false


text:
function q(ex, module′)
    @switch ex begin
    @case Expr(:macrocall, a, b, args...)  && let new_args = Any[a, b] end ||
          Expr(head, args...) && let new_args = [] end
        for i in eachindex(args)
            @switch args[i] begin
            @case ::LineNumberNode
            @case _
                push!(new_args, q(args[i], module′))
            end
        end
        if length(ex.args) !== length(new_args)
            ex.args = new_args
        end
        nothing
    @case _
    end
    ex
end

macro q(ex)
    esc(Expr(:quote, q(ex, __module__)))
end

macro codegen(ex)
    ex = __module__.eval(ex)
    if DEBUG
        @info macroexpand(__module__, ex)
    end
    __module__.eval(ex)
end

@inline function index(xs, x)
    @inbounds for i = eachindex(xs)
        if xs[i] === x
            return i
        end
    end
    nothing
end

@generated function GEP(a::Ptr{T}, ::Val{s}) where {T, s}
    off = fieldoffset(T, index(fieldnames(T), s))
    R = fieldtype(T, s)
    :(reinterpret(Ptr{$R}, a + $off))
end

function _pointer_access(@nospecialize(ex::Expr))
    @match ex begin
        :($a.$(s::Symbol) :: $A) =>
            :(reinterpret(Ptr{$A}, $GEP($(_pointer_access(a)),  $(Val(s)))))
        :($a.$(s::Symbol)) =>
            :($GEP($(_pointer_access(a)),  $(Val(s))))
        ex =>
            begin
                for i in eachindex(ex.args)
                    ex.args[i] = _pointer_access(ex.args[i])
                end
                ex
            end
    end
end

_pointer_access(a) = a

macro pointer_access(ex)
    esc(_pointer_access(ex))
end