text:
import AtomicLevels: spectroscopic_label

# * Plain hydrogenics

"""
    hydrogenic!(atom[; find_lowest=false, find_lowest_ℓmax=Inf, kwargs...])

Initialize the radial orbitals of `atom` to their unscreened
hydrogenic values. This is done via simple diagonalization of the
one-body Hamiltonian for each angular symmetry. If `find_lowest` is
`true`, only the orbital(s) with the lowest energy is kept (out of
those with `ℓ≤find_lowest_ℓmax`). The `kwargs` are passed on to
[`diagonalize_one_body`](@ref) and can be used to influence how the
diagonalization is performed.
"""
function hydrogenic!(atom::Atom{T,B,O,TC,C,P}; verbosity=0, find_lowest::Bool=false, find_lowest_ℓmax=Inf, kwargs...) where {T,B,O,TC,C,P}
    # This is only for non-relativistic atoms

    print_block() do io
        verbosity > 0 && println(io, "Hydrogenic initialization of the orbitals of $(atom)")

        R,Φ = atom.radial_orbitals.args
        Φ .= zero(T)

        # Orbitals are grouped by ℓ, since they have the same one-body
        # Hamiltonian. They should rather be grouped by
        # `symmetry(orb)`, where `symmetry` dispatches on the type of
        # the orbital. This would allow for pseudo-potentials that are
        # ℓ- and s-dependent.
        orbitals = Dict{Int,Vector{O}}()
        for orb in atom.orbitals
            orbitals[getℓ(orb)] = vcat(get(orbitals, getℓ(orb), O[]), orb)
        end

        Z = charge(atom.potential)
        Q = num_electrons(outsidecoremodel(first(atom.configurations),
                                           atom.potential))
        # This is a heuristic to get energy shifts that work for both
        # point charge nuclei and pseudopotentials; the latter screen
        # the bare charge but not fully, and thus we admix with a
        # little bit of the unscreened charge.
        hydrogen_like_eng = n -> -(0.3Z+0.7Q)^2/(2n^2)

        eng_fmt = FormatSpec("+10.7f")
        Δeng_fmt = FormatSpec("+10.3e")
        fmt_eng_vec(fspec, v) = "["*join(fmt.(Ref(fspec), v), ", ")*"]"

        λmin = Inf

        for ℓ in keys(orbitals)
            ℓ > find_lowest_ℓmax && continue
            print_block(io) do io
                min_n,max_n,nev = if any(!isbound, orbitals[ℓ])
                    ℓ+1,0,1
                else
                    min_n = minimum(getn, orbitals[ℓ])
                    max_n = maximum(getn, orbitals[ℓ])

                    nev = max_n - (min_n-1)
                    nev ≥ 1 || throw(BoundsError("Trying to find $(nev) eigenvalues"))
                    min_n,max_n,nev
                end

                # For shift-and-invert type of diagonalization, we target an
                # eigenvalue slightly below the lowest energy in partial wave
                # ℓ of a hydrogen-like system with point charge Z.
                Iₚℓ = hydrogen_like_eng(min_n)
                σ = 2Iₚℓ

                verbosity > 1 &&
                    println(io, "Diagonalizing symmetry ℓ = $(spectroscopic_label(ℓ)), maximum n = $(max_n) => $(nev) eigenvalues required")
                verbosity > 2 &&
                    println(io, "Target eigenvalue: ≤ $(Iₚℓ) Ha")

                H = one_body_hamiltonian(atom, first(orbitals[ℓ]))
                λᴴ,Φᴴ = diagonalize_one_body(H, nev; σ=σ, io=io, verbosity=verbosity, kwargs...)
                if find_lowest
                    if first(λᴴ) < λmin
                        λmin = first(λᴴ)
                        Φ .= zero(T)
                        for orb in orbitals[ℓ]
                            j = findfirst(isequal(orb), atom.orbitals)
                            copyto!(view(Φ, :, j), view(Φᴴ, :, 1))
                        end
                    end
                else
                    for orb in orbitals[ℓ]
                        j = findfirst(isequal(orb), atom.orbitals)
                        copyto!(view(Φ, :, j), view(Φᴴ, :, nev > 1 ? getn(orb)-(min_n-1) : 1))
                    end
                end

                if verbosity > 2 && max_n > 0
                    λᴴₐ = hydrogen_like_eng.(min_n:max_n)
                    println(io, "Hydrogenic energies $(fmt_eng_vec(eng_fmt, λᴴ)) Ha")
                    println(io, "Analytic energies   $(fmt_eng_vec(eng_fmt, λᴴₐ)) Ha")
                    println(io, "Δ                   $(fmt_eng_vec(Δeng_fmt, λᴴ-λᴴₐ)) Ha")
                end
            end
        end

        print_block(io) do io
            ml = maximum(length.(string.(atom.orbitals)))
            linefmt = FormatExpr("{1:$(ml)s} {2:<9.7f} {3:12.5e}")
            verbosity > 2 && printfmtln(io, "{1:$(ml)s} {2:9s}  {3:11s}", "", "Initial n", "1-n")
            for j in eachindex(atom.orbitals)
                n₀ = norm(atom[j])
                iszero(n₀) || norm_rot!(view(atom, j))

                if verbosity > 2
                    n = norm(atom[j])
                    printfmtln(io, linefmt, atom.orbitals[j], n₀, 1-n)
                end
            end
        end
        if verbosity > 3
            print_block(io) do io
                nconfigs = min(10, length(atom.configurations))
                ml = maximum(length.(string.(atom.configurations)))
                configfmt = "{1:<$(ml+3)s}"
                linefmt = FormatExpr("$(configfmt) {2:7.5f} {3:12.5e}")
                N = Q
                printfmtln(io, "$(configfmt) {2:7s}  {3:11s} ", "Cfg", "√N", "$(N)-N²")
                for (i,config) in enumerate(atom.configurations)
                    i > nconfigs && break
                    n = norm(atom, configuration=i)
                    printfmtln(io, linefmt, config, n, N-n^2)
                end
                length(atom.configurations) > nconfigs && println(io, "⋮")
            end
        end
    end

    atom
end

# * Screened hydrogenics

"""
    screened_hydrogenic!(atom[; kwargs...])

Initialize the radial orbitals of `atom` to their screened hydrogenic
values. This is done via simple diagonalization of the one-body
Hamiltonian for each orbital with screening computed from all the
other orbitals of the first configuration of `atom`. The `kwargs` are
passed on to [`diagonalize_one_body`](@ref) and can be used to
influence how the diagonalization is performed.
"""
function screened_hydrogenic!(atom::Atom{T,B,O,TC,C,P}; verbosity=0, kwargs...) where {T,B,O,TC,C,P}
    # This is only for non-relativistic atoms

    print_block() do io
        verbosity > 0 && println(io, "Hydrogenic initialization of the orbitals of $(atom)")

        R,Φ = atom.radial_orbitals.args
        Φ .= zero(T)

        cfg = first(atom.configurations)

        Z = charge(atom.potential)
        Q = num_electrons(outsidecoremodel(cfg, atom.potential))
        # This is a heuristic to get energy shifts that work for both
        # point charge nuclei and pseudopotentials; the latter screen
        # the bare charge but not fully, and thus we admix with a
        # little bit of the unscreened charge.
        hydrogen_like_eng = σ -> -(0.3Z+0.7Q-σ)^2/2

        eng_fmt = FormatSpec("+10.7f")
        Δeng_fmt = FormatSpec("+10.3e")

        for (i,o) in enumerate(atom.orbitals)
            isbound(o) || continue
            print_block(io) do io
                # For shift-and-invert type of diagonalization, we
                # target an eigenvalue slightly below the lowest
                # energy of a hydrogen-like system with point charge
                # Z - sc.
                sc = screening(getspatialorb(o), cfg)
                Vsc = operator(potential_matrix(r -> sc/r, R), R)
                Iₚ = hydrogen_like_eng(0)
                σ = 2Iₚ
                nev = getn(o)-getℓ(o)

                verbosity > 1 &&
                    println(io, "Diagonalizing orbital $o with screening $sc")
                verbosity > 2 &&
                    println(io, "Target eigenvalue: ≤ $(Iₚ) Ha")

                H = one_body_hamiltonian(atom, o) + Vsc
                λᴴ,Φᴴ = diagonalize_one_body(H, nev; σ=σ, io=io, verbosity=verbosity, kwargs...)
                copyto!(view(Φ, :, i), view(Φᴴ, :, nev))

                if verbosity > 2
                    println(io, "Hydrogenic energy $(fmt(eng_fmt, λᴴ[1])) Ha")
                    println(io, "Target energy     $(fmt(eng_fmt, Iₚ)) Ha")
                    println(io, "Δ                 $(fmt(Δeng_fmt, λᴴ[1]-Iₚ)) Ha")
                end
            end
        end
    end

    atom
end

# ** Screening

@doc raw"""
    screening(i, j)

Compute the amount of screening of orbital `i` due to orbital `j`,
according to the formula

```math
\sigma_{ij} = \left\{
1 + \left[
\frac{3n_j^2 - \ell_j(\ell_j+1)}{3n_i^2 - \ell_i(\ell_i+1)}
\right]^2
\right\}^{-3/2}
```

taken from Eq. (10) of

    Bessis, N., & Bessis, G. (1981). Analytic Atomic Shielding
    Parameters. The Journal of Chemical Physics, 74(6),
    3628–3630. http://dx.doi.org/10.1063/1.441475
"""
function screening(i::Orbital, j::Orbital)
    r = (3j.n^2 - j.ℓ*(j.ℓ+1))/(3i.n^2 - i.ℓ*(i.ℓ+1))
    (1 + r^2)^(-3/2)
end

@doc raw"""
    screening(i, c)

Compute the screening of orbital `i` due to all _other_ orbitals of
configuration `c`:

```math
\sigma_i = \sum_j (w_j-\delta_{ij})\sigma_{ij}
```

where ``w_j`` is the occupancy of orbital `j`.

# Examples

```jldoctest
julia> Atoms.screening(o"1s", c"1s2 2s2")
0.3820869935387247
```

The `1s` orbital is only slightly screened by the other `1s` electron
and the 2 `2s` electrons, whereas

```jldoctest
julia> Atoms.screening(o"2s", c"1s2 2s2")
2.179703979102134
```

shows that the `2s` electron is screened by both the `1s` electrons
and a little bit of the the other `2s` electron.
"""
function screening(i::Orbital, c::Configuration)
    σᵢ = 0.0
    for (j,wⱼ,_) in c
        σᵢ += (wⱼ - (i==j))*screening(i, getspatialorb(j))
    end
    σᵢ
end

export hydrogenic!


text:
using Test, Qsosed

@testset "test parse parameters" begin
    config_path = String(@__DIR__) * "/../configs/config_example.yaml"
    parameters = Parameters(config_path)
    @test parameters.M == 1e8
    @test parameters.mdot == 0.5
    @test parameters.mdot == 0.5
    @test parameters.spin == 0.0
    @test parameters.mu_nucleon == 0.61
    @test parameters.reprocessing == true
    @test parameters.hard_xray_fraction == 0.02
    @test parameters.corona_electron_energy == 100
    @test parameters.warm_electron_energy == 0.2
    @test parameters.warm_photon_index == 2.5
    @test parameters.reflection_albedo == 0.3
    @test parameters.min_energy == 1e-4 # kev
    @test parameters.max_energy == 200
    @test parameters.uv_min_energy == 0.00387
    @test parameters.uv_max_energy == 0.06
    @test parameters.xray_min_energy == 0.06
    @test parameters.xray_max_energy == 1
end

ext: jl


text:
using MonteCarloMeasurements: Particles

@testset "InferenceData" begin
    data = load_arviz_data("centered_eight")

    @testset "construction" begin
        pydata = PyObject(data)
        @test InferenceData(pydata) isa InferenceData
        @test PyObject(InferenceData(pydata)) === pydata
        @test InferenceData(data) === data
        @test_throws ArgumentError InferenceData(py"PyNullObject()")
        data2 = InferenceData(; posterior=data.posterior)
        @test data2 isa InferenceData
        @test :posterior in propertynames(data2)
        @test hash(data) == hash(pydata)
    end

    @testset "properties" begin
        @test :posterior in propertynames(data)
        @test length(propertynames(data)) > 1
        data3 = InferenceData(; posterior=data.posterior, prior=data.prior)
        @test :prior in propertynames(data3)
        delete!(data3, :prior)
        @test :prior ∉ propertynames(data3)
    end

    @testset "groups" begin
        data4 = InferenceData(; posterior=data.posterior)
        @test ArviZ.groupnames(data4) == [:posterior]
        g = ArviZ.groups(data4)
        @test g isa Dict
        @test :posterior in keys(g)
    end

    @testset "isempty" begin
        @test !isempty(data)
        @test isempty(InferenceData())
    end

    @testset "conversion" begin
        @test pyisinstance(PyObject(data), ArviZ.arviz.InferenceData)
        data4 = convert(InferenceData, PyObject(data))
        @test data4 isa InferenceData
        @test PyObject(data4) === PyObject(data)

        # TODO: improve this test
        @test convert(InferenceData, [1.0, 2.0, 3.0, 4.0]) isa InferenceData
    end

    @testset "show" begin
        @testset "plain" begin
            text = sprint(show, data)
            @test startswith(text, "InferenceData with groups:")
            rest = split(PyObject(data).__str__(), '\n'; limit=2)[2]
            @test split(text, '\n'; limit=2)[2] == rest
        end

        @testset "html" begin
            text = repr(MIME("text/html"), data)
            @test text isa String
            @test !occursin("<xarray.Dataset>", text)
            @test !occursin("<xarray.Dataset>", text)
            @test !occursin("arviz.InferenceData", text)
            @test occursin("Dataset (xarray.Dataset)", text)
            @test occursin("InferenceData", text)
        end
    end
end

@testset "+(::InferenceData, ::InferenceData)" begin
    rng = MersenneTwister(42)
    idata1 = from_dict(;
        posterior=Dict("A" => randn(rng, 2, 10, 2), "B" => randn(rng, 2, 10, 5, 2))
    )
    idata2 = from_dict(;
        prior=Dict("C" => randn(rng, 2, 10, 2), "D" => randn(rng, 2, 10, 5, 2))
    )

    new_idata = idata1 + idata2
    @test new_idata isa InferenceData
    @test check_multiple_attrs(
        Dict(:posterior => ["A", "B"], :prior => ["C", "D"]), new_idata
    ) == []
end

@testset "concat" begin
    rng = MersenneTwister(42)
    idata1 = from_dict(;
        posterior=Dict("A" => randn(rng, 2, 10, 2), "B" => randn(rng, 2, 10, 5, 2))
    )
    idata2 = from_dict(;
        prior=Dict("C" => randn(rng, 2, 10, 2), "D" => randn(rng, 2, 10, 5, 2))
    )

    new_idata = concat(idata1, idata2)
    @test new_idata isa InferenceData
    @test check_multiple_attrs(
        Dict(:posterior => ["A", "B"], :prior => ["C", "D"]), new_idata
    ) == []

    new_idata = concat!(idata1, idata2)
    @test new_idata === idata1
    @test check_multiple_attrs(
        Dict(:posterior => ["A", "B"], :prior => ["C", "D"]), idata1
    ) == []

    idata3 = concat!()
    @test idata3 isa InferenceData
    @test isempty(idata3)

    idata4 = InferenceData()
    new_idata = concat!(idata4, idata1)
    @test new_idata === idata4
    @test check_multiple_attrs(
        Dict(:posterior => ["A", "B"], :prior => ["C", "D"]), idata4
    ) == []
end

@testset "convert_to_inference_data" begin
    rng = MersenneTwister(42)

    @testset "convert_to_inference_data(::Dict)" begin
        dataset = Dict("A" => randn(rng, 2, 10, 2), "B" => randn(rng, 2, 10, 5, 2))
        idata1 = convert_to_inference_data(dataset)
        @test idata1 isa InferenceData
        @test check_multiple_attrs(Dict(:posterior => ["A", "B"]), idata1) == []
    end

    @testset "convert_to_inference_data(::Array)" begin
        arr = randn(rng, 2, 10, 2)
        idata2 = convert_to_inference_data(arr)
        @test check_multiple_attrs(Dict(:posterior => ["x"]), idata2) == []
    end

    @testset "convert_to_inference_data(::Nothing)" begin
        idata3 = convert_to_inference_data(nothing)
        @test idata3 isa InferenceData
        @test isempty(idata3)
    end

    @testset "convert_to_inference_data(::Particles)" begin
        p = Particles(randn(rng, 10))
        idata4 = convert_to_inference_data(p)
        @test check_multiple_attrs(Dict(:posterior => ["x"]), idata4) == []
    end

    @testset "convert_to_inference_data(::Vector{Particles})" begin
        p = [Particles(randn(rng, 10)) for _ in 1:4]
        idata5 = convert_to_inference_data(p)
        @test check_multiple_attrs(Dict(:posterior => ["x"]), idata5) == []
    end
    @testset "convert_to_inference_data(::Vector{Array{Particles}})" begin
        p = [Particles(randn(rng, 10, 3)) for _ in 1:4]
        idata6 = convert_to_inference_data(p)
        @test check_multiple_attrs(Dict(:posterior => ["x"]), idata6) == []
    end
end

@testset "from_dict" begin
    rng = MersenneTwister(42)

    posterior = Dict("A" => randn(rng, 2, 10, 2), "B" => randn(rng, 2, 10, 5, 2))
    prior = Dict("C" => randn(rng, 2, 10, 2), "D" => randn(rng, 2, 10, 5, 2))

    idata = from_dict(posterior; prior=prior)
    @test check_multiple_attrs(
        Dict(:posterior => ["A", "B"], :prior => ["C", "D"]), idata
    ) == []

    idata2 = from_dict(; prior=prior)
    @test check_multiple_attrs(Dict(:prior => ["C", "D"]), idata2) == []
end

@testset "netcdf roundtrip" begin
    data = load_arviz_data("centered_eight")
    mktempdir() do path
        filename = joinpath(path, "tmp.nc")
        to_netcdf(data, filename)
        data2 = from_netcdf(filename)
        @test data2 isa InferenceData
        @test propertynames(data) == propertynames(data2)
        return nothing
    end
end

@testset "load_arviz_data" begin
    data = load_arviz_data("centered_eight")
    @test data isa InferenceData

    datasets = load_arviz_data()
    @test datasets isa Dict
end


text:
# ## `base.jl`
#
# First, we define the abstract types for each concept.

"""
    abstract type AbstractItem

Abstract supertype for all items. To implement items, subtype
either [`Item`](#) to create a new item or [`ItemWrapper`](#)
to wrap an existing item.
"""
abstract type AbstractItem end

"""
    abstract type Item

Abstract supertype of concrete items.

Subtype if you want to create a new item. If you want to wrap
an existing item, see [`ItemWrapper`](#).
"""
abstract type Item <: AbstractItem end

"""
    showitem(item)

Visualize `item`. Should return an image.
"""
function showitem end

function showitem(items::Tuple)
    return blend.(showitem.(items)...)
end

"""
    abstract type Transform

Abstract supertype for all transformations.
"""
abstract type Transform end

# Plus the default implementations of [`itemdata`](#) and
# [`getrandstate`](#).

"""
    itemdata(item)
    itemdata(items)

Access the data wrapped in `item` or a tuple of items.
"""
itemdata(item::Item) = item.data
itemdata(items::Tuple) = itemdata.(items)

"""
    getrandstate(transform)

Generates random state for stochastic transformations.
Calling `apply(tfm, item)` is equivalent to
`apply(tfm, item; randstate = getrandstate(tfm))`. It
defaults to `nothing`, so you it only needs to be implemented
for stochastic `Transform`s.
"""
getrandstate(::Transform) = nothing

# The following 2 methods are convenience for `apply` so that
# you don't have to explicitly handle the random state. Passing
# a tuple of `Item`s to `apply` will also use the same random
# state for each item.

"""
    apply(tfm, item[; randstate])
    apply(tfm, items[; randstate])

Apply `tfm` to an `item` or a tuple `items`.

"""
apply(tfm::Transform, items) = apply(tfm, items; randstate = getrandstate(tfm))

function apply(tfm::Transform, items::Tuple; randstate = getrandstate(tfm))
    map(item -> apply(tfm, item; randstate = randstate), items)
end

#

"""
    compose(transforms...)

Compose tranformations. Use `|>` as an alias.

Defaults to creating a [`Sequence`](#) of transformations,
but smarter behavior can be implemented.
For example, `MapElem(f) |> MapElem(g) == MapElem(g ∘ f)`.
"""
function compose end
compose(tfm) = tfm
compose(tfms...) = compose(compose(tfms[1], tfms[2]), tfms[3:end]...)
Base.:(|>)(tfm1::Transform, tfm2::Transform) = compose(tfm1, tfm2)

# [`Identity`](#) is the identity transformation.

"""
    Identity()

The identity transformation.
"""
struct Identity <: Transform end

apply(::Identity, item::Item; randstate = nothing) = item

compose(::Identity, ::Identity) = Identity()
compose(tfm::Transform, ::Identity) = tfm
compose(::Identity, tfm::Transform) = tfm

# Lastly, [`setdata`](#) provides a convenient way to create a copy
# of an item, replacing only the wrapped data. This relies on the
# wrapped data field being named `data`, though.

function setdata(item::Item, data)
    item = Setfield.@set item.data = data
    return item
end


text:
## Splitting Strategies
## ====================

# This splits along all dimensions, ignores the depth
function weighted_mid_split(o::Omega{Float64}, depth::Int)
  splitted = mid_split(o::Omega{Float64})
  transitionprob = 1/length(splitted)
  [(event, transitionprob) for event in splitted]
end

# This splits along the 1st dimension in for the first split
# The second dimension for the second split, then cycles
function weighted_partial_split(o::Omega{Float64}, depth::Int)
  dimindices = collect(keys(o.intervals))
  ndims = length(dimindices)

  if (depth % length(dimindices)) + 1 > length(dimindices)
    @show depth
    @show dimindices
  end

  dimtosplit = dimindices[(depth % length(dimindices)) + 1]
  splitted = mid_partial_split(o::Omega{Float64}, [dimtosplit])
  transitionprob = 1/length(splitted)
  [(event, transitionprob) for event in splitted]
end

# Randomly select a dimension
function rand_partial_split(o::Omega{Float64}, depth::Int; ndims = 4)
  dimindices = collect(keys(o.intervals))
  randindices = unique([rand_select(dimindices) for i = 1:ndims])
  splitted = mid_partial_split(o::Omega{Float64}, randindices)
  transitionprob = 1/length(splitted)
  [(event, transitionprob) for event in splitted]
end


text:
include("Class_point.jl")
include("Class_iterate.jl")
include("misc.jl")
include("eval.jl")
include("Class_cutest.jl")
include("autotune.jl")


text:
using JLD2, Plots

cd(@__DIR__)
begin
    @load "results/collocation.jld2" x sol1 sol2
    x0 = deepcopy(x)
    sol_collo = [sol1, sol2]

    @load "results/l2.jld2" x sol1 sol2
    sol_l2 = [sol1, sol2]

    @load "results/l2_apt.jld2" x sol1 sol2
    sol_apt = [sol1, sol2]

    @load "results/lasso.jld2" x sol1 sol2
    sol_lasso = [sol1, sol2]
end

#--- case 1 ---#
plot(x, sol_lasso[1][:, 1, 1], lw=2, label="Lasso", xlabel="x", ylabel="density")
plot!(x, sol_l2[1][:, 1, 1], lw=2, label="Adaptive L²", line=:dash)
plot!(x0, sol_collo[1][:, 1, 1], label="Reference", line=:dashdot, lw=2, color=:gray27)
savefig("sod1_density_mean.pdf")

plot(x, sol_lasso[1][:, 2, 1], lw=2, label="Lasso", xlabel="x", ylabel="velocity", legend=:topleft)
plot!(x, sol_l2[1][:, 2, 1], lw=2, label="Adaptive L²", line=:dash)
plot!(x0, sol_collo[1][:, 2, 1], label="Reference", line=:dashdot, lw=2, color=:gray27)
savefig("sod1_velocity_mean.pdf")

plot(x, sol_lasso[1][:, 3, 1], lw=2, label="Lasso", xlabel="x", ylabel="temperature", legend=:topleft)
plot!(x, sol_l2[1][:, 3, 1], lw=2, label="Adaptive L²", line=:dash)
plot!(x0, sol_collo[1][:, 3, 1], label="Reference", line=:dashdot, lw=2, color=:gray27)
savefig("sod1_temperature_mean.pdf")

plot(x, sol_lasso[1][:, 1, 2], lw=2, label="Lasso", xlabel="x", ylabel="density")
plot!(x, sol_l2[1][:, 1, 2], lw=2, label="Adaptive L²", line=:dash)
plot!(x0, sol_collo[1][:, 1, 2], label="Reference", line=:dashdot, lw=2, color=:gray27)
savefig("sod1_density_std.pdf")

plot(x, sol_lasso[1][:, 2, 2], lw=2, label="Lasso", legend=:topleft, xlabel="x", ylabel="velocity")
plot!(x, sol_l2[1][:, 2, 2], lw=2, label="Adaptive L²", line=:dash)
plot!(x0, sol_collo[1][:, 2, 2], label="Reference", line=:dashdot, lw=2, color=:gray27)
savefig("sod1_velocity_std.pdf")

plot(x, sol_lasso[1][:, 3, 2], lw=2, label="Lasso", xlabel="x", ylabel="temperature", legend=:topleft)
plot!(x, sol_l2[1][:, 3, 2], lw=2, label="Adaptive L²", line=:dash)
plot!(x0, sol_collo[1][:, 3, 2], label="Reference", line=:dashdot, lw=2, color=:gray27)
savefig("sod1_temperature_std.pdf")

#--- case 2 ---#
plot(x, sol_lasso[2][:, 1, 1], lw=2, label="Lasso", xlabel="x", ylabel="density")
plot!(x, sol_l2[2][:, 1, 1], lw=2, label="Adaptive L²", line=:dash)
plot!(x0, sol_collo[2][:, 1, 1], label="Reference", line=:dashdot, lw=2, color=:gray27)
savefig("sod2_density_mean.pdf")

plot(x, sol_lasso[2][:, 2, 1], lw=2, label="Lasso", xlabel="x", ylabel="velocity")
plot!(x, sol_l2[2][:, 2, 1], lw=2, label="Adaptive L²", line=:dash)
plot!(x0, sol_collo[2][:, 2, 1], label="Reference", line=:dashdot, lw=2, color=:gray27)
savefig("sod2_velocity_mean.pdf")

plot(x, sol_lasso[2][:, 3, 1], lw=2, label="Lasso", xlabel="x", ylabel="temperature", legend=:topleft)
plot!(x, sol_l2[2][:, 3, 1], lw=2, label="Adaptive L²", line=:dash)
plot!(x0, sol_collo[2][:, 3, 1], label="Reference", line=:dashdot, lw=2, color=:gray27)
savefig("sod2_temperature_mean.pdf")

plot(x, sol_lasso[2][:, 1, 2], lw=2, label="Lasso", xlabel="x", ylabel="density")
plot!(x, sol_l2[2][:, 1, 2], lw=2, label="Adaptive L²", line=:dash)
plot!(x0, sol_collo[2][:, 1, 2], label="Reference", line=:dashdot, lw=2, color=:gray27)
savefig("sod2_density_std.pdf")

plot(x, sol_lasso[2][:, 2, 2], lw=2, legend=:topleft, label="Lasso", xlabel="x", ylabel="density")
plot!(x, sol_l2[2][:, 2, 2], lw=2, label="Adaptive L²", line=:dash)
plot!(x0, sol_collo[2][:, 2, 2], label="Reference", line=:dashdot, lw=2, color=:gray27)
savefig("sod2_velocity_std.pdf")

plot(x, sol_lasso[2][:, 3, 2], lw=2, label="Lasso", xlabel="x", ylabel="density")
plot!(x, sol_l2[2][:, 3, 2], lw=2, label="Adaptive L²", line=:dash)
plot!(x0, sol_collo[2][:, 3, 2], label="Reference", line=:dashdot, lw=2, color=:gray27)
savefig("sod2_temperature_std.pdf")

max_stars_repo_name: CSMMLab/FRSG

max_issues_repo_name: CSMMLab/FRSG

max_forks_repo_name: CSMMLab/FRSG


text:
"""
    struct SparseMatrixCSR{Bi,Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}

Matrix type for storing Bi-based sparse matrices
in the Compressed Sparse Row format. The standard
way of constructing SparseMatrixCSR is through the
[`sparsecsr`](@ref) function.
"""
struct SparseMatrixCSR{Bi,Tv,Ti} <: AbstractSparseMatrix{Tv,Ti}
    offset :: Int             # Bi-index offset (Bi-1)
    m      :: Int             # Number of columns
    n      :: Int             # Number of rows
    rowptr :: Vector{Ti}      # Row i is in (rowptr[i]+offset):(rowptr[i+1]-Bi)
    colval :: Vector{Ti}      # Col indices (with base Bi) of stored values
    nzval  :: Vector{Tv}      # Stored values, typically nonzeros

    function SparseMatrixCSR{Bi}(m::Integer, n::Integer, rowptr::Vector{Ti}, colval::Vector{Ti},
                                    nzval::Vector{Tv}) where {Bi,Tv,Ti<:Integer}
        @noinline throwsz(str, lbl, k) =
            throw(ArgumentError("number of $str ($lbl) must be ≥ 0, got $k"))
        m < 0 && throwsz("rows", 'm', m)
        n < 0 && throwsz("columns", 'n', n)
        rowptr[m+1] == length(nzval) && throwsz("nnz", "rowptr[m+1]", length(nzval))
        offset = Bi-1
        rowptr .+= offset
        colval .+= offset
        new{Bi,Tv,Ti}(Int(offset), Int(m), Int(n), rowptr, colval, nzval)
    end
end

SparseMatrixCSR(transpose::SparseMatrixCSC) where {Tv,Ti} =
        SparseMatrixCSR{1}(transpose.n, transpose.m, transpose.colptr, transpose.rowval, transpose.nzval)
SparseMatrixCSR{Bi}(transpose::SparseMatrixCSC{Tv,Ti}) where {Bi,Tv,Ti} =
        SparseMatrixCSR{Bi}(transpose.n, transpose.m, transpose.colptr, transpose.rowval, transpose.nzval)

size(S::SparseMatrixCSR) = (S.m, S.n)
show(io::IO, S::SparseMatrixCSR) = show(convert(IOContext, io), S::SparseMatrixCSR)

hasrowmajororder(::Type{<:SparseMatrixCSR}) = true
hasrowmajororder(a::SparseMatrixCSR) = hasrowmajororder(SparseMatrixCSR)
hascolmajororder(::Type{<:SparseMatrixCSR}) = false
hascolmajororder(a::SparseMatrixCSR) = hascolmajororder(SparseMatrixCSR)

getptr(S::SparseMatrixCSR) = S.rowptr
getindices(S::SparseMatrixCSR) = colvals(S)
nnz(S::SparseMatrixCSR{Bi}) where {Bi} = length(nonzeros(S))
count(pred, S::SparseMatrixCSR) = count(pred, nonzeros(S))
nonzeros(S::SparseMatrixCSR) = S.nzval
nzrange(S::SparseMatrixCSR{Bi}, row::Integer) where {Bi} = S.rowptr[row]-S.offset:S.rowptr[row+1]-Bi
rowvals(S::SparseMatrixCSR) = error("CSR sparse matrix does not contain raw row values")
colvals(S::SparseMatrixCSR) = S.colval

"""
    sparsecsr(args...)

Same args than `sparse`.
"""
sparsecsr(I,J,V,m,n,args...) =
        SparseMatrixCSR(sparse(J,I,V,n,m,args...))
sparsecsr(::Type{<:SparseMatrixCSR},I,J,args...) =
        sparsecsr(I,J,args...)
sparsecsr(::Type{<:SparseMatrixCSR{Bi}},I,J,V,m,n,args...) where {Bi} =
        SparseMatrixCSR{Bi}(sparse(J,I,V,n,m,args...))
sparsecsr(::Type{<:SparseMatrixCSR{Bi,Tv,Ti}},I::Vector{Ti},J::Vector{Ti},V::Vector{Tv},m,n,args...) where {Bi,Tv,Ti} =
        SparseMatrixCSR{Bi}(sparse(J,I,V,n,m,args...))

*(A::SparseMatrixCSR, v::Vector) = (y = similar(v,A.n);mul!(y,A,v))

convert(::Type{<:SparseMatrixCSR}, x::AbstractSparseMatrix)  = convert(SparseMatrixCSR{1}, x)

function show(io::IO, ::MIME"text/plain", S::SparseMatrixCSR)
    xnnz = nnz(S)
    print(io, S.m, "×", S.n, " ", typeof(S), " with ", xnnz, " stored ",
              xnnz == 1 ? "entry" : "entries")
    if xnnz != 0
        print(io, ":")
        show(IOContext(io, :typeinfo => eltype(S)), S)
    end
end

function show(io::IOContext, S::SparseMatrixCSR{Bi}) where{Bi}
    nnz(S) == 0 && return show(io, MIME("text/plain"), S)

    ioc = IOContext(io, :compact => true)
    function _format_line(r, col, padr, padc)
        print(ioc, "\n  [", rpad(col, padr), ", ", lpad(S.colval[r], padc), "]  =  ")
        if isassigned(S.nzval, Int(r))
            show(ioc, S.nzval[r])
        else
            print(ioc, Base.undef_ref_str)
        end
    end

    function _get_cols(from, to)
        idx = eltype(S.rowptr)[]
        c = searchsortedlast(S.rowptr, from)
        for i = from:to
            while i == S.rowptr[c+1]
                c +=1
            end
            push!(idx, c)
        end
        idx
    end

    rows = displaysize(io)[1] - 4 # -4 from [Prompt, header, newline after elements, new prompt]
    if !get(io, :limit, false) || rows >= nnz(S) # Will the whole matrix fit when printed?
        cols = _get_cols(Bi, nnz(S)+S.offset)
        padr, padc = ndigits.((maximum(S.colval[1:nnz(S)]), cols[end]))
        _format_line.(1:nnz(S), cols.+S.offset, padr, padc)
    else
        if rows <= 2
            print(io, "\n  \u22ee")
            return
        end
        s1, e1 = 1, div(rows - 1, 2) # -1 accounts for \vdots
        s2, e2 = nnz(S) - (rows - 1 - e1) + 1, nnz(S)
        cols1, cols2 = _get_cols(s1+S.offset, e1+S.offset), _get_cols(s2+S.offset, e2+S.offset)
        padr = ndigits(max(maximum(S.colval[s1:e1]), maximum(S.colval[s2:e2])))
        padc = ndigits(cols2[end])
        _format_line.(s1:e1, cols1, padr, padc)
        print(io, "\n  \u22ee")
        _format_line.(s2:e2, cols2, padr, padc)
    end
    return
end

function getindex(S::SparseMatrixCSR{Bi,Tv,Ti}, i0::Integer, i1::Integer) where {Bi,Tv,Ti}
    if !(1 <= i0 <= S.m && 1 <= i1 <= S.n); throw(BoundsError()); end
    r1 = Int(S.rowptr[i0]-S.offset)
    r2 = Int(S.rowptr[i0+1]-Bi)
    (r1 > r2) && return zero(Tv)
    i1_= i1+S.offset
    r1 = searchsortedfirst(S.colval, i1_, r1, r2, Base.Order.Forward)
    ((r1 > r2) || (S.colval[r1] != i1_)) ? zero(Tv) : S.nzval[r1]
end

function findnz(S::SparseMatrixCSR{Bi,Tv,Ti}) where {Bi,Tv,Ti}
    numnz = nnz(S)
    I = Vector{Ti}(undef, numnz)
    J = Vector{Ti}(undef, numnz)
    V = Vector{Tv}(undef, numnz)

    count = 1
    @inbounds for row in 1:S.m
        @inbounds for k in nzrange(S,row)
            I[count] = row
            J[count] = S.colval[k]-S.offset
            V[count] = S.nzval[k]
            count += 1
        end
    end

    return (I, J, V)
end

function sparse_from_coo(M::Type{<:SparseMatrixCSR}, I,J,V,m,n)
  sparsecsr(M, I,J,V,m,n)
end

function push_coo!(::Type{<:SparseMatrixCSR},
    I::Vector,J::Vector,V::Vector,ik::Integer,jk::Integer,vk::Number) where {Bi}

    push!(I, ik)
    push!(J, jk)
    push!(V, vk)
    nothing
end

@inline function is_entry_stored(::Type{<:SparseMatrixCSR},i::Integer,j::Integer)
  true
end

function finalize_coo!(::Type{<:SparseMatrixCSR},
    I::Vector,J::Vector,V::Vector,m::Integer,n::Integer)
end

function add_entry!(A::SparseMatrixCSR{Bi,Tv,Ti},v::Number,i::Integer,j::Integer,combine::Function=+) where {Bi,Tv,Ti<:Integer}
    cv = convert(Tv, v)
    ci = convert(Ti, i)
    cj = convert(Ti, j+A.offset)
    m,n = size(A)
    if !((1 <= i <= m) & (1 <= j <= n))
        throw(BoundsError(A, (i,j)))
    end
    rows = getptr(A)
    indices = getindices(A)
    vals = nonzeros(A)
    rowifirstk = Int(rows[ci]-A.offset)
    rowilastk = Int(rows[ci+1]-Bi)
    searchk = searchsortedfirst(indices, cj, rowifirstk, rowilastk, Base.Order.Forward)
    if searchk <= rowilastk-A.offset && indices[searchk] == cj
        # Column j contains entry A[i,j]. Update
        vold = vals[searchk]
        vals[searchk] = combine(vold, cv)
    else
        throw(ArgumentError("Entry ($i,$j) does not exists in the sparsity pattern"))
    end
    return A
end

function copy_entries!(a::SparseMatrixCSR,b::SparseMatrixCSR)
  _copy_entries_sparse!(a,b)
end

function mul!(y::AbstractVector,A::SparseMatrixCSR,v::AbstractVector{T}) where {T}
    A.n == size(v, 1) || throw(DimensionMismatch())
    A.m == size(y, 1) || throw(DimensionMismatch())

    y .= zero(T)
    for row = 1:size(y, 1)
        @inbounds for nz in nzrange(A,row)
            col = A.colval[nz]-A.offset
            y[row] += A.nzval[nz]*v[col]
        end
    end
    return y
end

function convert(::Type{<:SparseMatrixCSR{Bi}}, x::SparseMatrixCSR{xBi}) where {Bi,xBi}
    if Bi == xBi
        return x
    else
        return SparseMatrixCSR{Bi}( x.m,
                                    x.n,
                                    copy(getptr(x)).-x.offset,
                                    copy(getindices(x)).-x.offset,
                                    copy(nonzeros(x)))
    end
end

function convert(::Type{<:SparseMatrixCSR{Bi,Tv,Ti}}, x::SparseMatrixCSR{xBi,xTv,xTi}) where {Bi,Tv,Ti,xBi,xTv,xTi}
    if (Bi,Tv,Ti) == (xBi,xTv,xTi)
        return x
    else
        return SparseMatrixCSR{Bi}( x.m,
                                    x.n,
                                    convert(Vector{Ti}, copy(getptr(x)).-x.offset),
                                    convert(Vector{Ti}, copy(getindices(x)).-x.offset),
                                    convert(Vector{Tv}, copy(nonzeros(x))))
    end
end

function convert(::Type{<:SparseMatrixCSR{Bi}}, x::SparseMatrixCSC) where {Bi}
    A = sparse(transpose(x))
    (m, n) = size(A)
    return SparseMatrixCSR{Bi}(m, n, getptr(A), getindices(A), nonzeros(A))
end

function convert(::Type{<:SparseMatrixCSR{Bi,Tv,Ti}}, x::SparseMatrixCSC) where {Bi,Tv,Ti}
    A = sparse(transpose(x))
    return SparseMatrixCSR{Bi}( A.m,
                                A.n,
                                convert(Vector{Ti}, getptr(A)),
                                convert(Vector{Ti}, getindices(A)),
                                convert(Vector{Tv}, nonzeros(A)))
end

function convert(::Type{SparseMatrixCSC{Tv,Ti}}, x::SparseMatrixCSR{xBi,xTv,xTi}) where {Tv,Ti,xBi,xTv,xTi}
    A = sparse(transpose(x))
    return SparseMatrixCSR{Bi}( A.m,
                                A.n,
                                convert(Vector{Ti}, getptr(A)),
                                convert(Vector{Ti}, getindices(A)),
                                convert(Vector{Tv}, nonzeros(A)))
end


text:
# Tests for NEP types functionality

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test

@bench @testset "logger (ErrorLogger)" begin
    nep=nep_gallery("dep0")
    logger=ErrorLogger(1,100,0);
    λ1,x1 = newton(nep,logger=logger,v=ones(size(nep,1)),λ=0,tol=eps()*10)
    E=logger.errs[:,1];
    EE=E[(!isnan).(E)];
    # Check "essentially" quad convergence
    @test abs(log10(EE[end-1]^2)-log10(EE[end]))<2
    @test length(EE) > 5 #Has done more than 5 iterations
    @test EE[end] < eps()*10 # Error measure fulfills stopping criteria

end


text:
@doc raw"""
# Lotka-Volterra Model in 3D

The Lotka–Volterra model in 3D is an example of a Hamiltonian system with degenerate Poisson structure.

The equations read
```math
\begin{aligned}
\dot{q}_{1} &= q_{1} (           -  a_{2} q_{2} + a_{3} q_{3} - b_{2} + b_{3} ) , \\
\dot{q}_{2} &= q_{2} ( \hphantom{-} a_{1} q_{1} - a_{3} q_{3} + b_{1} - b_{3} ) , \\
\dot{q}_{3} &= q_{3} (           -  a_{1} q_{1} + a_{2} q_{2} - b_{1} + b_{2} ) , \\
\end{aligned}
```
which can be written in Poisson-form as
```math
\dot{q} = P(q) \nabla H(q) ,
```
with Poisson matrix
```math
P(q) = \begin{pmatrix}
            0 & - q_{1} q_{2} & \hphantom{-} q_{1} q_{3} \\
\hphantom{-} q_{1} q_{2} &             0 & - q_{2} q_{3} \\
- q_{1} q_{3} & \hphantom{-} q_{2} q_{3} & 0             \\
\end{pmatrix} ,
```
and Hamiltonian
```math
H(q) = a_{1} q_{1} + a_{2} q_{2} + a_{3} q_{3} + b_{1} \ln q_{1} + b_{2} \ln q_{2} + b_{3} \ln q_{3} .
```

References:
* A. M. Perelomov. Selected topics on classical integrable systems,
  Troisième cycle de la physique, expanded version of lectures
  delivered in May 1995.

* Yuri B. Suris. Integrable discretizations for lattice systems: local
  equations of motion and their Hamiltonian properties,
  Rev. Math. Phys. 11, pp. 727–822, 1999.
"""
module LotkaVolterra3d

    using GeometricBase
    using GeometricEquations

    export lotka_volterra_3d_ode

    export hamiltonian, casimir
    export compute_energy_error, compute_casimir_error

    Δt = 0.01
    nt = 1000

    const A1=+1.0
    const A2=+1.0
    const A3=+1.0

    const B1= 0.0
    const B2=+1.0
    const B3=+1.0

    const X0=1.0
    const Y0=1.0
    const Z0=2.0

    function v₁(t, q)
        q[1] * ( - A2 * q[2] + A3 * q[3] + B2 - B3)
    end

    function v₂(t, q)
        q[2] * ( + A1 * q[1] - A3 * q[3] - B1 + B3)
    end

    function v₃(t, q)
        q[3] * ( - A1 * q[1] + A2 * q[2] + B1 - B2)
    end

    const q₀=[X0, Y0, Z0]
    const v₀=[v₁(0, q₀), v₂(0, q₀), v₃(0, q₀)]

    function hamiltonian(t, q)
        A1*q[1] + A2*q[2] + A3*q[3] - B1*log(q[1]) - B2*log(q[2]) - B3*log(q[3])
    end

    hamiltonian_iode(t, q, v) = hamiltonian(t, q)

    function casimir(t, q)
        log(q[1]) + log(q[2]) + log(q[3])
    end

    function lotka_volterra_3d_v(t, q, v)
        v[1] = v₁(t, q)
        v[2] = v₂(t, q)
        v[3] = v₃(t, q)
        nothing
    end

    function lotka_volterra_3d_ode(q₀=q₀)
        ODE(lotka_volterra_3d_v, q₀; invariants=(h=hamiltonian,))
    end

    function compute_energy_error(t, q::AbstractDataSeries{<:AbstractVector{T}}) where {T}
        h = DataSeries(T, ntime(q))
        e = DataSeries(T, ntime(q))

        for i in axes(q,1)
            h[i] = hamiltonian(t[i], q[i])
            e[i] = (h[i] - h[0]) / h[0]
        end

        (h, e)
    end

    function compute_casimir_error(t, q::AbstractDataSeries{<:AbstractVector{T}}) where {T}
        c = DataSeries(T, ntime(q))
        e = DataSeries(T, ntime(q))

        for i in axes(q,1)
            c[i] = casimir(t[i], q[i])
            e[i] = (c[i] - c[0]) / c[0]
        end

        (c, e)
    end

end


text:
# GFF3 Record
# ===========

mutable struct Record
    # data and filled range
    data::Vector{UInt8}
    filled::UnitRange{Int}
    # :undefiend | :feature | :directive | :comment
    kind::Symbol
    # indexes
    seqid::UnitRange{Int}
    source::UnitRange{Int}
    type_::UnitRange{Int}
    start::UnitRange{Int}
    end_::UnitRange{Int}
    score::UnitRange{Int}
    strand::Int
    phase::Int
    attribute_keys::Vector{UnitRange{Int}}
end

"""
    GFF3.Record()

Create an unfilled GFF3 record.
"""
function Record()
    return Record(
        UInt8[], 1:0, :undefiend,
        # seqid-end_
        1:0, 1:0, 1:0, 1:0, 1:0,
        # score-attribute_keys
        1:0, 0, 0, UnitRange{Int}[])
end

"""
    GFF3.Record(data::Vector{UInt8})

Create a GFF3 record object from `data`.
This function verifies and indexes fields for accessors.
Note that the ownership of `data` is transferred to a new record object.
"""
function Record(data::Vector{UInt8})
    return convert(Record, data)
end

function Base.convert(::Type{Record}, data::Vector{UInt8})
    record = Record(
        data, 1:0, :undefiend,
        # seqid-end_
        1:0, 1:0, 1:0, 1:0, 1:0,
        # score-attribute_keys
        1:0, 0, 0, UnitRange{Int}[])
    index!(record)
    return record
end

"""
    GFF3.Record(str::AbstractString)

Create a GFF3 record object from `str`.
This function verifies and indexes fields for accessors.
"""
function Record(str::AbstractString)
    return convert(Record, str)
end

function Base.convert(::Type{Record}, str::AbstractString)
    return Record(Vector{UInt8}(str))
end

function Base.empty!(record::Record)
    record.filled = 1:0
    record.kind = :undefiend
    record.seqid = 1:0
    record.source = 1:0
    record.type_ = 1:0
    record.start = 1:0
    record.end_ = 1:0
    record.score = 1:0
    record.strand = 0
    record.phase = 0
    empty!(record.attribute_keys)
    return record
end

function Base.convert(::Type{String}, record::Record)
    return String(record.data[datarange(record)])
end

function Base.convert(::Type{Interval}, record::Record)
    name = BioGenerics.seqname(record)
    lpos = BioGenerics.leftposition(record)
    rpos = BioGenerics.rightposition(record)
    strd = hasstrand(record) ? GenomicFeatures.strand(record) : GenomicFeatures.STRAND_BOTH
    return Interval(name, lpos, rpos, strd, record)
end

function Base.convert(::Type{Interval{Record}}, record::Record)
    return convert(Interval, record)
end

function BioGenerics.isfilled(record::Record)
    return !isempty(record.filled)
end

function datarange(record::Record)
    return record.filled
end

function Base.:(==)(record1::Record, record2::Record)
    if isfilled(record1) == isfilled(record2) == true
        r1 = datarange(record1)
        r2 = datarange(record2)
        return length(r1) == length(r2) && memcmp(pointer(record1.data, first(r1)), pointer(record2.data, first(r2)), length(r1)) == 0
    end
    return isfilled(record1) == isfilled(record2) == false
end

function Base.copy(record::Record)
    return Record(
        record.data[record.filled],
        record.filled,
        record.kind,
        record.seqid,
        record.source,
        record.type_,
        record.start,
        record.end_,
        record.score,
        record.strand,
        record.phase,
        copy(record.attribute_keys))
end

function Base.write(io::IO, record::Record)
    return unsafe_write(io, pointer(record.data, first(record.filled)), length(record.filled))
end

function Base.print(io::IO, record::Record)
    write(io, record)
    return nothing
end

function Base.show(io::IO, record::Record)
    print(io, summary(record), ':')
    if isfilled(record)
        if isfeature(record)
            println(io)
            println(io, "       seqid: ", hasseqid(record) ? seqid(record) : "<missing>")
            println(io, "      source: ", hassource(record) ? source(record) : "<missing>")
            println(io, "        type: ", hasfeaturetype(record) ? featuretype(record) : "<missing>")
            println(io, "       start: ", hasseqstart(record) ? seqstart(record) : "<missing>")
            println(io, "         end: ", hasseqend(record) ? seqend(record) : "<missing>")
            println(io, "       score: ", hasscore(record) ? score(record) : "<missing>")
            println(io, "      strand: ", hasstrand(record) ? strand(record) : "<missing>")
            println(io, "       phase: ", hasphase(record) ? phase(record) : "<missing>")
              print(io, "  attributes:")
            for (key, vals) in attributes(record)
                print(io, ' ', key, '=', join(vals, ','))
            end
        elseif isdirective(record)
            println(io)
              print(io, "   directive: ", String(record))
        elseif iscomment(record)
            println(io)
              print(io, "     comment: ", String(record))
        else
            @assert false
        end
    else
        @assert record.kind == :undefiend
        print(io, " <not filled>")
    end
end

# Accessor functions
# ------------------

"""
    isfeature(record::Record)::Bool

Test if `record` is a feature record.
"""
function isfeature(record::Record)
    return record.kind == :feature
end

"""
    isdirective(record::Record)::Bool

Test if `record` is a directive record.
"""
function isdirective(record::Record)
    return record.kind == :directive
end

"""
    iscomment(record::Record)::Bool

Test if `record` is a comment record.
"""
function iscomment(record::Record)
    return record.kind == :comment
end

function checkfilled(record)
    if !isfilled(record)
        throw(ArgumentError("unfilled GFF3 record"))
    end
end

function checkkind(record::Record, kind::Symbol)
    if record.kind != kind
        throw(ArgumentError("$(kind) is expected but $(record.kind) is given"))
    end
end

"""
    seqid(record::Record)::String

Get the sequence ID of `record`.
"""
function seqid(record::Record)
    checkfilled(record)
    checkkind(record, :feature)
    if ismissing(record, record.seqid)
        missingerror(:seqname)
    end
    return decode(String(record.data[record.seqid]))
end

function hasseqid(record::Record)
    return record.kind == :feature && !ismissing(record, record.seqid)
end

function BioGenerics.seqname(record::Record)
    return seqid(record)
end

function BioGenerics.hasseqname(record::Record)
    return hasseqid(record)
end

"""
    source(record::Record)::String

Get the source of `record`.
"""
function source(record::Record)
    checkfilled(record)
    checkkind(record, :feature)
    if ismissing(record, record.source)
        missingerror(:source)
    end
    return decode(String(record.data[record.source]))
end

function hassource(record::Record)
    return record.kind == :feature && !ismissing(record, record.source)
end

"""
    featuretype(record::Record)::String

Get the type of `record`.
"""
function featuretype(record::Record)
    checkfilled(record)
    checkkind(record, :feature)
    if ismissing(record, record.type_)
        missingerror(:featuretype)
    end
    return decode(String(record.data[record.type_]))
end

function hasfeaturetype(record::Record)
    return record.kind == :feature && !ismissing(record, record.type_)
end

"""
    seqstart(record::Record)::Int

Get the start coordinate of `record`.
"""
function seqstart(record::Record)
    checkfilled(record)
    checkkind(record, :feature)
    if ismissing(record, record.start)
        missingerror(:start)
    end
    return unsafe_parse_decimal(Int, record.data, record.start)
end

function hasseqstart(record::Record)
    return record.kind == :feature && !ismissing(record, record.start)
end

function BioGenerics.leftposition(record::Record)
    return seqstart(record)
end

function BioGenerics.hasleftposition(record::Record)
    return hasseqstart(record)
end

"""
    seqend(record::Record)::Int

Get the end coordinate of `record`.
"""
function seqend(record::Record)
    checkfilled(record)
    checkkind(record, :feature)
    if ismissing(record, record.end_)
        missingerror(:end)
    end
    return unsafe_parse_decimal(Int, record.data, record.end_)
end

function hasseqend(record::Record)
    return record.kind == :feature && !ismissing(record, record.end_)
end

function BioGenerics.rightposition(record::Record)
    return seqend(record)
end

function BioGenerics.hasrightposition(record::Record)
    return hasseqend(record)
end

"""
    score(record::Record)::Float64

Get the score of `record`.
"""
function score(record::Record)
    checkfilled(record)
    checkkind(record, :feature)
    if ismissing(record, record.score)
        missingerror(:score)
    end
    # TODO: non-copy parsing
    return parse(Float64, String(record.data[record.score]))
end

function hasscore(record::Record)
    return record.kind == :feature && !ismissing(record, record.score)
end

"""
    strand(record::Record)::GenomicFeatures.Strand

Get the strand of `record`.
"""
function strand(record::Record)
    checkfilled(record)
    checkkind(record, :feature)
    return convert(GenomicFeatures.Strand, Char(record.data[record.strand]))
end

function hasstrand(record::Record)
    return record.kind == :feature && !ismissing(record, record.strand)
end

function GenomicFeatures.strand(record::Record)
    return strand(record)
end

"""
    phase(record::Record)::Int

Get the phase of `record`.
"""
function phase(record::Record)
    checkfilled(record)
    checkkind(record, :feature)
    if ismissing(record, record.phase)
        missingerror(:phase)
    end
    return unsafe_parse_decimal(Int, record.data, record.phase:record.phase)
end

function hasphase(record::Record)
    return record.kind == :feature && !ismissing(record, record.phase)
end

"""
    attributes(record::Record)::Vector{Pair{String,Vector{String}}}

Get the attributes of `record`.
"""
function attributes(record::Record)
    checkfilled(record)
    checkkind(record, :feature)
    ret = Pair{String,Vector{String}}[]
    for (i, key) in enumerate(record.attribute_keys)
        push!(ret, decode(String(record.data[key])) => attrvals(record, i))
    end
    return ret
end

function attrvals(record::Record, i::Int)
    key = record.attribute_keys[i]
    @assert record.data[last(key)+1] == UInt8('=')
    if i == lastindex(record.attribute_keys)
        valsrange = last(key)+2:last(datarange(record))
    else
        nextkey = record.attribute_keys[i+1]
        @assert record.data[first(nextkey)-1] == UInt8(';')
        valsrange = last(key)+2:first(nextkey)-2
    end
    stop = last(valsrange)
    ps = first(valsrange)
    pe = searchend(record.data, UInt8(','), ps, stop)
    vals = [decode(String(record.data[ps:pe-1]))]
    while pe ≤ stop
        ps = pe + 1
        pe = searchend(record.data, UInt8(','), ps, stop)
        push!(vals, decode(String(record.data[ps:pe-1])))
    end
    return vals
end

"""
    attributes(record::Record, key::String)::Vector{String}

Get the attributes of `record` with `key`.
"""
function attributes(record::Record, key::String)
    checkfilled(record)
    checkkind(record, :feature)
    for (i, k) in enumerate(record.attribute_keys)
        if isequaldata(key, record.data, k)
            return attrvals(record, i)
        end
    end
    throw(KeyError(key))
end

function searchend(data::Vector{UInt8}, b::UInt8, start::Int, stop::Int)
    p = start
    while p ≤ stop
        if data[p] == b
            return p
        end
        p += 1
    end
    return p
end

"""
    content(record::Record)::String

Get the content of `record`. Leading '#' letters are removed.
"""
function content(record::Record)
    checkfilled(record)
    lo = first(datarange(record))
    hi = last(datarange(record))
    if isfeature(record)
        return decode(String(record.data[lo:hi]))
    end
    if isdirective(record)
        return decode(String(record.data[lo+2:hi]))
    end
    if iscomment(record)
        return decode(String(record.data[lo+1:hi]))
    end
    @assert false
end

function is_fasta_directive(record::Record)
    return isdirective(record) && isequaldata("##FASTA", record.data, datarange(record))
end

function ismissing(record::Record, index::Int)
    return record.data[index] == UInt8('.')
end

function ismissing(record::Record, range::UnitRange{Int})
    return length(range) == 1 && record.data[first(range)] == UInt8('.')
end

# r"[-+]?[0-9]+" must match `data[range]`.
function unsafe_parse_decimal(::Type{T}, data::Vector{UInt8}, range::UnitRange{Int}) where T<:Signed
    lo = first(range)
    if data[lo] == UInt8('-')
        sign = T(-1)
        lo += 1
    elseif data[lo] == UInt8('+')
        sign = T(+1)
        lo += 1
    else
        sign = T(+1)
    end
    x = zero(T)
    @inbounds for i in lo:last(range)
        x = Base.Checked.checked_mul(x, 10 % T)
        x = Base.Checked.checked_add(x, (data[i] - UInt8('0')) % T)
    end
    return sign * x
end

# Decode the percent-encoded string if needed.
function decode(str::String)
    if '%' ∈ str
        return URIParser.unescape(str)
    end
    return str
end

# Check if `str == data[range]`
function isequaldata(str::String, data::Vector{UInt8}, range::UnitRange{Int})
    rlen = length(range)
    return rlen == sizeof(str) && memcmp(pointer(data, first(range)), pointer(str), rlen) == 0
end

function memcmp(p1::Ptr, p2::Ptr, n::Integer)
    return ccall(:memcmp, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), p1, p2, n)
end


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

module DataLakeStoreAccountManagementClient

using Random
using Dates
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_AccountsApi.jl")
include("api_FirewallRulesApi.jl")
include("api_LocationUsageApi.jl")
include("api_LocationsApi.jl")
include("api_OperationsApi.jl")
include("api_TrustedIdProvidersApi.jl")
include("api_VirtualNetworkRulesApi.jl")

# export models
export convert, CapabilityInformation
export convert, CheckNameAvailabilityParameters
export convert, CreateDataLakeStoreAccountParameters
export convert, CreateDataLakeStoreAccountProperties
export convert, CreateFirewallRuleWithAccountParameters
export convert, CreateOrUpdateFirewallRuleParameters
export convert, CreateOrUpdateFirewallRuleProperties
export convert, CreateOrUpdateTrustedIdProviderParameters
export convert, CreateOrUpdateTrustedIdProviderProperties
export convert, CreateOrUpdateVirtualNetworkRuleParameters
export convert, CreateOrUpdateVirtualNetworkRuleProperties
export convert, CreateTrustedIdProviderWithAccountParameters
export convert, CreateVirtualNetworkRuleWithAccountParameters
export convert, DataLakeStoreAccount
export convert, DataLakeStoreAccountBasic
export convert, DataLakeStoreAccountListResult
export convert, DataLakeStoreAccountProperties
export convert, DataLakeStoreAccountPropertiesBasic
export convert, EncryptionConfig
export convert, EncryptionIdentity
export convert, FirewallRule
export convert, FirewallRuleListResult
export convert, FirewallRuleProperties
export convert, KeyVaultMetaInfo
export convert, NameAvailabilityInformation
export convert, Operation
export convert, OperationDisplay
export convert, OperationListResult
export convert, Resource
export convert, SubResource
export convert, TrustedIdProvider
export convert, TrustedIdProviderListResult
export convert, TrustedIdProviderProperties
export convert, UpdateDataLakeStoreAccountParameters
export convert, UpdateDataLakeStoreAccountProperties
export convert, UpdateEncryptionConfig
export convert, UpdateFirewallRuleParameters
export convert, UpdateFirewallRuleProperties
export convert, UpdateFirewallRuleWithAccountParameters
export convert, UpdateKeyVaultMetaInfo
export convert, UpdateTrustedIdProviderParameters
export convert, UpdateTrustedIdProviderProperties
export convert, UpdateTrustedIdProviderWithAccountParameters
export convert, UpdateVirtualNetworkRuleParameters
export convert, UpdateVirtualNetworkRuleProperties
export convert, UpdateVirtualNetworkRuleWithAccountParameters
export convert, Usage
export convert, UsageListResult
export convert, UsageName
export convert, VirtualNetworkRule
export convert, VirtualNetworkRuleListResult
export convert, VirtualNetworkRuleProperties

# export operations
export convert, AccountsApi, FirewallRulesApi, LocationUsageApi, LocationsApi, OperationsApi, TrustedIdProvidersApi, VirtualNetworkRulesApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
"""
NCDatasets is a module to read and write NetCDF files.
This is minimal example how to read a NetCDF file.

```julia
using NCDatasets
# open the file and show its metadata (if called in the REPL without ending semicolon)
ds = NCDataset("filename.nc","r")
# load all data of the variable temperature
v = ds["temperature"][:,:]
# load the attribute units
unit = v.attrib["units"]
# close the file
close(ds)
```

More information is available at https://github.com/Alexander-Barth/NCDatasets.jl .
"""
module NCDatasets

using NetCDF_jll
using Dates
using Printf

using Base
using DataStructures: OrderedDict
import Base.convert

import Base: close
import Base: Array

using CFTime
export CFTime
export daysinmonth, daysinyear, yearmonthday, yearmonth, monthday
export dayofyear, firstdayofyear
export DateTimeStandard, DateTimeJulian, DateTimeProlepticGregorian,
    DateTimeAllLeap, DateTimeNoLeap, DateTime360Day, AbstractCFDateTime

const default_timeunits = "days since 1900-00-00 00:00:00"

const SymbolOrString = Union{Symbol, AbstractString}

include("CatArrays.jl")
include("types.jl")
include("colors.jl")
include("errorhandling.jl")
include("netcdf_c.jl")
include("dataset.jl")
include("attributes.jl")
include("dimensions.jl")
include("groupes.jl")
include("variable.jl")
include("cfvariable.jl")
include("cfconventions.jl")
include("defer.jl")
include("multifile.jl")
include("ncgen.jl")
include("precompile.jl")

export CatArrays

end # module


text:
using Quadrature , Distributions , Test

μ = [0.00 , 0.00]
Σ = [0.4 0.0 ; 0.00 0.4]
d = MvNormal(μ , Σ)
m(x , p) = pdf(d, x)
prob = QuadratureProblem(m, [-Inf , -Inf] , [Inf , Inf])
sol = solve(prob,HCubatureJL(),reltol=1e-3,abstol=1e-3)
@test  (1.00 - sol.u)^2  < 1e-6

μ = [0.00 , 0.00]
Σ = [0.4 0.0 ; 0.00 0.4]
d = MvNormal(μ , Σ)
m(x , p) = pdf(d, x)
prob = QuadratureProblem(m, [-Inf , 0.00] , [Inf , Inf])
sol = solve(prob,HCubatureJL(),reltol=1e-3,abstol=1e-3)
@test  (0.500 - sol.u)^2  < 1e-6

f(x , p) = pdf(Normal(0.00 , 1.00) , x)
prob = QuadratureProblem(f, -Inf , Inf)
sol = solve(prob,HCubatureJL(),reltol=1e-3,abstol=1e-3)
@test  (1.00 - sol.u)^2  < 1e-6

f(x , p) = pdf(Normal(0.00 , 1.00) , x)
prob = QuadratureProblem(f, -Inf , 0.00)
sol = solve(prob,HCubatureJL(),reltol=1e-3,abstol=1e-3)
@test  (0.50 - sol.u)^2  < 1e-6

f(x , p) = (1 / (x^2 + 1))
prob = QuadratureProblem(f, 0 , Inf)
sol = solve(prob,HCubatureJL(),reltol=1e-3,abstol=1e-3)
@test  (pi/2 - sol.u)^2  < 1e-6


text:
using LinearAlgebra  # ; similar to `import numpy` in Python
using Random
# import LinearAlgebra  # different from `using LinearAlgebra`

"""
Julia functions receive two types of arguments: arguments (args) and keyword arguments (kwargs)
"""
function args_and_kwargs(args...; kwargs...)  # args... is similar to `*args` in Python
    @show args
    @show kwargs
end

"""
`Random` and `LinearAlgebra` are basic packages provided by Julia.
Run REPL (Julia session) by e.g. `julia -q`.
Type `include("test/julia-tutorials/basic.jl")` in the REPL to read this script.
To run `main`, type `main()` in the REPL.

To reflect the updated codes without re-`include`,
use `includet` with [Revise.jl](https://github.com/timholy/Revise.jl).
I highly recommend you to read Revise.jl's documentation.
"""
function main(; seed=2021)
    Random.seed!(seed)  # how to control the random seed
    args_and_kwargs(1, 2; a=1, b=2)  # args = (1, 2), kwargs = Base.Pairs(:a => 1, :b => 2)

    x = [1, 2, 3]
    y = [4, 5, 6]
    @show dot(x, y)  # LinearAlgebra.dot; 1*4 + 2*5 + 3*6 = 32
    nothing  # similar to `None` in Python; the output of the last line is automatically returned if there is no `return` command
end


text:
"""
	grd2xyz(cmd0::String="", arg1=nothing, kwargs...)

Reads one 2-D grid and returns xyz-triplets.

Full option list at [`grd2xyz`]($(GMTdoc)grd2xyz.html)

Parameters
----------

- $(GMT.opt_J)
- **C** | **row_col** :: [Type => Bool]

    Replace the x- and y-coordinates on output with the corresponding column and row numbers.
    ($(GMTdoc)grd2xyz.html#c)
- $(GMT.opt_R)
- $(GMT.opt_V)
- **W** | **weight** :: [Type => Str]           `Arg = [a|weight]`

    Write out x,y,z,w, where w is the supplied weight (or 1 if not supplied) [Default writes x,y,z only].
    ($(GMTdoc)grd2xyz.html#w)
- **Z** | **flags** :: [Type => Str]

    Write a 1-column table. Output will be organized according to the specified ordering
    convention contained in ``flags``.
    ($(GMTdoc)grd2xyz.html#z)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_bo)
- $(GMT.opt_d)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_o)
- $(GMT.opt_s)
"""
function grd2xyz(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("grd2xyz", cmd0, arg1)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode

	cmd, = parse_common_opts(d, "", [:R :V_params :bo :d :f :h :o :s])
    cmd  = parse_these_opts(cmd, d, [[:C :row_col], [:W :weight], [:Z :flags]])

	common_grd(d, cmd0, cmd, "grd2xyz ", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
grd2xyz(arg1; kw...) = grd2xyz("", arg1; kw...)


text:
# Automatically generated using Clang.jl

# Skipping MacroDefinition: P2EST_QUADRANT_INIT ( q ) ( ( void ) memset ( ( q ) , - 1 , sizeof ( p2est_quadrant_t ) ) )
# Skipping MacroDefinition: P6EST_COLUMN_GET_RANGE ( q , f , l ) do { * ( f ) = ( size_t ) ( q ) -> p . piggy3 . local_num ; * ( l ) = * ( f ) + ( size_t ) ( q ) -> p . piggy3 . which_tree ; } while ( 0 ) ;
# Skipping MacroDefinition: P6EST_COLUMN_SET_RANGE ( q , f , l ) do { ( q ) -> p . piggy3 . local_num = ( p4est_locidx_t ) ( f ) ; ( q ) -> p . piggy3 . which_tree = ( p4est_topidx_t ) ( ( l ) - ( f ) ) ; } while ( 0 ) ;

struct p6est_connectivity
    conn4::Ptr{p4est_connectivity_t}
    top_vertices::Ptr{Cdouble}
    height::NTuple{3, Cdouble}
end

const p6est_connectivity_t = p6est_connectivity

struct p2est_quadrant
    z::p4est_qcoord_t
    level::Int8
    pad8::Int8
    pad16::Int16
    p::p6est_quadrant_data
end

const p2est_quadrant_t = p2est_quadrant

struct p6est
    mpicomm::MPI_Comm
    mpisize::Cint
    mpirank::Cint
    mpicomm_owned::Cint
    data_size::Csize_t
    user_pointer::Ptr{Cvoid}
    connectivity::Ptr{p6est_connectivity_t}
    columns::Ptr{p4est_t}
    layers::Ptr{sc_array_t}
    user_data_pool::Ptr{sc_mempool_t}
    layer_pool::Ptr{sc_mempool_t}
    global_first_layer::Ptr{p4est_gloidx_t}
    root_len::p4est_qcoord_t
end

const p6est_t = p6est
const p6est_init_t = Ptr{Cvoid}
const p6est_replace_t = Ptr{Cvoid}
const p6est_refine_column_t = Ptr{Cvoid}
const p6est_refine_layer_t = Ptr{Cvoid}
const p6est_coarsen_column_t = Ptr{Cvoid}
const p6est_coarsen_layer_t = Ptr{Cvoid}
const p6est_weight_t = Ptr{Cvoid}

@cenum p6est_comm_tag_t::UInt32 begin
    P6EST_COMM_PARTITION = 1
    P6EST_COMM_GHOST = 2
    P6EST_COMM_BALANCE = 3
end

struct p6est_ghost
    mpisize::Cint
    num_trees::p4est_topidx_t
    btype::p4est_connect_type_t
    column_ghost::Ptr{p4est_ghost_t}
    column_layer_offsets::Ptr{sc_array_t}
    ghosts::sc_array_t
    tree_offsets::Ptr{p4est_locidx_t}
    proc_offsets::Ptr{p4est_locidx_t}
    mirrors::sc_array_t
    mirror_tree_offsets::Ptr{p4est_locidx_t}
    mirror_proc_mirrors::Ptr{p4est_locidx_t}
    mirror_proc_offsets::Ptr{p4est_locidx_t}
    mirror_proc_fronts::Ptr{p4est_locidx_t}
    mirror_proc_front_offsets::Ptr{p4est_locidx_t}
end

const p6est_ghost_t = p6est_ghost


text:
"""
    @test_formany [ntests = 100] [maxtests = 1000] [logto = ""] argument_data, prop

Checks `prop` on randomly generated samples values specified in `argument_data`.
Returns `Pass` if `prop` holds for all values in the sample,
`Fail` if `prop` does not hold for some sampled value,
and `Error` if a part of the code could not be evaluated
or if insufficient non-vacuous tests are generated.

## Arguments

* `argument_data`:
  a sequence of comma-separated expressions,
  where each expression specifies
  a dummy variable, a universe from which to sample it,
  and a restricted comprehension of that universe.
  In particular, expressions must have one of the following forms:

  * `a<x::T<b` or `a>x::T>b`, where `x` is any symbol, treated as a dummy
    variable; `T` is a type; `a` & `b` are arbitrary expressions to constrain
    the range of values for `x` (provided that `a<x::T<b` is meaningful
    given `T` and the types of `a` and `b`);
  * `a<<x::T<b`, `a<x::T<<b`, `a>>x::T>b` or `a>x::T>>b`,
    as above, but where using `<<` will instruct `@test_formany`
    to generate test cases *away* from that boundary, i.e. closer
    to the other boundary;
  * `a<x<b` etc., where the type `T` is omitted,
    and therefore defaults to `Float64`

* `prop`: the property to be tested. Two types of properties are allowed:

  * any boolean-valued expression, referencing the dummy variables defined
    in `argument_data`, or
  * a property of the form `antecedent --> consequent`,
    where `antecedent` and `consequent` are themselves
    boolean-valued expressions that may reference dummy variables
    defined in `argument_data`. This form is useful for
    expressing "if-then" relationships without thinking through
    how to generate test cases that satisfy the antecedent,
    as `@test_formany` will do this for you; see below for more detail.

Optional keyword arguments include the following:

* `ntests = 100`: the number of tests that `test_formany` will perform
before returning `Pass`.

* `maxtests = 1000`: the maximal number of times to attempt to generate
test cases that satisfy `antecedent`, when the `antecedent --> consequent`
form is used.

* `logto = "path_to_file"`: the path to a file to which
a log of the actual test cases used will be recorded.
Since `@test_formany` is pseudo-random and might pass a test that you
expect to fail (or vice versa), this can be useful
for debugging to see which actual test cases were tested.

## Details

`@test_formany` generates a sample of test cases based on `argument_data`,
and tests `prop` on these values, returning `Pass` only if every test passes.

### Basic Use

When `prop` is just a boolean-valued expression,
`@test_formany` simply checks whether `prop` holds for the values
generated according to `argument_data`.

### Advanced Use: Conditional Properties and Vacuity Avoidance

When using `ntests = [some integer]`, with a `prop` of the form `antecedent --> consequent`,
`@test_formany` will attempt to test `ntests` cases where the
conditional is *not* vacuous, i.e. `ntests` cases where `antecedent`
evaluates to true.
However, to prevent the test from looping infinitely,
`@test_formany` uses `maxtests = [some integer]`
to limit the number of total attempts it will make.

Consider the following example: you might wish to test
this statement:

    ∀ 0 < x < 1, x^2 < 0.0001 --> x < 0.01

If you simply tested this statement with 100 randomly sampled
values between 0 and 1, you would be unlikely actually ever
to evaluate the consequent of the conditional,
because most of the time x^2 < 0.0001 would not be true.
So, you could use

    @test_formany ntests=100 0 < x::Float64 < 1, x^2 < 0.0001 --> x < 0.01

to guarantee 100 tests where `x^2 < 0.0001` was true.

However, by the same reasoning, you might actually be testing
an expression for which the antecedent is itself virtually always false.
Then you don't want the tests to continue forever---you want to specify
an upper bound on how many tests `@test_formany` should run before
it gives up.
For example, using `maxtests` in this case
prevents `@test_formany` from looping infinitely
looking for a value that satisfies an impossible antecedent:

    @test_formany ntests=100 maxtests=10_000 0 < x::Float64 < 1, (x < 1 && x > 2) --> x < 0.01

In this case, if the number of attempts to satisfy the antecedent
exceeds `maxtests` before `ntests` tests *do* satisfy the antecedent,
then `@test_formany` will return `Base.Test.Error`.

### Keyword Argument Defaults

- If no keywords are specified, `@test_formany` assumes `ntests = 100, maxtests = 1000`.
- If only `ntests` is specified, `@test_formany` assumes `maxtests = 10*ntests`.
- By default, `@test_formany` does *not* log output

Keyword arguments are order-invariant.

### Exhaustion

`@test_formany` tests a random sample of values specified
by `argument_data`, which in typical use will imply that it
tests a random sample of values drawn from a mathematically infinite
set.
Consequently, `@test_formany` is not exhaustive,
and can therefore generate false positives,
i.e. cases that generate `Pass` even when a counter-example
exists.

A similar and more complicated concern applies when using
conditional properties.
Here, the vacuity-avoidance arguments `ntests` and `maxtests`
may result in false positives and false negatives:
a conditional that is true, but with an unlikely-to-satisfy antecedent,
mail fail to generate enough test cases to test the entire
conditional, and therefore return `Error` even when the
property is true for the entire universe.
On the other hand, a counter-example may still
exists even when `@test_formany` returns `Pass`
on a conditional property after testing many non-vacuous antecedents.

## Examples

```julia
julia> @test_formany 100 > x::Float64 > 10, x+10 < y::Float64 << 1000, y-5 < z::Float64 < Inf, z > x+5
Test Passed
  Expression: (:((100 > x::Float64 > 10,x + 10 < y::Float64 << 1000,y - 5 < z::Float64 < Inf,z > x + 5)),)

julia> @test_formany ntests = 1000 10 < x::Float64 < 100, x+10 < y::Float64 << 1000, y-5 < z::Float64 < Inf, z > x+6
String["x = 25.18875813550991","y = 35.192306691353075","z = 31.183953098939906"]
Test Failed
  Expression: (:(ntests = 1000),:((100 > x::Float64 > 10,x + 10 < y::Float64 << 1000,y - 5 < z::Float64 < Inf,z > x + 6)))

## Failing to find enough non-vacuous test cases:
julia> @test_formany -1000 < x::Float64 < 1000, x > 999 --> x+1 > 1000
Error During Test
  Test threw an exception of type ErrorException
  Expression: (:((-1000 < x::Float64 < 1000,$(Expr(:-->, :(x > 999), :(x + 1 > 1000))))),)
  Found only 1/100 values satisfying given antecedent.

## Test that `log` is an increasing function
julia>  @test_formany ntests=1000 Inf > x::Float64 > 0, Inf > y::Float64 > 0, x < y --> log(x) < log(y)
Test Passed
  Expression: (:((Inf > x::Float64 > 0,Inf > y::Float64 > 0,$(Expr(:-->, :(x < y), :(log(x) < log(y)))))),:(ntests = 1000))

## Test that f(x) = x^3 is convex on (-100,100)
## Fails, because f(x) = x^3 is not convex on that interval!
julia> @test_formany 0 < a < 1,-100 < x < 100,-100 < y < 100,(a*x + (1-a)*y)^3 < a*x^3 + (1-a)*y^3
String["a = 0.39535367198058546","x = -13.538004569422625","y = 0.8504731053549079"]
Test Failed
  Expression: (:((0 < a < 1,-100 < x < 100,-100 < y < 100,(a * x + (1 - a) * y) ^ 3 < a * x ^ 3 + (1 - a) * y ^ 3)),)

## Test that f(x) = x^3 is convex on (0,100)
julia> @test_formany 0<a<1,0<x<100,0<y<100, (a*x+(1-a)*y)^3<=a*x^3+(1-a)*y^3
Test Passed
  Expression: (:((0 < a < 1,0 < x < 100,0 < y < 100,(a * x + (1 - a) * y) ^ 3 <= a * x ^ 3 + (1 - a) * y ^ 3)),)

## Test that Cobb-Douglas Utility is concave:
## First, define a helper function for concavity testing
julia> function respects_concavity(f,t,x1,x2)
                  f(t*x1+(1-t)*x2) > t * f(x1) + (1-t) * f(x2)
              end
## Same for weak concavity
julia> function respects_weak_concavity(f,t,x1,x2)
                  f(t*x1+(1-t)*x2) >= t * f(x1) + (1-t) * f(x2)
              end
## And one for Cobb-Douglas
julia> function cobb_douglas(a,x)
                  prod(x.^a)
              end

## The test below passes the @test_formany,
## since the probability of discovering a counter-example is zero.
julia> @test_formany 0<t<1, 0<a<1, 1<x1<10, 1<x2<10, 1<y1<10, 1<y2<10, respects_concavity(x -> cobb_douglas( [a,1-a], x), t, [x1,y1], [x2,y2])
Test Passed
  Expression: (:((0 < t < 1,0 < a < 1,1 < x1 < 10,1 < x2 < 10,1 < y1 < 10,1 < y2 < 10,respects_concavity((x->begin  # REPL[5], line 1:
                    cobb_douglas([a,1 - a],x)
                end),t,[x1,y1],[x2,y2]))),)

julia> @test_formany 0 < t < 1, 0 < a < 1, 1 < x1 < 10, 1 < x2 < 10, 1 < y1 < 10, 1 < y2 < 10,respects_weak_concavity(x -> cobb_douglas( [a,1-a], x), t, [x1,y1], [x2,y2])
Test Passed
  Expression: (:((0 < t < 1,0 < a < 1,1 < x1 < 10,1 < x2 < 10,1 < y1 < 10,1 < y2 < 10,respects_weak_concavity((x->begin  # REPL[6], line 1:
                    cobb_douglas([a,1 - a],x)
                end),t,[x1,y1],[x2,y2]))),)

## Handling Errors
julia> @test_formany -1 < x < 1, abs( log(x) ) > 0
Error During Test
  Test threw an exception of type DomainError
  Expression: (:((-1 < x < 1,abs(log(x)) > 0)),)
  DomainError
```
"""
macro test_formany(exprs...)
    outex = Expr(:macrocall, Symbol("@test_cases"))
    for ex in exprs
        push!(outex.args, esc(ex))
    end
    ex1 = Expr(:(=), :mode, :test_formany)
    push!(outex.args, esc(ex1))
    return outex
end


text:
function neural_ode(model,x,tspan,
                    args...;kwargs...)
  p = destructure(model)
  dudt_(du,u::TrackedArray,p,t) = du .= restructure(model,p)(u)
  dudt_(du,u::AbstractArray,p,t) = du .= Flux.data(restructure(model,p)(u))
  prob = ODEProblem(dudt_,x,tspan,p)
  return diffeq_adjoint(p,prob,args...;kwargs...)
end

function neural_ode_rd(model,x,tspan,
                    args...;kwargs...)
  Tracker.istracked(x) && error("u0 is not currently differentiable.")
  p = destructure(model)
  dudt_(u::TrackedArray,p,t) = restructure(model,p)(u)
  dudt_(u::AbstractArray,p,t) = Flux.data(restructure(model,p)(u))
  prob = ODEProblem(dudt_,x,tspan,p)
  return Flux.Tracker.collect(diffeq_rd(p,prob,args...;kwargs...))
end

neural_msde(x,model,mp,tspan,args...;kwargs...) = neural_msde(x,model,mp,tspan,
                                                         diffeq_fd,
                                                         args...;kwargs...)
function neural_dmsde(model,x,mp,tspan,
                      args...;kwargs...)
  Tracker.istracked(x) && error("u0 is not currently differentiable.")
  p = destructure(model)
  dudt_(u::TrackedArray,p,t) = restructure(model,p)(u)
  dudt_(u::AbstractArray,p,t) = Flux.data(restructure(model,p)(u))
  g(u,p,t) = mp.*u
  prob = SDEProblem(dudt_,g,x,tspan,p)
  return Flux.Tracker.collect(diffeq_rd(p,prob,args...;kwargs...))
end


text:
export isinvertible

@inline order(I::AlgAssAbsOrdIdl) = I.order

iszero(I::AlgAssAbsOrdIdl) = (I.iszero == 1)
isone(I::AlgAssAbsOrdIdl) = isone(abs(det(basis_mat(I, copy = false))))

function Base.one(S::AlgAssAbsOrdIdlSet)
  O = order(S)
  M = identity_matrix(FlintZZ, degree(O))
  return ideal(O, M, :twosided, true)
end

function one(I::AlgAssAbsOrdIdl)
  return ideal(order(I), one(order(I)))
end

function Base.copy(I::AlgAssAbsOrdIdl)
  return I
end

###############################################################################
#
#  String I/O
#
###############################################################################

function show(io::IO, a::AlgAssAbsOrdIdl)
  print(io, "Ideal of ")
  show(IOContext(io, :compact => true), order(a))
  println(io, " with basis matrix")
  print(io, a.basis_mat)
end

################################################################################
#
#  Deepcopy
#
################################################################################

function Base.deepcopy_internal(a::AlgAssAbsOrdIdl, dict::IdDict)
  b = typeof(a)(order(a))
  for i in fieldnames(typeof(a))
    if isdefined(a, i)
      if i != :order
        setfield!(b, i, Base.deepcopy_internal(getfield(a, i), dict))
      end
    end
  end
  return b
end

################################################################################
#
#  Ideal Set
#
################################################################################

function show(io::IO, a::AlgAssAbsOrdIdlSet)
  print(io, "Set of ideals of $(order(a))\n")
end

order(a::AlgAssAbsOrdIdlSet) = a.order

parent(I::AlgAssAbsOrdIdl) = AlgAssAbsOrdIdlSet(order(I))

function IdealSet(O::AlgAssAbsOrd)
   return AlgAssAbsOrdIdlSet(O)
end

elem_type(::Type{AlgAssAbsOrdIdlSet{S, T}}) where {S, T} = AlgAssAbsOrdIdl{S, T}

elem_type(::AlgAssAbsOrdIdlSet{S, T}) where {S, T} = AlgAssAbsOrdIdl{S, T}

parent_type(::Type{AlgAssAbsOrdIdl{S, T}}) where {S, T} = AlgAssAbsOrdIdlSet{S, T}

################################################################################
#
#  Basis (matrices)
#
################################################################################

function assure_has_basis(a::AlgAssAbsOrdIdl{S, T}) where {S, T}
  if isdefined(a, :basis)
    return nothing
  end

  O = order(a)
  a.basis = Vector{AlgAssAbsOrdElem{S, T}}(undef, degree(O))
  for i = 1:degree(O)
    a.basis[i] = elem_from_mat_row(O, basis_mat(a, copy = false), i)
  end
  return nothing
end

function assure_has_basis_mat_inv(a::AlgAssAbsOrdIdl)
  if isdefined(a, :basis_mat_inv)
    return nothing
  else
    a.basis_mat_inv = FakeFmpqMat(pseudo_inv(basis_mat(a, copy = false)))
    return nothing
  end
end

function basis(a::AlgAssAbsOrdIdl; copy::Bool = true)
  assure_has_basis(a)
  if copy
    return deepcopy(a.basis)
  else
    return a.basis
  end
end

function basis_mat(a::AlgAssAbsOrdIdl; copy::Bool = true)
  if copy
    return deepcopy(a.basis_mat)
  else
    return a.basis_mat
  end
end

function basis_mat_inv(a::AlgAssAbsOrdIdl; copy::Bool = true)
  assure_has_basis_mat_inv(a)
  if copy
    return deepcopy(a.basis_mat_inv)
  else
    return a.basis_mat_inv
  end
end

function det_of_basis_mat(a::AlgAssAbsOrdIdl)
  d = fmpz(1)
  M = basis_mat(a, copy = false)
  for i = 1:nrows(M)
    d = mul!(d, d, M[i, i])
  end
  return d
end

################################################################################
#
#  Arithmetic
#
################################################################################

function +(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where {S, T}
  if iszero(a)
    return deepcopy(b)
  elseif iszero(b)
    return deepcopy(a)
  end

  d = degree(order(a))
  M = vcat(basis_mat(a), basis_mat(b))
  deta = det_of_basis_mat(a)
  detb = det_of_basis_mat(b)
  M = hnf_modular_eldiv!(M, gcd(deta, detb), :lowerleft)
  M = sub(M, (d + 1):2*d, 1:d)
  return ideal(order(a), M, :nothing, true)
end

function *(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where {S, T}
  if iszero(a)
    return deepcopy(a)
  elseif iszero(b)
    return deepcopy(b)
  end

  d = degree(order(a))
  ba = basis(a, copy = false)
  bb = basis(b, copy = false)
  d2 = d^2

  M = zero_matrix(FlintZZ, d2, d)
  t = one(order(a))
  for i = 1:d
    for j = 1:d
      t = mul!(t, ba[i], bb[j])
      for k = 1:d
        M[(i - 1)*d + j, k] = coordinates(t; copy = false)[k]
      end
    end
  end

  deta = det_of_basis_mat(a)
  detb = det_of_basis_mat(b)
  M = hnf_modular_eldiv!(M, deta*detb, :lowerleft)
  M = sub(M, (d2 - d + 1):d2, 1:d)
  return ideal(order(a), M, :nothing, true)
end

^(A::AlgAssAbsOrdIdl, e::Int) = Base.power_by_squaring(A, e)
^(A::AlgAssAbsOrdIdl, e::fmpz) = Base.power_by_squaring(A, BigInt(e))

function intersect(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where {S, T}
  d = degree(order(a))
  H = vcat(basis_mat(a), basis_mat(b))
  K = left_kernel(H)[2]
  return ideal(order(a), _hnf(view(K, 1:d, 1:d)*basis_mat(a, copy = false), :lowerleft), :nothing, true)
end

################################################################################
#
#  Ad hoc multiplication
#
################################################################################

*(x::fmpz, a::AlgAssAbsOrdIdl) = ideal(order(a), x*basis_mat(a, copy = false), :nothing, true)
*(a::AlgAssAbsOrdIdl, x::fmpz) = ideal(order(a), basis_mat(a, copy = false)*x, :nothing, true)

function *(x::AlgAssAbsOrdElem, a::AlgAssAbsOrdIdl)
  @assert parent(x) === order(a)
  @assert isleft_ideal(a)
  O = order(a)
  return ideal(O, x)*a
end

function *(a::AlgAssAbsOrdIdl, x::AlgAssAbsOrdElem)
  @assert parent(x) === order(a)
  @assert isright_ideal(a)
  O = order(a)
  return a*ideal(O, x)
end

################################################################################
#
#  Division
#
################################################################################

function divexact(a::AlgAssAbsOrdIdl, x::fmpz)
  O = order(a)
  M = FakeFmpqMat(basis_mat(a, copy = false), deepcopy(x))
  if M.den != 1
    error("Ideal not divisible by x")
  end
  return ideal(O, M.num)
end

################################################################################
#
#  Construction
#
################################################################################

function ideal(O::AlgAssAbsOrd{S, T}, M::fmpz_mat, side::Symbol = :nothing, M_in_hnf::Bool = false) where {S, T}
  !M_in_hnf ? M = _hnf(M, :lowerleft) : nothing
  a = AlgAssAbsOrdIdl{S, T}(O, M)
  _set_sidedness(a, side)
  return a
end

function ideal(O::AlgAssAbsOrd{S, T}, x::AlgAssAbsOrdElem{S, T}) where {S, T}
  @assert parent(x) == O

  if iscommutative(O)
    a = ideal(O, x, :left)
    a.isright = 1
    return a
  end

  t1 = O()
  t2 = O()
  M = zero_matrix(FlintZZ, degree(O)^2, degree(O))
  b = basis(O, copy = false)
  for i = 1:degree(O)
    t1 = mul!(t1, b[i], x)
    ii = (i - 1)*degree(O)
    for j = 1:degree(O)
      t2 = mul!(t2, t1, b[j])
      elem_to_mat_row!(M, ii + j, t2)
    end
  end
  M = sub(_hnf(M, :lowerleft), nrows(M) - degree(O) + 1:nrows(M), 1:ncols(M))

  return ideal(O, M, :twosided, true)
end

function ideal(O::AlgAssAbsOrd{S, T}, x::AlgAssAbsOrdElem{S, T}, action::Symbol) where {S, T}
  @assert parent(x) == O
  M = representation_matrix(x, action)
  a = ideal(O, M)
  if iszero(x)
    a.iszero = 1
  end

  if action == :left
    a.isright = 1
  elseif action == :right
    a.isleft = 1
  end

  return a
end

*(O::AlgAssAbsOrd{S, T}, x::AlgAssAbsOrdElem{S, T}) where {S, T} = ideal(O, x, :left)
*(x::AlgAssAbsOrdElem{S, T}, O::AlgAssAbsOrd{S, T}) where {S, T} = ideal(O, x, :right)

function ideal_from_z_gens(O::AlgAssAbsOrd, b::Vector{T}, side::Symbol = :nothing) where { T <: AlgAssAbsOrdElem }
  d = degree(O)
  @assert length(b) >= d

  M = zero_matrix(FlintZZ, length(b), d)
  for i = 1:length(b)
    elem_to_mat_row!(M, i, b[i])
  end
  M = _hnf(M, :lowerleft)
  if d < length(b)
    M = sub(M, (length(b) - d + 1):length(b), 1:d)
  end
  return ideal(O, M, side, true)
end

################################################################################
#
#  Extend/contract
#
################################################################################

function extend(A::AlgAssAbsOrdIdl, O::AlgAssAbsOrd, action::Symbol = :twosided)
  # Assumes order(A) \subseteq O

  d = degree(O)
  X = map(O, basis(A, copy = false))
  Y = basis(O, copy = false)
  if action == :left || action == :twosided
    M = zero_matrix(FlintZZ, d^2, d)
    t = O()
    for i = 1:d
      ii = (i - 1)*d
      for j = 1:d
        t = Y[i]*X[j]
        elem_to_mat_row!(M, ii + j, t)
      end
    end
    M = sub(_hnf(M, :lowerleft), nrows(M) - d + 1: nrows(M), 1:d)
    if action == :left
      return ideal(O, M, :left, true)
    end
    for i = 1:d
      X[i] = elem_from_mat_row(O, M, i)
    end
  end
  if action == :right || action == :twosided
    M = zero_matrix(FlintZZ, d^2, d)
    t = O()
    for i = 1:d
      ii = (i - 1)*d
      for j = 1:d
        t = X[i]*Y[j]
        elem_to_mat_row!(M, ii + j, t)
      end
    end
    M = sub(_hnf(M, :lowerleft), nrows(M) - d + 1: nrows(M), 1:d)
    return ideal(O, M, action, true)
  end
  error("action should be :left, :right or :twosided")
end

*(A::AlgAssAbsOrdIdl, O::AlgAssAbsOrd) = extend(A, O, :right)
*(O::AlgAssAbsOrd, A::AlgAssAbsOrdIdl) = extend(A, O, :left)

function contract(A::AlgAssAbsOrdIdl, O::AlgAssAbsOrd)
  # Assumes O \subseteq order(A)

  d = degree(O)
  M = basis_mat(O, copy = false)*basis_mat_inv(order(A), copy = false)
  @assert M.den == 1
  H = vcat(basis_mat(A), M.num)
  K = left_kernel(H)[2]
  M = sub(K, 1:d, 1:d)*basis_mat(A, copy = false)
  M = M*basis_mat(order(A), copy = false)*basis_mat_inv(O, copy = false)
  @assert M.den == 1
  M = _hnf(M.num, :lowerleft)
  return ideal(O, M, :nothing, true)
end

intersect(O::AlgAssAbsOrd, A::AlgAssAbsOrdIdl) = contract(A, O)
intersect(A::AlgAssAbsOrdIdl, O::AlgAssAbsOrd) = contract(A, O)

################################################################################
#
#  Move ideals to another algebra
#
################################################################################

# Let m: A -> B be an injection and I an ideal in B. Then this computes the
# ideal O_A\cap I, where O_A is the maximal order of A.
function _as_ideal_of_smaller_algebra(m::AbsAlgAssMor, I::AlgAssAbsOrdIdl)
  A = domain(m)
  B = codomain(m)
  @assert dim(A) <= dim(B)
  @assert algebra(order(I)) == B
  OA = maximal_order(A)
  OB = maximal_order(B)
  IB = extend(I, OB)
  # Transport OA to B
  M = zero_matrix(FlintZZ, dim(B), dim(B))
  for i = 1:dim(A)
    t = OB(m(elem_in_algebra(basis(OA, copy = false)[i])))
    elem_to_mat_row!(M, i, t)
  end
  # Compute the intersection of M and IB
  H = vcat(M, basis_mat(IB))
  K = left_kernel(H)[2]
  N = view(K, 1:dim(B), 1:dim(B))*M
  # Map the basis to A
  basis_in_A = Vector{elem_type(OA)}(undef, dim(B))
  for i = 1:dim(B)
    t = elem_from_mat_row(OB, N, i)
    b, s = haspreimage(m, elem_in_algebra(t, copy = false))
    @assert b
    basis_in_A[i] = OA(s)
  end
  J = ideal_from_z_gens(OA, basis_in_A)
  return J
end

# Let m: A -> B be an injection and O an order in B. Then this computes the
# order O_A\cap O, where O_A is the maximal order of A.
function _as_order_of_smaller_algebra(m::AbsAlgAssMor, O::AlgAssAbsOrd)
  A = domain(m)
  B = codomain(m)
  @assert dim(A) <= dim(B)
  @assert algebra(O) == B
  OA = maximal_order(A)
  OB = maximal_order(B)
  # Transport OA to B
  M = zero_matrix(FlintZZ, dim(B), dim(B))
  for i = 1:dim(A)
    t = OB(m(elem_in_algebra(basis(OA, copy = false)[i])))
    elem_to_mat_row!(M, i, t)
  end
  # Compute the intersection of M and O (in OB)
  N = basis_mat(O, copy = false)*basis_mat_inv(OB, copy = false)
  @assert N.den == 1
  H = vcat(M, N.num)
  K = left_kernel(H)[2]
  N = sub(K, 1:dim(B), 1:dim(B))*M
  # Map the basis to A
  basis_in_A = Vector{elem_type(A)}(undef, dim(B))
  for i = 1:dim(B)
    t = elem_from_mat_row(OB, N, i)
    b, s = haspreimage(m, elem_in_algebra(t, copy = false))
    @assert b
    basis_in_A[i] = s
  end
  M = zero_matrix(base_ring(A), length(basis_in_A), dim(A))
  for i = 1:length(basis_in_A)
    elem_to_mat_row!(M, i, basis_in_A[i])
  end
  M = FakeFmpqMat(M)
  MM = sub(hnf!(M, :lowerleft), (dim(B) - dim(A) + 1):dim(B), 1:dim(A))
  OO = Order(A, MM)
  return OO
end

# Let m: A -> B be an injection, I an ideal in A and O an order of B. Then this
# computes I*O.
# WE ASSUME THAT m(A) LIES IN THE CENTRE OF B!
function _as_ideal_of_larger_algebra(m::AbsAlgAssMor, I::AlgAssAbsOrdIdl, O::AlgAssAbsOrd)
  A = domain(m)
  B = codomain(m)
  @assert dim(A) <= dim(B)
  @assert algebra(order(I)) == A
  @assert algebra(O) == B

  M = zero_matrix(FlintZZ, dim(A)*dim(B), dim(B))
  X = Vector{elem_type(O)}(undef, dim(A))
  for i = 1:dim(A)
    t = elem_in_algebra(basis(I, copy = false)[i])
    X[i] = O(m(t))
  end
  for i = 1:dim(A)
    for j = 1:dim(B)
      t = X[i]*basis(O, copy = false)[j]
      elem_to_mat_row!(M, (i - 1)*dim(B) + j, t)
    end
  end
  M = sub(_hnf(M, :lowerleft), dim(B)*(dim(A) - 1) + 1:dim(A)*dim(B), 1:dim(B))
  return ideal(O, M, :nothing, true)
end

################################################################################
#
#  Inclusion of elements in ideals
#
################################################################################

function in(x::AlgAssAbsOrdElem, I::AlgAssAbsOrdIdl)
  el = coordinates(x)
  y = matrix(FlintZZ, 1, length(el), el)
  M1, d =pseudo_inv(I.basis_mat)
  if FakeFmpqMat(y*M1, d).den==1
    return true
  else
    return false
  end
end

################################################################################
#
#  Equality
#
################################################################################

function ==(A::AlgAssAbsOrdIdl, B::AlgAssAbsOrdIdl)
  order(A) !== order(B) && return false
  return basis_mat(A, copy = false) == basis_mat(B, copy = false)
end

################################################################################
#
#  Hashing
#
################################################################################

function Base.hash(A::AlgAssAbsOrdIdl, h::UInt)
  return Base.hash(basis_mat(A, copy = false), h)
end

################################################################################
#
#  Idempotents
#
################################################################################

# Algorithm 1.3.2 in Cohen "Advanced Topics in Computational Number Theory"
function idempotents(a::AlgAssAbsOrdIdl, b::AlgAssAbsOrdIdl)
  !(order(a) === order(b)) && error("Parent mismatch")
  O = order(a)
  d = degree(O)
  A = basis_mat(a, copy = false)
  B = basis_mat(b, copy = false)

  C = vcat(A, B)
  H, U = hnf_with_transform(C)

  if H != vcat(identity_matrix(FlintZZ, d), zero_matrix(FlintZZ, d, d))
    error("Ideals are not coprime")
  end

  X = sub(U, 1:1, 1:d)
  XA = X*A
  x = O([ XA[1, i] for i = 1:d ])
  @assert x in a
  y = one(O) - x
  @assert y in b
  return x, y
end

################################################################################
#
#  From algebra to number field
#
################################################################################

function _as_ideal_of_number_field(I::AlgAssAbsOrdIdl, m::AbsAlgAssToNfAbsMor)
  @assert algebra(order(I)) == domain(m)
  K = codomain(m)
  OK = maximal_order(K)

  b = Vector{elem_type(OK)}()
  for i = 1:dim(domain(m))
    push!(b, OK(m(elem_in_algebra(basis(I, copy = false)[i]))))
  end
  return ideal_from_z_gens(OK, b)
end

function _as_ideal_of_number_field(I::FacElem{<: AlgAssAbsOrdIdl, <: AlgAssAbsOrdIdlSet}, m::AbsAlgAssToNfAbsMor)
  @assert algebra(order(base_ring(parent(I)))) == domain(m)
  K = codomain(m)
  OK = maximal_order(K)

  if isempty(I)
    return FacElemMon(IdealSet(OK))()
  end

  bases = Vector{ideal_type(OK)}()
  exps = Vector{fmpz}()
  for (b, e) in I
    push!(bases, _as_ideal_of_number_field(b, m))
    push!(exps, e)
  end
  return FacElem(bases, exps)
end

function _as_ideal_of_algebra(I::NfAbsOrdIdl, i::Int, O::AlgAssAbsOrd, one_ideals::Vector{Vector{T}}) where { T <: AlgAssAbsOrdElem }
  A = algebra(O)
  fields_and_maps = as_number_fields(A)
  b = Vector{elem_type(O)}()
  for j = 1:length(fields_and_maps)
    K, AtoK = fields_and_maps[j]
    if j == i
      @assert nf(order(I)) == K
      append!(b, [ O(AtoK\K(bb)) for bb in basis(I, copy = false) ])
    else
      append!(b, one_ideals[j])
    end
  end
  return ideal_from_z_gens(O, b)
end

function _as_ideal_of_algebra(I::FacElem{NfOrdIdl, NfOrdIdlSet}, i::Int, O::AlgAssAbsOrd, one_ideals::Vector{Vector{T}}) where { T <: AlgAssAbsOrdElem }
  if isempty(I)
    return FacElemMon(IdealSet(O))()
  end

  base = Vector{ideal_type(O)}()
  exp = Vector{fmpz}()
  for (b, e) in I
    push!(base, _as_ideal_of_algebra(b, i, O, one_ideals))
    push!(exp, e)
  end
  return FacElem(base, exp)
end

# Returns an array of bases of the ideals O_i(1)*O_i lifted to O, where O_i
# are the maximal orders of the number fields in the decomposition of
# algebra(O).
function _lift_one_ideals(O::AlgAssAbsOrd)
  A = algebra(O)
  fields_and_maps = as_number_fields(A)

  one_ideals = Vector{Vector{elem_type(O)}}()
  for i = 1:length(fields_and_maps)
    K, AtoK = fields_and_maps[i]
    OK = maximal_order(K)
    a = OK(1)*OK
    b = Vector{elem_type(O)}()
    for i = 1:degree(K)
      push!(b, O(AtoK\K(basis(a, copy = false)[i])))
    end
    push!(one_ideals, b)
  end
  return one_ideals
end

################################################################################
#
#  Factorization
#
################################################################################

function factor(I::AlgAssAbsOrdIdl)
  O = order(I)
  A = algebra(O)
  fields_and_maps = as_number_fields(A)

  one_ideals = _lift_one_ideals(O)

  fac = Dict{ideal_type(O), Int}()
  for i = 1:length(fields_and_maps)
    K, AtoK = fields_and_maps[i]
    J = _as_ideal_of_number_field(I, AtoK)
    fac2 = factor(J)
    for (p, e) in fac2
      P = _as_ideal_of_algebra(p, i, O, one_ideals)
      fac[P] = e
    end
  end
  return fac
end

################################################################################
#
#  Inverse
#
################################################################################

# If I is not coprime to the conductor of O in the maximal order, then this might
# not be an inverse.
function inv(a::AlgAssAbsOrdIdl)
  @assert !iszero(a)
  O = order(a)
  return colon(ideal(O, one(O)), a)
end

# Mostly the same as colon in NfOrd/Ideal/Ideal.jl
function colon(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where {S, T}
  O = order(a)
  n = degree(O)
  B = basis(b)

  bmatinv = basis_mat_inv(a, copy = false)

  n = FakeFmpqMat(representation_matrix(B[1]), FlintZZ(1))*bmatinv
  m = numerator(n)
  d = denominator(n)
  for i in 2:length(B)
    n = FakeFmpqMat(representation_matrix(B[i]), FlintZZ(1))*bmatinv
    l = lcm(denominator(n), d)
    if l == d
      m = hcat(m, n.num)
    else
      m = hcat(m*div(l, d), n.num*div(l, denominator(n)))
      d = l
    end
  end
  m = hnf(transpose(m))
  # n is upper right HNF
  m = transpose(sub(m, 1:degree(O), 1:degree(O)))
  b, l = pseudo_inv(m)
  return frac_ideal_type(O)(O, ideal(O, b), l)
end

function isinvertible(a::AlgAssAbsOrdIdl)
  if iszero(a)
    return false, a
  end

  bmata = basis_mat(a, copy = false)

  if any(iszero(bmata[i, i]) for i in 1:degree(order(a)))
    return false, a
  end

  if order(a).ismaximal == 1
    return true, inv(a)
  end

  b = inv(a)
  c = a*b
  return isone(c), b
end


text:
#=
A regular log-λ grid is characterized by:
log10(λ_i+1) - log10(λ_i)  =  log10(λ_i+1 / λ_i)  =  costant step

hence:
λ_i+1 / λ_i  =  cost.
λ_i+1 / λ_i - 1  =  cost.
(λ_i+1 - λ_i) / λ_i  =  Δλ / λ  =  1/R   is constant ∀i

i.e. it has a constant spectral resolution given by:
R = c / σ_kms

The step of the grid is:
log10(σ_kms / c + 1)
=#

function conv_gauss(λ, data, σ_kms; oversampling::Int=1)
    @assert oversampling >= 1

    # Interpolate input data on the regular log-λ grid with proper resolution
    # check with: extrema([(grid_λ[i+1] - grid_λ[i]) / grid_λ[i] * 3e5 for i in 1:length(grid_λ)-1]) .* oversampling, σ_kms
    step = log10(σ_kms / 3e5 + 1) / oversampling
    grid_λ = 10. .^ range(log10(minimum(λ)), log10(maximum(λ)), step=step)
    grid_y = Spline1D(λ, data, k=1)(grid_λ)

    # Convolution with Gaussian kernel and extraction of proper subset
    kernel = gauss.(-5:1. / oversampling:5, 0., 1.)
    i1 = div(length(kernel), 2) + 1
    i2 = i1 + length(grid_y) - 1
    out = conv(grid_y, kernel)[i1:i2] ./ oversampling

    # Interpolate back to original domain
    out = Spline1D(grid_λ, out, k=1, bc="extrapolate")(λ)

    # Replace data close to the edges
    ee = (1 + 2σ_kms / 3e5) * minimum(λ)
    i = findall(λ .< ee)
    out[i] = Spline1D(grid_λ, grid_y, k=1, bc="extrapolate")(λ[i])

    ee = (1 - 2σ_kms / 3e5) * maximum(λ)
    i = findall(λ .> ee)
    out[i] = Spline1D(grid_λ, grid_y, k=1, bc="extrapolate")(λ[i])

    # @gp    λ data "w lp t 'Data'"
    # @gp :- λ out  "w lp t 'Convolution'"
    return out
end

function delta_conv_gauss(λ, data, σ_kms; oversampling::Int=1)
    step = log10(σ_kms / 3e5 + 1) / oversampling
    grid_λ = 10. .^ range(log10(minimum(λ) * (1 - 5 * σ_kms/3e5)),
                          log10(maximum(λ) * (1 + 5 * σ_kms/3e5)),
                          step=step)
    out = grid_λ .* 0.

    for ii in 1:length(data)
        out .+= data[ii] .* gauss(grid_λ, λ[ii], σ_kms / 3e5 * λ[ii])
    end

    return grid_λ, out
end

#=

x, y = QSFit.delta_conv_gauss([1e3, 3e3, 1e4, 3e4], [0., 1., 1., 0.], 3000.)
QSFit.int_tabulated(x, y)  # should be ~ 2
@gp xlog=true x y "w lp"

y2 = QSFit.conv_gauss(x, y, 10000.)
@gp :- x y2 "w lp"
QSFit.int_tabulated(x, y)  # still ~ 2

=#