text:
module Lar

	using NearestNeighbors
	using DataStructures
	using IntervalTrees
	#using Triangle
	using SparseArrays
	using LinearAlgebra
	#using Distributed
	using LinearAlgebraicRepresentation
	using LightXML

	LAR = LinearAlgebraicRepresentation

	"""
		Points = Array{Number,2}

	Alias declation of LAR-specific data structure.
	Dense `Array{Number,2,1}` ``M x N`` to store the position of *vertices* (0-cells)
	of a *cellular complex*. The number of rows ``M`` is the dimension
	of the embedding space. The number of columns ``N`` is the number of vertices.
	"""
	const Points = Matrix

	"""
		Cells = Array{Array{Int,1},1}

	Alias declation of LAR-specific data structure.
	Dense `Array` to store the indices of vertices of `P-cells`
	of a cellular complex.
	The linear space of `P-chains` is generated by `Cells` as a basis.
	Simplicial `P-chains` have ``P+1`` vertex indices for `cell` element in `Cells` array.
	Cuboidal `P-chains` have ``2^P`` vertex indices for `cell` element in `Cells` array.
	Other types of chain spaces may have different numbers of vertex indices for `cell`
	element in `Cells` array.
	"""
	const Cells = Array{Array{Int,1},1}

	const Cell = SparseVector{Int8, Int}

	"""
		Chain = SparseArrays.SparseVector{Int8,Int}

	Alias declation of LAR-specific data structure.
	Binary `SparseVector` to store the coordinates of a `chain` of `N-cells`. It is
	`nnz=1` with `value=1` for the coordinates of an *elementary N-chain*, constituted by
	a single *N-chain*.
	"""
	const Chain = SparseArrays.SparseVector{Int8,Int}

	"""
		ChainOp = SparseArrays.SparseMatrixCSC{Int8,Int}

	Alias declation of LAR-specific data structure.
	`SparseMatrix` in *Compressed Sparse Column* format, contains the coordinate
	representation of an operator between linear spaces of `P-chains`.
	Operators ``P-Boundary : P-Chain -> (P-1)-Chain``
	and ``P-Coboundary : P-Chain -> (P+1)-Chain`` are typically stored as
	`ChainOp` with elements in ``{-1,0,1}`` or in ``{0,1}``, for
	*signed* and *unsigned* operators, respectively.
	"""
	const ChainOp = SparseArrays.SparseMatrixCSC{Int8,Int}

	"""
		ChainComplex = Array{ChainOp,1}

	Alias declation of LAR-specific data structure. It is a
	1-dimensional `Array` of `ChainOp` that provides storage for either the
	*chain of boundaries* (from `D` to `0`) or the transposed *chain of coboundaries*
	(from `0` to `D`), with `D` the dimension of the embedding space, which may be either
	``R^2`` or ``R^3``.
	"""
	const ChainComplex = Array{ChainOp,1}

	include("./model.jl")
	include("./arrangement/arrangement.jl")
	include("lar2LAR.jl")
	include("file_extraction.jl")
	#include("./interface.jl")
    #include("./utilities.jl")
    #include("./simplexn.jl")
    #include("./largrid.jl")
    #include("./mapper.jl")
    #include("./struct.jl")
    #include("./integr.jl")
    #include("./refactoring.jl")
    #include("./graphalgorithms.jl")
    #include("./fragface.jl")
    #include("./svg2lar.jl")
    #include("./bool3d.jl")
    #include("./bool2d.jl")
end


text:
export binhist

probabilities(x) = _non0hist(x)
function _non0hist(x)
    L = length(x)

    hist = Vector{Float64}()
    # Reserve enough space for histogram:
    sizehint!(hist, L)

    sx = sort(x, alg = QuickSort)

    # Fill the histogram by counting consecutive equal values:
    prev_val, count = sx[1], 0
    for val in sx
        if val == prev_val
            count += 1
        else
            push!(hist, count)
            prev_val = val
            count = 1
        end
    end
    push!(hist, count)

    # Shrink histogram capacity to fit its size:
    sizehint!(hist, length(hist))
    return Probabilities(hist ./ L)
end
_non0hist(x::AbstractDataset) = _non0hist(x.data)
probabilities(x::AbstractDataset) = _non0hist(x.data)

ext: jl


text:
module Rendering
export merge_and_render_update
import WhereTraits
using WhereTraits: CONFIG
using WhereTraits.Utils
using WhereTraits.InternalState

using ExprParsers
using Markdown
using Setfield

abstract type RenderType end
struct RenderOuterFunc{Signature} <: RenderType
  outer::InternalState.DefOuterFunc{Signature}
end
struct RenderOuterAndInnerFuncs{Signature} <: RenderType
  outer::InternalState.DefOuterFunc{Signature}
  inners::InternalState.InnerFuncs
end
struct RenderInnerFunc{Signature} <: RenderType
  # we need the outerfunc to construct the innerfunc, as it depends on the ordering of the traits functions
  # which are defined in the outerfunc
  outer::InternalState.DefOuterFunc{Signature}
  inner::InternalState.DefInnerFunc
end

struct RenderDoc{Signature}
  deftraitsfunction::InternalState.DefTraitsFunction{Signature}
  inner::InternalState.DefInnerFunc
end

"""
    merge_and_render_update(env::MacroEnv, outfunc, innerfunc; [doc = true])

Merges the new traits definition (given by `outerfunc` and `innerfunc`) and renders the update.

Returns a list of Base.Expr which can be evaled or returned in macro.
"""
function merge_and_render_update(
    env::MacroEnv,
    outerfunc::InternalState.DefOuterFunc,
    innerfunc::InternalState.DefInnerFunc;
    doc = true)

  store = get_traitsstore(env.mod, outerfunc.fixed.name, outerfunc.fixed.signature)
  store_new, func_rendering, doc_rendering = _merge(env, store, outerfunc, innerfunc)

  exprs = [render(env, store_new, func_rendering)]
  if doc && CONFIG.auto_documentation
    push!(exprs, render_doc(env, store_new, doc_rendering))
  end
  push!(exprs, render_store_update(env, store_new))
  exprs
end

function render_store_update(env::MacroEnv, store::InternalState.TraitsStore)
  name = store.original_function
  if env.mod === name.mod
    # if we are rendering code for the same module, we need to drop the module information
    # this is needed for defining the function the very first time as `MyModule.func(...) = ...` is invalid syntax
    # for the initial definition
    name = name.name
  end
  :(function $(to_expr(name))(::WhereTraits.InternalState.TraitsStoreSingleton, ::Type{$(store.definitions.outer.fixed.signature)})
    $store
  end)
end

"""
merge the new traits information into the given traitsstore
and return whatever needs to be rendered for a correct update of the traits
"""
function _merge(env::MacroEnv,
                store::Union{Nothing, InternalState.TraitsStore},
                outerfunc::InternalState.DefOuterFunc,
                innerfunc::InternalState.DefInnerFunc)
  signature = outerfunc.fixed.signature
  if isnothing(store)
    mod_original, funcname_original = normalize_mod_and_name(env.mod, outerfunc.fixed.name)
    innerfuncs_new = InternalState.InnerFuncs(innerfunc.fixed => innerfunc.nonfixed)
    store_new = InternalState.TraitsStore(
      original_function = InternalState.Reference(mod_original, funcname_original),
      definitions = InternalState.DefTraitsFunction(outerfunc, innerfuncs_new)
    )
    func_rendering = RenderOuterAndInnerFuncs(outerfunc, innerfuncs_new)

  else
    outerfunc_old, innerfuncs = store.definitions.outer, store.definitions.inners

    innerfuncs_new = copy(innerfuncs)
    innerfuncs_new[innerfunc.fixed] = innerfunc.nonfixed

    outerfunc_new_nonfixed = InternalState.DefOuterFuncNonFixedPart(
      # we aggregate all unique traits and ensure order
      innerargs_traits = sortexpr(unique([outerfunc_old.nonfixed.innerargs_traits; outerfunc.nonfixed.innerargs_traits])),
    )

    if outerfunc_old.nonfixed == outerfunc_new_nonfixed  # if same WhereTraits, only the inner function needs to be rendered
      store_new = @set store.definitions = InternalState.DefTraitsFunction(outerfunc_old, innerfuncs_new)
      func_rendering = RenderInnerFunc(outerfunc_old, innerfunc)
    else
      outerfunc_new = InternalState.DefOuterFunc(
        # outerfunc.fixed == outerfunc_old.fixed, because of same signature
        fixed = outerfunc_old.fixed,
        nonfixed = outerfunc_new_nonfixed,
      )
      store_new = @set store.definitions = InternalState.DefTraitsFunction(outerfunc_new, innerfuncs_new)
      func_rendering = RenderOuterAndInnerFuncs(outerfunc_new, innerfuncs_new)
    end
  end

  # also return all the information about the state after this merge within an update variable
  doc_rendering = RenderDoc(store_new.definitions, innerfunc)
  store_new, func_rendering, doc_rendering
end

# Render
# ======

# we use special Singletons as separators to distinguish different kinds of parameters
struct _BetweenTypeVarsAndTraits end
struct _BetweenArgsAndTypeVars end

"""
render a whole TraitsStore

for debugging purposes only
"""
function render(env::MacroEnv, store::InternalState.TraitsStore)
  exprs = []
  for (outerfunc, innerfuncs) in values(store)
    push!(exprs, render(env, store, RenderOuterFunc(outerfunc)))
    for (fixed, nonfixed) in innerfuncs
      innerfunc = (fixed = fixed, nonfixed = nonfixed)
      push!(exprs, render(env, store, RenderInnerFunc(outerfunc, innerfunc)))
    end
  end
  flatten_blocks(Expr(:block, exprs...))
end

"""
rerender one single outerfunc and respective innerfuncs
"""
function render(env::MacroEnv, store::InternalState.TraitsStore, torender::RenderOuterAndInnerFuncs)
  outerfunc, innerfuncs = torender.outer, torender.inners
  exprs = []
  push!(exprs, render(env, store, RenderOuterFunc(outerfunc)))
  for (fixed, nonfixed) in innerfuncs
    innerfunc = InternalState.DefInnerFunc(fixed = fixed, nonfixed = nonfixed)
    push!(exprs, render(env, store, RenderInnerFunc(outerfunc, innerfunc)))
  end
  flatten_blocks(Expr(:block, exprs...))
end

function _map_args(new_to_old, innerargs)
  map(innerargs) do a
    # pure `_` is currently buggy, see https://github.com/JuliaLang/julia/issues/32727
    # hence we use ::Any instead
    get(new_to_old, a, Expr(:(::), Symbol("'", a, "'"), :(Any)))
  end
end

"""
render innerfunction
(this is only possible with informations from outerfunc)
"""
function render(env::MacroEnv, store::InternalState.TraitsStore, torender::RenderInnerFunc)
  outerfunc, innerfunc = torender.outer, torender.inner
  args = [
    :(::$(InternalState.TraitsDefSingleton));
    # we need to dispatch on the signature so that different outerfuncs don't
    # overwrite each other's innerfunc
    :(::$(Type{outerfunc.fixed.signature}));
    _map_args(innerfunc.fixed.args_mapping, outerfunc.fixed.innerargs_args);
    :(::$_BetweenArgsAndTypeVars);
    _map_args(innerfunc.fixed.typevars_mapping, outerfunc.fixed.innerargs_typevars);
    :(::$_BetweenTypeVarsAndTraits);
    _map_args(innerfunc.fixed.traits_mapping, outerfunc.nonfixed.innerargs_traits);
  ]
  # if we are rendering code for the same module, we need to drop the module information
  # this is needed for defining the function the very first time as `MyModule.func(...) = ...` is invalid syntax
  # for the initial definition
  name = store.original_function
  if env.mod === name.mod
    name = name.name
  end

  innerfunc_parsed = EP.Function_Parsed(
    name = name,
    curlies = [],
    args = args,
    kwargs = innerfunc.nonfixed.kwargs,
    wheres = [],
    body = innerfunc.nonfixed.body
  )
  to_expr(innerfunc_parsed)
end

"""
render outer function
"""
function render(env::MacroEnv, store::InternalState.TraitsStore, torender::RenderOuterFunc)
  outerfunc = torender.outer

  innerargs = [
    InternalState.traitsdefsingleton;
    # we need to dispatch on the signature so that different outerfuncs don't
    # overwrite each other's innerfunc
    outerfunc.fixed.signature;
    outerfunc.fixed.innerargs_args;
    _BetweenArgsAndTypeVars();
    outerfunc.fixed.innerargs_typevars;
    _BetweenTypeVarsAndTraits();
    outerfunc.nonfixed.innerargs_traits;
  ]
  innerfunc_call = EP.Call_Parsed(
    name = store.original_function,
    curlies = [],
    args = innerargs,
    kwargs = [:(kwargs...)],
  )
  # add LineNumberNode for debugging purposes
  body = Expr(:block, env.source, innerfunc_call)

  outerfunc_parsed = EP.Function_Parsed(
    name = outerfunc.fixed.name,
    curlies = outerfunc.fixed.curlies,
    args = outerfunc.fixed.args,
    kwargs = [:(kwargs...)],
    wheres = outerfunc.fixed.wheres,
    body = body,
  )
  to_expr(outerfunc_parsed)
end

"""
render documentation

extra effort needs to be done to properly document the outer function by referring
to innerfunctions
"""
function render_doc(env::MacroEnv, store::InternalState.TraitsStore, torender::RenderDoc)
  outerfunc = torender.deftraitsfunction.outer
  innerfuncs = torender.deftraitsfunction.inners
  innerfunc = torender.inner

  signature = to_expr(EP.Signature_Parsed(
    name = outerfunc.fixed.name,
    curlies = outerfunc.fixed.curlies,
    args = outerfunc.fixed.args,
    kwargs = [:(kwargs...)],
    wheres = outerfunc.fixed.wheres,
  ))

  # start documentation with autosignature of outer function
  header = Markdown.parse("""
  ```
  $signature
  ```
  ------ Original @traits definitions follow ------

  """)
  separator = Markdown.parse("- - -\n")

  doc_exprs = Any[header]
  for (fixed, nonfixed) in innerfuncs
    # automatic signature string of inner function
    signature_original = Markdown.parse("```julia\n$(nonfixed.expr_original.args[1])\n```")  # TODO this assumes that expr_original is a function, can we do this?
    push!(doc_exprs, signature_original)
    # manual doc string of respective inner function
    push!(doc_exprs, :(WhereTraits.Utils.DocsHelper.mygetdoc(
      $(to_expr(store.original_function)),
      Tuple{WhereTraits.InternalState.TraitsDocSingleton,
            Type{$(outerfunc.fixed.signature)},
            Type{$(innerfunc_fixed_to_doctype(fixed))}}
    )))
    # automatic doc string of inner function definition
    expr_original = Markdown.parse("Original @traits definition:\n```julia\n$(nonfixed.expr_original)\n```")
    push!(doc_exprs, expr_original)
    # better visual separation
    push!(doc_exprs, separator)
  end
  # get rid of last separator
  deleteat!(doc_exprs, lastindex(doc_exprs))

  # if we are rendering code for the same module, we need to drop the module information
  # this is needed for defining the function the very first time as `MyModule.func(...) = ...` is invalid syntax
  # for the initial definition
  name = store.original_function
  if env.mod === name.mod
    name = name.name
  end
  name = to_expr(name)

  quote
    # first the documentation of the inner function as this needs to be updated BEFORE the outer doc-string
    # is updated below
    Base.@__doc__ function $name(::$(InternalState.TraitsDocSingleton),
                                 ::Type{$(outerfunc.fixed.signature)},
                                 ::Type{$(innerfunc_fixed_to_doctype(innerfunc.fixed))}) end

    # documentation of outer function (we need to manually ignore nothing docs)
    let docstring = Base.Docs.catdoc(filter(!isnothing, [$(doc_exprs...)])...)
      WhereTraits.Utils.DocsHelper.@doc_signature docstring ($signature)
    end
  end
end

struct InnerFuncFixedDocSig{Args, TypeVars, WhereTraits} end
function innerfunc_fixed_to_doctype(innerfunc_fixed)
  dicts = [innerfunc_fixed.args_mapping, innerfunc_fixed.typevars_mapping, innerfunc_fixed.traits_mapping]
  types = Dict_to_normalizedType.(dicts)
  InnerFuncFixedDocSig{types...}
end

end # module


text:
@testset "buildSquare" begin
    mesh = buildSquare(0, 1, 0, 1, 2, 2)

    # Element 1
    @test mesh.elements[1].nodes[1].coords == Vec{2}((0.,0.))
    @test mesh.elements[1].nodes[2].coords == Vec{2}((0.5,0.))
    @test mesh.elements[1].nodes[3].coords == Vec{2}((0.5,0.5))
    @test mesh.elements[1].nodes[4].coords == Vec{2}((0.,0.5))

    # Element 4
    @test mesh.elements[4].nodes[1].coords == Vec{2}((0.5,0.5))
    @test mesh.elements[4].nodes[2].coords == Vec{2}((1.,0.5))
    @test mesh.elements[4].nodes[3].coords == Vec{2}((1.,1.))
    @test mesh.elements[4].nodes[4].coords == Vec{2}((.5,1.))

    # Boundary Info
    boundary_info = mesh.boundary_info

    @test boundary_info.sidesets == [1,2,3,4]
    @test boundary_info.nodesets == [1,2,3,4]

    # Sideset 1
    @test length(boundary_info.side_list[1]) == 2
    @test boundary_info.side_list[1][1].element.id == 1
    @test boundary_info.side_list[1][1].side == 1
    @test boundary_info.side_list[1][2].element.id == 2
    @test boundary_info.side_list[1][2].side == 1

    # Sideset 2
    @test length(boundary_info.side_list[2]) == 2
    @test boundary_info.side_list[2][1].element.id == 2
    @test boundary_info.side_list[2][1].side == 2
    @test boundary_info.side_list[2][2].element.id == 4
    @test boundary_info.side_list[2][2].side == 2

    # Sideset 3
    @test length(boundary_info.side_list[3]) == 2
    @test boundary_info.side_list[3][1].element.id == 3
    @test boundary_info.side_list[3][1].side == 3
    @test boundary_info.side_list[3][2].element.id == 4
    @test boundary_info.side_list[3][2].side == 3

    # Sideset 4
    @test length(boundary_info.side_list[4]) == 2
    @test boundary_info.side_list[4][1].element.id == 1
    @test boundary_info.side_list[4][1].side == 4
    @test boundary_info.side_list[4][2].element.id == 3
    @test boundary_info.side_list[4][2].side == 4
end


text:
function iterate!(ph::AbstractProgressiveHedgingSolver, ::AbstractExecution)
    # Resolve all subproblems at the current optimal solution
    Q = resolve_subproblems!(ph)
    if Q == Inf
        return :Infeasible
    elseif Q == -Inf
        return :Unbounded
    end
    ph.data.Q = Q
    # Update iterate
    update_iterate!(ph)
    # Update subproblems
    update_subproblems!(ph)
    # Get dual gap
    update_dual_gap!(ph)
    # Update penalty (if applicable)
    update_penalty!(ph)
    # Update progress
    @unpack δ₁, δ₂ = ph.data
    ph.data.δ = sqrt(δ₁ + δ₂)/(1e-10+norm(ph.ξ,2))
    # Log progress
    log!(ph)
    # Check optimality
    if check_optimality(ph)
        # Optimal
        return :Optimal
    end
    # Just return a valid status for this iteration
    return :Valid
end

function start_workers!(::AbstractProgressiveHedgingSolver, ::AbstractExecution)
    return nothing
end

function close_workers!(::AbstractProgressiveHedgingSolver, ::AbstractExecution)
    return nothing
end


text:
function choose_block_size(C, A, B, ::Nothing)
    if (*)(length(C) |> Int128, length(A) |> Int128, length(B) |> Int128) >= ((3default_block_size()) >>> 1)^6
        default_block_size()
    else
        32
    end
end

choose_block_size(C, A, B, block_size::Integer) = block_size


text:
struct DiffusionMatrix
    M
    coeff
    Δx
end
(deriv::DiffusionMatrix)(x::AbstractVector) = deriv.coeff*(deriv.M*x)

function derivative_matrix_first_order(n, xl, xu)
    Δx = (xu - xl)/(n - 1)
    r8fΔx = 1/(40320Δx)

    first_block = [
        -109584.0 +322560 -564480 +752640 -705600 +451584 -188160 +46080 -5040
        -5040.0 -64224 +141120 -141120 +117600 -70560 +28224 -6720 +720
        +720.0 -11520 -38304 +80640 -50400 +26880 -10080 +2304 -240
        -240.0 +2880 -20160 -18144 +50400 -20160 +6720 -1440 +144
        ]
    first_block = [first_block zeros(Float64, 4, n-9)]

    i_th_row = [+144.0 -1536 +8064 -32256 +0 +32256 -8064 +1536 -144]

    last_block = [
        -144.0 +1440 -6720 +20160 -50400 +18144 +20160 -2880 +240
        +240.0 -2304 +10080 -26880 +50400 -80640 +38304 +11520 -720
        -720.0 +6720 -28224 +70560 -117600 +141120 -141120 +64224 +5040
        +5040.0 -46080 +188160 -451584 +705600 -752640 +564480 -322560 +109584
        ]
    last_block = [zeros(Float64, 4, n-9) last_block]

    mid_block = zeros(Float64, n-8, n)
    @inbounds for (row_index, i) in enumerate(5:n-4)
        @inbounds mid_block[row_index,row_index:row_index+8] = i_th_row
    end

    deriv_matrix = [
        first_block
        mid_block
        last_block
        ]

    return DiffusionMatrix(deriv_matrix, r8fΔx, Δx)
end

function derivative_matrix_second_order(n, xl, xu)
    Δx = (xu - xl)/(n - 1)

    r12Δxs = 1/(12Δx^2)

    mid_block = zeros(Float64, n - 4, n)

    first_block = [
        -415/6 +96 -36 +32/3 -3/2 0
        +10.0 -15 -4 +14 -6 +1
        ]
    first_block = [first_block zeros(Float64, 2, n-6)]

    i_th_row = [-1.0 +16 -30 +16 -1]

    last_block = [
        +1.0 -6 +14 -4 -15 +10
        0.0 -3/2 +32/3 -36 +96 -415/6
        ]
    last_block = [zeros(Float64, 2, n-6) last_block]

    @inbounds for (row_index, i) in enumerate(3:n-2)
        @inbounds mid_block[row_index, row_index:row_index+4] = i_th_row
    end

    deriv_matrix = [
        first_block
        mid_block
        last_block
        ]

    return DiffusionMatrix(deriv_matrix, r12Δxs, Δx)
end

function derivative_matrices_first_and_second_order(N, xl=0, xu=1)
    """
    Combines the first and second order derivative matrices into a dictionary
    """
    out = Dict{Int64,DiffusionMatrix}()
    out[1] = derivative_matrix_first_order(N, xl, xu)
    out[2] = derivative_matrix_second_order(N, xl, xu)

    return out
end

function block_matrix_maker(p, X, Y, Z)
    """
    Finite volume stencil for the electrolyte grid (cathode, separator, anode)
    """
    single_block(x) = Matrix(Tridiagonal{eltype(x)}(-x[1:end-1],x .+ [0;x[1:end-1]],-x[1:end-1]))

    ind_p = (1:p.N.p)
    ind_s = (1:p.N.s) .+ (p.N.p)
    ind_n = (1:p.N.n) .+ (p.N.p+p.N.s)

    A_tot = zeros(eltype(X), ind_n[end], ind_n[end])

    A_tot[ind_p,ind_p] .= single_block(X)
    A_tot[ind_s,ind_s] .= single_block(Y)
    A_tot[ind_n,ind_n] .= single_block(Z)

    return A_tot
end

function interpolate_electrolyte_conductivities(Keff_p, Keff_s, Keff_n, p::AbstractModel)
    """
    interpolate_electrolyte_conductivities interpolates electrolyte conductivities at the edges of control volumes using harmonic mean.
    """

    Δx = Δx_values(p.N)

    Keff_i_medio(β_i, Keff_i) = Keff_i[1:end-1].*Keff_i[2:end]./(β_i*Keff_i[2:end]+(1-β_i)*Keff_i[1:end-1])
    β_i_j(Δx_i, l_i, Δx_j, l_j) = Δx_i*l_i/2 /(Δx_j*l_j/2+Δx_i*l_i/2)
    Keff_i_j_interface(β_i_j, Keff_i, Keff_j) = Keff_i[end]*Keff_j[1] / (β_i_j*Keff_j[1] + (1-β_i_j)*Keff_i[end])

    ## Positive electrode mean conductivity
    β_p = 0.5
    Keff_p_medio = Keff_i_medio(β_p, Keff_p)

    # The last element of Keff_p_medio will be the harmonic mean of the
    # elements at the interface positive-separator

    β_p_s = β_i_j(Δx.p, p.θ[:l_p], Δx.s, p.θ[:l_s])

    Keff_p_s_interface = Keff_i_j_interface(β_p_s, Keff_p, Keff_s)

    append!(Keff_p_medio, Keff_p_s_interface)

    ## Separator mean conductivity
    # Compute the harmonic mean values for the separator
    β_s = 0.5
    Keff_s_medio = Keff_i_medio(β_s, Keff_s)

    # The last element of Keff_s_medio will be the harmonic mean of the
    # elements at the interface separator-negative

    β_s_n = β_i_j(Δx.s, p.θ[:l_s], Δx.n, p.θ[:l_n])

    Keff_s_n_interface = Keff_i_j_interface(β_s_n, Keff_s, Keff_n)

    append!(Keff_s_medio, Keff_s_n_interface)

    ## Negative electrode mean conductivity
    # Compute the harmonic mean values for the negative electrode

    β_n = 0.5
    Keff_n_medio = Keff_i_medio(β_n, Keff_n)

    append!(Keff_n_medio, 0.0) # the cc interface is not used. The zero is only to match the dimensions

    return Keff_p_medio, Keff_s_medio, Keff_n_medio

end

harmonic_mean(β, x₁, x₂) = @. x₁*x₂/(β*x₂ + (1.0 - β)*x₁)

function interpolate_electrolyte_concentration(c_e, p::AbstractModel)
    """
    interpolate_electrolyte_concentration interpolates the value of electrolyte concentration at the edges of control volumes using harmonic mean
    """

    Δx = Δx_values(p.N)

    ## Electrolyte concentration interpolation

    # Interpolation within the positive electrode
    β_c_e_p = 0.5
    @inbounds @views c̄_e_p = harmonic_mean(β_c_e_p, c_e[1:p.N.p-1], c_e[2:p.N.p])

    # Interpolation on the interface between separator & positive electrode
    @inbounds @views β_c_e_ps = Δx.p*p.θ[:l_p]/2 / (Δx.p*p.θ[:l_p]/2 + Δx.s*p.θ[:l_s]/2)
    @inbounds @views c̄_e_ps = harmonic_mean(β_c_e_ps, c_e[p.N.p], c_e[p.N.p+1])

    # Interpolation within the separator
    β_c_e_s = 0.5
    @inbounds @views c̄_e_s = harmonic_mean(β_c_e_s, c_e[p.N.p+1:p.N.p+p.N.s-1], c_e[p.N.p+2:p.N.p+p.N.s])

    # Interpolation on the interface between separator & negative electrode
    @inbounds @views β_c_e_sn = Δx.s*p.θ[:l_s]/2 / (Δx.n*p.θ[:l_n]/2 + Δx.s*p.θ[:l_s]/2)
    @inbounds @views c̄_e_sn = harmonic_mean(β_c_e_sn, c_e[p.N.p+p.N.s], c_e[p.N.p+p.N.s+1])

    # Interpolation within the negative electrode
    β_c_e_n = 0.5
    @inbounds @views c̄_e_n = harmonic_mean(β_c_e_n, c_e[p.N.p+p.N.s+1:end-1], c_e[p.N.p+p.N.s+2:end])

    return [c̄_e_p; c̄_e_ps], [c̄_e_s; c̄_e_sn], c̄_e_n

end

interpolate_temperature(T, p::AbstractModel) = interpolate_electrolyte_concentration((@views @inbounds T[p.N.a+1:end-p.N.z]),p)

function interpolate_electrolyte_concetration_fluxes(c_e, p::AbstractModel)
    """
    interpolate_electrolyte_concetration_fluxes interpolates the electrolyte concentration flux at the edges of control volumes using harmonic mean
    """
    Δx = Δx_values(p.N)

    # Fluxes within the positive electrode
    @inbounds @views ∂ₓc_e_p = (c_e[2:p.N.p] .- c_e[1:p.N.p-1])/(Δx.p*p.θ[:l_p])

    # Fluxes at the separator-positive interface
    @inbounds @views ∂ₓc_e_ps = (c_e[p.N.p+1] .- c_e[p.N.p]) / ((Δx.p*p.θ[:l_p]/2+Δx.s*p.θ[:l_s]/2))

    # Fluxes within the separator
    @inbounds @views ∂ₓc_e_s = (c_e[p.N.p+2:p.N.p+p.N.s] .- c_e[p.N.p+1:p.N.p+p.N.s-1])/(Δx.s*p.θ[:l_s])

    # Fluxes at the separator-negative interface
    @inbounds @views ∂ₓc_e_sn = (c_e[p.N.p+p.N.s+1] .- c_e[p.N.p+p.N.s]) / ((Δx.n*p.θ[:l_n]/2+Δx.s*p.θ[:l_s]/2))

    # Fluxes within the negative electrode
    @inbounds @views ∂ₓc_e_n = (c_e[p.N.p+p.N.s+2:end] .- c_e[p.N.p+p.N.s+1:end-1])/(Δx.n*p.θ[:l_n])

    return [∂ₓc_e_p; ∂ₓc_e_ps], [∂ₓc_e_s; ∂ₓc_e_sn], ∂ₓc_e_n
end

Δx_values(N) = (p=1/N.p, s=1/N.s, n=1/N.n, a=1/N.a, z=1/N.z)


text:
module Molecules

using ForwardDiff
using LinearAlgebra

export @Orbital_str, @orbital_str,
    chemicalpotential, hamiltonian, makeaxes, molecularorbitals, overlap,
    selfenergy, slaterkoster, total_energy, vibrationalcoupling

const ⊗ = kron
const σ0 = I(2)
const σ = (
    x = [0.0 1.0; 1.0 0.0],
    y = [0.0 -1.0im; 1.0im 0.0],
    z = [1.0 0.0; 0.0 -1.0]
)

struct Orbital{N,L,M} end
macro Orbital_str(str::String)
    reg = r"^([0-9\*]+)([a-zA-Z]+)(?:_(?:\{(.+)\}|([a-zA-Z\*]+)))?$"
    m = match(reg, str)
    m === nothing && error("Invalid orbital")
    N = m[1] === nothing ? error("Invalid orbital") : m[1] == "*" ? nothing : parse(Int, m[1])
    L = m[2] === nothing ? error("Invalid orbital") : Symbol(m[2])
    M = m[3] === nothing ? m[4] === nothing ? Nothing : Symbol(m[4]) : Symbol(m[3])
    M === Symbol("*") && (M = nothing)

    N === nothing && M === nothing && return Orbital{<:Any,L,<:Any}
    N === nothing && return Orbital{<:Any,L,M}
    M === nothing && return Orbital{N,L,<:Any}
    return Orbital{N,L,M}
end

macro orbital_str(str::String)
    reg = r"^([0-9]+)([a-zA-Z]+)(?:_(?:\{(.+)\}|([a-zA-Z]+)))?$"
    m = match(reg, str)
    m === nothing && error("Invalid orbital")
    N = m[1] === nothing ? error("Invalid orbital") : parse(Int, m[1])
    L = m[2] === nothing ? error("Invalid orbital") : Symbol(m[2])
    M = m[3] === nothing ? m[4] === nothing ? Nothing : Symbol(m[4]) : Symbol(m[3])
    return Orbital{N,L,M}()
end

function Base.show(io::IO, ::Orbital{N,L,M}) where {N,L,M}
    if M === Nothing
        print(io, "orbital\"$N$L\"")
    elseif match(r"^[a-zA-Z]$", string(M)) !== nothing
        print(io, "orbital\"$N$(L)_$M\"")
    else
        print(io, "orbital\"$N$(L)_{$M}\"")
    end
end

function shortform(::Orbital{N,L,M}) where {N,L,M}
    l = L
    m = M == Nothing ? "" : replace(string(M), r"[_^+\-{}]" => s"")
    return Symbol("$l$m")
end

mutable struct Atom{SYMBOL,Z,N}
    position
    axes
    valence::Union{Int,Missing}
    orbitals

    function Atom{SYMBOL,Z,N}(position, axes, valence, orbitals::Orbital...) where {SYMBOL,Z,N}
        return new{SYMBOL,Z,N}(
            position,
            axes,
            valence,
            orbitals
        )
    end
end

function Atom{SYMBOL,Z,N}(position, orbitals::Orbital...; valence=missing, axes=makeaxes(0., 0.)) where {SYMBOL,Z,N}
    return Atom{SYMBOL,Z,N}(position, axes, valence, orbitals...)
end

Base.:(==)(atom1::Atom, atom2::Atom) = false
function Base.:(==)(atom1::T, atom2::T) where {T <: Atom}
    return (
        atom1.position == atom2.position &&
        atom1.orbitals == atom2.orbitals &&
        isequal(atom1.valence, atom2.valence)
    )
end

Base.getindex(atom::Atom, i::Integer) = collect(atom.orbitals)[i]

position(atom::Atom) = atom.position
atomic_positions(atoms::AbstractVector{<:Atom}) = position.(atoms)

struct Bond{A1<:Atom,A2<:Atom}
    atom1::A1
    atom2::A2

    function Bond{A1,A2}(atom1::A1, atom2::A2) where {A1<:Atom,A2<:Atom}
        if hash(atom1) > hash(atom2)
            return new{A1,A2}(atom1, atom2)
        else
            return new{A2,A1}(atom2, atom1)
        end
    end
end
Bond(atom1::A1, atom2::A2) where {A1, A2} = Bond{A1,A2}(atom1, atom2)

mutable struct Molecule
    position
    axes
    atoms::Vector{Atom}
    bonds::Set{Bond}
end

Molecule(atoms=Vector{Atom}(), bonds=Set{Bond}()) = Molecule(
    [0., 0., 0.],
    makeaxes(0., 0.),
    collect(Atom, atoms),
    Set{Bond}(bonds))

function Base.:(==)(mol1::Molecule, mol2::Molecule)
    return (
        mol1.position == mol2.position &&
        length(mol1.atoms) == length(mol2.atoms) &&
        all(atom1 -> any(atom2 -> atom1 == atom2, mol2.atoms), mol1.atoms) &&
        all(atom2 -> any(atom1 -> atom1 == atom2, mol1.atoms), mol2.atoms)
    )
end

Base.push!(mol::Molecule, atoms::Atom...) = push!(mol.atoms, atoms...)
Base.push!(mol::Molecule, bonds::Bond...) = push!(mol.bonds, bonds...)

Base.eachindex(mol::Molecule) = eachindex(collect(mol.atoms))
Base.iterate(mol::Molecule) = iterate(collect(mol.atoms))
Base.iterate(mol::Molecule, i) = iterate(collect(mol.atoms), i)
Base.length(mol::Molecule) = length(mol.atoms)
Base.getindex(mol::Molecule, i::Integer) = collect(mol.atoms)[i]
Base.lastindex(mol::Molecule) = lastindex(collect(mol.atoms))

function indices(mol::Molecule, atom::Atom)
    N0 = 1
    for atom′ in mol.atoms
        N1 = N0 + countorbitals(atom′)
        if atom === atom′
            return N0:(N1-1)
        end
        N0 = N1
    end
    return nothing
end

function indices(mol::Molecule, atom::Atom, orbital::Orbital)
    N0 = 1
    for atom′ in mol.atoms
        N1 = N0 + countorbitals(atom′)
        if atom === atom′
            for (i, orbital′) in enumerate(atom.orbitals)
                if orbital == orbital′
                    return (N0+i-1):(N0+i-1)
                end
            end
            return nothing
        end
        N0 = N1
    end
    return nothing
end

Base.Vector{T}(mol::Molecule) where {T} = normalize!(fill(oneunit(T), countorbitals(mol)))

function Base.Vector{T}(mol::Molecule, atom::Atom) where {T}
    u = fill(zero(T), countorbitals(mol))
    N = countorbitals(atom)
    u[indices(mol, atom)] .= 1 / sqrt(N)
    return u
end

function Base.Vector{T}(mol::Molecule, n::Integer) where {T}
    atom = mol[n]
    return Vector{T}(mol, atom)
end

function Base.Vector{T}(mol::Molecule, atom::Atom, orbital::Orbital) where {T}
    u = fill(zero(T), countorbitals(mol))
    u[indices(mol, atom, orbital)] .= 1
    return u
end

function Base.Vector{T}(mol::Molecule, atom::Atom, s::Integer) where {T}
    orbital = atom[s]
    return Vector{T}(mol, atom, orbital)
end

function Base.Vector{T}(mol::Molecule, n::Integer, s::Integer) where {T}
    atom = mol[n]
    orbital = atom[s]
    return Vector{T}(mol, atom, orbital)
end

function centre_of_mass(f, mol::Molecule)
    return sum(mass(f, atom)*position(atom) for atom in mol) / mass(f, mol)
end

function balance_position!(f, mol::Molecule)
    p0 = centre_of_mass(f, mol)
    for atom in mol.atoms
        atom.position -= p0
    end
    return mol
end

##################################################
# Chemical properties
##################################################

protonnumber(::Atom{SYMBOL,Z,N}) where {SYMBOL,Z,N} = Z
neutronnumber(::Atom{SYMBOL,Z,N}) where {SYMBOL,Z,N} = N
chemical_symbol(::Atom{SYMBOL,Z,N}) where {SYMBOL,Z,N} = SYMBOL
chemical_symbols(atoms::AbstractVector{<:Atom}) = chemical_symbol.(atoms)

countorbitals(atom::Atom)::Int = length(atom.orbitals)
countorbitals(mol::Molecule)::Int = sum(map(countorbitals, mol.atoms))

mass(f, atom::Atom) = f(chemical_symbol(atom))
mass(f, mol::Molecule) = sum(map(atom -> mass(f, atom), mol))

function electronvalence(atom::Atom)
    !ismissing(atom.valence) && return atom.valence
    n = protonnumber(atom)
    n <= 2 && return n
    n <= 10 && return n - 2
    n <= 18 && return n - 10 # Not very generic ...
    n > 18 && error("atom numbers greater than 18 are currently not supported")
end
"""
    valencecount(mol)

Return the sum of electron valences of each atom in `mol`.

"""
valencecount(mol::Molecule, oxidation=0) = sum(electronvalence(atom) for atom in mol.atoms) - oxidation

function molecularorbitals(H)
    vals, vecs = eigen(H)
    s = sortperm(vals)
    return vals[s], vecs[:, s]
end

function chemicalpotential(mol::Molecule, levels, oxidation=0)
    n = valencecount(mol, oxidation)
    μ = real(levels[ceil(Int, n/2)] + levels[ceil(Int, (n + 1)/2)])/2
    Δ = real(levels[ceil(Int, (n + 1)/2)] - levels[ceil(Int, n/2)])/2
    return μ, Δ
end

bond_length(bond::Bond) = norm(position(bond.atom2) - position(bond.atom1))

function total_energy(mol::Molecule, levels; oxidation = 0)
    n = valencecount(mol, oxidation)
    return sum(i -> real(levels[i]), 1:ceil(Int, n/2))
end

##################################################
# Hamiltonians
##################################################

include("hamiltonian.jl")
include("overlap.jl")

include("unitcell.jl")

##################################################
# Hamiltonian/overlap schemes
##################################################

include("slaterkoster.jl")

##################################################
# Elements
##################################################

chemical_symbol_to_atom_type(x::Symbol) = chemical_symbol_to_atom_type(Val(x))

const Hydrogen = Atom{:H, 1, 0}
Hydrogen(position; valence=1, kwargs...) = Hydrogen(position, orbital"1s"; valence=valence, kwargs...)
chemical_symbol_to_atom_type(::Val{:H}) = Hydrogen

const Deuterium = Atom{:D, 1, 1}
const Tritium = Atom{:T, 1, 2}

const Carbon = Atom{:C, 6, 6}
Carbon(position; valence=4, kwargs...) = Carbon(position, orbital"2s", orbital"2p_x", orbital"2p_y", orbital"2p_z"; valence=valence, kwargs...)
chemical_symbol_to_atom_type(::Val{:C}) = Carbon

const Carbon13 = Atom{:C, 6, 7}
const Nitrogen = Atom{:N, 7, 7}
const Oxygen = Atom{:O, 8, 8}

##################################################
# Coordinates
##################################################

function makeaxes(θ = 0.0, ϕ = 0.0, ρ = 0.0)
    R1 = LinearAlgebra.Givens(2, 1, cos(ρ), sin(ρ))
    R2 = LinearAlgebra.Givens(3, 1, cos(ϕ), sin(ϕ))
    R3 = LinearAlgebra.Givens(2, 1, cos(θ), sin(θ))
    return R3*(R2*(R1*Matrix(I, 3, 3)))
end

##################################################
# Basis change
##################################################

localbasis(atom::Atom) = localbasis(atom.axes, atom.orbitals...)

function localbasis(axes::AbstractMatrix{T}, ::Orbital"*s") where {T}
    return fill(one(T), 1, 1)
end

function localbasis(axes::AbstractMatrix{T}, ::Orbital"*s", ::Orbital"*p_x", ::Orbital"*p_y", ::Orbital"*p_z") where {T}
    Z = zero(T) / oneunit(T)
    return [one(T) fill(Z, 1, 3)
            fill(Z, 3, 1) axes]
end

function localbasis(axes::AbstractMatrix{T}, ::Orbital"*p_x", ::Orbital"*p_y") where {T}
    Z = zero(T) / oneunit(T)
    return axes[1:2,1:2]
end

##################################################
# Position operator
##################################################

_position(::Val{:x}, atom::Atom) = atom.position[1]
_position(::Val{:y}, atom::Atom) = atom.position[2]
_position(::Val{:z}, atom::Atom) = atom.position[3]

position(x::Symbol, args...) = position(Val(x), args...)

position(::Val, s1::Orbital, s2::Orbital) = s1 == s2 ? 1.0 : 0.0

function position(val::Val, atom::Atom)
    N = countorbitals(atom)
    U = localbasis(atom)
    X = fill(0.0im, N, N)
    for (j, s2) in enumerate(atom.orbitals), (i, s1) in enumerate(atom.orbitals)
        X[i, j] = _position(val, atom) * position(val, s1, s2)
    end
    return U'X*U
end

function position(val::Val, mol::Molecule)
    N = countorbitals(mol)
    X = fill(0.0im, N, N)

    i = 1
    for atom in mol.atoms
        Ni = countorbitals(atom)
        X[i:i+Ni-1, i:i+Ni-1] = position(val, atom)
        i += Ni
    end

    return X
end

# function position(val::Val, mol::Molecule, n)
#     N = countorbitals(mol)
#     X = fill(0.0im, N, N)
#     inds = indices(mol, mol[n])
#     X[inds, inds] = position(val, mol[n])
#     return X
# end

##################################################
# Angular momentum
##################################################

angularmomentum(vec::AbstractVector, args...) =
    vec[1]*angularmomentum(:x, args...) +
    vec[2]*angularmomentum(:y, args...) +
    vec[3]*angularmomentum(:z, args...)

angularmomentum(x::Symbol, args...) = angularmomentum(Val(x), args...)

angularmomentum(::Val, ::Orbital"*s", ::Orbital"*s") = 0.0
angularmomentum(::Val, ::Orbital"*p_{*}", ::Orbital"*s") = 0.0
angularmomentum(::Val, ::Orbital"*s", ::Orbital"*p_{*}") = 0.0
angularmomentum(::Val, ::Orbital"*p_{*}", ::Orbital"*p_{*}") = 0.0
angularmomentum(::Val{:x}, ::Orbital"*p_y", ::Orbital"*p_z") = -im
angularmomentum(::Val{:x}, ::Orbital"*p_z", ::Orbital"*p_y") = im
angularmomentum(::Val{:y}, ::Orbital"*p_z", ::Orbital"*p_x") = -im
angularmomentum(::Val{:y}, ::Orbital"*p_x", ::Orbital"*p_z") = im
angularmomentum(::Val{:z}, ::Orbital"*p_x", ::Orbital"*p_y") = -im
angularmomentum(::Val{:z}, ::Orbital"*p_y", ::Orbital"*p_x") = im

function angularmomentum(val::Val, atom::Atom)
    N = countorbitals(atom)
    U = localbasis(atom)
    L = fill(0.0im, N, N)
    for (j, s2) in enumerate(atom.orbitals), (i, s1) in enumerate(atom.orbitals)
        L[i, j] = angularmomentum(val, s1, s2)
    end
    return U'L*U
end

function angularmomentum(val::Val, mol::Molecule)
    N = countorbitals(mol)
    L = fill(0.0im, N, N)

    i = 1
    for atom in mol.atoms
        Ni = countorbitals(atom)
        L[i:i+Ni-1, i:i+Ni-1] = angularmomentum(val, atom)
        i += Ni
    end

    return L
end

function angularmomentum(val::Val, mol::Molecule, atom::Atom)
    N = countorbitals(mol)
    L = fill(0.0im, N, N)
    inds = indices(mol, atom)
    L[inds, inds] = angularmomentum(val, atom)
    return L
end

function angularmomentum(val::Val, mol::Molecule, n::Int)
    N = countorbitals(mol)
    L = fill(0.0im, N, N)
    inds = indices(mol, mol[n])
    L[inds, inds] = angularmomentum(val, mol[n])
    return L
end

##################################################
# Symmetries
##################################################

abstract type AbstractSymmetry end
struct Mirror{S} <: AbstractSymmetry end

parity(::Mirror{:X}, ::Orbital"*s") = 1.0
parity(::Mirror{:X}, ::Orbital"*p_*") = 1.0
parity(::Mirror{:X}, ::Orbital"*p_x") = -1.0

parity(::Mirror{:Y}, ::Orbital"*s") = 1.0
parity(::Mirror{:Y}, ::Orbital"*p_*") = 1.0
parity(::Mirror{:Y}, ::Orbital"*p_y") = -1.0

parity(::Mirror{:Z}, ::Orbital"*s") = 1.0
parity(::Mirror{:Z}, ::Orbital"*p_*") = 1.0
parity(::Mirror{:Z}, ::Orbital"*p_z") = -1.0

(s::Mirror)(a::T) where {T <: Atom} = T(s(a.position), a.orbitals...; valence=a.valence, axes=a.axes)

(s::Mirror{:X})(v::T) where {T <: AbstractVector} = convert(T, [-v[1], v[2], v[3]])

##################################################
# Distribution functions
##################################################

# In atomic units, room temperature is 300K / 315774K = 9.5e-4
nB(E, T = 9.5e-4) = 1/(exp(E/T) - 1)
nF(E, T = 9.5e-4) = 1/(exp(E/T) + 1)

##################################################
# Vibrations
##################################################

"""
    vibrationalcoupling(mol, atom, ω0, dr)

Calculate the numerical derivative of the `hamiltonian(mol)` w.r.t. a
displacement along `dr` of the position of `atom`.  `ω0` is the
characteristic frequency of the atomic vibration.

"""
function vibrationalcoupling(mol::Molecule, atom::Atom, ω0, dr)
    m = Molecules.mass(atom)
    H0 = hamiltonian(mol)
    atom.position += dr
    H′ = hamiltonian(mol)
    atom.position -= dr
    return (H′ - H0) / norm(dr) * √(1 / (2 * m * ω0))
end

##################################################
# Plots recipes
##################################################

include("recipe.jl")
include("blender.jl")
include("xyz.jl")

end #module


text:
using Transformers
using Transformers.Basic

# BertAbs transformer model for abstractive summarization
# from the "Text Summarization with Pretrained Encoders" paper
# by Liu et al. (2019) as described on page 6.
function BertAbs(
    bert_model::TransformerModel{<:CompositeEmbedding,<:Bert,<:Any},
    vocab_size::Integer
)::Translator
    return Translator(
        TransformersModel(
            AbsEmbed(768, vocab_size),
            bert_model.transformers,
            Decoder(768, 8, 96, 2048, 6, pdrop=0.1),
            Generator(768, vocab_size)
        ),
        AbsSearch()
    )
    return
end

# TransformerAbs transformer model for abstractive summarization
# from the "Text Summarization with Pretrained Encoders" paper
# by Liu et al. (2019) as described on pages 6-7.
function TransformerAbs(vocab_size::Integer)::Translator
    return Translator(
        TransformersModel(
            AbsEmbed(768, vocab_size),
            # FIXME The paper mentions a hidden size of 512.
            768, 8, 96, 2048, 6,
            vocab_size,
            pdrop=0.1,
        ),
        AbsSearch()
    )
end

# Tiny version of `TransformerAbs` for testing training on devices with
# less GPU memory.
function TransformerAbsTiny(vocab_size::Integer)::Translator
    return Translator(
        TransformersModel(
            AbsEmbed(8, vocab_size),
            8, 4, 2, 16, 2,
            vocab_size,
            pdrop=0.1,
        ),
        AbsSearch()
    )
end

# Beam search used in the abstractive summarization models
# from the "Text Summarization with Pretrained Encoders" paper
# by Liu et al. (2019) as described on page 6.
function AbsSearch()
    return BeamSearch(
        5,
        # This parameter is not described in the paper.
        32,
        # The length normalization parameter α should be
        # tuned on the development set from 0.6 to 1.0.
        length_normalization=0.8
    )
end

# Word embeddings used in the abstractive summarization models
# from the "Text Summarization with Pretrained Encoders" paper
# by Liu et al. (2019) as described on page 4.
# TODO It's not clear how the random embedding is
# "added to" BERT embeddings.
function AbsEmbed(size::Integer, vocab_size::Integer)
    return WordPositionEmbed(
        size,
        vocab_size,
        # Not mentioned in the paper, but CNN and Daily Mail
        # documents are limited to 2000 tokens, XSum documents
        # are shorter then CNN / Daily Mail on average.
        # Our guess is that no document would be longer
        # than 4096 tokens.
        4096,
        trainable=true
    )
end


text:
# Autogenerated wrapper script for ZITSOL_1_jll for i686-w64-mingw32-libgfortran4
export libZITSOL_1

using LAPACK_jll
using OpenBLAS_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("ZITSOL_1")
JLLWrappers.@declare_library_product(libZITSOL_1, "libZITSOL_1.dll")
function __init__()
    JLLWrappers.@generate_init_header(LAPACK_jll, OpenBLAS_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libZITSOL_1,
        "bin\\libZITSOL_1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    LinearKernel(; c = 0.0)

The linear kernel is a Mercer kernel given by
```
    κ(x,y) = xᵀy + c
```
Where `c` is a real number
"""
struct LinearKernel{Tc<:Real} <: SimpleKernel
    c::Vector{Tc}
    function LinearKernel(;c::T=0.0) where {T}
        new{T}([c])
    end
end

kappa(κ::LinearKernel, xᵀy::Real) = xᵀy + first(κ.c)

metric(::LinearKernel) = DotProduct()

Base.show(io::IO, κ::LinearKernel) = print(io, "Linear Kernel (c = ", first(κ.c), ")")

"""
    PolynomialKernel(; d = 2.0, c = 0.0)

The polynomial kernel is a Mercer kernel given by
```
    κ(x,y) = (xᵀy + c)^d
```
Where `c` is a real number, and `d` is a shape parameter bigger than 1. For `d = 1` see [`LinearKernel`](@ref)
"""
struct PolynomialKernel{Td<:Real, Tc<:Real} <: SimpleKernel
    d::Vector{Td}
    c::Vector{Tc}
    function PolynomialKernel(; d::Td=2.0, c::Tc=0.0) where {Td<:Real, Tc<:Real}
        @check_args(PolynomialKernel, d, d >= one(Td), "d >= 1")
        return new{Td, Tc}([d], [c])
    end
end

kappa(κ::PolynomialKernel, xᵀy::Real) = (xᵀy + first(κ.c))^(first(κ.d))

metric(::PolynomialKernel) = DotProduct()

Base.show(io::IO, κ::PolynomialKernel) = print(io, "Polynomial Kernel (c = ", first(κ.c), ", d = ", first(κ.d), ")")


text:
@testset "periodic_transform" begin
    @testset "compare to periodic exponentiated quadratic" begin
        rng = MersenneTwister(123456)
        f = rand(rng) + 2.0
        x = collect(range(0.0, 3.0 / f; length=1_000))

        # Construct in the usual way.
        k_eq_periodic = PeriodicKernel(; r=[sqrt(0.25)]) ∘ ScaleTransform(f)

        # Construct using the peridic transform.
        k_eq_transform = SqExponentialKernel() ∘ PeriodicTransform(f)

        @test kernelmatrix(k_eq_periodic, x) ≈ kernelmatrix(k_eq_transform, x)
        # TODO - add interface_tests once #159 is merged.
    end
end


text:
@testset "Testing abrkga.jl" begin

    @testset "Testing abrkga.jl with toi file" begin

        instance = read_dataset("datasets/movie.toi")
        objective, bucket_order, statistics = execute_abrkga(instance,seed=0,local_search=true,debug=false)

        @test objective == 72
        @test bucket_order == " 1 | 2 , 3 | 4 | 5 , 6 "
        @test statistics.time_best < statistics.time_total

    end;

end

max_stars_repo_name: LuizHNLorena/AbrkgaObop

max_issues_repo_name: LuizHNLorena/AbrkgaObop

max_forks_repo_name: LuizHNLorena/AbrkgaObop


text:
export PassManager,
       add!, dispose

# subtypes are expected to have a 'ref::API.LLVMPassManagerRef' field
abstract type PassManager end

Base.unsafe_convert(::Type{API.LLVMPassManagerRef}, pm::PassManager) = pm.ref

function add!(pm::PassManager, pass::Pass)
    push!(pm.roots, pass)
    API.LLVMAddPass(pm, pass)
end

dispose(pm::PassManager) = API.LLVMDisposePassManager(pm)

#
# Module pass manager
#

export ModulePassManager, run!

@checked struct ModulePassManager <: PassManager
    ref::API.LLVMPassManagerRef
    roots::Vector{Any}
end

ModulePassManager() = ModulePassManager(API.LLVMCreatePassManager(), [])

function ModulePassManager(f::Core.Function, args...; kwargs...)
    mpm = ModulePassManager(args...; kwargs...)
    try
        f(mpm)
    finally
        dispose(mpm)
    end
end

run!(mpm::ModulePassManager, mod::Module) =
    convert(Core.Bool, API.LLVMRunPassManager(mpm, mod))

#
# Function pass manager
#

export FunctionPassManager,
       initialize!, finalize!, run!

@checked struct FunctionPassManager <: PassManager
    ref::API.LLVMPassManagerRef
    roots::Vector{Any}
end

FunctionPassManager(mod::Module) =
    FunctionPassManager(API.LLVMCreateFunctionPassManagerForModule(mod), [])

function FunctionPassManager(f::Core.Function, args...; kwargs...)
    fpm = FunctionPassManager(args...; kwargs...)
    try
        f(fpm)
    finally
        dispose(fpm)
    end
end

initialize!(fpm::FunctionPassManager) =
    convert(Core.Bool, API.LLVMInitializeFunctionPassManager(fpm))
finalize!(fpm::FunctionPassManager) =
    convert(Core.Bool, API.LLVMFinalizeFunctionPassManager(fpm))

run!(fpm::FunctionPassManager, f::Function) =
    convert(Core.Bool, API.LLVMRunFunctionPassManager(fpm, f))


text:
# Given two particle data sets collected on the same sample but at different
# positions and rotations, find the offset and rotation which best aligns
# datasets.
using Random
using GeometryBasics: Rect2
using StaticArrays
using Rotations
using NearestNeighbors
using LinearAlgebra
using StatsBase
using Distributions
using CoordinateTransformations
using LsqFit
using InvertedIndices: Not

"""
generate_ground_truth(rnd::Random, bounds::Rect2, npart::Integer, inside::Function=p->p in bounds)

Generate a simulated ground-truth particle data set.  This represents the "true" location of `npart` particles that are randomly
distributed within `bounds` and for which `inside(pt)` evaluates true.
"""
function generate_ground_truth(rnd::AbstractRNG, bounds::Rect2, npart::Integer, inside::Function = p -> p in bounds)
    @assert minimum(bounds.widths) > 0.0 "Rectangle widths must be positive."
    sc = maximum(bounds.widths)
    res = map(Base.OneTo(npart)) do _
        pt = SVector{2}(sc .* rand(rnd, Float64, 2) .+ bounds.origin)
        while !inside(pt)
            pt = SVector{2}(sc .* rand(rnd, Float64, 2) .+ bounds.origin)
        end
        pt
    end
    return res
end

"""
    measure_particles(rnd::AbstractRNG, pd::Vector{<:SVector{2}}, offset::SVector{2}, θ::Float64, frac::Float64, eps::Float64)

Sample a ground-truth particle data set in `pd`.  A randomized fraction `frac` of the data set `pd` is sampled and then the
positions are transformed according to an `offset` and a rotation by `θ`. In addition, a measurement jitter the scale of
which is determiend by `eps` is added to each measured coordinate.
"""
function measure_particles(rnd::AbstractRNG, pd::AbstractVector{<:SVector{2}}, offset::SVector{2}, θ::Float64, frac::Float64, eps::Float64)
    @assert frac > 0.0 "The fraction measured `frac` must be larger than 0.0."
    @assert frac <= 1.0 "The fraction measured `frac` must be less than or equal to 1.0."
    rot = RotMatrix{2}(θ)
    n = Distributions.Normal(0.0, eps)
    return shuffle!(rnd, map(filter(f -> rand(rnd) < frac, pd)) do part
        rot * part + offset + rand(rnd, n, 2)
    end)
end

"""
    rough_align(ps1::AbstractVector{StaticVector{2, T}}, ps2::AbstractVector{StaticVector{2, T}}; groupsize=3, tol=0.001, corrtol=10.0) where { T<: AbstractFloat }

Perform a rough alignment of two sets of measured particle coordinates.  The assumption is that the coordinates come from two difference measurements
of the same particle data set which may vary in sample orientation and offset.  The task is to find the affine transformations that
will transform ps1 and ps2 to overlap at the origin.

The particle data sets should include "most" of the same particles or at least a significant overlapping region.  The distance metric should be the
same in both x & y and also between data sets.

The algorithm uses a nearest neighbor method to identify triplets (or groupsize-lets) of particles.  The triplets are matched from `ps1` with the closest
match in `ps2`.  Then the distances between matched triplets is compared to identify triplets that appear similar in particle spacing and in spacing
between center-of-gravities.  These pairs are then used to compute the rotation angle and the center of gravity for each data set.

  * `groupsize` is the number of nearest neighbors to consider (nominally 3)
  * `tol` is the similarity tolerance for separation edge length (<<1.0). `tol` should be approximately the uncertainty in each
  component of the measured positions.

A description of the algorithm is available [here](https://docs.google.com/presentation/d/1XvkeB1AY75s0X-f992S1Qf2H895_9MoI/edit?usp=sharing&ouid=113647574378052266554&rtpof=true&sd=true).

Example:

    julia> (ct1, ct2) = rough_align(ps1, ps2)
    julia> ct1.(ps1) # Translates ps1 towards the origin (no rotation)
    julia> ct2.(ps2) # Translates ps2 towards the origin and rotates it to overlay `ct1.(ps1)`
    julia> (inv(ct1)∘ct2).(ps2) # Transforms ps2 to overlay ps1
    julia> (inv(ct2)∘ct1).(ps1) # Transforms ps1 to overlay ps2
"""
function rough_align(
    ps1::AbstractVector{<:StaticVector{2,T}}, #
    ps2::AbstractVector{<:StaticVector{2,T}}; #
    groupsize = 3, #
    tol = 0.001 #
)::NTuple{2,AffineMap} where {T<:AbstractFloat}
    # For each data set find the (groupsize-1) nearest neighbors for each particle
    idx1, _ = knn(KDTree(ps1), ps1, groupsize, true)
    idx2, _ = knn(KDTree(ps2), ps2, groupsize, true)
    # Compute the separations between all pairs of nearest neighbors
    function compute_separations(idxs, pts)
        # add a sign to indicate clockwise vs counter-clockwise
        signednorm(a, b)::T = (a[1] * b[2] > a[2] * b[1] ? -one(T) : one(T)) * norm(b)
        map(idxs) do idx
            # DOF = groupsize particles times 2 dimensions - (rotate, x_off, y_off)
            dnn12 = pts[idx[1]] - pts[idx[2]]
            SA[collect(signednorm(dnn12, pts[idx[i]] - pts[idx[j]]) for i in Base.OneTo(groupsize) for j in i+1:groupsize)...]
        end # separations
    end
    # Transform the particle data into "nearest neighbor separation space"
    si1 = compute_separations(idx1, ps1)
    si2 = compute_separations(idx2, ps2)
    # Match particle groupings between ps1 and ps2 by matching the ones that are most similar in shape and size
    correspondences = collect(zip(eachindex(idx1), nn(KDTree(si2), si1)...)) # index in idx1, index in idx2, distance
    # Order the correspondences by similarity
    sort!(correspondences, lt = (a, b) -> a[3] < b[3])
    # These functions return the equivalent i-th triplet of particle coordinates in either sp1 or sp2
    triplet1(i::Integer) = ps1[idx1[correspondences[i][1]]]
    triplet2(i::Integer) = ps2[idx2[correspondences[i][2]]]
    # Compute the "center-of-gravity" for a series of coordinates
    cog(coords) = mean(coords)
    # Now perform a second filter.  Look at the center-of-gravity for the matching triplets.
    # Build a list of pairs for which the distance between the cog for triplet1(i) and triplet2(i) are similar.
    res = Tuple{Int,Int,Float64,Float64}[] # index into tripletX(), index into tripletX(), difference in separation, angle between
    clast = min(length(correspondences), 100)
    for ii in 1:clast, jj in ii+1:clast
        dc1, dc2 = cog(triplet1(ii)) - cog(triplet1(jj)), cog(triplet2(ii)) - cog(triplet2(jj))
        ndc1, ndc2 = norm(dc1), norm(dc2)
        # Use determinant to diffentiate 0 to π from π to 2π
        dets = dc1[1] * dc2[2] > dc1[2] * dc2[1] ? -one(T) : one(T)
        # if the norms are exactly equal to zero then they likely represent the same triplets but starting with a different seed particle
        if abs(ndc1 - ndc2) < 6 * tol && ((ndc1 > 1.0e-12) || (ndc2 > 1.0e-12))
            push!(res, (ii, jj, abs(ndc1 - ndc2), mod(dets * acos(clamp(dc1 ⋅ dc2 / (ndc1 * ndc2), -1.0, 1.0)), 2π)))
        end
    end
    # res now contains the indices into ps1 & ps2 for triplets that match in shape and separation
    sort!(res, lt = (a, b) -> a[3] < b[3]) # sort by difference in length
    # Handles cyclic boundaries on angle measurement
    function between_angles(ϕ, ϕmin, ϕmax)
        @assert ϕmin <= ϕmax
        ϕ, ϕmin, ϕmax = mod(ϕ, 2π), mod(ϕmin, 2π), mod(ϕmax, 2π)
        return (ϕmin != ϕmax) && (ϕmin < ϕmax ? (ϕ > ϕmin) && (ϕ < ϕmax) : (ϕ > ϕmin) || (ϕ < ϕmax))
    end
    # Construct a histogram to determine the most probable angle...
    bins = 0.0:π/180.0:2π
    h = fit(Histogram, map(r -> r[4], res), bins)
    θi = findmax(h.weights)[2] # most probably rotation angle
    # Find the indices of corresponding pairs of groups close to this angle
    idxs = filter(i -> between_angles(res[i][4], bins[θi] - π / 180, bins[θi+1] + π / 180), eachindex(res))
    @assert length(idxs) > 0 "θi = $(findmax(h.weights))"
    good = res[idxs]
    # Compute the best estimate of the rotation angle
    θ = mean(r -> r[4], good)
    goodidxs = unique(mapreduce(g -> [g[1], g[2]], append!, good)) # g[1] & g[2] are indices of corresponding particle groups
    # Use these indices to compute the best estimates of the center-of-gravity for ps1 and ps2
    cog1 = cog(map(i -> cog(triplet1(i)), goodidxs))
    cog2 = cog(map(i -> cog(triplet2(i)), goodidxs))
    # Convert this into two affine transformations.
    return (
        LinearMap(Angle2d(zero(T))) ∘ inv(Translation(cog1)),
        LinearMap(Angle2d(T(θ))) ∘ inv(Translation(cog2))
    )
end

"""
    correspondences(ps1::AbstractVector{<:StaticVector{2, T}}, ps2::AbstractVector{<:StaticVector{2, T}}; tol=0.01, invert=false) where { T<: AbstractFloat }

Returns a pair of indices `(idx1, idx2)` into `ps1` and `ps2` respectively such that `ps1[idx1]` are likely the same particle as `ps2[idx2]`.

`invert=true` inverts the meaning so returns particles which don't match well with another.
"""
function correspondences(ps1::AbstractVector{<:StaticVector{2,T}}, ps2::AbstractVector{<:StaticVector{2,T}}; tol = 0.01, invert = false) where {T<:AbstractFloat}
    (idx1, _) = nn(KDTree(ps2), ps1) # idx1 is an index into ps2
    (idx2, dist2) = nn(KDTree(ps1), ps2) # idx2 is an index into ps1
    keep2 = map(enumerate(idx2)) do (i2, i1)
        # Within tolerance and closest alternative in either direction
        (dist2[i2] < tol) && (i2==idx1[i1])
    end
    c1 = idx2[keep2]
    c2 = idx1[c1]
    return invert ? ( Not(c1), Not(c2) ) : ( c1, c2 )
end

function orthogonal_procrustes_alignment(ps1::AbstractVector{<:StaticVector{2,T}}, ps2::AbstractVector{<:StaticVector{2,T}}; tol = 0.01) where {T<:AbstractFloat}
    c1, c2 = correspondences(ps1, ps2; tol = tol, invert = false)
    cps1, cps2 = ps1[c1], ps2[c2]
    com1, com2 = mean(cps1), mean(cps2)
    cpst1, cpst2 = Translation(-com1).(cps1), Translation(-com2).(cps2)
    m2(vs) = transpose(reshape(collect(Iterators.flatten(vs)), length(vs[1]), length(vs)))
    # Solve the orthogonal Procrustes problem
    f = svd(m2(cpst1) * transpose(m2(cpst2)))
    Ω = f.U * f.Vt
    # Now translate and rotated the reordered data towards `pts1`
    return cps1, Translation(com1).(Ω * cpst2)
end

"""
    refined_alignment(ps1::Vector{<:StaticVector{2,T}}, ps2::Vector{<:StaticVector{2,T}}) where { T <: AbstractFloat }

Uses a non-linear algorithm to refine a rough alignment by adjusting the rotation and offset to produce the smallest
mean-square distance between particles in `ps1` and `ps2`.  Return the `AffineMap` that best transforms `ps2` to match
`ps1`.

The inputs `ps1` and `ps2` must correspond index-by-index with each other. Meaning they represent the exact same particles
in the same order.  This is where the `correspondences(...)` function comes in handy.

Example:

    julia> ct1, ct2 = rough_align(ps1, ps2)                  # Get the alignment transforms
    julia  rs1, rs2 = ct1.(ps1),ct2.(ps2)                    # Rough align the data sets
    julia> ci1, ci2 = correspondences(rs1, rs2)              # Identify corresponding rough-aligned particles by index
    julia> cs1, cs2 = rs1[ci1], rs2[ci2]                     # Extract these particles
    julia> n(ds1, ds2) = sum(x->x^2, norm.(ds1 .- ds2))       # Metric function
    julia> n(rs1[ci1], rs2[ci2])                             # Pre-Error measure
    julia> ct3 = refined_alignment(cs1, cs2)                 # Get the refined alignment transform
    julia> n(rs1, ct3.(rs2))                                 # Post-Error measure
    julia> n(rs1, (ct3∘ct2).(cs2));
    julia> aps1, aps2 = ct1.(ps1), (ct3∘ct2).(ps2)           # Transform all points
"""
function refined_alignment(ps1::AbstractVector{<:StaticVector{2,T}}, ps2::AbstractVector{<:StaticVector{2,T}}) where {T<:AbstractFloat}
    @assert length(ps1) == length(ps2)
    function flatten(ps)
        res = Array{T}(undef, 2 * length(ps))
        foreach(i -> res[2i-1:2i] .= ps[i], eachindex(ps))
        return res
    end
    # Compute the function
    function f(res, ps, param)
        A, B, C, D = cos(param[1]), sin(param[1]), param[2], param[3]
        for i in 1:2:length(ps)
            x, y = ps[i], ps[i+1]
            res[i] = A * (x + C) - B * (y + D)
            res[i+1] = B * (x + C) + A * (y + D)
        end
    end
    # Compute the Jacobian
    function jac(res, ps, param)
        A, B, C, D = cos(param[1]), sin(param[1]), param[2], param[3]
        for i in 1:2:length(ps)
            x, y = ps[i], ps[i+1]
            res[i, :] .= (-B * (x + C) - A * (y + D), A, -B)
            res[i+1, :] .= (-B * (y + D) + (x + C) * A, B, A)
        end
        return res
    end
    fit = curve_fit(f, jac, flatten(ps2), flatten(ps1), [0.0, 0.0, 0.0]; inplace = true)
    return LinearMap(Angle2d(fit.param[1])) ∘ Translation(fit.param[2:3])
end

"""
    align(ps1::Vector{<:StaticVector{2,T}}, ps2::Vector{<:StaticVector{2,T}}; tol=0.001, finealign=true)  where { T <: AbstractFloat }

Perform a `rough_align(...)` followed by a "refined_alignment(...)" on `ps1` and `ps2`.  Return the `AffineMap`
transformations `(ct1, ct2) such that `ct1.(ps1)` and `ct2.(ps2)` are registered.  If `finealign=true` then a
second-stage non-linear optimization of corresponding points is performed.

Example:

    julia> (ct1, ct2) = align(ps1, ps2)
    julia> ct1.(ps1) # Translates ps1 towards the origin (no rotation)
    julia> ct2.(ps2) # Translates ps2 towards the origin and rotates it to overlay `ct1.(ps1)`
    julia> (inv(ct1)∘ct2).(ps2) # Transforms ps2 to overlay ps1
    julia> (inv(ct2)∘ct1).(ps1) # Transforms ps1 to overlay ps2
"""
function align(ps1::AbstractVector{<:StaticVector{2,T}}, ps2::AbstractVector{<:StaticVector{2,T}}; tol = 0.001, finealign = true) where {T<:AbstractFloat}
    ct1, ct2 = rough_align(ps1, ps2; tol = tol)
    if finealign
        rs1, rs2 = ct1.(ps1), ct2.(ps2)
        ci1, ci2 = correspondences(rs1, rs2)
        cs1, cs2 = rs1[ci1], rs2[ci2]
        ct3 = refined_alignment(cs1, cs2)
        ct2 = ct3 ∘ ct2
    end
    return (ct1, ct2)
end

"""
    identify(pss::AbstractVector{<:AbstractVector{<:StaticVector{2,T}}}; tol=0.001, ctol=0.01 )::DataFrame where { T <: AbstractFloat}

Takes a vector of X-Y coordinate vectors.  Aligns the X-Y coordinate vectors to the first one.  Then it labels each unique particle and
identifies the index (if any) in which it appears in each data set.  Thus we are able to track particle identities through the data sets.
Each row in the returned `DataFrame` represents a "unique particle".  The `:PS?` columns represent the index of the particle in that
X-Y coordinate vectors. The `:COUNT` column represents the number of data sets in which the particle was measured.  The `:FIRST` column
indicates the first data set in which the particle was observed.
"""
function identify(pss::AbstractVector{<:AbstractVector{<:StaticVector{2,T}}}; tol=0.001, ctol=0.01 )::DataFrame where { T <: AbstractFloat}
    # Align all pss relative to pss[1]
    als = map(pss) do psj
        if psj==pss[1]
            pss[1]
        else
            cti, ctj = align(pss[1], psj, tol=tol, finealign=true)
            (inv(cti)∘ctj).(psj)
        end
    end
    df = DataFrame( (Symbol("PS$i")=> Union{Missing,Int}[] for i in eachindex(als))...)
    for i in eachindex(als)
        alsi = als[i]
        # Determine the indices of the previously unseen particles in als[i]
        unseen = fill(true, length(alsi))
        foreach(i-> (!ismissing(i)) && (unseen[i]=false), df[:,i])
        idxs=collect(eachindex(alsi)[unseen])
        # Get the coordinates of these particle
        uni=collect(alsi[idxs])
        # Find correspondences for these particles with later data sets
        cijs = map(i+1:length(als)) do j
            ci, cj = correspondences(uni, als[j], tol=ctol)
            tmp = zeros(Int, length(uni))
            tmp[ci] .= cj
            tmp
        end
        ztom(x) = x==zero(typeof(x)) ? missing : x
        for (ii, ind) in enumerate(idxs)
            # uni[ii]=alsi[ind]
            push!(df, [ ( missing for jj in 1:i-1)..., ind, ( ztom(cij[ii]) for cij in cijs )...])
        end
    end
    f = CategoricalArray(map(eachrow(df)) do r
        findfirst(c->!ismissing(c), Tuple(r))
    end, ordered=true)
    c = CategoricalArray(map(eachrow(df)) do r
        count(c->!ismissing(c), Tuple(r))
    end, ordered=true)
    x = map(eachrow(df)) do r
        mean(skipmissing(map(i->ismissing(r[i]) ? missing : als[i][r[i]][1], 1:length(r))))
    end
    y = map(eachrow(df)) do r
        mean(skipmissing(map(i->ismissing(r[i]) ? missing : als[i][r[i]][2], 1:length(r))))
    end
    insertcols!(df, :APPEARS => f, :COUNT => c, :XABS => x, :YABS => y )
    sort!(df, [ :APPEARS, order(:COUNT,rev=true) ], alg=MergeSort) # Maintains order
    df
end

max_stars_repo_name: NicholasWMRitchie/NeXLParticle

max_issues_repo_name: NicholasWMRitchie/NeXLParticle

max_forks_repo_name: NicholasWMRitchie/NeXLParticle


text:
using LinearAlgebra
# using SpecialFunctions

"""
    isdependent returns true if eq is dependent on x
"""
isdependent(eq, x) = !isequal(expand_derivatives(Differential(x)(eq)), 0)

# this is the main heurisctic used to find the test fragments
function generate_basis(eq, x)
    Δeq = expand_derivatives(Differential(x)(eq))
    kers = expand(eq + Δeq)
    return [one(x); candidates(kers, x)]
end

"""
    candidates returns a list of candidate expressions to form the integration
    basis
"""
candidates(eq, x) = isdependent(eq,x) ? [eq] : []

candidates(eq::Num, x) = candidates(value(eq), x)

# the candidates of an Add is the union of the candidates of the terms
# ∫ Σᵢ fᵢ(x) dx = Σᵢ ∫ fᵢ(x) dx
candidates(eq::SymbolicUtils.Add, x) = unique(∪([candidates(t,x) for t in arguments(eq)]...))

# the candidates of a Mul is the outer product of the candidates of the terms
# d(uv)/dx = u dv/dx + v + du/dx
function candidates(eq::SymbolicUtils.Mul, x)
    terms = [candidates(q,x) for q in arguments(eq)]
    n = length(terms)

    l = Any[one(x)]

    for j = 1:n
        m = length(l)
        for t in terms[j]
            for k = 1:m
                push!(l, l[k]*t)
            end
        end
    end

    unique(l[2:end])    # removing the initial 1
end

# the candidates of a Pow encode different integration rules
function candidates(eq::SymbolicUtils.Pow, x)
    if !isdependent(eq,x) return [one(x)] end

    p = arguments(eq)[1]    # eq = p ^ k
    k = arguments(eq)[2]

    if k == -1
        return candidate_pow_minus_one(p, x)
    elseif k ≈ 0.5 || k ≈ -0.5
        # ∫ √f(x) dx = ... + c * log(df/dx + √f) if deg(f) == 2
        Δ = expand_derivatives(Differential(x)(p))
        return [[p^k, p^(k+1)]; log(abs(0.5*Δ + sqrt(p)))]
    end

    # ∫ p^k dp = c * p^(k+1)
    return [p^k, p^(k+1)]
end

function candidate_pow_minus_one(p, x)
    q = to_poly(p, x)   # q is a DynamicPolynomials Polynomial
                        # the reason is to be able to use find_roots
    if q == nothing return [p, p^-1, log(p)] end

    q = 0 + 1*q         # take care of monomials and constants
    r, s = find_roots(q, var(q))
    s = s[1:2:end]
    r = nice_parameters(r)
    s = Complex.(nice_parameters(real.(s)), nice_parameters(imag(s)))

    # ∫ 1 / ((x-z₁)(x-z₂)) dx = ... + c₁ * log(x-z₁) + c₂ * log(x-z₂)
    q = sum(log(x - u) for u in r; init=zero(x)) +
        sum(atan((x - real(u))/imag(u)) for u in s; init=zero(x)) +
        sum(log(x^2 - 2*real(u)*x + abs2(u)) for u in s; init=zero(x))

    return [[p, p^-1]; candidates(q, x)]
end

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

"""
    integrate is the main entry point

    input:
    ------
    eq: a Symbolics expression to integrate
    abstol: the desired tolerance
    num_steps: the number of different steps with expanding basis to be tried
    num_trials: the number of trials in each step (no changes to the basis)
    lo and hi: the range used to generate random values of x (the independent variable)
    show_basis: if true, the basis is printed

    output:
    -------
    solved, unsolved

    a pair of expressions, solved is the solved integral and unsolved is the residual unsolved
    portion of the input
"""
function integrate(eq; abstol=1e-6, num_steps=3, num_trials=7, lo=-5.0, hi=5.0, show_basis=false)
    x = var(eq)
    if x == nothing
        @syms 𝑥
        return 𝑥 * eq, 0
    end
    integrate(eq, x; abstol, num_trials, num_steps, lo, hi, show_basis)
end

function integrate(eq::SymbolicUtils.Add, x; abstol=1e-6, num_steps=3, num_trials=7, lo=-5.0, hi=5.0, show_basis=false)
    solved = 0
    unsolved = 0

    for p in arguments(eq)
        s, u = integrate(p, x; abstol, num_steps, num_trials, lo, hi, show_basis)
        solved += s
        unsolved += u
    end

    solved, unsolved
end

function integrate(eq, x; abstol=1e-6, num_steps=3, num_trials=7, lo=-5.0, hi=5.0, show_basis=false)
    eq₁ = apply_integration_rules(expand(eq))
    basis = generate_basis(eq₁, x)

    if show_basis println(basis) end

    D = Differential(x)

    for i = 1:num_steps
        basis = unique([basis; basis*x])
        Δbasis = [expand_derivatives(D(f)) for f in basis]

        for j = 1:num_trials
            y, ϵ = try_integrate(eq, x, basis, Δbasis; abstol, lo, hi)
            if ϵ < abstol return y, 0 end
        end
    end
    # @warn "no solution is found"
    0, eq
end

rms(x) = sqrt(sum(x.^2) / length(x))

"""
    the core of the randomized parameter-fitting algorithm

    `try_integrate` tries to find a linear combination of the basis, whose
    derivative is equal to eq

    output
    -------
    integral, error
"""
function try_integrate(eq, x, basis, Δbasis; abstol=1e-6, lo=-5.0, hi=5.0, attemp_ratio=5)
    n = length(basis)
    # A is an nxn matrix holding the values of the fragments at n random points
    # b hold the value of the input function at those points
    A = zeros(n, n)
    b = zeros(n)

    i = 1
    k = 1

    while i <= n
        x₀ = rand()*(hi - lo) + lo
        d = Dict(x => x₀)
        try
            for j = 1:n
                A[i, j] = Float64(substitute(Δbasis[j], d))
            end
            b[i] = Float64(substitute(eq, d))
            i += 1
        catch e
            # println("exclude ", x₀)
        end
        if k > attemp_ratio*n return nothing, 1e6 end
        k += 1
    end

    # find a linearly independent subset of the basis
    l = find_independent_subset(A; abstol)
    A, b, basis = A[l,l], b[l], basis[l]

    if det(A) ≈ 0 return nothing, 1e6 end

    q₀ = A \ b
    q = nice_parameters(q₀)
    ϵ = rms(A * q - b)
    sum(q[i]*basis[i] for i = 1:length(basis) if q[i] != 0; init=zero(x)), ϵ
end

"""
    returns a list of the indices of a linearly independent subset of the columns of A
"""
function find_independent_subset(A; abstol=1e-5)
    _, R = qr(A)
    abs.(diag(R)) .> abstol
end

"""
    converts float to int or small rational numbers
"""
function nice_parameters(p; abstol=1e-3)
    c = lcm(collect(1:10)...)
    n = length(p)
    q = Array{Any}(undef, n)
    for i = 1:n
        den = 1
        while den < 10
            if abs(round(p[i]*den) - p[i]*den) < abstol
                a = round(Int, p[i]*den) // den
                q[i] = (denominator(a) == 1 ? numerator(a) : a)
                den = 10
            else
                q[i] = p[i]
            end
            den += 1
        end
    end
    q
end

########################## Transformation Rules ###############################

trig_rule1 = @rule tan(~x) => sin(~x) / cos(~x)
trig_rule2 = @rule sec(~x) => one(~x) / cos(~x)
trig_rule3 = @rule csc(~x) => one(~x) / sin(~x)
trig_rule4 = @rule cot(~x) => cos(~x) / sin(~x)

trig_rules = [trig_rule1, trig_rule2, trig_rule3, trig_rule4]

hyper_rule1 = @rule tanh(~x) => sinh(~x) / cosh(~x)
hyper_rule2 = @rule sech(~x) => one(~x) / cosh(~x)
hyper_rule3 = @rule csch(~x) => one(~x) / sinh(~x)
hyper_rule4 = @rule coth(~x) => cosh(~x) / sinh(~x)

hyper_rules = [hyper_rule1, hyper_rule2, hyper_rule3, hyper_rule4]

misc_rule1 = @rule sqrt(~x) => ^(~x, 0.5)

misc_rules = [misc_rule1]

int_rules = [trig_rules; hyper_rules; misc_rules]

apply_integration_rules(eq) = Fixpoint(Prewalk(PassThrough(Chain(int_rules))))(value(eq))

########################## Convert to a Polynomial? #############################

@polyvar 𝑦

is_number(x::T) where T<:Integer = true
is_number(x::T) where T<:Real = true
is_number(x::T) where T<:Complex = true
is_number(x::T) where T<:Rational = true
is_number(x) = false

function is_poly(eq, x)
    p = collect_powers(value(eq), x)
    all(is_number, values(p))
end

function to_poly(p::SymbolicUtils.Add, x)
    l = [to_poly(t,x) for t in arguments(p)]
    if any(x->x==nothing, l) return nothing end
    sum(l; init=0)
end

function to_poly(p::SymbolicUtils.Mul, x)
    l = [to_poly(t,x) for t in arguments(p)]
    if any(x->x==nothing, l) return nothing end
    prod(l; init=1)
end

function to_poly(p::SymbolicUtils.Pow, x)
    if !isequal(arguments(p)[1], x) return nothing end
    return 𝑦^arguments(p)[2]
end

to_poly(p::SymbolicUtils.Sym, x) = isequal(p, x) ? 𝑦 : nothing
to_poly(p::SymbolicUtils.Term, x) = nothing
to_poly(p, x) = is_number(p) ? p : nothing

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

@syms x

"""
    a list of basic standard integral tests
    based on http://integral-table.com/ with modifications
"""
basic_integrals = [
# Basic Forms
    1,
    x^2,
    4x^3,
# Integrals of Rational Functions
    1/x,
    1/(2x + 5),
    1/(x + 1)^2,
    (x + 3)^3,
    x*(x - 2)^4,
    1/(1 + x^2),
    1/(9 + x^2),
    x/(4 + x^2),
    x^2/(16 + x^2),
    x^3/(1 + x^2),
    1/(x^2 - 5x + 6),
    1/(x^2 + x + 1),
    x/(x + 4)^2,
    x/(x^2 + x + 1),
# Integrals with Roots
    sqrt(x - 2),
    1 / sqrt(x - 1),
    1 / sqrt(x + 1),
    1 / sqrt(4 - x),
    x * sqrt(x - 3),
    sqrt(2x + 5),
    (3x - 1)^1.5,
    x / sqrt(x - 1),
    x / sqrt(x + 1),
    sqrt(x / (4 - x)),
    sqrt(x / (4 + x)),
    x * sqrt(2x + 3),
    sqrt(x*(x+2)),
    sqrt(x^3*(x+3)),
    sqrt(x^2 + 4),
    sqrt(x^2 - 4),
    sqrt(4 - x^2),
    x * sqrt(x^2 + 9),
    x * sqrt(x^2 - 9),
    1 / sqrt(x^2 + 4),
    1 / sqrt(x^2 - 4),
    1 / sqrt(4 - x^2),
    x / sqrt(x^2 + 4),
    x / sqrt(x^2 - 4),
    x / sqrt(4 - x^2),
    x^2 / sqrt(x^2 + 4),
    x^2 / sqrt(x^2 - 4),
    sqrt(x^2 - 5x + 6),
    x * sqrt(x^2 - 5x + 6),
    1 / sqrt(x^2 - 5x + 6),
    1 / (4 + x^2)^1.5,
# Integrals with Logarithms
    log(x),
    x * log(x),
    x^2 * log(x),
    log(2x) / x,
    log(x) / x^2,
    log(2x + 1),
    log(x^2 + 4),
    log(x^2 - 4),
    log(x^2 - 5x + 6),
    x * log(x + 2),
    x * log(9 - 4x^2),
    log(x)^2,
    log(x)^3,
    x * log(x)^2,
    x^2 * log(x)^2,
# Integrals with Exponentials
    exp(x),
    sqrt(x) * exp(x),
    x * exp(x),
    x * exp(3x),
    x^2 * exp(x),
    x^2 * exp(5x),
    x^3 * exp(x),
    x^3 * exp(2x),
    exp(x^2),
    x * exp(x^2),
# Integrals with Trigonometric Functions
    sin(4x),
    sin(x)^2,
    sin(x)^3,
    cos(3x),
    cos(x)^2,
    cos(2x)^3,
    sin(x) * cos(x),
    sin(3x) * cos(5x),
    sin(x)^2 * cos(x),
    sin(3x)^2 * cos(x),
    sin(x) * cos(x)^2,
    sin(x) * cos(5x)^2,
    sin(x)^2 * cos(x),
    sin(x)^2 * cos(x)^2,
    sin(4x)^2 * cos(4x^2),
    tan(x),
    tan(7x),
    tan(x)^2,
    tan(x)^3,
    sec(x),
    sec(x)^2,
    sec(x)^3,
    sec(x) * tan(x),
    sec(x)^2 * tan(x),
    csc(x),
    csc(x)^2,
    csc(x)^3,
    sec(x) * csc(x),
# Products of Trigonometric Functions and Monomials
    x * cos(x),
    x * cos(3x),
    x^2 * cos(x),
    x^2 * cos(5x),
    x * sin(x),
    x * sin(3x),
    x^2 * sin(x),
    x^2 * sin(5x),
    x * cos(x)^2,
    x * sin(x)^2,
    x * tan(x)^2,
    x * sec(x)^2,
# Products of Trigonometric Functions and Exponentials
    exp(x) * sin(x),
    exp(3x) * sin(2x),
    exp(x) * cos(x),
    exp(2x) * cos(7x),
    x * exp(x) * sin(x),
    x * exp(x) * cos(x),
# Integrals of Hyperbolic Functions
    cosh(x),
    exp(x) * cosh(x),
    sinh(3x),
    exp(2x) * sinh(3x),
    tanh(x),
    exp(x) * tanh(x),
    cos(x) * cosh(x),
    cos(x) * sinh(x),
    sin(x) * cosh(x),
    sin(x) * sinh(x),
    sinh(x) * cosh(x),
    sinh(3x) * cosh(5x),
]

function test_integrals()
    for eq in basic_integrals
        printstyled(eq, " =>\t"; color=:green)
        solved, unsolved = integrate(eq)
        printstyled(solved; color=:white)
        if isequal(unsolved, 0)
            println()
        else
            printstyled(" + ∫ ", unsolved, '\n'; color=:red)
        end
    end
end

##################### Special Functions ######################################

# """
#     logarithmic integral
# """
# function li(x; n=10)
#     z = log(abs(x))
#     s = sum(z^k / (factorial(k) * k) for k = 1:n)
#     return SpecialFunctions.γ + log(z) + s
# end


text:
@testset "Non_negative" begin

    p, p_inv = Non_negative()

    x = abs.(randn((10, 10)))
    x2 = 100 .*randn((10, 10))

    @test x ≈ p(p_inv(x))
    @test x ≈ p_inv(p(x))

    @test all(p(x2) .>= 0)
end

@testset "Map_0_1" begin

    p, p_inv = Map_0_1()

    x = abs.(randn((10, 10)))
    x2 = 100 .*randn((10, 10))

    @test x ≈ p_inv(p(x))

    @test all(p(x2) .>= 0)
    @test all(p(x2) .<= 1)
end

@testset "Pow4_positive" begin
    p, p_inv = Pow4_positive()

    x = abs.(randn((10, 10)))
    x2 = 100 .*randn((10, 10))

    @test x ≈ p(p_inv(x))
    @test x ≈ p_inv(p(x))

    @test all(p(x2) .>= 0)
end

@testset "Piecewise_positive" begin
    p, p_inv = Piecewise_positive()

    x = abs.(randn((10, 10)))
    x2 = 100 .*randn((10, 10))

    @test x ≈ p(p_inv(x))
    @test x ≈ p_inv(p(x))

    @test all(p(x2) .>= 0)

    function f(x)
        @test all(.≈(Zygote.gradient(x -> sum(p(x)), x)[1], (p(x .+ 1e-8) .- p(x))./1e-8, rtol=1e-4))
        @test Zygote.gradient(x -> sum(p(x)), x)[1] ≈ DeconvOptim.f_pw_pos_grad(x)
    end

    f([1.1, 12312.2, -10.123, 22.2, -123.23, 0])
end

@testset "Abs_positive" begin
    p, p_inv = Abs_positive()

    x = abs.(randn((10, 10)))
    x2 = 100 .*randn((10, 10))

    @test x ≈ p(p_inv(x))
    @test x ≈ p_inv(p(x))

    @test all(p(x2) .>= 0)
end


text:
ENV["JULIA_DATASETS_PATH"] = ""

using DataSets
using Test
using UUIDs
using ResourceContexts

using DataSets: FileSystemRoot

#-------------------------------------------------------------------------------
@testset "DataSet config" begin
    proj = DataSets.load_project("Data.toml")

    ds = dataset(proj, "a_text_file")
    @test ds.uuid == UUID("b498f769-a7f6-4f67-8d74-40b770398f26")
    @test ds.name == "a_text_file"
    @test ds.description == "A text file"
    @test ds.storage["driver"] == "FileSystem"
end

@testset "DataSet config from Dict" begin
    config = Dict(
        "data_config_version"=>0,
        "datasets"=>[Dict(
            "description"=>"A text file",
            "name"=>"a_text_file",
            "uuid"=>"b498f769-a7f6-4f67-8d74-40b770398f26",

            "storage"=>Dict(
                "driver"=>"FileSystem",
                "type"=>"Blob",
                "path"=>joinpath(@__DIR__, "data", "file.txt")
               )
           )]
       )

    proj = DataSets.load_project(config)

    ds = dataset(proj, "a_text_file")
    @test ds.uuid == UUID("b498f769-a7f6-4f67-8d74-40b770398f26")
end

@testset "open() for DataSet" begin
    proj = DataSets.load_project("Data.toml")

    text_data = dataset(proj, "a_text_file")
    @test open(text_data) isa Blob
    @test read(open(text_data), String) == "Hello world!\n"
    @context begin
        @test read(@!(open(text_data)), String) == "Hello world!\n"
    end

    tree_data = dataset(proj, "a_tree_example")
    @test open(tree_data) isa BlobTree
    @context begin
        @test @!(open(tree_data)) isa BlobTree
        tree = @! open(tree_data)
        @test readdir(tree) == ["1.csv", "2.csv"]
    end
end

#-------------------------------------------------------------------------------
@testset "open() for Blob and BlobTree" begin
    blob = Blob(FileSystemRoot("data/file.txt"))
    @test        open(identity, String, blob)         == "Hello world!\n"
    @test String(open(identity, Vector{UInt8}, blob)) == "Hello world!\n"
    @test open(io->read(io,String), IO, blob)         == "Hello world!\n"
    @test open(identity, Blob, blob) === blob
    # Unscoped forms
    @test open(String, blob)                == "Hello world!\n"
    @test String(open(Vector{UInt8}, blob)) == "Hello world!\n"
    @test read(open(IO, blob), String)      == "Hello world!\n"

    tree = BlobTree(FileSystemRoot("data"))
    @test open(identity, BlobTree, tree) === tree

    # Context-based forms
    @context begin
        @test @!(open(String, blob))               == "Hello world!\n"
        @test String(@! open(Vector{UInt8}, blob)) == "Hello world!\n"
        @test read(@!(open(IO, blob)), String)     == "Hello world!\n"
        @test @!(open(Blob, blob))                 === blob
        @test @!(open(BlobTree, tree))             === tree
    end
end

#-------------------------------------------------------------------------------
@testset "Data set name parsing" begin
    # Valid names
    @test DataSets.check_dataset_name("a_b") === nothing
    @test DataSets.check_dataset_name("a1") === nothing
    @test DataSets.check_dataset_name("δεδομένα") === nothing
    @test DataSets.check_dataset_name("a/b") === nothing
    @test DataSets.check_dataset_name("a/b/c") === nothing
    # Invalid names
    @test_throws ErrorException("DataSet name \"a?b\" is invalid. DataSet names must start with a letter and can contain only letters, numbers, `_` or `/`.") DataSets.check_dataset_name("a?b")
    @test_throws ErrorException DataSets.check_dataset_name("1")
    @test_throws ErrorException DataSets.check_dataset_name("a b")
    @test_throws ErrorException DataSets.check_dataset_name("a.b")
    @test_throws ErrorException DataSets.check_dataset_name("a/b/")
    @test_throws ErrorException DataSets.check_dataset_name("/a/b")
end

#-------------------------------------------------------------------------------
# Trees
@testset "Temporary trees" begin
    function write_dir(j)
        d = newdir()
        for i=1:2
            d["hi_$i.txt"] = newfile() do io
                println(io, "hi $j $i")
            end
        end
        return d
    end

    temptree = newdir()
    for j=1:3
        temptree["d$j"] = write_dir(j)
    end
    @test open(io->read(io,String), IO, temptree["d1"]["hi_2.txt"]) == "hi 1 2\n"
    @test open(io->read(io,String), IO, temptree["d3"]["hi_1.txt"]) == "hi 3 1\n"
    @test isfile(DataSets.sys_abspath(temptree["d1"]["hi_2.txt"]))
end

include("projects.jl")
include("entrypoint.jl")
include("repl.jl")
include("DataTomlStorage.jl")
include("backend_compat.jl")
include("driver_autoload.jl")


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

mutable struct TroubleshootingProperties <: SwaggerModel
    storageId::Any # spec type: Union{ Nothing, String } # spec name: storageId
    storagePath::Any # spec type: Union{ Nothing, String } # spec name: storagePath

    function TroubleshootingProperties(;storageId=nothing, storagePath=nothing)
        o = new()
        validate_property(TroubleshootingProperties, Symbol("storageId"), storageId)
        setfield!(o, Symbol("storageId"), storageId)
        validate_property(TroubleshootingProperties, Symbol("storagePath"), storagePath)
        setfield!(o, Symbol("storagePath"), storagePath)
        o
    end
end # type TroubleshootingProperties

const _property_map_TroubleshootingProperties = Dict{Symbol,Symbol}(Symbol("storageId")=>Symbol("storageId"), Symbol("storagePath")=>Symbol("storagePath"))
const _property_types_TroubleshootingProperties = Dict{Symbol,String}(Symbol("storageId")=>"String", Symbol("storagePath")=>"String")
Base.propertynames(::Type{ TroubleshootingProperties }) = collect(keys(_property_map_TroubleshootingProperties))
Swagger.property_type(::Type{ TroubleshootingProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_TroubleshootingProperties[name]))}
Swagger.field_name(::Type{ TroubleshootingProperties }, property_name::Symbol) =  _property_map_TroubleshootingProperties[property_name]

function check_required(o::TroubleshootingProperties)
    (getproperty(o, Symbol("storageId")) === nothing) && (return false)
    (getproperty(o, Symbol("storagePath")) === nothing) && (return false)
    true
end

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


text:
using Test

@testset "MAFE" begin
    include("types_test.jl")
    include("stata_files.jl")
    include("histories_test.jl")
end

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


text:
export distproxy
function distproxy end

import MonteCarloMeasurements

PROXIES = Dict(
    :Distributions => [
        :mean
        :std
        :entropy
        :cdf
        ],
    :MonteCarloMeasurements => [
        :Particles
    ]
)

for m in keys(PROXIES)
    for f in PROXIES[m]
        @eval begin
            import $m: $f
            export $f
            $m.$f(d::AbstractMeasure, args...) = $m.$f(MeasureTheory.distproxy(d), args...)
        end
    end
end

MonteCarloMeasurements.Particles(N::Int, d::AbstractMeasure) = MonteCarloMeasurements.Particles(N, distproxy(d))

# using MonteCaroMeasurements

# MonteCaroMeasurementsPROXIES = [
#     :Particles
# ]

# for f in DistributionsPROXIES
#     @eval begin
#         import Distributions: $f
#         export $f
#         Distributions.$f(d::AbstractMeasure) = Distributions.$f(MeasureTheory.distproxy(d))
#     end
# end


text:
using DoubleExponentialFormulas: quaddeo
using LinearAlgebra
using Test

# Test functions in
# Ooura, T., Mori, M., 1991.
# The double exponential formula for oscillatory functions over the half infinite interval.
# Journal of Computational and Applied Mathematics 38, 353–360.
# https://doi.org/10.1016/0377-0427(91)90181-I

# 1
let f(x) = exp(-x)*cos(x), expect = 1/2
    I, E = quaddeo(f, 1.0, π/2, 0.0, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

# 2
let f(x) = x*sin(x)/(1 + x^2), expect = π/(2*exp(1))
    I, E = quaddeo(f, 1.0, 0.0, 0.0, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

# 3
let f(x) = cos(x)/(1 + x^2), expect = π/(2*exp(1))
    I, E = quaddeo(f, 1.0, π/2, 0.0, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

# 4
let f(x) = log((x^2 + 4)/(x^2 + 1))*cos(x), expect = (exp(-1) - exp(-2))*π
    I, E = quaddeo(f, 1.0, π/2, 0.0, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

# 5 (Dirichlet integral)
let f(x) = sin(x)/x, expect = π/2
    I, E = quaddeo(f, 1.0, 0.0, 0.0, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

# 6
let f(x) = 1/sqrt(x)*sin(x), expect = sqrt(π/2)
    I, E = quaddeo(f, 1.0, 0.0, 0.0, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

# 7
let f(x) = cos(x)/sqrt(x), expect = sqrt(π/2)
    I, E = quaddeo(f, 1.0, π/2, 0.0, Inf)
    @test I isa Float64
    @test_broken I ≈ expect
    @test abs(I - expect) < 1e-6
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

# 8
let f(x) = sin(x)*log(x), expect = -Base.MathConstants.eulergamma
    I, E = quaddeo(f, 1.0, 0.0, 0.0, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

# Various interval (-∞, ∞)
let f(x) = sin(x)/x, expect = π
    I, E = quaddeo(f, 1.0, 0.0, -Inf, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

# Various interval (0, ∞)
let f(x) = sin(x)/x, expect = π/2
    I, E = quaddeo(f, 1.0, 0.0, 0.0, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

# Various interval (a, ∞)
let f(x) = sin(x)/x, expect = π/2 - BigFloat("1.851937051982466170361053370157991363345809728981154909804")
    I, E = quaddeo(f, 1.0, 0.0, π, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

# Various interval (-∞, 0)
let f(x) = sin(x)/x, expect = π/2
    I, E = quaddeo(f, 1.0, 0.0, -Inf, 0.0)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

# Various interval (-∞, b)
let f(x) = sin(x)/x, expect = π/2 - BigFloat("1.851937051982466170361053370157991363345809728981154909804")
    I, E = quaddeo(f, 1.0, 0.0, -Inf, -π)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

# Various interval (a, b)
let f(x) = sin(x)/x, expect = BigFloat("-0.43378547584983772011527320785824193420335637968613863796")
    I, E = quaddeo(f, 1.0, 0.0, π, 2π)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end
let f(x) = sin(x)/x, expect = BigFloat("0.43378547584983772011527320785824193420335637968613863796")
    I, E = quaddeo(f, 1.0, 0.0, -π, -2π)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end
let f(x) = sin(x)/x, expect = BigFloat("3.703874103964932340722106740315982726691619457962309819609")
    I, E = quaddeo(f, 1.0, 0.0, -π, π)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

# Corner case (a == b)
let f(x) = sin(x)/x, expect = π/2
    I, E = quaddeo(f, 1.0, 0.0, 0.0, 0.0)
    @test I isa Float64
    @test I ≈ 0
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

let f(x) = sin(x)/x, expect = π/2
    I, E = quaddeo(f, 1.0, 0.0, 1.0, 1.0)
    @test I isa Float64
    @test I ≈ 0
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

# Corner case (a > b)
let f(x) = sin(x)/x, expect = π/2
    I, E = quaddeo(f, 1.0, 0.0, Inf, 0.0)
    @test I isa Float64
    @test I ≈ -π/2
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

# Multiple integral intervals
let f(x) = sin(x)/x, expect = π/2
    I, E = quaddeo(f, 1.0, 0.0, 0.0, 1.0, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

let f(x) = sin(x)/x, expect = π/2
    I, E = quaddeo(f, 1.0, 0.0, 0.0, 1.0, 2.0, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I))*norm(I))
end

# Vectorized numerical integration
let f(x) = [sin(x)/x, exp(-x)*cos(x)], expect = [π/2, 1/2]
    I, E = quaddeo(f, 1.0, 0.0, 0.0, Inf)
    @test eltype(I) == Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I))*norm(I))
end

# Type handling
let f(x) = exp(-x)*cos(x), expect = BigFloat("0.5")
    ω = BigFloat(1.0)
    θ = BigFloat(π)/2
    a = BigFloat(0.0)
    b = BigFloat(Inf)

    for T in [Float32, Float64, BigFloat]
        ω_T, θ_T, a_T, b_T = T.((ω, θ, a, b))
        I, E = quaddeo(f, ω_T, θ_T, a_T, b_T)
        @test I isa T
        @test I ≈ expect
        @test E ≤ sqrt(eps(typeof(I))*norm(I))
    end
end


text:
Pkg.add("DataFrames")
Pkg.add("Gadfly")
Pkg.add("Cairo")
Pkg.add("Fontconfig")
using DataFrames
using Gadfly
using Cairo
using Fontconfig

# Plot all validities

style = Theme(
    major_label_font_size=35pt,
    minor_label_font_size=30pt,
    background_color="white",
    major_label_font="CMU Serif",
    minor_label_font="CMU Serif",
    bar_spacing=150px
)

numbers = ["377", "11177", "127413", "140054", "7506"];

for number in numbers
    println("Creating " , number , "...")
    table = readtable(string("output/" , number , "_answers.csv"));
    output = plot(
        table,
        x=:answer,
        y=:count,
        Geom.bar(position=:dodge),
        Coord.Cartesian(),
        Guide.xlabel("Answer"),
        Guide.ylabel("# of answers"),
        style
    );
    draw(PDF(string("graphs/" , number , "-answers.pdf"), 1600px, 800px), output);
end


text:
module TokyoVariantAnalysis4A
using Logging
using CSV
using DataFrames
using Dates
using LinearAlgebra: dot
using SpecialFunctions: beta_inc
using LsqFit: curve_fit, coef

using ..TokyoVariantAnalysis4

export likelihoods_df, likelihood_max

@info "module", @__MODULE__

@info "========"
@info "特定変異オッズ 最尤値探索"

# モデル
lh_model(t; α, β) = exp(α + β * t)

# 累積分布関数
# 不完全ベータ関数により二項分布の累積分布を連続値関数として表す
function Φ(z; n, zm)
    @assert z ≥ 0.0
    @assert n > 0
    @assert zm ≥ 0.0
    q = zm / (zm + 1.0)
    y = n * z / (z + 1.0)
    return beta_inc(n - y + 0.5, y + 0.5, 1.0 - q)[1]
end
# 確率密度関数（数値微分による）
const h_float64 = 2^-12
function φ(z; n, zm)
    @assert z ≥ 0.0
    @assert n > 0
    @assert zm ≥ 0.0
    h::Float64 = max(z, 1.0) * h_float64
    if (h ≤ z)
        return (Φ(z + h; n=n, zm=zm) - Φ(z - h; n=n, zm=zm)) / (2.0 * h)
    else
        return (Φ(z + h; n=n, zm=zm) - Φ(z; n=n, zm=zm)) / h
    end
end

const lh_ps_niter       =    5
const lh_ps_hwidth_α    =  5.0
const lh_ps_hwidth_β    =  1.0
const lh_ps_nslice_α    =  200
const lh_ps_nslice_β    =  200
const lh_nslice_α       =  500
const lh_nslice_β       = 1000
const lh_logL_threshold = -6.0
@info "予備探索回数", lh_ps_niter
@info "予備探索パラメータ初期探索幅", lh_ps_hwidth_α, lh_ps_hwidth_β
@info "予備探索パラメータ・スライス数", lh_ps_nslice_α, lh_ps_nslice_β
@info "本探索パラメータ・スライス数", lh_nslice_α, lh_nslice_β
@info "打ち切り尤度閾値", lh_logL_threshold

function search_odd_likelihoods(df, varname, period)::DataFrame
    oddname = varname*"_odds"
    (sdate_start, sdate_end) = period
    @info oddname
    @info sdate_start, sdate_end
    @info "--------"
    @info "データ点収集"
    @info "対象行のみを取り出した DataFrame を作成する"
    sample_df = filter(
        (row) -> (begin
            var_n = row[varname]
            odds = row[oddname]
            (sdate_start ≤ row.date_end) && (row.date_start ≤ sdate_end) &&
            !ismissing(var_n) && !ismissing(odds) && !isinf(odds) && !isnan(odds)
        end),
        df
    )
    @info size(sample_df)
    @info "--------"
    @info "サンプルサイズ"
    sample_size = nrow(sample_df)
    @info sample_size
    @info "--------"
    @info "基準日算出"
    date_epoch = round(dot(sample_df.date_value, sample_df.n) / sum(sample_df.n))
    @info date_epoch
    (t_start, t_end) = date_to_value.(period) .- date_epoch
    @info t_start, t_end
    sample_df.t = sample_df.date_value .- date_epoch
    @info "--------"
    @info "非線形フィッティングの予備的実行"
    @info "初期パラメーターを定めるため非線形最小二乗フィッティング"
    lsq_α = nothing
    lsq_β = nothing
    for once in 1:1
        cils = sample_df[!,oddname*"_cil"]
        cius = sample_df[!,oddname*"_ciu"]
        ws = 1.0 ./ (cius .- cils).^2
        βs0 = [0.0, 0.0]
        cf = curve_fit(
            (t, βs) -> lh_model.(t; α=βs[1], β=βs[2]),
            sample_df.t, sample_df[!,oddname], ws, βs0
        )
        (lsq_α, lsq_β) = coef(cf)
    end
    @info lsq_α, lsq_β
    @info "--------"
    @info "尤度予備探索"
    @info "計算すべきパラメーターの範囲を絞り込む"
    mn_α_range = nothing
    mn_β_range = nothing
    for once in 1:1
        α_range = range(lsq_α - lh_ps_hwidth_α, lsq_α + lh_ps_hwidth_α; length=lh_ps_nslice_α)
        β_range = range(lsq_β - lh_ps_hwidth_β, lsq_β + lh_ps_hwidth_β; length=lh_ps_nslice_β)
        for iter in 1:lh_ps_niter
            @info "予備探索 $(iter) 回目"
            @info α_range
            @info β_range
            ps_df = DataFrame(α=Float64[], β=Float64[], logL=Float64[])
            logL_max::Union{Float64,Nothing} = nothing
            for α in α_range
                for β in β_range
                    logL = 0.0
                    for row in eachrow(sample_df)
                        zm = lh_model(row.t; α=α, β=β)
                        logL += log(φ(row[oddname]; n=row.n, zm=zm))
                    end
                    if isnothing(logL_max) || (logL_max < logL)
                        logL_max = logL
                    end
                    if logL_max + lh_logL_threshold ≤ logL
                        push!(ps_df, [α, β, logL])
                    end
                end
            end
            @info logL_max
            α_min = last(α_range)
            α_max = first(α_range)
            β_min = last(β_range)
            β_max = first(β_range)
            for row in eachrow(ps_df)
                if logL_max + lh_logL_threshold ≤ row.logL
                    if row.α < α_min; α_min = row.α end
                    if α_max < row.α; α_max = row.α end
                    if row.β < β_min; β_min = row.β end
                    if β_max < row.β; β_max = row.β end
                end
            end
            α_torr = 1.5 * step(α_range)
            β_torr = 1.5 * step(β_range)
            α_range = range(α_min - α_torr, α_max + α_torr; length=lh_ps_nslice_α)
            β_range = range(β_min - β_torr, β_max + β_torr; length=lh_ps_nslice_β)
        end
        mn_α_range = range(first(α_range), last(α_range); length=lh_nslice_α)
        mn_β_range = range(first(β_range), last(β_range); length=lh_nslice_β)
    end
    @info "--------"
    @info "尤度本探索"
    @info "パラメーターの範囲を分割し、各パラメーターにおける尤度を計算"
    @info mn_α_range
    @info mn_β_range
    lh_df::Union{DataFrame,Nothing} = nothing
    for once in 1:1
        mn_df = DataFrame(α=Float64[], β=Float64[], logL=Float64[])
        logL_max::Union{Float64,Nothing} = nothing
        for (iα, α) in enumerate(mn_α_range)
            if iα % 10 == 0
                print(stderr, "\r$(round(Int, 100.0 * (iα / lh_nslice_α)))%")
            end
            for β in mn_β_range
                logL = 0.0
                for row in eachrow(sample_df)
                    zm = lh_model(row.t; α=α, β=β)
                    logL += log(φ(row[oddname]; n=row.n, zm=zm))
                end
                if isnothing(logL_max) || (logL_max < logL)
                    logL_max = logL
                end
                if logL_max + lh_logL_threshold ≤ logL
                    push!(mn_df, [α, β, logL])
                end
            end
        end
        print(stderr, "\r100%\n")
        lh_df = filter((row) -> logL_max + lh_logL_threshold ≤ row.logL, mn_df)
        for i in 1:nrow(lh_df)
            lh_df.α[i] = lh_df.α[i] - lh_df.β[i] * date_epoch
        end
        sort!(lh_df, :logL; rev=true)
        @info size(lh_df)
    end
    return lh_df
end

lh_varname = "L452R"
lh_sdate = Date("2021-06-01")
lh_edate = Date("2021-08-31")
lh_sdate_value = date_to_value(lh_sdate)
lh_edate_value = date_to_value(lh_edate)

likelihoods_df = search_odd_likelihoods(variants_df, lh_varname, (lh_sdate, lh_edate))
@info "size(likelihoods_df)", size(likelihoods_df)
likelihood_max = likelihoods_df[1,:]
@info "likelihood_max", likelihood_max

@info "--------"
@info "尤度データ書き出し"
const lh_output_csv_filename = "tokyo_variant_analysis4_likelihoods.csv"
@info "lh_output_csv_filename", lh_output_csv_filename
CSV.write(lh_output_csv_filename, likelihoods_df; missingstring="NA")

end #module

max_stars_repo_name: muimuiz/TokyoVariant

max_issues_repo_name: muimuiz/TokyoVariant

max_forks_repo_name: muimuiz/TokyoVariant