text:
struct Stockholm <: FileFormat end

@inline function _fill_with_sequence_line!(IDS, SEQS, line)
    if !startswith(line,'#') && !startswith(line,"//")
        words = get_n_words(line, 2)
        @inbounds id = words[1]
        if id in IDS
            # It's useful when sequences are split into several lines
            # It can be a problem with duplicated IDs
            i = something(findfirst(isequal(id), IDS), 0)
            SEQS[i] = SEQS[i] * words[2]
        else
            push!(IDS, id)
            push!(SEQS, words[2])
        end
    end
end

function _fill_with_line!(IDS, SEQS, GF, GS, GC, GR, line)
    if startswith(line,"#=GF")
        words = get_n_words(line, 3)
        id = words[2]
        if id in keys(GF)
            GF[ id ] = GF[ id ] * "\n" * words[3]
        else
            GF[ id ] = words[3]
        end
    elseif startswith(line,"#=GS")
        words = get_n_words(line, 4)
        idtuple = (words[2],words[3])
        if idtuple in keys(GS)
            GS[ idtuple ] = GS[ idtuple ] * "\n" * words[4]
        else
            GS[ idtuple ] = words[4]
        end
    elseif startswith(line,"#=GC")
        words = get_n_words(line, 3)
        GC[words[2]] = words[3]
    elseif startswith(line,"#=GR")
        words = get_n_words(line, 4)
        GR[(words[2],words[3])] = words[4]
    else
        _fill_with_sequence_line!(IDS, SEQS, line)
    end
end

function _pre_readstockholm(io::Union{IO, AbstractString})
    IDS  = OrderedSet{String}()
    SEQS = String[]
    GF = OrderedDict{String,String}()
    GC = Dict{String,String}()
    GS = Dict{Tuple{String,String},String}()
    GR = Dict{Tuple{String,String},String}()

    @inbounds for line::String in lineiterator(io)
        _fill_with_line!(IDS, SEQS, GF, GS, GC, GR, line)
        if startswith(line, "//")
           break
        end
    end

    _check_seq_len(IDS, SEQS)

    GF = sizehint!(GF, length(GF))
    GC = sizehint!(GC, length(GC))
    GS = sizehint!(GS, length(GS))
    GR = sizehint!(GR, length(GR))
    (IDS, SEQS, GF, GS, GC, GR)
end

function _pre_readstockholm_sequences(io::Union{IO, AbstractString})
    IDS  = OrderedSet{String}()
    SEQS = String[]
    @inbounds for line::String in lineiterator(io)
        _fill_with_sequence_line!(IDS, SEQS, line)
        if startswith(line, "//")
           break
        end
    end
    _check_seq_len(IDS, SEQS)
    (IDS, SEQS)
end

function Base.parse(io::Union{IO, AbstractString},
                   format::Type{Stockholm},
                   output::Type{AnnotatedMultipleSequenceAlignment};
                   generatemapping::Bool=false,
                   useidcoordinates::Bool=false,
                   deletefullgaps::Bool=true,
                   keepinserts::Bool=false)::AnnotatedMultipleSequenceAlignment
    IDS, SEQS, GF, GS, GC, GR = _pre_readstockholm(io)
    annot = Annotations(GF, GS, GC, GR)
    _generate_annotated_msa(annot, collect(IDS), SEQS, keepinserts,
                            generatemapping, useidcoordinates, deletefullgaps)
end

function Base.parse(io::Union{IO, AbstractString},
                   format::Type{Stockholm},
                   output::Type{NamedResidueMatrix{Array{Residue,2}}};
                   deletefullgaps::Bool=true)::NamedResidueMatrix{Array{Residue,2}}
    IDS, SEQS = _pre_readstockholm_sequences(io)
    msa = _generate_named_array(SEQS, IDS)
    if deletefullgaps
        return deletefullgapcolumns(msa)
    end
    msa
end

function Base.parse(io::Union{IO, AbstractString},
                   format::Type{Stockholm},
                   output::Type{MultipleSequenceAlignment};
                   deletefullgaps::Bool=true)::MultipleSequenceAlignment
    msa = parse(io,
                format,
                NamedResidueMatrix{Array{Residue,2}},
                deletefullgaps=deletefullgaps)
    MultipleSequenceAlignment(msa)
end

function Base.parse(io::Union{IO,AbstractString},
                   format::Type{Stockholm},
                   output::Type{Matrix{Residue}};
                   deletefullgaps::Bool=true)::Matrix{Residue}
    IDS, SEQS = _pre_readstockholm_sequences(io)
    _strings_to_matrix_residue_unsafe(SEQS, deletefullgaps)
end

function Base.parse(io, format::Type{Stockholm};
                    generatemapping::Bool=false,
                    useidcoordinates::Bool=false,
                    deletefullgaps::Bool=true,
                    keepinserts::Bool=false)::AnnotatedMultipleSequenceAlignment
    parse(io, Stockholm, AnnotatedMultipleSequenceAlignment,
          generatemapping=generatemapping,
          useidcoordinates=useidcoordinates,
          deletefullgaps=deletefullgaps,
          keepinserts=keepinserts)
end

# Print Pfam
# ==========

function _to_sequence_dict(annotation::Dict{Tuple{String,String},String})
    seq_dict = Dict{String,Vector{String}}()
    for (key, value) in annotation
        seq_id = key[1]
        if haskey(seq_dict, seq_id)
            push!(seq_dict[seq_id], string(seq_id, '\t', key[2], '\t', value))
        else
            seq_dict[seq_id] = [ string(seq_id, '\t', key[2], '\t', value) ]
        end
    end
    sizehint!(seq_dict, length(seq_dict))
end

function Base.print(io::IO, msa::AnnotatedMultipleSequenceAlignment,
                    format::Type{Stockholm})
    _printfileannotations(io, msa.annotations)
    _printsequencesannotations(io, msa.annotations)
    res_annotations = _to_sequence_dict(msa.annotations.residues)
    seqnames = sequencenames(msa)
    for i in 1:nsequences(msa)
        id = seqnames[i]
        seq = stringsequence(msa, i)
        println(io, id, "\t\t\t", seq)
        if id in keys(res_annotations)
            for line in res_annotations[id]
                println(io, "#=GR ", line)
            end
        end
    end
    _printcolumnsannotations(io, msa.annotations)
    println(io, "//")
end

function Base.print(io::IO, msa::AbstractMatrix{Residue}, format::Type{Stockholm})
    seqnames = sequencenames(msa)
    for i in 1:nsequences(msa)
        println(io, seqnames[i], "\t\t\t", stringsequence(msa, i))
    end
    println(io, "//")
end

Base.print(msa::AnnotatedMultipleSequenceAlignment) = print(stdout, msa, Stockholm)


text:
using Base: @deprecate

@deprecate ci(args...) confint(args...)

ext: jl


text:
function interactiveWorkaround(window::Gtk.GtkWindow)
    # If used outside of REPL this keeps the program running until the window is closed
    if !isinteractive()
        c = Condition()
        signal_connect(window, :destroy) do widget
            notify(c)
        end
        wait(c)
    end
end


text:
# Oil production and petro states
export total_oil_prod, oil_gdp_of_nordics, per_capita_pumped_oil, oil_percent_gdp, oil_rent
using CairoMakie, ElectronDisplay

noto_sans = assetpath("fonts", "NotoSans-Regular.ttf")
noto_sans_bold = assetpath("fonts", "NotoSans-Bold.ttf")

function oil_gdp_of_nordics()
    petrostates_path = joinpath(datadir, "petro-states-data.csv")
    petrostates = CSV.File(petrostates_path) |> DataFrame
    petrostates.OilGDP = petrostates.GDP .* (petrostates.OilPercent ./ 100)

    nordics = ["Norway", "Denmark", "Sweden", "Finland", "Iceland"]
    nordics = DataFrame(Country=nordics)
    nordics = innerjoin(petrostates, nordics, on=:Country)

    (sum(nordics.OilGDP) / sum(nordics.GDP)) * 100
end

function oil_percent_gdp()
    petrostates_path = joinpath(datadir, "petro-states-data.csv")
    petrostates = CSV.File(petrostates_path) |> DataFrame

    countries = ["Norway", "Nordics", "United States", "Canada", "Russia", "Saudi Arabia"]

    filter!(petrostates) do row
       row.Country in countries
    end

    fig = Figure(resolution = (1024, 600), font = noto_sans)
    xs = 1:nrow(petrostates)
    ys = petrostates.OilPercent

    countries = replace(countries, "United States" => "USA")
    ax = Axis(fig[1, 1],
        xticks = (xs, countries),
        ylabel="Percent",
        title="GDP based on Oil")
    barplot!(ax, xs, ys, color=xs)

    fig
end

function oil_rent()
    petrostates_path = joinpath(datadir, "petro-states-data.csv")
    petrostates = CSV.File(petrostates_path) |> DataFrame

    countries = ["Norway", "Kuwait", "Iraq", "UAE", "Saudi Arabia"]

    filter!(petrostates) do row
       row.Country in countries
    end

    fig = Figure(resolution = (1024, 600), font = noto_sans)
    xs = 1:nrow(petrostates)
    ys = petrostates.OilRent

    ax = Axis(fig[1, 1],
        xticks = (xs, countries),
        ylabel="Percent",
        title="Oil Rent")
    barplot!(ax, xs, ys, color=xs)

    fig
end

# Plot total oil production for various petro states from 1900 to present time
function total_oil_prod()
    oilprod_path = joinpath(datadir, "oil-production-by-country.csv")
    oilprod = CSV.File(oilprod_path) |> DataFrame

    petrostates_path = joinpath(datadir, "petro-states-data.csv")

    select!(oilprod, :Entity => :Country, :Year, Symbol("Oil production (TWh)") => :Oil)
    petrostates = CSV.File(petrostates_path) |> DataFrame
    sample = select(petrostates, :Country)
    oilprod = innerjoin(oilprod, sample, on=:Country)

    filter!(oilprod) do row
        # row.Year > 2015 &&
        row.Country ∉ ["Finland", "Sweden", "Denmark", "Iceland", "Nordic"]
    end

    oilgroups = groupby(oilprod, :Country)
    totaloilprod = combine(oilgroups, :Oil => sum => :Oil)

    fig = Figure(resolution = (1024, 600), font = noto_sans)
    xs = 1:length(totaloilprod.Country)
    ys = totaloilprod.Oil
    countries = replace(totaloilprod.Country, "United States" => "USA")
    ax = Axis(fig[1, 1], xticks = (xs, countries), ylabel="Oil Prod (TWh)")
    barplot!(ax, xs, ys, color=xs)

    fig
end

function per_capita_pumped_oil()
    oilprod_path = joinpath(datadir, "oil-prod-per-capita.csv")
    oilprod = CSV.File(oilprod_path) |> DataFrame

    petrostates_path = joinpath(datadir, "petro-states-data.csv")

    select!(oilprod, :Entity => :Country, :Year, Symbol("Oil production per capita (kWh)") => :Oil)
    petrostates = CSV.File(petrostates_path) |> DataFrame
    sample = select(petrostates, :Country)
    oilprod = innerjoin(oilprod, sample, on=:Country)

    endyear = 1989

    filter!(oilprod) do row
        row.Year <= endyear &&
        row.Country ∉ ["Finland", "Sweden", "Denmark", "Iceland", "Kuwait", "Saudi Arabia", "Nordic"]
    end

    oilgroups = groupby(oilprod, :Country)
    totaloilprod = combine(oilgroups, :Oil => sum => :Oil)

    fig = Figure(resolution = (1024, 600), font = noto_sans)
    xs = 1:length(totaloilprod.Country)

    # Get average per capita oil production from 1900 to 1990
    ys = totaloilprod.Oil ./ (endyear - 1900)

    countries = replace(totaloilprod.Country, "United States" => "USA")
    ax = Axis(fig[1, 1],
        xticks = (xs, countries),
        ylabel="per capita oil production (kWh)",
        title = "Average oil production per capita (1900 to $endyear)")
    barplot!(ax, xs, ys, color=xs)

    fig
end


text:
# the name of a distribution
#
#   Generally, this should be just the type name, e.g. Normal.
#   Under certain circumstances, one may want to specialize
#   this function to provide a name that is easier to read,
#   especially when the type is parametric.
#
distrname(d::Distribution) = string(typeof(d))

show(io::IO, d::Distribution) = show(io, d, typeof(d).names)

# For some distributions, the fields may contain internal details,
# which we don't want to show, this function allows one to
# specify which fields to show.
#
function show(io::IO, d::Distribution, pnames::(Symbol...))
    uml, namevals = _use_multline_show(d, pnames)
    uml ? show_multline(io, d, namevals) : show_oneline(io, d, namevals)
end

function _use_multline_show(d::Distribution, pnames::(Symbol...))
    # decide whether to use one-line or multi-line format
    #
    # Criteria: if total number of values is greater than 8, or
    # there are matrix-valued params, we use multi-line format
    #
    namevals = (Symbol, Any)[]
    multline = false
    tlen = 0
    for (i, p) in enumerate(pnames)
        pv = d.(p)
        if !(isa(pv, Number) || isa(pv, (Number...)) || isa(pv, AbstractVector))
            multline = true
        end
        tlen += length(pv)
        push!(namevals, (p, pv))
    end
    if tlen > 8
        multline = true
    end
    return (multline, namevals)
end
_use_multline_show(d::Distribution) = _use_multline_show(d, typeof(d).names)

function show_oneline(io::IO, d::Distribution, namevals)
    print(io, distrname(d))
    np = length(namevals)
    print(io, '(')
    for (i, nv) in enumerate(namevals)
        (p, pv) = nv
        print(io, p)
        print(io, '=')
        show(io, pv)
        if i < np
            print(io, ", ")
        end
    end
    print(io, ')')
end

function show_multline(io::IO, d::Distribution, namevals)
    print(io, distrname(d))
    println(io, "(")
    for (p, pv) in namevals
        print(io, p)
        print(io, ": ")
        println(io, pv)
    end
    println(io, ")")
end


text:
mutable struct FilledExtrapolation{T,N,ITP<:AbstractInterpolation,IT,FT} <: AbstractExtrapolation{T,N,ITP,IT}
    itp::ITP
    fillvalue::FT
end

function FilledExtrapolation(itp::AbstractInterpolation{T,N,IT}, fillvalue) where {T,N,IT}
    Te = promote_type(T,typeof(fillvalue))
    FilledExtrapolation{Te,N,typeof(itp),IT,typeof(fillvalue)}(itp, fillvalue)
end

@noinline FilledExtrapolation(itp::AbstractInterpolation{T,N,IT}, ::Type{F}) where {T,N,IT,F} =
    throw(ArgumentError("cannot create a filled extrapolation with a type $F, use a value of this type instead (e.g., $F(0))"))
@noinline FilledExtrapolation(itp::AbstractInterpolation{T,N,IT}, ::Type{F}) where {T,N,IT,F<:BoundaryCondition} =
    throw(ArgumentError("cannot create a filled extrapolation with a type $F, use a value of this type instead (e.g., $F())"))

Base.parent(A::FilledExtrapolation) = A.itp
etpflag(A::FilledExtrapolation) = A.fillvalue
itpflag(A::FilledExtrapolation) = itpflag(A.itp)

"""
`extrapolate(itp, fillvalue)` creates an extrapolation object that returns the `fillvalue` any time the indexes in `itp(x1,x2,...)` are out-of-bounds.
"""
extrapolate(itp::AbstractInterpolation{T,N,IT}, fillvalue) where {T,N,IT} = FilledExtrapolation(itp, fillvalue)

@inline function (etp::FilledExtrapolation{T,N})(x::Vararg{Number,N}) where {T,N}
    itp = parent(etp)
    coefs = coefficients(itp)
    Tret = typeof(prod(x) * zero(eltype(coefs)))
    if checkbounds(Bool, itp, x...)
        @inbounds itp(x...)
    else
        convert(Tret, etp.fillvalue)
    end
end
@inline function (etp::FilledExtrapolation{T,N})(args::Vararg{Number,M}) where {T,M,N}
    inds, trailing = Base.IteratorsMD.split(args, Val(N))
    @boundscheck all(x->x==1, trailing) || Base.throw_boundserror(etp, args)
    @assert length(inds) == N
    etp(inds...)
end
@inline function (etp::FilledExtrapolation{T,N})(x::Vararg{Union{Number,AbstractVector},N}) where {T,N}
    itp = parent(etp)
    Tret = typeof(lispyprod(zero(T), x...))
    ret = fill(convert(Tret, etp.fillvalue), shape(x...))
    axsib = inbounds(itp, x...)
    any(isempty, axsib) && return ret
    xib = getindex.(x, axsib)
    getindex!(view(ret, keepvectors(axsib...)...), itp, xib...)
    return ret
end

expand_index_resid_etp(deg, fillvalue, (l, u), x, etp::FilledExtrapolation, xN) =
    (l <= x <= u || Base.throw_boundserror(etp, xN))

# expand_etp_valueE(fv::FT, etp::FilledExtrapolation{T,N,ITP,IT,FT}, x) where {T,N,ITP,IT,FT} = fv
# expand_etp_gradientE(fv::FT, etp::FilledExtrapolation{T,N,ITP,IT,FT}, x) where {T,N,ITP,IT,FT} =
#     zero(SVector{N,FT})
# expand_etp_hessianE(fv::FT, etp::FilledExtrapolation{T,N,ITP,IT,FT}, x) where {T,N,ITP,IT,FT} =
#     zero(Matrix{N,N,FT})


text:
import ..Resampling: resample
import ..UncertainValues: UncertainValue
import ..UncertainDatasets: UncertainValueDataset

include("add_uncertainvaluedataset.jl")
include("subtract_uncertainvaluedataset.jl")
include("multiply_uncertainvaluedataset.jl")
include("divide_uncertainvaluedataset.jl")
include("exponentiate_uncertainvaluedataset.jl")


text:
using SymPy.SpecialFuncs
using SymPy: Sym, sqrt, conjugate, symbols, PI, simplify,
             expand_func, rewrite, N, gamma
if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

@testset "SpecialFuns" begin
    a, b, x, y = symbols("a, b, x, y")
    n, m, θ, ϕ = symbols("n, m, theta, phi")
    ν = symbols("nu", integer=true)

    @test fresnels(Sym(0)) == 0
    @test fresnels(Sym(oo)) == Sym(1)/2
    @test diff(fresnels(x), x) == sin(PI*x^2/2)
    #@test evalf(fresnels(Sym(2)), 30) == Sym(parse(BigFloat, "0.343415678363698242195300815958"))

    @test  diff(fresnelc(x), x) == cos(PI*x^2/2)

    @test Ei(Sym(-1)) == Ei(exp(1im*PI))
    @test diff(Ei(x), x) == exp(x)/x
    @test diff(Si(x), x) == sin(x)/x
    @test diff(Ci(x), x) == cos(x)/x

    @test airyai(Sym(0)) == 3^(Sym(1)/3)/(3*gamma(Sym(2)/3))
    @test airybi(Sym(0)) == 3^(Sym(5)/6)/(3*gamma(Sym(2)/3))

    @test jacobi(0, a, b, x) == 1
    @test jacobi(1, a, b, x) == a/2 - b/2 + x*(a/2 + b/2 + 1)
    @test jacobi(n, 0, 0, x) == legendre(n, x)
    @test jacobi(n, a, b, -x) == (-1)^n*jacobi(n, b, a, x)
    @test conjugate(jacobi(n, a, b, x)) == jacobi(n, conjugate(a), conjugate(b), conjugate(x))

    @test gegenbauer(0, a, x) == 1
    @test gegenbauer(1, a, x) == 2a*x
    @test gegenbauer(2, a, x) == -a + x^2*(2a^2 + 2a)
    @test conjugate(gegenbauer(n, a, x)) == gegenbauer(n, conjugate(a), conjugate(x))
    @test diff(gegenbauer(n, a, x), x) == 2a*gegenbauer(n-1, a+1, x)

    @test chebyshevt(0, x) == 1
    @test chebyshevt(1, x) == x
    @test chebyshevt(2, x) == 2x^2-1

    @test chebyshevu(0, x) == 1
    @test chebyshevu(1, x) == 2x
    @test chebyshevu(2, x) == 4x^2 - 1
    @test chebyshevu(n, 0) == cos(π*n/2)
    @test chebyshevu(n, 1) == n + 1

    @test legendre(0, x) == 1
    @test legendre(1, x) == x
    @test legendre(2, x) == 3x^2/2 - 1/2
    @test legendre(n, x) == legendre(n, x)
    @test diff(legendre(n,x), x) == n*(x*legendre(n, x) - legendre(n - 1, x))/(x^2 - 1)

    @test assoc_legendre(0,0, x) == 1
    @test assoc_legendre(1,0, x) == x
    @test assoc_legendre(1,1, x) == -sqrt(-x^2 + 1)
    #@test assoc_legendre(n,m,x) == assoc_legendre(n, m, x)

    @test hermite(0, x) == 1
    @test hermite(1, x) == 2x
    @test hermite(2, x) == 4x^2 - 2
    @test hermite(n, x) == hermite(n, x)
    @test diff(hermite(n,x), x) == 2n*hermite(n - 1, x)
    @test hermite(n, -x) == (-1)^n*hermite(n, x)

    @test laguerre(0, x) == 1
    @test laguerre(1, x) == -x + 1
    @test laguerre(2, x) == x^2/2 - 2x + 1
    @test laguerre(3, x) == -x^3/6 + 3x^2/2 - 3x + 1
    @test diff(laguerre(n, x), x) == -assoc_laguerre(n - 1, 1, x)

    @test assoc_laguerre(0, a, x) == 1
    @test assoc_laguerre(1, a, x) == a - x + 1
    @test assoc_laguerre(2, a, x) == a^2/2 + 3a/2 + x^2/2 + x*(-a - 2) + 1
    @test assoc_laguerre(n, 0, x) == laguerre(n, x)
    @test diff(assoc_laguerre(n, a, x), x) == -assoc_laguerre(n - 1, a + 1, x)

    @test Ynm(n, m, θ, ϕ) == Ynm(n, m, θ, ϕ)
    @test Ynm(n, -m, θ, ϕ) == (-1)^m*exp(-2im*m*ϕ)*Ynm(n, m, θ, ϕ)
    @test Ynm(n, m, -θ, ϕ) == Ynm(n, m, θ, ϕ)
    @test Ynm(n, m, θ, -ϕ) == exp(-2im*m*ϕ)*Ynm(n, m, θ, ϕ)
    @test expand(simplify(Ynm(0, 0, θ, ϕ)), func=true) == 1/(2*sqrt(PI))

    @test Ynm_c(n, m, θ, ϕ) == conjugate(Ynm(n, m, θ, ϕ))

    @test diff(hankel1(n, x), x) == hankel1(n - 1, x)/2 - hankel1(n + 1, x)/2
    @test diff(hankel2(n, x), x) == hankel2(n - 1, x)/2 - hankel2(n + 1, x)/2

    @test expand_func(jn(0, x)) == sin(x)/x
    @test expand_func(jn(1, x)) == sin(x)/x^2 - cos(x)/x
    @test rewrite(jn(ν, x), "besselj") == sqrt(2PI/x)*besselj(ν + Sym(1)/2, x)/2
    VERSION < v"0.6.0-dev" && @test rewrite(jn(ν, x), "bessely") == (-1)^ν*sqrt(2PI/x)*bessely(-ν - Sym(1)/2, x)/2
    u = N(jn(2, 5.2+0.3im), 20)
    @test norm(real(u) - 0.099419756723640344491) <= 1e-15 && norm(imag(u) + 0.054525080242173562897) <= 1e-15

    @test expand_func(yn(0, x)) == -cos(x)/x
    @test expand_func(yn(1, x)) == -cos(x)/x^2-sin(x)/x
    VERSION < v"0.6.0-dev" && @test rewrite(yn(ν, x), "besselj") == (-1)^(ν + 1)*sqrt(2PI/x)*besselj(-ν - Sym(1)/2, x)/2
    @test rewrite(yn(ν, x), "bessely") == sqrt(2PI/x)*bessely(ν + Sym(1)/2, x)/2
    @test N(yn(2, 5.2+0.3im)) ≈ 0.18525034196069722536 + 0.014895573969924817587im

    # gamma, beta and related functions
    @test gamma(Sym(1)) == 1
    @test gamma(Sym(3)/2) == sqrt(PI)/2

    @test diff(polygamma(Sym(0), x), x) == polygamma(Sym(1), x)
    @test diff(polygamma(Sym(0), x), x, 2) == polygamma(Sym(2), x)

    @test diff(beta(x, y), x) == (polygamma(Sym(0), x) - polygamma(Sym(0), x + y)) * beta(x, y)
    @test diff(beta(x, y), y) == (polygamma(Sym(0), y) - polygamma(Sym(0), x + y)) * beta(x, y)

    # test numerical consistency with Julia functions
    @test N(gamma(Sym(4.1))) ≈ gamma(4.1)
    @test N(polygamma(Sym(2), Sym(3.2))) ≈ polygamma(2, 3.2)
    VERSION >= v"0.5.0" && @test N(beta(Sym(1)+1im, Sym(1)+1im)) ≈ beta(1.0+1im, 1.0+1im)

    # Elliptic-type functions

    @test elliptic_k(Sym(0)) == PI/2
    @test N(elliptic_k(Sym(1.0 + im))) ≈ 1.50923695405127 + 0.625146415202697im

    @test N(elliptic_f(Sym(3.0 + im/2), Sym(1.0 + im))) ≈ 2.909449841483 + 1.74720545502474im

    @test elliptic_e(Sym(0)) == PI/2
    @test N(elliptic_e(Sym(2.0 - im))) ≈ 0.991052601328069 + 0.81879421395609im

    @test elliptic_pi(Sym(0), Sym(0)) == PI/2
    @test N(elliptic_pi(Sym(1.0 - im/3), Sym(2.0 + im))) ≈ 3.29136443417283 + 0.32555634906645im

    # Bessel-type functions
    @test diff(besselj(n, x), x) == (besselj(n - 1, x) - besselj(n + 1, x))/2
    @test rewrite(besselj(n, x), "jn") == sqrt(2x/PI)*jn(n - 1/2, x)

    @test diff(bessely(n, x), x) == (bessely(n - 1, x) - bessely(n + 1, x))/2
    @test rewrite(bessely(n, x), "yn") == sqrt(2x/PI)*yn(n - 1/2, x)

    @test diff(besseli(n, x), x) == (besseli(n - 1, x) + besseli(n + 1, x))/2

    @test diff(besselk(n, x), x) == -(besselk(n - 1, x) + besselk(n + 1, x))/2

    # test numerical consistency with Julia functions
    if VERSION < v"0.6.0-dev"
      @test N(besselj(3.2, Sym(1.5))) ≈ besselj(3.2, 1.5)
      @test N(bessely(3.2, Sym(1.5))) ≈ bessely(3.2, 1.5)
      @test N(besseli(3.2, Sym(1.5))) ≈ besseli(3.2, 1.5)
      @test N(besselk(3.2, Sym(1.5))) ≈ besselk(3.2, 1.5)
    end

    @test expand_func(x*hyper([1, 1], [2], -x)) == log(x + 1)
end


text:
using Metamath
using Base.Test

getenvsig(env) = map(f->length(getfield(env,f)),fieldnames(env))
mmexample(fname) = joinpath(Pkg.dir("Metamath"),"data",fname)

env = Metamath.main(mmexample("demo.mm"))
@test getenvsig(env) == [1,0,0,6,5,3,12,1,0,0,1]

@test env == Metamath.globalenv

empty!(env)
@test getenvsig(env) == [0,0,0,0,0,0,0,1,0,0,1]

mmverify!(env,mmexample("hol.mm"))
@test getenvsig(env) == [1,0,0,31,289,18,209,1,0,0,1]

empty!(env)
mmverify!(env,mmexample("miu.mm"))
@test getenvsig(env) == [1,0,0,5,6,2,11,1,0,0,1]


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

"""
    CompositeTask(tasks, [adjacency])

A learning task composed of multiple sub-`tasks`.
The tasks are assumed to be independent, unless
an `adjacency` matrix is provided with the list
of dependencies.

## Examples

A composite task composed of three independent sub-tasks
`t₁`, `t₂` and `t₃`:

```julia
julia> CompositeTask([t₁,t₂,t₃])
```

A composite sequential task `t₁ → t₂ → t₃`:

```julia
julia> CompositeTask([t₁,t₂,t₃], [0 1 0; 0 0 1; 0 0 0])
```

### Notes

Only directed acyclic graphs (DAGs) are permitted.
"""
struct CompositeTask <: AbstractLearningTask
  tasks::Vector{AbstractLearningTask}
  adjacency::Matrix{Int}

  function CompositeTask(tasks, adjacency)
    @assert all(length(tasks) .== size(adjacency)) "invalid adjacency for tasks"
    @assert isacyclic(adjacency) "dependencies must be acyclic"
    new(tasks, adjacency)
  end
end

CompositeTask(tasks::AbstractVector) =
  CompositeTask(tasks, zeros(Int, length(tasks), length(tasks)))

iscomposite(task::CompositeTask) = true

"""
    isacyclic(adjacency)

Check whether or not the `adjacency` matrix is acyclic
using Kahn's topological sort.
"""
function isacyclic(adjacency::AbstractMatrix{Int})
  # copy input to avoid side effects
  A = copy(adjacency)

  # find root nodes
  roots = Vector{Int}()
  for j in 1:size(A, 2)
    if all(A[:,j] .== 0)
      push!(roots, j)
    end
  end

  # Kahn's algorithm
  sorted = Vector{Int}()
  while !isempty(roots)
    i = pop!(roots)
    push!(sorted, i)

    # for all edges i → j
    for j in findall(A[i,:] .== 1)
      # remove edge i → j
      A[i,j] = 0

      # if j has no other incoming edge
      if all(A[:,j] .== 0)
        push!(sorted, j)
      end
    end
  end

  any(A .== 1) ? false : true
end


text:
using OffsetArrays

struct BernsteinBasis{T,N} <: PolynomialBasis{T,N}
    x::OffsetArray{T,1,Array{T,1}}

    function BernsteinBasis{T,N}(x) where {T,N}
        @assert length(x) == N
        ox  = OffsetArray{T}(undef, 0:N-1)
        for i in eachindex(x)
            ox[i-1] = x[i]
        end
        new(ox)
    end
end

function BernsteinBasis(x::Vector{T}) where {T}
    BernsteinBasis{T,length(x)}(x)
end

nodes(b::BernsteinBasis{T,N})  where {T,N} = b.x.parent

Base.hash(b::BernsteinBasis, h::UInt) = hash(b.c, h)

Base.:(==)(b1::BernsteinBasis, b2::BernsteinBasis) = (b1.x == b2.x)

Base.isequal(b1::BernsteinBasis{T1,N1}, b2::BernsteinBasis{T2,N2}) where {T1,N1,T2,N2} = (b1 == b2 && T1 == T2 && N1 == N2)

function bernstein(b::BernsteinBasis{T,N}, i::Int, n::Int, x::T) where {T,N}
    if i < 0 || i > n
        return zero(T)
    else
        if n == 0
            return one(T)
        else
            return bernstein(b, i, n-1, x) * (1-x) + bernstein(b, i-1, n-1, x) * x
        end
    end
end

function eval_basis(b::BernsteinBasis{T,N}, i::Int, x::T) where {T,N}
    @assert i ≥ 1 && i ≤ N
    bernstein(b, i-1, N-1, x)
end

function deriv_basis(b::BernsteinBasis{T,N}, i::Int, x::T) where {T,N}
    @assert i ≥ 1 && i ≤ N
    (N-1) * ( bernstein(b, i-2, N-2, x) - bernstein(b, i-1, N-2, x) )
end

deriv_basis(b::BernsteinBasis, i::Int, j::Int) = derivative(b, i, b.x[j])

# function int_basis(b::BernsteinBasis{T,N}, i::Int, x::T) where {T,N}
# end
#
# int_basis(b::BernsteinBasis, i::Int, j::Int) = integral(b, i, b.x[j])


text:
using Clang
using Clang.Generators
using Telescope_jll

cd(@__DIR__)

include_dir = joinpath(Telescope_jll.artifact_dir, "include") |> normpath
telescope_h = joinpath(include_dir, "telescope.h")

options = load_options(joinpath(@__DIR__, "generator.toml"))

args = get_default_args()
push!(args, "-I$include_dir")

ctx = create_context([telescope_h], args, options)

build!(ctx)


text:
function log_one_node(com, nvars, back_idx, step_nr)
    tree_log_node = TreeLogNode()
    tree_log_node.step_nr = step_nr
    tree_log_node.id = back_idx
    parent_idx = 0
    if length(com.backtrack_vec) > 0
        tree_log_node.status = com.backtrack_vec[back_idx].status
        tree_log_node.var_idx = com.backtrack_vec[back_idx].variable_idx
        tree_log_node.set_val = com.backtrack_vec[back_idx].pval
        tree_log_node.best_bound = com.backtrack_vec[back_idx].best_bound

        parent_idx = com.backtrack_vec[back_idx].parent_idx
    else # for initial solve
        tree_log_node.status = :Closed
        tree_log_node.var_idx = 0
        tree_log_node.set_val = 0
        tree_log_node.best_bound = com.best_bound
    end
    tree_log_node.var_changes = Dict{Int64,Vector{Tuple{Symbol, Int64, Int64, Int64}}}()
    tree_log_node.var_states = Dict{Int64,Vector{Int64}}()
    if tree_log_node.status == :Closed
        for var in com.search_space
            if length(var.changes[back_idx]) > 0
                tree_log_node.var_states[var.idx] = sort!(values(var))
                tree_log_node.var_changes[var.idx] = var.changes[back_idx]
            end
        end
    end
    tree_log_node.children = Vector{TreeLogNode}()
    if parent_idx > 0
        changed = false
        for (i,child) in enumerate(com.logs[parent_idx].children)
            if child.id == back_idx
                com.logs[parent_idx].children[i] = tree_log_node
                changed = true
                break
            end
        end
        if !changed
            push!(com.logs[parent_idx].children, tree_log_node)
        end
    end
    return tree_log_node
end

"""
    save_logs(com::CS.CoM, filepath)

Save the tree structure and some additional problem information in a json file `filepath`.
Can be only used if `keep_logs` is set to `true` in the [`solve!`](@ref) call.
"""
function save_logs(com::CS.CoM, filepath)
    logs = Dict{Symbol, Any}()
    nvars = length(com.search_space)

    logs[:init] = Vector{Vector{Int64}}(undef, nvars)
    for var in com.init_search_space
        logs[:init][var.idx] = values(var)
    end

    logs[:tree] = com.logs[1]
    write(filepath, JSON.json(logs))
end


text:
#!/cvmfs/soft.computecanada.ca/easybuild/software/2020/avx2/Core/julia/1.6.0/bin/julia
#SBATCH --mail-user=chunlok@ualberta.ca
#SBATCH --mail-type=ALL
#SBATCH -o emphatic_control.out # Standard output
#SBATCH -e emphatic_control.err # Standard error
#SBATCH --mem-per-cpu=4000M # Memory request of 3 GB
#SBATCH --time=09:30:00 #
#SBATCH --ntasks=164
#SBATCH --account=def-whitem

using Pkg
Pkg.activate(".")

include(joinpath(ENV["SLURM_SUBMIT_DIR"], "parallel/parallel_config.jl"))
reproduce_config_experiment("configs/interestTB_test/experiment_4_gpi.toml", save_path="~/scratch/curiosity")

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
type RsvgHandle
	ptr::Ptr{Void}

    function RsvgHandle(ptr::Ptr{Void})
        self = new(ptr)
        #finalizer(self, destroy)
        self
    end
end

function destroy(handle::RsvgHandle)
    if handle.ptr == C_NULL
        return
    end
    Gtk.GLib.gc_unref(handle)
    handle.ptr = C_NULL
    nothing
end

GError = Gtk.GLib.GError

type RsvgDimensionData
    width::Int32
    height::Int32
    em::Float64
    ex::Float64
end

RsvgDimensionData() = RsvgDimensionData(0,0,0,0)


text:
module BLAS

using LinearAlgebra, SparseArrays
using LinearAlgebra: BlasInt, checksquare

# For testing purposes:
global const __counter = Ref(0)

# TODO: BLAS1

# BLAS and BLAS3
include(joinpath("level_2_3", "matdescra.jl"))
include(joinpath("level_2_3", "generator.jl"))
include(joinpath("level_2_3", "matmul.jl"))

end


text:
struct SieveLayer{T₁ <: Number, T₂ <: KDEComparator}
    input_weights::Matrix{T₁}
    comparators::Vector{T₂}

    function SieveLayer(input_weights::Matrix{T₁},
                        comparators::Vector{T₂}) where {T₁ <: Number,
                                                        T₂ <: KDEComparator}
        new{T₁,T₂}(input_weights, comparators)
    end
end

abstract type Sieve end

struct DataPassingSieve{T <: Number} <: Sieve
    n_neurons::Int
    layers::Vector{SieveLayer{T}}
    output_weights::Matrix{T}
    consensus_threshold::Real
    training_samples_per_layer::Vector{Int}
    depth_per_sample::Vector{Int}

    function DataPassingSieve{T₁}(samples::T₂,
                                  targets::T₃,
                                  n_neurons::Int;
                                  initial_boundary_offset::Int = 1,
                                  subsequent_boundary_offset::Int = 1,
                                  consensus_threshold::Real = 0.9,
                                  max_layers::Int = 3,
                                  min_split::Int = 2) where {T₁ <: Number,
                                                             T₂ <: AbstractMatrix,
                                                             T₃ <: AbstractVector}
        X = samples
        t = targets
        L = n_neurons
        N = last(size(X))
        H = zeros(T₁, L * max_layers, N)
        layers = Vector{SieveLayer{T₁}}()
        training_samples_per_layer = Vector{Int}()
        depth_per_sample = zeros(Int, N)

        X₀ = X
        t₀ = t
        active_samples = Set{Int}(1:N)
        indices = get_indices(active_samples, N)

        for depth in 1:max_layers
            if length(active_samples) < min_split || all(t .== 1) || all(t .!= 1)
                break
            end

            push!(training_samples_per_layer, length(active_samples))
            depth_per_sample[indices] .= depth
            boundary_offset = depth == 1 ? initial_boundary_offset : subsequent_boundary_offset

            i₀ = (depth-1) * L + 1
            i₁ = i₀ + L - 1

            H₀ = @view H[i₀:i₁,indices]
            D, N₀ = size(X₀)

            W, fs = train_kde_comparators(T₁, X₀, t₀, L, boundary_offset = boundary_offset)
            H₀ .= project(X₀, W, fs)
            push!(layers, SieveLayer(W, fs))

            to_remove = Set{Int}()
            for n in active_samples
                if consensus(H[i₀:i₁,n], consensus_threshold)
                    push!(to_remove, n)
                end
            end
            setdiff!(active_samples, to_remove)
            indices = get_indices(active_samples, N)

            X₀ = @view X[:,indices]
            t₀ = @view t[indices]
        end

        H = [H; ones(eltype(H), 1, N)]

        Ψ = LinearAlgebra.Diagonal(calculate_sample_weights(t))
        H = H * Ψ
        T = reshape(t, 1, :) * Ψ

        B = (T * H') * LinearAlgebra.pinv(H * H')

        new{T₁}(L, layers, B, consensus_threshold, training_samples_per_layer, depth_per_sample)
    end
end

function predict(sieve::T₁,
                 samples::T₂) where {T₁ <: DataPassingSieve,
                                     T₂ <: AbstractMatrix}
    param(::DataPassingSieve{T}) where {T} = T

    X = samples
    L = sieve.n_neurons
    N = last(size(X))
    H = zeros(param(sieve), L * length(sieve.layers), N)

    X₀ = X
    active_samples = Set{Int}(1:N)
    indices = get_indices(active_samples, N)

    for (depth, layer) in enumerate(sieve.layers)
        i₀ = (depth-1) * L + 1
        i₁ = i₀ + L - 1

        H₀ = @view H[i₀:i₁,indices]
        D, N₀ = size(X₀)

        W = layer.input_weights
        fs = layer.comparators
        H₀ .= project(X₀, W, fs)

        to_remove = Set{Int}()
        for n in active_samples
            if consensus(H[i₀:i₁,n], sieve.consensus_threshold)
                push!(to_remove, n)
            end
        end
        setdiff!(active_samples, to_remove)
        indices = get_indices(active_samples, N)

        X₀ = @view X[:,indices]
    end

    H = [H; ones(eltype(H), 1, N)]
    B = sieve.output_weights
    y = vec(B * H)
end

function predict(sieve::T₁,
                 sample::T₂) where {T₁ <: DataPassingSieve,
                                    T₂ <: AbstractVector}
    y = predict(sieve, reshape(sample, :, 1))
    first(y)
end

struct ProjectionPassingSieve{T <: Number} <: Sieve
    n_neurons::Int
    layers::Vector{SieveLayer{T}}
    output_weights::Matrix{T}
    consensus_threshold::Real
    training_samples_per_layer::Vector{Int}
    depth_per_sample::Vector{Int}

    function ProjectionPassingSieve{T₁}(samples::T₂,
                                        targets::T₃,
                                        n_neurons::Int;
                                        initial_boundary_offset::Int = 1,
                                        subsequent_boundary_offset::Int = 1,
                                        consensus_threshold::Real = 0.9,
                                        max_layers::Int = 3,
                                        min_split::Int = 2) where {T₁ <: Number,
                                                                   T₂ <: AbstractMatrix,
                                                                   T₃ <: AbstractVector}
        X = samples
        t = targets
        L = n_neurons
        N = last(size(X))
        H = zeros(T₁, L * max_layers, N)
        layers = Vector{SieveLayer{T₁}}()

        X₀ = X
        t₀ = t
        active_samples = Set{Int}(1:N)
        indices = get_indices(active_samples, N)

        for depth in 1:max_layers
            if length(active_samples) < min_split || all(t .== 1) || all(t .!= 1)
                break
            end

            push!(training_samples_per_layer, length(active_samples))
            depth_per_sample[indices] .= depth
            boundary_offset = depth == 1 ? initial_boundary_offset : subsequent_boundary_offset

            i₀ = (depth-1) * L + 1
            i₁ = i₀ + L - 1

            H₀ = @view H[i₀:i₁,indices]
            D, N₀ = size(X₀)

            W, fs = train_kde_comparators(T₁, X₀, t₀, L, boundary_offset = boundary_offset)
            H₀ .= project(X₀, W, fs)
            push!(layers, SieveLayer(W, fs))

            to_remove = Set{Int}()
            for n in active_samples
                if consensus(H[i₀:i₁,n], consensus_threshold)
                    push!(to_remove, n)
                end
            end
            setdiff!(active_samples, to_remove)
            indices = get_indices(active_samples, N)

            X₀ = @view H[i₀:i₁,indices]
            t₀ = @view t[indices]
        end

        H = [H; ones(eltype(H), 1, N)]

        Ψ = LinearAlgebra.Diagonal(calculate_sample_weights(t))
        H = H * Ψ
        T = reshape(t, 1, :) * Ψ

        B = (T * H') * LinearAlgebra.pinv(H * H')

        new{T₁}(L, layers, B, consensus_threshold, training_samples_per_layer, depth_per_sample)
    end
end

function predict(sieve::T₁,
                 samples::T₂) where {T₁ <: ProjectionPassingSieve,
                                     T₂ <: AbstractMatrix}
    param(::ProjectionPassingSieve{T}) where {T} = T

    X = samples
    L = sieve.n_neurons
    N = last(size(X))
    H = zeros(param(sieve), L * length(sieve.layers), N)

    X₀ = X
    active_samples = Set{Int}(1:N)
    indices = get_indices(active_samples, N)

    for (depth, layer) in enumerate(sieve.layers)
        i₀ = (depth-1) * L + 1
        i₁ = i₀ + L - 1

        H₀ = @view H[i₀:i₁,indices]
        D, N₀ = size(X₀)

        W = layer.input_weights
        fs = layer.comparators
        H₀ .= project(X₀, W, fs)

        to_remove = Set{Int}()
        for n in active_samples
            if consensus(H[i₀:i₁,n], sieve.consensus_threshold)
                push!(to_remove, n)
            end
        end
        setdiff!(active_samples, to_remove)
        indices = get_indices(active_samples, N)

        X₀ = @view H[i₀:i₁,indices]
    end

    H = [H; ones(eltype(H), 1, N)]
    B = sieve.output_weights
    y = vec(B * H)
end

function predict(sieve::T₁,
                 sample::T₂) where {T₁ <: ProjectionPassingSieve,
                                    T₂ <: AbstractVector}
    y = predict(sieve, reshape(sample, :, 1))
    first(y)
end

function consensus(samples::T₁,
                   consensus_threshold::Real) where {T₁ <: AbstractVector}
    y = sign.(samples)
    th = consensus_threshold
    f(x) = x == 0 || x == 1
    g(x) = x == 0 || x == -1

    all(f.(y)) || all(g.(y))
end

function get_indices(active_samples::Set{Int}, N::Int)
    [n in active_samples for n in 1:N]
end


text:
function PhaseIonosphericDelayModifier(arg0::IonosphericModel, arg1::jdouble)
    return PhaseIonosphericDelayModifier((IonosphericModel, jdouble), arg0, arg1)
end

function get_parameters_drivers(obj::PhaseIonosphericDelayModifier)
    return jcall(obj, "getParametersDrivers", List, ())
end

function modify(obj::PhaseIonosphericDelayModifier, arg0::EstimatedMeasurement)
    return jcall(obj, "modify", void, (EstimatedMeasurement,), arg0)
end


text:
module UserExceptionModule

struct UserError <: Exception end

end

max_stars_repo_name: ThomasRiddick/DynamicHD

max_issues_repo_name: ThomasRiddick/DynamicHD

max_forks_repo_name: ThomasRiddick/DynamicHD


text:
using Symbolics: scalarize
@variables x[1:3] u[1:2]
@parameters p[1:2] t

g = [x[1] * p[1] + p[2] * x[2]; x[2] * u[1]; u[2] * x[3] + exp(-t)]

b = Basis(g, x, parameters = p, controls = u, iv = t)

# Vector call
x0 = randn(3)
p0 = randn(2)
t0 = 0.0
u0 = randn(2)

true_res(x, p, t, u) =
    [sum(x[1:2] .* p); x[2] .* u[1]; u[2] .* x[3] .+ exp.(-t)]
true_res_(x, p, t) =
    scalarize.(hcat([true_res(x[:, i], p, t[i], zeros(2)) for i = 1:100]...))
true_res_(x, p, t, u) =
    scalarize.(hcat([true_res(x[:, i], p, t[i], u[:, i]) for i = 1:100]...))

@test isequal(b(x0), DataDrivenDiffEq.get_f(b)(x0, p, t, zeros(2)))
@test isequal(b(x0, p0), DataDrivenDiffEq.get_f(b)(x0, p0, t, zeros(2)))
@test isequal(b(x0, p0, t0), DataDrivenDiffEq.get_f(b)(x0, p0, t0, zeros(2)))
@test isequal(b(x0, p0, t0, u0), DataDrivenDiffEq.get_f(b)(x0, p0, t0, u0))

# Array call
x0 = randn(3, 100)
p0 = randn(2)
t0 = randn(100)
u0 = randn(2, 100)

# These first two fail, since exp(-t) != exp(getindex(t,1))
@test isequal(b(x0), true_res_(x0, p, [t for i = 1:100]))
@test isequal(b(x0, p0), true_res_(x0, p0, [t for i = 1:100]))
@test isequal(b(x0, p0, t0), true_res_(x0, p0, t0))
@test isequal(b(x0, p0, t0, u0), true_res_(x0, p0, t0, u0))

@parameters w[1:2] t
@variables u[1:3](t)

h = [u[1]; u[2]; cos(w[1] * u[2] + w[2] * u[3]); 5*u[3] + u[2]]
h_not_unique = [u[1]; u[1]; u[1]^1; h; 1]
basis = Basis(h_not_unique, u, parameters = w, iv = t)
basis_2 =
    Basis(h_not_unique, u, parameters = w, iv = t, linear_independent = true)
# Check getters
@test isequal(states(basis), u)
@test isequal(parameters(basis), w)
@test isequal(ModelingToolkit.get_iv(basis), t)
@test isequal(controls(basis), [])

# Check free parameter calculation
@test free_parameters(basis, operations = [+]) == 7
@test free_parameters(basis_2) == 6
@test free_parameters(basis, operations = [+, cos]) == 8
@test free_parameters(basis_2, operations = [+, cos]) == 7
@test DataDrivenDiffEq.count_operation(
    (1 + cos(u[2]) * sin(u[1]))^3,
    [+, cos, ^, *],
) == 4

# Check array functionalities
basis_2 = unique(basis)
@test isequal(basis, basis_2)
@test size(basis) == (5,)
@test basis([1.0; 2.0; π], [0.0; 1.0]) ≈ [1.0; 2.0; -1.0; 5*π + 2.0; 1.0]
@test size(basis) == size(basis_2)
push!(basis_2, sin(u[2]))
@test size(basis_2)[1] == length(h) + 2
basis_3 = merge(basis, basis_2)
@test size(basis_3) == (6,)
@test isequal(states(basis_3), states(basis_2))
@test isequal(parameters(basis_3), parameters(basis_2))
merge!(basis_3, basis)
@test basis_3 == basis_2
push!(basis, 5*u[3] + u[2])
unique!(basis) # Does not remove
@test size(basis) == (5,)

# Further callables
@variables a
g = [u[1]; u[3]; a]
basis = Basis(g, [u; a])
@test basis([1; 2; 3; 4]) == [1; 3; 4]
g = [u[1]; u[3]; u[2]]
basis = Basis(g, u, parameters = [])
X = ones(Float64, 3, 10)
X[1, :] .= 3 * X[1, :]
X[3, :] .= 5 * X[3, :]
# Check the array evaluation
@test basis(X) ≈ [1.0 0.0 0.0; 0.0 0.0 1.0; 0.0 1.0 0.0] * X
Y = similar(X)
basis(Y, X)
@test Y ≈ [1.0 0.0 0.0; 0.0 0.0 1.0; 0.0 1.0 0.0] * X
f = DataDrivenDiffEq.jacobian(basis)
@test f([1; 1; 1], [0.0; 0.0], 0.0) ≈ [1.0 0.0 0.0; 0.0 0.0 1.0; 0.0 1.0 0.0]
@test_nowarn [xi for xi in basis]
@test_nowarn basis[2:end];

@variables u[1:3] t

g = [u[2]; -sin(u[1]) * exp(-t); u[2] + u[3]]
basis = Basis(g, u, iv = t)

f_(u, p, t) = [u[3]; u[2] * u[1]; p[1] * sin(u[1]) * u[2]; p[2] * t]
b = Basis(f_, u, parameters = w, iv = t)
@test f_([1; 2; 3], [2; 0], 3.0) ≈ b([1; 2; 3], [2; 0], 3.0)


text:
#=
Core type definitions for the `NCDataset` struct,
the `Attributes` and `Group` parts as well.
and the `Attributes` part of it.
Helper functions defined for these as well.

High-level interface is at the "High-level" section about actually
loading/reading/making datasets.
=#

############################################################
# Type definitions
############################################################
# -----------------------------------------------------
# base type of attributes list
# concrete types are Attributes (single NetCDF file) and
# MFAttributes (multiple NetCDF files)

abstract type BaseAttributes
end

abstract type AbstractDataset
end

abstract type AbstractDimensions
end

abstract type AbstractGroups
end

# -----------------------------------------------------
# List of attributes (for a single NetCDF file)
# all ids should be Cint

mutable struct Attributes <: BaseAttributes
    ncid::Cint
    varid::Cint
    isdefmode::Vector{Bool}
end

mutable struct Groups <: AbstractGroups
    ncid::Cint
    isdefmode::Vector{Bool}
end

mutable struct Dimensions <: AbstractDimensions
    ncid::Cint
    isdefmode::Vector{Bool}
end

mutable struct NCDataset <: AbstractDataset
    ncid::Cint
    # true of the NetCDF is in define mode (i.e. metadata can be added, but not data)
    # need to be an array, so that it is copied by reference
    isdefmode::Vector{Bool}
    attrib::Attributes
    dim::Dimensions
    group::Groups
end

"Alias to `NCDataset`"
const Dataset = NCDataset

const NCIterable = Union{BaseAttributes,AbstractDimensions,AbstractDataset,AbstractGroups}
Base.length(a::NCIterable) = length(keys(a))

############################################################
# Mappings
############################################################

# Mapping between NetCDF types and Julia types
const jlType = Dict(
    NC_BYTE   => Int8,
    NC_UBYTE  => UInt8,
    NC_SHORT  => Int16,
    NC_USHORT => UInt16,
    NC_INT    => Int32,
    NC_UINT   => UInt32,
    NC_INT64  => Int64,
    NC_UINT64 => UInt64,
    NC_FLOAT  => Float32,
    NC_DOUBLE => Float64,
    NC_CHAR   => Char,
    NC_STRING => String
)

# Inverse mapping
const ncType = Dict(value => key for (key, value) in jlType)

"""
    fillvalue(::Type{Int8})
    fillvalue(::Type{UInt8})
    fillvalue(::Type{Int16})
    fillvalue(::Type{UInt16})
    fillvalue(::Type{Int32})
    fillvalue(::Type{UInt32})
    fillvalue(::Type{Int64})
    fillvalue(::Type{UInt64})
    fillvalue(::Type{Float32})
    fillvalue(::Type{Float64})
    fillvalue(::Type{Char})
    fillvalue(::Type{String})

Default fill-value for the given type.
"""
@inline fillvalue(::Type{Int8})    = NC_FILL_BYTE
@inline fillvalue(::Type{UInt8})   = NC_FILL_UBYTE
@inline fillvalue(::Type{Int16})   = NC_FILL_SHORT
@inline fillvalue(::Type{UInt16})  = NC_FILL_USHORT
@inline fillvalue(::Type{Int32})   = NC_FILL_INT
@inline fillvalue(::Type{UInt32})  = NC_FILL_UINT
@inline fillvalue(::Type{Int64})   = NC_FILL_INT64
@inline fillvalue(::Type{UInt64})  = NC_FILL_UINT64
@inline fillvalue(::Type{Float32}) = NC_FILL_FLOAT
@inline fillvalue(::Type{Float64}) = NC_FILL_DOUBLE
@inline fillvalue(::Type{Char})    = NC_FILL_CHAR
@inline fillvalue(::Type{String})  = NC_FILL_STRING

"Make sure that a dataset is in data mode"
function datamode(ncid,isdefmode::Vector{Bool})
    if isdefmode[1]
        nc_enddef(ncid)
        isdefmode[1] = false
    end
end

"Make sure that a dataset is in define mode"
function defmode(ncid,isdefmode::Vector{Bool})
    if !isdefmode[1]
        nc_redef(ncid)
        isdefmode[1] = true
    end
end

############################################################
# High-level
############################################################

"""
    NCDataset(filename::AbstractString, mode = "r";
            format::Symbol = :netcdf4, attrib = [])

Load, create, or even overwrite a NetCDF file at `filename`, depending on `mode`:

* `"r"` (default) : open an existing netCDF file or OPeNDAP URL
   in read-only mode.
* `"c"` : create a new NetCDF file at `filename` (an existing file with the same
  name will be overwritten).
* `"a"` : open `filename` into append mode (i.e. existing data in the netCDF
  file is not overwritten and a variable can be added).

Notice that this does not close the dataset, use `close` on the
result (or see below the `do`-block).

The optional parameter `attrib` is an iterable of attribute name and attribute
value pairs, for example a `Dict`, `DataStructures.OrderedDict` or simply a
vector of pairs (see example below).

# Supported `format` values:

* `:netcdf4` (default): HDF5-based NetCDF format.
* `:netcdf4_classic`: Only netCDF 3 compatible API features will be used.
* `:netcdf3_classic`: classic netCDF format supporting only files smaller than 2GB.
* `:netcdf3_64bit_offset`: improved netCDF format supporting files larger than 2GB.

Files can also be open and automatically closed with a `do` block.

```julia
NCDataset("file.nc") do ds
    data = ds["temperature"][:,:]
end
```

Here is an attribute example:
```julia
using DataStructures
NCDataset("file.nc", "c", attrib = OrderedDict("title" => "my first netCDF file")) do ds
   defVar(ds,"temp",[10.,20.,30.],("time",))
end;
```

`NCDataset` is an alias to `NCDataset`.
"""
function NCDataset(filename::AbstractString,
                 mode::AbstractString = "r";
                 format::Symbol = :netcdf4,
                 diskless = false,
                 persist = false,
                 attrib = [])

    ncid = -1
    isdefmode = [false]

    if (mode == "r") || (mode == "a")
        ncmode =
            if (mode == "r")
                NC_NOWRITE
            else
                NC_WRITE
            end

        if diskless
            ncmode = ncmode | NC_DISKLESS
        end

        ncid = nc_open(filename,ncmode)
    elseif mode == "c"
        ncmode  = NC_CLOBBER

        if format == :netcdf3_64bit_offset
            ncmode = ncmode | NC_64BIT_OFFSET
        elseif format == :netcdf4_classic
            ncmode = ncmode | NC_NETCDF4 | NC_CLASSIC_MODEL
        elseif format == :netcdf4
            ncmode = ncmode | NC_NETCDF4
        elseif format == :netcdf3_classic
            # do nothing
        else
            throw(NetCDFError(-1, "Unkown format '$(format)' for filename '$(filename)'"))
        end

        if diskless
            ncmode = ncmode | NC_DISKLESS

            if persist
                ncmode = ncmode | NC_PERSIST
            end
        end

        ncid = nc_create(filename,ncmode)
        isdefmode[1] = true
    else
        throw(NetCDFError(-1, "Unsupported mode '$(mode)' for filename '$(filename)'"))
    end

    ds = NCDataset(ncid,isdefmode)

    # set global attributes
    for (attname,attval) in attrib
        ds.attrib[attname] = attval
    end

    return ds
end

function NCDataset(ncid::Integer,
                 isdefmode::Vector{Bool})
    attrib = Attributes(ncid,NC_GLOBAL,isdefmode)
    dim = Dimensions(ncid,isdefmode)
    group = Groups(ncid,isdefmode)
    return NCDataset(ncid,isdefmode,attrib,dim,group)
end

function NCDataset(f::Function,args...; kwargs...)
    ds = NCDataset(args...; kwargs...)
    try
        f(ds)
    finally
        #@debug "closing netCDF NCDataset $(NCDatasets.path(ds))"
        close(ds)
    end
end

export NCDataset, Dataset

############################################################
# High-level: user convenience
############################################################
"""
    keys(ds::NCDataset)

Return a list of all variables names in NCDataset `ds`.
"""
Base.keys(ds::NCDataset) = listVar(ds.ncid)

"""
    path(ds::NCDataset)
Return the file path (or the opendap URL) of the NCDataset `ds`
"""
path(ds::NCDataset) = nc_inq_path(ds.ncid)
export path

"""
    sync(ds::NCDataset)

Write all changes in NCDataset `ds` to the disk.
"""
function sync(ds::NCDataset)
    datamode(ds.ncid,ds.isdefmode)
    nc_sync(ds.ncid)
end
export sync

"""
    close(ds::NCDataset)

Close the NCDataset `ds`. All pending changes will be written
to the disk.
"""
Base.close(ds::NCDataset) = nc_close(ds.ncid)
export close

############################################################
# Common methods
############################################################
function Base.iterate(a::NCIterable, state = keys(a))
    if length(state) == 0
        return nothing
    end

    return (state[1] => a[popfirst!(state)], state)
end

"""
    varbyattrib(ds, attname = attval)

Returns a list of variable(s) which has the attribute `attname` matching the value `attval`
in the dataset `ds`.
The list is empty if the none of the variables has the match.
The output is a list of `CFVariable`s.

# Examples

Load all the data of the first variable with standard name "longitude" from the
NetCDF file `results.nc`.

```julia-repl
julia> ds = NCDataset("results.nc", "r");
julia> data = varbyattrib(ds, standard_name = "longitude")[1][:]
```

"""
function varbyattrib(ds::NCDataset; kwargs...)
    # Start with an empty list of variables
    varlist = []

    # Loop on the variables
    for v in keys(ds)
        var = ds[v]

        matchall = true

        for (attsym,attval) in kwargs
            attname = String(attsym)

            # Check if the variable has the desired attribute
            if haskey(var.attrib, attname)
                # Check if the attribute value is the selected one
                if var.attrib[attname] != attval
                    matchall = false
                    break
                end
            else
                matchall = false
                break
            end
        end

        if matchall
            push!(varlist, var)
        end
    end

    return varlist
end
export varbyattrib

"""
    haskey(ds::NCDataset,name)
    haskey(d::Dimensions,name)
    haskey(ds::Attributes,name)

Return true if the NCDataset `ds` (or dimension/attribute list) has a variable (dimension/attribute) with the name `name`.
For example:

```julia
ds = NCDataset("/tmp/test.nc","r")
if haskey(ds,"temperature")
    println("The file has a variable 'temperature'")
end

if haskey(ds.dim,"lon")
    println("The file has a dimension 'lon'")
end
```

This example checks if the file `/tmp/test.nc` has a variable with the
name `temperature` and a dimension with the name `lon`.
"""
Base.haskey(a::NCIterable,name::AbstractString) = name in keys(a)
Base.in(name::AbstractString,a::NCIterable) = name in keys(a)

function Base.show(io::IO,ds::AbstractDataset; indent="")
    try
        dspath = path(ds)
        printstyled(io, indent, "NCDataset: ",dspath,"\n", color=:red)
    catch err
        if isa(err,NetCDFError)
            if err.code == NC_EBADID
                print(io,"closed NetCDF NCDataset")
                return
            end
        end
        rethrow
    end

    print(io,indent,"Group: ",groupname(ds),"\n")
    print(io,"\n")

    dims = collect(ds.dim)

    if length(dims) > 0
        printstyled(io, indent, "Dimensions\n",color=:red)

        for (dimname,dimlen) in dims
            print(io,indent,"   $(dimname) = $(dimlen)\n")
        end
        print(io,"\n")
    end

    varnames = keys(ds)

    if length(varnames) > 0

        printstyled(io, indent, "Variables\n",color=:red)

        for name in varnames
            show(io,variable(ds,name); indent = "$(indent)  ")
            print(io,"\n")
        end
    end

    # global attribues
    if length(ds.attrib) > 0
        printstyled(io, indent, "Global attributes\n",color=:red)
        show(io,ds.attrib; indent = "$(indent)  ");
    end

    # groups
    groupnames = keys(ds.group)

    if length(groupnames) > 0
        printstyled(io, indent, "Groups\n",color = :red)
        for groupname in groupnames
            show(io,group(ds,groupname); indent = "  ")
        end
    end

end

"""
    write(dest_filename::AbstractString, src::AbstractDataset; include = keys(src), exclude = [])
    write(dest::NCDataset, src::AbstractDataset; include = keys(src), exclude = [])

Write the variables of `src` dataset into an empty `dest` dataset (which must be opened in mode `"a"` or `"c"`).
The keywords `include` and `exclude` configure which variable of `src` should be included
(by default all), or which should be `excluded` (by default none).

If the first argument is a file name, then the dataset is open in create mode (`c`).

This function is useful when you want to save the dataset from a multi-file dataset.
"""
function Base.write(dest::NCDataset, src::AbstractDataset;
                     include = keys(src), exclude = String[])

    unlimited_dims = unlimited(src.dim)

    for (dimname,dimlength) in src.dim
        isunlimited = dimname in unlimited_dims

        # if haskey(dest.dim,dimname)
        #     # check length
        #     if (dest.dim[dimname] !== src.dim[dimname]) && !isunlimited
        #         throw(DimensionMismatch("length of the dimensions $dimname are inconstitent in files $(path(dest)) and $(path(src))"))
        #     end
        # else
            if isunlimited
                defDim(dest, dimname, Inf)
            else
                defDim(dest, dimname, dimlength)
            end
        # end
    end

    # loop over variables
    for varname in include
        (varname ∈ exclude) && continue
        @debug "Writing variable $varname..."
        cfvar = src[varname]
        defVar(dest, varname, Array(cfvar), dimnames(cfvar); attrib = cfvar.attrib)
    end

    # loop over all global attributes
    for (attribname,attribval) in src.attrib
        dest.attrib[attribname] = attribval
    end

    # loop over all groups
    for (groupname,groupsrc) in src.group
        groupdest = defGroup(dest,groupname)
        write(groupdest,groupsrc)
    end
    return dest
end

function Base.write(dest_filename::AbstractString, src::AbstractDataset; kwargs...)
    NCDataset(dest_filename,"c") do dest
        write(dest,src, kwargs...)
    end
    return nothing
end


text:
module AddressEvent

""" Event timestamp in us. """
function event_time end

""" Raw event coordinates (x,y). """
function event_coord end

""" Event location (i,j) in standard image coordinates """
function event_location end

""" Event polarity. """
function event_pol end

export event_time, event_coord, event_pol, event_location

function read_events end
function read_events! end

""" Groups `events` into frames of duration `period_us`. """
function iter_frames(f::Function, events, period_us::Integer)
    E = eltype(events)
    frame = E[]
    t0 = -1
    for event in events
        if t0 == -1
            t0 = Int(event_time(event))
        end
        while (event_time(event) - t0) >= period_us
            f(frame)
            frame = E[]
            t0 += Int(period_us)
        end
        push!(frame, event)
    end
    f(frame)
end

function iter_frames(events, period_us::Integer)
    out = Channel{Vector{eltype(events)}}(0)
    @schedule try
        iter_frames(events, period_us) do frame
            put!(out, frame)
        end
    finally
        close(out)
    end
    out
end

function read_frames(events, period_us::Integer)
    frames = Vector{eltype(events)}[]
    iter_frames(events, period_us) do frame
        push!(frames, frame)
    end
    frames
end

export read_events, read_events!, iter_frames, read_frames

""" Standard image size for a particular event type """
function image_size(E)
    error("No default image size for $E")
end

""" Draw events at their location on a 2D canvas """
function draw_events!(events, canvas::Matrix)
    for e in events
        i, j = event_location(e, size(canvas))
        canvas[i,j] += event_pol(e) ? 1 : -1
    end
end

function draw_events(events, imsize)
    canvas = zeros(imsize)
    draw_events!(events, canvas)
    canvas
end

function draw_events(events::Vector{E}) where {E}
    draw_events(events, image_size(E))
end

""" Draws a film strip with the supplied frames """
function film_strip(frames::Vector{<:Vector{<:E}}, imsize=image_size(E)) where {E}
    imgs = hcat([draw_events(events, zeros(imsize)) for events in frames]...)
end

using Images
using StatsBase

""" Turns rasterized events into an image, ignoring hot pixels """
gray_img(img) = Gray.(0.5 .+ img./max(1,2*quantile(img[:],0.99)))

""" Exports events as a video. """
function export_video(filename::String, events, period_us::Number, imsize=image_size(events))
    open(`ffmpeg -f image2pipe -vcodec png -r $(1000000/period_us) -i - -vcodec h264 -pix_fmt yuv420p $filename`, "w") do io
        for frame in iter_frames(events, period_us)
            img = gray_img(draw_events(frame, imsize))
            show(io, MIME("image/png"), img)
        end
    end
end

export image_size, draw_events, draw_events!, film_strip, export_video

include("aer_dat.jl")
include("ros_dvs.jl")

end # module


text:
using Oceananigans.Architectures
using Oceananigans.BoundaryConditions
using Oceananigans.Utils

"""
    HorizontalAverage{F, R, P, I, Ω} <: AbstractDiagnostic

A diagnostic for computing horizontal average of a field.
"""
mutable struct HorizontalAverage{F, R, P, I, Ω, G} <: AbstractDiagnostic
          field :: F
         result :: P
      frequency :: Ω
       interval :: I
       previous :: Float64
    return_type :: R
           grid :: G
end

"""
    HorizontalAverage(model, field; frequency=nothing, interval=nothing, return_type=Array)

Construct a `HorizontalAverage` of `field`.

After the horizontal average is computed it will be stored in the `result` property.

The `HorizontalAverage` can be used as a callable object that computes and returns the
horizontal average.

A `frequency` or `interval` (or both) can be passed to indicate how often to run this
diagnostic if it is part of `model.diagnostics`. `frequency` is a number of iterations
while `interval` is a time interval in units of `model.clock.time`.

A `return_type` can be used to specify the type returned when the `HorizontalAverage` is
used as a callable object. The default `return_type=Array` is useful when running a GPU
model and you want to save the output to disk by passing it to an output writer.
"""
function HorizontalAverage(field; frequency=nothing, interval=nothing, return_type=Array)
    arch = architecture(field)
    result = zeros(arch, field.grid, 1, 1, field.grid.Tz)
    return HorizontalAverage(field, result, frequency, interval, 0.0, return_type, field.grid)
end

# Normalize a horizontal sum to get the horizontal average.
normalize_horizontal_sum!(havg, grid) = havg.result /= (grid.Nx * grid.Ny)

"""
    run_diagnostic(model, havg::HorizontalAverage{NTuple{1}})

Compute the horizontal average of `havg.field` and store the result in `havg.result`.
"""
function run_diagnostic(model, havg::HorizontalAverage)
    zero_halo_regions!(parent(havg.field), model.grid)
    sum!(havg.result, parent(havg.field))
    normalize_horizontal_sum!(havg, model.grid)
    return nothing
end

function (havg::HorizontalAverage{F, Nothing})(model) where F
    run_diagnostic(model, havg)
    return havg.result
end

function (havg::HorizontalAverage)(model)
    run_diagnostic(model, havg)
    return havg.return_type(havg.result)
end


text:
module MatrixProductStates

# Basic functionality for contracting tensors in an efficient manner
include("tensorfunctions.jl")
export contract_tensors

# Functionality to provide some basic spin Hamiltonians as MPOs
include("operators.jl")
export getPauliMatrices, getIsingMPO, getHeisenbergMPO, getTotalSpinMPO, getIdentityMPO

# Methods implementing the MPS algorithm with open boundary conditions
include("MPS_OBC.jl")
export Site, MPS, Operator, MPO, random_mps_obc, basis_state_obc, calculate_overlap, expectation_value, gaugeMPS, gaugeMPS!, apply_operator, apply_operator!, sum_states, sum_operators, contract_virtual_indices, find_groundstate, approximate_mps, compute_entropy, sample_from_mps!, sample_from_mps, svd_compress_mps, decompose_into_mpo

end # module


text:
using KNITRO
using Test

@testset "License manager test" begin
    m = KNITRO.Model()
    KNITRO.KN_free(m)
    # create license manager context
    @show lm = KNITRO.LMcontext()

    # we create 10 KNITRO instances with the license manager
    n_instances = 10
    kcs = [KNITRO.KN_new_lm(lm) for i in 1:n_instances]

    # then we free the instances
    for kc in kcs
        KNITRO.KN_free(kc)
    end

    # and release the license
    KNITRO.KN_release_license(lm)
    @test lm.ptr_lmcontext == C_NULL
end

@testset "User params in newpoint callback with license manager (issue #151)" begin
    function callbackNewPoint(kc, x, duals, userParams)
        # Test that inputs are of valid type.
        @test isa(kc, KNITRO.Model)
        @test isa(x, Vector{Cdouble})
        @test isa(duals, Vector{Cdouble})
        return 0
    end
    # Instantiate license manager context.
    lm = KNITRO.LMcontext()
    # Test different values for userparams: nothing is the default value
    # so allows to test default Knitro behavior. "myuserparam" tries
    # setting a string userparam in Knitro.
    for userparam in [nothing, "myuserparam"]
        kc = KNITRO.KN_new_lm(lm)
        KNITRO.KN_set_param(kc, "outlev", 0)
        # At first, `newpoint_user` is nothing
        @test isnothing(kc.newpoint_user)
        KNITRO.KN_add_vars(kc, 1)
        KNITRO.KN_set_var_lobnd(kc, Cint(0), 2.0)
        KNITRO.KN_set_var_primal_init_values_all(kc, [0.0])
        KNITRO.KN_set_obj_goal(kc, KNITRO.KN_OBJGOAL_MINIMIZE)
        KNITRO.KN_add_obj_quadratic_struct(kc, Cint[0], Cint[0], [1.0])
        KNITRO.KN_set_newpt_callback(kc, callbackNewPoint, userparam)
        # Test that userparam is correctly set
        @test kc.newpoint_user == userparam
        nstatus = KNITRO.KN_solve(kc)
        @test nstatus == KNITRO.KN_RC_OPTIMAL
        KNITRO.KN_free(kc)
    end
    KNITRO.KN_release_license(lm)
end


text:
using Test
using GeneralTools

println("Running tests:")

(testfiles, allnames) = targetlist(r"^(?!runtests).*(\.jl)$", "./");
# pwd() here should be "./test/"
for f in testfiles
	include(f);
end


text:
################################################################################
#
#  Base field
#
################################################################################

@doc Markdown.doc"""
    base_field(L::NumField) -> NumField

Given a number field $L/K$ this function returns the base field $K$.
For absolute extensions this returns $\mathbf{Q}$.
"""
base_field(::NumField)

_base_ring(K::NumField) = base_field(K)

_base_ring(::FlintRationalField) = FlintQQ

################################################################################
#
#  Predicates
#
################################################################################

export isabsolute

@doc Markdown.doc"""
    isabsolute(L::NumField) -> Bool

Returns whether $L$ is an absolute extension, that is, whether the base field
of $L$ is $\mathbf{Q}$.
"""
isabsolute(::NumField)

isabsolute(::NumField) = false

isabsolute(::NumField{fmpq}) = true

@doc Markdown.doc"""
    elem_type(L::NumField) -> Type

Returns the type of the elements of $L$.
"""
elem_type(::NumField)

################################################################################
#
#  Degree
#
################################################################################

@doc Markdown.doc"""
    degree(L::NumField) -> Int

Given a number field $L/K$, this function returns the degree of $L$ over $K$.
"""
degree(A::NumField)

dim(K::NumField) = degree(K)

################################################################################
#
#  Absolute degree
#
################################################################################

@doc Markdown.doc"""
    absolute_degree(L::NumField) -> Int

Given a number field $L/K$, this function returns the degree of $L$ over
$\mathbf Q$.
"""
function absolute_degree(A::NumField)
  return absolute_degree(base_field(A)) * degree(A)
end

################################################################################
#
#  Is simple extension
#
################################################################################

@doc Markdown.doc"""
    issimple(L::NumField) -> Bool

Given a number field $L/K$ this function returns whether $L$ is simple, that is,
whether $L/K$ is defined by a univariate polynomial$.
"""
issimple(a::NumField)

################################################################################
#
#  Number field creation
#
################################################################################

@doc Markdown.doc"""
    NumberField(f::Poly{NumFieldElem}, s::String;
                cached::Bool = false, check::Bool = false) -> NumField, NumFieldElem

Given an irreducible polynomial $f \in K[x]$ over some number field $K$, this
function creates the simple number field $L = K[x]/(x)$ and returns $(L, b)$,
where $b$ is the class of $x$ in $L$. The string `s` is used only for printing
the primitive element $b$.

Testing that $f$ is irreducible can be disabled by setting the keyword argument
`check` to `false`.
"""
NumberField(f::PolyElem{<:NumFieldElem}, s::String;
            cached::Bool = false, check::Bool = false)

################################################################################
#
#  Is commutative
#
################################################################################

iscommutative(K::NumField) = true

################################################################################
#
#  Pure and Kummer extensions
#
################################################################################

@doc Markdown.doc"""
    isradical_extension(L::SimpleNumField) -> Bool

Tests if $L/K$ is pure, that is, if the defining polynomial is of the form
$x^n - b$ for some $b \in K$.
"""
function isradical_extension(K::SimpleNumField)
  if !ismonic(K.pol)
    return false
  end
  return all(i -> iszero(coeff(K.pol, i)), 1:degree(K)-1)
end

@doc Markdown.doc"""
    iskummer_extension(L::SimpleNumField) -> Bool

Tests if $L/K$ is a Kummer extension, that is, if the defining polynomial is
of the form $x^n - b$ for some $b \in K$ and if $K$ contains the $n$-th roots
of unity.
"""
function iskummer_extension(K::SimpleNumField)
  if !isradical_extension(K)
    return false
  end

  k = base_field(K)
  Zk = maximal_order(k)
  _, o = torsion_units_gen_order(Zk)
  if o % degree(K) != 0
    return false
  end
  return true
end

function iskummer_extension(K::AnticNumberField)
  if degree(K) != 2
    return false
  end
  return isradical_extension(K)
end

function radical_extension(n::Int, a::FacElem;
                        cached::Bool = true, check::Bool = true)
  return radical_extension(n, evaluate(a), cached = cached, check = check)
end

@doc Markdown.doc"""
    radical_extension(n::Int, a::NumFieldElem; s = "_$",
                   check = true, cached = true) -> NumField, NumFieldElem

Given an element $a$ of a number field $K$ and an integer $n$, create the simple
extension of $K$ with the defining polynomial $x^n - a$.
"""
function radical_extension(n::Int, a::NumFieldElem; s::String = "_\$",
                        cached::Bool = true, check::Bool = true)
  k = parent(a)
  kx, x = PolynomialRing(k, cached = false)
  return number_field(x^n - a, s, check = check, cached = cached)
end

function radical_extension(n::Int, a::fmpq; s::String = "_\$",
                        cached::Bool = true, check::Bool = true)
  k = parent(a)
  kx, x = PolynomialRing(k, cached = false)
  return number_field(x^n - a, s, check = check, cached = cached)
end

## Non-simple

@doc Markdown.doc"""
    NumberField(f::Vector{PolyElem{<:NumFieldElem}}, s::String="_\$")
                                              -> NumField, Vector{NumFieldElem}

Given a list $f$ of univariate polynomials $f_1, \ldots, f_n \in K[x]$ over
some number field $K$, constructs the extension $K[x_1, \ldots, x_n]/(f_1(x_1),
\ldots, f_n(x_n))$.

The extensions $K[x]/(f_i)$ must be linearly disjoint or equivalently the ideal
$(f_1(x_1),\dotsc,f_n(x_n))$ must be maximal (although this is not tested).
"""
NumberField(::Vector{PolyElem{<:Union{NumFieldElem, fmpq}}}, ::String)

## Missing

@doc Markdown.doc"""
    basis(L::SimpleNumField) -> Vector{NumFieldElem}

Returns the canonical basis of a simple extension $L/K$, that is, the elements
$1,a,\dotsc,a^{d - 1}$, where $d$ is the degree of $K$ and $a$ the primitive
element.
"""
basis(::SimpleNumField)

@doc Markdown.doc"""
    basis(L::NonSimpleNumField) -> Vector{NumFieldElem}

Returns the canonical basis of a non-simple extension $L/K$. If $L = K(a_1,\dotsc,a_n)$
where each $a_i$ has degree $d_i$, then the basis will be $a_1^{i_1}\dotsm a_d^{i_d}$
with $0 \leq i_j d_j - 1$ for $1 \leq j \leq n$.
"""
basis(::NonSimpleNumField)

@doc Markdown.doc"""
    simple_extension(L::NonSimpleNumField) -> SimpleNumField, Map

Given a non-simple extension $L/K$, this function computes a simple extension $M/K$
and an $K$-linear isomorphism $M \to L$.
"""
simple_extension(::NonSimpleNumField)

@doc Markdown.doc"""
    absolute_field(L::NumField) -> NumField, Map

Given a number field $L$, this function returns an absolute simple number field $M/\mathbf{Q}$
together with a $\mathbf{Q}$-linear isomorphism $M \to K$.
"""
absolute_field(::NumField)

################################################################################
#
#  Defining polynomial
#
################################################################################

export defining_polynomial

@doc Markdown.doc"""
    defining_polynomial(L::SimpleNumField) -> PolyElem

Given a simple number field $L/K$, constructed as $L = K[x]/(f)$, this function
returns $f$.
"""
defining_polynomial(::SimpleNumField)

defining_polynomial(K::NfRel) = K.pol

defining_polynomial(K::AnticNumberField) = K.pol

@doc Markdown.doc"""
    discriminant(L::SimpleNumField) -> NumFieldElem

The discriminant of the defining polynomial of $L$, *not* the discriminant of
the maximal order of $L$.
"""
function discriminant(K::SimpleNumField)
  return discriminant(defining_polynomial(K))
end

@doc Markdown.doc"""
    absolute_discriminant(L::SimpleNumField, QQ) -> fmpq

The absolute discriminant of the defining polynomial of $L$, *not* the
discriminant of the maximal order of $L$. This is the norm of the discriminant
times the $d$-th power of the discriminant of the base field, where $d$ is the
degree of $L$.
"""
absolute_discriminant(::SimpleNumField)

function absolute_discriminant(K::AnticNumberField)
  return discriminant(K)
end

function absolute_discriminant(K::NfRel)
  d = norm(discriminant(K)) * absolute_discriminant(base_field(K))^degree(K)
  return d
end

function discriminant(K::FlintRationalField)
  return one(K)
end

##

@doc Markdown.doc"""
    issubfield(K::SimpleNumField, L::SimpleNumField) -> Bool, Map

Returns `true` and an injection from $K$ to $L$ if $K$ is a subfield of $L$.
Otherwise the function returns `false` and a morphism mapping everything to
$0$.
"""
issubfield(::SimpleNumField, ::SimpleNumField)

@doc Markdown.doc"""
    isisomorphic(K::SimpleNumField, L::SimpleNumField) -> Bool, Map

Returns `true` and an isomorphism from $K$ to $L$ if $K$ and $L$ are isomorphic.
Otherwise the function returns `false` and a morphism mapping everything to $0$.
"""
isisomorphic(::SimpleNumField, ::SimpleNumField)

################################################################################
#
#  Linear disjointness
#
################################################################################

# TODO (easy): Do this for Non-Simple number fields
@doc Markdown.doc"""
    islinear_disjoint(K::SimpleNumField, L::SimpleNumField) -> Bool

Given two number fields $K$ and $L$ with the same base field $k$, this function
returns whether $K$ and $L$ are linear disjoint over $k$.
"""
islinear_disjoint(K::SimpleNumField, L::SimpleNumField)

################################################################################
#
#  Normal basis
#
################################################################################

@doc Markdown.doc"""
    normal_basis(L::NumField) -> NumFieldElem

Given a normal number field $L/K$, this function returns an element $a$ of $L$,
such that the orbit of $a$ under the Galois group of $L/K$ is an $K$-basis
of $L$.
"""
function normal_basis(L::NumField)
  n = degree(L)
  K = base_field(L)
  Aut = automorphisms(L)

  length(Aut) != degree(L) && error("The field is not normal over the rationals!")

  A = zero_matrix(K, n, n)
  r = one(L)
  while true
    r = rand(basis(L), -n:n)
    for i = 1:n
      y = Aut[i](r)
      for j = 1:n
        A[i,j] = coeff(y, j - 1)
      end
    end
    if !iszero(det(A))
      break
    end
  end
  return r
end

################################################################################
#
#  Check linearly disjointness
#
################################################################################

function check_consistency(K::NonSimpleNumField)
  QQz, z = PolynomialRing(base_field(K), "z")
  for i = 1:length(K.pol)
    v = [zero(QQz) for i in 1:length(K.pol)]
    v[i] = z
    p = evaluate(K.pol[i], v)
    if !isirreducible(p)
      return false
    end
  end
  lg = gens(K)
  el = lg[1]
  f = minpoly(el)
  d = degree(f)
  for i = 2:length(lg)
    el += lg[i]
    f = minpoly(el)
    while !issquarefree(f)
      el += lg[i]
      f = minpoly(el)
    end
    if degree(f) != prod(degree(K.pol[j], j) for j = 1:i)
      return false
    end
    if !isirreducible(f)
      return false
    end
  end
  return true
end


text:
import FractionalCalculus.FracDiffAlg

"""
Riemann-Liouville sense fractional derivative algorithms, please refer to [Riemann-Liouville derivative](https://en.wikipedia.org/wiki/Fractional_calculus#Riemann%E2%80%93Liouville_fractional_derivative)
"""
abstract type RLDiff <: FracDiffAlg end

#RLDiff_LinearInterp maybe??
"""
# Riemann Liouville sense derivative approximation

    fracdiff(f, α, end_point, h, RLDiff_Approx())

Using Linear interpolation to approximate fractional derivative in Riemann Liouville  fractional derivative sense.

### Example

```julia-repl
julia> fracdiff(x->x^5, 0.5, 2.5, 0.0001, RLDiff_Approx())
```

!!! warning
    The RLDiff_Approx algorithm only support for 0 < α < 1.
"""
struct RLDiff_Approx <: RLDiff end

"""
# Riemann Liouville sense derivative using Triangular Strip Matrix to discrete and compute.

    fracdiff(f, α, end_point, h, RLDiff_Matrix())

Using [Triangular Strip Matrix](https://en.wikipedia.org/wiki/Triangle_strip) to approximate fractional derivative.

### Example

```julia-repl
julia> fracdiff(x->x^5, 0.5, 2.5, 0.0001, RLInt_Matrix())
```

!!! info
    Triangular Strip Matrix method returns the derivative in the interval ``[0, T]`` in ```Vector```

!!! tip
    With the advancing Triangular Strip Matrix method, you can not only compute fractional derivatives, integer order, higher order derivative is also supported!!
Try to set α as an integer, arbitrary integer of course! I promise you would enjoy it😏

```tex
@article{2009,
title={Matrix approach to discrete fractional calculus II: Partial fractional differential equations},
DOI={10.1016/j.jcp.2009.01.014},
author={Podlubny, Igor and Chechkin, Aleksei and Skovranek, Tomas and Chen, YangQuan and Vinagre Jara, Blas M.},
}
```
"""
struct RLDiff_Matrix <: RLDiff end

"""
Using linear spline interpolation method to approximate the Riemann Liouville fractional derivative.
"""
struct RL_Linear_Spline_Interp <: RLDiff end

################################################################
###                    Type defination done                  ###
################################################################

function fracdiff(f::Union{Number, Function}, α, end_point, h, ::RLDiff_Approx)
    #checks(f, α, 0, end_point)

    summation = 0
    n = Int64(floor(end_point/h))

    @fastmath @inbounds @simd for i ∈ 0:n-1
        summation += (f(end_point-i*h) - f(end_point-(i+1)*h))*((i+1)^(1-α) - i^(1-α))
    end

    result = ((1-α)*f(0)/n^α+summation)*end_point^(-α)*n^α/gamma(2-α)
    return result
end

function fracdiff(f::Union{Number, Function}, α::Float64, end_point::AbstractArray, h, ::RLDiff_Approx)::Vector
    result = map(x->fracdiff(f, α, x, h, RLDiff_Approx()), end_point)
    return result
end

function fracdiff(f, α, end_point, h::Float64, ::RLDiff_Matrix)
    N = Int64(end_point/h+1)
    @views tspan = collect(0:h:end_point)
    return B(N, α, h)*f.(tspan)
end

#Compute the eliminator matrix Sₖ by omiting n-th row
function eliminator(n, row)
    temp = zeros(n, n) + I
    return @views temp[Not(row), :]
end

#Generate elements in Matrix.
function omega(n, p)
    omega = zeros(n+1)

    omega[1]=1
    @fastmath @inbounds @simd for i ∈ 1:n
        omega[i+1]=(1-(p+1)/i)*omega[i]
    end

    return omega
end
function B(N, p, h::Float64)
    result = zeros(N, N)
    temp = omega(N, p)

    @inbounds @simd for i ∈ 1:N
        @views result[i, 1:i]=reverse(temp[1:i])
    end

    return h^(-p)*result
end
# Multiple dispatch for not assigning step size *h*.
function B(N, p)
    result=zeros(N, N)
    temp=omega(N, p)

    @inbounds @simd for i ∈ 1:N
        @views result[i, 1:i]=reverse(temp[1:i])
    end

    return result
end

#=
Numerical methods for fractional calculus by Li, Changpin
Page 57

Linear Spline Interpolation
=#
function fracdiff(f, α, x, h, ::RL_Linear_Spline_Interp)
    N = Int64(floor(x/h))

    result = 0

    for k = 0:(N+1)
        result += z̄ₘₖ(N, k, α)*f(k*h)
    end

    return 1/(gamma(4-α)h^α)*result

end

function z̄ₘₖ(m, k, α)
    if k ≤ m-1
        return c̄ⱼₖ(m-1, k, α)-2*c̄ⱼₖ(m, k, α)+c̄ⱼₖ(m+1, k, α)
    elseif k == m
        return -2*c̄ⱼₖ(m, k, α)+c̄ⱼₖ(m+1, k, α)
    elseif k == m+1
        return c̄ⱼₖ(m+1, k, α)
    elseif k > m+1
        return 0
    end
end

function c̄ⱼₖ(j, k, α)
    if k==0
        return (j-1)^(3-α)-j^(2-α)*(j-3+α)
    elseif 1 ≤ k ≤ j-1
        return (j-k+1)^(3-α)-2*(j-k)^(3-α)+(j-k-1)^(3-α)
    elseif k == j
        return 1
    end
end

function fracdiff(f::Union{Number, Function}, α::Float64, end_point::AbstractArray, h, ::RL_Linear_Spline_Interp)::Vector
    result = map(x->fracdiff(f, α, x, h, RL_Linear_Spline_Interp()), end_point)
    return result
end


text:
# This file is a part of SolidStateDetectors.jl, licensed under the MIT License (MIT).
"""
    struct RectangularCuboid{T} <: AbstractGeometry{T, 3}

...
"""
struct RectangularCuboid{T} <: AbstractGeometry{T, 3}
    org::CartesianPoint{T}
    e1::CartesianVector{T}
    e2::CartesianVector{T}
    e3::CartesianVector{T}
end

function in(pt::CartesianPoint{T}, rc::RectangularCuboid{T})::Bool where {T <: Real}
    shift::CartesianVector{T} = pt - rc.org
    return  shift ⋅ rc.e1 >= 0 &&
            shift ⋅ rc.e2 >= 0 &&
            shift ⋅ rc.e3 >= 0 &&
            (shift - rc.e1) ⋅ rc.e1 <= 0 &&
            (shift - rc.e2) ⋅ rc.e2 <= 0 &&
            (shift - rc.e3) ⋅ rc.e3 <= 0
end

function vertices(rc::RectangularCuboid{T})::Vector{CartesianPoint{T}} where {T}
    v::Vector{CartesianPoint{T}} = CartesianPoint{T}[
        rc.org,
        rc.org + rc.e1,
        rc.org + rc.e1 + rc.e2,
        rc.org + rc.e2,
        rc.org + rc.e3,
        rc.org + rc.e3 + rc.e1,
        rc.org + rc.e3 + rc.e1 + rc.e2,
        rc.org + rc.e3 + rc.e2,
    ]
    return v
end

function LineSegments(rc::RectangularCuboid{T})::Vector{LineSegment{T, 3, :cartesian}} where {T}
    return LineSegment{T, 3, :cartesian}[
        LineSegment{T, 3, :cartesian}( rc.org, rc.e1),
        LineSegment{T, 3, :cartesian}( rc.org + rc.e1,  rc.e2),
        LineSegment{T, 3, :cartesian}( rc.org + rc.e1 + rc.e2, -rc.e1),
        LineSegment{T, 3, :cartesian}( rc.org + rc.e2, -rc.e2),
        LineSegment{T, 3, :cartesian}( rc.org, rc.e3),
        LineSegment{T, 3, :cartesian}( rc.org + rc.e1,  rc.e3),
        LineSegment{T, 3, :cartesian}( rc.org + rc.e1 + rc.e2, rc.e3),
        LineSegment{T, 3, :cartesian}( rc.org + rc.e2, rc.e3),
        LineSegment{T, 3, :cartesian}( rc.org + rc.e3, rc.e1),
        LineSegment{T, 3, :cartesian}( rc.org + rc.e3 + rc.e1,  rc.e2),
        LineSegment{T, 3, :cartesian}( rc.org + rc.e3 + rc.e1 + rc.e2, - rc.e1),
        LineSegment{T, 3, :cartesian}( rc.org + rc.e3 + rc.e2, -rc.e2)
    ]
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule nlminb_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("nlminb")
JLLWrappers.@generate_main_file("nlminb", UUID("c04a2911-6870-53d6-aa47-4d7d5aff1028"))
end  # module nlminb_jll


text:
function _update_till_convergence!( pssrb::PotentialSimulationSetupRB{T, S, 3},
                                    convergence_limit::T,
                                    device_array_type::Type{DAT};
                                    n_iterations_between_checks = 500,
                                    depletion_handling::Val{depletion_handling_enabled} = Val{false}(),
                                    only2d::Val{only_2d} = Val{false}(),
                                    is_weighting_potential::Val{_is_weighting_potential} = Val{false}(),
                                    use_nthreads::Int = Base.Threads.nthreads(),
                                    max_n_iterations::Int = 10_000, # -1
                                    verbose::Bool = true
                                )::T where {T, S, depletion_handling_enabled, only_2d, _is_weighting_potential, DAT <: GPUArrays.AbstractGPUArray}
    device = get_device(DAT)
    N_grid_points = prod(size(pssrb.potential)[1:3] .- 2)
    kernel = get_sor_kernel(S, device)
    @showprogress for i in 1:max_n_iterations
        update_even_points = true
        wait(kernel(
            pssrb.potential, pssrb.point_types, pssrb.volume_weights, pssrb.q_eff_imp, pssrb.q_eff_fix, pssrb.ϵ_r,
            pssrb.geom_weights, pssrb.sor_const, update_even_points, depletion_handling_enabled, _is_weighting_potential, only_2d,
            ndrange = N_grid_points
        ))
        apply_boundary_conditions!(pssrb, Val(update_even_points), only2d)
        update_even_points = false
        wait(kernel(
            pssrb.potential, pssrb.point_types, pssrb.volume_weights, pssrb.q_eff_imp, pssrb.q_eff_fix, pssrb.ϵ_r,
            pssrb.geom_weights, pssrb.sor_const, update_even_points, depletion_handling_enabled, _is_weighting_potential, only_2d,
            ndrange = N_grid_points
        ))
        apply_boundary_conditions!(pssrb, Val(update_even_points), only2d)
    end
    return 0
end

get_sor_kernel(::Type{Cylindrical}, args...) = sor_cyl_gpu!(args...)
get_sor_kernel(::Type{Cartesian},   args...) = sor_car_gpu!(args...)

function get_device end


text:
using DashBootstrapComponents

nav = dbc_nav([
    dbc_navlink("Active", active = true, href = "#"),
    dbc_navlink("A link", href = "#"),
    dbc_navlink("Another link", href = "#"),
    dbc_navlink("Disabled", disabled = true, href = "#"),
]);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
Aij_indices::Vector{Tuple{Int64,Int64,Int64}}

indice1 = [1,1,1,1,1,1,1,2,2,4,4,4,5,5,6]
indice2 = [2,3,4,5,6,7,15,3,4,5,6,7,6,7,7]
indice3 = 1:length(indice1)
unique_indices = intersect([1,2,3,4,5,6,15],xsel)

Aij_indices =[(a1,a2,a3) for (a1,a2,a3) in
zip(indice1[unique_indices],indice2[unique_indices],indice3[unique_indices])]

for kk in 1:length(model.Aij_indices) # i are CartesianIndices

    i1,i2,i0 = model.Aij_indices[kk]
#    i0 = model.Aij_indices[kk]
 #   i1 = model.Aij_indices[kk]
  #  i2 = model.Aij_indices[kk]


text:
#=
I don't know if I should rescale by system size / unit cell
=#
function fidelity_susceptibility(state::Union{FiniteMPS,InfiniteMPS},H₀::T,Vs::AbstractVector{T},hpars = params(state,H₀);maxiter=Defaults.maxiter,tol=Defaults.tol) where T<:MPOHamiltonian
    init_v = rand_quasiparticle(state,excitation_space = oneunit(space(state.AC[1],1)));

    tangent_vecs = map(Vs) do V
        vpars = params(state,V)

        Tos = similar(init_v)
        for (i,ac) in enumerate(state.AC)
            temp = ac_prime(ac,i,state,vpars);
            help = Tensor(ones,oneunit(space(state.AC[1],1)))
            @tensor tor[-1 -2;-3 -4]:= temp[-1,-2,-4]*help[-3]
            Tos[i] = tor
        end

        (vec,convhist) = linsolve(Tos,Tos,GMRES(maxiter=maxiter,tol=tol)) do x
            effective_excitation_hamiltonian(H₀, x, params(x,H₀,hpars))
        end
        convhist.converged == 0 && @info "failed to converge $(convhist.normres)"

        vec
    end

    map(Iterators.product(tangent_vecs,tangent_vecs)) do (v,w)
        dot(v,w)
    end
end