text:
### GenePop parsing ###
"""
    genepop(infile::String; digits::Int64 = 3, popsep::Any = "POP", numpops::Int64)
Load a Genepop format file into memory as a PopObj object.
### Arguments
- `infile` : path to Genepop file
- `digits` : number of digits denoting each allele
- `popsep` : word that separates populations in `infile` (default: "POP")
- `numpops` : number of populations in `infile` (used for checking parser)
- `marker` : "snp" (default) or "msat"
### File must follow standard Genepop formatting:
- First line is a comment (and skipped)
- Loci are listed after first line as one-per-line without commas or in single comma-separated row
- A line with a particular keyword (default "POP") must delimit populations
- File is tab or space delimted

## Example

 `waspsNY = genepop("wasp_hive.gen", digits = 3, popsep = "POP", numpops = 2);`

### Genepop file example:
Wasp populations in New York \n
Locus1 \n
Locus2 \n
Locus3 \n
POP \n
Oneida_01,  250230 564568 110100  \n
Oneida_02,  252238 568558 100120  \n
Oneida_03,  254230 564558 090100  \n
POP \n
Newcomb_01,  254230 564558 080100 \n
Newcomb_02,  000230 564558 090080 \n
Newcomb_03,  254230 000000 090100 \n
Newcomb_04,  254230 564000 090120 \n
"""
function genepop(
    infile::String;
    digits::Int64 = 3,
    popsep::Any = "POP",
    numpops::Int64,
    marker = "snp",
)
    println("\n", "Input File : ", abspath(infile))
    if lowercase(marker) == "snp"
        geno_type = Int8
    else
        geno_type = Int16
    end
    gpop = split(open(readlines, infile)[2:end], popsep)
    if length(gpop) - 1 != numpops
        error("incorrect number of populations detected, see docstring for formatting
            expected : $numpops
            detected : $(length(gpop)-1) ")
    end
    if length(gpop[1]) == 1     # loci horizontally stacked
        locinames = strip.(split(gpop[1] |> join, ",") |> Array{String,1})
        replace!(locinames, "." => "_")
    else                        # loci vertically stacked
        locinames = replace(gpop[1], "." => "_")
    end
    d = Dict(string(i) => [] for i in locinames)
    popid = []
    indnames = []
    for i = 2:length(gpop)
        append!(popid, fill(i - 1, length(gpop[i])))
        for j = 1:length(gpop[i])
            #println(Base.Threads.threadid()) ########
            phasedloci = []
            push!(indnames, split(strip(gpop[i][j]), r"\,|\t")[1])
            unphasedloci = split(strip(gpop[i][j]), r"\s|\t")[2:end] |>
                           Array{String,1}
            replace!(unphasedloci, "-9" => "0"^digits) #just in case -9 = missing
            for locus in unphasedloci
                phasedlocus = parse.(
                    geno_type,
                    [join(i) for i in Iterators.partition(locus, digits)],
                ) |> sort |> Tuple
                push!(phasedloci, phasedlocus)
            end
            for (loc, geno) in zip(locinames, phasedloci)
                push!(d[loc], geno)
            end
        end
    end
    ploidy = length.(d[locinames[1]])   # lazy finding of ploidy from single locus
    for (loc, ploid) in zip(locinames, ploidy)
        miss_geno = fill(0, ploid) |> Tuple
        msat_miss_geno = ("0")
        replace!(d[loc], miss_geno => missing)
        replace!(d[loc], msat_miss_geno => missing)
    end
    # typesafe genotype DataFrame
    loci_df = DataFrame([Symbol(i) => Array{Union{Tuple,Missing},1}(d[i]) for i in locinames])
    samples_df = DataFrame(
        name = string.(indnames),
        population = string.(popid),
        ploidy = Int8.(ploidy),
        longitude = fill(missing, length(indnames)),
        latitude = fill(missing, length(indnames)),
    )
    PopObj(samples_df, loci_df)
end

# Alternative line-by-line reader
## NOT EXPORTED
#=
function gpop2(infile::String; digits::Int64 = 3, popsep::Any = "POP", numpops::Int64)
    println("\n", "Input File : ", abspath(infile))
    popid = []
    indnames = []
    locinames = []
    d = Dict()
    linenum = 1
    popcount = 0
    open(infile) do file
        for ln in eachline(file)
            if popcount == 0
                if linenum == 1
                    linenum += 1
                    continue
                elseif linenum == 2
                    if occursin(",", ln) == true  #loci are horizontally stacked
                        append!(locinames, strip.(split(ln |> join, ",") |> Array{String,1}))
                        replace!(locinames, "." => "_")
                        linenum += 1
                        continue
                    else        # loci are vertically stacked
                        push!(locinames, ln)
                        linenum += 1
                        continue
                    end
                else
                    if ln == popsep
                        popcount += 1
                        continue
                    else
                        push!(locinames, ln)
                        continue
                    end
                end
            end
            if ln == popsep
                popcount += 1
                continue
            else
                phasedloci = []
                push!(indnames, split(strip(ln), r"\,|\t")[1])
                push!(popid, popcount)
                unphasedloci = split(strip(ln), r"\s|\t")[2:end] |> Array{String,1}
                for locus in unphasedloci
                    phasedlocus = parse.(Int16,[join(i) for i in Iterators.partition(locus,digits)])  |> sort |> Tuple
                    push!(phasedloci, phasedlocus)
                end
                if locinames[1] ∉ keys(d)
                    [d[i] = [] for i in locinames]
                end
                for (loc,geno) in zip(locinames, phasedloci)
                    push!(d[loc], geno)
                end
            end
        end
    end
    ploidy = length.(d[locinames[1]])   # lazy finding of ploidy from single locus
    for (loc, ploid) in zip(locinames, ploidy)
        miss_geno = fill(0,ploid) |> Tuple
        replace!(d[loc], miss_geno => missing)
    end
    loci_df = DataFrame([i = Array{Union{Tuple, Missing},1}(d[i]) for i in locinames])
    names!(loci_df, Symbol.(locinames))
    samples_df = DataFrame(name = string.(indnames),
                           population = popid,
                           ploidy = Int8.(ploidy),
                           longitude = fill(missing,length(indnames)),
                           latitude = fill(missing,length(indnames)))
    PopObj(samples_df, loci_df)
end
=#

### CSV parsing ###

"""
    csv(infile::String; delim::Union{Char,String,Regex}, digits::Int64 = 3, location::Bool = false)
Load a CSV-type file into memory as a PopObj object
### Arguments
- `infile` : path to CSV file
- `delim` : delimiter characters. default is comma (","), can be space (" "), tab ("\\t"), etc.
- `digits` : number of digits denoting each allele
- `marker` : "snp" (default) or "msat"
- `location` : decimal degrees longitude/latitude provided as values 3/4
### File formatting:
- Loci names must be first row
- Individuals names must be first value in row
- Population ID's must be second value in row
- [Optional] longitude (x) values third value in row, latitude (y) fourth

## Example
`lizardsCA = Read.csv("CA_lizards.csv", delim = ",", digits = 3);`

### Formatting example:
Locus1,Locus2,Locus3   \n
sierra_01,mountain,001001,002002,001001   \n
sierra_02,mountain,001001,001001,001002   \n
snbarb_03,coast,001001,001001,001002 \n
snbarb_02,coast,001001,001001,001001 \n
snbarb_03,coast,001002,001001,001001 \n
"""
function csv(
    infile::String;
    delim::Union{Char,String,Regex} = ",",
    digits::Int = 3,
    marker = "snp",
    location::Bool = false,
)
    println("\n", "Input File : ", abspath(infile))
    popid = []
    indnames = []
    locx = []
    locy = []
    locinames = []
    d = Dict()
    linenum = 1
    if lowercase(marker) == "snp"
        geno_type = Int8
    else
        geno_type = Int16
    end
    file = open(infile, "r") #do file
        for ln in eachline(file)
            if linenum == 1
                loci_raw = split(ln, delim)
                loci_safe = replace(loci_raw, "." => "_")
                append!(locinames, loci_safe)
                [d[string(i)] = [] for i in locinames]
                linenum += 1
                continue
            end
            if location == false
                tmp = split(ln, delim) |> Array{String,1}
                # phase genotypes by ploidy
                phasedloci = []
                for locus in tmp[3:end]
                    phasedlocus = parse.(
                        geno_type,
                        [join(i) for i in Iterators.partition(locus, digits)],
                    ) |> sort |> Tuple
                    push!(phasedloci, phasedlocus)
                end
                for (loc, geno) in zip(locinames, phasedloci)
                    push!(d[loc], geno)
                end
                push!(indnames, tmp[1])
                push!(popid, tmp[2])
                push!(locx, missing)
                push!(locy, missing)
            else
                tmp = split(ln, delim) |> Array{String,1}
                replace!(tmp, "-9" => "0"^digits) #just in case -9 = missing
                phasedloci = []
                for locus in tmp[5:end]
                    phasedlocus = parse.(
                        Int16,
                        [join(i) for i in Iterators.partition(locus, digits)],
                    ) |> sort |> Tuple
                    push!(phasedloci, phasedlocus)
                end
                for (loc, geno) in zip(locinames, phasedloci)
                    push!(d[loc], geno)
                end
                push!(indnames, tmp[1])
                push!(popid, tmp[2])
                push!(locx, parse.(Float64, tmp[3]))
                push!(locy, parse.(Float64, tmp[4]))
            end
        end
    close(file)
    ploidy = length.(d[locinames[1]])   # lazy finding of ploidy from single locus
    for (loc, ploid) in zip(locinames, ploidy)
        miss_geno = fill(0, ploid) |> Tuple
        msat_miss_geno = ("0")
        replace!(d[loc], miss_geno => missing)
        replace!(d[loc], msat_miss_geno => missing)
    end
    # typesafe genotype DataFrame
    loci_df = DataFrame([Symbol(i) => Array{Union{Tuple,Missing},1}(d[i]) for i in locinames])
    #names!(loci_df, Symbol.(locinames))
    samples_df = DataFrame(
        name = string.(indnames),
        population = string.(popid),
        ploidy = Int8.(ploidy),
        longitude = locx,
        latitude = locy,
    )
    PopObj(samples_df, loci_df)
end

### VCF parsing ###

"""
    vcf(infile::String)
Load a VCF file into memory as a PopObj object. Population and [optional]
location information need to be provided separately.
- `infile` : path to VCF file
"""
function vcf(infile::String)
    vcf_file = VCF.Reader(open(infile, "r"))
    # get sample names from header
    sample_names = header(vcf_file).sampleID
    # fill in pop/lat/long with missing
    population = fill(missing, length(sample_names))
    lat = fill(missing, length(sample_names))
    long = fill(missing, length(sample_names))
    ## array of genotypes
    # get loci names
    locinames = []
    d = Dict()
    # get genotypes
    for record in vcf_file
        chr_safe = replace(VCF.chrom(record), "." => "_")
        chr_safer = replace(chr_safe, "|" => "_")
        pos = VCF.pos(record) |> string
        push!(locinames, chr_safer*"_"*pos)
        geno_raw = [split(i, ('/', '|')) for i in VCF.genotype(record, :, "GT")] |> sort
        # change missing data "." to "-1"
        geno_corr_miss = [replace(i, "." => "-1") for i in geno_raw]
        # convert everything to an integer
        geno_int = [parse.(Int8, i) for i in geno_corr_miss]
        # add 1 to shift genos so 0 is 1 and -1 is 0
        geno_shift = [i .+ Int8(1) for i in geno_int]
        geno_final = [replace(i, 0 => missing) for i in geno_shift]
        geno_tuple = [Tuple(i) for i in geno_final]
        d[locinames[end]] = geno_tuple
    end
    ploidy = length.(d[locinames[1]])
    loci_df = DataFrame([i = Array{Union{Tuple, Missing},1}(d[i]) for i in locinames])
    names!(loci_df, Symbol.(locinames))
    samples_df = DataFrame(name = sample_names,
                           population = population,
                           ploidy = Int8.(ploidy),
                           latitude = lat,
                           longitude = long)
    PopObj(samples_df, loci_df)
end


text:
import PyPlot as Plt
import Seaborn as Sns
using DataFrames

function plot_precision_recall(pr_dfs::T where T <: AbstractDict{String, DataFrame}; col::Symbol, color_per_label::Dict{String, String}, xlabel::String,
        ax=Plt.gca(), legend::Bool=true, bins=0.0:0.05:1.0, lw=2, alpha=0.1, legend_title::String="Target", legend_loc::String="best")
    for (k, df) in pr_dfs
        color = color_per_label[k]
        ax.hist(df[!,col], bins=bins, alpha=alpha, color=color);
        ax.hist(df[!,col], bins=bins, label=k, histtype="step", color=color, lw=lw);
    end

    ax.grid(true, alpha=0.1)
    ax.set_ylabel("Num. of cells");
    ax.set_xlim(0, 1);
    ax.set_xlabel(xlabel)
    if legend
        ax.legend(title=legend_title, labelspacing=0.2, handlelength=1, loc=legend_loc);
    end
    return ax
end

function plot_correlation_violins(datasets::NamedTuple, part_cor_key::Symbol, labels::Tuple{String, String}; color_per_label::Dict{String, String}, ax=Plt.gca())
    p_df = vcat([vcat([DataFrame(:Correlation => d[part_cor_key][i][1], :Type => d[:name], :Source => t) for d in datasets if (part_cor_key in keys(d)) && (d[:name] != "ISS")]...)
        for (i,t) in enumerate(labels)]...);
    p_df = @orderby(p_df, lowercase.(:Type))

    Sns.violinplot(x=p_df.Type, y=p_df.Correlation, hue=p_df.Source, split=true, inner="quart", palette=color_per_label, bw=0.2, scale="count", saturation=1, ax=ax)
    ax.legend(title="Source", loc="lower right", frameon=false, labelspacing=0.25)
    ax.set_xticklabels(ax.get_xticklabels(), fontsize=11)
    ax.set_ylim(-1.0, 1);
    ax.set_yticks(-1.0:0.5:1.0)
    ax.set_ylabel("Correlation between parts");
    return ax
end

function plot_cell_type_clustermap(expr_per_clust, gene_names::Vector{String}, gene_subs::Vector{Int}; title::String="",
        figsize::Tuple{<:Real, <:Real}=(4, 4), width_frac::Float64=0.5, kwargs...)
    Plt.figure(figsize=(figsize[1] * width_frac, figsize[2]))
    plt, cell_ord, gene_ord = clustermap(expr_per_clust, gene_names; gene_ord=gene_subs, transpose=true, cmap="OrRd", kwargs...);
    ax = plt.axes
    ax.tick_params(axis="both", which="both",length=0)
    ax.set_xticks([]); ax.set_yticklabels(ax.get_yticklabels(), fontsize=13, rotation=0, va="baseline");
    ax.set_title(title)
    return plt
end

ext: jl

max_stars_repo_name: hiraksarkar/BaysorAnalysis

max_issues_repo_name: hiraksarkar/BaysorAnalysis

max_forks_repo_name: hiraksarkar/BaysorAnalysis


text:
@doc raw"""
    Torus{N} <: AbstractPowerManifold

The n-dimensional torus is the $n$-dimensional product of the [`Circle`](@ref).

The [`Circle`](@ref) is stored internally within `M.manifold`, such that all functions of
[`AbstractPowerManifold`](@ref) can be used directly.
"""
struct Torus{N} <: AbstractPowerManifold{ℝ,Circle{ℝ},ArrayPowerRepresentation}
    manifold::Circle{ℝ}
end

Torus(n::Int) = Torus{n}(Circle())

Base.:^(M::Circle, n::Int) = Torus{n}(M)

@doc raw"""
    check_manifold_point(M::Torus{n},p)

Checks whether `p` is a valid point on the [`Torus`](@ref) `M`, i.e. each of
its entries is a valid point on the [`Circle`](@ref) and the length of `x` is `n`.
"""
check_manifold_point(::Torus, ::Any)
function check_manifold_point(M::Torus{N}, p; kwargs...) where {N}
    if length(p) != N
        return DomainError(
            length(p),
            "The number of elements in `p` ($(length(p))) does not match the dimension of the torus ($(N)).",
        )
    end
    return check_manifold_point(PowerManifold(M.manifold, N), p; kwargs...)
end
@doc raw"""
    check_tangent_vector(M::Torus{n}, p, X; check_base_point = true, kwargs...)

Checks whether `X` is a valid tangent vector to `p` on the [`Torus`](@ref) `M`.
This means, that `p` is valid, that `X` is of correct dimension and elementwise
a tangent vector to the elements of `p` on the [`Circle`](@ref).
The optional parameter `check_base_point` indicates, whether to call [`check_manifold_point`](@ref)  for `p`.
"""
function check_tangent_vector(
    M::Torus{N},
    p,
    X;
    check_base_point = true,
    kwargs...,
) where {N}
    if check_base_point && length(p) != N
        return DomainError(
            length(p),
            "The number of elements in `p` ($(length(p))) does not match the dimension of the torus ($(N)).",
        )
    end
    if length(X) != N
        return DomainError(
            length(X),
            "The number of elements in `X` ($(length(X))) does not match the dimension of the torus ($(N)).",
        )
    end
    return check_tangent_vector(
        PowerManifold(M.manifold, N),
        p,
        X;
        check_base_point = check_base_point,
        kwargs...,
    )
end

get_iterator(M::Torus{N}) where {N} = 1:N

@generated manifold_dimension(::Torus{N}) where {N} = N

@generated representation_size(::Torus{N}) where {N} = (N,)

Base.show(io::IO, ::Torus{N}) where {N} = print(io, "Torus($(N))")


text:
using LinearAlgebraicRepresentation, ViewerGL, SparseArrays
Lar = LinearAlgebraicRepresentation; GL = ViewerGL

using Revise, OhMyREPL

wall = Lar.svg2lar("test/svg/assembly/wall.svg", flag=false)
openings = Lar.svg2lar("test/svg/assembly/openings.svg", flag=false)
rectangle = Lar.svg2lar("test/svg/assembly/rectangle.svg", flag=false)
box = Lar.svg2lar("test/svg/assembly/box.svg", flag=false)

assembly = Lar.Struct([wall , openings, rectangle, box])
V,EV = Lar.struct2lar(assembly)
GL.VIEW([ GL.GLFrame2, GL.GLGrid(V,EV, GL.COLORS[1],1) ]);
W, copEV, copFE, boolmatrix = Lar.bool2d(assembly)
@show Matrix(boolmatrix)

innerpoints = Lar.internalpoints2d(W,copEV,copFE[1:end,:])
points = convert(Array{Float64,2}, hcat(innerpoints...)')
GL.VIEW([ GL.GLFrame2, GL.GLGrid(V,EV, GL.COLORS[1],1), GL.GLPoints(points) ]);

A = boolmatrix[:,1]
B = boolmatrix[:,2]
C = boolmatrix[:,3]
D = boolmatrix[:,4]
expr = ( D .& .!(.&(A, .!B, .!C)) )# A - B - C
formula = Matrix(copFE)' * Int.(expr) # coord vector of Faces
V = convert(Lar.Points,W')
EV = Lar.cop2lar(copEV)
EVmin = [ev for (k,ev) in enumerate(EV) if abs(formula[k])==1 ]
V = Lar.normalize(V)
V = Lar.normalize(V)
GL.VIEW([ GL.GLFrame2, GL.GLGrid(V,EVmin, GL.COLORS[1],1) ]);

VV = [[k] for k in 1:size(V,2)];
model = (V,[VV, EV])
GL.VIEW( GL.numbering(.1)( model,GL.COLORS[1],0.1 ) );


text:
@testset "8.6 Gamma functions" begin
    (a, b, c, d, e, m, n, p, x, ) = @variables a b c d e m n p x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integration*Problems*Involving*Gamma*Functions=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(d*x)^m*Gamma(n, b*x)=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*Gamma(n, b*x)=#

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

    @test_int [Gamma(0, a*x)*x^100, x, 1, (1/101)*x^101*Gamma(0, a*x) - Gamma(101, a*x)/(101*a^101)]

    @test_int [Gamma(0, a*x)*x^2, x, 1, (1/3)*x^3*Gamma(0, a*x) - Gamma(3, a*x)/(3*a^3)]
    @test_int [Gamma(0, a*x)*x^1, x, 1, (1/2)*x^2*Gamma(0, a*x) - Gamma(2, a*x)/(2*a^2)]
    @test_int [Gamma(0, a*x)*x^0, x, 1, -(1/(ℯ^(a*x)*a)) + x*Gamma(0, a*x)]
    @test_int [Gamma(0, a*x)/x^1, x, 1, a*x*HypergeometricPFQ([1, 1, 1], [2, 2, 2], (-a)*x) - EulerGamma*log(x) - (1/2)*log(a*x)^2]
    @test_int [Gamma(0, a*x)/x^2, x, 1, a*Gamma(-1, a*x) - Gamma(0, a*x)/x]
    @test_int [Gamma(0, a*x)/x^3, x, 1, (1/2)*a^2*Gamma(-2, a*x) - Gamma(0, a*x)/(2*x^2)]
    @test_int [Gamma(0, a*x)/x^4, x, 1, (1/3)*a^3*Gamma(-3, a*x) - Gamma(0, a*x)/(3*x^3)]

    @test_int [Gamma(1, a*x)*x^3, x, 4, -(6/(ℯ^(a*x)*a^4)) - (6*x)/(ℯ^(a*x)*a^3) - (3*x^2)/(ℯ^(a*x)*a^2) - x^3/(ℯ^(a*x)*a)]
    @test_int [Gamma(1, a*x)*x^2, x, 3, -(2/(ℯ^(a*x)*a^3)) - (2*x)/(ℯ^(a*x)*a^2) - x^2/(ℯ^(a*x)*a)]
    @test_int [Gamma(1, a*x)*x^1, x, 2, -(1/(ℯ^(a*x)*a^2)) - x/(ℯ^(a*x)*a)]
    @test_int [Gamma(1, a*x)*x^0, x, 1, -(1/(ℯ^(a*x)*a))]
    @test_int [Gamma(1, a*x)/x^1, x, 1, ExpIntegralEi((-a)*x)]
    @test_int [Gamma(1, a*x)/x^2, x, 2, -(1/(ℯ^(a*x)*x)) - a*ExpIntegralEi((-a)*x)]
    @test_int [Gamma(1, a*x)/x^3, x, 3, -(1/(ℯ^(a*x)*(2*x^2))) + a/(ℯ^(a*x)*(2*x)) + (1/2)*a^2*ExpIntegralEi((-a)*x)]
    @test_int [Gamma(1, a*x)/x^4, x, 4, -(1/(ℯ^(a*x)*(3*x^3))) + a/(ℯ^(a*x)*(6*x^2)) - a^2/(ℯ^(a*x)*(6*x)) - (1/6)*a^3*ExpIntegralEi((-a)*x)]

    @test_int [Gamma(2, a*x)*x^100, x, 1, (1/101)*x^101*Gamma(2, a*x) - Gamma(103, a*x)/(101*a^101)]

    @test_int [Gamma(2, a*x)*x^2, x, 1, (1/3)*x^3*Gamma(2, a*x) - Gamma(5, a*x)/(3*a^3)]
    @test_int [Gamma(2, a*x)*x^1, x, 1, (1/2)*x^2*Gamma(2, a*x) - Gamma(4, a*x)/(2*a^2)]
    @test_int [Gamma(2, a*x)*x^0, x, 1, x*Gamma(2, a*x) - Gamma(3, a*x)/a]
    @test_int [Gamma(2, a*x)/x^1, x, 2, -ℯ^(-a*x) + ExpIntegralEi(-a*x)]
    @test_int [Gamma(2, a*x)/x^2, x, 1, a/ℯ^(a*x) - Gamma(2, a*x)/x]
    @test_int [Gamma(2, a*x)/x^3, x, 1, (1/2)*a^2*Gamma(0, a*x) - Gamma(2, a*x)/(2*x^2)]
    @test_int [Gamma(2, a*x)/x^4, x, 1, (1/3)*a^3*Gamma(-1, a*x) - Gamma(2, a*x)/(3*x^3)]

    @test_int [Gamma(3, a*x)*x^100, x, 1, (1/101)*x^101*Gamma(3, a*x) - Gamma(104, a*x)/(101*a^101)]

    @test_int [Gamma(3, a*x)*x^2, x, 1, (1/3)*x^3*Gamma(3, a*x) - Gamma(6, a*x)/(3*a^3)]
    @test_int [Gamma(3, a*x)*x^1, x, 1, (1/2)*x^2*Gamma(3, a*x) - Gamma(5, a*x)/(2*a^2)]
    @test_int [Gamma(3, a*x)*x^0, x, 1, x*Gamma(3, a*x) - Gamma(4, a*x)/a]
    @test_int [Gamma(3, a*x)/x^1, x, 3, -2/ℯ^(a*x) + 2*ExpIntegralEi((-a)*x) - Gamma(2, a*x)]
    @test_int [Gamma(3, a*x)/x^2, x, 1, a*Gamma(2, a*x) - Gamma(3, a*x)/x]
    @test_int [Gamma(3, a*x)/x^3, x, 1, ((1/2)*a^2)/ℯ^(a*x) - Gamma(3, a*x)/(2*x^2)]
    @test_int [Gamma(3, a*x)/x^4, x, 1, (1/3)*a^3*Gamma(0, a*x) - Gamma(3, a*x)/(3*x^3)]

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

    @test_int [Gamma(-1, a*x)*x^100, x, 1, (1/101)*x^101*Gamma(-1, a*x) - Gamma(100, a*x)/(101*a^101)]

    @test_int [Gamma(-1, a*x)*x^3, x, 1, (1/4)*x^4*Gamma(-1, a*x) - Gamma(3, a*x)/(4*a^4)]
    @test_int [Gamma(-1, a*x)*x^2, x, 1, (1/3)*x^3*Gamma(-1, a*x) - Gamma(2, a*x)/(3*a^3)]
    @test_int [Gamma(-1, a*x)*x^1, x, 1, -(1/(ℯ^(a*x)*(2*a^2))) + (1/2)*x^2*Gamma(-1, a*x)]
    @test_int [Gamma(-1, a*x)*x^0, x, 1, x*Gamma(-1, a*x) - Gamma(0, a*x)/a]
    @test_int [Gamma(-1, a*x)/x^1, x, 2, -Gamma(-1, a*x) - a*x*HypergeometricPFQ([1, 1, 1], [2, 2, 2], (-a)*x) + EulerGamma*log(x) + (1/2)*log(a*x)^2]
    @test_int [Gamma(-1, a*x)/x^2, x, 1, a*Gamma(-2, a*x) - Gamma(-1, a*x)/x]
    @test_int [Gamma(-1, a*x)/x^3, x, 1, (1/2)*a^2*Gamma(-3, a*x) - Gamma(-1, a*x)/(2*x^2)]
    @test_int [Gamma(-1, a*x)/x^4, x, 1, (1/3)*a^3*Gamma(-4, a*x) - Gamma(-1, a*x)/(3*x^3)]

    @test_int [Gamma(-2, a*x)*x^100, x, 1, (1/101)*x^101*Gamma(-2, a*x) - Gamma(99, a*x)/(101*a^101)]

    @test_int [Gamma(-2, a*x)*x^3, x, 1, (1/4)*x^4*Gamma(-2, a*x) - Gamma(2, a*x)/(4*a^4)]
    @test_int [Gamma(-2, a*x)*x^2, x, 1, -(1/(ℯ^(a*x)*(3*a^3))) + (1/3)*x^3*Gamma(-2, a*x)]
    @test_int [Gamma(-2, a*x)*x^1, x, 1, (1/2)*x^2*Gamma(-2, a*x) - Gamma(0, a*x)/(2*a^2)]
    @test_int [Gamma(-2, a*x)*x^0, x, 1, x*Gamma(-2, a*x) - Gamma(-1, a*x)/a]
    @test_int [Gamma(-2, a*x)/x^1, x, 3, (-(1/2))*Gamma(-2, a*x) + (1/2)*Gamma(-1, a*x) + (1/2)*a*x*HypergeometricPFQ([1, 1, 1], [2, 2, 2], (-a)*x) - (1/2)*EulerGamma*log(x) - (1/4)*log(a*x)^2]
    @test_int [Gamma(-2, a*x)/x^2, x, 1, a*Gamma(-3, a*x) - Gamma(-2, a*x)/x]
    @test_int [Gamma(-2, a*x)/x^3, x, 1, (1/2)*a^2*Gamma(-4, a*x) - Gamma(-2, a*x)/(2*x^2)]
    @test_int [Gamma(-2, a*x)/x^4, x, 1, (1/3)*a^3*Gamma(-5, a*x) - Gamma(-2, a*x)/(3*x^3)]

    @test_int [Gamma(-3, a*x)*x^100, x, 1, (1/101)*x^101*Gamma(-3, a*x) - Gamma(98, a*x)/(101*a^101)]

    @test_int [Gamma(-3, a*x)*x^3, x, 1, -(1/(ℯ^(a*x)*(4*a^4))) + (1/4)*x^4*Gamma(-3, a*x)]
    @test_int [Gamma(-3, a*x)*x^2, x, 1, (1/3)*x^3*Gamma(-3, a*x) - Gamma(0, a*x)/(3*a^3)]
    @test_int [Gamma(-3, a*x)*x^1, x, 1, (1/2)*x^2*Gamma(-3, a*x) - Gamma(-1, a*x)/(2*a^2)]
    @test_int [Gamma(-3, a*x)*x^0, x, 1, x*Gamma(-3, a*x) - Gamma(-2, a*x)/a]
    @test_int [Gamma(-3, a*x)/x^1, x, 4, (-(1/3))*Gamma(-3, a*x) + (1/6)*Gamma(-2, a*x) - (1/6)*Gamma(-1, a*x) - (1/6)*a*x*HypergeometricPFQ([1, 1, 1], [2, 2, 2], (-a)*x) + (1/6)*EulerGamma*log(x) + (1/12)*log(a*x)^2]
    @test_int [Gamma(-3, a*x)/x^2, x, 1, a*Gamma(-4, a*x) - Gamma(-3, a*x)/x]
    @test_int [Gamma(-3, a*x)/x^3, x, 1, (1/2)*a^2*Gamma(-5, a*x) - Gamma(-3, a*x)/(2*x^2)]
    @test_int [Gamma(-3, a*x)/x^4, x, 1, (1/3)*a^3*Gamma(-6, a*x) - Gamma(-3, a*x)/(3*x^3)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*Gamma(n/2, b*x)=#

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

    @test_int [Gamma(1/2, a*x)*x^100, x, 1, (1/101)*x^101*Gamma(1/2, a*x) - Gamma(203/2, a*x)/(101*a^101)]

    @test_int [Gamma(1/2, a*x)*x^2, x, 1, (1/3)*x^3*Gamma(1/2, a*x) - Gamma(7/2, a*x)/(3*a^3)]
    @test_int [Gamma(1/2, a*x)*x^1, x, 1, (1/2)*x^2*Gamma(1/2, a*x) - Gamma(5/2, a*x)/(2*a^2)]
    @test_int [Gamma(1/2, a*x)*x^0, x, 1, x*Gamma(1/2, a*x) - Gamma(3/2, a*x)/a]
    @test_int [Gamma(1/2, a*x)/x^1, x, 1, -4*sqrt(a*x)*HypergeometricPFQ([1/2, 1/2], [3/2, 3/2], (-a)*x) + sqrt(pi)*log(x)]
    @test_int [Gamma(1/2, a*x)/x^2, x, 1, a*Gamma(-(1/2), a*x) - Gamma(1/2, a*x)/x]
    @test_int [Gamma(1/2, a*x)/x^3, x, 1, (1/2)*a^2*Gamma(-(3/2), a*x) - Gamma(1/2, a*x)/(2*x^2)]
    @test_int [Gamma(1/2, a*x)/x^4, x, 1, (1/3)*a^3*Gamma(-(5/2), a*x) - Gamma(1/2, a*x)/(3*x^3)]

    @test_int [Gamma(3/2, a*x)*x^100, x, 1, (1/101)*x^101*Gamma(3/2, a*x) - Gamma(205/2, a*x)/(101*a^101)]

    @test_int [Gamma(3/2, a*x)*x^2, x, 1, (1/3)*x^3*Gamma(3/2, a*x) - Gamma(9/2, a*x)/(3*a^3)]
    @test_int [Gamma(3/2, a*x)*x^1, x, 1, (1/2)*x^2*Gamma(3/2, a*x) - Gamma(7/2, a*x)/(2*a^2)]
    @test_int [Gamma(3/2, a*x)*x^0, x, 1, x*Gamma(3/2, a*x) - Gamma(5/2, a*x)/a]
    @test_int [Gamma(3/2, a*x)/x^1, x, 1, (-(4/9))*(a*x)^(3/2)*HypergeometricPFQ([3/2, 3/2], [5/2, 5/2], (-a)*x) + (1/2)*sqrt(pi)*log(x)]
    @test_int [Gamma(3/2, a*x)/x^2, x, 1, a*Gamma(1/2, a*x) - Gamma(3/2, a*x)/x]
    @test_int [Gamma(3/2, a*x)/x^3, x, 1, (1/2)*a^2*Gamma(-(1/2), a*x) - Gamma(3/2, a*x)/(2*x^2)]
    @test_int [Gamma(3/2, a*x)/x^4, x, 1, (1/3)*a^3*Gamma(-(3/2), a*x) - Gamma(3/2, a*x)/(3*x^3)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(d*x)^m*Gamma(n, b*x)*when*m*symbolic=#

    @test_int [(d*x)^m*Gamma(3, b*x), x, 1, ((d*x)^(1 + m)*Gamma(3, b*x))/(d*(1 + m)) - ((d*x)^m*Gamma(4 + m, b*x))/((b*x)^m*(b*(1 + m)))]
    @test_int [(d*x)^m*Gamma(2, b*x), x, 1, ((d*x)^(1 + m)*Gamma(2, b*x))/(d*(1 + m)) - ((d*x)^m*Gamma(3 + m, b*x))/((b*x)^m*(b*(1 + m)))]
    @test_int [(d*x)^m*Gamma(1, b*x), x, 1, -(((d*x)^m*Gamma(1 + m, b*x))/((b*x)^m*b))]
    @test_int [(d*x)^m*Gamma(0, b*x), x, 1, ((d*x)^(1 + m)*Gamma(0, b*x))/(d*(1 + m)) - ((d*x)^m*Gamma(1 + m, b*x))/((b*x)^m*(b*(1 + m)))]
    @test_int [(d*x)^m*Gamma(-1, b*x), x, 1, ((d*x)^(1 + m)*Gamma(-1, b*x))/(d*(1 + m)) - ((d*x)^m*Gamma(m, b*x))/((b*x)^m*(b*(1 + m)))]
    @test_int [(d*x)^m*Gamma(-2, b*x), x, 1, ((d*x)^(1 + m)*Gamma(-2, b*x))/(d*(1 + m)) - ((d*x)^m*Gamma(-1 + m, b*x))/((b*x)^m*(b*(1 + m)))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(d*x)^m*Gamma(n, b*x)*when*n*symbolic=#

    @test_int [x^m*Gamma(n, x), x, 1, (x^(1 + m)*Gamma(n, x))/(1 + m) - Gamma(1 + m + n, x)/(1 + m)]
    @test_int [x^m*Gamma(n, b*x), x, 1, (x^(1 + m)*Gamma(n, b*x))/(1 + m) - (x^m*Gamma(1 + m + n, b*x))/((b*x)^m*(b*(1 + m)))]
    @test_int [(d*x)^m*Gamma(n, x), x, 1, ((d*x)^(1 + m)*Gamma(n, x))/(d*(1 + m)) - ((d*x)^m*Gamma(1 + m + n, x))/(x^m*(1 + m))]
    @test_int [(d*x)^m*Gamma(n, b*x), x, 1, ((d*x)^(1 + m)*Gamma(n, b*x))/(d*(1 + m)) - ((d*x)^m*Gamma(1 + m + n, b*x))/((b*x)^m*(b*(1 + m)))]

    @test_int [Gamma(n, a*x)*x^100, x, 1, (1/101)*x^101*Gamma(n, a*x) - Gamma(101 + n, a*x)/(101*a^101)]

    @test_int [Gamma(n, a*x)*x^2, x, 1, (1/3)*x^3*Gamma(n, a*x) - Gamma(3 + n, a*x)/(3*a^3)]
    @test_int [Gamma(n, a*x)*x^1, x, 1, (1/2)*x^2*Gamma(n, a*x) - Gamma(2 + n, a*x)/(2*a^2)]
    @test_int [Gamma(n, a*x)*x^0, x, 1, x*Gamma(n, a*x) - Gamma(1 + n, a*x)/a]
    @test_int [Gamma(n, a*x)/x^1, x, 1, -(((a*x)^n*HypergeometricPFQ([n, n], [1 + n, 1 + n], (-a)*x))/n^2) + Gamma(n)*log(x)]
    @test_int [Gamma(n, a*x)/x^2, x, 1, a*Gamma(-1 + n, a*x) - Gamma(n, a*x)/x]
    @test_int [Gamma(n, a*x)/x^3, x, 1, (1/2)*a^2*Gamma(-2 + n, a*x) - Gamma(n, a*x)/(2*x^2)]
    @test_int [Gamma(n, a*x)/x^4, x, 1, (1/3)*a^3*Gamma(-3 + n, a*x) - Gamma(n, a*x)/(3*x^3)]

    @test_int [Gamma(n, 2*x)*x^100, x, 1, (1/101)*x^101*Gamma(n, 2*x) - Gamma(101 + n, 2*x)/256065421246102339102334047485952]

    @test_int [Gamma(n, 2*x)*x^2, x, 1, (1/3)*x^3*Gamma(n, 2*x) - (1/24)*Gamma(3 + n, 2*x)]
    @test_int [Gamma(n, 2*x)*x^1, x, 1, (1/2)*x^2*Gamma(n, 2*x) - (1/8)*Gamma(2 + n, 2*x)]
    @test_int [Gamma(n, 2*x)*x^0, x, 1, x*Gamma(n, 2*x) - (1/2)*Gamma(1 + n, 2*x)]
    @test_int [Gamma(n, 2*x)/x^1, x, 1, -((2^n*x^n*HypergeometricPFQ([n, n], [1 + n, 1 + n], -2*x))/n^2) + Gamma(n)*log(x)]
    @test_int [Gamma(n, 2*x)/x^2, x, 1, 2*Gamma(-1 + n, 2*x) - Gamma(n, 2*x)/x]
    @test_int [Gamma(n, 2*x)/x^3, x, 1, 2*Gamma(-2 + n, 2*x) - Gamma(n, 2*x)/(2*x^2)]
    @test_int [Gamma(n, 2*x)/x^4, x, 1, (8/3)*Gamma(-3 + n, 2*x) - Gamma(n, 2*x)/(3*x^3)]

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

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

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

    @test_int [(c+ d*x)^3*Gamma(0, a + b*x), x, 8, -(((b*c - a*d)^3*ℯ^(-a - b*x))/(4*b^4)) - ((b*c - a*d)^4*Gamma(0, a + b*x))/(4*b^4*d) + ((c + d*x)^4*Gamma(0, a + b*x))/(4*d) - (d*(b*c - a*d)^2*ℯ^(-a + (b*c)/d)*Gamma(2, (b*(c + d*x))/d))/(4*b^4) - (d^2*(b*c - a*d)*ℯ^(-a + (b*c)/d)*Gamma(3, (b*(c + d*x))/d))/(4*b^4) - (d^3*ℯ^(-a + (b*c)/d)*Gamma(4, (b*(c + d*x))/d))/(4*b^4)]
    @test_int [(c+ d*x)^2*Gamma(0, a + b*x), x, 7, -(((b*c - a*d)^2*ℯ^(-a - b*x))/(3*b^3)) - ((b*c - a*d)^3*Gamma(0, a + b*x))/(3*b^3*d) + ((c + d*x)^3*Gamma(0, a + b*x))/(3*d) - (d*(b*c - a*d)*ℯ^(-a + (b*c)/d)*Gamma(2, (b*(c + d*x))/d))/(3*b^3) - (d^2*ℯ^(-a + (b*c)/d)*Gamma(3, (b*(c + d*x))/d))/(3*b^3)]
    @test_int [(c+ d*x)^1*Gamma(0, a + b*x), x, 6, -(((b*c - a*d)*ℯ^(-a - b*x))/(2*b^2)) - ((b*c - a*d)^2*Gamma(0, a + b*x))/(2*b^2*d) + ((c + d*x)^2*Gamma(0, a + b*x))/(2*d) - (d*ℯ^(-a + (b*c)/d)*Gamma(2, (b*(c + d*x))/d))/(2*b^2)]
    @test_int [(c+ d*x)^0*Gamma(0, a + b*x), x, 1, -(ℯ^(-a - b*x)/b) + ((a + b*x)*Gamma(0, a + b*x))/b]
    @test_int [Gamma(0, a + b*x)/(c+ d*x)^1, x, 0, Unintegrable(Gamma(0, a + b*x)/(c + d*x), x)]
    @test_int [Gamma(0, a + b*x)/(c+ d*x)^2, x, 5, (b*Gamma(0, a + b*x))/(d*(b*c - a*d)) - Gamma(0, a + b*x)/(d*(c + d*x)) - (b*ℯ^(-a + (b*c)/d)*Gamma(0, (b*(c + d*x))/d))/(d*(b*c - a*d))]
    @test_int [Gamma(0, a + b*x)/(c+ d*x)^3, x, 6, -((b^2*ℯ^(-a + (b*c)/d)*Gamma(-1, (b*(c + d*x))/d))/(2*d^2*(b*c - a*d))) + (b^2*Gamma(0, a + b*x))/(2*d*(b*c - a*d)^2) - Gamma(0, a + b*x)/(2*d*(c + d*x)^2) - (b^2*ℯ^(-a + (b*c)/d)*Gamma(0, (b*(c + d*x))/d))/(2*d*(b*c - a*d)^2)]
    @test_int [Gamma(0, a + b*x)/(c+ d*x)^4, x, 7, -((b^3*ℯ^(-a + (b*c)/d)*Gamma(-2, (b*(c + d*x))/d))/(3*d^3*(b*c - a*d))) - (b^3*ℯ^(-a + (b*c)/d)*Gamma(-1, (b*(c + d*x))/d))/(3*d^2*(b*c - a*d)^2) + (b^3*Gamma(0, a + b*x))/(3*d*(b*c - a*d)^3) - Gamma(0, a + b*x)/(3*d*(c + d*x)^3) - (b^3*ℯ^(-a + (b*c)/d)*Gamma(0, (b*(c + d*x))/d))/(3*d*(b*c - a*d)^3)]

    @test_int [(c+ d*x)^4*Gamma(1, a + b*x), x, 5, -((24*d^4*ℯ^(-a - b*x))/b^5) - (24*d^3*ℯ^(-a - b*x)*(c + d*x))/b^4 - (12*d^2*ℯ^(-a - b*x)*(c + d*x)^2)/b^3 - (4*d*ℯ^(-a - b*x)*(c + d*x)^3)/b^2 - (ℯ^(-a - b*x)*(c + d*x)^4)/b]
    @test_int [(c+ d*x)^3*Gamma(1, a + b*x), x, 4, -((6*d^3*ℯ^(-a - b*x))/b^4) - (6*d^2*ℯ^(-a - b*x)*(c + d*x))/b^3 - (3*d*ℯ^(-a - b*x)*(c + d*x)^2)/b^2 - (ℯ^(-a - b*x)*(c + d*x)^3)/b]
    @test_int [(c+ d*x)^2*Gamma(1, a + b*x), x, 3, -((2*d^2*ℯ^(-a - b*x))/b^3) - (2*d*ℯ^(-a - b*x)*(c + d*x))/b^2 - (ℯ^(-a - b*x)*(c + d*x)^2)/b]
    @test_int [(c+ d*x)^1*Gamma(1, a + b*x), x, 2, -((d*ℯ^(-a - b*x))/b^2) - (ℯ^(-a - b*x)*(c + d*x))/b]
    @test_int [(c+ d*x)^0*Gamma(1, a + b*x), x, 1, -(ℯ^(-a - b*x)/b)]
    @test_int [Gamma(1, a + b*x)/(c+ d*x)^1, x, 1, (ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/d]
    @test_int [Gamma(1, a + b*x)/(c+ d*x)^2, x, 2, -(ℯ^(-a - b*x)/(d*(c + d*x))) - (b*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/d^2]
    @test_int [Gamma(1, a + b*x)/(c+ d*x)^3, x, 3, -(ℯ^(-a - b*x)/(2*d*(c + d*x)^2)) + (b*ℯ^(-a - b*x))/(2*d^2*(c + d*x)) + (b^2*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/(2*d^3)]
    @test_int [Gamma(1, a + b*x)/(c+ d*x)^4, x, 4, -(ℯ^(-a - b*x)/(3*d*(c + d*x)^3)) + (b*ℯ^(-a - b*x))/(6*d^2*(c + d*x)^2) - (b^2*ℯ^(-a - b*x))/(6*d^3*(c + d*x)) - (b^3*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/(6*d^4)]

    @test_int [(c+ d*x)^3*Gamma(2, a + b*x), x, 5, ((c + d*x)^4*Gamma(2, a + b*x))/(4*d) + (d^2*(b*c - a*d)*ℯ^(-a + (b*c)/d)*Gamma(5, (b*(c + d*x))/d))/(4*b^4) - (d^3*ℯ^(-a + (b*c)/d)*Gamma(6, (b*(c + d*x))/d))/(4*b^4)]
    @test_int [(c+ d*x)^2*Gamma(2, a + b*x), x, 5, ((c + d*x)^3*Gamma(2, a + b*x))/(3*d) + (d*(b*c - a*d)*ℯ^(-a + (b*c)/d)*Gamma(4, (b*(c + d*x))/d))/(3*b^3) - (d^2*ℯ^(-a + (b*c)/d)*Gamma(5, (b*(c + d*x))/d))/(3*b^3)]
    @test_int [(c+ d*x)^1*Gamma(2, a + b*x), x, 5, ((c + d*x)^2*Gamma(2, a + b*x))/(2*d) + ((b*c - a*d)*ℯ^(-a + (b*c)/d)*Gamma(3, (b*(c + d*x))/d))/(2*b^2) - (d*ℯ^(-a + (b*c)/d)*Gamma(4, (b*(c + d*x))/d))/(2*b^2)]
    @test_int [(c+ d*x)^0*Gamma(2, a + b*x), x, 1, ((a + b*x)*Gamma(2, a + b*x))/b - Gamma(3, a + b*x)/b]
    @test_int [Gamma(2, a + b*x)/(c+ d*x)^1, x, 6, -(ℯ^(-a - b*x)/d) + (ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/d - ((b*c - a*d)*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/d^2]
    @test_int [Gamma(2, a + b*x)/(c+ d*x)^2, x, 5, (b*ℯ^(-a - b*x))/d^2 - (b*(b*c - a*d)*ℯ^(-a + (b*c)/d)*Gamma(0, (b*(c + d*x))/d))/d^3 - Gamma(2, a + b*x)/(d*(c + d*x))]
    @test_int [Gamma(2, a + b*x)/(c+ d*x)^3, x, 5, -((b^2*(b*c - a*d)*ℯ^(-a + (b*c)/d)*Gamma(-1, (b*(c + d*x))/d))/(2*d^4)) + (b^2*ℯ^(-a + (b*c)/d)*Gamma(0, (b*(c + d*x))/d))/(2*d^3) - Gamma(2, a + b*x)/(2*d*(c + d*x)^2)]
    @test_int [Gamma(2, a + b*x)/(c+ d*x)^4, x, 5, -((b^3*(b*c - a*d)*ℯ^(-a + (b*c)/d)*Gamma(-2, (b*(c + d*x))/d))/(3*d^5)) + (b^3*ℯ^(-a + (b*c)/d)*Gamma(-1, (b*(c + d*x))/d))/(3*d^4) - Gamma(2, a + b*x)/(3*d*(c + d*x)^3)]
    @test_int [Gamma(2, a + b*x)/(c+ d*x)^5, x, 5, -((b^4*(b*c - a*d)*ℯ^(-a + (b*c)/d)*Gamma(-3, (b*(c + d*x))/d))/(4*d^6)) + (b^4*ℯ^(-a + (b*c)/d)*Gamma(-2, (b*(c + d*x))/d))/(4*d^5) - Gamma(2, a + b*x)/(4*d*(c + d*x)^4)]

    @test_int [(c+ d*x)^3*Gamma(3, a + b*x), x, 6, ((c + d*x)^4*Gamma(3, a + b*x))/(4*d) - (d*(b*c - a*d)^2*ℯ^(-a + (b*c)/d)*Gamma(5, (b*(c + d*x))/d))/(4*b^4) + (d^2*(b*c - a*d)*ℯ^(-a + (b*c)/d)*Gamma(6, (b*(c + d*x))/d))/(2*b^4) - (d^3*ℯ^(-a + (b*c)/d)*Gamma(7, (b*(c + d*x))/d))/(4*b^4)]
    @test_int [(c+ d*x)^2*Gamma(3, a + b*x), x, 6, ((c + d*x)^3*Gamma(3, a + b*x))/(3*d) - ((b*c - a*d)^2*ℯ^(-a + (b*c)/d)*Gamma(4, (b*(c + d*x))/d))/(3*b^3) + (2*d*(b*c - a*d)*ℯ^(-a + (b*c)/d)*Gamma(5, (b*(c + d*x))/d))/(3*b^3) - (d^2*ℯ^(-a + (b*c)/d)*Gamma(6, (b*(c + d*x))/d))/(3*b^3)]
    @test_int [(c+ d*x)^1*Gamma(3, a + b*x), x, 6, -(((b*c - a*d)^2*Gamma(3, a + b*x))/(2*b^2*d)) + ((c + d*x)^2*Gamma(3, a + b*x))/(2*d) - ((b*c - a*d)*Gamma(4, a + b*x))/b^2 - (d*Gamma(5, a + b*x))/(2*b^2)]
    @test_int [(c+ d*x)^0*Gamma(3, a + b*x), x, 1, ((a + b*x)*Gamma(3, a + b*x))/b - Gamma(4, a + b*x)/b]
    @test_int [Gamma(3, a + b*x)/(c+ d*x)^1, x, 13, -((3*ℯ^(-a - b*x))/d) + ((b*c - a*d)*ℯ^(-a - b*x))/d^2 - (ℯ^(-a - b*x)*(a + b*x))/d + (2*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/d - (2*(b*c - a*d)*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/d^2 + ((b*c - a*d)^2*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/d^3]
    @test_int [Gamma(3, a + b*x)/(c+ d*x)^2, x, 6, -((b*(b*c - a*d)*ℯ^(-a - b*x))/d^3) + (b*(b*c - a*d)^2*ℯ^(-a + (b*c)/d)*Gamma(0, (b*(c + d*x))/d))/d^4 + (b*Gamma(2, a + b*x))/d^2 - Gamma(3, a + b*x)/(d*(c + d*x))]
    @test_int [Gamma(3, a + b*x)/(c+ d*x)^3, x, 6, (b^2*ℯ^(-a - b*x))/(2*d^3) + (b^2*(b*c - a*d)^2*ℯ^(-a + (b*c)/d)*Gamma(-1, (b*(c + d*x))/d))/(2*d^5) - (b^2*(b*c - a*d)*ℯ^(-a + (b*c)/d)*Gamma(0, (b*(c + d*x))/d))/d^4 - Gamma(3, a + b*x)/(2*d*(c + d*x)^2)]
    @test_int [Gamma(3, a + b*x)/(c+ d*x)^4, x, 6, (b^3*(b*c - a*d)^2*ℯ^(-a + (b*c)/d)*Gamma(-2, (b*(c + d*x))/d))/(3*d^6) - (2*b^3*(b*c - a*d)*ℯ^(-a + (b*c)/d)*Gamma(-1, (b*(c + d*x))/d))/(3*d^5) + (b^3*ℯ^(-a + (b*c)/d)*Gamma(0, (b*(c + d*x))/d))/(3*d^4) - Gamma(3, a + b*x)/(3*d*(c + d*x)^3)]
    @test_int [Gamma(3, a + b*x)/(c+ d*x)^5, x, 6, (b^4*(b*c - a*d)^2*ℯ^(-a + (b*c)/d)*Gamma(-3, (b*(c + d*x))/d))/(4*d^7) - (b^4*(b*c - a*d)*ℯ^(-a + (b*c)/d)*Gamma(-2, (b*(c + d*x))/d))/(2*d^6) + (b^4*ℯ^(-a + (b*c)/d)*Gamma(-1, (b*(c + d*x))/d))/(4*d^5) - Gamma(3, a + b*x)/(4*d*(c + d*x)^4)]

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

    @test_int [(c+ d*x)^3*Gamma(-1, a + b*x), x, 8, -((3*d*(b*c - a*d)^2*ℯ^(-a - b*x))/(2*b^4)) - ((b*c - a*d)^4*Gamma(-1, a + b*x))/(4*b^4*d) + ((c + d*x)^4*Gamma(-1, a + b*x))/(4*d) - ((b*c - a*d)^3*Gamma(0, a + b*x))/b^4 - (d^2*(b*c - a*d)*Gamma(2, a + b*x))/b^4 - (d^3*Gamma(3, a + b*x))/(4*b^4)]
    @test_int [(c+ d*x)^2*Gamma(-1, a + b*x), x, 7, -((d*(3*b*c - 2*a*d)*ℯ^(-a - b*x))/(3*b^3)) - ((b*c - a*d)^3*Gamma(-1, a + b*x))/(3*b^3*d) + ((c + d*x)^3*Gamma(-1, a + b*x))/(3*d) - ((b*c - a*d)^2*Gamma(0, a + b*x))/b^3 - (d^2*Gamma(2, b*x))/(ℯ^a*(3*b^3))]
    @test_int [(c+ d*x)^1*Gamma(-1, a + b*x), x, 6, -((d*ℯ^(-a - b*x))/(2*b^2)) - ((b*c - a*d)^2*Gamma(-1, a + b*x))/(2*b^2*d) + ((c + d*x)^2*Gamma(-1, a + b*x))/(2*d) - ((b*c - a*d)*Gamma(0, a + b*x))/b^2]
    @test_int [(c+ d*x)^0*Gamma(-1, a + b*x), x, 1, ((a + b*x)*Gamma(-1, a + b*x))/b - Gamma(0, a + b*x)/b]
    @test_int [Gamma(-1, a + b*x)/(c+ d*x)^1, x, 0, Unintegrable(Gamma(-1, a + b*x)/(c + d*x), x)]
    @test_int [Gamma(-1, a + b*x)/(c+ d*x)^2, x, 6, (b*Gamma(-1, a + b*x))/(d*(b*c - a*d)) - Gamma(-1, a + b*x)/(d*(c + d*x)) - (b*Gamma(0, a + b*x))/(b*c - a*d)^2 + (b*ℯ^(-a + (b*c)/d)*Gamma(0, (b*(c + d*x))/d))/(b*c - a*d)^2]
    @test_int [Gamma(-1, a + b*x)/(c+ d*x)^3, x, 7, (b^2*Gamma(-1, a + b*x))/(2*d*(b*c - a*d)^2) - Gamma(-1, a + b*x)/(2*d*(c + d*x)^2) + (b^2*ℯ^(-a + (b*c)/d)*Gamma(-1, (b*(c + d*x))/d))/(2*d*(b*c - a*d)^2) - (b^2*Gamma(0, a + b*x))/(b*c - a*d)^3 + (b^2*ℯ^(-a + (b*c)/d)*Gamma(0, (b*(c + d*x))/d))/(b*c - a*d)^3]
    @test_int [Gamma(-1, a + b*x)/(c+ d*x)^4, x, 8, (b^3*ℯ^(-a + (b*c)/d)*Gamma(-2, (b*(c + d*x))/d))/(3*d^2*(b*c - a*d)^2) + (b^3*Gamma(-1, a + b*x))/(3*d*(b*c - a*d)^3) - Gamma(-1, a + b*x)/(3*d*(c + d*x)^3) + (2*b^3*ℯ^(-a + (b*c)/d)*Gamma(-1, (b*(c + d*x))/d))/(3*d*(b*c - a*d)^3) - (b^3*Gamma(0, a + b*x))/(b*c - a*d)^4 + (b^3*ℯ^(-a + (b*c)/d)*Gamma(0, (b*(c + d*x))/d))/(b*c - a*d)^4]

    @test_int [(c+ d*x)^3*Gamma(-2, a + b*x), x, 8, -((d^2*(4*b*c - 3*a*d)*ℯ^(-a - b*x))/(4*b^4)) - ((b*c - a*d)^4*Gamma(-2, a + b*x))/(4*b^4*d) + ((c + d*x)^4*Gamma(-2, a + b*x))/(4*d) - ((b*c - a*d)^3*Gamma(-1, a + b*x))/b^4 - (3*d*(b*c - a*d)^2*Gamma(0, a + b*x))/(2*b^4) - (d^3*Gamma(2, b*x))/(ℯ^a*(4*b^4))]
    @test_int [(c+ d*x)^2*Gamma(-2, a + b*x), x, 7, -((d^2*ℯ^(-a - b*x))/(3*b^3)) - ((b*c - a*d)^3*Gamma(-2, a + b*x))/(3*b^3*d) + ((c + d*x)^3*Gamma(-2, a + b*x))/(3*d) - ((b*c - a*d)^2*Gamma(-1, a + b*x))/b^3 - (d*(b*c - a*d)*Gamma(0, a + b*x))/b^3]
    @test_int [(c+ d*x)^1*Gamma(-2, a + b*x), x, 6, -(((b*c - a*d)^2*Gamma(-2, a + b*x))/(2*b^2*d)) + ((c + d*x)^2*Gamma(-2, a + b*x))/(2*d) - ((b*c - a*d)*Gamma(-1, a + b*x))/b^2 - (d*Gamma(0, a + b*x))/(2*b^2)]
    @test_int [(c+ d*x)^0*Gamma(-2, a + b*x), x, 1, ((a + b*x)*Gamma(-2, a + b*x))/b - Gamma(-1, a + b*x)/b]
    @test_int [Gamma(-2, a + b*x)/(c+ d*x)^1, x, 0, Unintegrable(Gamma(-2, a + b*x)/(c + d*x), x)]
    @test_int [Gamma(-2, a + b*x)/(c+ d*x)^2, x, 7, (b*Gamma(-2, a + b*x))/(d*(b*c - a*d)) - Gamma(-2, a + b*x)/(d*(c + d*x)) - (b*Gamma(-1, a + b*x))/(b*c - a*d)^2 + (b*d*Gamma(0, a + b*x))/(b*c - a*d)^3 - (b*d*ℯ^(-a + (b*c)/d)*Gamma(0, (b*(c + d*x))/d))/(b*c - a*d)^3]
    @test_int [Gamma(-2, a + b*x)/(c+ d*x)^3, x, 8, (b^2*Gamma(-2, a + b*x))/(2*d*(b*c - a*d)^2) - Gamma(-2, a + b*x)/(2*d*(c + d*x)^2) - (b^2*Gamma(-1, a + b*x))/(b*c - a*d)^3 - (b^2*ℯ^(-a + (b*c)/d)*Gamma(-1, (b*(c + d*x))/d))/(2*(b*c - a*d)^3) + (3*b^2*d*Gamma(0, a + b*x))/(2*(b*c - a*d)^4) - (3*b^2*d*ℯ^(-a + (b*c)/d)*Gamma(0, (b*(c + d*x))/d))/(2*(b*c - a*d)^4)]
    @test_int [Gamma(-2, a + b*x)/(c+ d*x)^4, x, 9, (b^3*Gamma(-2, a + b*x))/(3*d*(b*c - a*d)^3) - Gamma(-2, a + b*x)/(3*d*(c + d*x)^3) - (b^3*ℯ^(-a + (b*c)/d)*Gamma(-2, (b*(c + d*x))/d))/(3*d*(b*c - a*d)^3) - (b^3*Gamma(-1, a + b*x))/(b*c - a*d)^4 - (b^3*ℯ^(-a + (b*c)/d)*Gamma(-1, (b*(c + d*x))/d))/(b*c - a*d)^4 + (2*b^3*d*Gamma(0, a + b*x))/(b*c - a*d)^5 - (2*b^3*d*ℯ^(-a + (b*c)/d)*Gamma(0, (b*(c + d*x))/d))/(b*c - a*d)^5]

    @test_int [(c+ d*x)^3*Gamma(-3, a + b*x), x, 8, -((d^3*ℯ^(-a - b*x))/(4*b^4)) - ((b*c - a*d)^4*Gamma(-3, a + b*x))/(4*b^4*d) + ((c + d*x)^4*Gamma(-3, a + b*x))/(4*d) - ((b*c - a*d)^3*Gamma(-2, a + b*x))/b^4 - (3*d*(b*c - a*d)^2*Gamma(-1, a + b*x))/(2*b^4) - (d^2*(b*c - a*d)*Gamma(0, a + b*x))/b^4]
    @test_int [(c+ d*x)^2*Gamma(-3, a + b*x), x, 7, -(((b*c - a*d)^3*Gamma(-3, a + b*x))/(3*b^3*d)) + ((c + d*x)^3*Gamma(-3, a + b*x))/(3*d) - ((b*c - a*d)^2*Gamma(-2, a + b*x))/b^3 - (d*(b*c - a*d)*Gamma(-1, a + b*x))/b^3 - (d^2*Gamma(0, a + b*x))/(3*b^3)]
    @test_int [(c+ d*x)^1*Gamma(-3, a + b*x), x, 6, -(((b*c - a*d)^2*Gamma(-3, a + b*x))/(2*b^2*d)) + ((c + d*x)^2*Gamma(-3, a + b*x))/(2*d) - ((b*c - a*d)*Gamma(-2, a + b*x))/b^2 - (d*Gamma(-1, a + b*x))/(2*b^2)]
    @test_int [(c+ d*x)^0*Gamma(-3, a + b*x), x, 1, ((a + b*x)*Gamma(-3, a + b*x))/b - Gamma(-2, a + b*x)/b]
    @test_int [Gamma(-3, a + b*x)/(c+ d*x)^1, x, 0, Unintegrable(Gamma(-3, a + b*x)/(c + d*x), x)]
    @test_int [Gamma(-3, a + b*x)/(c+ d*x)^2, x, 8, (b*Gamma(-3, a + b*x))/(d*(b*c - a*d)) - Gamma(-3, a + b*x)/(d*(c + d*x)) - (b*Gamma(-2, a + b*x))/(b*c - a*d)^2 + (b*d*Gamma(-1, a + b*x))/(b*c - a*d)^3 - (b*d^2*Gamma(0, a + b*x))/(b*c - a*d)^4 + (b*d^2*ℯ^(-a + (b*c)/d)*Gamma(0, (b*(c + d*x))/d))/(b*c - a*d)^4]
    @test_int [Gamma(-3, a + b*x)/(c+ d*x)^3, x, 9, (b^2*Gamma(-3, a + b*x))/(2*d*(b*c - a*d)^2) - Gamma(-3, a + b*x)/(2*d*(c + d*x)^2) - (b^2*Gamma(-2, a + b*x))/(b*c - a*d)^3 + (3*b^2*d*Gamma(-1, a + b*x))/(2*(b*c - a*d)^4) + (b^2*d*ℯ^(-a + (b*c)/d)*Gamma(-1, (b*(c + d*x))/d))/(2*(b*c - a*d)^4) - (2*b^2*d^2*Gamma(0, a + b*x))/(b*c - a*d)^5 + (2*b^2*d^2*ℯ^(-a + (b*c)/d)*Gamma(0, (b*(c + d*x))/d))/(b*c - a*d)^5]
    @test_int [Gamma(-3, a + b*x)/(c+ d*x)^4, x, 10, (b^3*Gamma(-3, a + b*x))/(3*d*(b*c - a*d)^3) - Gamma(-3, a + b*x)/(3*d*(c + d*x)^3) - (b^3*Gamma(-2, a + b*x))/(b*c - a*d)^4 + (b^3*ℯ^(-a + (b*c)/d)*Gamma(-2, (b*(c + d*x))/d))/(3*(b*c - a*d)^4) + (2*b^3*d*Gamma(-1, a + b*x))/(b*c - a*d)^5 + (4*b^3*d*ℯ^(-a + (b*c)/d)*Gamma(-1, (b*(c + d*x))/d))/(3*(b*c - a*d)^5) - (10*b^3*d^2*Gamma(0, a + b*x))/(3*(b*c - a*d)^6) + (10*b^3*d^2*ℯ^(-a + (b*c)/d)*Gamma(0, (b*(c + d*x))/d))/(3*(b*c - a*d)^6)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^(m/2)*Gamma(n, a+b*x)=#

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

    @test_int [x^(5/2)*Gamma(2, a + b*x), x, 5, (2/7)*x^(7/2)*Gamma(2, a + b*x) - (2*a*sqrt(x)*Gamma(9/2, b*x))/(ℯ^a*(7*b^3*sqrt(b*x))) - (2*sqrt(x)*Gamma(11/2, b*x))/(ℯ^a*(7*b^3*sqrt(b*x)))]
    @test_int [x^(3/2)*Gamma(2, a + b*x), x, 5, (2/5)*x^(5/2)*Gamma(2, a + b*x) - (2*a*sqrt(x)*Gamma(7/2, b*x))/(ℯ^a*(5*b^2*sqrt(b*x))) - (2*sqrt(x)*Gamma(9/2, b*x))/(ℯ^a*(5*b^2*sqrt(b*x)))]
    @test_int [x^(1/2)*Gamma(2, a + b*x), x, 5, (2/3)*x^(3/2)*Gamma(2, a + b*x) - (2*a*sqrt(x)*Gamma(5/2, b*x))/(ℯ^a*(3*b*sqrt(b*x))) - (2*sqrt(x)*Gamma(7/2, b*x))/(ℯ^a*(3*b*sqrt(b*x)))]
    @test_int [Gamma(2, a + b*x)/x^(1/2), x, 5, -((2*a*sqrt(x)*Gamma(3/2, b*x))/(ℯ^a*sqrt(b*x))) + 2*sqrt(x)*Gamma(2, a + b*x) - (2*sqrt(x)*Gamma(5/2, b*x))/(ℯ^a*sqrt(b*x))]
    @test_int [Gamma(2, a + b*x)/x^(3/2), x, 5, (2*a*sqrt(b*x)*Gamma(1/2, b*x))/(ℯ^a*sqrt(x)) + (2*b*sqrt(x)*Gamma(3/2, b*x))/(ℯ^a*sqrt(b*x)) - (2*Gamma(2, a + b*x))/sqrt(x)]
    @test_int [Gamma(2, a + b*x)/x^(5/2), x, 5, (2*a*b*sqrt(b*x)*Gamma(-(1/2), b*x))/(ℯ^a*(3*sqrt(x))) + (2*b*sqrt(b*x)*Gamma(1/2, b*x))/(ℯ^a*(3*sqrt(x))) - (2*Gamma(2, a + b*x))/(3*x^(3/2))]
    @test_int [Gamma(2, a + b*x)/x^(7/2), x, 5, (2*a*b^2*sqrt(b*x)*Gamma(-(3/2), b*x))/(ℯ^a*(5*sqrt(x))) + (2*b^2*sqrt(b*x)*Gamma(-(1/2), b*x))/(ℯ^a*(5*sqrt(x))) - (2*Gamma(2, a + b*x))/(5*x^(5/2))]
    @test_int [Gamma(2, a + b*x)/x^(9/2), x, 5, (2*a*b^3*sqrt(b*x)*Gamma(-(5/2), b*x))/(ℯ^a*(7*sqrt(x))) + (2*b^3*sqrt(b*x)*Gamma(-(3/2), b*x))/(ℯ^a*(7*sqrt(x))) - (2*Gamma(2, a + b*x))/(7*x^(7/2))]

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

    @test_int [x^(3/2)*Gamma(-2, a + b*x), x, 0, Unintegrable(x^(3/2)*Gamma(-2, a + b*x), x)]
    @test_int [x^(1/2)*Gamma(-2, a + b*x), x, 0, Unintegrable(sqrt(x)*Gamma(-2, a + b*x), x)]
    @test_int [Gamma(-2, a + b*x)/x^(1/2), x, 0, Unintegrable(Gamma(-2, a + b*x)/sqrt(x), x)]
    @test_int [Gamma(-2, a + b*x)/x^(3/2), x, 0, Unintegrable(Gamma(-2, a + b*x)/x^(3/2), x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^(m/3)*Gamma(n, a+b*x)=#

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

    @test_int [x^(4/3)*Gamma(2, a + b*x), x, 5, (3/7)*x^(7/3)*Gamma(2, a + b*x) - (3*a*x^(1/3)*Gamma(10/3, b*x))/(ℯ^a*(7*b^2*(b*x)^(1/3))) - (3*x^(1/3)*Gamma(13/3, b*x))/(ℯ^a*(7*b^2*(b*x)^(1/3)))]
    @test_int [x^(2/3)*Gamma(2, a + b*x), x, 5, (3/5)*x^(5/3)*Gamma(2, a + b*x) - (3*a*x^(2/3)*Gamma(8/3, b*x))/(ℯ^a*(5*b*(b*x)^(2/3))) - (3*x^(2/3)*Gamma(11/3, b*x))/(ℯ^a*(5*b*(b*x)^(2/3)))]
    @test_int [x^(1/3)*Gamma(2, a + b*x), x, 5, (3/4)*x^(4/3)*Gamma(2, a + b*x) - (3*a*x^(1/3)*Gamma(7/3, b*x))/(ℯ^a*(4*b*(b*x)^(1/3))) - (3*x^(1/3)*Gamma(10/3, b*x))/(ℯ^a*(4*b*(b*x)^(1/3)))]
    @test_int [Gamma(2, a + b*x)/x^(1/3), x, 5, -((3*a*x^(2/3)*Gamma(5/3, b*x))/(ℯ^a*(2*(b*x)^(2/3)))) + (3/2)*x^(2/3)*Gamma(2, a + b*x) - (3*x^(2/3)*Gamma(8/3, b*x))/(ℯ^a*(2*(b*x)^(2/3)))]
    @test_int [Gamma(2, a + b*x)/x^(2/3), x, 5, -((3*a*x^(1/3)*Gamma(4/3, b*x))/(ℯ^a*(b*x)^(1/3))) + 3*x^(1/3)*Gamma(2, a + b*x) - (3*x^(1/3)*Gamma(7/3, b*x))/(ℯ^a*(b*x)^(1/3))]
    @test_int [Gamma(2, a + b*x)/x^(4/3), x, 5, (3*a*(b*x)^(1/3)*Gamma(2/3, b*x))/(ℯ^a*x^(1/3)) + (3*b*x^(2/3)*Gamma(5/3, b*x))/(ℯ^a*(b*x)^(2/3)) - (3*Gamma(2, a + b*x))/x^(1/3)]

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

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*Gamma(n, a+b*x)*when*m*symbolic=#

    @test_int [(c + d*x)^m*Gamma(3, a + b*x), x, 6, ((c + d*x)^(1 + m)*Gamma(3, a + b*x))/(d*(1 + m)) - ((b*c - a*d)^2*ℯ^(-a + (b*c)/d)*(c + d*x)^m*Gamma(2 + m, (b*(c + d*x))/d))/(((b*(c + d*x))/d)^m*(b*d^2*(1 + m))) + (2*(b*c - a*d)*ℯ^(-a + (b*c)/d)*(c + d*x)^m*Gamma(3 + m, (b*(c + d*x))/d))/(((b*(c + d*x))/d)^m*(b*d*(1 + m))) - (ℯ^(-a + (b*c)/d)*(c + d*x)^m*Gamma(4 + m, (b*(c + d*x))/d))/(((b*(c + d*x))/d)^m*(b*(1 + m)))]
    @test_int [(c + d*x)^m*Gamma(2, a + b*x), x, 5, ((c + d*x)^(1 + m)*Gamma(2, a + b*x))/(d*(1 + m)) + ((b*c - a*d)*ℯ^(-a + (b*c)/d)*(c + d*x)^m*Gamma(2 + m, (b*(c + d*x))/d))/(((b*(c + d*x))/d)^m*(b*d*(1 + m))) - (ℯ^(-a + (b*c)/d)*(c + d*x)^m*Gamma(3 + m, (b*(c + d*x))/d))/(((b*(c + d*x))/d)^m*(b*(1 + m)))]
    @test_int [(c + d*x)^m*Gamma(1, a + b*x), x, 1, -((ℯ^(-a + (b*c)/d)*(c + d*x)^m*Gamma(1 + m, (b*(c + d*x))/d))/(((b*(c + d*x))/d)^m*b))]
    @test_int [(c + d*x)^m*Gamma(0, a + b*x), x, 0, Unintegrable((c + d*x)^m*Gamma(0, a + b*x), x)]
    @test_int [(c + d*x)^m*Gamma(-1, a + b*x), x, 0, Unintegrable((c + d*x)^m*Gamma(-1, a + b*x), x)]
    @test_int [(c + d*x)^m*Gamma(-2, a + b*x), x, 0, Unintegrable((c + d*x)^m*Gamma(-2, a + b*x), x)]
    @test_int [(c + d*x)^m*Gamma(-3, a + b*x), x, 0, Unintegrable((c + d*x)^m*Gamma(-3, a + b*x), x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*Gamma(n, a+b*x)*when*n*symbolic=#

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

    @test_int [(c + d*x)^4*Gamma(n, a + b*x), x, 9, -(((b*c - a*d)^5*Gamma(n, a + b*x))/(5*b^5*d)) + ((c + d*x)^5*Gamma(n, a + b*x))/(5*d) - ((b*c - a*d)^4*Gamma(1 + n, a + b*x))/b^5 - (2*d*(b*c - a*d)^3*Gamma(2 + n, a + b*x))/b^5 - (2*d^2*(b*c - a*d)^2*Gamma(3 + n, a + b*x))/b^5 - (d^3*(b*c - a*d)*Gamma(4 + n, a + b*x))/b^5 - (d^4*Gamma(5 + n, a + b*x))/(5*b^5)]
    @test_int [(c + d*x)^3*Gamma(n, a + b*x), x, 8, -(((b*c - a*d)^4*Gamma(n, a + b*x))/(4*b^4*d)) + ((c + d*x)^4*Gamma(n, a + b*x))/(4*d) - ((b*c - a*d)^3*Gamma(1 + n, a + b*x))/b^4 - (3*d*(b*c - a*d)^2*Gamma(2 + n, a + b*x))/(2*b^4) - (d^2*(b*c - a*d)*Gamma(3 + n, a + b*x))/b^4 - (d^3*Gamma(4 + n, a + b*x))/(4*b^4)]
    @test_int [(c + d*x)^2*Gamma(n, a + b*x), x, 7, -(((b*c - a*d)^3*Gamma(n, a + b*x))/(3*b^3*d)) + ((c + d*x)^3*Gamma(n, a + b*x))/(3*d) - ((b*c - a*d)^2*Gamma(1 + n, a + b*x))/b^3 - (d*(b*c - a*d)*Gamma(2 + n, a + b*x))/b^3 - (d^2*Gamma(3 + n, a + b*x))/(3*b^3)]
    @test_int [(c + d*x)^1*Gamma(n, a + b*x), x, 6, -(((b*c - a*d)^2*Gamma(n, a + b*x))/(2*b^2*d)) + ((c + d*x)^2*Gamma(n, a + b*x))/(2*d) - ((b*c - a*d)*Gamma(1 + n, a + b*x))/b^2 - (d*Gamma(2 + n, a + b*x))/(2*b^2)]
    @test_int [(c + d*x)^0*Gamma(n, a + b*x), x, 1, ((a + b*x)*Gamma(n, a + b*x))/b - Gamma(1 + n, a + b*x)/b]
    @test_int [Gamma(n, a + b*x)/(c + d*x)^1, x, 0, Unintegrable(Gamma(n, a + b*x)/(c + d*x), x)]
    @test_int [Gamma(n, a + b*x)/(c + d*x)^2, x, 0, Unintegrable(Gamma(n, a + b*x)/(c + d*x)^2, x)]
    @test_int [Gamma(n, a + b*x)/(c + d*x)^3, x, 0, Unintegrable(Gamma(n, a + b*x)/(c + d*x)^3, x)]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*Gamma(p, d*(a+b*log(c*x^n)))=#

    @test_int [x^2*Gamma(p, d*(a + b*log(c*x^n))), x, 4, (1/3)*x^3*Gamma(p, d*(a + b*log(c*x^n))) - ((1/3)*x^3*Gamma(p, -(((3 - b*d*n)*(a + b*log(c*x^n)))/(b*n)))*(d*(a + b*log(c*x^n)))^p)/(ℯ^((3*a)/(b*n))*(c*x^n)^(3/n)*(-(((3 - b*d*n)*(a + b*log(c*x^n)))/(b*n)))^p)]
    @test_int [x^1*Gamma(p, d*(a + b*log(c*x^n))), x, 4, (1/2)*x^2*Gamma(p, d*(a + b*log(c*x^n))) - ((1/2)*x^2*Gamma(p, -(((2 - b*d*n)*(a + b*log(c*x^n)))/(b*n)))*(d*(a + b*log(c*x^n)))^p)/(ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)*(-(((2 - b*d*n)*(a + b*log(c*x^n)))/(b*n)))^p)]
    @test_int [x^0*Gamma(p, d*(a + b*log(c*x^n))), x, 5, x*Gamma(p, d*(a + b*log(c*x^n))) - (x*Gamma(p, -(((1 - b*d*n)*(a + b*log(c*x^n)))/(b*n)))*(d*(a + b*log(c*x^n)))^p)/(ℯ^(a/(b*n))*(c*x^n)^n^(-1)*(-(((1 - b*d*n)*(a + b*log(c*x^n)))/(b*n)))^p)]
    @test_int [Gamma(p, d*(a + b*log(c*x^n)))/x^1, x, 3, -(Gamma(1 + p, a*d + b*d*log(c*x^n))/(b*d*n)) + (Gamma(p, a*d + b*d*log(c*x^n))*(a + b*log(c*x^n)))/(b*n)]
    @test_int [Gamma(p, d*(a + b*log(c*x^n)))/x^2, x, 4, -(Gamma(p, d*(a + b*log(c*x^n)))/x) + (ℯ^(a/(b*n))*(c*x^n)^(1/n)*Gamma(p, ((1 + b*d*n)*(a + b*log(c*x^n)))/(b*n))*(d*(a + b*log(c*x^n)))^p)/((((1 + b*d*n)*(a + b*log(c*x^n)))/(b*n))^p*x)]
    @test_int [Gamma(p, d*(a + b*log(c*x^n)))/x^3, x, 4, -(Gamma(p, d*(a + b*log(c*x^n)))/(2*x^2)) + (ℯ^((2*a)/(b*n))*(c*x^n)^(2/n)*Gamma(p, ((2 + b*d*n)*(a + b*log(c*x^n)))/(b*n))*(d*(a + b*log(c*x^n)))^p)/((((2 + b*d*n)*(a + b*log(c*x^n)))/(b*n))^p*(2*x^2))]

    @test_int [(e*x)^m*Gamma(p, d*(a + b*log(c*x^n))), x, 4, ((e*x)^(1 + m)*Gamma(p, d*(a + b*log(c*x^n))))/(e*(1 + m)) - ((e*x)^(1 + m)*Gamma(p, -(((1 + m - b*d*n)*(a + b*log(c*x^n)))/(b*n)))*(d*(a + b*log(c*x^n)))^p)/(ℯ^((a*(1 + m))/(b*n))*(c*x^n)^((1 + m)/n)*(-(((1 + m - b*d*n)*(a + b*log(c*x^n)))/(b*n)))^p*(e*(1 + m)))]

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

    @test_int [(c + d*x)^3*LogGamma(a + b*x), x, 4, -((6*d^3*PolyGamma(-5, a + b*x))/b^4) + (6*d^2*(c + d*x)*PolyGamma(-4, a + b*x))/b^3 - (3*d*(c + d*x)^2*PolyGamma(-3, a + b*x))/b^2 + ((c + d*x)^3*PolyGamma(-2, a + b*x))/b]
    @test_int [(c + d*x)^2*LogGamma(a + b*x), x, 3, (2*d^2*PolyGamma(-4, a + b*x))/b^3 - (2*d*(c + d*x)*PolyGamma(-3, a + b*x))/b^2 + ((c + d*x)^2*PolyGamma(-2, a + b*x))/b]
    @test_int [(c + d*x)^1*LogGamma(a + b*x), x, 2, -((d*PolyGamma(-3, a + b*x))/b^2) + ((c + d*x)*PolyGamma(-2, a + b*x))/b]
    @test_int [(c + d*x)^0*LogGamma(a + b*x), x, 1, PolyGamma(-2, a + b*x)/b]
    @test_int [LogGamma(a + b*x)/(c + d*x)^1, x, 0, Unintegrable(LogGamma(a + b*x)/(c + d*x), x)]
    @test_int [LogGamma(a + b*x)/(c + d*x)^2, x, 0, Unintegrable(LogGamma(a + b*x)/(c + d*x)^2, x)]

    @test_int [(c + d*x)^(3/2)*LogGamma(a + b*x), x, 0, Unintegrable((c + d*x)^(3/2)*LogGamma(a + b*x), x)]
    @test_int [(c + d*x)^(1/2)*LogGamma(a + b*x), x, 0, Unintegrable(sqrt(c + d*x)*LogGamma(a + b*x), x)]
    @test_int [LogGamma(a + b*x)/(c + d*x)^(1/2), x, 0, Unintegrable(LogGamma(a + b*x)/sqrt(c + d*x), x)]

    @test_int [(c + d*x)^2*log(Gamma(a + b*x)), x, 6, ((c + d*x)^3*log(Gamma(a + b*x)))/(3*d) - ((c + d*x)^3*LogGamma(a + b*x))/(3*d) + (2*d^2*PolyGamma(-4, a + b*x))/b^3 - (2*d*(c + d*x)*PolyGamma(-3, a + b*x))/b^2 + ((c + d*x)^2*PolyGamma(-2, a + b*x))/b]
    @test_int [(c + d*x)^1*log(Gamma(a + b*x)), x, 5, ((c + d*x)^2*log(Gamma(a + b*x)))/(2*d) - ((c + d*x)^2*LogGamma(a + b*x))/(2*d) - (d*PolyGamma(-3, a + b*x))/b^2 + ((c + d*x)*PolyGamma(-2, a + b*x))/b]
    @test_int [(c + d*x)^0*log(Gamma(a + b*x)), x, 4, x*log(Gamma(a + b*x)) - x*LogGamma(a + b*x) + PolyGamma(-2, a + b*x)/b]
    @test_int [log(Gamma(a + b*x))/(c + d*x)^1, x, 2, (log(c + d*x)*(log(Gamma(a + b*x)) - LogGamma(a + b*x)))/d + Unintegrable(LogGamma(a + b*x)/(c + d*x), x)]
    @test_int [log(Gamma(a + b*x))/(c + d*x)^2, x, 2, -(log(Gamma(a + b*x))/(d*(c + d*x))) + (b*Unintegrable(PolyGamma(0, a + b*x)/(c + d*x), x))/d]

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

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

    @test_int [(c + d*x)^3*PolyGamma(n, a + b*x), x, 4, -((6*d^3*PolyGamma(-4 + n, a + b*x))/b^4) + (6*d^2*(c + d*x)*PolyGamma(-3 + n, a + b*x))/b^3 - (3*d*(c + d*x)^2*PolyGamma(-2 + n, a + b*x))/b^2 + ((c + d*x)^3*PolyGamma(-1 + n, a + b*x))/b]
    @test_int [(c + d*x)^2*PolyGamma(n, a + b*x), x, 3, (2*d^2*PolyGamma(-3 + n, a + b*x))/b^3 - (2*d*(c + d*x)*PolyGamma(-2 + n, a + b*x))/b^2 + ((c + d*x)^2*PolyGamma(-1 + n, a + b*x))/b]
    @test_int [(c + d*x)^1*PolyGamma(n, a + b*x), x, 2, -((d*PolyGamma(-2 + n, a + b*x))/b^2) + ((c + d*x)*PolyGamma(-1 + n, a + b*x))/b]
    @test_int [(c + d*x)^0*PolyGamma(n, a + b*x), x, 1, PolyGamma(-1 + n, a + b*x)/b]
    @test_int [PolyGamma(n, a + b*x)/(c + d*x)^1, x, 0, Unintegrable(PolyGamma(n, a + b*x)/(c + d*x), x)]
    @test_int [PolyGamma(n, a + b*x)/(c + d*x)^2, x, 1, -(PolyGamma(n, a + b*x)/(d*(c + d*x))) + (b*Unintegrable(PolyGamma(1 + n, a + b*x)/(c + d*x), x))/d]
    @test_int [PolyGamma(n, a + b*x)/(c + d*x)^3, x, 2, -(PolyGamma(n, a + b*x)/(2*d*(c + d*x)^2)) - (b*PolyGamma(1 + n, a + b*x))/(2*d^2*(c + d*x)) + (b^2*Unintegrable(PolyGamma(2 + n, a + b*x)/(c + d*x), x))/(2*d^2)]

    @test_int [(c + d*x)^(3/2)*PolyGamma(n, a + b*x), x, 2, -((3*d*sqrt(c + d*x)*PolyGamma(-2 + n, a + b*x))/(2*b^2)) + ((c + d*x)^(3/2)*PolyGamma(-1 + n, a + b*x))/b + (3*d^2*Unintegrable(PolyGamma(-2 + n, a + b*x)/sqrt(c + d*x), x))/(4*b^2)]
    @test_int [(c + d*x)^(1/2)*PolyGamma(n, a + b*x), x, 1, (sqrt(c + d*x)*PolyGamma(-1 + n, a + b*x))/b - (d*Unintegrable(PolyGamma(-1 + n, a + b*x)/sqrt(c + d*x), x))/(2*b)]
    @test_int [PolyGamma(n, a + b*x)/(c + d*x)^(1/2), x, 0, Unintegrable(PolyGamma(n, a + b*x)/sqrt(c + d*x), x)]
    @test_int [PolyGamma(n, a + b*x)/(c + d*x)^(3/2), x, 1, -((2*PolyGamma(n, a + b*x))/(d*sqrt(c + d*x))) + (2*b*Unintegrable(PolyGamma(1 + n, a + b*x)/sqrt(c + d*x), x))/d]

    @test_int [x^2*PolyGamma(1, a + b*x), x, 3, -((2*x*LogGamma(a + b*x))/b^2) + (2*PolyGamma(-2, a + b*x))/b^3 + (x^2*PolyGamma(0, a + b*x))/b]

    @test_int [PolyGamma(1, a + b*x)/x^2 - (b*PolyGamma(2, a + b*x))/x, x, 2, -(PolyGamma(1, a + b*x)/x)]
    @test_int [PolyGamma(n, a + b*x)/x^2 - (b*PolyGamma(1 + n, a + b*x))/x, x, 2, -(PolyGamma(n, a + b*x)/x)]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*Gamma(c+dx)^m*PolyGamma(n, a+b*x)=#

    @test_int [Gamma(a + b*x)^n*PolyGamma(0, a + b*x), x, 1, Gamma(a + b*x)^n/(b*n)]
    @test_int [factorial(a + b*x)^n*PolyGamma(0, 1 + a + b*x), x, 1, factorial(a + b*x)^n/(b*n)]
    end

max_stars_repo_path: test/8 Special functions/8.6 Gamma functions.jl


text:
# This file is a part of EncodedArrays.jl, licensed under the MIT License (MIT).

"""
    abstract type AbstractArrayCodec <: Codecs.Codec end

Abstract type for arrays codecs.

Subtypes must implement the [`AbstractEncodedArray`](@ref) API.
Most coded should use [`EncodedArray`](@ref) as the concrete subtype of
`AbstractArrayCodec`. Codecs that use a custom subtype of
`AbstractEncodedArray` must implement

    EncodedArrays.encarraytype(::Type{<:AbstractArrayCodec},::Type{<:AbstractArray{T,N}})::Type{<:AbstractEncodedArray{T,N}}
"""
abstract type AbstractArrayCodec end
export AbstractArrayCodec

import Base.|>

"""
    ¦>(A::AbstractArray{T}, codec::AbstractArrayCodec)::AbstractEncodedArray

Encode `A` using `codec` and return an [`AbstractEncodedArray`](@ref). The
default implementation returns an [`EncodedArray`](@ref).
"""
function |>(A::AbstractArray{T}, codec::AbstractArrayCodec) where T
    encoded = Vector{UInt8}()
    encode_data!(encoded, codec, A)
    EncodedArray{T}(codec, size(A), encoded)
end

# Make AbstractArrayCodec behave as a Scalar for broadcasting
@inline Base.Broadcast.broadcastable(codec::AbstractArrayCodec) = (codec,)

"""
    encode_data!(encoded::AbstractVector{UInt8}, codec::AbstractArrayCodec, data::AbstractArray)

Will resize `encoded` as necessary to fit the encoded data.

Returns `encoded`.
"""
function encode_data! end

"""
    decode_data!(data::AbstractArray, codec::AbstractArrayCodec, encoded::AbstractVector{UInt8})

Depending on `codec`, may or may not resize `decoded` to fit the size of the
decoded data. Codecs may require `decoded` to be of correct size (e.g. to
improved performance or when the size/shape of the decoded data cannot be
easily inferred from the encoded data.

Returns `data`.
"""
function decode_data! end

"""
    AbstractEncodedArray{T,N} <: AbstractArray{T,N}

Abstract type for arrays that store their elements in encoded/compressed form.

In addition to the standard `AbstractArray` API, an `AbstractEncodedArray`
must support the functions

* `EncodedArrays.getcodec(A::EncodedArray)`: Returns the codec.
* `Base.codeunits(A::EncodedArray)`: Returns the internal encoded data
  representation.

Encoded arrays will typically be created via

    A_enc = (codec::AbstractArrayCodec)(A::AbstractArray)

or

    A_enc = AbstractEncodedArray(undef, codec::AbstractArrayCodec)
    append!(A_enc, B::AbstractArray)

Decoding happens via standard array conversion or assignment:

    A_dec = Array(A)
    A_dec = convert(Array,A)
    A_dec = A[:]

    A_dec = Array{T,N}(undef, size(A_enc)...)
    A_dec[:] = A_enc
"""
abstract type AbstractEncodedArray{T,N} <: AbstractArray{T,N} end
export AbstractEncodedArray

import Base.==
==(A::AbstractArray, B::AbstractEncodedArray) = A == Array(B)
==(A::AbstractEncodedArray, B::AbstractArray) = Array(A) == Array(B)
==(A::AbstractEncodedArray, B::AbstractEncodedArray) = Array(A) == Array(B)

"""
    EncodedArrays.getcodec(A::AbstractEncodedArray)::AbstractArrayCodec

Returns the codec used to encode/compress A.
"""
function getcodec end

"""
    EncodedArray{T,N,C,DV} <: AbstractEncodedArray{T,N}

Concrete type for [`AbstractEncodedArray`](@ref)s.

Constructor:

```julia
EncodedArray{T}(
    codec::AbstractArrayCodec,
    size::NTuple{N,Integer},
    encoded::AbstractVector{UInt8}
)
```

Codecs using `EncodedArray` only need to implement
[`EncodedArrays.encode_data!`](@ref) and [`EncodedArrays.decode_data!`](@ref).

If length of the decoded data can be inferred from the encoded data,
a constructor

    EncodedArray{T,N}(codec::MyCodec,encoded::AbstractVector{UInt8})

should also be defined. By default, two `EncodedArray`s that have the same
codec and size are assumed to be equal if and only if their code units are
equal.

Generic methods for the rest of the [`AbstractEncodedArray`](@ref) API are
already provided for `EncodedArray`.
"""
struct EncodedArray{T,N,C<:AbstractArrayCodec,DV<:AbstractVector{UInt8}} <: AbstractEncodedArray{T,N}
    codec::C
    size::NTuple{N,Int}
    encoded::DV
end
export EncodedArray

EncodedArray{T}(
    codec::AbstractArrayCodec,
    size::NTuple{N,Integer},
    encoded::AbstractVector{UInt8}
) where {T,N} = EncodedArray{T, N, typeof(codec),typeof(encoded)}(codec, size, encoded)

EncodedArray{T}(
    codec::AbstractArrayCodec,
    length::Integer,
    encoded::AbstractVector{UInt8}
) where {T} = EncodedArray{T, typeof(codec),typeof(encoded)}(codec, (len,), encoded)

EncodedArray{T,N,C,DV}(A::EncodedArray{T,N,C}) where {T,N,C,DV} = EncodedArray{T,N,C,DV}(A.codec, A.size, A.encoded)
Base.convert(::Type{EncodedArray{T,N,C,DV}}, A::EncodedArray{T,N,C}) where {T,N,C,DV} = EncodedArray{T,N,C,DV}(A)

@inline Base.size(A::EncodedArray) = A.size
@inline getcodec(A::EncodedArray) = A.codec
@inline Base.codeunits(A::EncodedArray) = A.encoded

# ToDo: Base.iscontiguous

function Base.Array{T,N}(A::EncodedArray{U,N}) where {T,N,U}
    B = Array{T,N}(undef, size(A)...)
    decode_data!(B, getcodec(A), codeunits(A))
end

Base.Array{T}(A::EncodedArray{U,N}) where {T,N,U} = Array{T,N}(A)
Base.Array(A::EncodedArray{T,N}) where {T,N} = Array{T,N}(A)

Base.Vector(A::EncodedArray{T,1}) where {T} = Array{T,1}(A)
Base.Matrix(A::EncodedArray{T,2}) where {T} = Array{T,2}(A)

Base.convert(::Type{Array{T,N}}, A::EncodedArray) where {T,N} = Array{T,N}(A)
Base.convert(::Type{Array{T}}, A::EncodedArray) where {T} = Array{T}(A)
Base.convert(::Type{Array}, A::EncodedArray) = Array(A)

Base.convert(::Type{Vector}, A::EncodedArray) = Vector(A)
Base.convert(::Type{Matrix}, A::EncodedArray) = Matrix(A)

Base.IndexStyle(A::EncodedArray) = IndexLinear()

function _getindex(A::EncodedArray, idxs::AbstractVector{Int})
    B = collect(A)
    if idxs == eachindex(IndexLinear(), A)
        B
    else
        B[idxs]
    end
end

_getindex(A::EncodedArray, i::Int) = collect(A)[i]

Base.@propagate_inbounds Base.getindex(A::EncodedArray, idxs) =
    _getindex(A, Base.to_indices(A, (idxs,))...)

@inline function _setindex!(A::AbstractArray, B::EncodedArray, idxs::AbstractVector{Int})
    @boundscheck let n = length(idxs), len_B = length(eachindex(B))
        n == len_B || Base.throw_setindex_mismatch(B, (n,))
    end

    if idxs == eachindex(A) || idxs == axes(A)
        decode_data!(A, getcodec(B), codeunits(B))
    else
        decode_data!(view(A, idxs), getcodec(B), codeunits(B))
    end

    A
end

Base.@propagate_inbounds function Base.setindex!(A::AbstractArray, B::EncodedArray, idxs::Colon)
    _setindex!(A, B, Base.to_indices(A, (idxs,))...)
end

@inline Base.@propagate_inbounds function Base.setindex!(A::Array, B::EncodedArray, idxs::AbstractVector{Int})
    @boundscheck checkbounds(A, idxs)
    _setindex!(A, B, Base.to_indices(A, (idxs,))...)
end

function _append!(A::AbstractVector, B::EncodedArray)
    n = length(eachindex(B))
    from = lastindex(A) + 1
    to = lastindex(A) + n
    resize!(A, to + 1 - firstindex(A))
    A[from:to] = B
    A
end

Base.append!(A::AbstractVector, B::EncodedArray) = _append!(A, B)
Base.append!(A::Vector, B::EncodedArray) = _append!(A, B)

# # ToDo (compatible with ElasticArrays.ElasticArray):
# Base.append!(A::AbstractArray{T,N}, B::EncodedArray) where {T,N} = ...

@inline function Base.copyto!(dest::AbstractArray, src::EncodedArray)
    @boundscheck if length(eachindex(dest)) < length(eachindex(src))
        throw(BoundsError())
    end
    decode_data!(dest, getcodec(src), codeunits(src))
end

# # ToDo:
# Base.copyto!(dest::AbstractArray, destoffs, src::EncodedArray, srcoffs, N) = ...

import Base.==
function ==(A::EncodedArray, B::EncodedArray)
    if getcodec(A) == getcodec(B) && size(A) == size(B)
        codeunits(A) == codeunits(B)
    else
        Array(A) == Array(B)
    end
end

"""
    const VectorOfEncodedArrays{T,N,...} = StructArray{EncodedArray{...},...}

Alias for vectors of encoded arrays that store the code units of all elements
in contiguous fashion using `StructArrays.StructArray`
and `ArraysOfArray.VectorOfArrays`.
"""
const VectorOfEncodedArrays{
    T,N,
    C<:AbstractArrayCodec,
    VC <: StructArray{C,1},
    VS <: AbstractVector{<:NTuple{N,Integer}},
    VOA <: VectorOfArrays
} = StructArray{
    EncodedArray{T,N,C,Array{UInt8,1}},
    1,
    NamedTuple{
        (:codec, :size, :encoded),
        Tuple{VC,VS,VOA}
    }
}

export VectorOfEncodedArrays

# Specialize getindex to properly support ArraysOfArrays, preventing
# conversion to exact element type:
@inline Base.getindex(A::StructArray{<:EncodedArray{T}}, I::Int...) where T =
    EncodedArray{T}(A.codec[I...], A.size[I...], A.encoded[I...])

const BroadcastedEncodeVectorOfArrays{T,N,C<:AbstractArrayCodec} = Base.Broadcast.Broadcasted{
    <:Base.Broadcast.AbstractArrayStyle{1},
    Tuple{Base.OneTo{Int64}},
    typeof(|>),
    <:Tuple{
        AbstractVector{<:AbstractArray{T,N}},
        Union{Tuple{C},AbstractVector{<:C}}
    }
}

@inline _get_1st_or_ith(A, i::Int) = (length(A) == 1) ? A[1] : A[i]

function _bcast_enc_impl(::Type{T}, ::Val{N}, ::Type{C}, data_arg, codec_arg) where {T,N,C}
    idxs_tuple = Base.Broadcast.combine_axes(data_arg, codec_arg)
    @assert length(idxs_tuple) == 1
    idxs = idxs_tuple[1]

    n = length(idxs)
    codec_vec = StructArray{C,1}(undef, n)
    size_vec = Vector{NTuple{N,Int}}(undef, n)
    encoded_vec = VectorOfVectors{UInt8}()

    sizehint!(encoded_vec.elem_ptr, n + 1)
    sizehint!(encoded_vec.kernel_size, n)

    for i in idxs
        data = _get_1st_or_ith(data_arg, i)
        codec = _get_1st_or_ith(codec_arg, i)

        codec_vec[i] = codec
        size_vec[i] = size(data)

        # ToDo: Improve, eliminate temporary memory allocation:
        tmp_encoded = encode_data!(Vector{UInt8}(), codec, data)
        push!(encoded_vec, tmp_encoded)
    end

    result = StructArray{EncodedArray{T,1,C,Vector{UInt8}}}((
        codec_vec,
        size_vec,
        encoded_vec
    ))

    @assert result isa VectorOfEncodedArrays{T,1,C}
    result
end

function Base.copy(instance::BroadcastedEncodeVectorOfArrays{T,N,C}) where {T,N,C}
    data_arg = instance.args[1]
    codec_arg = instance.args[2]
    _bcast_enc_impl(T, Val{N}(), C, data_arg, codec_arg)
end

const BroadcastedDecodeVectorOfArrays{T,N,C<:AbstractArrayCodec} = Base.Broadcast.Broadcasted{
    Base.Broadcast.DefaultArrayStyle{1},
    Tuple{Base.OneTo{Int64}},
    typeof(collect),
    <:Tuple{AbstractVector{<:EncodedArray{T,N,C}}}
}

function _bcast_dec_impl(::Type{T}, ::Val{N}, ::Type{C}, encoded_data) where {T,N,C}
    result = VectorOfArrays{T,N}()
    @inbounds for i in eachindex(encoded_data)
        x = encoded_data[i]
        push!(result, Fill(typemax(T), length(x)))
        copyto!(last(result), x)
    end
    result
end

function Base.copy(instance::BroadcastedDecodeVectorOfArrays{T,N,C}) where {T,N,C}
    _bcast_dec_impl(T, Val{N}(), C, instance.args[1])
end

# ToDo: SerialArrayCodec with decode_next, encode_next!, pos_type(codec),
#       finalize_codeunits!

# ToDo: Custom broadcasting over encoded array.


text:
srand(1000)

using SpikingNetworks;@unit
@codegen
parseParameters"""
area = 20000*um^2
Cm = 1*ufarad*cm^(-2)*area
gl = 5e-5*siemens*cm^(-2)*area
El = -65*mV
EK = -90*mV
ENa = 50*mV
g_na = 100*msiemens*cm^(-2)*area
g_kd = 30*msiemens*cm^(-2)*area
VT = -63*mV
taue = 5*ms
taui = 10*ms
Ee = 0*mV
Ei = -80*mV"""

parseParameters"""
N=4000;Ne=4*N÷5;Ni=N÷5;
we=spzeros(N,N);wi=spzeros(N,N)
we[1:N,1:Ne]=sprand(N,Ne,0.02)*6*nS
wi[1:N,Ne+1:N]=sprand(N,Ni,0.02)*67*nS"""

parseEquations"""
dv/dt=(gl*(El-v)+ge*(Ee-v)+gi*(Ei-v)-g_na*m^3*h*(v-ENa)-g_kd*n^4*(v-EK))/Cm
dm/dt=0.32*(13-v+VT)/(exp((13-v+VT)/4)-1)*(1-m)-0.28*(v-VT-40)/(exp((v-VT-40)/5)-1)*m
dn/dt=0.032*(15-v+VT)/(exp((15-v+VT)/5)-1)*(1-n)-0.5*exp((10-v+VT)/40)*n
dh/dt=0.128*exp((17-v+VT)/18)*(1-h)-4/(1+exp((40-v+VT)/5))*h
dge/dt=-ge/taue
dgi/dt=-gi/taui"""

net=init"""
N=N
v=El+5*(randn(N)-1)
m=0*ones(N)
n=0*ones(N)
h=1*ones(N)
ge=(randn(N)*1.5+4)*10*nS
gi=(randn(N)*12+20)*10*nS
we=we
wi=wi"""

spike"""
v.>-20
3"""

synapse"""
ge+=we
gi+=wi"""

stateMonitor"[1,10,100]"

net.dt=0.06
@time run!(net,1000/net.dt)
using MatlabPlot
mplot(net.vRecord')
# workspace()


text:
using Test
using Forecast

@testset "boxcox" begin

    x = rand(100)
    bc = boxcox(x)
    @test iboxcox(bc[:x],bc[:λ]) ≈ x

    x = rand(100) .- 0.5
    bc = boxcox(x)
    @test iboxcox(bc[:x],bc[:λ]) ≈ x

    x = rand(100)
    bc = boxcox(x)
    @test iboxcox(bc[:x],bc[:λ]) ≈ x

    x =  Int64.(round.(100*rand(100)))
    bc = boxcox(x)
    @test iboxcox(bc[:x],bc[:λ]) ≈ x

    x =  Int64.(round.(100*rand(100) .- 50))
    bc = boxcox(x)
    @test iboxcox(bc[:x],bc[:λ]) ≈ x

end


text:
### Indexing

## Helper functions
## locate runs, get index of run corresponding to the i'th value in the expanded runs
function ind2run(rle::RLEVector, i::Integer)
  re = rle.runends
  n = length(re)
  run = searchsortedfirst(re,i,1,n)
  run <= n || throw(BoundsError())  # Can't be < 1
  return(run)
end

function ind2run(rle::RLEVector,i::UnitRange)
  re = rle.runends
  n = length(re)
  left_run = searchsortedfirst(re,first(i),1,n)
  right_run = searchsortedfirst(re,last(i),left_run,n)
  right_run <= n || throw(BoundsError())  # Can't be < 1
  return( left_run:right_run )
end

function ind2run(rle::RLEVector, i::AbstractArray)
  re = rle.runends
  n = length(re)
  runs = searchsortedfirst(re,i,1,n)
  maximum(runs) <= n || throw(BoundsError())  # Can't be < 1
  return(runs)
end

# get index of the run corresponding to the i'th value in the expanded runs, index in run and remainder of run
#  (runindex, index_in_run, values_in_run_after_i)
function ind2runcontext(rle::RLEVector, i::Integer)
  run = ind2run(rle, i)
  runend = rle.runends[run]
  ind_in_run = run == 1 ? i : i - rle.runends[run-1]
  (run, ind_in_run, runend - i)
end

function ind2runcontext(rle::RLEVector, i::UnitRange)
  s = start(i)
  e = last(i)
  n = length(rle)
  runs = ind2run(rle, i)
  left_run = start(runs)
  right_run = last(runs)
  runend = rle.runends[right_run]
  ind_in_run = left_run == 1 ? s : s - rle.runends[left_run-1]
  (left_run, right_run, ind_in_run, runend - e)
end

function setrun!(rle::RLEVector, value, i::Int)
  run = ind2run(rle,i)
  rle.runvalues[run] = value
  return(rle)
end

## Just enough for AbstractArray
Base.linearindexing{T<:RLEVector}(::Type{T}) = Base.LinearFast()

endof(rle::RLEVector) = length(rle)

function Base.getindex(rle::RLEVector, i::Int)
    run = ind2run(rle,i)
    return( rle.runvalues[run] )
end

function Base.setindex!(rle::RLEVector, value, i::Int)
  run = ind2run(rle,i)
  runvalue = rle.runvalues[run]
  runend = rle.runends[run]
  value == runvalue && return rle # replace with same value, no-op
  previous_run = run - 1
  next_run = run + 1
  at_start_of_run = (previous_run > 0 && i == rle.runends[previous_run] + 1) || i == 1
  at_end_of_run = i == runend
  match_left = run > 1 && rle.runvalues[previous_run] == value
  match_right = run < nrun(rle) && rle.runvalues[next_run] == value
  if at_end_of_run
    if at_start_of_run # in a run of length 1
      if match_right && match_left
        splice!(rle.runvalues, previous_run:run)
        splice!(rle.runends, previous_run:run)
      elseif match_right
        splice!(rle.runvalues,run)
        splice!(rle.runends,run)
      elseif match_left
        splice!(rle.runvalues,run)
        splice!(rle.runends,previous_run)
      else
        rle.runvalues[run] = value
      end
    else # end of a run longer than 1
      if match_right
        rle.runends[run] = runend - 1
      else
        insert!(rle.runvalues, next_run, value)
        insert!(rle.runends, run, runend - 1)
      end
    end
  elseif at_start_of_run
    if match_left
      rle.runends[previous_run] = rle.runends[previous_run] + 1
    else
      insert!(rle.runvalues, run, value)
      insert!(rle.runends, run, i)
    end
 else # middle of a run, average case
    splice!(rle.runvalues, run, [runvalue,value,runvalue])
    splice!(rle.runends, run, [i-1,i,runend])
  end
  return(rle)
end

function Base.getindex(rle::RLEVector, ind::Array{Bool, 1})
    rle[ find(ind) ]
end

function Base.setindex!(rle::RLEVector, value::AbstractArray, ind::Array{Bool, 1})
    rle[ find(ind) ] = value
end

## Indexing optimizations
function Base.getindex(rle::RLEVector, ind::UnitRange)
    runs = ind2run(rle,collect(ind)) # OPTIMIZE ME: ind2runcontext
    rv = rle.runvalues[runs]
    return(rv)
end

#function Base.slice(rle::RLEVector, indices::UnitRange)
#    runs = ind2run(rle,indices)
#    rv = rle.runvalues[runs]
#    nrun = length(runs)
#    re = similar(rle.runends, nrun)
#    offset = start(indices) - 1
#    @inbounds for (i,r) in enumerate( start(runs):(last(runs) - 1) )
#        re[i] = rle.runends[r] - offset
#    end
#    re[nrun] = last(indices) - offset # length(indices)
#    rval = RLEVector{eltype(rv),eltype(re)}(rv,re)
#    return(rval)
#end

function Base.getindex(rle::RLEVector, ind::AbstractVector)
    run_indices = ind2run(rle, ind)
    return( RLEVector( rle.runvalues[ run_indices ] ) )
end

function Base.setindex!{T1,T2}(rle::RLEVector{T1,T2}, value::T1, indices::UnitRange)
  runs = ind2run(rle,indices)
  left_run = first(runs)
  right_run = last(runs)
  left_runvalue = rle.runvalues[left_run]
  right_runvalue = rle.runvalues[right_run]
  left_runend = rle.runends[left_run]
  right_runend = rle.runends[right_run]
  left_i = start(indices)
  right_i = last(indices)
  previous_run = left_run - 1
  next_run = right_run + 1
  at_start_of_run = (previous_run > 0 && left_i == rle.runends[previous_run] + 1) || left_i == 1
  at_end_of_run = right_i == right_runend
  match_left = left_run > 1 && rle.runvalues[previous_run] == value
  match_right = right_run < nrun(rle) && rle.runvalues[next_run] == value
  adjusted_runvalues = Array(typeof(left_runvalue),0)
  adjusted_runends = Array(typeof(left_runend),0)
  if at_end_of_run
    if at_start_of_run # in a run of length 1
      if match_right && match_left
        left_run = previous_run
      elseif match_right
        # do nothing
      elseif match_left
        rle.runends[previous_run] = right_runend
      else
        rle.runvalues[right_run] = value
        right_run = right_run - 1
      end
    else
      if match_right
        rle.runends[left_run] = left_i - 1
        left_run = left_run + 1
      else
        adjusted_runvalues = [left_runvalue,value]
        adjusted_runends = [left_i-1,right_i]
      end
    end
  elseif at_start_of_run
    if match_left
      rle.runends[previous_run] = last(indices)
      right_run = right_run - 1
    else
      adjusted_runvalues = [value,right_runvalue]
      adjusted_runends = [right_i,right_runend]
    end
  else # middle of a run, average case
    adjusted_runvalues = [left_runvalue,value,right_runvalue]
    adjusted_runends = [left_i-1,right_i,right_runend]
  end
  adjusted_runs = left_run:right_run
  splice!(rle.runvalues,adjusted_runs,adjusted_runvalues)
  splice!(rle.runends,adjusted_runs,adjusted_runends)
  return(rle)
end

## Getter shortcuts
function head(x::RLEVector,l::Integer=6)
    collect(x[ 1:l ])
end

function tail(x::RLEVector,l::Integer=6)
    collect( x[ length(x)-(l-1):end ] )
end

## Iterators
function start(rle::RLEVector)
  (1,1)
end

function next(rle::RLEVector, state)
  if state[2] == rle.runends[ state[1] ]
    newstate = (state[1] + 1, state[2] + 1)
  else
    newstate = (state[1],state[2] + 1)
  end
  return( (rle.runvalues[state[1]], newstate) )
end

function done(rle::RLEVector, state)
  state[1] > nrun(rle)
end

# Iterator for ranges based on RLE e.g. (value, start, end)
immutable RLEEachIterator{T1,T2}
    rle::RLEVector{T1,T2}
end
each(x::RLEVector) = RLEEachIterator(x)

function start(x::RLEEachIterator)
    1
end

function next(x::RLEEachIterator, state)
    newstate = state + 1
    first = rfirst(x.rle,state)
    last = rlast(x.rle)[state]
    ( (rvalue(x.rle)[state], first:last ), newstate )
end

function done(x::RLEEachIterator, state)
    state > nrun(x.rle)
end


text:
# This file is a part of EmpiricalDistributions.jl, licensed under the MIT License (MIT).

__precompile__(true)

module EmpiricalDistributions

using LinearAlgebra
using Random
using Statistics

using Adapt
using ArraysOfArrays
using Distributions
using StatsBase

include("utils.jl")
include("uv_binned_dist.jl")
include("mv_binned_dist.jl")

end # module


text:
module AWS4AuthRequest

using ..Base64
using ..Dates
using MbedTLS: digest, MD_SHA256, MD_MD5
import ..Layer, ..request, ..Headers
using ..URIs
using ..Pairs: getkv, setkv, rmkv
import ..@debug, ..DEBUG_LEVEL

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

Add a [AWS Signature Version 4](http://docs.aws.amazon.com/general/latest/gr/signature-version-4.html)
`Authorization` header to a `Request`.

Credentials are read from environment variables `AWS_ACCESS_KEY_ID`,
`AWS_SECRET_ACCESS_KEY` and `AWS_SESSION_TOKEN`.
"""
abstract type AWS4AuthLayer{Next <: Layer} <: Layer end
export AWS4AuthLayer

function request(::Type{AWS4AuthLayer{Next}},
                 url::URI, req, body; kw...) where Next

    if !haskey(kw, :aws_access_key_id) &&
       !haskey(ENV, "AWS_ACCESS_KEY_ID")
        kw = merge(dot_aws_credentials(), kw)
    end

    sign_aws4!(req.method, url, req.headers, req.body; kw...)

    return request(Next, url, req, body; kw...)
end

# Normalize whitespace to the form required in the canonical headers.
# Note that the expected format for multiline headers seems not to be explicitly
# documented, but Amazon provides a test case for it, so we'll match that behavior.
# We replace each `\n` with a `,` and remove all whitespace around the newlines,
# then any remaining contiguous whitespace is replaced with a single space.
function _normalize_ws(s::AbstractString)
    if any(isequal('\n'), s)
        join(map(_normalize_ws, split(s, '\n')), ',')
    else
        replace(strip(s), r"\s+" => " ")
    end
end

function sign_aws4!(method::String,
                    url::URI,
                    headers::Headers,
                    body::Vector{UInt8};
                    body_sha256::Vector{UInt8}=digest(MD_SHA256, body),
                    body_md5::Vector{UInt8}=digest(MD_MD5, body),
                    t::Union{DateTime,Nothing}=nothing,
                    timestamp::DateTime=now(Dates.UTC),
                    aws_service::String=String(split(url.host, ".")[1]),
                    aws_region::String=String(split(url.host, ".")[2]),
                    aws_access_key_id::String=ENV["AWS_ACCESS_KEY_ID"],
                    aws_secret_access_key::String=ENV["AWS_SECRET_ACCESS_KEY"],
                    aws_session_token::String=get(ENV, "AWS_SESSION_TOKEN", ""),
                    token_in_signature=true,
                    include_md5=true,
                    include_sha256=true,
                    kw...)
    if t !== nothing
        Base.depwarn("The `t` keyword argument to `sign_aws4!` is deprecated; use " *
                     "`timestamp` instead.", :sign_aws4!)
        timestamp = t
    end

    # ISO8601 date/time strings for time of request...
    date = Dates.format(timestamp, dateformat"yyyymmdd")
    datetime = Dates.format(timestamp, dateformat"yyyymmddTHHMMSS\Z")

    # Authentication scope...
    scope = [date, aws_region, aws_service, "aws4_request"]

    # Signing key generated from today's scope string...
    signing_key = string("AWS4", aws_secret_access_key)
    for element in scope
        signing_key = digest(MD_SHA256, element, signing_key)
    end

    # Authentication scope string...
    scope = join(scope, "/")

    # SHA256 hash of content...
    content_hash = bytes2hex(body_sha256)

    # HTTP headers...
    rmkv(headers, "Authorization")
    setkv(headers, "host", url.host)
    setkv(headers, "x-amz-date", datetime)
    include_md5 && setkv(headers, "Content-MD5", base64encode(body_md5))
    if (aws_service == "s3" && method == "PUT") || include_sha256
        # This header is required for S3 PUT requests. See the documentation at
        # https://docs.aws.amazon.com/AmazonS3/latest/API/sig-v4-header-based-auth.html
        setkv(headers, "x-amz-content-sha256", content_hash)
    end
    if aws_session_token != ""
        setkv(headers, "x-amz-security-token", aws_session_token)
    end

    # Sort and lowercase() Headers to produce canonical form...
    unique_header_keys = Vector{String}()
    normalized_headers = Dict{String,Vector{String}}()
    for (k, v) in sort!([lowercase(k) => v for (k, v) in headers], by=first)
        # Some services want the token included as part of the signature
        if k == "x-amz-security-token" && !token_in_signature
            continue
        end
        # In Amazon's examples, they exclude Content-Length from signing. This does not
        # appear to be addressed in the documentation, so we'll just mimic the example.
        if k == "content-length"
            continue
        end
        if !haskey(normalized_headers, k)
            normalized_headers[k] = Vector{String}()
            push!(unique_header_keys, k)
        end
        push!(normalized_headers[k], _normalize_ws(v))
    end
    canonical_headers = map(unique_header_keys) do k
        string(k, ':', join(normalized_headers[k], ','))
    end
    signed_headers = join(unique_header_keys, ';')

    # Sort Query String...
    query = sort!(collect(queryparams(url.query)), by=first)

    # Paths for requests to S3 should be escaped but not normalized. See
    # http://docs.aws.amazon.com/AmazonS3/latest/API/sig-v4-header-based-auth.html#canonical-request
    # Note that escapepath escapes ~ per RFC 1738, but Amazon includes an example in their
    # signature v4 test suite where ~ remains unescaped. We follow the spec here and thus
    # deviate from Amazon's example in this case.
    path = escapepath(aws_service == "s3" ? url.path : URIs.normpath(url.path))

    # Create hash of canonical request...
    canonical_form = join([method,
                           path,
                           escapeuri(query),
                           join(canonical_headers, "\n"),
                           "",
                           signed_headers,
                           content_hash], "\n")
    @debug 3 "AWS4 canonical_form: $canonical_form"

    canonical_hash = bytes2hex(digest(MD_SHA256, canonical_form))

    # Create and sign "String to Sign"...
    string_to_sign = "AWS4-HMAC-SHA256\n$datetime\n$scope\n$canonical_hash"
    signature = bytes2hex(digest(MD_SHA256, string_to_sign, signing_key))

    @debug 3 "AWS4 string_to_sign: $string_to_sign"
    @debug 3 "AWS4 signature: $signature"

    # Append Authorization header...
    setkv(headers, "Authorization", string(
        "AWS4-HMAC-SHA256 ",
        "Credential=$aws_access_key_id/$scope, ",
        "SignedHeaders=$signed_headers, ",
        "Signature=$signature"
    ))
end

using IniFile

credentials = NamedTuple()

"""
Load Credentials from [AWS CLI ~/.aws/credentials file]
(http://docs.aws.amazon.com/cli/latest/userguide/cli-config-files.html).
"""
function dot_aws_credentials()::NamedTuple

    global credentials
    if !isempty(credentials)
        return credentials
    end

    f = get(ENV, "AWS_CONFIG_FILE", joinpath(homedir(), ".aws", "credentials"))
    p = get(ENV, "AWS_DEFAULT_PROFILE", get(ENV, "AWS_PROFILE", "default"))

    if !isfile(f)
        return NamedTuple()
    end

    ini = read(Inifile(), f)

    credentials = (
        aws_access_key_id = String(get(ini, p, "aws_access_key_id")),
        aws_secret_access_key = String(get(ini, p, "aws_secret_access_key")))
end

end # module AWS4AuthRequest


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

struct DecisionIndex <: MOI.AbstractConstraintAttribute end
MOI.is_set_by_optimize(::DecisionIndex) = true

struct DecisionCoefficientChange{T} <: MOI.AbstractFunctionModification
    decision::MOI.VariableIndex
    new_coefficient::T
end

struct DecisionMultirowChange{T} <: MOI.AbstractFunctionModification
    decision::MOI.VariableIndex
    new_coefficients::Vector{Tuple{Int64, T}}
end

struct DecisionStateChange <: MOI.AbstractFunctionModification
    index::Int
    new_state::DecisionState
end

struct KnownValuesChange <: MOI.AbstractFunctionModification end

const DecisionModification = Union{DecisionCoefficientChange, DecisionStateChange, KnownValuesChange}
const VectorDecisionModification = Union{DecisionMultirowChange, DecisionStateChange, KnownValuesChange}


text:
import Mustache

template = """
#=
This file is auto-generated. Do not edit.
=#
\"\"\"
    mutable struct {{struct_name}}{{#parametric}}{T <: {{parametric}}}{{/parametric}} <: {{supertype}}
        {{#parameters}}
        {{name}}::{{{data_type}}}
        {{/parameters}}
    end

{{#docstring}}{{{docstring}}}{{/docstring}}

# Arguments
{{#parameters}}
- `{{name}}::{{{data_type}}}`{{#comment}}: {{{comment}}}{{/comment}}{{#valid_range}}, validation range: `{{{valid_range}}}`{{/valid_range}}{{#validation_action}}, action if invalid: `{{{validation_action}}}`{{/validation_action}}
{{/parameters}}
\"\"\"
mutable struct {{struct_name}}{{#parametric}}{T <: {{parametric}}}{{/parametric}} <: {{supertype}}
    {{#parameters}}
    {{#comment}}"{{{comment}}}"\n    {{/comment}}{{name}}::{{{data_type}}}
    {{/parameters}}
    {{#inner_constructor_check}}

    function {{struct_name}}({{#parameters}}{{name}}, {{/parameters}})
        ({{#parameters}}{{name}}, {{/parameters}}) = {{inner_constructor_check}}(
            {{#parameters}}
            {{name}},
            {{/parameters}}
        )
        new({{#parameters}}{{name}}, {{/parameters}})
    end
    {{/inner_constructor_check}}
end

{{#needs_positional_constructor}}
function {{constructor_func}}({{#parameters}}{{^internal_default}}{{name}}{{#default}}={{default}}{{/default}}, {{/internal_default}}{{/parameters}}){{{closing_constructor_text}}}
    {{constructor_func}}({{#parameters}}{{^internal_default}}{{name}}, {{/internal_default}}{{/parameters}}{{#parameters}}{{#internal_default}}{{{internal_default}}}, {{/internal_default}}{{/parameters}})
end
{{/needs_positional_constructor}}

function {{constructor_func}}(; {{#parameters}}{{name}}{{#kwarg_value}}{{{kwarg_value}}}{{/kwarg_value}}, {{/parameters}}){{{closing_constructor_text}}}
    {{constructor_func}}({{#parameters}}{{name}}, {{/parameters}})
end

{{#has_null_values}}
# Constructor for demo purposes; non-functional.
function {{constructor_func}}(::Nothing){{{closing_constructor_text}}}
    {{constructor_func}}(;
        {{#parameters}}
        {{^internal_default}}
        {{name}}={{#quotes}}"{{null_value}}"{{/quotes}}{{^quotes}}{{null_value}}{{/quotes}},
        {{/internal_default}}
        {{/parameters}}
    )
end

{{/has_null_values}}
{{#accessors}}
{{#create_docstring}}\"\"\"Get [`{{struct_name}}`](@ref) `{{name}}`.\"\"\"{{/create_docstring}}
{{accessor}}(value::{{struct_name}}) = {{#needs_conversion}}get_value(value, value.{{name}}){{/needs_conversion}}{{^needs_conversion}}value.{{name}}{{/needs_conversion}}
{{/accessors}}

{{#setters}}
{{#create_docstring}}\"\"\"Set [`{{struct_name}}`](@ref) `{{name}}`.\"\"\"{{/create_docstring}}
{{setter}}(value::{{struct_name}}, val) = value.{{name}} = {{#needs_conversion}}set_value(value, val){{/needs_conversion}}{{^needs_conversion}}val{{/needs_conversion}}
{{/setters}}

{{#custom_code}}
{{{custom_code}}}
{{/custom_code}}
"""

function read_json_data(filename::String)
    return open(filename) do io
        data = JSON3.read(io, Dict)
        if data isa Array
            return data
        elseif data isa Dict && haskey(data, "auto_generated_structs")
            return data["auto_generated_structs"]
        else
            throw(DataFormatError("{filename} has invalid format"))
        end
    end
end

function generate_structs(directory, data::Vector; print_results=true)
    struct_names = Vector{String}()
    unique_accessor_functions = Set{String}()
    unique_setter_functions = Set{String}()

    for item in data
        has_internal = false
        accessors = Vector{Dict}()
        setters = Vector{Dict}()
        item["has_null_values"] = true
        has_non_default_values = false

        item["constructor_func"] = item["struct_name"]
        item["closing_constructor_text"] = ""
        if haskey(item, "parametric")
            item["constructor_func"] *= "{T}"
            item["closing_constructor_text"] = " where T <: $(item["parametric"])"
        end

        parameters = Vector{Dict}()
        for field in item["fields"]
            param = field
            param["struct_name"] = item["struct_name"]
            if haskey(param, "valid_range")
                if typeof(param["valid_range"]) == Dict{String, Any}
                    min = param["valid_range"]["min"]
                    max = param["valid_range"]["max"]
                    param["valid_range"] = "($min, $max)"
                elseif typeof(param["valid_range"]) == String
                    param["valid_range"] = param["valid_range"]
                end
            end
            push!(parameters, param)

            # Allow accessor functions to be re-implemented from another module.
            # If this key is defined then the accessor function will not be exported.
            # Example:  get_name is defined in InfrastructureSystems and re-implemented in
            # PowerSystems.
            if haskey(param, "accessor_module")
                accessor_module = param["accessor_module"] * "."
                create_docstring = false
            else
                accessor_module = ""
                create_docstring = true
            end
            accessor_name = accessor_module * "get_" * param["name"]
            setter_name = accessor_module * "set_" * param["name"] * "!"
            push!(
                accessors,
                Dict(
                    "name" => param["name"],
                    "accessor" => accessor_name,
                    "create_docstring" => create_docstring,
                    "needs_conversion" => get(param, "needs_conversion", false),
                ),
            )
            include_setter = !get(param, "exclude_setter", false)
            if include_setter
                push!(
                    setters,
                    Dict(
                        "name" => param["name"],
                        "setter" => setter_name,
                        "data_type" => param["data_type"],
                        "create_docstring" => create_docstring,
                        "needs_conversion" => get(param, "needs_conversion", false),
                    ),
                )
            end
            if field["name"] != "internal" && accessor_module == ""
                push!(unique_accessor_functions, accessor_name)
                push!(unique_setter_functions, setter_name)
            end

            param["kwarg_value"] = ""
            if haskey(param, "default")
                param["kwarg_value"] = "=" * param["default"]
            elseif haskey(param, "internal_default")
                param["kwarg_value"] = "=" * string(param["internal_default"])
                has_internal = true
                continue
            else
                has_non_default_values = true
            end

            # This controls whether a demo constructor will be generated.
            if !haskey(param, "null_value") && !haskey(param, "default")
                item["has_null_values"] = false
            else
                if !haskey(param, "null_value")
                    item["null_value"] = param["default"]
                end
                if param["data_type"] == "String"
                    param["quotes"] = true
                end
            end
        end

        item["parameters"] = parameters
        item["accessors"] = accessors
        item["setters"] = setters
        # If all parameters have defaults then the positional constructor will
        # collide with the kwarg constructor.
        item["needs_positional_constructor"] = has_internal && has_non_default_values

        filename = joinpath(directory, item["struct_name"] * ".jl")
        open(filename, "w") do io
            write(io, Mustache.render(template, item))
            push!(struct_names, item["struct_name"])
        end

        if print_results
            println("Wrote $filename")
        end
    end

    accessors = sort!(collect(unique_accessor_functions))
    setters = sort!(collect(unique_setter_functions))
    filename = joinpath(directory, "includes.jl")
    open(filename, "w") do io
        for name in struct_names
            write(io, "include(\"$name.jl\")\n")
        end
        write(io, "\n")

        for accessor in accessors
            write(io, "export $accessor\n")
        end
        for setter in setters
            write(io, "export $setter\n")
        end
        if print_results
            println("Wrote $filename")
        end
    end
end

function generate_structs(
    input_file::AbstractString,
    output_directory::AbstractString;
    print_results=true,
)
    # Include each generated file.
    if !isdir(output_directory)
        mkdir(output_directory)
    end

    data = read_json_data(input_file)
    generate_structs(output_directory, data, print_results=print_results)
    return
end

"""
Return true if the structs defined in existing_dir match structs freshly-generated
from descriptor_file.
"""
function test_generated_structs(descriptor_file, existing_dir)
    output_dir = "tmp-test-generated-structs"
    if isdir(output_dir)
        rm(output_dir; recursive=true)
    end
    mkdir(output_dir)

    generate_structs(descriptor_file, output_dir; print_results=false)

    matched = true
    for (file1, file2) in zip(readdir(output_dir), readdir(existing_dir))
        path1 = joinpath(output_dir, file1)
        path2 = joinpath(existing_dir, file2)
        for (line1, line2) in zip(readlines(path1), readlines(path2))
            # Note: must strip the line endings.
            line1 = strip(line1)
            line2 = strip(line2)
            if line1 != line2
                @error "Generated structs do not match descriptor file" file1 line1 line2
                matched = false
                # Every line will now fail. Trying to use system utilities like diff didn't
                # work well across platforms.
                break
            end
        end
    end

    rm(output_dir; recursive=true)
    return matched
end