text:
"""
    TemporalHistory{M}

Linear event histories with marks of type `M`.

# Fields

- `times::Vector{Float64}`: vector of event times
- `marks::Vector{M}`: vector of event marks
- `tmin::Float64`: start time
- `tmax::Float64`: end time
"""
mutable struct TemporalHistory{M} <: AbstractHistory{Float64,M}
    times::Vector{Float64}
    marks::Vector{M}
    tmin::Float64
    tmax::Float64
end

event_times(h::TemporalHistory) = h.times

event_marks(h::TemporalHistory) = h.marks

"""
    min_time(h)

Return the starting time of `h` (not the same as the first event time).
"""
min_time(h::TemporalHistory) = h.tmin

"""
    max_time(h)

Return the end time of `h` (not the same as the last event time).
"""
max_time(h::TemporalHistory) = h.tmax

"""
    nb_events(h::TemporalHistory, tmin=-Inf, tmax=Inf)

Count events in `h` during the interval `[tmin, tmax)`.
"""
function nb_events(h::TemporalHistory, tmin = -Inf, tmax = Inf)
    i_min = searchsortedfirst(event_times(h), tmin)
    i_max = searchsortedlast(event_times(h), tmax - eps(tmax))
    return i_max - i_min + 1
end

"""
    has_events(h::TemporalHistory, tmin=-Inf, tmax=Inf)

Check the presence of events in `h` during the interval `[tmin, tmax)`.
"""
function has_events(h::TemporalHistory, tmin = -Inf, tmax = Inf)
    return nb_events(h, tmin, tmax) > 0
end

"""
    duration(h::TemporalHistory)

Compute the difference `h.tmax - h.tmin`.
"""
function duration(h::TemporalHistory)
    return max_time(h) - min_time(h)
end

"""
    push!(h::TemporalHistory{M}, t::Float64, m::M)

Add event `(t, m)` at the end of history `h`.
"""
function Base.push!(h::TemporalHistory, t, m)
    @assert h.tmin <= t < h.tmax
    push!(h.times, t)
    push!(h.marks, m)
    return nothing
end

"""
    append!(h1::TemporalHistory, h2::TemporalHistory)

Add all the events of `h2` at the end of `h1`.
"""
function Base.append!(h1::TemporalHistory, h2::TemporalHistory)
    if has_events(h1) && has_events(h2)
        @assert maximum(event_times(h1)) < minimum(event_times(h2))
    end
    for (t, m) in zip(event_times(h2), event_marks(h2))
        push!(h1, t, m)
    end
    return nothing
end

@doc raw"""
    time_change(h, Λ)

Apply the time rescaling $t \mapsto \Lambda(t)$ to history `h`.
"""
function time_change(h::TemporalHistory, Λ)
    new_times = Λ.(event_times(h))
    new_marks = copy(event_marks(h))
    new_tmin = Λ(min_time(h))
    new_tmax = Λ(max_time(h))
    return TemporalHistory(new_times, new_marks, new_tmin, new_tmax)
end


text:
using easyGit
gitpush("easyGit")

ext: jl


text:
#### Boundary conditions

export TurbConvBC,
    NoTurbConvBC,
    turbconv_bcs,
    turbconv_boundary_state!,
    turbconv_normal_boundary_flux_second_order!

abstract type TurbConvBC end

"""
    NoTurbConvBC <: TurbConvBC

Boundary conditions are not applied
"""
struct NoTurbConvBC <: TurbConvBC end

turbconv_bcs(::NoTurbConv) = (NoTurbConvBC(), NoTurbConvBC())

function turbconv_boundary_state!(nf, bc_turbulence::NoTurbConvBC, bl, _...)
    nothing
end

function turbconv_normal_boundary_flux_second_order!(
    nf,
    bc_turbulence::NoTurbConvBC,
    bl,
    _...,
)
    nothing
end


text:
using TestPkg
using Test

@testset "TestPkg.jl" begin
  @test f1(5) == 10
end


text:
# convert the vpv,vph,vsv,vsh to bulk_c_v and bulk_c_h, which is used for the PSF test.
include("../specfem_gll.jl/src/utils/readfiles.jl")

max_stars_repo_name: ziyixi/seisflow

max_issues_repo_name: ziyixi/seisflow

max_forks_repo_name: ziyixi/seisflow


text:
function LevenbergMarquardtOptimizer()
    return LevenbergMarquardtOptimizer(())
end

function LevenbergMarquardtOptimizer(arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble)
    return LevenbergMarquardtOptimizer((jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function get_cost_relative_tolerance(obj::LevenbergMarquardtOptimizer)
    return jcall(obj, "getCostRelativeTolerance", jdouble, ())
end

function get_initial_step_bound_factor(obj::LevenbergMarquardtOptimizer)
    return jcall(obj, "getInitialStepBoundFactor", jdouble, ())
end

function get_ortho_tolerance(obj::LevenbergMarquardtOptimizer)
    return jcall(obj, "getOrthoTolerance", jdouble, ())
end

function get_parameter_relative_tolerance(obj::LevenbergMarquardtOptimizer)
    return jcall(obj, "getParameterRelativeTolerance", jdouble, ())
end

function get_ranking_threshold(obj::LevenbergMarquardtOptimizer)
    return jcall(obj, "getRankingThreshold", jdouble, ())
end

function optimize(obj::LevenbergMarquardtOptimizer, arg0::LeastSquaresProblem)
    return jcall(obj, "optimize", LeastSquaresOptimizer_Optimum, (LeastSquaresProblem,), arg0)
end

function with_cost_relative_tolerance(obj::LevenbergMarquardtOptimizer, arg0::jdouble)
    return jcall(obj, "withCostRelativeTolerance", LevenbergMarquardtOptimizer, (jdouble,), arg0)
end

function with_initial_step_bound_factor(obj::LevenbergMarquardtOptimizer, arg0::jdouble)
    return jcall(obj, "withInitialStepBoundFactor", LevenbergMarquardtOptimizer, (jdouble,), arg0)
end

function with_ortho_tolerance(obj::LevenbergMarquardtOptimizer, arg0::jdouble)
    return jcall(obj, "withOrthoTolerance", LevenbergMarquardtOptimizer, (jdouble,), arg0)
end

function with_parameter_relative_tolerance(obj::LevenbergMarquardtOptimizer, arg0::jdouble)
    return jcall(obj, "withParameterRelativeTolerance", LevenbergMarquardtOptimizer, (jdouble,), arg0)
end

function with_ranking_threshold(obj::LevenbergMarquardtOptimizer, arg0::jdouble)
    return jcall(obj, "withRankingThreshold", LevenbergMarquardtOptimizer, (jdouble,), arg0)
end


text:
const SimpleDiGraphEdge = SimpleEdge

"""
    SimpleDiGraph{T}

A type representing a directed graph.
"""
mutable struct SimpleDiGraph{T <: Integer} <: AbstractSimpleGraph{T}
    ne::Int
    fadjlist::Vector{Vector{T}} # [src]: (dst, dst, dst)
    badjlist::Vector{Vector{T}} # [dst]: (src, src, src)
end

eltype(x::SimpleDiGraph{T}) where T = T

# DiGraph{UInt8}(6), DiGraph{Int16}(7), DiGraph{Int8}()
"""
    SimpleDiGraph{T}(n=0)

Construct a `SimpleDiGraph{T}` with `n` vertices and 0 edges.
If not specified, the element type `T` is the type of `n`.

## Examples
```jldoctest
julia> SimpleDiGraph(UInt8(10))
{10, 0} directed simple UInt8 graph
```
"""
function SimpleDiGraph{T}(n::Integer=0) where T <: Integer
    fadjlist = [Vector{T}() for _ = one(T):n]
    badjlist = [Vector{T}() for _ = one(T):n]
    return SimpleDiGraph(0, fadjlist, badjlist)
end

# SimpleDiGraph(6), SimpleDiGraph(0x5)
SimpleDiGraph(n::T) where T <: Integer = SimpleDiGraph{T}(n)

# SimpleDiGraph()
SimpleDiGraph() = SimpleDiGraph{Int}()

# SimpleDiGraph(UInt8)
"""
    SimpleDiGraph(::Type{T})

Construct an empty `SimpleDiGraph{T}` with 0 vertices and 0 edges.

## Examples
```jldoctest
julia> SimpleDiGraph(UInt8)
{0, 0} directed simple UInt8 graph
```
"""
SimpleDiGraph(::Type{T}) where T <: Integer = SimpleDiGraph{T}(zero(T))

# SimpleDiGraph(adjmx)
"""
    SimpleDiGraph{T}(adjm::AbstractMatrix)

Construct a `SimpleDiGraph{T}` from the adjacency matrix `adjm`.
If `adjm[i][j] != 0`, an edge `(i, j)` is inserted. `adjm` must be a square matrix.
The element type `T` can be omitted.

## Examples
```jldoctest
julia> A1 = [false true; false false]
julia> SimpleDiGraph(A1)
{2, 1} directed simple Int64 graph

julia> A2 = [2 7; 5 0]
julia> SimpleDiGraph{Int16}(A2)
{2, 3} directed simple Int16 graph
```
"""
SimpleDiGraph(adjmx::AbstractMatrix) = SimpleDiGraph{Int}(adjmx)

# sparse adjacency matrix constructor: SimpleDiGraph(adjmx)
function SimpleDiGraph{T}(adjmx::SparseMatrixCSC{U}) where T <: Integer where U <: Real
    dima, dimb = size(adjmx)
    isequal(dima, dimb) || throw(ArgumentError("Adjacency / distance matrices must be square"))

    g = SimpleDiGraph(T(dima))
    maxc = length(adjmx.colptr)
    @inbounds for c = 1:(maxc - 1)
        for rind = adjmx.colptr[c]:(adjmx.colptr[c + 1] - 1)
            isnz = (adjmx.nzval[rind] != zero(U))
            if isnz
                r = adjmx.rowval[rind]
                add_edge!(g, r, c)
            end
        end
    end
    return g
end

# dense adjacency matrix constructor: DiGraph{UInt8}(adjmx)
function SimpleDiGraph{T}(adjmx::AbstractMatrix{U}) where T <: Integer where U <: Real
    dima, dimb = size(adjmx)
    isequal(dima, dimb) || throw(ArgumentError("Adjacency / distance matrices must be square"))

    g = SimpleDiGraph(T(dima))
    @inbounds for i in findall(adjmx .!= zero(U))
        add_edge!(g, i[1], i[2])
    end
    return g
end

# converts DiGraph{Int} to DiGraph{Int32}
"""
    SimpleDiGraph{T}(g::SimpleDiGraph)

Construct a copy of g.
If the element type `T` is specified, the vertices of `g` are converted to this type.
Otherwise the element type is the same as for `g`.

## Examples
```jldoctest
julia> g = CompleteDiGraph(5)
julia> SimpleDiGraph{UInt8}(g)
{5, 20} directed simple UInt8 graph
```
"""
function SimpleDiGraph{T}(g::SimpleDiGraph) where T <: Integer
    h_fadj = [Vector{T}(x) for x in fadj(g)]
    h_badj = [Vector{T}(x) for x in badj(g)]
    return SimpleDiGraph(ne(g), h_fadj, h_badj)
end

SimpleDiGraph(g::SimpleDiGraph) = copy(g)

# constructor from abstract graph: SimpleDiGraph(graph)
"""
    SimpleDiGraph(g::AbstractSimpleGraph)

Construct an directed `SimpleDiGraph` from a graph `g`.
The element type is the same as for `g`.

## Examples
```jldoctest
julia> g = PathGraph(Int8(5))
julia> SimpleDiGraph(g)
{5, 8} directed simple Int8 graph
```
"""
function SimpleDiGraph(g::AbstractSimpleGraph)
    h = SimpleDiGraph(nv(g))
    h.ne = ne(g) * 2 - num_self_loops(g)
    h.fadjlist = deepcopy(fadj(g))
    h.badjlist = deepcopy(badj(g))
    return h
end

@inbounds function cleanupedges!(fadjlist::Vector{Vector{T}},
                                 badjlist::Vector{Vector{T}}) where T <: Integer
    neg = 0
    for v in 1:length(fadjlist)
        if !issorted(fadjlist[v])
            sort!(fadjlist[v])
        end
        if !issorted(badjlist[v])
            sort!(badjlist[v])
        end
        unique!(fadjlist[v])
        unique!(badjlist[v])
        neg += length(fadjlist[v])
    end
    return neg
end

"""
    SimpleDiGraph(edge_list::Vector)

Construct a `SimpleDiGraph` from a vector of edges.
The element type is taken from the edges in `edge_list`.
The number of vertices is the highest that is used in an edge in `edge_list`.

### Implementation Notes
This constructor works the fastest when `edge_list` is sorted
by the lexical ordering and does not contain any duplicates.

### See also
[`SimpleDiGraphFromIterator`](@ref)

## Examples
```jldoctest

julia> el = Edge.([ (1, 3), (1, 5), (3, 1) ])
julia> SimpleDiGraph(el)
{5, 3} directed simple Int64 graph
```
"""
function SimpleDiGraph(edge_list::Vector{SimpleDiGraphEdge{T}}) where T <: Integer
    nvg = zero(T)
    @inbounds(
    for e in edge_list
        nvg = max(nvg, src(e), dst(e))
    end)

    list_sizes_out = ones(Int, nvg)
    list_sizes_in = ones(Int, nvg)
    degs_out = zeros(Int, nvg)
    degs_in = zeros(Int, nvg)
    @inbounds(
    for e in edge_list
        s, d = src(e), dst(e)
        (s >= 1 && d >= 1) || continue
        degs_out[s] += 1
        degs_in[d] += 1
    end)

    fadjlist = Vector{Vector{T}}(undef, nvg)
    badjlist = Vector{Vector{T}}(undef, nvg)
    @inbounds(
    for v in 1:nvg
        fadjlist[v] = Vector{T}(undef, degs_out[v])
        badjlist[v] = Vector{T}(undef, degs_in[v])
    end)

    @inbounds(
    for e in edge_list
        s, d = src(e), dst(e)
        (s >= 1 && d >= 1) || continue
        fadjlist[s][list_sizes_out[s]] = d
        list_sizes_out[s] += 1
        badjlist[d][list_sizes_in[d]] = s
        list_sizes_in[d] += 1
    end)

    neg = cleanupedges!(fadjlist, badjlist)
    g = SimpleDiGraph{T}()
    g.fadjlist = fadjlist
    g.badjlist = badjlist
    g.ne = neg

    return g
end

@inbounds function add_to_lists!(fadjlist::Vector{Vector{T}},
                                 badjlist::Vector{Vector{T}}, s::T, d::T) where T <: Integer
    nvg = length(fadjlist)
    nvg_new = max(nvg, s, d)
    for v = (nvg + 1):nvg_new
        push!(fadjlist, Vector{T}())
        push!(badjlist, Vector{T}())
    end

    push!(fadjlist[s], d)
    push!(badjlist[d], s)
end

function _SimpleDiGraphFromIterator(iter)::SimpleDiGraph
    T = Union{}
    fadjlist = Vector{Vector{T}}()
    badjlist = Vector{Vector{T}}()
    @inbounds(
    for e in iter
        typeof(e) <: SimpleDiGraphEdge ||
                        throw(ArgumentError("iter must be an iterator over SimpleEdge"))
        s, d = src(e), dst(e)
        (s >= 1 && d >= 1) || continue
        if T != eltype(e)
            T = typejoin(T, eltype(e))
            fadjlist = convert(Vector{Vector{T}}, fadjlist)
            badjlist = convert(Vector{Vector{T}}, badjlist)
        end
        add_to_lists!(fadjlist, badjlist, s, d)
    end)

    T == Union{} && return SimpleDiGraph(0)
    neg  = cleanupedges!(fadjlist, badjlist)
    g = SimpleDiGraph{T}()
    g.fadjlist = fadjlist
    g.badjlist = badjlist
    g.ne = neg

    return g
end

function _SimpleDiGraphFromIterator(iter, ::Type{SimpleDiGraphEdge{T}}) where T <: Integer
    fadjlist = Vector{Vector{T}}()
    badjlist = Vector{Vector{T}}()
    @inbounds(
    for e in iter
        s, d = src(e), dst(e)
        (s >= 1 && d >= 1) || continue
        add_to_lists!(fadjlist, badjlist, s, d)
    end)

    g = SimpleDiGraph{T}()
    neg  = cleanupedges!(fadjlist, badjlist)
    g.fadjlist = fadjlist
    g.badjlist = badjlist
    g.ne = neg

    return g
end

"""
    SimpleDiGraphFromIterator(iter)

Create a `SimpleDiGraph` from an iterator `iter`. The elements in `iter` must
be of `type <: SimpleEdge`.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = SimpleDiGraph(2);

julia> add_edge!(g, 1, 2);

julia> add_edge!(g, 2, 1);

julia> h = SimpleDiGraphFromIterator(edges(g))
{2, 2} directed simple Int64 graph

julia> collect(edges(h))
2-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 1 => 2
 Edge 2 => 1
```
"""
function SimpleDiGraphFromIterator(iter)::SimpleDiGraph
    if Base.IteratorEltype(iter) == Base.EltypeUnknown()
        return _SimpleDiGraphFromIterator(iter)
    end
    if eltype(iter) <: SimpleGraphEdge && isconcretetype(eltype(iter))
        return _SimpleDiGraphFromIterator(iter, eltype(iter))
    end
    return _SimpleDiGraphFromIterator(iter)
end

edgetype(::SimpleDiGraph{T}) where T <: Integer = SimpleGraphEdge{T}

badj(g::SimpleDiGraph) = g.badjlist
badj(g::SimpleDiGraph, v::Integer) = badj(g)[v]

copy(g::SimpleDiGraph{T}) where T <: Integer =
SimpleDiGraph{T}(g.ne, deepcopy(g.fadjlist), deepcopy(g.badjlist))

==(g::SimpleDiGraph, h::SimpleDiGraph) =
vertices(g) == vertices(h) &&
ne(g) == ne(h) &&
fadj(g) == fadj(h) &&
badj(g) == badj(h)

is_directed(g::SimpleDiGraph) = true
is_directed(::Type{SimpleDiGraph}) = true
is_directed(::Type{SimpleDiGraph{T}}) where T = true

function has_edge(g::SimpleDiGraph{T}, s, d) where T
    verts = vertices(g)
    (s in verts && d in verts) || return false  # edge out of bounds
    @inbounds list = g.fadjlist[s]
    @inbounds list_backedge = g.badjlist[d]
    if length(list) > length(list_backedge)
        d = s
        list = list_backedge
    end
    return insorted(d, list)
end

function has_edge(g::SimpleDiGraph{T}, e::SimpleDiGraphEdge{T}) where T
    s, d = T.(Tuple(e))
    return has_edge(g, s, d)
end

function add_edge!(g::SimpleDiGraph{T}, e::SimpleDiGraphEdge{T}) where T
    s, d = T.(Tuple(e))
    verts = vertices(g)
    (s in verts && d in verts) || return false  # edge out of bounds
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) && return false  # edge already in graph
    insert!(list, index, d)

    g.ne += 1

    @inbounds list = g.badjlist[d]
    index = searchsortedfirst(list, s)
    insert!(list, index, s)
    return true  # edge successfully added
end

function rem_edge!(g::SimpleDiGraph{T}, e::SimpleDiGraphEdge{T}) where T
    s, d = T.(Tuple(e))
    verts = vertices(g)
    (s in verts && d in verts) || return false  # edge out of bounds
    @inbounds list = g.fadjlist[s]
    index = searchsortedfirst(list, d)
    @inbounds (index <= length(list) && list[index] == d) || return false   # edge not in graph
    deleteat!(list, index)

    g.ne -= 1

    @inbounds list = g.badjlist[d]
    index = searchsortedfirst(list, s)
    deleteat!(list, index)
    return true # edge successfully removed
end

function add_vertex!(g::SimpleDiGraph{T}) where T
    (nv(g) + one(T) <= nv(g)) && return false       # test for overflow
    push!(g.badjlist, Vector{T}())
    push!(g.fadjlist, Vector{T}())

    return true
end


text:
typealias SamIOStream BamIOStream

@doc """ sam file name must endwith bam
""" ->
function sam_open(fname::AbstractString, mode::AbstractString)
    fname = realpath(fname)
    split(fname,".")[end] == "sam" ||
    split(fname,".")[end] == "bam" || info("file name doesn't endwith sam/bam")

    sam_fl = hts_open(fname,mode)
    if mode == "r"
        phdr = sam_hdr_read(sam_fl) #header info
    else
        phdr = bam_hdr_init()
    end
    prec = bam_init1()
    kstr = KString(0,0,C_NULL)
    pkstr = convert(Ptr{KString}, pointer_from_objref(kstr))
    SamIOStream(fname,sam_fl,phdr,prec,pkstr)
end

function sam_hdr_parse(str::AbstractString)
    l_text = Int32(length(str))
    str = string(str,'\0')
    text = convert(Cstring, pointer(str))
    sam_hdr_parse(l_text,text)
end

@doc """ read a record from a sam iostream
""" ->
function readline(sios::SamIOStream)
    handle = sios.handle
    phdr = sios.phdr
    prec = sios.precord
    ret = sam_read!(handle,phdr,prec)
    if ret == 0
        sios.eof = false
    elseif ret == -1
        sios.eof = true
    end
    pkstr = sios.pkstr
    sam_format!(phdr,prec,pkstr)

    strptr(pkstr)
end

@doc """ eof(sios) -> Bool
         Tests whether an Bam I/O stream is at end-of-file.
""" ->
function eof(sios::SamIOStream)
    sios.eof
end

@doc """ read all the records from a bam file
""" ->
function readlines(fname::AbstractString)
    sios = sam_open(fname,"r")
    data = readlines(sios)
    close(sios)
    data
end

@doc """ read all the records from a SamIOStream
""" ->
function readlines(sios::SamIOStream)
    size = 256
    data = Array{ASCIIString,1}(size)
    i = 0
    while !eof(sios)
        i += 1
        if i > size
            resize!(data,2*size)
            size = 2*size
        end
        line = readline(sios)
        data[i] = line
    end
    data[1:i]
end

@doc """ write a record to a SamIOStream
""" ->
function writeline(sios::SamIOStream, str::ASCIIString)
#    error("Wait for hdr function finished")
    kstr = KString(str)
    pkstr = convert(Ptr{KString}, pointer_from_objref(kstr))
    handle = sios.handle
    phdr = sios.phdr
    prec = sios.precord
    sam_parse1(pkstr, phdr, prec)
    sam_write1(handle,phdr,prec)
end

@doc """ write all the records in strs to a SamIOStream
""" ->
function writelines(sios::SamIOStream, strs::Array{ASCIIString,1})
    sam_hdr_write(sios.handle,sios.phdr)
    for str in strs
        writeline(sios, str)
    end
end

@doc """ close a SamIOStream
""" ->
function close(sios::SamIOStream)
    bam_destroy1(sios.precord)
#    bam_hdr_destroy(sios.phdr) # not yet to find why freed already
    sam_close(sios.handle)
end


text:
module Formatting

    import Base.show

    export
        FormatSpec, FormatExpr,
        printfmt, printfmtln, fmt, format,
        sprintf1, generate_formatter

    using Compat

    include("cformat.jl" )
    include("fmtspec.jl")
    include("fmtcore.jl")
    include("formatexpr.jl")

end # module


text:
int_rules_1_2_2_8 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.2.2.8*P(x)*(d+e*x)^q*(a+b*x^2+c*x^4)^p =#
    @apply_utils Antiderivative(1 / ((~d + ~(e') * ~x) * sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4)), ~x) => ~d * Antiderivative(1 / (((~d) ^ 2 - (~e) ^ 2 * (~x) ^ 2) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)), ~x) - ~e * Antiderivative(~x / (((~d) ^ 2 - (~e) ^ 2 * (~x) ^ 2) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x)
    @apply_utils Antiderivative(1 / ((~d + ~(e') * ~x) * sqrt(~a + ~(c') * (~x) ^ 4)), ~x) => ~d * Antiderivative(1 / (((~d) ^ 2 - (~e) ^ 2 * (~x) ^ 2) * sqrt(~a + ~c * (~x) ^ 4)), ~x) - ~e * Antiderivative(~x / (((~d) ^ 2 - (~e) ^ 2 * (~x) ^ 2) * sqrt(~a + ~c * (~x) ^ 4)), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x)
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~q / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => ((~e) ^ 3 * (~d + ~e * ~x) ^ (~q + 1) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)) / ((~q + 1) * (~c * (~d) ^ 4 + ~b * (~d) ^ 2 * (~e) ^ 2 + ~a * (~e) ^ 4)) + (1 / ((~q + 1) * (~c * (~d) ^ 4 + ~b * (~d) ^ 2 * (~e) ^ 2 + ~a * (~e) ^ 4))) * Antiderivative(((~d + ~e * ~x) ^ (~q + 1) / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)) * Simp(((~d * (~q + 1) * (~c * (~d) ^ 2 + ~b * (~e) ^ 2) - ~e * (~c * (~d) ^ 2 * (~q + 1) + ~b * (~e) ^ 2 * (~q + 2)) * ~x) + ~c * ~d * (~e) ^ 2 * (~q + 1) * (~x) ^ 2) - ~c * (~e) ^ 3 * (~q + 3) * (~x) ^ 3, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ(~c * (~d) ^ 4 + ~b * (~d) ^ 2 * (~e) ^ 2 + ~a * (~e) ^ 4, 0) && ILtQ(~q, -1))
    @apply_utils Antiderivative((~d + ~(e') * ~x) ^ ~q / sqrt(~a + ~(c') * (~x) ^ 4), ~x) => ((~e) ^ 3 * (~d + ~e * ~x) ^ (~q + 1) * sqrt(~a + ~c * (~x) ^ 4)) / ((~q + 1) * (~c * (~d) ^ 4 + ~a * (~e) ^ 4)) + (~c / ((~q + 1) * (~c * (~d) ^ 4 + ~a * (~e) ^ 4))) * Antiderivative(((~d + ~e * ~x) ^ (~q + 1) / sqrt(~a + ~c * (~x) ^ 4)) * Simp((((~d) ^ 3 * (~q + 1) - (~d) ^ 2 * ~e * (~q + 1) * ~x) + ~d * (~e) ^ 2 * (~q + 1) * (~x) ^ 2) - (~e) ^ 3 * (~q + 3) * (~x) ^ 3, ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NeQ(~c * (~d) ^ 4 + ~a * (~e) ^ 4, 0) && ILtQ(~q, -1))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~(p') / (~d + ~(e') * ~x), ~x) => ~d * Antiderivative((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p / ((~d) ^ 2 - (~e) ^ 2 * (~x) ^ 2), ~x) - ~e * Antiderivative((~x * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p) / ((~d) ^ 2 - (~e) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && IntegerQ(~p + 1 / 2)
    @apply_utils Antiderivative((~a + ~(c') * (~x) ^ 4) ^ ~(p') / (~d + ~(e') * ~x), ~x) => ~d * Antiderivative((~a + ~c * (~x) ^ 4) ^ ~p / ((~d) ^ 2 - (~e) ^ 2 * (~x) ^ 2), ~x) - ~e * Antiderivative((~x * (~a + ~c * (~x) ^ 4) ^ ~p) / ((~d) ^ 2 - (~e) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && IntegerQ(~p + 1 / 2)
    @apply_utils Antiderivative(~Px * (~d + ~(e') * ~x) ^ ~(q') * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => Antiderivative(PolynomialQuotient(~Px, ~d + ~e * ~x, ~x) * (~d + ~e * ~x) ^ (~q + 1) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p, ~q], ~x) && (PolyQ(~Px, ~x) && EqQ(PolynomialRemainder(~Px, ~d + ~e * ~x, ~x), 0))
    @apply_utils Antiderivative(~Px * (~d + ~(e') * ~x) ^ ~(q') * (~a + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => Antiderivative(PolynomialQuotient(~Px, ~d + ~e * ~x, ~x) * (~d + ~e * ~x) ^ (~q + 1) * (~a + ~c * (~x) ^ 4) ^ ~p, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~p, ~q], ~x) && (PolyQ(~Px, ~x) && EqQ(PolynomialRemainder(~Px, ~d + ~e * ~x, ~x), 0))
    @apply_utils Antiderivative(~Px * (~d + ~(e') * ~x) ^ ~(q') * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => Antiderivative(PolynomialQuotient(~Px, ~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4, ~x) * (~d + ~e * ~x) ^ ~q * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p, ~q], ~x) && (PolyQ(~Px, ~x) && EqQ(PolynomialRemainder(~Px, ~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4, ~x), 0))
    @apply_utils Antiderivative(~Px * (~d + ~(e') * ~x) ^ ~(q') * (~a + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => Antiderivative(PolynomialQuotient(~Px, ~a + ~c * (~x) ^ 4, ~x) * (~d + ~e * ~x) ^ ~q * (~a + ~c * (~x) ^ 4) ^ (~p + 1), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~p, ~q], ~x) && (PolyQ(~Px, ~x) && EqQ(PolynomialRemainder(~Px, ~a + ~c * (~x) ^ 4, ~x), 0))
    @apply_utils Antiderivative((~Px * (~d + ~(e') * ~x) ^ ~q) / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([A = Coeff(~Px, ~x, 0), B = Coeff(~Px, ~x, 1), C = Coeff(~Px, ~x, 2), D = Coeff(~Px, ~x, 3)], Antiderivative(((~d + ~e * ~x) ^ (~q - 1) * (A * ~d + (B * ~d + A * ~e) * ~x + (C * ~d + B * ~e) * (~x) ^ 2 + C * ~e * (~x) ^ 3)) / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (PolyQ(~Px, ~x) && (LeQ(Expon(~Px, ~x), 2) && (NeQ(~c * (~d) ^ 4 + ~b * (~d) ^ 2 * (~e) ^ 2 + ~a * (~e) ^ 4, 0) && GtQ(~q, 0))))
    @apply_utils Antiderivative((~Px * (~d + ~(e') * ~x) ^ ~q) / sqrt(~a + ~(c') * (~x) ^ 4), ~x) => With([A = Coeff(~Px, ~x, 0), B = Coeff(~Px, ~x, 1), C = Coeff(~Px, ~x, 2), D = Coeff(~Px, ~x, 3)], Antiderivative(((~d + ~e * ~x) ^ (~q - 1) * (A * ~d + (B * ~d + A * ~e) * ~x + (C * ~d + B * ~e) * (~x) ^ 2 + C * ~e * (~x) ^ 3)) / sqrt(~a + ~c * (~x) ^ 4), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (PolyQ(~Px, ~x) && (LeQ(Expon(~Px, ~x), 2) && (NeQ(~c * (~d) ^ 4 + ~a * (~e) ^ 4, 0) && GtQ(~q, 0))))
    @apply_utils Antiderivative((~Px * (~d + ~(e') * ~x) ^ ~q) / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([A = Coeff(~Px, ~x, 0), B = Coeff(~Px, ~x, 1), C = Coeff(~Px, ~x, 2), D = Coeff(~Px, ~x, 3)], (D * (~d + ~e * ~x) ^ ~q * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)) / (~c * (~q + 2)) - (1 / (~c * (~q + 2))) * Antiderivative(((~d + ~e * ~x) ^ (~q - 1) / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)) * Simp(((~a * D * ~e * ~q - A * ~c * ~d * (~q + 2)) + ((~b * ~d * D - B * ~c * ~d * (~q + 2)) - A * ~c * ~e * (~q + 2)) * ~x + (~b * D * ~e * (~q + 1) - ~c * (C * ~d + B * ~e) * (~q + 2)) * (~x) ^ 2) - ~c * (~d * D * ~q + C * ~e * (~q + 2)) * (~x) ^ 3, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (PolyQ(~Px, ~x, 3) && (NeQ(~c * (~d) ^ 4 + ~b * (~d) ^ 2 * (~e) ^ 2 + ~a * (~e) ^ 4, 0) && GtQ(~q, 0)))
    @apply_utils Antiderivative((~Px * (~d + ~(e') * ~x) ^ ~q) / sqrt(~a + ~(c') * (~x) ^ 4), ~x) => With([A = Coeff(~Px, ~x, 0), B = Coeff(~Px, ~x, 1), C = Coeff(~Px, ~x, 2), D = Coeff(~Px, ~x, 3)], (D * (~d + ~e * ~x) ^ ~q * sqrt(~a + ~c * (~x) ^ 4)) / (~c * (~q + 2)) - (1 / (~c * (~q + 2))) * Antiderivative(((~d + ~e * ~x) ^ (~q - 1) / sqrt(~a + ~c * (~x) ^ 4)) * Simp((((~a * D * ~e * ~q - A * ~c * ~d * (~q + 2)) - ~c * (B * ~d * (~q + 2) + A * ~e * (~q + 2)) * ~x) - ~c * (C * ~d + B * ~e) * (~q + 2) * (~x) ^ 2) - ~c * (~d * D * ~q + C * ~e * (~q + 2)) * (~x) ^ 3, ~x), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (PolyQ(~Px, ~x, 3) && (NeQ(~c * (~d) ^ 4 + ~a * (~e) ^ 4, 0) && GtQ(~q, 0)))
    @apply_utils Antiderivative((~Px * (~d + ~(e') * ~x) ^ ~q) / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([A = Coeff(~Px, ~x, 0), B = Coeff(~Px, ~x, 1), C = Coeff(~Px, ~x, 2), D = Coeff(~Px, ~x, 3)], (-(((((~d) ^ 3 * D - C * (~d) ^ 2 * ~e) + B * ~d * (~e) ^ 2) - A * (~e) ^ 3)) * (~d + ~e * ~x) ^ (~q + 1) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)) / ((~q + 1) * (~c * (~d) ^ 4 + ~b * (~d) ^ 2 * (~e) ^ 2 + ~a * (~e) ^ 4)) + (1 / ((~q + 1) * (~c * (~d) ^ 4 + ~b * (~d) ^ 2 * (~e) ^ 2 + ~a * (~e) ^ 4))) * Antiderivative(((~d + ~e * ~x) ^ (~q + 1) / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)) * Simp(((~q + 1) * (~a * ~e * (((~d) ^ 2 * D - C * ~d * ~e) + B * (~e) ^ 2) + A * ~d * (~c * (~d) ^ 2 + ~b * (~e) ^ 2)) - ((~e * (~q + 1) * (A * ~c * (~d) ^ 2 + ~a * ~e * (~d * D - C * ~e)) - B * ~d * (~c * (~d) ^ 2 * (~q + 1) + ~b * (~e) ^ 2 * (~q + 2))) - ~b * (((~d) ^ 3 * D - C * (~d) ^ 2 * ~e) - A * (~e) ^ 3 * (~q + 2))) * ~x) + (~q + 1) * (D * ~e * (~b * (~d) ^ 2 + ~a * (~e) ^ 2) + ~c * ~d * (C * (~d) ^ 2 - ~e * (B * ~d - A * ~e))) * (~x) ^ 2 + ~c * (~q + 3) * ((((~d) ^ 3 * D - C * (~d) ^ 2 * ~e) + B * ~d * (~e) ^ 2) - A * (~e) ^ 3) * (~x) ^ 3, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (PolyQ(~Px, ~x) && (LeQ(Expon(~Px, ~x), 3) && (NeQ(~c * (~d) ^ 4 + ~b * (~d) ^ 2 * (~e) ^ 2 + ~a * (~e) ^ 4, 0) && LtQ(~q, -1))))
    @apply_utils Antiderivative((~Px * (~d + ~(e') * ~x) ^ ~q) / sqrt(~a + ~(c') * (~x) ^ 4), ~x) => With([A = Coeff(~Px, ~x, 0), B = Coeff(~Px, ~x, 1), C = Coeff(~Px, ~x, 2), D = Coeff(~Px, ~x, 3)], (-(((((~d) ^ 3 * D - C * (~d) ^ 2 * ~e) + B * ~d * (~e) ^ 2) - A * (~e) ^ 3)) * (~d + ~e * ~x) ^ (~q + 1) * sqrt(~a + ~c * (~x) ^ 4)) / ((~q + 1) * (~c * (~d) ^ 4 + ~a * (~e) ^ 4)) + (1 / ((~q + 1) * (~c * (~d) ^ 4 + ~a * (~e) ^ 4))) * Antiderivative(((~d + ~e * ~x) ^ (~q + 1) / sqrt(~a + ~c * (~x) ^ 4)) * Simp(((~q + 1) * (~a * ~e * (((~d) ^ 2 * D - C * ~d * ~e) + B * (~e) ^ 2) + A * ~d * (~c * (~d) ^ 2)) - (~e * (~q + 1) * (A * ~c * (~d) ^ 2 + ~a * ~e * (~d * D - C * ~e)) - B * ~d * (~c * (~d) ^ 2 * (~q + 1))) * ~x) + (~q + 1) * (D * ~e * (~a * (~e) ^ 2) + ~c * ~d * (C * (~d) ^ 2 - ~e * (B * ~d - A * ~e))) * (~x) ^ 2 + ~c * (~q + 3) * ((((~d) ^ 3 * D - C * (~d) ^ 2 * ~e) + B * ~d * (~e) ^ 2) - A * (~e) ^ 3) * (~x) ^ 3, ~x), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (PolyQ(~Px, ~x) && (LeQ(Expon(~Px, ~x), 3) && (NeQ(~c * (~d) ^ 4 + ~a * (~e) ^ 4, 0) && LtQ(~q, -1))))
    @apply_utils Antiderivative((~A + ~(B') * ~x) / ((~d + ~(e') * ~x) * sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4)), ~x) => ((-((~A) ^ 2) * (~B * ~d + ~A * ~e)) / ~e) * Subst(Antiderivative(1 / ((6 * (~A) ^ 3 * ~B * ~d + 3 * (~A) ^ 4 * ~e) - ~a * ~e * (~x) ^ 2), ~x), ~x, (~A + ~B * ~x) ^ 2 / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B], ~x) && (NeQ(~B * ~d - ~A * ~e, 0) && (EqQ((~c) ^ 2 * (~d) ^ 6 + ~a * (~e) ^ 4 * (13 * ~c * (~d) ^ 2 + ~b * (~e) ^ 2), 0) && (EqQ((~b) ^ 2 * (~e) ^ 4 - 12 * ~c * (~d) ^ 2 * (~c * (~d) ^ 2 - ~b * (~e) ^ 2), 0) && EqQ(4 * ~A * ~c * ~d * ~e + ~B * (2 * ~c * (~d) ^ 2 - ~b * (~e) ^ 2), 0))))
    @apply_utils Antiderivative(~Px / ((~d + ~(e') * ~x) * sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4)), ~x) => With([A = Coeff(~Px, ~x, 0), B = Coeff(~Px, ~x, 1), C = Coeff(~Px, ~x, 2), D = Coeff(~Px, ~x, 3)], Antiderivative((~x * ((B * ~d - A * ~e) + (~d * D - C * ~e) * (~x) ^ 2)) / (((~d) ^ 2 - (~e) ^ 2 * (~x) ^ 2) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)), ~x) + Antiderivative(((A * ~d + (C * ~d - B * ~e) * (~x) ^ 2) - D * ~e * (~x) ^ 4) / (((~d) ^ 2 - (~e) ^ 2 * (~x) ^ 2) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (PolyQ(~Px, ~x) && (LeQ(Expon(~Px, ~x), 3) && NeQ(~c * (~d) ^ 4 + ~b * (~d) ^ 2 * (~e) ^ 2 + ~a * (~e) ^ 4, 0)))
    @apply_utils Antiderivative(~Px / ((~d + ~(e') * ~x) * sqrt(~a + ~(c') * (~x) ^ 4)), ~x) => With([A = Coeff(~Px, ~x, 0), B = Coeff(~Px, ~x, 1), C = Coeff(~Px, ~x, 2), D = Coeff(~Px, ~x, 3)], Antiderivative((~x * ((B * ~d - A * ~e) + (~d * D - C * ~e) * (~x) ^ 2)) / (((~d) ^ 2 - (~e) ^ 2 * (~x) ^ 2) * sqrt(~a + ~c * (~x) ^ 4)), ~x) + Antiderivative(((A * ~d + (C * ~d - B * ~e) * (~x) ^ 2) - D * ~e * (~x) ^ 4) / (((~d) ^ 2 - (~e) ^ 2 * (~x) ^ 2) * sqrt(~a + ~c * (~x) ^ 4)), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (PolyQ(~Px, ~x) && (LeQ(Expon(~Px, ~x), 3) && NeQ(~c * (~d) ^ 4 + ~a * (~e) ^ 4, 0)))
    @apply_utils Antiderivative((~Px * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~(p')) / (~d + ~(e') * ~x), ~x) => ~d * Antiderivative((~Px * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p) / ((~d) ^ 2 - (~e) ^ 2 * (~x) ^ 2), ~x) - ~e * Antiderivative((~x * ~Px * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p) / ((~d) ^ 2 - (~e) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (PolyQ(~Px, ~x) && IntegerQ(~p + 1 / 2))
    @apply_utils Antiderivative((~Px * (~a + ~(c') * (~x) ^ 4) ^ ~(p')) / (~d + ~(e') * ~x), ~x) => ~d * Antiderivative((~Px * (~a + ~c * (~x) ^ 4) ^ ~p) / ((~d) ^ 2 - (~e) ^ 2 * (~x) ^ 2), ~x) - ~e * Antiderivative((~x * ~Px * (~a + ~c * (~x) ^ 4) ^ ~p) / ((~d) ^ 2 - (~e) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (PolyQ(~Px, ~x) && IntegerQ(~p + 1 / 2))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.2 Trinomial products/1.2.2 Quartic/1.2.2.8 P(x) (d+e x)^q (a+b x^2+c x^4)^p.jl


text:
using Documenter, BridgeSDEInference

makedocs(
    modules = [BridgeSDEInference],
    format = Documenter.HTML(
        #prettyurls = get(ENV, "CI", nothing) == "true"
        mathengine = Documenter.MathJax(Dict(:TeX => Dict(
    :equationNumbers => Dict(:autoNumber => "AMS"),
    :Macros => Dict(
        :ket => ["|#1\\rangle", 1],
        :bra => ["\\langle#1|", 1],
    ),
    )))

    ),
    sitename="BridgeSDEInference.jl",
    pages = [
        "Home" => "index.md",
        "Tutorial" => Any[
            "Workflow" => joinpath("man", "overview.md"),
            "Defining a diffusion" => joinpath("man", "model_definition.md"),
            "Data generation" => joinpath("man", "generate_data.md"),
            "Diffusion Setup" => joinpath("man", "diffusion_setup.md"),
            "MCMC Setup" => joinpath("man", "mcmc_setup.md"),
            "Running the sampler" => joinpath("man", "run.md"),
            "Querying the results" => joinpath("man", "query.md"),
            ],
        "Examples" => Any[
            "FitzHugh-Nagumo model" => joinpath("examples", "fitzhugh_nagumo.md"),
            "Jansen-Rit model" => joinpath("examples", "jansen_rit.md"),
            "Lorenz63 system" => joinpath("examples", "lorenz63.md"),
            "Lorenz96 system" => joinpath("examples", "lorenz96.md"),
            "Prokaryotic autoregulatory gene network" => joinpath("examples", "prokaryote.md"),
            "Sine diffusion" => joinpath("examples", "sine.md")
        ],
        "Visualisation tools" => joinpath("vis", "visualisation.md"),
        "Extras" => Any[
            "Generic MCMC" => joinpath("extras", "generic_mcmc.md"),
            "First passage times" => joinpath("extras", "first_passage_times.md"),
            "Blocking" => joinpath("extras", "blocking.md")
            ],
        "Adaptations" => Any[
            "Adaptive proposals" => joinpath("adaptations", "adaptive_proposals.md"),
            "Fusion" => joinpath("adaptations", "fusion.md")
        ]
    ],
)

deploydocs(
    repo = "github.com/mmider/BridgeSDEInference.jl.git",
)


text:
include("test_common.jl")

@testset "FastBEAST.jl" begin
    include("test_boundingbox.jl")
    include("test_boxtree.jl")
    include("test_kmeans.jl")
    include("test_skeletons.jl")
    include("test_aca.jl")
    include("test_hmatrix.jl")
    include("test_beast.jl")
    include("test_beast_efie.jl")
end

max_stars_repo_name: sbadrian/FastBEAST

max_issues_repo_name: sbadrian/FastBEAST

max_forks_repo_name: sbadrian/FastBEAST


text:
#############################################################################
#############################################################################
#
# This file implements polynomial division
#
#############################################################################
#############################################################################

"""  Modular algorithm.
f divide by g

f = q*g + r

p is a prime
"""
function divide(num::PolynomialModP, den::PolynomialModP)
    @assert num.prime == den.prime "Both polynomials must share the same prime divisor"
    f, g = mod(num.terms,num.prime), mod(den.terms,den.prime)
    degree(f) < degree(num.terms) && return nothing
    iszero(g) && throw(DivideError())
    q = Polynomial()
    prev_degree = degree(f)
    while degree(f) ≥ degree(g)
        h = Polynomial( (leading(f) ÷ leading(g))(num.prime) )  #syzergy
        f = mod((f - h*g), num.prime)
        q = mod((q + h), num.prime)
        prev_degree == degree(f) && break
        prev_degree = degree(f)
    end
    @assert iszero( mod((num.terms  - (q*g + f)),num.prime))
    return q, f
end

"""
The quotient from polynomial division. Returns a function of an integer.
"""

÷(num::PolynomialModP, den::PolynomialModP)  = first(divide(num,den))

"""
The remainder from polynomial division. Returns a function of an integer.
"""

rem(num::PolynomialModP, den::PolynomialModP)  = last(divide(num,den))


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    GeoStats

`ScientificTypes.jl` convention adopted in the GeoStats.jl project.
"""
struct GeoStats <: Convention end

# Basic values
scitype(::Integer,        ::GeoStats) = Count
scitype(::AbstractFloat,  ::GeoStats) = Continuous
scitype(::AbstractString, ::GeoStats) = Textual

# Categorical values
function scitype(c::CategoricalValue, ::GeoStats)
  L = length(levels(c))
  ifelse(isordered(pool(c)), OrderedFactor{L}, Multiclass{L})
end
function scitype(A::CategoricalArray{T,N}, ::GeoStats) where {T,N}
  L = length(levels(A))
  S = ifelse(isordered(A), OrderedFactor{L}, Multiclass{L})
  T >: Missing && (S = Union{S,Missing})
  AbstractArray{S,N}
end

# Scitype for fast array broadcasting
Scitype(::Type{<:Integer},        ::GeoStats) = Count
Scitype(::Type{<:AbstractFloat},  ::GeoStats) = Continuous
Scitype(::Type{<:AbstractString}, ::GeoStats) = Textual


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

using Markdown
using InteractiveUtils

# ╔═╡ e3311240-8431-11eb-1d8d-61736f938fc4
md"""

# Pluto.jl + fastpages!

> converting pluto.jl notebook to jupyter and cleaning them to host as a fastpages blog post

- toc: true
- categories: [blog post, fastai, fastpages]
- comments: true
"""

# ╔═╡ 17224d40-8431-11eb-07db-d5e5f18760cb
md"""

# Pluto.jl + Fastpages

Hi all! Today's post is for people who want to use Julia to write blog posts. I expect that you must have heard about this cool thing by [fast.ai](https://github.com/fastai) folks called [fastpages](https://fastpages.fast.ai/)

## What is fastpages and why to use it?

> An easy to use blogging platform, with support for Jupyter notebooks, Word docs, and Markdown.

Whole of my blog is made using fastpages. It takes just a few minutes to setup a whole blog using fastpages. Now, the best feature of fastpages is that you can use Jupyter notebooks to write a post and fastpages will automatically convert them into a blog post on your blog. Now this is seriously cool! Not just that you can write posts using regular markdown files, and also word document. fastpages will automatically render them as posts. You can refer their repo for more details: [fastpages](https://github.com/fastai/fastpages).

## Julia + fastpages

I have been trying julia out for few days now and instead of Python I wanted to use Julia to write blog posts. You can install the IJulia package and use Jupyter notebooks to write blog posts just like you do it in python. However, I would like to bring to your notice another very nice tool which I found about during JuliaCon 2020. I have used it a lot now and I can tell you that it addresses a lot of pitfalls of Jupyter Notebooks. It kind of changes the game. Currently, it only supports Julia and not Python. It's called [Pluto.jl](https://github.com/fonsp/Pluto.jl).

## Pluto.jl

Pluto.jl is:

> Simple reactive notebooks for Julia

Let me explain it to you what **"reactivity"** is. In Jupyter notebook, say you have used a function which you haven't imported. Naturally, as you run the cell it would raise an error. To fix this error, you go all the way to the top where your imports are, import the particular function and re-run all the cells where the function is used. This is where Pluto.jl shines over Jupyter notebook. In Pluto as you import the function all the cells which depend on it are re-run automatically. It's pretty smart in figuring out dependencies. This is a very powerful feature.

## Pluto.jl + fastpages!

Finally we look at how we can use Pluto.jl to write blog posts. Now, fastpages requires the notebooks to be Ipython notebooks,to be able to render them as posts. So, what do we do? We will convert the Pluto.jl notebooks to Ipython notebook first and then use them.

To do this, there's this tool made by the creator of Pluto.jl we can use -> [pluto-jl-jupyter-conversion](https://observablehq.com/@olivier_plas/pluto-jl-jupyter-conversion). Pluto saves your notebooks as .jl files. All you have to do is upload your Pluto.jl notebook and it will convert it into Jupyter notebook after which u can download it.

There's only one issue. The converted notebook has all the cells as code cells. Even the cells we used as markdown cells are converted to code cells. We need to fix this!
I have made a simple utility to solve this issue -> [pluto-jl-jupyter-conversion-cleaner](https://github.com/Abhiswain97/pluto-jl-jupyter-conversion-cleaner)

So, finally these are the steps you need to follow to convert your Pluto.jl notebook into a Jupyter notebook you can use with fastpages

### Steps

1. Use [pluto-jl-jupyter-conversion](https://observablehq.com/@olivier_plas/pluto-jl-jupyter-conversion) to convert the Pluto.jl notebook to a Ipython notebook.

2. Install all the requirements:

   For windows: `run.bat` \
   For Linux: `sh run.sh` \

3. Clone repo and use it:
   ```
   git clone https://github.com/Abhiswain97/pluto-jl-jupyter-conversion-cleaner.git
   cd pluto-jl-jupyter-conversion-cleaner
   python cleaner.py --path <your converted .ipynb file>
   ```

The converted file is saved as file-name-cleaned.ipynb in the same folder as your original Jupyter notebook.

Now, all you got to do is make it in the format fastpages expects, it has to be saved as: **yyyy-mm-dd-name.ipynb**. There are various other options you can specify. Have a look at their docs at: [fastpages](https://fastpages.fast.ai/)

Yay! we finally have a Jupyter notebook we can use with fastpages!

If you made it this far, then let me tell you something, the post you're reading right now is written using Pluto.jl & fastpages!
"""

# ╔═╡ Cell order:
# ╠═e3311240-8431-11eb-1d8d-61736f938fc4
# ╠═17224d40-8431-11eb-07db-d5e5f18760cb


text:
function ITensors.OpSum(::Model"hubbard", n1, n2; t, U, V)
  opsum = OpSum()
  opsum += -t, "Cdagup", n1, "Cup", n2
  opsum += -t, "Cdagup", n2, "Cup", n1
  opsum += -t, "Cdagdn", n1, "Cdn", n2
  opsum += -t, "Cdagdn", n2, "Cdn", n1
  if U ≠ 0
    opsum += U, "Nupdn", n1
  end
  if V ≠ 0
    opsum += V, "Ntot", n1, "Ntot", n2
  end
  return opsum
end

function ITensors.MPO(::Model"hubbard", s; t, U, V)
  N = length(s)
  opsum = OpSum()
  for n in 1:(N - 1)
    n1, n2 = n, n + 1
    opsum .+= -t, "Cdagup", n1, "Cup", n2
    opsum .+= -t, "Cdagup", n2, "Cup", n1
    opsum .+= -t, "Cdagdn", n1, "Cdn", n2
    opsum .+= -t, "Cdagdn", n2, "Cdn", n1
    if V ≠ 0
      opsum .+= V, "Ntot", n1, "Ntot", n2
    end
  end
  if U ≠ 0
    for n in 1:N
      opsum .+= U, "Nupdn", n
    end
  end
  return splitblocks(linkinds, MPO(opsum, s))
end

function ITensors.ITensor(model::Model"hubbard", s; kwargs...)
  return prod(MPO(model, s; kwargs...))
end

"""
@article{PhysRevB.6.930,
  title = {Magnetic Susceptibility at Zero Temperature for the One-Dimensional Hubbard Model},
  author = {Shiba, H.},
  journal = {Phys. Rev. B},
  volume = {6},
  issue = {3},
  pages = {930--938},
  numpages = {0},
  year = {1972},
  month = {Aug},
  publisher = {American Physical Society},
  doi = {10.1103/PhysRevB.6.930},
  url = {https://link.aps.org/doi/10.1103/PhysRevB.6.930}
}
"""
function reference(::Model"hubbard", ::Observable"energy"; U, Npoints=50)
  f(x) = 1 / π * sum(n -> (-1)^(n + 1) * (2n) / (x^2 + (2n)^2), 1:10000)
  function m(i, j)
    return (8π) / U * (2π) / Npoints * f(
      4 / U * (sin(-π + (2(i - 1) * π) / Npoints) - sin(-π + (2(j - 1) * π) / Npoints))
    )
  end
  matrix = [m(i, j) for i in 1:Npoints, j in 1:Npoints]
  for i in 1:Npoints
    matrix[i, :] *= cos(-π + (i - 1) * (2π) / Npoints)
  end
  ones_vec = ones(Npoints)
  d = Diagonal(2π .* ones_vec)
  v, _ = linsolve(d - matrix, ones_vec)
  g(i) = -2 * (2π) / Npoints * cos(-π + (i - 1) / Npoints * 2π) * v[i]
  return sum(g, 1:Npoints)
end


text:
###############
# Arithmetics
###############

# We can not specialize (*/∘/...)(ops::AbstractOperator...), since too general
(*)(op1::AbstractOperator, op2::AbstractOperator) = compose(op2, op1)
(∘)(op1::AbstractOperator, op2::AbstractOperator) = (*)(op1, op2)

# make times (*) a synonym for applying the operator
(*)(op::AbstractOperator, object) = apply(op, object)

mul!(y::AbstractVector, A::DictionaryOperator, x::AbstractVector) = apply!(A, tocoefficientformat(y,dest(A)), tocoefficientformat(x,src(A)))
mul!(y::AbstractMatrix, A::DictionaryOperator, x::AbstractMatrix) = apply_multiple!(A, y, x)

# LazyArrays does not work as well as composite operator.
# include("../util/lazyarrays.jl")
# (∘)(op1::ArrayOperator, op2::ArrayOperator) = ArrayOperator(Mul(unsafe_matrix(op1), unsafe_matrix(op2)), src(op2), dest(op1))
# (*)(op1::ArrayOperator, op2::ArrayOperator) = ArrayOperator(_checked_mul(unsafe_matrix(op1), unsafe_matrix(op2)), src(op2), dest(op1))
# (+)(op1::ArrayOperator, op2::ArrayOperator) = ArrayOperator(_checked_sum(unsafe_matrix(op1), unsafe_matrix(op2)), src(op2), dest(op1))
# (-)(op1::ArrayOperator, op2::ArrayOperator) = ArrayOperator(_checked_sum(unsafe_matrix(op1), -unsafe_matrix(op2)), src(op2), dest(op1))

# Hacks for some operators

for (OP) in (:DiagonalOperator, :CirculantOperator, :ScalingOperator)
    @eval begin
        (*)(op1::$OP, op2::$OP) = ArrayOperator(_checked_mul(unsafe_matrix(op1), unsafe_matrix(op2)), src(op2), dest(op1))
        (+)(op1::$OP, op2::$OP) = ArrayOperator(_checked_sum(unsafe_matrix(op1), unsafe_matrix(op2)), src(op2), dest(op1))
        (-)(op1::$OP, op2::$OP) = ArrayOperator(_checked_sum(unsafe_matrix(op1), -unsafe_matrix(op2)), src(op2), dest(op1))
    end
end

(*)(a::Number, op::ArrayOperator) = ArrayOperator(a*unsafe_matrix(op), src(op), dest(op))
(*)(op::ArrayOperator, a::Number) = ArrayOperator(a*unsafe_matrix(op), src(op), dest(op))

(*)(a::Number, op::DictionaryOperator) = ScalingOperator(dest(op), a) * op
(*)(op::DictionaryOperator, a::Number) = op * ScalingOperator(src(op), a)

for op in (:unchecked_mul,:unchecked_sum, :∘)
    # copied from base such that is is sufficient to implement the two argument problem
    @eval ($op)(a, b, c, xs...) = Base.afoldl($op, ($op)(($op)(a,b),c), xs...)
end

# Default implementation, not typestable...
function _checked_mul(a1::AbstractArray, a2::AbstractArray)
    a = a1*a2
    if isefficient(a1) && isefficient(a2) && !isefficient(a) && !(ignore_mul_message(a1, a2, a))
            @warn "The multiplication, $(typeof(a)), of efficient $(typeof(a1)) and efficient $(typeof(a2)) is not efficient.\n
        Consider implementing *(a1::$(typeof(a1)), a2::$(typeof(a2)))"
    end
    a
end

function _checked_sum(a1::AbstractArray, a2::AbstractArray)
    a = a1+a2
    if isefficient(a1) && isefficient(a2) && !isefficient(a) && !(ignore_sum_message(a1, a2, a))
        @warn "The sum, $(typeof(a)), of efficient $(typeof(a1)) and efficient $(typeof(a2)) is not efficient.\n
        Consider implementing +(a1::$(typeof(a1)), a2::$(typeof(a2)))"
    end
    a
end

ignore_mul_message(a, b, c) = false
ignore_sum_message(a, b, c) = false


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: iot_data_plane
using AWS.Compat
using AWS.UUIDs

"""
    delete_thing_shadow(thing_name)
    delete_thing_shadow(thing_name, params::Dict{String,<:Any})

Deletes the shadow for the specified thing. For more information, see DeleteThingShadow in
the AWS IoT Developer Guide.

# Arguments
- `thing_name`: The name of the thing.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"name"`: The name of the shadow.
"""
delete_thing_shadow(thingName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_data_plane("DELETE", "/things/$(thingName)/shadow"; aws_config=aws_config)
delete_thing_shadow(thingName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_data_plane("DELETE", "/things/$(thingName)/shadow", params; aws_config=aws_config)

"""
    get_thing_shadow(thing_name)
    get_thing_shadow(thing_name, params::Dict{String,<:Any})

Gets the shadow for the specified thing. For more information, see GetThingShadow in the
AWS IoT Developer Guide.

# Arguments
- `thing_name`: The name of the thing.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"name"`: The name of the shadow.
"""
get_thing_shadow(thingName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_data_plane("GET", "/things/$(thingName)/shadow"; aws_config=aws_config)
get_thing_shadow(thingName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_data_plane("GET", "/things/$(thingName)/shadow", params; aws_config=aws_config)

"""
    list_named_shadows_for_thing(thing_name)
    list_named_shadows_for_thing(thing_name, params::Dict{String,<:Any})

Lists the shadows for the specified thing.

# Arguments
- `thing_name`: The name of the thing.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"nextToken"`: The token to retrieve the next set of results.
- `"pageSize"`: The result page size.
"""
list_named_shadows_for_thing(thingName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_data_plane("GET", "/api/things/shadow/ListNamedShadowsForThing/$(thingName)"; aws_config=aws_config)
list_named_shadows_for_thing(thingName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_data_plane("GET", "/api/things/shadow/ListNamedShadowsForThing/$(thingName)", params; aws_config=aws_config)

"""
    publish(topic)
    publish(topic, params::Dict{String,<:Any})

Publishes state information. For more information, see HTTP Protocol in the AWS IoT
Developer Guide.

# Arguments
- `topic`: The name of the MQTT topic.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"payload"`: The state information, in JSON format.
- `"qos"`: The Quality of Service (QoS) level.
"""
publish(topic; aws_config::AbstractAWSConfig=global_aws_config()) = iot_data_plane("POST", "/topics/$(topic)"; aws_config=aws_config)
publish(topic, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_data_plane("POST", "/topics/$(topic)", params; aws_config=aws_config)

"""
    update_thing_shadow(payload, thing_name)
    update_thing_shadow(payload, thing_name, params::Dict{String,<:Any})

Updates the shadow for the specified thing. For more information, see UpdateThingShadow in
the AWS IoT Developer Guide.

# Arguments
- `payload`: The state information, in JSON format.
- `thing_name`: The name of the thing.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"name"`: The name of the shadow.
"""
update_thing_shadow(payload, thingName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_data_plane("POST", "/things/$(thingName)/shadow", Dict{String, Any}("payload"=>payload); aws_config=aws_config)
update_thing_shadow(payload, thingName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_data_plane("POST", "/things/$(thingName)/shadow", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("payload"=>payload), params)); aws_config=aws_config)


text:
function processinterval(   𝑖::Int,
                            𝓟::Vector{T},
                            sample_position_template::Vector{T},
                            𝑒::Function,
                            L::Int,
                            verbose_flag) where T
    @assert 1 < 𝑖 <= length(𝓟)

    a = 𝓟[𝑖-1]
    b = 𝓟[𝑖]

    # build approximation of f over current interval.
    X = samplepoints( sample_position_template, a, b )

    𝑒_itp_X = xx->Barycentric2nditp(xx, X, 𝑒.(X))

    # fit Chebyshevinterpolator over this approximation.
    d, ν = getChebyshevinterpolator(𝑒_itp_X, L, a, b)

    𝑒_extrema_positions_X = scalefromChebyinterval.(findChebyitpextremawithchk(d), a, b)

    # d and X can be used for plotting or debugging purposes.
    return 𝑒_extrema_positions_X, d, X
end

# sample_position_template is between [-1,1]. For each interval, rescale it to the interval at hand.
function getextremafromroots(  𝑒::Function,
                                L::Int,
                                𝓟_in::Vector{T},
                                sample_position_template::Vector{T},
                                tol_𝓧_spacing, # was defaulted to convert(T,1e-6/length(𝓟_in)
                                tol_derivative_zero,
                                verbose_flag = true) where T
    𝓟 = Vector{T}(undef,length(𝓟_in)+2)
    𝓟[1] = -one(T)
    𝓟[2:end-1] = 𝓟_in
    𝓟[end] = one(T)

    makenotclose!(𝓟, convert(T,1e-14/length(𝓟_in))) # just to be sure positions are spread out enough.

    # prepare outputs.
    out = Vector{T}(undef,0)
    d_set = Vector{Vector{T}}(undef,length(𝓟)-1)
    X_set = Vector{Vector{T}}(undef,length(𝓟)-1)

    # work on intervals.
    for i = 2:length(𝓟)

        # build approximation of f over current interval.
        X = samplepoints( sample_position_template,𝓟[i-1],𝓟[i] )

        𝑒_itp_X = xx->Barycentric2nditp(xx, X, 𝑒.(X))

        # fit Chebyshevinterpolator over this approximation.
        d, ν = getChebyshevinterpolator(𝑒_itp_X, L, 𝓟[i-1], 𝓟[i])

        if verbose_flag
            Printf.@printf("working on %d-th interval\n",i-1)
        end
        𝑒_extrema_positions_X = scalefromChebyinterval.(findChebyitpextremawithchk(d), 𝓟[i-1], 𝓟[i])

        push!(out, 𝑒_extrema_positions_X...)

        # for plotting or debugging.
        d_set[i-1] = d
        X_set[i-1] = X
    end

    # remove extrema search results that are too close to zero.
    #   These are probably incorrect results, since the local extremas
    #   should be greater or equal to the level error, which is non-zero
    #   unless the target function is in the space of Chebyshev approximations.
    while sign(𝑒(out[1])) == zero(T) && !isempty(out)
        popfirst!(out)
    end

    return out, d_set, X_set
end

function getextremafromrootsparallel(𝑒::Function,
                                L::Int,
                                𝓧::Vector{T},
                                sample_position_template::Vector{T},
                                tol_𝓧_spacing, # was defaulted to convert(T,1e-6/length(𝓟_in)
                                tol_derivative_zero,
                                verbose_flag = true) where T

    ## construct the boundary point sets that will dictate the piece-wise proxy intervals.
    𝓑 = copy(𝓧)
    if !isnumericallyclose(-one(T), 𝓧[1])
        insert!(𝓑,1,-one(T))
    end
    if !isnumericallyclose(one(T), 𝓧[end])
        push!(𝓑,one(T))
    end

    # work on intervals.
    workerfunc = xx->processinterval(xx, 𝓑, sample_position_template, 𝑒, L, verbose_flag)
    sol = pmap(workerfunc, collect(2:length(𝓑)) )

    # package up outputs.
    A = sol
    k = 1
    len = sum( length(A[i][k]) for i = 1:length(A))
    out = Vector{T}(undef,len)

    st = 1
    for i = 1:length(A)
        fin = st + length(A[i][k]) -1
        out[st:fin] = A[i][k]

        st = fin + 1
    end

    d_set = collect( A[i][2] for i = 1:length(A) )
    X_set = collect( A[i][3] for i = 1:length(A) )

    return out, d_set, X_set
end

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
using Interpolations: degree,
                      getknots,
                      coefficients,
                      itpflag,
                      DimSpec,
                      AbstractLanczos,
                      lbounds,
                      ubounds

@testset "Lanczos" begin

@testset "Lanczos($N)" for N in 2:4
    X = 1:100
    X = [X; reverse(X)[2:end]]
    itp = interpolate(X, Lanczos(N))

    # properties
    @test getknots(itp) == axes(itp) == axes(X)
    @test coefficients(itp)[axes(itp)...] == X
    @test itpflag(itp) isa Lanczos{N}
    @test size(itp) == size(X)
    @test lbounds(itp) == (firstindex(X),)
    @test ubounds(itp) == (lastindex(X),)

    # recovers input points exactly
    @test itp.(X) == X

    @test 1 < itp(1.5) < 2
    @test 99 < itp(99.5) < 100

    # symmetry check
    interpolant = itp.(X)
    @test interpolant ≈ reverse(interpolant)

end

@testset "Lanczos OpenCV4" begin
    X = 1:100
    X = [X; reverse(X)[2:end]]
    itp = interpolate(X, Lanczos4OpenCV())

    # properties
    @test getknots(itp) == axes(itp) == axes(X)
    @test coefficients(itp)[axes(itp)...] == X
    @test itpflag(itp) isa Lanczos4OpenCV
    @test size(itp) == size(X)
    @test lbounds(itp) == (firstindex(X),)
    @test ubounds(itp) == (lastindex(X),)

    # recovers input points exactly
    @test itp.(X) == X

    @test 1 < itp(1.5) < 2
    @test 99 < itp(99.5) < 100

    # symmetry check
    interpolant = itp.(X)
    @test interpolant ≈ reverse(interpolant)

end

end


text:
import GeometryBasics
import Observables
import Makie

export Triad
export update

"""
    Represents a triad.

fields: q, length, width, color and observable.
"""
mutable struct Triad <: Sprite
    q::Biquaternion
    length::Float64
    width::Int
    color::Array{Symbol}
    observable::Observables.Observable
end

"""
    Triad(q, scene)

Construct a Triad with the given configuration `q`, `scene`, and the optional arguments:
`length`, `width` and `color`.
"""
function Triad(q::Biquaternion,
               scene::Makie.Scene;
               length::Float64 = 1.0,
               width::Int = 5,
               color::Array{Symbol} = [:red, :green, :blue])
    triad = constructtriad(q, length = length)
    observable = Observables.Observable(triad)
    array = Makie.@lift begin
        s = map(x -> GeometryBasics.Point3f(vec(x)), $observable)
        [s[1] => s[2], s[3] => s[4], s[5] => s[6]]
    end
    Makie.linesegments!(scene,
                                   array,
                                   color = color,
                                   linewidth = width)
    Triad(q, length, width, color, observable)
end

"""
    update(triad, q)

Update a Triad by changing its observable with the given `triad` and configuration `q`.
"""
function update(triad::Triad, q::Biquaternion)
    triad.q = q
    triad.observable[] = constructtriad(triad.q, length = triad.length)
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
function transform(input::Dict)
    output = []
    for i = input
        for j = map(lowercase, i[2])
            pair = (j, i[1])
            push!(output, pair)
        end
    end
    Dict(output)
end


text:
using PhaseSlopeIndex
using Test

@testset "PhaseSlopeIndex.jl" begin
    # Write your tests here.
end


text:
export Hermite,GaussWeight

#TODO: Add general lines

"""
`Hemite(L)` represents `H_k(sqrt(L) * x)` where `H_k` are Hermite polynomials.
`Hermite()` is equivalent to `Hermite(1.0)`.
"""
struct Hermite{T} <: PolynomialSpace{Line{false,Float64},Float64}
    L::T
end
Hermite() = Hermite(1.0)

domain(::Hermite) = Line()
canonicalspace(H::Hermite) = Hermite()
spacescompatible(::Hermite,::Hermite) = true #TODO:L
canonicaldomain(H::Hermite) = Line()
tocanonical(H::Hermite,x) = x

#####
# recα/β/γ are given by
#       x p_{n-1} =γ_n p_{n-2} + α_n p_{n-1} +  p_n β_n
#####

recα(::Type,::Hermite,k) = 0;
recβ(::Type,H::Hermite,k) = 0.5 / sqrt(H.L);
recγ(::Type,H::Hermite,k) = (k-1) / sqrt(H.L);
recA(::Type,H::Hermite,k) = 2 * sqrt(H.L);
recB(::Type,::Hermite,k) = 0;
recC(::Type,::Hermite,k) = 2k;

Derivative(H::Hermite,order) = ConcreteDerivative(H,order)

bandwidths(D::ConcreteDerivative{H}) where {H<:Hermite} = 0,D.order
rangespace(D::ConcreteDerivative{H}) where {H<:Hermite} = domainspace(D)
getindex(D::ConcreteDerivative{H},k::Integer,j::Integer) where {H<:Hermite} =
        j==k+D.order ? one(eltype(D))*2^D.order*pochhammer(k,D.order) : zero(eltype(D))

function hermitep(r::AbstractRange,x::Number)
    n = r[end] + 1

    T = typeof(x)
    H = Hermite()
    if n ≤ 2
        v = [1.,2x]
    else
        v=Array{promote_type(Float64,typeof(x))}(undef, n)  # x may be complex
        v[1] = 1.
        v[2] = 2x

        for k = 2:n-1
            v[k + 1] = ((x - recα(T,H,k)) * v[k] - recγ(T,H,k) * v[k - 1]) / recβ(T,H,k)
        end
    end
    v[r.+1]
end
hermitep(n::Integer, v::Number) = hermitep(n:n,v)[1]

# exp(-Lx^2)
struct GaussWeight{S,T} <: WeightSpace{S,Line{Float64},Float64}
    space::S
    L::T
end

GaussWeight(H::Hermite) = GaussWeight(H,H.L)
GaussWeight() = GaussWeight(Hermite())

"""
`GaussWeight(Hermite(L), L)` is a space spanned by `exp(-Lx²) * H_k(sqrt(L) * x)`
where `H_k(x)`'s are Hermite polynomials.

`GaussWeight()` is equivalent to `GaussWeight(Hermite(), 1.0)` by default.
"""
Fun(::typeof(identity), sp::Hermite) = Fun(sp,[0.,0.5])
Fun(::typeof(identity), sp::GaussWeight) = Fun(identity, sp.space)

spacescompatible(a::GaussWeight,b::GaussWeight)=spacescompatible(a.space,b.space)&&isapprox(a.L,b.L)

function Derivative(sp::GaussWeight,k::Integer)
   if k == 1
        x = Multiplication(Fun(identity,sp.space),sp.space)
        D = Derivative(sp.space)
        D2 = D - (2sp.L)x
        DerivativeWrapper(SpaceOperator(D2,sp,GaussWeight(rangespace(D2),sp.L)),1)
    else
        D = Derivative(sp)
        DerivativeWrapper(TimesOperator(Derivative(rangespace(D),k-1).op,D.op),k)
    end
end

weight(H::GaussWeight,x) = exp(-H.L * x^2)

function Base.sum(f::Fun{GaussWeight{H,T}}) where {H<:Hermite,T}
    @assert space(f).space.L==space(f).L  # only implemented with matching weight
    f.coefficients[1]*sqrt(convert(T,π))/sqrt(space(f).L)
end

include("hermitetransform.jl")

function Multiplication(f::Fun{H}, S::GaussWeight{H}) where H<:Hermite
    M = Multiplication(f, S.space)
    rs = rangespace(M)
    MultiplicationWrapper(f, SpaceOperator(M, S, GaussWeight(rs, S.L)))
end

function Multiplication(f::Fun{GaussWeight{H,T}}, S::Hermite) where {H<:Hermite,T}
    M = Multiplication(Fun(space(f).space, f.coefficients), S)
    rs = rangespace(M)
    MultiplicationWrapper(f, SpaceOperator(M, S, GaussWeight(rs, space(f).L)))
end

function integrate(f::Fun{GW}) where GW<:GaussWeight{H} where H<:Hermite
    n = length(f.coefficients);
    L = space(f).L;
    if L == 0
        return Fun(Hermite(space(f).space.L), [0; f.coefficients] ./ [1; 2:2:2n])
    elseif space(f).space.L != L
        throw(ArgumentError("`integrate` is applicable if only parameters of GaussWeight and Hermite are equal."))
    else
        if n == 0
            return Fun(0)
        elseif f.coefficients[1] == 0
            return Fun(GaussWeight(Hermite(L), L), f.coefficients[2:end] / -sqrt(L))
        else
            g = Fun(GaussWeight(Hermite(L), L), f.coefficients[2:end] / -sqrt(L));
            f₀ = Fun(GaussWeight(Hermite(L), L), [f.coefficients[1]]);
            f₀ = Fun(f₀, Chebyshev(Line()));
            g₀ = integrate(f₀);
            g₀ = g₀ - last(g₀);
            return g + g₀
        end
    end
end


text:
using OrdinaryDiffEq, DiffEqBase
using DynamicalSystemsBase: DEFAULT_DIFFEQ_KWARGS, _get_solver
using Roots: find_zero, A42, AlefeldPotraShi, Brent
export poincaresos, produce_orbitdiagram

const ROOTS_ALG = A42()

#####################################################################################
#                               Poincare Section                                    #
#####################################################################################
"""
    poincaresos(ds::ContinuousDynamicalSystem, plane, tfinal = 1000.0; kwargs...)
Calculate the Poincaré surface of section (also called Poincaré map) [1, 2]
of the given system with the given `plane`.
The system is evolved for total time of `tfinal`.

If the state of the system is ``\\mathbf{u} = (u_1, \\ldots, u_D)`` then the
equation intersecting the hyperplane is
```math
a_1u_1 + \\dots + a_Du_D = \\mathbf{a}\\cdot\\mathbf{u}=b
```
where ``\\mathbf{a}, b`` are the parameters that define the hyperplane.

In code, `plane` can be either:

* A `Tuple{Int, <: Number}`, like `(j, r)` : the planecrossing is defined
  as when the `j` variable of the system crosses the value `r`.
* An `AbstractVector` of length `D+1`. The first `D` elements of the
  vector correspond to ``\\mathbf{a}`` while the last element is ``b``.

Returns a [`Dataset`](@ref) of the points that are on the surface of section.

## Keyword Arguments
* `direction = 1` : Only crossings with `sign(direction)` are considered to belong to
  the surface of section. Positive direction means going from less than ``b``
  to greater than ``b``.
* `idxs = 1:dimension(ds)` : Optionally you can choose which variables to save.
  Defaults to the entire state.
* `Ttr = 0.0` : Transient time to evolve the system before starting
  to compute the PSOS.
* `warning = true` : Throw a warning if the Poincaré section was empty.
* `rootkw = (xrtol = 1e-6, atol = 1e-6)` : A `NamedTuple` of keyword arguments
  passed to `find_zero` from [Roots.jl](https://github.com/JuliaMath/Roots.jl).
* `diffeq...` : All other extra keyword arguments are propagated into `init`
  of DifferentialEquations.jl. See [`trajectory`](@ref) for examples.

## References
[1] : H. Poincaré, *Les Methods Nouvelles de la Mécanique Celeste*,
Paris: Gauthier-Villars (1892)

[2] : M. Tabor, *Chaos and Integrability in Nonlinear Dynamics: An Introduction*,
§4.1, in pp. 118-126, New York: Wiley (1989)

See also [`orbitdiagram`](@ref), [`produce_orbitdiagram`](@ref).
"""
function poincaresos(ds::CDS{IIP, S, D}, plane, tfinal = 1000.0;
    direction = +1, Ttr::Real = 0.0, warning = true, idxs = 1:D,
    rootkw = (xrtol = 1e-6, atol = 1e-6), diffeq...) where {IIP, S, D}

    _check_plane(plane, D)
    integ = integrator(ds; diffeq...)
    planecrossing = PlaneCrossing{D}(plane, direction > 0 )
    f = (t) -> planecrossing(integ(t))

    i = typeof(idxs) <: Int ? i : SVector{length(idxs), Int}(idxs...)

    data = _initialize_output(get_state(ds), i)

    poincare_cross!(data, integ,
                     f, planecrossing, tfinal, Ttr, i, rootkw)
    warning && length(data) == 0 && @warn PSOS_ERROR

    return Dataset(data)
end

function _initialize_output(u::S, i::Int) where {S}
    output = eltype(S)[]
end
function _initialize_output(u::S, i::SVector) where {S}
    s = u[i]
    output = typeof(s)[]
end

const PSOS_ERROR =
"the Poincaré surface of section did not have any points!"

struct PlaneCrossing{D, P}
    plane::P
    dir::Bool
end
PlaneCrossing{D}(p::P, b) where {P, D} = PlaneCrossing{D, P}(p, b)
function (hp::PlaneCrossing{D, P})(u::AbstractVector) where {D, P<:Tuple}
    @inbounds x = u[hp.plane[1]] - hp.plane[2]
    hp.dir ? x : -x
end
function (hp::PlaneCrossing{D, P})(u::AbstractVector) where {D, P<:AbstractVector}
    x = zero(eltype(u))
    @inbounds for i in 1:D
        x += u[i]*hp.plane[i]
    end
    @inbounds x -= hp.plane[D+1]
    hp.dir ? x : -x
end

function poincare_cross!(data, integ,
                          f, planecrossing, tfinal, Ttr, j, rootkw)

    Ttr != 0 && step!(integ, Ttr)

    side = planecrossing(integ.u)

    while integ.t < tfinal + Ttr
        while side < 0
            integ.t > tfinal + Ttr && break
            step!(integ)
            side = planecrossing(integ.u)
        end
        while side ≥ 0
            integ.t > tfinal + Ttr && break
            step!(integ)
            side = planecrossing(integ.u)
        end
        integ.t > tfinal + Ttr && break

        # I am now guaranteed to have `t` in negative and `tprev` in positive
        tcross = find_zero(f, (integ.tprev, integ.t), ROOTS_ALG; rootkw...)
        ucross = integ(tcross)

        push!(data, ucross[j])
    end
    return data
end

function _check_plane(plane, D)
    P = typeof(plane)
    L = length(plane)
    if P <: AbstractVector
        if L != D + 1
            throw(ArgumentError(
            "The plane for the `poincaresos` must be either a 2-Tuple or a vector of "*
            "length D+1 with D the dimension of the system."
            ))
        end
    elseif P <: Tuple
        if !(P <: Tuple{Int, Number})
            throw(ArgumentError(
            "If the plane for the `poincaresos` is a 2-Tuple then "*
            "it must be subtype of `Tuple{Int, Number}`."
            ))
        end
    else
        throw(ArgumentError(
        "Unrecognized type for the `plane` argument."
        ))
    end
end

#####################################################################################
#                            Produce Orbit Diagram                                  #
#####################################################################################
"""
    produce_orbitdiagram(ds::ContinuousDynamicalSystem, plane, i::Int,
                         p_index, pvalues; kwargs...)
Produce an orbit diagram (also called bifurcation diagram)
for the `i` variable(s) of the given continuous
system by computing Poincaré surfaces of section using `plane`
for the given parameter values (see [`poincaresos`](@ref)).

`i` can be `Int` or `AbstractVector{Int}`.
If `i` is `Int`, returns a vector of vectors. Else
it returns vectors of vectors of vectors.
Each entry are the points at each parameter value.

## Keyword Arguments
* `printparams::Bool = false` : Whether to print the parameter used
  during computation in order to keep track of running time.
* `direction, warning, Ttr, rootkw, diffeq...` :
  Propagated into [`poincaresos`](@ref).
* `u0 = get_state(ds)` : Initial condition. Besides a vector you can also give
  a vector of vectors such that `length(u0) == length(pvalues)`. Then each parameter
  has a different initial condition.

## Description
For each parameter, a PSOS reduces the system from a flow to a map. This then allows
the formal computation of an "orbit diagram" for the `i` variable
of the system, just like it is done in [`orbitdiagram`](@ref).

The parameter change is done as `p[p_index] = value` taking values from `pvalues`
and thus you must use a parameter container that supports this
(either `Array`, `LMArray`, dictionary or other).

See also [`poincaresos`](@ref), [`orbitdiagram`](@ref).
"""
function produce_orbitdiagram(
    ds::CDS{IIP, S, D},
    plane,
    idxs,
    p_index,
    pvalues;
    tfinal::Real = 100.0,
    direction = +1,
    printparams = false,
    warning = true,
    Ttr = 0.0,
    u0 = get_state(ds),
    rootkw = (xrtol = 1e-6, atol = 1e-6),
    diffeq...) where {IIP, S, D}

    i = typeof(idxs) <: Int ? idxs : SVector{length(idxs), Int}(idxs...)

    _check_plane(plane, D)
    typeof(u0) <: Vector{<:AbstractVector} && @assert length(u0)==length(p)

    integ = integrator(ds; diffeq...)
    planecrossing = PlaneCrossing{D}(plane, direction > 0 )
    f = (t) -> planecrossing(integ(t))
    p0 = ds.p[p_index]

    output = [_initialize_output(ds.u0, i) for k in 1:length(pvalues)]

    for (n, p) in enumerate(pvalues)
        integ.p[p_index] = p
        printparams && println("parameter = $p")

        if typeof(u0) <: Vector{<:AbstractVector}
            st = u0[n]
        else
            st = u0
        end

        reinit!(integ, st)
        poincare_cross!(output[n], integ,
                         f, planecrossing, tfinal, Ttr, i, rootkw)

        warning && length(output[n]) == 0 && @warn "For parameter $p $PSOS_ERROR"

    end
    # Reset the parameter of the system:
    ds.p[p_index] = p0
    return output
end


text:
using Conda

function conda_setup()
    Conda.add_channel("conda-forge")
    Conda.add("selenium")
end

function build()
    conda_setup()
end

build()

max_stars_repo_name: stivenroytman/Selenium

max_issues_repo_name: stivenroytman/Selenium

max_forks_repo_name: stivenroytman/Selenium


text:
TRUE = Node("TRUE", ASingletonStatement)
FALSE = Node("FALSE", ASingletonStatement)
UNKNOWN = Node("UNKNOWN", ASingletonStatement)

macro singleton_set(parent::Symbol, abstract_type::Symbol, name::Symbol)
    return quote
        abstract type $(esc(abstract_type)) <: $(esc(parent)) end
        $(esc(name)) = Node($(String(name)), $(esc(abstract_type)))
    end
end

@singleton_set ASingletonSet AΩ Ω
@singleton_set AΩ AStatements Statements
@singleton_set AΩ ASets Sets
# @singleton_set AΩ AFunctions Functions
abstract type AFunctions{T} <: AΩ where T end
Functions(set) = Node("Functions", AFunctions{typeof(set)})  # the set of functions whose ouput is an element of set
@singleton_set AΩ ANumbers Numbers

@singleton_set AΩ A∅ ∅
@singleton_set AΩ Aℝ ℝ

quick_range(::Type{<:ABasic}) = AΩ
quick_range(::Type{<:ALogicOperator}) = AStatements
quick_range(::Type{<:ASetOperator}) = ASets
quick_range(::Type{<:AExpressionOperator}) = ANumbers
quick_range(x::Type{<:AUnknown}) = x.parameters[1]
quick_range(::Type{<:ASingletonSet}) = ASets
quick_range(::Type{<:ASingletonStatement}) = AStatements

quick_range(x::Node) = quick_range(class(x))


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "ElTopoBuilder"
version = v"0.1.3"

# Collection of sources required to build ElTopoBuilder
sources = [
    "https://github.com/tysonbrochu/eltopo.git" =>
    "14b1d7cbd45def90cfce04d381e92c4fefc5fab7",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/eltopo/eltopo3d/

cat <<EOF > Makefile.local_defs
INCLUDE_PATH = -I. -I../common -I../common/newsparse -I../common/meshes -I../common/tunicate
DEPEND = g++ -D__LITTLE_ENDIAN__ -DUSE_FORTRAN_BLAS -DNO_GUI
CC = g++ -Wall -D__LITTLE_ENDIAN__ -DUSE_FORTRAN_BLAS -DNO_GUI -fPIC
RELEASE_FLAGS = -O3 -funroll-loops
DEBUG_FLAGS = -g
LINK = g++
LINK_LIBS = -lGL -lGLU -lglut destdir/lib/libopenblas.a
EOF

make depend
make release

cd $WORKSPACE
(cd srcdir/eltopo/eltopo3d && find . -name '*.h' -print | tar --create --files-from -) | (cd $WORKSPACE/destdir/include && tar xvfp -)
(cd srcdir/eltopo/common && find . -name '*.h' -print | tar --create --files-from -) | (cd $WORKSPACE/destdir/include && tar xvfp -)

if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
g++ srcdir/eltopo/eltopo3d/obj/*.o destdir/lib/libopenblas.a -o destdir/lib/eltopo.dll -fPIC -shared -lm -static-libgfortran -static-libstdc++ -fno-exceptions
else
g++ srcdir/eltopo/eltopo3d/obj/*.o destdir/lib/libopenblas.a -o destdir/lib/eltopo.so -fPIC -shared -lm -static-libgfortran -static-libstdc++
fi
"""
# -Wl,-Bstatic -lgfortran

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Windows(:x86_64),
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "eltopo", :eltopo)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BinaryBuilder.InlineBuildDependency(join(readlines("build_OpenBLAS.v0.3.0.jl"),"\n"))
]

# Build the tarballs, and possibly a `build.jl` as well.
# platforms = expand_gcc_versions(platforms)
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
mutable struct CRT
    data::Ptr{Cvoid}

    function CRT()
        c = new()
        c.data = Libc.malloc(1000)  # 552
        ccall((:mbedtls_x509_crt_init, libmbedx509), Cvoid, (Ptr{Cvoid},), c.data)
        finalizer(c->begin
            ccall((:mbedtls_x509_crt_free, libmbedx509), Cvoid, (Ptr{Cvoid},), c.data)
            Libc.free(c.data)
        end, c)
        c
    end
end

show(io::IO, crt::CRT) = print(io, crt_info(crt))

function crt_info(crt::CRT)
    buf = zeros(UInt8, 1000)
    ccall((:mbedtls_x509_crt_info, libmbedx509), Cint,
        (Ptr{Cvoid}, Csize_t, Cstring, Ptr{Cvoid}),
        buf, 1000, "", crt.data)
    GC.@preserve buf unsafe_string(pointer(buf))
end

function crt_parse!(chain, buf::String)
    ret = ccall((:mbedtls_x509_crt_parse, libmbedx509), Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Csize_t),
        chain.data, buf, sizeof(buf)+1)
    ret == 0 || mbed_err(ret)
    chain
end

crt_parse!(chain, buf::IOStream) = crt_parse!(chain, String(read(buf)))

crt_parse_file(path) = crt_parse(String(read(path)))

function crt_parse(buf)
    crt = CRT()
    crt_parse!(crt, String(buf))
    crt
end


text:
abstract type AbstractMonomer end
abstract type AbstractGenomicData end
abstract type AbstractAssimilation end
abstract type AbstractMetabolism end
abstract type AbstractSoil end

@units @description struct BaseMonomer <: AbstractMonomer
    D_S::Array{Float64,1}     | m^2/s   | "Aqueous diffusivity"
    N_C::Array{Float64,1}               | "Number of C atoms"
    N_N::Array{Float64,1}               | "Number of N atoms"
    y_DE::Array{Float64,1}    | mol/mol | "Assimilation yield"
end

@units @description mutable struct PAssimilation{T<:AbstractGenomicData, D<:AbstractMonomer, G<:AbstractSoil, H<:AbstractMetabolism} <: AbstractAssimilation
   Isolate::T                           | "Isolate defined by BaseGenome"
   Monomer::D                           | "Monomer defined by BaseMonomer"
   Soil::G                              | "Soil properties defined by BaseSoil"
   PMetabolismC::H                      | "C metabolism"
   ρ_p::Array{Float64,2}                | "Transporter density on cell surface"
   N_SB::Array{Float64,2}  |mol/mol     | "Number of uptake sites"
   K_D_0::Array{Float64,2} |mol/m^3     | "Reference affinity constant"
   K_D::Array{Float64,2}   |mol/m^3     | "Affinity constant"
   #
   function PAssimilation{T, D, G, H}(Isolate::T, Monomer::D, Soil::G, PMetabolismC::H) where {T<:AbstractGenomicData, D<:AbstractMonomer, G<:AbstractSoil, H<:AbstractMetabolism}
      ρ_p = transporter_density(genome_size_to_cell_volume(Isolate.L_DNA), Isolate.min_gentime, PMetabolismC.k_E, PMetabolismC.α, PMetabolismC.y_EM, PMetabolismC.y_EV, Isolate.transporter_distr)
      N_SB = monomer_uptake_sites(genome_size_to_cell_volume(Isolate.L_DNA), ρ_p)
      K_D_0 = specific_reference_affinity(ρ_p, genome_size_to_cell_volume(Isolate.L_DNA), Monomer.D_S)
      K_D = affinity_constant(ρ_p, genome_size_to_cell_volume(Isolate.L_DNA), Monomer.D_S, Soil.pct_sand, Soil.pct_clay, Soil.s_sat, Soil.N_cell)
      new(Isolate, Monomer, Soil, PMetabolismC, ρ_p, N_SB, K_D_0, K_D)
   end
end

PAssimilation(Isolate::T, Monomer::D, Soil::G, PMetabolismC::H) where {T<:AbstractGenomicData, D<:AbstractMonomer, G<:AbstractSoil, H<:AbstractMetabolism} = PAssimilation{T, D, G, H}(Isolate::T, Monomer::D, Soil::G, PMetabolismC::H)

function specific_reference_affinity(ρ_p::Array{Float64,2}, V_c::Array{Float64,1}, D_S::Array{Float64,1})
    V_S = 100.0                   # s
    N_A = 6.022e23                # Avogrado constant [1/mol]
    r_p = 1e-9                    # 1 μm
    r_c = (3*V_c/(4*pi)).^(1/3)
    #
    p_inv = zeros(size(ρ_p,1), size(ρ_p,2))
    for i in 1:size(ρ_p,2)
        p_inv[:,i] = @. (4*ρ_p[:,i]*r_c[i]^2/r_p + pi*r_c[i])/(4*ρ_p[:,i]*r_c[i]^2/r_p)
    end
    #
    K_SC_0 = zeros(size(ρ_p,1), size(ρ_p,2))
    for i in 1:size(ρ_p,2)
        K_SC_0[:,i] = @. (V_S*ρ_p[:,i]*r_c[i])/(pi*D_S*r_p^2*N_A)*p_inv[:,i]
    end
    #
    any(x->x==true, isnan.(K_SC_0)) ? throw(DomainError("NaN in DEBJulia.ECA_kinetics!"))  : return K_SC_0
end

function affinity_constant(ρ_p::Array{Float64,2}, V_c::Array{Float64,1}, D_S::Array{Float64,1}, pct_sand::Float64, pct_clay::Float64, s_sat::Float64, N_cell::Array{Float64,1})
    V_S = 100.0
    N_A = 6.022e23
    ϕ_w, τ_g, τ_w, δ  = soil_affinity_properties(pct_sand, pct_clay, s_sat)
    K_S_0 = specific_reference_affinity(ρ_p, V_c, D_S)
    r_c = (3*V_c/(4*pi)).^(1/3)
    r_p = 1e-9
    k1_s = zeros(size(ρ_p,1), size(ρ_p,2))
    for i in 1:size(ρ_p,2)
        k1_s[:,i] = @. V_S*ρ_p[:,i]*4*r_c[i]^2/(K_S_0[:,i]*r_p^2)
    end
    #k1_s  = @. V_S*ρ_p*4*r_c^2/(K_S_0*r_p^2)
    r_m = @. r_c*(80.0*N_cell)^(1/3)
    v_m = @. 4/3*pi*r_m^3
    #
    kapvsi = zeros(size(ρ_p,1), size(ρ_p,2))
    for i in 1:size(ρ_p,2)
        kapvsi[:,i] = @. (δ/(D_S*r_m[i]*(r_m[i]+δ)) + 1/(D_S*τ_w*ϕ_w*(r_m[i]+δ)+1.e-20))*v_m[i]/(4.0*pi)
    end
    #
    K = zeros(size(ρ_p,1), size(ρ_p,2))
    for i in 1:size(ρ_p, 2)
        K[:,i] =  @. K_S_0[:,i]*(1 + (kapvsi[:,i])*k1_s[:,i]*N_cell[i]/(N_A*v_m[i]))
    end
    any(x->x==true, isnan.(K)) ? throw(DomainError("NaN in DEBJulia.ECA_kinetics!"))  : return K
end

function monomer_uptake_sites(V_c::Array{Float64,1}, ρ_p::Array{Float64,2})
    n_0 = cell_volume_to_cell_number_density(V_c)
    r_p = 1e-9                    # 1 μm
    r_c = (3*V_c/(4*pi)).^(1/3)
    N_SB = zeros(size(ρ_p,1), size(ρ_p,2))
    for i in 1:size(V_c,1)
        N_SB[:,i] = n_0[i]*4*r_c[i]^2/(r_p^2)*ρ_p[:,i]
    end
    N_SB[N_SB.<1e-16].=0.0
    return N_SB
end

function transporter_density(V_c::Array{Float64,1}, N_p::Array{Float64,2})
    r_p = 1e-9                    # 1 μm
    r_c = (3*V_c/(4*pi)).^(1/3)
    ρ_p = zeros(size(N_p,1), size(N_p,2))
    for i in 1:size(V_c,1)
        ρ_p[:,i] = N_p[:,i]*r_p^2/(4*r_c[i]^2)
    end
    return ρ_p
end

function transporter_density(V_c::Array{Float64,1}, min_gentime::Array{Float64,1}, k_E::Array{Float64,1}, α::Array{Float64,1}, y_EM::Array{Float64,1}, y_EV::Array{Float64,1}, genome_distr::Array{Float64,2})
    N_C = 10
    V_S = 100.0
    n_0 = cell_volume_to_cell_number_density(V_c)
    min_gentime = min_gentime*60*60   # in seconds
    gmax = log(2)./(min_gentime) # max. specific growth rate [1/s]
    k_M = cell_volume_to_specific_maintenance_rate(V_c)
    N_p = @. (gmax*(1+α)*y_EV + y_EM*k_M)/(N_C*V_S*n_0 - gmax*V_S*n_0*N_C/k_E)
    #
    closure   = zeros(size(genome_distr))
    for i in 1:size(closure,2)
        closure[:,i] = N_p[i]*genome_distr[:,i]./max(1e-8,sum(genome_distr[:,i]))
    end
    ceil_closure = ceil.(closure)
    ceil_closure[ceil_closure.==0.0].=1e-8
    #
    r_p = 1e-9                    # 1 μm
    r_c = (3*V_c/(4*pi)).^(1/3)
    ρ_p = zeros(size(ceil_closure,1), size(ceil_closure,2))
    for i in 1:size(V_c,1)
        ρ_p[:,i] = ceil_closure[:,i]*r_p^2/(4*r_c[i]^2)
    end
    return ρ_p
end


text:
"""
    unsafe_createfielddefn(name::AbstractString, etype::OGRFieldType)

Create a new field definition.

By default, fields have no width, precision, are nullable and not ignored.
"""
unsafe_createfielddefn(name::AbstractString, etype::OGRFieldType)::FieldDefn =
    FieldDefn(GDAL.ogr_fld_create(name, etype))

"Destroy a field definition."
function destroy(fielddefn::FieldDefn)::Nothing
    GDAL.ogr_fld_destroy(fielddefn.ptr)
    fielddefn.ptr = C_NULL
    return nothing
end

function destroy(fielddefn::IFieldDefnView)::Nothing
    fielddefn.ptr = C_NULL
    return nothing
end

"Set the name of this field."
function setname!(fielddefn::FieldDefn, name::AbstractString)::FieldDefn
    GDAL.ogr_fld_setname(fielddefn.ptr, name)
    return fielddefn
end

"Fetch the name of this field."
getname(fielddefn::AbstractFieldDefn)::String =
    GDAL.ogr_fld_getnameref(fielddefn.ptr)

"Fetch the type of this field."
gettype(fielddefn::AbstractFieldDefn)::OGRFieldType =
    GDAL.ogr_fld_gettype(fielddefn.ptr)

"Set the type of this field."
function settype!(fielddefn::FieldDefn, etype::OGRFieldType)::FieldDefn
    GDAL.ogr_fld_settype(fielddefn.ptr, etype)
    return fielddefn
end

"""
    getsubtype(fielddefn::AbstractFieldDefn)

Fetch subtype of this field.

### Parameters
* `fielddefn`: handle to the field definition to get subtype from.

### Returns
field subtype.
"""
getsubtype(fielddefn::AbstractFieldDefn)::OGRFieldSubType =
    GDAL.ogr_fld_getsubtype(fielddefn.ptr)

"""
    setsubtype!(fielddefn::FieldDefn, subtype::OGRFieldSubType)

Set the subtype of this field.

This should never be done to an OGRFieldDefn that is already part of an
OGRFeatureDefn.

### Parameters
* `fielddefn`: handle to the field definition to set type to.
* `subtype`: the new field subtype.
"""
function setsubtype!(fielddefn::FieldDefn, subtype::OGRFieldSubType)::FieldDefn
    GDAL.ogr_fld_setsubtype(fielddefn.ptr, subtype)
    return fielddefn
end

"""
    getjustify(fielddefn::AbstractFieldDefn)

Get the justification for this field.

Note: no driver is know to use the concept of field justification.
"""
getjustify(fielddefn::AbstractFieldDefn)::OGRJustification =
    GDAL.ogr_fld_getjustify(fielddefn.ptr)

"""
    setjustify!(fielddefn::FieldDefn, ejustify::OGRJustification)

Set the justification for this field.

Note: no driver is know to use the concept of field justification.
"""
function setjustify!(
    fielddefn::FieldDefn,
    ejustify::OGRJustification,
)::FieldDefn
    GDAL.ogr_fld_setjustify(fielddefn.ptr, ejustify)
    return fielddefn
end

"""
    getwidth(fielddefn::AbstractFieldDefn)

Get the formatting width for this field.

### Returns
the width, zero means no specified width.
"""
getwidth(fielddefn::AbstractFieldDefn)::Integer =
    GDAL.ogr_fld_getwidth(fielddefn.ptr)

"""
    setwidth!(fielddefn::FieldDefn, width::Integer)

Set the formatting width for this field in characters.

This should never be done to an OGRFieldDefn that is already part of an
OGRFeatureDefn.
"""
function setwidth!(fielddefn::FieldDefn, width::Integer)::FieldDefn
    GDAL.ogr_fld_setwidth(fielddefn.ptr, width)
    return fielddefn
end

"""
    getprecision(fielddefn::AbstractFieldDefn)

Get the formatting precision for this field.

This should normally be zero for fields of types other than OFTReal.
"""
getprecision(fielddefn::AbstractFieldDefn)::Integer =
    GDAL.ogr_fld_getprecision(fielddefn.ptr)

"""
    setprecision!(fielddefn::FieldDefn, precision::Integer)

Set the formatting precision for this field in characters.

This should normally be zero for fields of types other than OFTReal.
"""
function setprecision!(fielddefn::FieldDefn, precision::Integer)::FieldDefn
    GDAL.ogr_fld_setprecision(fielddefn.ptr, precision)
    return fielddefn
end

"""
    setparams!(fielddefn, name, etype, [nwidth, [nprecision, [justify]]])

Set defining parameters for a field in one call.

### Parameters
* `fielddefn`:  the field definition to set to.
* `name`:       the new name to assign.
* `etype`:      the new type (one of the OFT values like OFTInteger).
* `nwidth`:     the preferred formatting width. 0 (default) indicates undefined.
* `nprecision`: number of decimals for formatting. 0 (default) for undefined.
* `justify`:    the formatting justification ([OJUndefined], OJLeft or OJRight)
"""
function setparams!(
    fielddefn::FieldDefn,
    name::AbstractString,
    etype::OGRFieldType;
    nwidth::Integer = 0,
    nprecision::Integer = 0,
    justify::OGRJustification = OJUndefined,
)::FieldDefn
    GDAL.ogr_fld_set(fielddefn.ptr, name, etype, nwidth, nprecision, justify)
    return fielddefn
end

"""
    isignored(fielddefn::AbstractFieldDefn)

Return whether this field should be omitted when fetching features.
"""
isignored(fielddefn::AbstractFieldDefn)::Bool =
    Bool(GDAL.ogr_fld_isignored(fielddefn.ptr))

"""
    setignored!(fielddefn::FieldDefn, ignore::Bool)

Set whether this field should be omitted when fetching features.
"""
function setignored!(fielddefn::FieldDefn, ignore::Bool)::FieldDefn
    GDAL.ogr_fld_setignored(fielddefn.ptr, ignore)
    return fielddefn
end

"""
    isnullable(fielddefn::AbstractFieldDefn)

Return whether this field can receive null values.

By default, fields are nullable.

Even if this method returns `false` (i.e not-nullable field), it doesn't mean
that OGRFeature::IsFieldSet() will necessarily return `true`, as fields can be
temporarily unset and null/not-null validation is usually done when
OGRLayer::CreateFeature()/SetFeature() is called.
"""
isnullable(fielddefn::AbstractFieldDefn)::Bool =
    Bool(GDAL.ogr_fld_isnullable(fielddefn.ptr))

"""
    setnullable!(fielddefn::FieldDefn, nullable::Bool)

Set whether this field can receive null values.

By default, fields are nullable, so this method is generally called with `false`
to set a not-null constraint.

Drivers that support writing not-null constraint will advertize the
GDAL_DCAP_NOTNULL_FIELDS driver metadata item.
"""
function setnullable!(fielddefn::FieldDefn, nullable::Bool)::FieldDefn
    GDAL.ogr_fld_setnullable(fielddefn.ptr, nullable)
    return fielddefn
end

"""
    getdefault(fielddefn::AbstractFieldDefn)

Get default field value
"""
function getdefault(fielddefn::AbstractFieldDefn)::String
    result =
        @gdal(OGR_Fld_GetDefault::Cstring, fielddefn.ptr::GDAL.OGRFieldDefnH)
    return if result == C_NULL
        ""
    else
        unsafe_string(result)
    end
end

"""
    setdefault!(fielddefn::AbstractFieldDefn, default)

Set default field value.

The default field value is taken into account by drivers (generally those with
a SQL interface) that support it at field creation time. OGR will generally not
automatically set the default field value to null fields by itself when calling
OGRFeature::CreateFeature() / OGRFeature::SetFeature(), but will let the
low-level layers to do the job. So retrieving the feature from the layer is
recommended.

The accepted values are NULL, a numeric value, a literal value enclosed between
single quote characters (and inner single quote characters escaped by repetition
of the single quote character), CURRENT_TIMESTAMP, CURRENT_TIME, CURRENT_DATE or
a driver specific expression (that might be ignored by other drivers). For a
datetime literal value, format should be 'YYYY/MM/DD HH:MM:SS[.sss]'
(considered as UTC time).

Drivers that support writing DEFAULT clauses will advertize the
GDAL_DCAP_DEFAULT_FIELDS driver metadata item.
"""
function setdefault!(fielddefn::T, default)::T where {T<:AbstractFieldDefn}
    GDAL.ogr_fld_setdefault(fielddefn.ptr, default)
    return fielddefn
end

"""
    isdefaultdriverspecific(fielddefn::AbstractFieldDefn)

Returns whether the default value is driver specific.

Driver specific default values are those that are not NULL, a numeric value, a
literal value enclosed between single quote characters, CURRENT_TIMESTAMP,
CURRENT_TIME, CURRENT_DATE or datetime literal value.
"""
isdefaultdriverspecific(fielddefn::AbstractFieldDefn)::Bool =
    Bool(GDAL.ogr_fld_isdefaultdriverspecific(fielddefn.ptr))

"""
    unsafe_creategeomdefn(name::AbstractString, etype::OGRwkbGeometryType)

Create a new field geometry definition.
"""
function unsafe_creategeomdefn(
    name::AbstractString,
    etype::OGRwkbGeometryType,
)::GeomFieldDefn
    return GeomFieldDefn(GDAL.ogr_gfld_create(name, etype))
end

"Destroy a geometry field definition."
function destroy(geomdefn::GeomFieldDefn)::Nothing
    GDAL.ogr_gfld_destroy(geomdefn.ptr)
    geomdefn.ptr = C_NULL
    geomdefn.spatialref = SpatialRef()
    return nothing
end

"Destroy a geometry field definition."
function destroy(geomdefn::IGeomFieldDefnView)::Nothing
    geomdefn.ptr = C_NULL
    return nothing
end

"Set the name of this field."
function setname!(geomdefn::GeomFieldDefn, name::AbstractString)::GeomFieldDefn
    GDAL.ogr_gfld_setname(geomdefn.ptr, name)
    return geomdefn
end

"Fetch name of this field."
getname(geomdefn::AbstractGeomFieldDefn)::String =
    GDAL.ogr_gfld_getnameref(geomdefn.ptr)

"Fetch geometry type of this field."
gettype(geomdefn::AbstractGeomFieldDefn)::OGRwkbGeometryType =
    GDAL.ogr_gfld_gettype(geomdefn.ptr)

"Set the geometry type of this field."
function settype!(
    geomdefn::GeomFieldDefn,
    etype::OGRwkbGeometryType,
)::GeomFieldDefn
    GDAL.ogr_gfld_settype(geomdefn.ptr, etype)
    return geomdefn
end

"""
    getspatialref(geomdefn::AbstractGeomFieldDefn)

Returns a clone of the spatial reference system for this field. May be NULL.
"""
function getspatialref(geomdefn::AbstractGeomFieldDefn)::ISpatialRef
    result = GDAL.ogr_gfld_getspatialref(geomdefn.ptr)
    if result == C_NULL
        return ISpatialRef()
    else
        # NOTE(yeesian): we make a clone here so that the spatialref does not
        # depend on the GeomFieldDefn/Dataset.
        return ISpatialRef(GDAL.osrclone(result))
    end
end

function unsafe_getspatialref(geomdefn::AbstractGeomFieldDefn)::SpatialRef
    result = GDAL.ogr_gfld_getspatialref(geomdefn.ptr)
    if result == C_NULL
        return SpatialRef()
    else
        # NOTE(yeesian): we make a clone here so that the spatialref does not
        # depend on the GeomFieldDefn/Dataset.
        return SpatialRef(GDAL.osrclone(result))
    end
end

"""
    setspatialref!(geomdefn::GeomFieldDefn, spatialref::AbstractSpatialRef)

Set the spatial reference of this field.

This function drops the reference of the previously set SRS object and acquires
a new reference on the passed object (if non-NULL).
"""
function setspatialref!(
    geomdefn::GeomFieldDefn,
    spatialref::AbstractSpatialRef,
)::GeomFieldDefn
    clonespatialref = clone(spatialref)
    GDAL.ogr_gfld_setspatialref(geomdefn.ptr, clonespatialref.ptr)
    geomdefn.spatialref = clonespatialref
    return geomdefn
end

"""
    isnullable(geomdefn::AbstractGeomFieldDefn)

Return whether this geometry field can receive null values.

By default, fields are nullable.

Even if this method returns `false` (i.e not-nullable field), it doesn't mean
that OGRFeature::IsFieldSet() will necessary return `true`, as fields can be
temporarily unset and null/not-null validation is usually done when
OGRLayer::CreateFeature()/SetFeature() is called.

Note that not-nullable geometry fields might also contain 'empty' geometries.
"""
isnullable(geomdefn::AbstractGeomFieldDefn)::Bool =
    Bool(GDAL.ogr_gfld_isnullable(geomdefn.ptr))

"""
    setnullable!(geomdefn::GeomFieldDefn, nullable::Bool)

Set whether this geometry field can receive null values.

By default, fields are nullable, so this method is generally called with `false`
to set a not-null constraint.

Drivers that support writing not-null constraint will advertize the
GDAL_DCAP_NOTNULL_GEOMFIELDS driver metadata item.
"""
function setnullable!(geomdefn::GeomFieldDefn, nullable::Bool)::GeomFieldDefn
    GDAL.ogr_gfld_setnullable(geomdefn.ptr, nullable)
    return geomdefn
end

"""
    isignored(geomdefn::AbstractGeomFieldDefn)

Return whether this field should be omitted when fetching features.
"""
isignored(geomdefn::AbstractGeomFieldDefn)::Bool =
    Bool(GDAL.ogr_gfld_isignored(geomdefn.ptr))

"""
    setignored!(geomdefn::GeomFieldDefn, ignore::Bool)

Set whether this field should be omitted when fetching features.
"""
function setignored!(geomdefn::GeomFieldDefn, ignore::Bool)::GeomFieldDefn
    GDAL.ogr_gfld_setignored(geomdefn.ptr, ignore)
    return geomdefn
end


text:
module ConnectionRequest

import ..Layer, ..request
using ..URIs
using ..Messages
using ..IOExtras
using ..ConnectionPool
using MbedTLS.SSLContext
import ..@debug, ..DEBUG_LEVEL

"""
    request(ConnectionPoolLayer, ::URI, ::Request, body) -> HTTP.Response

Retrieve an `IO` connection from the [`ConnectionPool`](@ref).

Close the connection if the request throws an exception.
Otherwise leave it open so that it can be reused.

`IO` related exceptions from `Base` are wrapped in `HTTP.IOError`.
See [`isioerror`](@ref).
"""

abstract type ConnectionPoolLayer{Next <: Layer} <: Layer end
export ConnectionPoolLayer

function request(::Type{ConnectionPoolLayer{Next}}, url::URI, req, body;
                 socket_type::Type=TCPSocket, kw...) where Next

    IOType = ConnectionPool.Transaction{sockettype(url, socket_type)}
    io = getconnection(IOType, url.host, url.port; kw...)

    try
        return request(Next, io, req, body; kw...)
    catch e
        @debug 1 "❗  ConnectionLayer $e. Closing: $io"
        close(io)
        rethrow(isioerror(e) ? IOError(e, "during request($url)") : e)
    end
end

sockettype(url::URI, default) = url.scheme in ("wss", "https") ? SSLContext :
                                                                 default

end # module ConnectionRequest