text:
"""
    timeevolution.master_h(tspan, rho0, H, J; <keyword arguments>)

Integrate the master equation with dmaster_h as derivative function.

Further information can be found at [`master`](@ref).
"""
function master_h(tspan, rho0::Operator, H::AbstractOperator, J;
                rates=nothing,
                Jdagger=dagger.(J),
                fout=nothing,
                kwargs...)
    check_master(rho0, H, J, Jdagger, rates)
    tmp = copy(rho0)
    dmaster_(t, rho, drho) = dmaster_h!(drho, H, J, Jdagger, rates, rho, tmp)
    integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
end

"""
    timeevolution.master_nh(tspan, rho0, H, J; <keyword arguments>)

Integrate the master equation with dmaster_nh as derivative function.

In this case the given Hamiltonian is assumed to be the non-hermitian version:
```math
H_{nh} = H - \\frac{i}{2} \\sum_k J^†_k J_k
```
Further information can be found at [`master`](@ref).
"""
function master_nh(tspan, rho0::Operator, Hnh::AbstractOperator, J;
                rates=nothing,
                Hnhdagger::AbstractOperator=dagger(Hnh),
                Jdagger=dagger.(J),
                fout=nothing,
                kwargs...)
    check_master(rho0, Hnh, J, Jdagger, rates)
    tmp = copy(rho0)
    dmaster_(t, rho, drho) = dmaster_nh!(drho, Hnh, Hnhdagger, J, Jdagger, rates, rho, tmp)
    integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
end

"""
    timeevolution.master(tspan, rho0, H, J; <keyword arguments>)

Time-evolution according to a master equation.

There are two implementations for integrating the master equation:

* [`master_h`](@ref): Usual formulation of the master equation.
* [`master_nh`](@ref): Variant with non-hermitian Hamiltonian.

For dense arguments the `master` function calculates the
non-hermitian Hamiltonian and then calls master_nh which is slightly faster.

# Arguments
* `tspan`: Vector specifying the points of time for which output should
        be displayed.
* `rho0`: Initial density operator. Can also be a state vector which is
        automatically converted into a density operator.
* `H`: Arbitrary operator specifying the Hamiltonian.
* `J`: Vector containing all jump operators which can be of any arbitrary
        operator type.
* `rates=nothing`: Vector or matrix specifying the coefficients (decay rates)
        for the jump operators. If nothing is specified all rates are assumed
        to be 1.
* `Jdagger=dagger.(J)`: Vector containing the hermitian conjugates of the jump
        operators. If they are not given they are calculated automatically.
* `fout=nothing`: If given, this function `fout(t, rho)` is called every time
        an output should be displayed. ATTENTION: The given state rho is not
        permanent! It is still in use by the ode solver and therefore must not
        be changed.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master(tspan, rho0::Operator, H::AbstractOperator, J;
                rates=nothing,
                Jdagger=dagger.(J),
                fout=nothing,
                kwargs...)
    isreducible = check_master(rho0, H, J, Jdagger, rates)
    if !isreducible
        tmp = copy(rho0)
        dmaster_h_(t, rho, drho) = dmaster_h!(drho, H, J, Jdagger, rates, rho, tmp)
        return integrate_master(tspan, dmaster_h_, rho0, fout; kwargs...)
    else
        Hnh = nh_hamiltonian(H,J,Jdagger,rates)
        Hnhdagger = dagger(Hnh)
        tmp = copy(rho0)
        dmaster_nh_(t, rho, drho) = dmaster_nh!(drho, Hnh, Hnhdagger, J, Jdagger, rates, rho, tmp)
        return integrate_master(tspan, dmaster_nh_, rho0, fout; kwargs...)
    end
end

"""
    timeevolution.master(tspan, rho0, L; <keyword arguments>)

Time-evolution according to a master equation with a Liouvillian superoperator `L`.
# Arguments
* `tspan`: Vector specifying the points of time for which output should
        be displayed.
* `rho0`: Initial density operator. Can also be a state vector which is
        automatically converted into a density operator.
* `L`: Superoperator defining the right-hand-side of the master equation.
* `fout=nothing`: If given, this function `fout(t, rho)` is called every time
        an output should be displayed. ATTENTION: The given state rho is not
        permanent! It is still in use by the ode solver and therefore must not
        be changed.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master(tspan, rho0::Operator, L::SuperOperator; fout=nothing, kwargs...)
    # Rewrite rho as Ket and L as Operator
    dim = length(rho0.basis_l)*length(rho0.basis_r)
    b = GenericBasis(dim)
    rho_ = Ket(b,reshape(rho0.data, dim))
    L_ = Operator(b,b,L.data)
    dmaster_(t,rho,drho) = dmaster_liouville!(drho,L_,rho)

    # Rewrite into density matrix when saving
    tmp = copy(rho0)
    if fout===nothing
        fout_ = function(t,rho)
            tmp.data[:] = rho.data
            return copy(tmp)
        end
    else
        fout_ = function(t,rho)
            tmp.data[:] = rho.data
            return fout(t,tmp)
        end
    end

    # Solve
    return integrate_master(tspan, dmaster_, rho_, fout_; kwargs...)
end

"""
    timeevolution.master_dynamic(tspan, rho0, f; <keyword arguments>)

Time-evolution according to a master equation with a dynamic non-hermitian Hamiltonian and J.

In this case the given Hamiltonian is assumed to be the non-hermitian version.
```math
H_{nh} = H - \\frac{i}{2} \\sum_k J^†_k J_k
```
The given function can either be of the form `f(t, rho) -> (Hnh, Hnhdagger, J, Jdagger)`
or `f(t, rho) -> (Hnh, Hnhdagger, J, Jdagger, rates)` For further information look
at [`master_dynamic`](@ref).
"""
function master_nh_dynamic(tspan, rho0::Operator, f;
                rates=nothing,
                fout=nothing,
                kwargs...)
    tmp = copy(rho0)
    dmaster_(t, rho, drho) = dmaster_nh_dynamic!(drho, f, rates, rho, tmp, t)
    integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
end

"""
    timeevolution.master_dynamic(tspan, rho0, f; <keyword arguments>)

Time-evolution according to a master equation with a dynamic Hamiltonian and J.

There are two implementations for integrating the master equation with dynamic
operators:

* [`master_dynamic`](@ref): Usual formulation of the master equation.
* [`master_nh_dynamic`](@ref): Variant with non-hermitian Hamiltonian.

# Arguments
* `tspan`: Vector specifying the points of time for which output should be displayed.
* `rho0`: Initial density operator. Can also be a state vector which is
        automatically converted into a density operator.
* `f`: Function `f(t, rho) -> (H, J, Jdagger)` or `f(t, rho) -> (H, J, Jdagger, rates)`
* `rates=nothing`: Vector or matrix specifying the coefficients (decay rates)
        for the jump operators. If nothing is specified all rates are assumed
        to be 1.
* `fout=nothing`: If given, this function `fout(t, rho)` is called every time
        an output should be displayed. ATTENTION: The given state rho is not
        permanent! It is still in use by the ode solver and therefore must not
        be changed.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master_dynamic(tspan, rho0::Operator, f;
                rates=nothing,
                fout=nothing,
                kwargs...)
    tmp = copy(rho0)
    dmaster_(t, rho, drho) = dmaster_h_dynamic!(drho, f, rates, rho, tmp, t)
    integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
end

# Automatically convert Ket states to density operators
for f ∈ [:master,:master_h,:master_nh,:master_dynamic,:master_nh_dynamic]
    @eval $f(tspan,psi0::Ket,args...;kwargs...) = $f(tspan,dm(psi0),args...;kwargs...)
end

# Non-hermitian Hamiltonian
function nh_hamiltonian(H,J,Jdagger,::Nothing)
    Hnh = copy(H)
    for i=1:length(J)
        Hnh -= complex(float(eltype(H)))(0.5im)*Jdagger[i]*J[i]
    end
    return Hnh
end
function nh_hamiltonian(H,J,Jdagger,rates::AbstractVector)
    Hnh = copy(H)
    for i=1:length(J)
        Hnh -= complex(float(eltype(H)))(0.5im*rates[i])*Jdagger[i]*J[i]
    end
    return Hnh
end
function nh_hamiltonian(H,J,Jdagger,rates::AbstractMatrix)
    Hnh = copy(H)
    for i=1:length(J), j=1:length(J)
        Hnh -= complex(float(eltype(H)))(0.5im*rates[i,j])*Jdagger[i]*J[j]
    end
    return Hnh
end

# Recasting needed for the ODE solver is just providing the underlying data
function recast!(rho::Operator{B,B,T},x::T) where {B,T}
    rho.data = x
end
recast!(x::T,rho::Operator{B,B,T}) where {B,T} = nothing

function integrate_master(tspan, df, rho0, fout; kwargs...)
    x0 = rho0.data
    state = deepcopy(rho0)
    dstate = deepcopy(rho0)
    integrate(tspan, df, x0, state, dstate, fout; kwargs...)
end

# Time derivative functions
#   * dmaster_h
#   * dmaster_nh
#   * dmaster_h_dynamic -> callback(t, rho) -> dmaster_h
#   * dmaster_nh_dynamic -> callback(t, rho) -> dmaster_nh
# dmaster_h and dmaster_nh provide specialized implementations depending on
# the type of the given decay rate object which can either be nothing, a vector
# or a matrix.

"""
    dmaster_h!(drho, H, J, Jdagger, rates, rho, drho_cache)

Update `drho` according to a master equation given in standard Lindblad form.
A cached copy `drho_cache` of `drho` is used as a temporary saving step.

See also: [`master`](@ref), [`dmaster_nh!`](@ref), [`dmaster_h_dynamic!`](@ref),
    [`dmaster_nh_dynamic!`](@ref), [`dmaster_liouville!`](@ref)
"""
function dmaster_h!(drho, H, J, Jdagger, rates::Nothing, rho, drho_cache)
    QuantumOpticsBase.mul!(drho,H,rho,-eltype(rho)(im),zero(eltype(rho)))
    QuantumOpticsBase.mul!(drho,rho,H,eltype(rho)(im),one(eltype(rho)))
    for i=1:length(J)
        QuantumOpticsBase.mul!(drho_cache,J[i],rho)
        QuantumOpticsBase.mul!(drho,drho_cache,Jdagger[i],true,true)

        QuantumOpticsBase.mul!(drho,Jdagger[i],drho_cache,eltype(rho)(-0.5),one(eltype(rho)))

        QuantumOpticsBase.mul!(drho_cache,rho,Jdagger[i],true,false)
        QuantumOpticsBase.mul!(drho,drho_cache,J[i],eltype(rho)(-0.5),one(eltype(rho)))
    end
    return drho
end

function dmaster_h!(drho, H, J, Jdagger, rates::AbstractVector, rho, drho_cache)
    QuantumOpticsBase.mul!(drho,H,rho,-eltype(rho)(im),zero(eltype(rho)))
    QuantumOpticsBase.mul!(drho,rho,H,eltype(rho)(im),one(eltype(rho)))
    for i=1:length(J)
        QuantumOpticsBase.mul!(drho_cache,J[i],rho,eltype(rho)(rates[i]),zero(eltype(rho)))
        QuantumOpticsBase.mul!(drho,drho_cache,Jdagger[i],true,true)

        QuantumOpticsBase.mul!(drho,Jdagger[i],drho_cache,eltype(rho)(-0.5),one(eltype(rho)))

        QuantumOpticsBase.mul!(drho_cache,rho,Jdagger[i],eltype(rho)(rates[i]),zero(eltype(rho)))
        QuantumOpticsBase.mul!(drho,drho_cache,J[i],eltype(rho)(-0.5),one(eltype(rho)))
    end
    return drho
end

function dmaster_h!(drho, H, J, Jdagger, rates::AbstractMatrix, rho, drho_cache)
    QuantumOpticsBase.mul!(drho,H,rho,-eltype(rho)(im),zero(eltype(rho)))
    QuantumOpticsBase.mul!(drho,rho,H,eltype(rho)(im),one(eltype(rho)))
    for j=1:length(J), i=1:length(J)
        QuantumOpticsBase.mul!(drho_cache,J[i],rho,eltype(rho)(rates[i,j]),zero(eltype(rho)))
        QuantumOpticsBase.mul!(drho,drho_cache,Jdagger[j],true,true)

        QuantumOpticsBase.mul!(drho,Jdagger[j],drho_cache,eltype(rho)(-0.5),one(eltype(rho)))

        QuantumOpticsBase.mul!(drho_cache,rho,Jdagger[j],eltype(rho)(rates[i,j]),zero(eltype(rho)))
        QuantumOpticsBase.mul!(drho,drho_cache,J[i],eltype(rho)(-0.5),one(eltype(rho)))
    end
    return drho
end

"""
    dmaster_nh!(drho, Hnh, Hnh_dagger, J, Jdagger, rates, rho, drho_cache)

Updates `drho` according to a master equation given in standard Lindblad form.
The part of the Liuovillian which can be written as a commutator should be
contained in `Hnh` and `Hnh_dagger`. This allows to skip a number of matrix
multiplications making it slightly faster than [`dmaster_h!`](@ref).

See also: [`master`](@ref), [`dmaster_h!`](@ref), [`dmaster_h_dynamic!`](@ref),
    [`dmaster_nh_dynamic!`](@ref), [`dmaster_liouville!`](@ref)
"""
function dmaster_nh!(drho, Hnh, Hnh_dagger, J, Jdagger, rates::Nothing, rho, drho_cache)
    QuantumOpticsBase.mul!(drho,Hnh,rho,-eltype(rho)(im),zero(eltype(rho)))
    QuantumOpticsBase.mul!(drho,rho,Hnh_dagger,eltype(rho)(im),one(eltype(rho)))
    for i=1:length(J)
        QuantumOpticsBase.mul!(drho_cache,J[i],rho)
        QuantumOpticsBase.mul!(drho,drho_cache,Jdagger[i],true,true)
    end
    return drho
end

function dmaster_nh!(drho, Hnh, Hnh_dagger, J, Jdagger, rates::AbstractVector, rho, drho_cache)
    QuantumOpticsBase.mul!(drho,Hnh,rho,-eltype(rho)(im),zero(eltype(rho)))
    QuantumOpticsBase.mul!(drho,rho,Hnh_dagger,eltype(rho)(im),one(eltype(rho)))
    for i=1:length(J)
        QuantumOpticsBase.mul!(drho_cache,J[i],rho,eltype(rho)(rates[i]),zero(eltype(rho)))
        QuantumOpticsBase.mul!(drho,drho_cache,Jdagger[i],true,true)
    end
    return drho
end

function dmaster_nh!(drho, Hnh, Hnh_dagger, J, Jdagger, rates::AbstractMatrix, rho, drho_cache)
    QuantumOpticsBase.mul!(drho,Hnh,rho,-eltype(rho)(im),zero(eltype(rho)))
    QuantumOpticsBase.mul!(drho,rho,Hnh_dagger,eltype(rho)(im),one(eltype(rho)))
    for j=1:length(J), i=1:length(J)
        QuantumOpticsBase.mul!(drho_cache,J[i],rho,eltype(rho)(rates[i,j]),zero(eltype(rho)))
        QuantumOpticsBase.mul!(drho,drho_cache,Jdagger[j],true,true)
    end
    return drho
end

"""
    dmaster_liouville!(drho,L,rho)

Update `drho` according to a master equation as `L*rho`, where `L` is an arbitrary
(super-)operator.

See also: [`master`](@ref), [`dmaster_h!`](@ref), [`dmaster_nh!`](@ref),
    [`dmaster_h_dynamic!`](@ref), [`dmaster_nh_dynamic!`](@ref)
"""
function dmaster_liouville!(drho,L,rho)
    mul!(drho,L,rho)
    return drho
end

"""
    dmaster_h_dynamic!(drho, f, rates, rho, drho_cache, t)

Computes the Hamiltonian and jump operators as `H,J,Jdagger=f(t,rho)` and
update `drho` according to a master equation. Optionally, rates can also be
returned from `f`.

See also: [`master_dynamic`](@ref), [`dmaster_h!`](@ref), [`dmaster_nh!`](@ref),
    [`dmaster_nh_dynamic!`](@ref)
"""
function dmaster_h_dynamic!(drho, f, rates, rho, drho_cache, t)
    result = f(t, rho)
    QO_CHECKS[] && @assert 3 <= length(result) <= 4
    if length(result) == 3
        H, J, Jdagger = result
        rates_ = rates
    else
        H, J, Jdagger, rates_ = result
    end
    QO_CHECKS[] && check_master(rho, H, J, Jdagger, rates_)
    dmaster_h!(drho, H, J, Jdagger, rates_, rho, drho_cache)
end

"""
    dmaster_nh_dynamic!(drho, f, rates, rho, drho_cache, t)

Computes the non-hermitian Hamiltonian and jump operators as
`Hnh,Hnh_dagger,J,Jdagger=f(t,rho)` and update `drho` according to a master
equation. Optionally, rates can also be returned from `f`.

See also: [`master_dynamic`](@ref), [`dmaster_h!`](@ref), [`dmaster_nh!`](@ref),
    [`dmaster_h_dynamic!`](@ref)
"""
function dmaster_nh_dynamic!(drho, f, rates, rho, drho_cache, t)
    result = f(t, rho)
    QO_CHECKS[] && @assert 4 <= length(result) <= 5
    if length(result) == 4
        Hnh, Hnh_dagger, J, Jdagger = result
        rates_ = rates
    else
        Hnh, Hnh_dagger, J, Jdagger, rates_ = result
    end
    QO_CHECKS[] && check_master(rho, Hnh, J, Jdagger, rates_)
    dmaster_nh!(drho, Hnh, Hnh_dagger, J, Jdagger, rates_, rho, drho_cache)
end

function check_master(rho0, H, J, Jdagger, rates)
    isreducible = true # test if all operators are sparse or dense
    if !(isa(H, DenseOpType) || isa(H, SparseOpType))
        isreducible = false
    end
    for j=J
        @assert isa(j, AbstractOperator)
        if !(isa(j, DenseOpType) || isa(j, SparseOpType))
            isreducible = false
        end
        check_samebases(rho0, j)
    end
    for j=Jdagger
        @assert isa(j, AbstractOperator)
        if !(isa(j, DenseOpType) || isa(j, SparseOpType))
            isreducible = false
        end
        check_samebases(rho0, j)
    end
    @assert length(J)==length(Jdagger)
    if isa(rates,AbstractMatrix)
        @assert size(rates, 1) == size(rates, 2) == length(J)
    elseif isa(rates,AbstractVector)
        @assert length(rates) == length(J)
    end
    isreducible
end


text:
using NBodySimulator, StaticArrays, LinearAlgebra, StochasticDiffEq
using Test

include("lennard_jones_test.jl")
include("electrostatics_test.jl")
include("gravitational_test.jl")
include("custom_potential_test.jl")
include("magnetostaic_test.jl")
include("thermostat_test.jl")
include("water_test.jl")

ext: jl


text:
"""
    psconvert(cmd0::String="", arg1=[]; kwargs...)

Place images or EPS files on maps.

Full option list at [`psconvert`](http://gmt.soest.hawaii.edu/doc/latest/psconvert.html)

Parameters
----------

- **A** : **adjust** : -- Str or Number --

    Adjust the BoundingBox and HiResBoundingBox to the minimum required by the image content.
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/psconvert.html#a)
- **C** : **gs_option** : -- Str or Array os strings --

    Specify a single, or an araay of, custom option that will be passed on to GhostScript as is.
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/psconvert.html#c)
- **D** : **out_dir*** : *output_dir** : -- Str --

    Sets an alternative output directory (which must exist) [Default is the same directory
    as the PS files].
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/psconvert.html#d)
- **E** : **dpi** : -- Number --

    Set raster resolution in dpi [default = 720 for PDF, 300 for others].
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/psconvert.html#e)
- **F** : **:out_name** : **output_name** : -- Str --

    Force the output file name.
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/psconvert.html#f)
- **G** : **ghost_path** : -- Bool or [] --

    Full path to your GhostScript executable.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/psconvert.html#g)
- **I** : **icc_gray** : -- Bool or [] --

    Enforce gray-shades by using ICC profiles.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/psconvert.html#i)
- **in_memory** : -- Bool or [] --

	Process a in memory PS file. No other input file should be provided.
	Currently works on Windows only.
- **L** : **list_file** : -- Str --

    The listfile is an ASCII file with the names of the PostScript files to be converted.
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/psconvert.html#l)
- **Q** : **anti_aliasing** : -- Str --

	Set the anti-aliasing options for graphics or text. Append the size of the subsample box
	(1, 2, or 4) [4]. This option is set by default.
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/psconvert.html#q)
- **S** : **gs_command** : -- Bool or [] --

    Print to standard error the GhostScript command after it has been executed.
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/psconvert.html#s)
- **T** : **format** : -- Str --

    b|e|E|f|F|j|g|G|m|s|t Sets the output format, where b = BMP, e = EPS, E = EPS with PageSize command,
    f = PDF, F = multi-page PDF, j = JPEG, g = PNG, G = transparent PNG (untouched regions are
	transparent), m = PPM,  and t = TIFF [default is JPEG].
	Alternatively, the format may be set with the *fmt* keyword, e.g. *fmt="png"*.
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/psconvert.html#t)
- **W** : **world_file** : -- Str --

    Write a ESRI type world file suitable to make (e.g) .tif files be recognized as geotiff by
    software that know how to do it.
    [`-W`](http://gmt.soest.hawaii.edu/doc/latest/psconvert.html#w)
- **kml** -- Str or [] --

    Create a minimalist KML file that allows loading the image in GoogleEarth.
    [`-W`](http://gmt.soest.hawaii.edu/doc/latest/psconvert.html#w)
- **Z** : **del_input_ps** : -- Bool or [] --

    Remove the input PostScript file(s) after the conversion.
    [`-Z`](http://gmt.soest.hawaii.edu/doc/latest/psconvert.html#z)
- $(GMT.opt_V)
"""
function psconvert(cmd0::String="", arg1=[]; kwargs...)

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

	if (!isempty(cmd0)) && isempty_(arg1)  arg1 = cmd0  end

	d = KW(kwargs)
	cmd = add_opt("", 'A', d, [:A :adjust])
	if (cmd == " -A")  cmd = cmd * "1p"  end		# If just -A default to -A1p
	cmd = add_opt(cmd, 'D', d, [:D :out_dir :output_dir])
	cmd = add_opt(cmd, 'E', d, [:E :dpi])
	cmd = add_opt(cmd, 'F', d, [:F :out_name :output_name])
	cmd = add_opt(cmd, 'G', d, [:G :ghost_path])
	cmd = add_opt(cmd, 'I', d, [:I :icc_gray])
	cmd = add_opt(cmd, 'L', d, [:L :list_file])
	cmd = add_opt(cmd, 'Q', d, [:Q :anti_aliasing])
	cmd = add_opt(cmd, 'S', d, [:S :gs_command])
	cmd = add_opt(cmd, 'T', d, [:T :format])
	cmd = add_opt(cmd, 'Z', d, [:Z :del_input_ps])
	cmd = parse_V_params(cmd, d)

	fmt = ""
	if (haskey(d, :fmt))
		fmt = isa(d[:fmt], Symbol) ? string(d[:fmt]) : d[:fmt]
	end

	if (!occursin("-T", cmd) && !isempty(fmt) && fmt != "ps")	# Must convert the FMT into a -T opt
		if (fmt == "pdf")      cmd = cmd * " -Tf"
		elseif (fmt == "eps")  cmd = cmd * " -Te"
		elseif (fmt == "png")  cmd = cmd * " -Tg"
		elseif (fmt == "jpg")  cmd = cmd * " -Tj"
		elseif (fmt == "tif")  cmd = cmd * " -Tt"
		end
	end

	for sym in [:C :gs_option]
		if (haskey(d, sym))
			if (isa(d[sym], String))
				cmd = cmd * " -C" * d[sym]
			elseif (isa(d[sym], Array{Any}))
				for k = 1:length(d[sym])
					cmd = cmd * " -C" * d[sym][k]
				end
			end
			break
		end
	end

	cmd = add_opt(cmd, 'W', d, [:W :world_file])
	if (haskey(d, :kml))
		cmd = cmd * " -W+k" * d[:kml]
	end

	if (haskey(d, :in_memory))
		cmd = cmd * " ="
		if (!isempty_(arg1))
			@warn("The IN_MEMORY option is imcompatible with passing an input file name.
			      Dropping this one.")
			arg1 = []
		end
	end

	# In case DATA holds a file name, copy it into cmd.
	#cmd, arg1, = read_data(d, cmd0, cmd, arg1)
	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)		# Find how data was transmitted

	if (isempty(cmd))          cmd = "-A1p -Tj -Qg4 -Qt4"  end 	# Means no options were used. Allowed case
	if (!occursin("-Q", cmd))  cmd = cmd * " -Qt4 -Qg4"    end	# We promised to have these as default

	is_struct_PS = false
	if (isa(arg1, GMTps))
		is_struct_PS = true
	end
	if (is_struct_PS)  return gmt("psconvert " * cmd, arg1)
	else               return gmt("psconvert " * cmd)
	end
end


text:
using Gadfly

plot(xgroup=["A", "A", "B", "B"],
     x=["X", "Y", "X", "Y"],
     y=rand(4),
     Geom.subplot_grid(Geom.bar))


text:
using Statistics
using StatsBase: sample
# using Revise
using EvoTrees

# prepare a dataset
features = rand(Int(1.25e6), 100)
# features = rand(100, 10)
X = features
Y = rand(size(X, 1))
𝑖 = collect(1:size(X,1))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace = false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1))+1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

#############################
# CPU - linear
#############################
params1 = EvoTreeRegressor(T=Float32,
    loss=:linear, metric=:none,
    nrounds=100,
    λ = 1.0, γ=0.1, η=0.1,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=0.5, nbins=64)

@time model, cache = init_evotree(params1, X_train, Y_train);
# model = fit_evotree(params1, X_train, Y_train);
@time model = fit_evotree(params1, X_train, Y_train);
@time model = fit_evotree(params1, X_train, Y_train);


text:
"""a is a vector of length M, b is a matrix shape M x N
"""
function encoding_distance(a, b)
    @tullio distances[n] := a[m] - b[m,n]
end

function link_nodes_by_Δ!(sg::SpectroscopicGraph; Δ=1)
    encodings = reduce(hcat, sg.levels.encoding)
    for node ∈ 1:length(sg.levels)
        distances = abs.(encoding_distance(encodings[:,node], encodings))
        for (i, d) in enumerate(distances)
            d != Δ ? continue : MetaGraphs.add_edge!(sg.graph, node, i, Dict(:Δ => d))
        end
    end
end


text:
module RationalRoots

export
    # Abstract types
    RationalRoot,

    # Functions
    signedroot,
    signedsquare

const IntOrRational{T} = Union{T,Rational{T}} where T<:Integer

"""
    RationalRoot{T} where {T<:Integer} <: AbstractIrrational

Type for representing the positive or negative square root of a positive `Rational{T}`.
"""
struct RationalRoot{T<:Integer} <: AbstractIrrational
    signedsquare::Rational{T}
    # Inner constructor that is only used to define signedroot(::Type{RationalRoot{T}}, x)
    RationalRoot{T}(::Val{:inner}, signedsquare) where T<:Integer = new(signedsquare)
end
RationalRoot{T}(x::Real) where T<:Integer = signedroot(RationalRoot{T}, signedsquare(x))
RationalRoot{T}(x::RationalRoot{T}) where T<:Integer = x

"""
    signedroot([R<:RationalRoot,] x)

Return `sign(x)*sqrt(abs(x)) == x/sqrt(abs(x))`. If the first argument is not present, this value is given as a floating point number if `x isa AbstractFloat`, and as a `RationalRoot` if `x isa Union{Integer,Rational}`. With the first argument present, the return type is specified as `RationalRoot` or a specific `RationalRoot{T}`, and for floating point numbers `x` the result will first be rationalized by `Base.rationalize{T}`.

# Examples

```jldoctest
julia> signedroot(3)
+√(3//1)

julia> signedroot(-4.0)
-2.0

julia> signedroot(RationalRoot, 5.2)
+√(26//5)

julia> signedroot(RationalRoot{Int8}, 8)
+√(8//1)
```
"""
signedroot(x::Real) = sign(x)*sqrt(abs(x))
signedroot(x::IntOrRational{T}) where {T<:Integer} = signedroot(RationalRoot{T}, x)
signedroot(::Type{RationalRoot}, x) = signedroot(RationalRoot, rationalize(x))
signedroot(::Type{RationalRoot}, x::T) where T<:Integer = signedroot(RationalRoot{T}, x)
signedroot(::Type{RationalRoot}, x::Rational{T}) where T<:Integer =
    signedroot(RationalRoot{T}, x)
signedroot(::Type{RationalRoot{T}}, x::AbstractFloat) where T<:Integer =
    signedroot(RationalRoot{T}, rationalize(T, x))
signedroot(::Type{RationalRoot{T}}, x) where T<:Integer = RationalRoot{T}(Val{:inner}(), x)

"""
    signedsquare(x)

Return `sign(x)*x*x = x*abs(x)`. If `x` is a `RationalRoot`, return an appropriate rational type.

```jldoctest
julia> signedsquare(-2)
-4

julia> signedsquare(4.0)
16.0

julia> signedsquare(-RationalRoot{Int}(3//2))
-9//4
```
"""
signedsquare(x) = x*abs(x)
signedsquare(x::RationalRoot) = x.signedsquare

Base.promote_rule(::Type{RationalRoot{T1}}, ::Type{T2}) where
    {T1<:Integer, T2<:Integer} = RationalRoot{promote_type(T1, T2)}
Base.promote_rule(::Type{RationalRoot{T1}}, ::Type{Rational{T2}}) where
    {T1<:Integer, T2<:Integer} = RationalRoot{promote_type(T1, T2)}
Base.promote_rule(::Type{RationalRoot{T1}}, ::Type{RationalRoot{T2}}) where
    {T1<:Integer, T2<:Integer} = RationalRoot{promote_type(T1, T2)}

RationalRoot(x::RationalRoot) = x
RationalRoot(x::Number) = signedroot(RationalRoot, signedsquare(x))

function Base.convert(T::Type{<:AbstractFloat}, x::RationalRoot)
    s = _convert(T, signedsquare(x))
    s < zero(s) ? -sqrt(-s) : sqrt(s)
end
function _convert(T::Type{<:Union{Float32,Float64}}, x::Rational{BigInt})
    n, d = numerator(x), denominator(x)
    if typemin(Int) <= n <= typemax(Int) && typemin(Int) <= d <= typemax(Int)
        # fast path, don't go via BigFloat
        convert(T, convert(Int, n)//convert(Int, d))
    else
        convert(T, x)
    end
end
_convert(T, x) = convert(T, x)

function Base.convert(T::Type{<:Integer}, x::RationalRoot)
    s = convert(Rational{T}, signedsquare(x))
    if denominator(s) == one(T)
        a = findsignedroot(numerator(s))
        if a !== nothing
            return a
        end
    end
    throw(InexactError(nameof(T), T, x))
end
function Base.convert(T::Type{<:Rational}, x::RationalRoot)
    s = convert(T, signedsquare(x))
    a = findsignedroot(s)
    if a !== nothing
        return a
    else
        throw(InexactError(nameof(T), T, x))
    end
end

Base.AbstractFloat(x::RationalRoot) = convert(AbstractFloat, x)
Base.Float32(x::RationalRoot) = convert(Float32, x)
Base.Float64(x::RationalRoot) = convert(Float64, x)
Base.BigFloat(x::RationalRoot) = convert(BigFloat, x)

function Base.hash(a::RationalRoot, h::UInt)
    s = signedsquare(a)
    x = findsignedroot(s)
    if x !== nothing
        return hash(x, h)
    else
        return hash(s, hash(s, h)) # something arbitrary
    end
end

for op in (:<, :≤, :(==))
    @eval Base.$op(x::RationalRoot, y::RationalRoot) = $op(signedsquare(x), signedsquare(y))

    for T in (Rational, Rational{BigInt}, Int8, Int16, Int32, Int64, BigInt, AbstractFloat, Float32, Float64)

        @eval Base.$op(x::RationalRoot, y::$T) = $op(signedsquare(x), signedsquare(y))
        @eval Base.$op(x::$T, y::RationalRoot) = $op(signedsquare(x), signedsquare(y))
    end
end

Base.:+(x::RationalRoot) = signedroot(+signedsquare(x))
Base.:-(x::RationalRoot) = signedroot(-signedsquare(x))
for op in (:*, :/, :\, ://)
    @eval Base.$op(x::RationalRoot, y::RationalRoot) =
        signedroot($op(signedsquare(x), signedsquare(y)))
    @eval Base.$op(x::RationalRoot, y::IntOrRational) =
        signedroot($op(signedsquare(x), signedsquare(y)))
    @eval Base.$op(x::IntOrRational, y::RationalRoot) =
        signedroot($op(signedsquare(x), signedsquare(y)))
end

Base.inv(x::RationalRoot) = signedroot(inv(signedsquare(x)))

Base.one(x::RationalRoot) = one(typeof(x))
Base.zero(x::RationalRoot) = zero(typeof(x))
Base.one(::Type{RationalRoot{T}}) where T<:Integer = signedroot(one(T))
Base.zero(::Type{RationalRoot{T}}) where T<:Integer = signedroot(zero(T))
Base.isone(x::RationalRoot) = isone(signedsquare(x))
Base.iszero(x::RationalRoot) = iszero(signedsquare(x))

Base.sign(x::RationalRoot) = sign(signedsquare(x))*one(x)

Base.big(x::RationalRoot) = convert(RationalRoot{BigInt}, x)

Base.widen(::Type{RationalRoot{T}}) where T<:Integer = RationalRoot{widen(T)}

Base.typemax(::Type{RationalRoot{T}}) where T<:Integer = signedroot(typemax(Rational{T}))
Base.typemin(::Type{RationalRoot{T}}) where T<:Integer = signedroot(typemin(Rational{T}))

function Base.show(io::IO, x::RationalRoot)
    signedsquare(x) < 0 ? print(io, "-√(") : print(io, "+√(")
    show(io, abs(signedsquare(x)))
    print(io, ")")
end

function findsignedroot(x::Integer)
    n = abs(x)
    k = isqrt(n)
    if k*k == n
        return sign(x)*k
    else
        return nothing
    end
end
function findsignedroot(x::Rational)
    n = findsignedroot(numerator(x))
    d = findsignedroot(denominator(x))
    if n === nothing || d === nothing
        return nothing
    else
        return n//d
    end
end

function Base.isinteger(x::RationalRoot)
    s = signedsquare(x)
    return isone(denominator(s)) && findsignedroot(numerator(s)) !== nothing
end

end


text:
using Filthy
using Test
using StaticArrays
using BenchmarkTools
using CSV

#=
Nile data
==========================#

nile = CSV.read("Nile.csv", allowmissing=:none)

nile = Matrix(nile)

# Values from Durbin and Koopman
# Zn = fill(1.0, (1,1))
# Hn = fill(15099., (1,1))
# Rn = fill(1.0, (1,1))
# Tn = fill(1.0, (1,1))
# Qn = fill(1469.1, (1,1))

# Values from MLE (inexact, P = 1e07)
Zn = fill(1.0, (1,1))
Hn = fill(15099.82, (1,1))
Rn = fill(1.0, (1,1))
Tn = fill(1.0, (1,1))
Qn = fill(1468.487, (1,1))

Zs = SMatrix{1,1}(Zn)
Hs = SMatrix{1,1}(Hn)
Rs = SMatrix{1,1}(Rn)
Ts = SMatrix{1,1}(Tn)
Qs = SMatrix{1,1}(Qn)

#=
Constructors
=----------=#
@testset "Constructors (inexact).........." begin
    @test isa(LinearStateSpace(1.0, 15099.82, 1.0, 1.0, 1469.1, 0.0, 10.0^7), LinearStateSpace)
    @test isa(LinearStateSpace(Zn, Hn, Tn, Rn, Qn, [0.0], fill(10.0^7, (1,1))), LinearStateSpace)
    @test isa(LinearStateSpace(Zs, Hs, Ts, Rs, Qs, SVector{1}(0.0), SMatrix{1,1}(fill(10.0^7, (1,1)))), LinearStateSpace)

    @test_throws AssertionError LinearStateSpace(Zn, Hn, Tn, Rn, Qn, 0.0, fill(10.0^7, (1,1)))
    @test_throws AssertionError LinearStateSpace(Zn, Hn, Tn, Rn, Qn, [0.0], 10.0^7)
    @test_throws AssertionError LinearStateSpace(Zn, Hn, Tn, Rn, Qn, 0.0, 10.0^7)

    @test_throws AssertionError LinearStateSpace(Zs, Hs, Ts, Rs, Qs, SMatrix{1,1}([0.0]), rand(1,1))
    @test_throws AssertionError LinearStateSpace(Zs, Hs, Ts, Rs, Qs, SMatrix{1,1}([0.0]), SMatrix{1,1}(fill(10.0^7, (1,1))))
    @test_throws AssertionError LinearStateSpace(Zs, Hs, Ts, Rs, Qs, SVector{1}(0.0), rand(1,1))

    @test_throws AssertionError LinearStateSpace(1.0, 15099.82, 1.0, 1.0, 1469.1, [0.0], 10.0^7)
    @test_throws AssertionError LinearStateSpace(1.0, 15099.82, 1.0, 1.0, 1469.1, 0.0, [10.0^7])
    @test_throws AssertionError LinearStateSpace(1.0, 15099.82, 1.0, 1.0, 1469.1, [0.0], [10.0^7])
end

@testset "Constructors (exact).........." begin
    @test isa(LinearStateSpace(1.0, 15099.82, 1.0, 1.0, 1468.487, 0.0, 10.0^7, P1inf = 1.0), LinearStateSpace)
    @test isa(LinearStateSpace(Zn, Hn, Tn, Rn, Qn, [0.0], fill(10.0^7, (1,1)), P1inf = fill(1.0, (1,1))), LinearStateSpace)
    @test isa(LinearStateSpace(Zs, Hs, Ts, Rs, Qs, SVector{1}(0.0), SMatrix{1,1}(fill(10.0^7, (1,1)))), LinearStateSpace)

    ss = LinearStateSpace(Zs, Hs, Ts, Rs, Qs, SVector{1}(0.0), SMatrix{1,1}(fill(10.0^7, (1,1))), P1inf = SMatrix{1,1}([1.0]))
    @test isa(ss.f.Pinf[1], StaticArrays.SArray{Tuple{1,1},Float64,2,1})
    @test isa(ss.f.Pstar[1], StaticArrays.SArray{Tuple{1,1},Float64,2,1})
    @test ss.f.Pinf[1][1] == 1.0
    @test ss.f.Pstar[1][1] == 0.0

    ss = LinearStateSpace(Zn, Hn, Tn, Rn, Qn, [0.0], fill(10.0^7, (1,1)), P1inf = fill(1.0, (1,1)))
    @test isa(ss.f.Pinf[1], Array{Float64, 2})
    @test isa(ss.f.Pstar[1], Array{Float64, 2})
    @test ss.f.Pinf[1][1] == 1.0
    @test ss.f.Pstar[1][1] == 0.0

    ss = LinearStateSpace(1.0, 15099.82, 1.0, 1.0, 1468.487, 0.0, 10.0^7, Pinf = 1.0)
    @test isa(ss.f.Pinf[1], Float64)
    @test isa(ss.f.Pstar[1], Float64)
    @test ss.f.Pinf[1][1] == 1.0
    @test ss.f.Pstar[1][1] == 0.0

    # @test_throws AssertionError LinearStateSpace(Zn, Hn, Tn, Rn, Qn, 0.0, fill(10.0^7, (1,1)))
    # @test_throws AssertionError LinearStateSpace(Zn, Hn, Tn, Rn, Qn, [0.0], 10.0^7)
    # @test_throws AssertionError LinearStateSpace(Zn, Hn, Tn, Rn, Qn, 0.0, 10.0^7)

    # @test_throws AssertionError LinearStateSpace(Zs, Hs, Ts, Rs, Qs, SMatrix{1,1}([0.0]), rand(1,1))
    # @test_throws AssertionError LinearStateSpace(Zs, Hs, Ts, Rs, Qs, SMatrix{1,1}([0.0]), SMatrix{1,1}(fill(10.0^7, (1,1))))
    # @test_throws AssertionError LinearStateSpace(Zs, Hs, Ts, Rs, Qs, SVector{1}(0.0), rand(1,1))

    # @test_throws AssertionError LinearStateSpace(1.0, 15099.82, 1.0, 1.0, 1468.487, [0.0], 10.^7.)
    # @test_throws AssertionError LinearStateSpace(1.0, 15099.82, 1.0, 1.0, 1468.487, 0.0, [10.^7.])
    # @test_throws AssertionError LinearStateSpace(1.0, 15099.82, 1.0, 1.0, 1468.487, [0.0], [10.^7.])
    # @test_throw  AssertionError LinearStateSpace(Zn, Hn, Tn, Rn, Qn, [0.0], fill(10.^7., (1,1)))
end

#=
Online version
=#

ss_scalar = LinearStateSpace(1.0, 15099.82, 1.0, 1.0, 1468.487, 0.0, 10.0^7)
for j in eachindex(nile)
    Filthy.onlinefilter!(ss_scalar, nile[j])
end

ss_matrix = LinearStateSpace(Zn, Hn, Tn, Rn, Qn, [0.0], fill(10.0^7, (1,1)))
for j in eachindex(nile)
    Filthy.onlinefilter!(ss_matrix, nile[j,:])
end

ss_static = LinearStateSpace(Zs, Hs, Ts, Rs, Qs, SVector{1}(0.0), SMatrix{1,1}(fill(10.^7., (1,1))))
for j in eachindex(nile)
    Filthy.onlinefilter!(ss_static, nile[j,:])
end

@testset "Filter correctness (online)....." begin
    @test isa(ss_matrix.f.a[1], Vector)
    @test isa(ss_matrix.f.P[1], Matrix)

    @test isa(ss_scalar.f.a[1], Float64)
    @test isa(ss_scalar.f.P[1], Float64)

    @test isa(ss_static.f.a[1], StaticArray)
    @test isa(ss_static.f.P[1], StaticMatrix)

    @test ss_scalar.f.a[1] == ss_matrix.f.a[1][1]
    @test ss_scalar.f.a[10] == ss_matrix.f.a[10][1]
    @test ss_scalar.f.a[101] == ss_matrix.f.a[101][1]

    @test ss_scalar.f.a[1] == ss_static.f.a[1][1]
    @test ss_scalar.f.a[10] == ss_static.f.a[10][1]
    @test ss_scalar.f.a[101] == ss_static.f.a[101][1]

    @test ss_scalar.c.F[1]   == ss_static.c.F[1][1]
    @test ss_scalar.c.F[10]  == ss_static.c.F[10][1]
    @test ss_scalar.c.F[101] == ss_static.c.F[101][1]

    @test ss_scalar.f.a[1] == 0.0
    @test ss_scalar.f.a[10] ≈ 1171.22228
    @test ss_scalar.f.a[101] ≈ 798.3871578

    @test ss_scalar.c.F[2] ≈ 10015099.822227377445
    @test ss_scalar.c.F[11] ≈ 20635.5495
    @test ss_scalar.c.F[101] ≈ 20599.87957

end

#=
Vectorized version
=#

ss_scalar = LinearStateSpace(1.0, 15099.82, 1.0, 1.0, 1468.487, 0.0, 10.0^7.)
Filthy.filter!(ss_scalar, vec(nile))

ss_matrix = LinearStateSpace(Zn, Hn, Tn, Rn, Qn, [0.0], fill(10.0^7, (1,1)))
Filthy.filter!(ss_matrix, nile')

ss_static = LinearStateSpace(Zs, Hs, Ts, Rs, Qs, SVector{1}(0.0), SMatrix{1,1}(fill(10.^7., (1,1))))
Filthy.filter!(ss_static, nile')

@testset "Filter correctness (offline) [inexact]..." begin
    @test ss_scalar_vec.f.a[1] == ss_scalar.f.a[1]
    @test ss_scalar_vec.f.a[10] == ss_scalar.f.a[10]
    @test ss_scalar_vec.f.a[100] == ss_scalar.f.a[100]
    @test ss_matrix_vec.f.a[1] == ss_matrix.f.a[1]
    @test ss_matrix_vec.f.a[10] == ss_matrix.f.a[10]
    @test ss_matrix_vec.f.a[100] == ss_matrix.f.a[100]
    @test ss_static_vec.f.a[1] == ss_static.f.a[1]
    @test ss_static_vec.f.a[10] == ss_static.f.a[10]
    @test ss_static_vec.f.a[100] == ss_static.f.a[100]
    @test ss_scalar_vec.f.a[1] == 0.0
    @test ss_scalar_vec.f.a[10] ≈ 1171.22228
    @test ss_scalar_vec.f.a[101] ≈ 798.3871578
    @test ss_scalar_vec.c.F[2] ≈ 10015099.822227377445
    @test ss_scalar_vec.c.F[11] ≈ 20635.5495
    @test ss_scalar_vec.c.F[101] ≈ 20599.87957
end

# Values from MLE (inexact, P = 1)

Zn = fill(1.0, (1,1))
Hn = fill(15098.654334841132368, (1,1))
Rn = fill(1.0, (1,1))
Tn = fill(1.0, (1,1))
Qn = fill(1469.1632513366273542, (1,1))

Zs = SMatrix{1,1}(Zn)
Hs = SMatrix{1,1}(Hn)
Rs = SMatrix{1,1}(Rn)
Ts = SMatrix{1,1}(Tn)
Qs = SMatrix{1,1}(Qn)

ss = LinearStateSpace(Zs, Hs, Ts, Rs, Qs, SVector{1}(0.0), SMatrix{1,1}(fill(0., (1,1))), P1inf = SMatrix{1,1}([1.0]))
Filthy.filter!(ss, nile')

@testset "Filter correctness (offline) [exact]....." begin
    @test ss.f.Pinf[1][1] == 1.0
    @test ss.f.P[1][1] == 0.0
    @test ss.f.Pstar[1][1] == 0.0
    @test ss.f.Pstar[2][1] ≈ 16568.307
    @test ss.f.P[2][1][1] ≈ 16568.307
    @test ss.f.d==2
    @test ss.f.hasdiffended[1] == true
    @test_throws BoundsError ss.f.Pstar[ss.f.d+2]
    @test ss.loglik[1] ≈ -633.46456
end

#=
Smoothing
---------=#
Filthy.smooth!(ss)
@testset "State smoother" begin
    ss.s.a[1]
end

# @btime begin
#     ss_scalar = LinearStateSpace(1.0, 15099.82, 1.0, 1.0, 1468.487, 0.0, 10.^7.)
#     Filthy.filter!(ss_scalar, vec(nile))
# end

# # stat_online = @btime begin
# #     cf_static = Filthy.CovarianceFilter(Zs, Hs, Ts, Rs, Qs, SMatrix{1,1}([0.0]), SMatrix{1,1}(fill(10.^7., (1,1))))
# #     for j in eachindex(nile)
# #         filter!(cf_static, [nile[j]])
# #     end
# # end

# # stat_vec = @btime begin
# #     cf_static = Filthy.CovarianceFilter(Zs, Hs, Ts, Rs, Qs, SMatrix{1,1}([0.0]), SMatrix{1,1}(fill(10.^7., (1,1))))
# #     filter!(cf_static, nile)
# # end

# # @btime begin
# #     cf_matrix = Filthy.CovarianceFilter(Zn, Hn, Tn, Rn, Qn, [0.0], fill(10.^7., (1,1)))
# #     for j in eachindex(nile)
# #         filter!(cf_matrix, [nile[j]])
# #     end
# # end

# # @test all(cf_scalar.s.att[:,1] .≈ cf_matrix.s.att[:,1])
# # @test all(cf_static.s.att[:,1] .≈ cf_matrix.s.att[:,1])

# # using Plots
# # gr()
# # Plots.plot(1:length(nile), cf_static.f.a[2:end,1], color = :darkblue)
# # plot!(nile, color = :black)
# # plot!(cf_static.s.a, color = :red)
# # plot!(cf.s.att[2:end,1]-1.64*sqrt.(cf.s.Ptt[2:end,1]), color = :darkblue, linestyle = :dot)
# # plot!(cf.s.att[2:end,1]+1.64*sqrt.(cf.s.Ptt[2:end,1]), color = :darkblue, linestyle = :dot)

# p = 2
# m = 2
# r = 2

# Z = SMatrix{p,p}([0.2 0.1; 0.4 0.15])
# H = SMatrix{p,p}([1.2 0.1; 0.1 1.15])
# T = SMatrix{m,m}([0.92 0.2; -.2 -0.9])
# R = I
# Q = SMatrix{m,m}(eye(m))
# a1 = SVector{2}(zeros(2))
# P1 = SMatrix{2,2}(1000.*eye(2))
# Pinf = SMatrix{2,2}(eye(2))

# ss = Filthy.LinearStateSpace(Z, H, T, R, Q, a1, P1, Pinf = Pinf)
# Filthy.filter!(ss, Y)

# ss2 = Filthy.LinearStateSpace(Z, H, T, R, Q, a1, P1)
# Filthy.filter!(ss2, Y)

# srand(12234)
# Y, a = simulate(ss, 200)

#writecsv("Y.csv", Y')

# # for i in 1:200
# #     filter!(cf, Y[:, i])
# # end

# # @test all(convert(Vector, cf.f.a[10]) .≈ [0.46446643122638098244; 2.13857431622820914896])
# # @test all(convert(Vector, cf.f.a[100]) .≈ [1.20622036437758350935; -0.39600752142009676415])

# # cf = CovarianceFilter(Zp, Hp, Tp, Rp, Qp, SVector{2}(zeros(2)), SMatrix{2,2}(eye(2)))
# # filter!(cf, Y')
# # @test all(convert(Vector, cf.f.a[10]) .≈ [0.46446643122638098244; 2.13857431622820914896])
# # @test all(convert(Vector, cf.f.a[100]) .≈ [1.20622036437758350935; -0.39600752142009676415])

# # tt = @btime begin
# #     cf = CovarianceFilter(Zp, Hp, Tp, Rp, Qp, SVector{2}(zeros(2)), SMatrix{2,2}(eye(2)))
# #     for i in 1:size(Y,2)
# #         filter!(cf, Y[:, i])
# #     end
# # end

# # vv = @btime begin
# #     cf = CovarianceFilter(Zp, Hp, Tp, Rp, Qp, SVector{2}(zeros(2)), SMatrix{2,2}(eye(2)))
# #     filter!(cf, Y')
# # end

# # # using Plots
# # # gr()
# # # p1 = Plots.plot(1:size(Y,2), cf.s.att[2:end,1], color = :darkblue)
# # # p2 = Plots.plot(1:size(Y,2), cf.s.att[2:end,2], color = :darkblue)

# # p1 = Plots.plot(1:size(Y,2), cf.f.a[2:end,1], color = :darkblue)
# # p2 = Plots.plot(1:size(Y,2), cf.f.a[2:end,2], color = :darkblue)

# # # plot!(p1, a[1,:])
# # # plot!(p2, a[2,:])
# # # plot!(nile, color = :black)
# # # plot!(cf.s.att[2:end,1]-1.64*sqrt.(cf.s.Ptt[2:end,1]), color = :darkblue, linestyle = :dot)
# # # plot!(cf.s.att[2:end,1]+1.64*sqrt.(cf.s.Ptt[2:end,1]), color = :darkblue, linestyle = :dot)

# # # using Plots
# # # gr()
# # # Plots.plot(1:length(nile), cf.s.att[2:end,1], color = :darkblue)
# # # plot!(nile, color = :black)
# # # plot!(cf.s.att[2:end,1]-1.64*sqrt.(cf.s.Ptt[2:end,1]), color = :darkblue, linestyle = :dot)
# # # plot!(cf.s.att[2:end,1]+1.64*sqrt.(cf.s.Ptt[2:end,1]), color = :darkblue, linestyle = :dot)

# # # Plots.plot(1:length(nile), cf.s.Ptt[2:end,1], color = :darkblue)
# # # Plots.

# # @btime begin
# #     cf = Filthy.CovarianceFilter(Z, H, T, R, Q, zeros(2), eye(2))
# #     for i in 1:200
# #         filter!(cf, Y[:, i])
# #     end
# # end

# # using Plots
# # gr()

# # Plots.plot(1:200, a[1,:])
# # plot!(kf.s.att[1:end,1])

# ## Nile Data

# ## In the Nile data example,

# ## y = α + ϵ
# ## α' = α + η

# ## In terms of Matrices

# ## Y    = Zt αt + ϵt, ϵt ∼ N(0, Ht)
# ## αt+1 = Tt αt + Rt ηt, ηt ~ N(0, Qt)

# # scalar_online = @btime begin
# #     cf_scalar = LinearStateSpace(1.0, 15099., 1.0, 1.0, 1469.1, 0.0, 10.^7.)
# #     for j in eachindex(nile)
# #         Filthy.onlinefilter!(cf_scalar, nile[j])
# #     end
# # end

# # stat_online = @btime begin
# #     cf_static = LinearStateSpace(Zs, Hs, Ts, Rs, Qs, SMatrix{1,1}([0.0]), SMatrix{1,1}(fill(10.^7., (1,1))))
# #     for j in eachindex(nile)
# #         Filthy.onlinefilter!(cf_static, nile[j,:])
# #     end
# # end

# # matrix_online = @btime begin
# #     ss_matrix = LinearStateSpace(Zn, Hn, Tn, Rn, Qn, [0.0], fill(10.^7., (1,1)))
# #     for j in eachindex(nile)
# #         Filthy.onlinefilter!(ss_matrix, nile[j,:])
# #     end
# # end


text:
@testset "4.5.10 (c+d x)^m (a+b sec)^n" begin
    (a, b, c, d, e, f, m, n, x, ) = @variables a b c d e f m n x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*sec(e+f*x))^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*sec(e+f*x))^n*with*a^2-b^2=0=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+a*sec(e+f*x))^n=#

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [(a + a*sec(e + f*x))*(c + d*x)^3, x, 11, (a*(c + d*x)^4)/(4*d) - (2*im*a*(c + d*x)^3*atan(ℯ^(im*(e + f*x))))/f + (3*im*a*d*(c + d*x)^2*PolyLog(2, (-im)*ℯ^(im*(e + f*x))))/f^2 - (3*im*a*d*(c + d*x)^2*PolyLog(2, im*ℯ^(im*(e + f*x))))/f^2 - (6*a*d^2*(c + d*x)*PolyLog(3, (-im)*ℯ^(im*(e + f*x))))/f^3 + (6*a*d^2*(c + d*x)*PolyLog(3, im*ℯ^(im*(e + f*x))))/f^3 - (6*im*a*d^3*PolyLog(4, (-im)*ℯ^(im*(e + f*x))))/f^4 + (6*im*a*d^3*PolyLog(4, im*ℯ^(im*(e + f*x))))/f^4]
    @test_int [(a + a*sec(e + f*x))*(c + d*x)^2, x, 9, (a*(c + d*x)^3)/(3*d) - (2*im*a*(c + d*x)^2*atan(ℯ^(im*(e + f*x))))/f + (2*im*a*d*(c + d*x)*PolyLog(2, (-im)*ℯ^(im*(e + f*x))))/f^2 - (2*im*a*d*(c + d*x)*PolyLog(2, im*ℯ^(im*(e + f*x))))/f^2 - (2*a*d^2*PolyLog(3, (-im)*ℯ^(im*(e + f*x))))/f^3 + (2*a*d^2*PolyLog(3, im*ℯ^(im*(e + f*x))))/f^3]
    @test_int [(a + a*sec(e + f*x))*(c + d*x)^1, x, 7, (a*(c + d*x)^2)/(2*d) - (2*im*a*(c + d*x)*atan(ℯ^(im*(e + f*x))))/f + (im*a*d*PolyLog(2, (-im)*ℯ^(im*(e + f*x))))/f^2 - (im*a*d*PolyLog(2, im*ℯ^(im*(e + f*x))))/f^2]
    @test_int [(a + a*sec(e + f*x))/(c + d*x)^1, x, 0, Unintegrable((a + a*sec(e + f*x))/(c + d*x), x)]
    @test_int [(a + a*sec(e + f*x))/(c + d*x)^2, x, 0, Unintegrable((a + a*sec(e + f*x))/(c + d*x)^2, x)]

    @test_int [(a + a*sec(e + f*x))^2*(c + d*x)^3, x, 17, -((im*a^2*(c + d*x)^3)/f) + (a^2*(c + d*x)^4)/(4*d) - (4*im*a^2*(c + d*x)^3*atan(ℯ^(im*(e + f*x))))/f + (3*a^2*d*(c + d*x)^2*log(1 + ℯ^(2*im*(e + f*x))))/f^2 + (6*im*a^2*d*(c + d*x)^2*PolyLog(2, (-im)*ℯ^(im*(e + f*x))))/f^2 - (6*im*a^2*d*(c + d*x)^2*PolyLog(2, im*ℯ^(im*(e + f*x))))/f^2 - (3*im*a^2*d^2*(c + d*x)*PolyLog(2, -ℯ^(2*im*(e + f*x))))/f^3 - (12*a^2*d^2*(c + d*x)*PolyLog(3, (-im)*ℯ^(im*(e + f*x))))/f^3 + (12*a^2*d^2*(c + d*x)*PolyLog(3, im*ℯ^(im*(e + f*x))))/f^3 + (3*a^2*d^3*PolyLog(3, -ℯ^(2*im*(e + f*x))))/(2*f^4) - (12*im*a^2*d^3*PolyLog(4, (-im)*ℯ^(im*(e + f*x))))/f^4 + (12*im*a^2*d^3*PolyLog(4, im*ℯ^(im*(e + f*x))))/f^4 + (a^2*(c + d*x)^3*tan(e + f*x))/f]
    @test_int [(a + a*sec(e + f*x))^2*(c + d*x)^2, x, 14, -((im*a^2*(c + d*x)^2)/f) + (a^2*(c + d*x)^3)/(3*d) - (4*im*a^2*(c + d*x)^2*atan(ℯ^(im*(e + f*x))))/f + (2*a^2*d*(c + d*x)*log(1 + ℯ^(2*im*(e + f*x))))/f^2 + (4*im*a^2*d*(c + d*x)*PolyLog(2, (-im)*ℯ^(im*(e + f*x))))/f^2 - (4*im*a^2*d*(c + d*x)*PolyLog(2, im*ℯ^(im*(e + f*x))))/f^2 - (im*a^2*d^2*PolyLog(2, -ℯ^(2*im*(e + f*x))))/f^3 - (4*a^2*d^2*PolyLog(3, (-im)*ℯ^(im*(e + f*x))))/f^3 + (4*a^2*d^2*PolyLog(3, im*ℯ^(im*(e + f*x))))/f^3 + (a^2*(c + d*x)^2*tan(e + f*x))/f]
    @test_int [(a + a*sec(e + f*x))^2*(c + d*x)^1, x, 9, (a^2*(c + d*x)^2)/(2*d) - (4*im*a^2*(c + d*x)*atan(ℯ^(im*(e + f*x))))/f + (a^2*d*log(cos(e + f*x)))/f^2 + (2*im*a^2*d*PolyLog(2, (-im)*ℯ^(im*(e + f*x))))/f^2 - (2*im*a^2*d*PolyLog(2, im*ℯ^(im*(e + f*x))))/f^2 + (a^2*(c + d*x)*tan(e + f*x))/f]
    @test_int [(a + a*sec(e + f*x))^2/(c + d*x)^1, x, 0, Unintegrable((a + a*sec(e + f*x))^2/(c + d*x), x)]
    @test_int [(a + a*sec(e + f*x))^2/(c + d*x)^2, x, 0, Unintegrable((a + a*sec(e + f*x))^2/(c + d*x)^2, x)]

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [1/(a + a*sec(e + f*x))*(c + d*x)^3, x, 9, (im*(c + d*x)^3)/(a*f) + (c + d*x)^4/(4*a*d) - (6*d*(c + d*x)^2*log(1 + ℯ^(im*(e + f*x))))/(a*f^2) + (12*im*d^2*(c + d*x)*PolyLog(2, -ℯ^(im*(e + f*x))))/(a*f^3) - (12*d^3*PolyLog(3, -ℯ^(im*(e + f*x))))/(a*f^4) - ((c + d*x)^3*tan(e/2 + (f*x)/2))/(a*f)]
    @test_int [1/(a + a*sec(e + f*x))*(c + d*x)^2, x, 8, (im*(c + d*x)^2)/(a*f) + (c + d*x)^3/(3*a*d) - (4*d*(c + d*x)*log(1 + ℯ^(im*(e + f*x))))/(a*f^2) + (4*im*d^2*PolyLog(2, -ℯ^(im*(e + f*x))))/(a*f^3) - ((c + d*x)^2*tan(e/2 + (f*x)/2))/(a*f)]
    @test_int [1/(a + a*sec(e + f*x))*(c + d*x)^1, x, 5, (c + d*x)^2/(2*a*d) - (2*d*log(cos(e/2 + (f*x)/2)))/(a*f^2) - ((c + d*x)*tan(e/2 + (f*x)/2))/(a*f)]
    @test_int [1/(a + a*sec(e + f*x))/(c + d*x)^1, x, 0, Unintegrable(1/((c + d*x)*(a + a*sec(e + f*x))), x)]
    @test_int [1/(a + a*sec(e + f*x))/(c + d*x)^2, x, 0, Unintegrable(1/((c + d*x)^2*(a + a*sec(e + f*x))), x)]

    @test_int [1/(a + a*sec(e + f*x))^2*(c + d*x)^3, x, 19, (5*im*(c + d*x)^3)/(3*a^2*f) + (c + d*x)^4/(4*a^2*d) - (10*d*(c + d*x)^2*log(1 + ℯ^(im*(e + f*x))))/(a^2*f^2) + (4*d^3*log(cos(e/2 + (f*x)/2)))/(a^2*f^4) + (20*im*d^2*(c + d*x)*PolyLog(2, -ℯ^(im*(e + f*x))))/(a^2*f^3) - (20*d^3*PolyLog(3, -ℯ^(im*(e + f*x))))/(a^2*f^4) - (d*(c + d*x)^2*sec(e/2 + (f*x)/2)^2)/(2*a^2*f^2) + (2*d^2*(c + d*x)*tan(e/2 + (f*x)/2))/(a^2*f^3) - (5*(c + d*x)^3*tan(e/2 + (f*x)/2))/(3*a^2*f) + ((c + d*x)^3*sec(e/2 + (f*x)/2)^2*tan(e/2 + (f*x)/2))/(6*a^2*f)]
    @test_int [1/(a + a*sec(e + f*x))^2*(c + d*x)^2, x, 17, (5*im*(c + d*x)^2)/(3*a^2*f) + (c + d*x)^3/(3*a^2*d) - (20*d*(c + d*x)*log(1 + ℯ^(im*(e + f*x))))/(3*a^2*f^2) + (20*im*d^2*PolyLog(2, -ℯ^(im*(e + f*x))))/(3*a^2*f^3) - (d*(c + d*x)*sec(e/2 + (f*x)/2)^2)/(3*a^2*f^2) + (2*d^2*tan(e/2 + (f*x)/2))/(3*a^2*f^3) - (5*(c + d*x)^2*tan(e/2 + (f*x)/2))/(3*a^2*f) + ((c + d*x)^2*sec(e/2 + (f*x)/2)^2*tan(e/2 + (f*x)/2))/(6*a^2*f)]
    @test_int [1/(a + a*sec(e + f*x))^2*(c + d*x)^1, x, 9, (c + d*x)^2/(2*a^2*d) - (10*d*log(cos(e/2 + (f*x)/2)))/(3*a^2*f^2) - (d*sec(e/2 + (f*x)/2)^2)/(6*a^2*f^2) - (5*(c + d*x)*tan(e/2 + (f*x)/2))/(3*a^2*f) + ((c + d*x)*sec(e/2 + (f*x)/2)^2*tan(e/2 + (f*x)/2))/(6*a^2*f)]
    @test_int [1/(a + a*sec(e + f*x))^2/(c + d*x)^1, x, 0, Unintegrable(1/((c + d*x)*(a + a*sec(e + f*x))^2), x)]
    @test_int [1/(a + a*sec(e + f*x))^2/(c + d*x)^2, x, 0, Unintegrable(1/((c + d*x)^2*(a + a*sec(e + f*x))^2), x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+a*sec(e+f*x))^n*with*m*symbolic=#

    @test_int [(c + d*x)^m*(a + a*sec(e + f*x))^n, x, 0, Unintegrable((c + d*x)^m*(a + a*sec(e + f*x))^n, x)]

    @test_int [(c + d*x)^m*(a + a*sec(e + f*x))^1, x, 0, Unintegrable((c + d*x)^m*(a + a*sec(e + f*x)), x)]
    @test_int [(c + d*x)^m/(a + a*sec(e + f*x))^1, x, 0, Unintegrable((c + d*x)^m/(a + a*sec(e + f*x)), x)]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*sec(e+f*x))^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*sec(e+f*x))^n=#

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [(a + b*sec(e + f*x))*(c + d*x)^3, x, 11, (a*(c + d*x)^4)/(4*d) - (2*im*b*(c + d*x)^3*atan(ℯ^(im*(e + f*x))))/f + (3*im*b*d*(c + d*x)^2*PolyLog(2, (-im)*ℯ^(im*(e + f*x))))/f^2 - (3*im*b*d*(c + d*x)^2*PolyLog(2, im*ℯ^(im*(e + f*x))))/f^2 - (6*b*d^2*(c + d*x)*PolyLog(3, (-im)*ℯ^(im*(e + f*x))))/f^3 + (6*b*d^2*(c + d*x)*PolyLog(3, im*ℯ^(im*(e + f*x))))/f^3 - (6*im*b*d^3*PolyLog(4, (-im)*ℯ^(im*(e + f*x))))/f^4 + (6*im*b*d^3*PolyLog(4, im*ℯ^(im*(e + f*x))))/f^4]
    @test_int [(a + b*sec(e + f*x))*(c + d*x)^2, x, 9, (a*(c + d*x)^3)/(3*d) - (2*im*b*(c + d*x)^2*atan(ℯ^(im*(e + f*x))))/f + (2*im*b*d*(c + d*x)*PolyLog(2, (-im)*ℯ^(im*(e + f*x))))/f^2 - (2*im*b*d*(c + d*x)*PolyLog(2, im*ℯ^(im*(e + f*x))))/f^2 - (2*b*d^2*PolyLog(3, (-im)*ℯ^(im*(e + f*x))))/f^3 + (2*b*d^2*PolyLog(3, im*ℯ^(im*(e + f*x))))/f^3]
    @test_int [(a + b*sec(e + f*x))*(c + d*x)^1, x, 7, (a*(c + d*x)^2)/(2*d) - (2*im*b*(c + d*x)*atan(ℯ^(im*(e + f*x))))/f + (im*b*d*PolyLog(2, (-im)*ℯ^(im*(e + f*x))))/f^2 - (im*b*d*PolyLog(2, im*ℯ^(im*(e + f*x))))/f^2]
    @test_int [(a + b*sec(e + f*x))/(c + d*x)^1, x, 0, Unintegrable((a + b*sec(e + f*x))/(c + d*x), x)]
    @test_int [(a + b*sec(e + f*x))/(c + d*x)^2, x, 0, Unintegrable((a + b*sec(e + f*x))/(c + d*x)^2, x)]

    @test_int [(a + b*sec(e + f*x))^2*(c + d*x)^3, x, 17, -((im*b^2*(c + d*x)^3)/f) + (a^2*(c + d*x)^4)/(4*d) - (4*im*a*b*(c + d*x)^3*atan(ℯ^(im*(e + f*x))))/f + (3*b^2*d*(c + d*x)^2*log(1 + ℯ^(2*im*(e + f*x))))/f^2 + (6*im*a*b*d*(c + d*x)^2*PolyLog(2, (-im)*ℯ^(im*(e + f*x))))/f^2 - (6*im*a*b*d*(c + d*x)^2*PolyLog(2, im*ℯ^(im*(e + f*x))))/f^2 - (3*im*b^2*d^2*(c + d*x)*PolyLog(2, -ℯ^(2*im*(e + f*x))))/f^3 - (12*a*b*d^2*(c + d*x)*PolyLog(3, (-im)*ℯ^(im*(e + f*x))))/f^3 + (12*a*b*d^2*(c + d*x)*PolyLog(3, im*ℯ^(im*(e + f*x))))/f^3 + (3*b^2*d^3*PolyLog(3, -ℯ^(2*im*(e + f*x))))/(2*f^4) - (12*im*a*b*d^3*PolyLog(4, (-im)*ℯ^(im*(e + f*x))))/f^4 + (12*im*a*b*d^3*PolyLog(4, im*ℯ^(im*(e + f*x))))/f^4 + (b^2*(c + d*x)^3*tan(e + f*x))/f]
    @test_int [(a + b*sec(e + f*x))^2*(c + d*x)^2, x, 14, -((im*b^2*(c + d*x)^2)/f) + (a^2*(c + d*x)^3)/(3*d) - (4*im*a*b*(c + d*x)^2*atan(ℯ^(im*(e + f*x))))/f + (2*b^2*d*(c + d*x)*log(1 + ℯ^(2*im*(e + f*x))))/f^2 + (4*im*a*b*d*(c + d*x)*PolyLog(2, (-im)*ℯ^(im*(e + f*x))))/f^2 - (4*im*a*b*d*(c + d*x)*PolyLog(2, im*ℯ^(im*(e + f*x))))/f^2 - (im*b^2*d^2*PolyLog(2, -ℯ^(2*im*(e + f*x))))/f^3 - (4*a*b*d^2*PolyLog(3, (-im)*ℯ^(im*(e + f*x))))/f^3 + (4*a*b*d^2*PolyLog(3, im*ℯ^(im*(e + f*x))))/f^3 + (b^2*(c + d*x)^2*tan(e + f*x))/f]
    @test_int [(a + b*sec(e + f*x))^2*(c + d*x)^1, x, 9, (a^2*(c + d*x)^2)/(2*d) - (4*im*a*b*(c + d*x)*atan(ℯ^(im*(e + f*x))))/f + (b^2*d*log(cos(e + f*x)))/f^2 + (2*im*a*b*d*PolyLog(2, (-im)*ℯ^(im*(e + f*x))))/f^2 - (2*im*a*b*d*PolyLog(2, im*ℯ^(im*(e + f*x))))/f^2 + (b^2*(c + d*x)*tan(e + f*x))/f]
    @test_int [(a + b*sec(e + f*x))^2/(c + d*x)^1, x, 0, Unintegrable((a + b*sec(e + f*x))^2/(c + d*x), x)]
    @test_int [(a + b*sec(e + f*x))^2/(c + d*x)^2, x, 0, Unintegrable((a + b*sec(e + f*x))^2/(c + d*x)^2, x)]

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [(c + d*x)^3/(a + b*sec(e + f*x)), x, 14, (c + d*x)^4/(4*a*d) + (im*b*(c + d*x)^3*log(1 + (a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*f) - (im*b*(c + d*x)^3*log(1 + (a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*f) + (3*b*d*(c + d*x)^2*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*f^2) - (3*b*d*(c + d*x)^2*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*f^2) + (6*im*b*d^2*(c + d*x)*PolyLog(3, -((a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*f^3) - (6*im*b*d^2*(c + d*x)*PolyLog(3, -((a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*f^3) - (6*b*d^3*PolyLog(4, -((a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*f^4) + (6*b*d^3*PolyLog(4, -((a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*f^4)]
    @test_int [(c + d*x)^2/(a + b*sec(e + f*x)), x, 12, (c + d*x)^3/(3*a*d) + (im*b*(c + d*x)^2*log(1 + (a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*f) - (im*b*(c + d*x)^2*log(1 + (a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*f) + (2*b*d*(c + d*x)*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*f^2) - (2*b*d*(c + d*x)*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*f^2) + (2*im*b*d^2*PolyLog(3, -((a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*f^3) - (2*im*b*d^2*PolyLog(3, -((a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*f^3)]
    @test_int [(c + d*x)^1/(a + b*sec(e + f*x)), x, 10, (c + d*x)^2/(2*a*d) + (im*b*(c + d*x)*log(1 + (a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*f) - (im*b*(c + d*x)*log(1 + (a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*f) + (b*d*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*f^2) - (b*d*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*f^2)]
    @test_int [1/((c + d*x)^1*(a + b*sec(e + f*x))), x, 0, Unintegrable(1/((c + d*x)*(a + b*sec(e + f*x))), x)]
    @test_int [1/((c + d*x)^2*(a + b*sec(e + f*x))), x, 0, Unintegrable(1/((c + d*x)^2*(a + b*sec(e + f*x))), x)]

    @test_int [(c + d*x)^3/(a + b*sec(e + f*x))^2, x, 36, -((im*b^2*(c + d*x)^3)/(a^2*(a^2 - b^2)*f)) + (c + d*x)^4/(4*a^2*d) + (3*b^2*d*(c + d*x)^2*log(1 + (a*ℯ^(im*(e + f*x)))/(b - im*sqrt(a^2 - b^2))))/(a^2*(a^2 - b^2)*f^2) + (3*b^2*d*(c + d*x)^2*log(1 + (a*ℯ^(im*(e + f*x)))/(b + im*sqrt(a^2 - b^2))))/(a^2*(a^2 - b^2)*f^2) - (im*b^3*(c + d*x)^3*log(1 + (a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*f) + (2*im*b*(c + d*x)^3*log(1 + (a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*f) + (im*b^3*(c + d*x)^3*log(1 + (a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*f) - (2*im*b*(c + d*x)^3*log(1 + (a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*f) - (6*im*b^2*d^2*(c + d*x)*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b - im*sqrt(a^2 - b^2)))))/(a^2*(a^2 - b^2)*f^3) - (6*im*b^2*d^2*(c + d*x)*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b + im*sqrt(a^2 - b^2)))))/(a^2*(a^2 - b^2)*f^3) - (3*b^3*d*(c + d*x)^2*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*f^2) + (6*b*d*(c + d*x)^2*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*f^2) + (3*b^3*d*(c + d*x)^2*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*f^2) - (6*b*d*(c + d*x)^2*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*f^2) + (6*b^2*d^3*PolyLog(3, -((a*ℯ^(im*(e + f*x)))/(b - im*sqrt(a^2 - b^2)))))/(a^2*(a^2 - b^2)*f^4) + (6*b^2*d^3*PolyLog(3, -((a*ℯ^(im*(e + f*x)))/(b + im*sqrt(a^2 - b^2)))))/(a^2*(a^2 - b^2)*f^4) - (6*im*b^3*d^2*(c + d*x)*PolyLog(3, -((a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*f^3) + (12*im*b*d^2*(c + d*x)*PolyLog(3, -((a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*f^3) + (6*im*b^3*d^2*(c + d*x)*PolyLog(3, -((a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*f^3) - (12*im*b*d^2*(c + d*x)*PolyLog(3, -((a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*f^3) + (6*b^3*d^3*PolyLog(4, -((a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*f^4) - (12*b*d^3*PolyLog(4, -((a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*f^4) - (6*b^3*d^3*PolyLog(4, -((a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*f^4) + (12*b*d^3*PolyLog(4, -((a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*f^4) + (b^2*(c + d*x)^3*sin(e + f*x))/(a*(a^2 - b^2)*f*(b + a*cos(e + f*x)))]
    @test_int [(c + d*x)^2/(a + b*sec(e + f*x))^2, x, 30, -((im*b^2*(c + d*x)^2)/(a^2*(a^2 - b^2)*f)) + (c + d*x)^3/(3*a^2*d) + (2*b^2*d*(c + d*x)*log(1 + (a*ℯ^(im*(e + f*x)))/(b - im*sqrt(a^2 - b^2))))/(a^2*(a^2 - b^2)*f^2) + (2*b^2*d*(c + d*x)*log(1 + (a*ℯ^(im*(e + f*x)))/(b + im*sqrt(a^2 - b^2))))/(a^2*(a^2 - b^2)*f^2) - (im*b^3*(c + d*x)^2*log(1 + (a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*f) + (2*im*b*(c + d*x)^2*log(1 + (a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*f) + (im*b^3*(c + d*x)^2*log(1 + (a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*f) - (2*im*b*(c + d*x)^2*log(1 + (a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*f) - (2*im*b^2*d^2*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b - im*sqrt(a^2 - b^2)))))/(a^2*(a^2 - b^2)*f^3) - (2*im*b^2*d^2*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b + im*sqrt(a^2 - b^2)))))/(a^2*(a^2 - b^2)*f^3) - (2*b^3*d*(c + d*x)*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*f^2) + (4*b*d*(c + d*x)*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*f^2) + (2*b^3*d*(c + d*x)*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*f^2) - (4*b*d*(c + d*x)*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*f^2) - (2*im*b^3*d^2*PolyLog(3, -((a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*f^3) + (4*im*b*d^2*PolyLog(3, -((a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*f^3) + (2*im*b^3*d^2*PolyLog(3, -((a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*f^3) - (4*im*b*d^2*PolyLog(3, -((a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*f^3) + (b^2*(c + d*x)^2*sin(e + f*x))/(a*(a^2 - b^2)*f*(b + a*cos(e + f*x)))]
    @test_int [(c + d*x)^1/(a + b*sec(e + f*x))^2, x, 21, (c + d*x)^2/(2*a^2*d) - (im*b^3*(c + d*x)*log(1 + (a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*f) + (2*im*b*(c + d*x)*log(1 + (a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*f) + (im*b^3*(c + d*x)*log(1 + (a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*f) - (2*im*b*(c + d*x)*log(1 + (a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*f) + (b^2*d*log(b + a*cos(e + f*x)))/(a^2*(a^2 - b^2)*f^2) - (b^3*d*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*f^2) + (2*b*d*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*f^2) + (b^3*d*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*f^2) - (2*b*d*PolyLog(2, -((a*ℯ^(im*(e + f*x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*f^2) + (b^2*(c + d*x)*sin(e + f*x))/(a*(a^2 - b^2)*f*(b + a*cos(e + f*x)))]
    @test_int [1/((c + d*x)^1*(a + b*sec(e + f*x))^2), x, 0, Unintegrable(1/((c + d*x)*(a + b*sec(e + f*x))^2), x)]
    @test_int [1/((c + d*x)^2*(a + b*sec(e + f*x))^2), x, 0, Unintegrable(1/((c + d*x)^2*(a + b*sec(e + f*x))^2), x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*sec(e+f*x))^n*with*m*symbolic=#

    @test_int [(c + d*x)^m*(a + b*sec(e + f*x))^n, x, 0, Unintegrable((c + d*x)^m*(a + b*sec(e + f*x))^n, x)]

    @test_int [(c + d*x)^m*(a + b*sec(e + f*x))^1, x, 0, Unintegrable((c + d*x)^m*(a + b*sec(e + f*x)), x)]
    @test_int [(c + d*x)^m/(a + b*sec(e + f*x))^1, x, 0, Unintegrable((c + d*x)^m/(a + b*sec(e + f*x)), x)]
    end

max_stars_repo_path: test/4 Trig functions/4.5 Secant/4.5.10 (c+d x)^m (a+b sec)^n.jl


text:
export PIBasis

# ---------------------- Implementation of the PIBasisSpec

"""
`struct PIBasisSpec`
"""
struct PIBasisSpec
   orders::Vector{Int}     # order (length) of ith basis function
   iAA2iA::Matrix{Int}     # where in A can we find the ith basis function
end

==(B1::PIBasisSpec, B2::PIBasisSpec) = _allfieldsequal(B1, B2)

Base.length(spec::PIBasisSpec) = length(spec.orders)

maxcorrorder(spec::PIBasisSpec) = size(spec.iAA2iA, 2)

function _get_pibfcn(spec0, Aspec, vv)
   vv1 = vv[2:end]
   vvnz = vv1[findall(vv1 .!= 0)]
   return (spec0[vv[1]], Aspec[vvnz])
end

function _get_pibfcn(Aspec, vv)
   vvnz = vv[findall(vv .!= 0)]
   return Aspec[vvnz]
end

# TODO: maybe instead of property == nothing, there should be a
#       generic property with no symmetry attached to it.

function PIBasisSpec( basis1p::OneParticleBasis,
                      symgrp::SymmetryGroup,
                      Bsel::DownsetBasisSelector;
                      property = nothing,
                      filterfun = _->true,
                      init1pbasis = true )

   # we initialize the 1p-basis here; to prevent this it must be manually
   # avoided by passing in init1pbasis = false
   if init1pbasis
      init1pspec!(basis1p, Bsel)
   end

   # get the basis spec of the one-particle basis
   #  Aspec[i] described the basis function that will get written into A[i]
   Aspec = get_spec(basis1p)

   # we assume that `Aspec` is sorted by degree, but best to double-check this
   # since the notion of degree used to construct `Aspec` might be different
   # from the one used to construct AAspec.
   if !issorted(Aspec; by = b -> level(b, Bsel, basis1p))
      error("""PIBasisSpec : AAspec construction failed because Aspec is not
               sorted by degree. This could e.g. happen if an incompatible
               notion of degree was used to construct the 1-p basis spec.""")
   end
   # An AA basis function is given by a tuple 𝒗 = vv. Each index 𝒗ᵢ = vv[i]
   # corresponds to the basis function Aspec[𝒗ᵢ] and the tuple
   # 𝒗 = (𝒗₁, ...) to a product basis function
   #   ∏ A_{vₐ}
   tup2b = vv -> _get_pibfcn(Aspec, vv)

   #  degree or level of a basis function ↦ is it admissible?
   admissible = bb -> (level(bb, Bsel, basis1p) <= maxlevel(Bsel, basis1p))

   if property != nothing
      filter1 = bb -> filterfun(bb) && filter(bb, Bsel, basis1p) && filter(property, symgrp, bb)
   else
      filter1 = bb -> filterfun(bb) && filter(bb, Bsel, basis1p)
   end

   # we can now construct the basis specification; the `ordered = true`
   # keyword signifies that this is a permutation-invariant basis
   maxord = maxorder(Bsel)
   AAspec = gensparse(; NU = maxorder(Bsel),
                        tup2b = tup2b,
                        admissible = admissible,
                        ordered = true,
                        maxvv = [length(Aspec) for _=1:maxord],
                        filter = filter1)

   return PIBasisSpec(AAspec)
end

function PIBasisSpec(AAspec)
   orders = zeros(Int, length(AAspec))
   iAA2iA = zeros(Int, (length(AAspec), length(AAspec[1])))
   for (iAA, vv) in enumerate(AAspec)
      # we use reverse because gensparse constructs the indices in
      # ascending order, but we want descending here.
      # (I don't remember why though)
      iAA2iA[iAA, :] .= reverse(vv)
      orders[iAA] = length( findall( vv .!= 0 ) )
   end
   return PIBasisSpec(orders, iAA2iA)
end

get_spec(AAspec::PIBasisSpec, i::Integer) = AAspec.iAA2iA[i, 1:AAspec.orders[i]]

# --------------------------------- PIBasis implementation

"""
`mutable struct PIBasis:` implementation of a permutation-invariant
basis based on the density projection trick.

The standard constructor is
```
PIBasis(basis1p, N, D, maxdeg)
```
* `basis1p` : a one-particle basis
* `N` : maximum interaction order
* `D` : an abstract degee specification, e.g., SparsePSHDegree
* `maxdeg` : the maximum polynomial degree as measured by `D`
"""
mutable struct PIBasis{BOP, REAL, TB, TA} <: ACEBasis
   basis1p::BOP             # a one-particle basis
   spec::PIBasisSpec
   real::REAL     # could be `real` or `identity` to keep AA complex
   # evaluator    # classic vs graph
   B_pool::VectorPool{TB}
   dAA_pool::VectorPool{TA}
end

cutoff(basis::PIBasis) = cutoff(basis.basis1p)

==(B1::PIBasis, B2::PIBasis) =
      ( (B1.basis1p == B2.basis1p) &&
        (B1.spec == B2.spec) &&
        (B1.real == B2.real) )

valtype(basis::PIBasis) = basis.real( valtype(basis.basis1p) )

valtype(basis::PIBasis, cfg::AbstractConfiguration) =
      basis.real( valtype(basis.basis1p, cfg) )

gradtype(basis::PIBasis, cfgorX) =
      basis.real( gradtype(basis.basis1p, cfgorX) )

Base.length(basis::PIBasis) = length(basis.spec)

# default symmetry group
PIBasis(basis1p, Bsel::AbstractBasisSelector; kwargs...) =
   PIBasis(basis1p, O3(), Bsel; kwargs...)

PIBasis(basis1p, symgrp, Bsel::AbstractBasisSelector;
        isreal = false, kwargs...) =
   PIBasis(basis1p,
           PIBasisSpec(basis1p, symgrp, Bsel; kwargs...),
           isreal ? Base.real : Base.identity )

function PIBasis(basis1p::OneParticleBasis, spec::PIBasisSpec, real)
   VT1 = valtype(basis1p)
   VT = real(VT1)  # default valtype
   B_pool = VectorPool{VT}()
   dAA_pool = VectorPool{VT1}()
   return PIBasis(basis1p, spec, real, B_pool, dAA_pool)
end

get_spec(pibasis::PIBasis) =
   [ get_spec(pibasis, i) for i = 1:length(pibasis) ]

get_spec(pibasis::PIBasis, i::Integer) =
      get_spec.( Ref(pibasis.basis1p), get_spec(pibasis.spec, i) )

setreal(basis::PIBasis, isreal::Bool) =
   PIBasis(basis.basis1p, basis.spec, isreal)

maxcorrorder(basis::PIBasis) = maxcorrorder(basis.spec)

# TODO: this is a hack; cf. #68
function scaling(pibasis::PIBasis, p)
   _absvaluep(x::Number) = abs(x)^p
   _absvaluep(x::Symbol) = 0
   ww = zeros(Float64, length(pibasis))
   bspec = get_spec(pibasis)
   for i = 1:length(pibasis)
      for b in bspec[i]
         # TODO: revisit how this should be implemented for a general basis
         ww[i] += sum(_absvaluep, b)  #  abs.(values(b)).^p
      end
   end
   return ww
end

# function scaling(pibasis::PIBasis, p)
#    ww = zeros(Float64, length(pibasis))
#    for iz0 = 1:numz(pibasis)
#       wwin = @view ww[pibasis.inner[iz0].AAindices]
#       for i = 1:length(pibasis.inner[iz0])
#          bspec = get_basis_spec(pibasis, iz0, i)
#          wwin[i] = scaling(bspec, p)
#       end
#    end
#    return ww
# end

# graphevaluator(basis::PIBasis) =
#    PIBasis(basis.basis1p, zlist(basis), basis.inner, DAGEvaluator())
#
# standardevaluator(basis::PIBasis) =
#    PIBasis(basis.basis1p, zlist(basis), basis.inner, StandardEvaluator())

# -------------------------------------------------
# FIO codes

write_dict(basis::PIBasis) =
   Dict(  "__id__" => "ACE_PIBasis",
         "basis1p" => write_dict(basis.basis1p),
            "spec" => write_dict(basis.spec),
            "real" => basis.real == Base.real ? true : false )

read_dict(::Val{:ACE_PIBasis}, D::Dict) =
   PIBasis( read_dict(D["basis1p"]),
            read_dict(D["spec"]),
            D["real"] ? Base.real : Base.identity )

write_dict(spec::PIBasisSpec) =
   Dict( "__id__" => "ACE_PIBasisSpec",
         "orders" => spec.orders,
         "iAA2iA" => write_dict(spec.iAA2iA) )

read_dict(::Val{:ACE_PIBasisSpec}, D::Dict) =
   PIBasisSpec( D["orders"], read_dict(D["iAA2iA"]) )

# -------------------------------------------------
# Evaluation codes

function evaluate!(AA, basis::PIBasis, config::AbstractConfiguration)
   A = acquire_B!(basis.basis1p, config)   #  THIS ALLOCATES!!!!
   evaluate!(A, basis.basis1p, config)
   fill!(AA, 1)
   for iAA = 1:length(basis)
      aa = one(eltype(A))
      for t = 1:basis.spec.orders[iAA]
         aa *= A[ basis.spec.iAA2iA[ iAA, t ] ]
      end
      AA[iAA] = basis.real(aa)
   end
   release_B!(basis.basis1p, A)
   return AA
end

# -------------------------------------------------
# gradients

function evaluate_ed!(AA, dAA, basis::PIBasis,
                      cfg::AbstractConfiguration, args...)
   A = acquire_B!(basis.basis1p, cfg)
   dA = acquire_dB!(basis.basis1p, cfg)   # TODO: THIS WILL ALLOCATE!!!!!
   evaluate_ed!(A, dA, basis.basis1p, cfg, args...)
   evaluate_ed!(AA, dAA, basis, A, dA)
   release_dB!(basis.basis1p, dA)
   release_B!(basis.basis1p, A)
   return AA, dAA
end

function _AA_local_adjoints!(dAAdA, A, iAA2iA, iAA, ord, _real)
   @assert length(dAAdA) >= ord
   # TODO - optimize a bit more? can move one operation out of the loop
   # Forward pass:
   @inbounds dAAdA[1] = 1
   @inbounds AAfwd = A[iAA2iA[iAA, 1]]
   @inbounds for a = 2:ord
      dAAdA[a] = AAfwd
      AAfwd *= A[iAA2iA[iAA, a]]
   end
   aa = _real(AAfwd)
   # backward pass
   @inbounds AAbwd = A[iAA2iA[iAA, ord]]
   @inbounds for a = ord-1:-1:1
      dAAdA[a] *= AAbwd
      AAbwd *= A[iAA2iA[iAA, a]]
   end

   return aa
end

_acquire_dAAdA!(basis::PIBasis) = acquire!(basis.dAA_pool, maxcorrorder(basis))

function evaluate_ed!(AA, dAA, basis::PIBasis,
                      A::AbstractVector, dA::AbstractMatrix)
   orders = basis.spec.orders
   iAA2iA = basis.spec.iAA2iA
   dAAdA = _acquire_dAAdA!(basis)

   # Must treat the constants separately. This is not so elegant and could
   # maybe be improved?
   if orders[1] == 0  # SHOULD BE THE ONLY ONE with ord=0!!
      iAAinit = 2
      AA[1] = 1.0
      dAA[1, :] .= Ref(zero(eltype(dAA)))
   else
      iAAinit = 1
   end

   for iAA = iAAinit:length(basis)
      ord = orders[iAA]

      # ----- compute the local adjoints dAA / dA
      # dAAdA[a] ← ∏_{t ≂̸ a} A_{v_t}
      AA[iAA] = _AA_local_adjoints!(dAAdA, A, iAA2iA, iAA, orders[iAA], basis.real)

      # ----- now convert them into dAA / dX
      for j = 1:size(dA, 2)
         dAA[iAA, j] = sum(dAAdA[a] * dA[iAA2iA[iAA, a], j]
                           for a = 1:ord) |> basis.real
      end
   end

   return AA, dAA
end


text:
using Clang
using Test

include("generators.jl")


text:
#!/usr/bin/env julia
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 1, Zadanie 1.5.
#

println("floatmin(Float32) = ", floatmin(Float32), ", floatmin(Float64) = ", floatmin(Float64))

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
#RSDeltaSigma: simulateDSM() algorithm
#-------------------------------------------------------------------------------

#==simulateDSM()
===============================================================================#
#Main algorithm for simulateDSM
function _simulateDSM(u, nq, nlev, x0, order, A, B, C, D1, savestate::Bool, trackmax::Bool)
	if isnan(x0)
		x0 = zeros(order,1)
	end

	N = length(u)
	v = zeros(nq,N)
	y = zeros(nq,N)

	xn = nothing; xmax = nothing
	if savestate #Need to store the state information
		xn = zeros(order,N)
	end
	if trackmax #Need to keep track of the state maxima
		xmax = abs.(x0)
	end

	for i=1:N
		y[:,i] = C*x0 .+ D1 .* u[:,i]
		v[:,i] = ds_quantize_DSM(y[:,i], nlev)
		x0 = A * x0 .+ B * [u[:,i]; v[:,i]]
		if savestate #Save the next state
			xn[:,i] = x0
		end
		if trackmax #Keep track of the state maxima
			xmax = max.(abs.(x0), xmax)
		end
	end

	result = (v=v, y=y)
	if savestate
		result = merge(result, (xn=xn,))
	end
	if trackmax
		result = merge(result, (xmax=xmax,))
	end

	return result
end

function simulateDSM(u, ntf::ZPKData; nlev=2, x0=NaN, savestate::Bool=false, trackmax::Bool=false)
	u = conv2seriesmatrix2D(u)
	nq = length(nlev) #Number of quantizers
	order = length(ntf.z)

	_ss = _zp2ss(ntf.p,ntf.z,-1) #realization of 1/H
	A, B2, C, D2 = _ss.A, _ss.B, _ss.C, _ss.D
	#Transform the realization so that C = [1 0 0 ...]
	Sinv = orth([C' eye(order)])/norm(C); S = inv(Sinv)
	C = C*Sinv
	if C[1]<0
		S = -S
		Sinv = -Sinv
	end
	A = S*A*Sinv; B2 = S*B2; C = [1 zeros(1,order-1)]; #C=C*Sinv
	D2 = 0
	#!!!! Assume stf=1
	B1 = -B2
	D1 = 1
	B = [B1 B2]

	return _simulateDSM(u, nq, nlev, x0, order, A, B, C, D1, savestate, trackmax)
end

function simulateDSM(u, ABCD::Array; nlev=2, x0=NaN, savestate::Bool=false, trackmax::Bool=false)
	u = conv2seriesmatrix2D(u)
	nu, nq, order = get_nu_nq_order(u, ABCD, nlev)
	A = ABCD[1:order, 1:order]
	B = ABCD[1:order, order+1:order+nu+nq]
	C = ABCD[order+1:order+nq, 1:order]
	D1= ABCD[order+1:order+nq, order+1:order+nu]

	return _simulateDSM(u, nq, nlev, x0, order, A, B, C, D1, savestate, trackmax)
end

"""`simresult = simulateDSM(u, [NTF|ABCD]; nlev=2, x0=NaN, savestate=false, trackmax=false))`

Compute the output of a general ΔΣ modulator.

# Inputs
 - `u[]`: Input signal (Multiple inputs are implied by the number of rows in `u`)
 - `[NTF|ABCD]`: NTF or ABCD (state-space)
 - `nlev`: Quantizer levels (Multiple quantizers are implied by making `nlev` an array)
 - `x0[]`: Initial state x0 (default zero).

The NTF is zpk object. (The STF is assumed to be 1.)
The structure that is simulated is the block-diagional structure used by
zp2ss.m.

# Returns `simresult::NamedTuple`:
 - .v: Output signal of ΔΣ modulator.
 - .xn: Internal state of ΔΣ modulator.
 - .xmax: Maximum internal state of ΔΣ modulator.
 - .y:
""" simulateDSM

@warn("Implements .m version of simulateDSM(). C code might be more efficient/validated.")

#==simulateQDSM
===============================================================================#
simulateQDSM(args...; kwargs...) = throw("simulateQDSM() not implemented")

#Last line


text:
@testset "Testing logic using strings" begin
    if VERSION < v"0.6.0-dev.1833" # added in Julia PR #19548
        Base.one(::Type{T}) where {T<:AbstractString} = ""
    end
    z = ["x ", "y ", "z "]

    want_1 = ["" "x " "y " "z "]
    @test complete_polynomial(reshape(z, 1, length(z)), 1) == want_1
    # Test vector methods to make sure generates same thing as matrix methods
    @test complete_polynomial(z, 1) == vec(want_1)

    want_2 = ["" "x " "x x " "x y " "x z " "y " "y y " "y z " "z " "z z "]
    @test complete_polynomial(reshape(z, 1, length(z)), 2) == want_2
    # Test vector methods to make sure generates same thing as matrix methods
    @test complete_polynomial(z, 2) == vec(want_2)

end

@testset "testing numerical routintes" begin
    z = rand(100, 10)
    the_ones = ones(size(z, 1))
    the_zeros = zeros(size(z, 1))
    for d in 1:5
        n_comp = n_complete(size(z, 2), d)
        # Test basis matrices
        buff_d = Array{Float64}(size(z, 1), n_comp)
        out_d = complete_polynomial(z, d)
        complete_polynomial!(buff_d, z, d)
        @test size(out_d, 1) == size(z, 1)
        @test size(out_d, 2) == n_comp
        @test all(isapprox.(out_d[:, 1], the_ones))
        @test all(isapprox.(out_d[:, 2], z[:, 1]))
        @test d == 1 ? true : all(isapprox.(out_d[:, 3], z[:, 1].*z[:, 1]))
        @test all(isapprox.(out_d[:, end], z[:, end].^d))
        @test all(isapprox.(out_d, buff_d))

        # Test derivatives
        buff_der_d = Array{Float64}(size(z, 1), n_comp)
        out_der_d = complete_polynomial(z, d, 1)
        complete_polynomial!(buff_der_d, z, d, 1)
        @test size(out_der_d, 1) == size(z, 1)
        @test size(out_der_d, 2) == n_comp
        @test all(isapprox.(out_der_d[:, 1], the_zeros))
        @test all(isapprox.(out_der_d[:, 2], the_ones))
        @test d ==1 ? true : all(isapprox.(out_der_d[:, 3], 2.*z[:, 1]))
        @test d <= 2 ? true : all(isapprox.(out_der_d[:, 4], 3.*z[:, 1].*z[:, 1]))
        @test all(isapprox.(out_der_d[:, end], the_zeros))
        @test all(isapprox.(out_der_d, buff_der_d))

    end

    z = rand(10, 3)
    z2 = z[1, :]
    all_zero = zeros(10)
    all_ones = ones(10)
    out_21 = complete_polynomial(z, 2)
    out_22 = complete_polynomial(z, BasisMatrices.Degree{2}())
    for out_2 in (out_21, out_22)
        # need to test columns size(z, 2) + 2:end
        @test all(isapprox.(out_2[:, 1], all_ones))
        @test all(isapprox.(out_2[:, 2], z[:, 1]))
        @test all(isapprox.(out_2[:, 3], z[:, 1].*z[:, 1]))
        @test all(isapprox.(out_2[:, 4], z[:, 1].*z[:, 2]))
        @test all(isapprox.(out_2[:, 5], z[:, 1].*z[:, 3]))
        @test all(isapprox.(out_2[:, 6], z[:, 2]))
        @test all(isapprox.(out_2[:, 7], z[:, 2].*z[:, 2]))
        @test all(isapprox.(out_2[:, 8], z[:, 2].*z[:, 3]))
        @test all(isapprox.(out_2[:, 9], z[:, 3]))
        @test all(isapprox.(out_2[:, 10], z[:, 3].*z[:, 3]))
        # Test vector methods to make sure generates same thing as matrix methods
        @test all(isapprox.(out_2[1, :], complete_polynomial(z2, 2)))
    end

    out_der_21 = complete_polynomial(z, 2, 2)
    out_der_22 = complete_polynomial(
        z, BasisMatrices.Degree{2}(), BasisMatrices.Derivative{2}()
    )
    for out_der_2 in (out_der_21, out_der_22)
        # need to test columns size(z, 2) + 2:end
        @test all(isapprox.(out_der_2[:, 1], all_zero))
        @test all(isapprox.(out_der_2[:, 2], all_zero))
        @test all(isapprox.(out_der_2[:, 3], all_zero))
        @test all(isapprox.(out_der_2[:, 4], z[:, 1]))
        @test all(isapprox.(out_der_2[:, 5], all_zero))
        @test all(isapprox.(out_der_2[:, 6], all_ones))
        @test all(isapprox.(out_der_2[:, 7], 2 .* z[:, 2]))
        @test all(isapprox.(out_der_2[:, 8], z[:, 3]))
        @test all(isapprox.(out_der_2[:, 9], all_zero))
        @test all(isapprox.(out_der_2[:, 10], all_zero))
    end

end


text:
module TestFinSets
using Test

using Catlab.Theories, Catlab.CategoricalAlgebra

sshow(args...) = sprint(show, args...)

# Finite sets
#############

# Finite ordinals.
set = FinSet(3)
@test length(set) == 3
@test collect(set) == 1:3
@test 2 ∈ set && 10 ∉ set
@test sshow(set) == "FinSet(3)"

# Collections as sets.
set = FinSet(Set(1:2:5))
@test length(set) == 3
@test sort!(collect(set)) == [1,3,5]
@test 3 ∈ set && 4 ∉ set
@test startswith(sshow(set), "FinSet(Set(")

# Tables as sets.
set = FinSet((x=[1,3,5], y=['a','b','c']))
@test length(set) == 3
@test map(NamedTuple, set) == [(x=1, y='a'), (x=3, y='b'), (x=5, y='c')]

# Functions between finite sets
###############################

f = FinFunction([1,3,4], 5)
g = FinFunction([1,1,2,2,3], 3)
h = FinFunction([3,1,2], 3)
@test (dom(f), codom(f)) == (FinSet(3), FinSet(5))
@test force(f) === f
@test codom(FinFunction([1,3,4])) == FinSet(4)

# Evaluation.
rot3(x) = (x % 3) + 1
@test map(f, 1:3) == [1,3,4]
@test map(FinFunction(rot3, 3, 3), 1:3) == [2,3,1]
@test map(id(FinSet(3)), 1:3) == [1,2,3]

# Composition.
@test compose(f,g) == FinFunction([1,2,2], 3)
@test compose(g,h) == FinFunction([3,3,1,1,2], 3)
@test compose(compose(f,g),h) == compose(f,compose(g,h))
@test compose(id(dom(f)), f) == f
@test compose(f, id(codom(f))) == f

# Indexing.
@test !is_indexed(f)
@test is_indexed(id(FinSet(3)))
@test preimage(id(FinSet(3)), 2) == [2]

f = FinFunction([1,2,1,3], 5, index=true)
@test is_indexed(f)
@test force(f) === f
@test (dom(f), codom(f)) == (FinSet(4), FinSet(5))
@test f(1) == 1
@test preimage(f, 1) == [1,3]
@test preimage(f, 3) == [4]
@test isempty(preimage(f, 4))

g = FinFunction(5:-1:1)
@test compose(f,g) == FinFunction([5,4,5,3])

# Pretty-print.
@test sshow(FinFunction(rot3, 3, 3)) ==
  "FinFunction(rot3, FinSet(3), FinSet(3))"
@test sshow(id(FinSet(3))) == "FinFunction(identity, FinSet(3))"
@test sshow(FinFunction([1,3,4], 5)) == "FinFunction($([1,3,4]), 3, 5)"
@test sshow(FinFunction([1,3,4], 5, index=true)) ==
  "FinFunction($([1,3,4]), 3, 5, index=true)"

# Functions out of finite sets
##############################

k = FinDomFunction([:a,:b,:c,:d,:e])
@test (dom(k), codom(k)) == (FinSet(5), TypeSet(Symbol))
@test k(3) == :c
@test collect(k) == [:a,:b,:c,:d,:e]
@test sshow(k) ==
  "FinDomFunction($([:a,:b,:c,:d,:e]), FinSet(5), TypeSet(Symbol))"

f = FinFunction([1,3,4], 5)
@test compose(f,k) == FinDomFunction([:a,:c,:d])

# Indexing.
@test !is_indexed(k)
@test preimage(k, :c) == [3]

k = FinDomFunction(5:10)
@test is_indexed(k)
@test preimage(k, 6) == [2]
@test isempty(preimage(k, 4))

k = FinDomFunction([:a,:b,:a,:c], index=true)
@test is_indexed(k)
@test (dom(k), codom(k)) == (FinSet(4), TypeSet(Symbol))
@test k(1) == :a
@test preimage(k, :a) == [1,3]
@test preimage(k, :c) == [4]
@test isempty(preimage(k, :d))
@test sshow(k) ==
  "FinDomFunction($([:a,:b,:a,:c]), FinSet(4), TypeSet(Symbol), index=true)"

f = FinFunction([1,3,2], 4)
@test compose(f,k) == FinDomFunction([:a,:a,:b])

# Limits
########

# Products
#---------

# Terminal object.
I = terminal(FinSet{Int})
@test ob(I) == FinSet(1)
@test force(delete(I, FinSet(3))) == FinFunction([1,1,1])

# Binary product.
lim = product(FinSet(2), FinSet(3))
@test ob(lim) == FinSet(6)
@test force(proj1(lim)) == FinFunction([1,2,1,2,1,2])
@test force(proj2(lim)) == FinFunction([1,1,2,2,3,3])

lim = product(FinSet(4), FinSet(3))
f, g = FinFunction([2,1,4]), FinFunction([1,3,2])
@test force(pair(lim,f,g) ⋅ proj1(lim)) == f
@test force(pair(lim,f,g) ⋅ proj2(lim)) == g

# N-ary product.
lim = product([FinSet(2), FinSet(3)])
@test ob(lim) == FinSet(6)
@test force.(legs(lim)) ==
  [FinFunction([1,2,1,2,1,2]), FinFunction([1,1,2,2,3,3])]
@test ob(product(FinSet{Int}[])) == FinSet(1)

lim = product([FinSet(4), FinSet(3)])
@test force(pair(lim,[f,g]) ⋅ first(legs(lim))) == f
@test force(pair(lim,[f,g]) ⋅ last(legs(lim))) == g

# Equalizers
#-----------

f, g = FinFunction([1,2,4,3]), FinFunction([3,2,4,1])
eq = equalizer(f,g)
@test incl(eq) == FinFunction([2,3], 4)
@test incl(equalizer([f,g])) == incl(eq)
@test factorize(eq, FinFunction([2,3,2])) == FinFunction([1,2,1])

# Equalizer of identical functions.
f = FinFunction([4,2,3,1], 5)
eq = equalizer(f,f)
@test incl(eq) == force(id(FinSet(4)))
@test incl(equalizer([f,f])) == incl(eq)
@test factorize(eq, FinFunction([2,1,3,3])) == FinFunction([2,1,3,3], 4)

# Equalizer matching nothing.
f, g = id(FinSet(5)), FinFunction([2,3,4,5,1], 5)
eq = equalizer(f,g)
@test incl(eq) == FinFunction(Int[], 5)
@test incl(equalizer([f,g])) == incl(eq)
@test factorize(eq, FinFunction(Int[])) == FinFunction(Int[])

# Equalizer of functions into non-finite set.
f = FinDomFunction([:a, :b, :d, :c])
g = FinDomFunction([:c, :b, :d, :a])
eq = equalizer(f,g)
@test incl(eq) == FinFunction([2,3], 4)

# Pullbacks
#----------

lim = pullback(FinFunction([1,1,3,2],4), FinFunction([1,1,4,2],4))
@test ob(lim) == FinSet(5)
@test force(proj1(lim)) == FinFunction([1,2,1,2,4], 4)
@test force(proj2(lim)) == FinFunction([1,1,2,2,4], 4)

f, g = FinFunction([1,2,4]), FinFunction([2,1,4])
@test force(pair(lim,f,g) ⋅ proj1(lim)) == f
@test force(pair(lim,f,g) ⋅ proj2(lim)) == g

# Pullback from a singleton set: the degenerate case of a product.
lim = pullback(FinFunction([1,1]), FinFunction([1,1,1]))
@test ob(lim) == FinSet(6)
@test force(proj1(lim)) == FinFunction([1,2,1,2,1,2])
@test force(proj2(lim)) == FinFunction([1,1,2,2,3,3])

f, g = FinFunction([1,1,2]), FinFunction([3,2,1])
@test force(pair(lim,f,g) ⋅ proj1(lim)) == f
@test force(pair(lim,f,g) ⋅ proj2(lim)) == g

# Pullback of a cospan into non-finite set.
f = FinDomFunction([:a, :a, :c, :b])
g = FinDomFunction([:a, :a, :d, :b])
π1, π2 = lim = pullback(f, g)
@test ob(lim) == FinSet(5)
@test force(π1) == FinFunction([1,2,1,2,4], 4)
@test force(π2) == FinFunction([1,1,2,2,4], 4)

# Pullback using different algorithms.
tuples(lim::AbstractLimit) =
  sort!([ Tuple(map(π -> π(i), legs(lim))) for i in ob(lim) ])

f, g = FinFunction([3,1,1,5,2],5), FinFunction([4,1,1,3,2],5)
for Alg in (NestedLoopJoin, SortMergeJoin, HashJoin)
  lim = pullback(f, g, alg=Alg())
  @test ob(lim) == FinSet(6)
  @test tuples(lim) == [(1,4), (2,2), (2,3), (3,2), (3,3), (5,5)]
end

# Ternary pullback using different algorithms.
f, g = FinDomFunction([:a,:b,:c]), FinDomFunction([:c,:b,:a])
h = FinDomFunction([:a,:a,:b,:b,:c,:c])
lim = limit(SMulticospan{3}(f, g, h), alg=ComposeProductEqualizer())
@test ob(lim) == FinSet(6)
reference_tuples = tuples(lim)

for Alg in (NestedLoopJoin, SortMergeJoin, HashJoin)
  lim = limit(SMulticospan{3}(f, g, h), alg=Alg())
  @test ob(lim) == FinSet(6)
  @test tuples(lim) == reference_tuples
end

# Pullback involving a constant, which should be handled specially.
lim = pullback(FinFunction([3], 4), FinFunction([1,3,4,2,3,3]), alg=SmartJoin())
@test ob(lim)== FinSet(3)
@test proj1(lim) == ConstantFunction(1, FinSet(3), FinSet(1))
@test proj2(lim) == FinFunction([2,5,6], 6)

# General limits
#---------------

# Pullback as limit of free diagram.
f, g = FinFunction([1,1,3,2],4), FinFunction([1,1,4,2],4)
lim = limit(FreeDiagram(Cospan(f, g)))
@test ob(lim) == FinSet(5)
π1, π2 = legs(lim)[1:2]
@test force(π1) == FinFunction([1,2,1,2,4], 4)
@test force(π2) == FinFunction([1,1,2,2,4], 4)

f′, g′ = FinFunction([1,2,4]), FinFunction([2,1,4])
h = universal(lim, Multispan([f′, g′, f′⋅f])) # f′⋅f == g′⋅g
@test force(h ⋅ π1) == f′
@test force(h ⋅ π2) == g′

# Pullback as limit of bipartite free diagram.
lim = limit(BipartiteFreeDiagram(Cospan(f, g)))
π1, π2 = legs(lim)
@test π1 == FinFunction([1,1,2,2,4], 4)
@test π2 == FinFunction([1,2,1,2,4], 4)
@test π1 ⋅ f == π2 ⋅ g

# Equalizer as limit of bipartite free diagram.
f, g = [FinDomFunction(x -> x % i, FinSet(100), TypeSet(Int)) for i in 2:3]
lim = (ι,) = limit(BipartiteFreeDiagram(ParallelPair(f, g)))
@test ι == incl(equalizer(f, g))

# Two pullbacks, which should be reduced to a single pullback by pairing.
f1, g1 = FinDomFunction([1,1,2,2,3,3]), FinDomFunction([:a,:a,:a,:b,:b,:b])
f2, g2 = FinDomFunction([1,2,3]), FinDomFunction([:a,:b,:c])
d = BipartiteFreeDiagram{SetOb,FinDomFunction{Int}}(2, 2)
d[:ob₁], d[:ob₂] = [FinSet(6), FinSet(3)], [TypeSet(Int), TypeSet(Symbol)]
add_edges!(d, [1,1,2,2], [1,2,1,2], hom=[f1,g1,f2,g2])
lim = π1, π2 = limit(d)
@test π1 == FinFunction([1,2,4], 6)
@test π2 == FinFunction([1,1,2], 3)

# Colimits
##########

# Coproducts
#-----------

# Initial object.
@test ob(initial(FinSet{Int})) == FinSet(0)
@test create(initial(FinSet{Int}), FinSet(3)) == FinFunction(Int[], 3)

# Binary coproduct.
colim = coproduct(FinSet(2), FinSet(3))
@test ob(colim) == FinSet(5)
@test coproj1(colim) == FinFunction([1,2], 5)
@test coproj2(colim) == FinFunction([3,4,5], 5)

f, g = FinFunction([3,5], 5), FinFunction([1,2,3], 5)
@test force(coproj1(colim) ⋅ copair(colim,f,g)) == f
@test force(coproj2(colim) ⋅ copair(colim,f,g)) == g

# N-ary coproduct.
colim = coproduct([FinSet(2), FinSet(3)])
@test ob(colim) == FinSet(5)
@test legs(colim) == [FinFunction([1,2], 5), FinFunction([3,4,5], 5)]
@test ob(coproduct(FinSet{Int}[])) == FinSet(0)

@test force(first(legs(colim)) ⋅ copair(colim,[f,g])) == f
@test force(last(legs(colim)) ⋅ copair(colim,[f,g])) == g

# Cocartesian monoidal structure.
@test FinSet(2)⊕FinSet(3) == FinSet(5)
@test oplus([FinSet(2), FinSet(3), FinSet(4)]) == FinSet(9)
@test f⊕g == FinFunction([3,5,6,7,8], 10)
@test mzero(FinSet{Int}) == FinSet(0)
@test swap(FinSet(2), FinSet(3)) == FinFunction([4,5,1,2,3])
ι1, ι2 = coproj1(FinSet(2),FinSet(3)), coproj2(FinSet(2),FinSet(3))
@test ι1 == FinFunction([1,2], 5)
@test ι2 == FinFunction([3,4,5], 5)

# Coequalizers
#-------------

# Coequalizer from a singleton set.
f, g = FinFunction([1], 3), FinFunction([3], 3)
coeq = coequalizer(f,g)
@test proj(coeq) == FinFunction([1,2,1], 2)
@test proj(coequalizer([f,g])) == proj(coeq)
@test factorize(coeq, FinFunction([4,1,4])) == FinFunction([4,1])
@test_throws AssertionError factorize(coeq, FinFunction([3,1,4]))

# Coequalizer in case of identical functions.
f = FinFunction([4,2,3,1], 5)
coeq = coequalizer(f,f)
@test proj(coeq) == force(id(FinSet(5)))
@test proj(coequalizer([f,f])) == proj(coeq)
@test factorize(coeq, FinFunction([2,1,3,3,4])) == FinFunction([2,1,3,3,4])

# Coequalizer identifying everything.
f, g = id(FinSet(5)), FinFunction([2,3,4,5,1], 5)
coeq = coequalizer(f,g)
@test proj(coeq) == FinFunction(fill(1,5))
@test proj(coequalizer([f,g])) == proj(coeq)
@test factorize(coeq, FinFunction(fill(3,5))) == FinFunction([3])

# Pushouts
#---------

# Pushout from the empty set: the degenerate case of the coproduct.
f, g = FinFunction(Int[], 2), FinFunction(Int[], 3)
colim = pushout(f,g)
@test ob(colim) == FinSet(5)
@test coproj1(colim) == FinFunction([1,2], 5)
@test coproj2(colim) == FinFunction([3,4,5], 5)

h, k = FinFunction([3,5], 5), FinFunction([1,2,3], 5)
@test force(coproj1(colim) ⋅ copair(colim,h,k)) == h
@test force(coproj2(colim) ⋅ copair(colim,h,k)) == k

# Pushout from a singleton set.
f, g = FinFunction([1], 2), FinFunction([2], 3)
colim = pushout(f,g)
@test ob(colim) == FinSet(4)
ι1, ι2 = colim
@test compose(f,ι1) == compose(g,ι2)
@test ι1 == FinFunction([1,2], 4)
@test ι2 == FinFunction([3,1,4], 4)

h, k = FinFunction([3,5]), FinFunction([1,3,5])
@test force(coproj1(colim) ⋅ copair(colim,h,k)) == h
@test force(coproj2(colim) ⋅ copair(colim,h,k)) == k
k = FinFunction([1,2,5])
@test_throws AssertionError copair(colim,h,k)

# Same thing as a colimit of a general free diagram.
diagram = FreeDiagram([FinSet(1),FinSet(2),FinSet(3)],[(f,1,2), (g,1,3)])
colim = colimit(diagram)
@test ob(colim) == FinSet(4)
_, ι1, ι2 = colim
@test compose(f,ι1) == compose(g,ι2)
@test ι1 == FinFunction([1,2], 4)
@test ι2 == FinFunction([3,1,4], 4)

h, k = FinFunction([3,5]), FinFunction([1,3,5])
ℓ = universal(colim, Multicospan([f⋅h, h, k])) # f⋅h == g⋅k
@test force(ι1 ⋅ ℓ) == h
@test force(ι2 ⋅ ℓ) == k

# Pushout from a two-element set, with non-injective legs.
f, g = FinFunction([1,1], 2), FinFunction([1,2], 3)
colim = pushout(f,g)
@test ob(colim) == FinSet(3)
ι1, ι2 = colim
@test compose(f,ι1) == compose(g,ι2)
@test ι1 == FinFunction([1,2], 3)
@test ι2 == FinFunction([1,1,3], 3)

# Same thing as a colimit of a general free diagram.
diagram = FreeDiagram([FinSet(2),FinSet(2),FinSet(3)],[(f,1,2),(g,1,3)])
colim = colimit(diagram)
@test ob(colim) == FinSet(3)
_, ι1, ι2 = colim
@test ι1 == FinFunction([1,2], 3)
@test ι2 == FinFunction([1,1,3], 3)

# Same thing as a colimit of a bipartite free diagram.
diagram = BipartiteFreeDiagram([FinSet(2)], [FinSet(2),FinSet(3)],
                               [(f,1,1),(g,1,2)])
colim = colimit(diagram)
@test ob(colim) == FinSet(3)
ι1, ι2 = colim
@test ι1 == FinFunction([1,2], 3)
@test ι2 == FinFunction([1,1,3], 3)

# Pushout complements
#--------------------

f = FinFunction([1,3], 4)
g = FinFunction([1,2,5,6], 6)
@test can_pushout_complement(f, g)
h, k = pushout_complement(f, g)
@test f⋅g == h⋅k
colim = pushout(f,h)
@test ob(colim) == FinSet(6)
@test allunique(collect(copair(colim, g, k)))

# Identification condition failure.
f = FinFunction([1,3], 4)
g = FinFunction([1,2,2,3], 3)
@test !can_pushout_complement(f, g)
@test_throws ErrorException pushout_complement(f, g)

# Subsets
#########

X = FinSet(10)
A, B = SubFinSet(X, [1,2,5,6,8,9]), SubFinSet(X, [2,3,5,7,8])
@test ob(A) == X
A_pred = SubFinSet(Bool[1,1,0,0,1,1,0,1,1,0])
@test hom(A) == hom(A_pred)
@test FinSets.predicate(A) == FinSets.predicate(A_pred)

# Lattice of subsets.
@test A ∧ B |> force == SubFinSet(X, [2,5,8])
@test A ∨ B |> force == SubFinSet(X, [1,2,3,5,6,7,8,9])
@test ⊤(X) |> force == SubFinSet(X, 1:10)
@test ⊥(X) |> force == SubFinSet(X, 1:0)

for alg in (SubOpBoolean(), SubOpWithLimits())
  @test meet(A, B, alg) |> sort == SubFinSet(X, [2,5,8])
  @test join(A, B, alg) |> sort == SubFinSet(X, [1,2,3,5,6,7,8,9])
  @test top(X, alg) |> force == SubFinSet(X, 1:10)
  @test bottom(X, alg) |> force == SubFinSet(X, 1:0)
end

end