text:
mutable struct ILDA
    K::Int          # topics
    D::Int          # documents
    I::Int          # features
    J::Vector{Int}  # values feature

    η::Vector{Float64}              # topic hyperparameter
    λ::Vector{Matrix{Float64}}      # topic variational parameters
    β::Vector{Matrix{Float64}}      # topic-feature distributions
    Elnβ::Vector{Matrix{Float64}}   # expectation of ln(β)

    α::Float64                  # doc-topic hyperparameter
    γ::Matrix{Float64}          # doc-topic variational parameters
    θ::Matrix{Float64}          # doc-topic distributions
    Elnθ::Matrix{Float64}       # expectation of ln(θ)
    ϕ::Vector{Matrix{Float64}}  # Z variational parameters

    features::Matrix{Int}
    X::Vector{Matrix{Int}}

    converged::Bool
    elbo::Float64
    ll::Float64

    function ILDA(k::Int, α::Float64, η::Vector{Float64},
            features::Matrix{Int}, X::Vector{Matrix{Int}})
        model = new()

        model.K = k
        model.α = α
        model.η = copy(η)
        model.X = X
        model.D = length(X)
        model.I = size(features, 2)
        model.J = vec(maximum(features, dims=1))
        model.features = features

        model.λ = [rand(1:100, model.J[i], model.K) for i in 1:model.I]
        model.Elnβ = [
            Array{Float64}(undef, model.J[i], model.K) for i in 1:model.I
        ]
        update_Elnβ!(model)

        model.γ = fill(1.0, model.K, model.D)
        model.Elnθ = Array{Float64}(undef, model.K, model.D)
        update_Elnθ!(model)

        model.ϕ = [
            fill(1.0 / model.K, model.K, size(model.X[d], 1))
            for d in 1:model.D
        ]

        model.converged = false

        return model
    end
end

function ILDA(k::Int, α::Float64, η::Float64, features::Matrix{Int},
        X::Vector{Matrix{Int}})
    I = size(features, 2)
    return ILDA(k, α, fill(η, I), features, X)
end

function update_ϕ!(model::ILDA)
    for d in 1:model.D
        for w in 1:size(model.X[d], 1)
            model.ϕ[d][:, w] .= model.Elnθ[:, d]

            v = model.X[d][w, 1]
            for i in 1:model.I
                j = model.features[v, i]
                model.ϕ[d][:, w] .+= model.Elnβ[i][j, :]
            end
        end

        model.ϕ[d] .= exp.(model.ϕ[d]) ./ sum(exp.(model.ϕ[d]), dims=1)
    end
end

function update_Elnθ!(model::ILDA)
    model.Elnθ .= digamma.(model.γ) .- digamma.(sum(model.γ, dims=1))
end

function update_γ!(model::ILDA)
    model.γ .= model.α

    for d in 1:model.D
        model.γ[:, d] .+= model.ϕ[d] * model.X[d][:, 2]
    end

    update_Elnθ!(model)
end

function update_θ!(model::ILDA)
    model.θ = model.γ ./ sum(model.γ, dims=1)
end

function update_Elnβ!(model::ILDA)
    for i in 1:model.I
        model.Elnβ[i] .= (
            digamma.(model.λ[i]) .- digamma.(sum(model.λ[i], dims=1))
        )
    end
end

function update_λ!(model::ILDA)
    for i in 1:model.I
        model.λ[i] .= model.η[i]
    end

    for d in 1:model.D
        Nϕ = model.ϕ[d]' .* model.X[d][:, 2]

        for w in 1:size(model.X[d], 1)
            v = model.X[d][w, 1]

            for i in 1:model.I
                j = model.features[v, i]
                model.λ[i][j, :] .+= Nϕ[w, :]
            end
        end
    end

    update_Elnβ!(model)
end

function update_β!(model::ILDA)
    model.β = [model.λ[i] ./ sum(model.λ[i], dims=1) for i in 1:model.I]
end

function calculate_ElnPβ(model::ILDA)
    lnp = 0.0
    for i in 1:model.I
        lnp += model.K * (
            lgamma(model.J[i] * model.η[i]) - model.J[i] * lgamma(model.η[i])
        )
        lnp += (model.η[i] - 1) * sum(model.Elnβ[i])
    end
    return lnp
end

function calculate_ElnPθ(model::ILDA)
    lnp = model.D * (lgamma(model.K * model.α) - model.K * lgamma(model.α))
    lnp += (model.α - 1) * sum(model.Elnθ)
    return lnp
end

function calculate_ElnPZ(model::ILDA)
    lnp = 0.0
    for d in 1:model.D
        lnp += sum(model.ϕ[d] .* model.Elnθ[:, d] .* model.X[d][:, 2]')
    end
    return lnp
end

function calculate_ElnPX(model::ILDA)
    lnp = 0.0
    for d in 1:model.D
        Nϕ = model.ϕ[d]' .* model.X[d][:, 2]

        for w in 1:size(model.X[d], 1)
            v = model.X[d][w, 1]

            for i in 1:model.I
                j = model.features[v, i]
                lnp += sum(Nϕ[w, :] .* model.Elnβ[i][j, :])
            end
        end
    end
    return lnp
end

function calculate_ElnQβ(model::ILDA)
    lnq = 0.0
    for i in 1:model.I
        lnq = sum(lgamma.(model.λ[i])) - sum(lgamma.(sum(model.λ[i], dims=1)))
        lnq -= sum((model.λ[i] .- 1) .* model.Elnβ[i])
    end
    return lnq
end

function calculate_ElnQθ(model::ILDA)
    lnq = sum(lgamma.(model.γ)) - sum(lgamma.(sum(model.γ, dims=1)))
    lnq -= sum((model.γ .- 1) .* model.Elnθ)
    return lnq
end

function calculate_ElnQZ(model::ILDA)
    lnq = 0.0
    for d in 1:model.D
        lnq += sum(log.(model.ϕ[d] .^ model.ϕ[d]))
    end
    return lnq
end

function calculate_elbo(model::ILDA)
    elbo = 0.0
    elbo += calculate_ElnPβ(model)
    elbo += calculate_ElnPθ(model)
    elbo += calculate_ElnPZ(model)
    elbo += calculate_ElnPX(model)
    elbo -= calculate_ElnQβ(model)
    elbo -= calculate_ElnQθ(model)
    elbo -= calculate_ElnQZ(model)
    return elbo
end

function calculate_loglikelihood(X::Vector{Matrix{Int}}, features::Matrix{Int},
                                 θ::Matrix{Float64}, β::Vector{Matrix{Float64}})
    ll = 0.0

    K = size(θ, 1)
    I = size(features, 2)
    N = 0
    for d in 1:length(X)
        N += sum(X[d][:, 2])

        for w in 1:size(X[d], 1)
            v = X[d][w, 1]

            pw = 0.0
            for k in 1:K
                tmp = θ[k, d]
                for i in 1:I
                    j = features[v, i]
                    tmp *= β[i][j, k]
                end
                pw += tmp
            end
            ll += X[d][w, 2] * log(pw)
        end
    end

    return ll / N
end

function calculate_loglikelihood(X::Vector{Matrix{Int}}, model::ILDA)
    return calculate_loglikelihood(X, model.features, model.θ, model.β)
end

function calculate_loglikelihood(model::ILDA)
    return calculate_loglikelihood(model.X, model.features, model.θ, model.β)
end

function fit!(model::ILDA; maxiter=1000, tol=1e-4, verbose=true)
    ll = Float64[]

    for iter in 1:maxiter
        update_γ!(model)
        update_ϕ!(model)
        update_λ!(model)

        update_β!(model)
        update_θ!(model)

        push!(ll, calculate_loglikelihood(model))

        if verbose
            println("$iter\tLog-likelihood: ", ll[end])
        end

        if length(ll) > 10 && check_convergence(ll, tol=tol)
            model.converged = true
            break
        end
    end
    model.elbo = calculate_elbo(model)
    model.ll = ll[end]

    return ll
end

function unsmoothed_update_ϕ!(model::ILDA)
    for d in 1:model.D
        for w in 1:size(model.X[d], 1)
            model.ϕ[d][:, w] .= exp.(model.Elnθ[:, d])

            v = model.X[d][w, 1]
            for i in 1:model.I
                j = model.features[v, i]
                model.ϕ[d][:, w] .*= model.β[i][j, :]
            end
        end

        model.ϕ[d] .= model.ϕ[d] ./ sum(model.ϕ[d], dims=1)
    end
end

function transform(model::ILDA, X::Vector{Matrix{Int}};
                   maxiter=1000, tol=1e-4, verbose=false)

    newmodel = LDA(model.K, model.α, model.η, X)
    newmodel.β = deepcopy(model.β)

    ll = Float64[]

    for iter in 1:maxiter
        update_γ!(newmodel)
        unsmoothed_update_ϕ!(newmodel)

        update_θ!(newmodel)

        push!(ll, calculate_loglikelihood(newmodel))

        if verbose
            println("$iter\tLog-likelihood: ", ll[end])
        end

        if length(ll) > 10 && check_convergence(ll, tol=tol)
            newmodel.converged = true
            break
        end
    end

    if !newmodel.converged
        warn("transform did not converge")
    end

    return newmodel.θ
end

function fit_heldout(Xheldout::Vector{Matrix{Int}}, model::ILDA;
        maxiter=100, verbose=false)

    heldout_model = ILDA(model.K, model.α, model.η, model.features, Xheldout)
    heldout_model.λ = deepcopy(model.λ)
    heldout_model.β = deepcopy(model.β)
    heldout_model.Elnβ = deepcopy(model.Elnβ)

    ll = Float64[]
    for iter in 1:maxiter
        update_γ!(heldout_model)
        update_ϕ!(heldout_model)

        update_θ!(heldout_model)

        push!(ll, calculate_loglikelihood(heldout_model))

        if verbose
            println("$iter\tLog-likelihood: ", ll[end])
        end

        if length(ll) > 10 && check_convergence(ll)
            heldout_model.converged = true
            break
        end
    end
    heldout_model.elbo = calculate_elbo(heldout_model)
    heldout_model.ll = ll[end]

    return heldout_model
end


text:
#==============================================================================
For a given measurement vector, this example computes a MAP estimate, all global
minimizers of the log-likelihood, posterior means, and posterior standard
deviations.

In addition, we visualize the posterior density and the entropy of the
conditional posterior of M as functions of w.

Additional packages may be needed to install for the visualization.

This is the second example presented in Section 4.2 of the article:

L. Mustonen, X. Gao, A. Santana, R. Mitchell, Y. Vigfusson, and L. Ruthotto,
A Bayesian framework for molecular strain identification from mixed diagnostic samples
https://arxiv.org/abs/1803.02916
==============================================================================#

using StrainRecon
using PyPlot
using jInv.Utils
using jInv.Mesh
using jInvVis
srand(0) # initialize the random number generator for reproducibility

# define the measurement vector
d = [.1,.3,.5,.6]

# dimensions of the problem (the number of classes is implicitly p=2)
m = length(d)
n = 3

R = 100 # inverse square root of the diagonal entry in the matrix Gamma
param = getStrainReconParam(d,n,R=R)

# generate all possible binary matrices of size (m,n)
Mbits = [[bit=='1' for bit in bin(k,m*n)] for k in 0:2^(m*n)-1]
Mall = true*reshape(hcat(Mbits...),(m,n,2^(m*n)))
N = size(Mall,3)

# compute minimizing frequency vectors and the corresponding misfits for each
# binary matrix
wmins = zeros(n,N)
misfits = zeros(N)
for idx = 1:N
    wmins[:,idx] = findw(param,Mall[:,:,idx])
    misfits[idx] = misfit(param,Mall[:,:,idx],wmins[:,idx])[1]
end

# sort misfits in ascending order
srt = sortperm(misfits)

# print all pairs (M,w) that correspond to the minimum misfit
for idx = 1:N
    if abs(misfits[srt[idx]] - misfits[srt[1]]) > 1e-8
        break # not a global minimum anymore
    end
    Base.showarray(STDOUT,1*Mall[:,:,srt[idx]],false)
    println()
    println(wmins[:,srt[idx]])
end
println()

# compute posterior means and variances
Mcm,wcm,dummy,Mvar,wvar = wrapperIntegratePosterior(param)

println("Posterior means:")
Base.showarray(STDOUT,Mcm,false)
println()
println(wcm)
println()

println("Posterior standard deviations:")
Base.showarray(STDOUT,sqrt.(Mvar),false)
println()
println(sqrt.(wvar))
println()

# visualize posterior (reduced misfit) and entropy
Ma = getRegularMesh([1/3 1 0 1/2 ],3*[80 80])
as = getCellCenteredGrid(Ma)
redmisfits = zeros(Ma.nc)
redmisfits[:] = NaN # white background
entropy = zeros(Ma.nc)
entropy[:] = NaN

for k=1:size(as,1)
    at = vec([ vec(as[k,:]); 1-sum(as[k,:])])
    if all(at.>=0) && sum(at)<=1 && all(diff(at).<=0)
        redmisfits[k] = 2*reducedMisfit(param,at)[1] # drop the coefficient 1/2
        entropy[k] = entropyM(param,at)
    end
end

subplot(1,2,1)
viewImage2D(redmisfits,Ma)
xlabel("w1")
ylabel("w2")
axis("image")
title("Reduced misfit")
colorbar()
set_cmap("Spectral")
clim(0,240)

subplot(1,2,2)
viewImage2D(entropy,Ma)
xlabel("w1")
ylabel("w2")
axis("image")
title("Entropy")
colorbar()
set_cmap("Spectral")
clim(0,3)

# fix w and compute matrix minimizers
w0 = [0.5; 0.3; 0.2]
misfits0 = zeros(N)
for idx = 1:N
    misfits0[idx] = misfit(param,Mall[:,:,idx],w0)[1]
end

# sort misfits in ascending order
srt0 = sortperm(misfits0)

# print all matrices that correspond to the minimum misfit for the given w
for idx = 1:N
    if abs(misfits0[srt0[idx]] - misfits0[srt0[1]]) > 1e-8
        break # not a minimizing matrix anymore
    end
    Base.showarray(STDOUT,1*Mall[:,:,srt0[idx]],false)
    println()
end

ext: jl


text:
struct ConsistencyTest{E<:CalibrationErrorEstimator,V} <: CalibrationTest
    """Calibration estimator."""
    estimator::E
    """Calibration estimate."""
    estimate::V
    """Bootstrap samples using consistency resampling."""
    straps::Vector{V}
end

function ConsistencyTest(estimator::CalibrationErrorEstimator,
                         data::Tuple{<:AbstractMatrix{<:Real},<:AbstractVector{<:Integer}};
                         nruns::Int = 1_000,
                         rng::AbstractRNG = Random.GLOBAL_RNG)
    # check if the number of predictions and labels is equal
    predictions, labels = CalibrationErrors.get_predictions_labels(data)

    # perform consistency resampling and evaluate calibration error
    estimate, straps = consistency_resampling(rng, estimator, predictions, labels, nruns)

    ConsistencyTest(estimator, estimate, straps)
end

HypothesisTests.testname(::ConsistencyTest) =
    "Calibration test based on consistency resampling"

function HypothesisTests.pvalue(test::ConsistencyTest)
    # count number of sampled estimates greater or equal to the observed estimate
    estimate = test.estimate
    n = 0
    for sample in test.straps
        sample ≥ estimate && (n += 1)
    end

    n / length(test.straps)
end

# TODO: move the following code upstream
function consistency_resampling(rng, estimator, predictions, labels, nruns)
    nclasses, nsamples = size(predictions)

    # use same heuristic as StatsBase to decide whether to sample predictions
    # directly or to build an alias table
    # TODO: needs proper benchmarking
    if nsamples < 40
        consistency_resampling_direct(rng, estimator, predictions, labels, nruns)
    else
        t = nsamples < 500 ? 64 : 32
        if nclasses < t
            consistency_resampling_direct(rng, estimator, predictions, labels, nruns)
        else
            consistency_resampling_alias(rng, estimator, predictions, labels, nruns)
        end
    end
end

# sample labels directly (without alias table)
function consistency_resampling_direct(rng::AbstractRNG, estimator, predictions, labels,
                                       nruns)
    nclasses, nsamples = size(predictions)

    # evaluate statistic
    estimate = calibrationerror(estimator, (predictions, labels))

    # create caches
    resampled_predictions = similar(predictions)
    resampled_labels = similar(labels)
    resampled_data = (resampled_predictions, resampled_labels)

    # create sampler
    sp = Random.Sampler(rng, Base.OneTo(nsamples))

    # create output
    straps = Vector{typeof(estimate)}(undef, nruns)

    # for each resampling step
    @inbounds for i in 1:nruns
        # resample data
        for j in 1:nsamples
            # resample predictions
            idx = rand(rng, sp)
            for k in axes(predictions, 1)
                resampled_predictions[k, j] = predictions[k, idx]
            end

            # resample labels
            p = rand(rng)
            cw = resampled_predictions[1, j]
            label = 1
            while cw < p && label < nclasses
                label += 1
                cw += resampled_predictions[label, j]
            end
            resampled_labels[j] = label
        end

        # evaluate calibration error
        straps[i] = calibrationerror(estimator, resampled_data)
    end

    estimate, straps
end

# sample labels with alias table
function consistency_resampling_alias(rng::AbstractRNG, estimator, predictions, labels,
                                      nruns)
    nclasses, nsamples = size(predictions)

    # evaluate statistic
    estimate = calibrationerror(estimator, (predictions, labels))

    # create alias table
    accept = [Vector{Float64}(undef, nclasses) for _ in 1:nsamples]
    alias = [Vector{Int}(undef, nclasses) for _ in 1:nsamples]
    @inbounds for i in axes(predictions, 2)
        StatsBase.make_alias_table!(view(predictions, :, i), 1.0, accept[i], alias[i])
    end

    # create sampler of labels
    splabels = Random.Sampler(rng, Base.OneTo(nclasses))

    # create caches
    resampled_predictions = similar(predictions)
    resampled_labels = similar(labels)
    resampled_data = (resampled_predictions, resampled_labels)

    # create sampler
    sp = Random.Sampler(rng, Base.OneTo(nsamples))

    # create output
    straps = Vector{typeof(estimate)}(undef, nruns)

    # for each resampling step
    @inbounds for i in 1:nruns
        # resample data
        for j in 1:nsamples
            # resample predictions
            idx = rand(rng, sp)
            for k in axes(predictions, 1)
                resampled_predictions[k, j] = predictions[k, idx]
            end

            # resample labels
            l = rand(rng, splabels)
            resampled_labels[j] = rand(rng) < accept[j][l] ? l : alias[j][l]
        end

        # evaluate calibration error
        straps[i] = calibrationerror(estimator, resampled_data)
    end

    estimate, straps
end


text:
using AtomicMiscellany
using Documenter, Literate
using Plots, Formatting

# Stop GR from trying to connect to X when producing plots.
# This will make the documentation build about x10 faster and gets rid of all the
# GR error messages from the build logs.
# https://discourse.julialang.org/t/generation-of-documentation-fails-qt-qpa-xcb-could-not-connect-to-display/60988/2
ENV["GKSwstype"] = "100"

DocMeta.setdocmeta!(AtomicMiscellany, :DocTestSetup, quote
    using AtomicMiscellany, AtomicMiscellany.NuclearModels
end; recursive=true)

Literate.markdown(
    joinpath(@__DIR__, "src", "implementation", "hydrogenic.jl"),
    joinpath(@__DIR__, "src", "implementation"),
    documenter = true,
)

makedocs(;
    modules=[AtomicMiscellany],
    authors="JuliaAtoms contributors",
    repo="https://github.com/JuliaAtoms/AtomicMiscellany.jl/blob/{commit}{path}#{line}",
    sitename="AtomicMiscellany.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaAtoms.github.io/AtomicMiscellany.jl",
        assets=String[],
        collapselevel = 1,
    ),
    pages=[
        "Home" => "index.md",
        "nuclearmodels.md",
        "hydrogenic.md",
        "Implementation notes" => [
            "Nuclear models" => "implementation/nuclearmodels.md",
            "Hydrogenic energies" => "implementation/hydrogenic.md",
        ]
    ],
)

deploydocs(;
    repo="github.com/JuliaAtoms/AtomicMiscellany.jl",
)


text:
#
#+t
# Blalbla
#
# blibli
function foo(x::T) where {T}
    # This is the foo function body
    2+2
end


text:
export stieltjes, lanczos, mcdiscretization

function removeZeroWeights(n::Vector{Float64},w::Vector{Float64})
    nw = [n w]
    inds = findall(w->abs(w)>=eps(), nw[:,2])
    # nw = sortrows(nw, by=x->x[1])
    nw = nw[inds,:]
    nw = sortslices(nw, dims=1, by=x->x[1])
    nw[:,1], nw[:,2]
end

"""
    stieltjes(N::Int64,nodes_::Vector{Float64},weights_::Vector{Float64};removezeroweights::Bool=true)
Stieltjes procedure---Given the nodes and weights the function
generates the first`N` recurrence coefficients of the
corresponding discrete orthogonal polynomials.

Set the Boolean `removezeroweights` to `true` if zero weights should be removed.
"""
function stieltjes(N::Int64,nodes_::Vector{Float64},weights_::Vector{Float64};removezeroweights::Bool=true)
    tiny = 10*floatmin()
    huge = 0.1*floatmax()
    α, β = zeros(Float64,N), zeros(Float64,N)
    nodes, weights = removezeroweights ? removeZeroWeights(nodes_,weights_) : (nodes_, weights_)
    Ncap = length(nodes)
    @assert N > 0 && N <= Ncap "N is out of range."
    s0::Float64 = sum(weights)
    @inbounds α[1] = dot(nodes,weights)/s0
    @inbounds β[1] = s0
    N == 1 && return α, β
    p0, p1, p2 = zeros(Float64,Ncap), zeros(Float64,Ncap), ones(Float64,Ncap)
    for k=1:N-1
        p0 = p1
        p1 = p2
        @inbounds p2 = (nodes .- α[k]) .* p1 - β[k]*p0
        s1 = dot(weights,p2.^2)
        s2 = dot(nodes, weights.*p2.^2)
        # s1 = sum( weights[ν]*p2[ν]^2 for ν=1:Ncap )
        # s2 = sum( weights[ν]*nodes[ν]*p2[ν]^2 for ν=1:Ncap )
        @assert abs(s1)>=tiny "Underflow in stieltjes() for k=$k; try using `removeZeroWeights`"
        @assert maximum(abs.(p2))<=huge && abs(s2)<=huge "Overflow in stieltjes for k=$k"
        @inbounds α[k+1] = s2/s1
        @inbounds β[k+1] = s1/s0
        s0 = s1
    end
    α, β
end

"""
    lanczos(N::Int64,nodes::Vector{Float64},weights::Vector{Float64};removezeroweights::Bool=true)
Lanczos procedure---given the nodes and weights the function
generates the first `N` recurrence coefficients of the
corresponding discrete orthogonal polynomials.

Set the Boolean `removezeroweights` to `true` if zero weights should be removed.

The script is adapted from the routine RKPW in
W.B. Gragg and W.J. Harrod, *The numerically stable
reconstruction of Jacobi matrices from spectral data*,
Numer. Math. 44 (1984), 317-335.
"""
function lanczos(N::Int64,nodes::Vector{Float64},weights::Vector{Float64};removezeroweights::Bool=true)
    @assert length(nodes) == length(weights) > 0 "inconsistent number of nodes and weights"
    nodes, weights = removezeroweights ? removeZeroWeights(nodes,weights) : (nodes, weights)
    Ncap = length(nodes)
    (N<=0 || N>Ncap) && error("$N out of range")
    p0 = copy(nodes)
    p1 = zeros(Float64,Ncap);
    @inbounds p1[1] = first(weights)
    for n = 1:Ncap-1
      @inbounds pn = weights[n+1]
      gam = 1.; sig = 0.; t = 0.;
      @inbounds xlam = nodes[n+1]
      for k = 1:n+1
        @inbounds ρ = p1[k] + pn
        tmp = gam*ρ
        tsig = sig
        gam, sig = ρ > 0 ? (p1[k]/ρ, pn/ρ) : (1.,0.)
        @inbounds tk = sig * (p0[k] - xlam) - gam * t
        @inbounds p0[k] -= (tk - t);
        t = tk
        @inbounds pn = sig <= 0 ? tsig*p1[k] : (t^2)/sig
        tsig = sig
        @inbounds p1[k] = tmp;
      end
    end
    @inbounds p0[Base.OneTo(N)], p1[Base.OneTo(N)]
end

"""
    mcdiscretization(N::Int64,quads::Vector{},discretemeasure::Matrix{Float64}=zeros(0,2);discretization::Function=stieltjes,Nmax::Integer=300,ε::Float64=1e-8,gaussquad::Bool=false)
This routine returns ``N`` recurrence coefficients of the polynomials that are
orthogonal relative to a weight function ``w`` that
is decomposed as a sum of ``m`` weights ``w_i`` with domains ``[a_i,b_i]`` for ``i=1,\\dots,m``,
```math
w(t) = \\sum_{i}^{m} w_i(t) \\quad \\text{with } \\operatorname{dom}(w_i) = [a_i, b_i].
```
For each weight ``w_i`` and its domain ``[a_i, b_i]`` the function `mcdiscretization()`
expects a quadrature rule of the form
    nodes::Vector{Float64}, weights::Vector{Float64} = my_quad_i(N::Int64)
all of which are stacked in the parameter `quad`
    quad = [ my_quad_1, ..., my_quad_m ]
If the weight function has a discrete part (specified by `discretemeasure`)
it is added on to the discretized continuous weight function.

The function `mcdiscretization()` performs a sequence of discretizations of the
given weight ``w(t)``, each discretization being
followed by an application of the Stieltjes or Lanczos procedure
(keyword `discretization in [stieltjes, lanczos]`) to produce approximations to the desired recurrence
coefficients.
The function applies to each subinterval ``i`` an `N`-point quadrature rule (the ``i``th entry of `quad`)
to discretize the weight function ``w_i`` on that subinterval.
If the procedure converges to within a prescribed accuracy
`ε` before `N` reaches its maximum allowed value `Nmax`.
If the function does not converge, the function prompts an error message.

The keyword `gaussquad` should be set to `true` if Gauss quadrature rules
are available *for all* ``m`` weights ``w_i(t)`` with ``i = 1, \\dots, m``.

For further information, please see W. Gautschi "Orthogonal Polynomials: Approximation
and Computation", Section 2.2.4.
"""
function mcdiscretization(N::Int64,quads::Vector{},discretemeasure::Matrix{Float64}=zeros(0,2);discretization::Function=stieltjes,Nmax::Integer=300,ε::Float64=1e-8,gaussquad::Bool=false,removezeroweights::Bool=false)
    @assert Nmax>0 && Nmax>N "invalid choice of Nmax=$Nmax."
    @assert ε>0 "invalid choice of ε = $ε"
    @assert discretization in [stieltjes, lanczos] "unknown discretization $discretization"
    @assert length(quads)>0 "no quadrature rule specified"
    δ = !gaussquad ? 1 : 2
    # δ::Int64=1
    # gaussquad ? δ=2 : ()
    mc, mp = size(quads,1), size(discretemeasure,1)
    Δ::Int64=1; kount::Int64=-1;
    α, β, b = zeros(Float64,N), zeros(Float64,N), ones(Float64,N)
    Mi = floor(Int64,(2*N-1)/δ);
    while any(abs.(β-b).>ε*abs.(β))
        b = copy(β);
        kount += 1
        if kount > 1 Δ=2^(floor(Int64,kount/5))*N end
        # kount>1 ? Δ=2^(floor(Int64,kount/5))*N : ()
        Mi+=Δ;
        Mi > Nmax && error("Mi=$Mi exceeds Nmax=$Nmax")
        Ntot = mc*Mi;
        xx, ww = zeros(Float64,Ntot), zeros(Float64,Ntot)
        for i=1:mc
            nn=(i-1)*Mi;
            @inbounds x, w = quads[i](Mi)
            @inbounds xx[nn+1:nn+Mi], ww[nn+1:nn+Mi] = x, w
        end
        if mp>0
            @inbounds xx[Ntot+1:Ntot+mp], ww[Ntot+1:Ntot+mp] = discretemeasure[:,1], discretemeasure[:,2]
        end
        α, β = discretization(N,xx,ww;removezeroweights=removezeroweights)
        end
    # printstyled("\nSuccess: ",color=:green)
    # print("converged after $kount iteration(s).\n\n")
    return α, β
end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(x^5*(a+b*arctan(c*x)), x) == :(1//6*x^6*(a+b*arctan(c*x))+-1//6*b*x*c^-5+-1//30*b*c^-1*x^5+1//6*b*c^-6*arctan(c*x)+1//18*b*c^-3*x^3)
@test integrate(x^4*(a+b*arctan(c*x)), x) == :(1//5*x^5*(a+b*arctan(c*x))+-1//10*b*c^-5*log(1+c^2*x^2)+-1//20*b*c^-1*x^4+1//10*b*c^-3*x^2)
@test integrate(x^3*(a+b*arctan(c*x)), x) == :(1//4*x^4*(a+b*arctan(c*x))+-1//4*b*c^-4*arctan(c*x)+-1//12*b*c^-1*x^3+1//4*b*x*c^-3)
@test integrate(x^2*(a+b*arctan(c*x)), x) == :(1//3*x^3*(a+b*arctan(c*x))+-1//6*b*c^-1*x^2+1//6*b*c^-3*log(1+c^2*x^2))
@test integrate(x*(a+b*arctan(c*x)), x) == :((1/2)*x^2*(a+b*arctan(c*x))+(1/2)*b*c^-2*arctan(c*x)+-1//2*b*x*c^-1)
@test integrate(a+b*arctan(c*x), x) == :(a*x+b*x*arctan(c*x)+-1//2*b*c^-1*log(1+c^2*x^2))
@test integrate(x^-1*(a+b*arctan(c*x)), x) == :(a*log(x)+(1/2)*im*b*Polylogarithms.polylog(2,-1*im*c*x)+-1//2*im*b*Polylogarithms.polylog(2,im*c*x))
@test integrate(x^-2*(a+b*arctan(c*x)), x) == :(-1*x^-1*(a+b*arctan(c*x))+b*c*log(x)+-1//2*b*c*log(1+c^2*x^2))
@test integrate(x^-3*(a+b*arctan(c*x)), x) == :(-1//2*x^-2*(a+b*arctan(c*x))+-1//2*b*c*x^-1+-1//2*b*c^2*arctan(c*x))
@test integrate(x^-4*(a+b*arctan(c*x)), x) == :(-1//3*x^-3*(a+b*arctan(c*x))+-1//3*b*c^3*log(x)+-1//6*b*c*x^-2+1//6*b*c^3*log(1+c^2*x^2))
@test integrate(x^-5*(a+b*arctan(c*x)), x) == :(-1//4*x^-4*(a+b*arctan(c*x))+-1//12*b*c*x^-3+1//4*b*c^3*x^-1+1//4*b*c^4*arctan(c*x))
@test integrate(x^-6*(a+b*arctan(c*x)), x) == :(-1//5*x^-5*(a+b*arctan(c*x))+-1//10*b*c^5*log(1+c^2*x^2)+-1//20*b*c*x^-4+1//5*b*c^5*log(x)+1//10*b*c^3*x^-2)
@test integrate(x^5*(a+b*arctan(c*x))^2, x) == :(1//6*c^-6*(a+b*arctan(c*x))^2+1//6*x^6*(a+b*arctan(c*x))^2+-4//45*b^2*c^-4*x^2+1//60*b^2*c^-2*x^4+23//90*b^2*c^-6*log(1+c^2*x^2)+-1//3*a*b*x*c^-5+-1//3*x*b^2*c^-5*arctan(c*x)+-1//15*b*c^-1*x^5*(a+b*arctan(c*x))+1//9*b*c^-3*x^3*(a+b*arctan(c*x)))
@test integrate(x^4*(a+b*arctan(c*x))^2, x) == :(1//5*x^5*(a+b*arctan(c*x))^2+-3//10*x*b^2*c^-4+1//5*im*c^-5*(a+b*arctan(c*x))^2+1//30*b^2*c^-2*x^3+3//10*b^2*c^-5*arctan(c*x)+-1//10*b*c^-1*x^4*(a+b*arctan(c*x))+1//5*im*b^2*c^-5*Polylogarithms.polylog(2,1+-2*(1+im*c*x)^-1)+1//5*b*c^-3*x^2*(a+b*arctan(c*x))+2//5*b*c^-5*(a+b*arctan(c*x))*log(2*(1+im*c*x)^-1))
@test integrate(x^3*(a+b*arctan(c*x))^2, x) == :(-1//4*c^-4*(a+b*arctan(c*x))^2+1//4*x^4*(a+b*arctan(c*x))^2+-1//3*b^2*c^-4*log(1+c^2*x^2)+1//12*b^2*c^-2*x^2+(1/2)*a*b*x*c^-3+(1/2)*x*b^2*c^-3*arctan(c*x)+-1//6*b*c^-1*x^3*(a+b*arctan(c*x)))
@test integrate(x^2*(a+b*arctan(c*x))^2, x) == :(1//3*x^3*(a+b*arctan(c*x))^2+-1//3*im*c^-3*(a+b*arctan(c*x))^2+-1//3*b^2*c^-3*arctan(c*x)+1//3*x*b^2*c^-2+-2//3*b*c^-3*(a+b*arctan(c*x))*log(2*(1+im*c*x)^-1)+-1//3*im*b^2*c^-3*Polylogarithms.polylog(2,1+-2*(1+im*c*x)^-1)+-1//3*b*c^-1*x^2*(a+b*arctan(c*x)))
@test integrate(x*(a+b*arctan(c*x))^2, x) == :((1/2)*c^-2*(a+b*arctan(c*x))^2+(1/2)*x^2*(a+b*arctan(c*x))^2+(1/2)*b^2*c^-2*log(1+c^2*x^2)+-1*a*b*x*c^-1+-1*x*b^2*c^-1*arctan(c*x))
@test integrate((a+b*arctan(c*x))^2, x) == :(x*(a+b*arctan(c*x))^2+im*c^-1*(a+b*arctan(c*x))^2+im*b^2*c^-1*Polylogarithms.polylog(2,1+-2*(1+im*c*x)^-1)+2*b*c^-1*(a+b*arctan(c*x))*log(2*(1+im*c*x)^-1))
@test integrate(x^-1*(a+b*arctan(c*x))^2, x) == :((1/2)*b^2*Polylogarithms.polylog(3,-1+2*(1+im*c*x)^-1)+2*(a+b*arctan(c*x))^2*arctanh(1+-2*(1+im*c*x)^-1)+-1//2*b^2*Polylogarithms.polylog(3,1+-2*(1+im*c*x)^-1)+im*b*(a+b*arctan(c*x))*Polylogarithms.polylog(2,-1+2*(1+im*c*x)^-1)+-1*im*b*(a+b*arctan(c*x))*Polylogarithms.polylog(2,1+-2*(1+im*c*x)^-1))
@test integrate(x^-2*(a+b*arctan(c*x))^2, x) == :(-1*x^-1*(a+b*arctan(c*x))^2+-1*im*c*(a+b*arctan(c*x))^2+-1*im*c*b^2*Polylogarithms.polylog(2,-1+2*(1+-1*im*c*x)^-1)+2*b*c*(a+b*arctan(c*x))*log(2+-2*(1+-1*im*c*x)^-1))
@test integrate(x^-3*(a+b*arctan(c*x))^2, x) == :(-1//2*c^2*(a+b*arctan(c*x))^2+-1//2*x^-2*(a+b*arctan(c*x))^2+b^2*c^2*log(x)+-1//2*b^2*c^2*log(1+c^2*x^2)+-1*b*c*x^-1*(a+b*arctan(c*x)))
@test integrate(x^-4*(a+b*arctan(c*x))^2, x) == :(-1//3*x^-3*(a+b*arctan(c*x))^2+-1//3*b^2*c^2*x^-1+-1//3*b^2*c^3*arctan(c*x)+1//3*im*c^3*(a+b*arctan(c*x))^2+-2//3*b*c^3*(a+b*arctan(c*x))*log(2+-2*(1+-1*im*c*x)^-1)+-1//3*b*c*x^-2*(a+b*arctan(c*x))+1//3*im*b^2*c^3*Polylogarithms.polylog(2,-1+2*(1+-1*im*c*x)^-1))
@test integrate(x^-5*(a+b*arctan(c*x))^2, x) == :(-1//4*x^-4*(a+b*arctan(c*x))^2+1//4*c^4*(a+b*arctan(c*x))^2+-2//3*b^2*c^4*log(x)+-1//12*b^2*c^2*x^-2+1//3*b^2*c^4*log(1+c^2*x^2)+(1/2)*b*c^3*x^-1*(a+b*arctan(c*x))+-1//6*b*c*x^-3*(a+b*arctan(c*x)))
@test integrate(x^5*(a+b*arctan(c*x))^3, x) == :(1//6*c^-6*(a+b*arctan(c*x))^3+1//6*x^6*(a+b*arctan(c*x))^3+-19//60*b^3*c^-6*arctan(c*x)+-1//60*b^3*c^-3*x^3+19//60*x*b^3*c^-5+-23//15*b^2*c^-6*(a+b*arctan(c*x))*log(2*(1+im*c*x)^-1)+-23//30*im*b*c^-6*(a+b*arctan(c*x))^2+-23//30*im*b^3*c^-6*Polylogarithms.polylog(2,1+-2*(1+im*c*x)^-1)+-4//15*b^2*c^-4*x^2*(a+b*arctan(c*x))+-1//2*b*x*c^-5*(a+b*arctan(c*x))^2+-1//10*b*c^-1*x^5*(a+b*arctan(c*x))^2+1//6*b*c^-3*x^3*(a+b*arctan(c*x))^2+1//20*b^2*c^-2*x^4*(a+b*arctan(c*x)))
@test integrate(x^4*(a+b*arctan(c*x))^3, x) == :(1//5*x^5*(a+b*arctan(c*x))^3+(1/2)*b^3*c^-5*log(1+c^2*x^2)+-1//20*b^3*c^-3*x^2+1//5*im*c^-5*(a+b*arctan(c*x))^3+3//10*b^3*c^-5*Polylogarithms.polylog(3,1+-2*(1+im*c*x)^-1)+9//20*b*c^-5*(a+b*arctan(c*x))^2+-9//10*a*x*b^2*c^-4+-9//10*x*b^3*c^-4*arctan(c*x)+-3//20*b*c^-1*x^4*(a+b*arctan(c*x))^2+1//10*b^2*c^-2*x^3*(a+b*arctan(c*x))+3//5*b*c^-5*(a+b*arctan(c*x))^2*log(2*(1+im*c*x)^-1)+3//10*b*c^-3*x^2*(a+b*arctan(c*x))^2+3//5*im*b^2*c^-5*(a+b*arctan(c*x))*Polylogarithms.polylog(2,1+-2*(1+im*c*x)^-1))
@test integrate(x^3*(a+b*arctan(c*x))^3, x) == :(-1//4*c^-4*(a+b*arctan(c*x))^3+1//4*x^4*(a+b*arctan(c*x))^3+-1//4*x*b^3*c^-3+1//4*b^3*c^-4*arctan(c*x)+im*b*c^-4*(a+b*arctan(c*x))^2+im*b^3*c^-4*Polylogarithms.polylog(2,1+-2*(1+im*c*x)^-1)+2*b^2*c^-4*(a+b*arctan(c*x))*log(2*(1+im*c*x)^-1)+-1//4*b*c^-1*x^3*(a+b*arctan(c*x))^2+1//4*b^2*c^-2*x^2*(a+b*arctan(c*x))+3//4*b*x*c^-3*(a+b*arctan(c*x))^2)
@test integrate(x^2*(a+b*arctan(c*x))^3, x) == :(1//3*x^3*(a+b*arctan(c*x))^3+-1//2*b*c^-3*(a+b*arctan(c*x))^2+-1//2*b^3*c^-3*Polylogarithms.polylog(3,1+-2*(1+im*c*x)^-1)+-1//2*b^3*c^-3*log(1+c^2*x^2)+-1//3*im*c^-3*(a+b*arctan(c*x))^3+a*x*b^2*c^-2+x*b^3*c^-2*arctan(c*x)+-1*b*c^-3*(a+b*arctan(c*x))^2*log(2*(1+im*c*x)^-1)+-1//2*b*c^-1*x^2*(a+b*arctan(c*x))^2+-1*im*b^2*c^-3*(a+b*arctan(c*x))*Polylogarithms.polylog(2,1+-2*(1+im*c*x)^-1))
@test integrate(x*(a+b*arctan(c*x))^3, x) == :((1/2)*c^-2*(a+b*arctan(c*x))^3+(1/2)*x^2*(a+b*arctan(c*x))^3+-3*b^2*c^-2*(a+b*arctan(c*x))*log(2*(1+im*c*x)^-1)+-3//2*im*b*c^-2*(a+b*arctan(c*x))^2+-3//2*im*b^3*c^-2*Polylogarithms.polylog(2,1+-2*(1+im*c*x)^-1)+-3//2*b*x*c^-1*(a+b*arctan(c*x))^2)
@test integrate((a+b*arctan(c*x))^3, x) == :(x*(a+b*arctan(c*x))^3+im*c^-1*(a+b*arctan(c*x))^3+3//2*b^3*c^-1*Polylogarithms.polylog(3,1+-2*(1+im*c*x)^-1)+3*b*c^-1*(a+b*arctan(c*x))^2*log(2*(1+im*c*x)^-1)+3*im*b^2*c^-1*(a+b*arctan(c*x))*Polylogarithms.polylog(2,1+-2*(1+im*c*x)^-1))
@test integrate(x^-1*(a+b*arctan(c*x))^3, x) == :(2*(a+b*arctan(c*x))^3*arctanh(1+-2*(1+im*c*x)^-1)+-3//2*b^2*(a+b*arctan(c*x))*Polylogarithms.polylog(3,1+-2*(1+im*c*x)^-1)+-3//4*im*b^3*Polylogarithms.polylog(4,-1+2*(1+im*c*x)^-1)+3//2*b^2*(a+b*arctan(c*x))*Polylogarithms.polylog(3,-1+2*(1+im*c*x)^-1)+3//4*im*b^3*Polylogarithms.polylog(4,1+-2*(1+im*c*x)^-1)+-3//2*im*b*(a+b*arctan(c*x))^2*Polylogarithms.polylog(2,1+-2*(1+im*c*x)^-1)+3//2*im*b*(a+b*arctan(c*x))^2*Polylogarithms.polylog(2,-1+2*(1+im*c*x)^-1))
@test integrate(x^-2*(a+b*arctan(c*x))^3, x) == :(-1*x^-1*(a+b*arctan(c*x))^3+-1*im*c*(a+b*arctan(c*x))^3+3//2*c*b^3*Polylogarithms.polylog(3,-1+2*(1+-1*im*c*x)^-1)+3*b*c*(a+b*arctan(c*x))^2*log(2+-2*(1+-1*im*c*x)^-1)+-3*im*c*b^2*(a+b*arctan(c*x))*Polylogarithms.polylog(2,-1+2*(1+-1*im*c*x)^-1))
@test integrate(x^-3*(a+b*arctan(c*x))^3, x) == :(-1//2*c^2*(a+b*arctan(c*x))^3+-1//2*x^-2*(a+b*arctan(c*x))^3+3*b^2*c^2*(a+b*arctan(c*x))*log(2+-2*(1+-1*im*c*x)^-1)+-3//2*im*b*c^2*(a+b*arctan(c*x))^2+-3//2*im*b^3*c^2*Polylogarithms.polylog(2,-1+2*(1+-1*im*c*x)^-1)+-3//2*b*c*x^-1*(a+b*arctan(c*x))^2)
@test integrate(x^-4*(a+b*arctan(c*x))^3, x) == :(-1//3*x^-3*(a+b*arctan(c*x))^3+b^3*c^3*log(x)+-1//2*b*c^3*(a+b*arctan(c*x))^2+-1//2*b^3*c^3*Polylogarithms.polylog(3,-1+2*(1+-1*im*c*x)^-1)+-1//2*b^3*c^3*log(1+c^2*x^2)+1//3*im*c^3*(a+b*arctan(c*x))^3+-1*b*c^3*(a+b*arctan(c*x))^2*log(2+-2*(1+-1*im*c*x)^-1)+-1*b^2*c^2*x^-1*(a+b*arctan(c*x))+-1//2*b*c*x^-2*(a+b*arctan(c*x))^2+im*b^2*c^3*(a+b*arctan(c*x))*Polylogarithms.polylog(2,-1+2*(1+-1*im*c*x)^-1))
@test integrate(x^-5*(a+b*arctan(c*x))^3, x) == :(-1//4*x^-4*(a+b*arctan(c*x))^3+1//4*c^4*(a+b*arctan(c*x))^3+-1//4*b^3*c^3*x^-1+-1//4*b^3*c^4*arctan(c*x)+im*b*c^4*(a+b*arctan(c*x))^2+im*b^3*c^4*Polylogarithms.polylog(2,-1+2*(1+-1*im*c*x)^-1)+-2*b^2*c^4*(a+b*arctan(c*x))*log(2+-2*(1+-1*im*c*x)^-1)+-1//4*b*c*x^-3*(a+b*arctan(c*x))^2+-1//4*b^2*c^2*x^-2*(a+b*arctan(c*x))+3//4*b*c^3*x^-1*(a+b*arctan(c*x))^2)
@test integrate(x^(1/2)*arctan(x), x) == :(-4//3*x^(1/2)+-1//3*2^(1/2)*arctan(1+-1*2^(1/2)*x^(1/2))+-1//6*2^(1/2)*log(1+x+-1*2^(1/2)*x^(1/2))+1//3*2^(1/2)*arctan(1+2^(1/2)*x^(1/2))+1//6*2^(1/2)*log(1+x+2^(1/2)*x^(1/2))+2//3*x^3//2*arctan(x))
@test integrate(x^7*(a+b*arctan(c*x^2)), x) == :(1//8*x^8*(a+b*arctan(c*x^2))+-1//8*b*c^-4*arctan(c*x^2)+-1//24*b*c^-1*x^6+1//8*b*c^-3*x^2)
@test integrate(x^5*(a+b*arctan(c*x^2)), x) == :(1//6*x^6*(a+b*arctan(c*x^2))+-1//12*b*c^-1*x^4+1//12*b*c^-3*log(1+c^2*x^4))
@test integrate(x^3*(a+b*arctan(c*x^2)), x) == :(1//4*x^4*(a+b*arctan(c*x^2))+-1//4*b*c^-1*x^2+1//4*b*c^-2*arctan(c*x^2))
@test integrate(x*(a+b*arctan(c*x^2)), x) == :((1/2)*x^2*(a+b*arctan(c*x^2))+-1//4*b*c^-1*log(1+c^2*x^4))
@test integrate(x^-1*(a+b*arctan(c*x^2)), x) == :(a*log(x)+-1//4*im*b*Polylogarithms.polylog(2,im*c*x^2)+1//4*im*b*Polylogarithms.polylog(2,-1*im*c*x^2))
@test integrate(x^-3*(a+b*arctan(c*x^2)), x) == :(-1//2*x^-2*(a+b*arctan(c*x^2))+b*c*log(x)+-1//4*b*c*log(1+c^2*x^4))
@test integrate(x^-5*(a+b*arctan(c*x^2)), x) == :(-1//4*x^-4*(a+b*arctan(c*x^2))+-1//4*b*c*x^-2+-1//4*b*c^2*arctan(c*x^2))
@test integrate(x^-7*(a+b*arctan(c*x^2)), x) == :(-1//6*x^-6*(a+b*arctan(c*x^2))+-1//3*b*c^3*log(x)+-1//12*b*c*x^-4+1//12*b*c^3*log(1+c^2*x^4))
@test integrate(x^4*(a+b*arctan(c*x^2)), x) == :(1//5*x^5*(a+b*arctan(c*x^2))+-2//15*b*c^-1*x^3+-1//10*b*2^(1/2)*c^-5//2*arctan(1+-1*x*2^(1/2)*c^(1/2))+-1//20*b*2^(1/2)*c^-5//2*log(1+c*x^2+x*2^(1/2)*c^(1/2))+1//10*b*2^(1/2)*c^-5//2*arctan(1+x*2^(1/2)*c^(1/2))+1//20*b*2^(1/2)*c^-5//2*log(1+c*x^2+-1*x*2^(1/2)*c^(1/2)))
@test integrate(x^2*(a+b*arctan(c*x^2)), x) == :(1//3*x^3*(a+b*arctan(c*x^2))+-2//3*b*x*c^-1+-1//6*b*2^(1/2)*c^-3//2*arctan(1+-1*x*2^(1/2)*c^(1/2))+-1//12*b*2^(1/2)*c^-3//2*log(1+c*x^2+-1*x*2^(1/2)*c^(1/2))+1//6*b*2^(1/2)*c^-3//2*arctan(1+x*2^(1/2)*c^(1/2))+1//12*b*2^(1/2)*c^-3//2*log(1+c*x^2+x*2^(1/2)*c^(1/2)))
@test integrate(a+b*arctan(c*x^2), x) == :(a*x+b*x*arctan(c*x^2)+(1/2)*b*2^(1/2)*c^-1//2*arctan(1+-1*x*2^(1/2)*c^(1/2))+-1//2*b*2^(1/2)*c^-1//2*arctan(1+x*2^(1/2)*c^(1/2))+-1//4*b*2^(1/2)*c^-1//2*log(1+c*x^2+-1*x*2^(1/2)*c^(1/2))+1//4*b*2^(1/2)*c^-1//2*log(1+c*x^2+x*2^(1/2)*c^(1/2)))
@test integrate(x^-2*(a+b*arctan(c*x^2)), x) == :(-1*x^-1*(a+b*arctan(c*x^2))+(1/2)*b*2^(1/2)*c^(1/2)*arctan(1+x*2^(1/2)*c^(1/2))+-1//2*b*2^(1/2)*c^(1/2)*arctan(1+-1*x*2^(1/2)*c^(1/2))+-1//4*b*2^(1/2)*c^(1/2)*log(1+c*x^2+-1*x*2^(1/2)*c^(1/2))+1//4*b*2^(1/2)*c^(1/2)*log(1+c*x^2+x*2^(1/2)*c^(1/2)))
@test integrate(x^-4*(a+b*arctan(c*x^2)), x) == :(-1//3*x^-3*(a+b*arctan(c*x^2))+-2//3*b*c*x^-1+-1//6*b*2^(1/2)*c^3//2*arctan(1+x*2^(1/2)*c^(1/2))+-1//12*b*2^(1/2)*c^3//2*log(1+c*x^2+-1*x*2^(1/2)*c^(1/2))+1//6*b*2^(1/2)*c^3//2*arctan(1+-1*x*2^(1/2)*c^(1/2))+1//12*b*2^(1/2)*c^3//2*log(1+c*x^2+x*2^(1/2)*c^(1/2)))
@test integrate(x^-6*(a+b*arctan(c*x^2)), x) == :(-1//5*x^-5*(a+b*arctan(c*x^2))+-2//15*b*c*x^-3+-1//10*b*2^(1/2)*c^5//2*arctan(1+x*2^(1/2)*c^(1/2))+-1//20*b*2^(1/2)*c^5//2*log(1+c*x^2+x*2^(1/2)*c^(1/2))+1//10*b*2^(1/2)*c^5//2*arctan(1+-1*x*2^(1/2)*c^(1/2))+1//20*b*2^(1/2)*c^5//2*log(1+c*x^2+-1*x*2^(1/2)*c^(1/2)))
@test integrate(x^7*(a+b*arctan(c*x^2))^2, x) == :(-1//8*c^-4*(a+b*arctan(c*x^2))^2+1//8*x^8*(a+b*arctan(c*x^2))^2+-1//6*b^2*c^-4*log(1+c^2*x^4)+1//24*b^2*c^-2*x^4+-1//12*b*c^-1*x^6*(a+b*arctan(c*x^2))+1//4*a*b*c^-3*x^2+1//4*b^2*c^-3*x^2*arctan(c*x^2))
@test integrate(x^5*(a+b*arctan(c*x^2))^2, x) == :(1//6*x^6*(a+b*arctan(c*x^2))^2+-1//6*im*c^-3*(a+b*arctan(c*x^2))^2+-1//6*b^2*c^-3*arctan(c*x^2)+1//6*b^2*c^-2*x^2+-1//3*b*c^-3*(a+b*arctan(c*x^2))*log(2*(1+im*c*x^2)^-1)+-1//6*im*b^2*c^-3*Polylogarithms.polylog(2,1+-2*(1+im*c*x^2)^-1)+-1//6*b*c^-1*x^4*(a+b*arctan(c*x^2)))
@test integrate(x^3*(a+b*arctan(c*x^2))^2, x) == :(1//4*c^-2*(a+b*arctan(c*x^2))^2+1//4*x^4*(a+b*arctan(c*x^2))^2+1//4*b^2*c^-2*log(1+c^2*x^4)+-1//2*a*b*c^-1*x^2+-1//2*b^2*c^-1*x^2*arctan(c*x^2))
@test integrate(x*(a+b*arctan(c*x^2))^2, x) == :((1/2)*x^2*(a+b*arctan(c*x^2))^2+(1/2)*im*c^-1*(a+b*arctan(c*x^2))^2+b*c^-1*(a+b*arctan(c*x^2))*log(2*(1+im*c*x^2)^-1)+(1/2)*im*b^2*c^-1*Polylogarithms.polylog(2,1+-2*(1+im*c*x^2)^-1))
@test integrate(x^-1*(a+b*arctan(c*x^2))^2, x) == :((a+b*arctan(c*x^2))^2*arctanh(1+-2*(1+im*c*x^2)^-1)+-1//4*b^2*Polylogarithms.polylog(3,1+-2*(1+im*c*x^2)^-1)+1//4*b^2*Polylogarithms.polylog(3,-1+2*(1+im*c*x^2)^-1)+(1/2)*im*b*(a+b*arctan(c*x^2))*Polylogarithms.polylog(2,-1+2*(1+im*c*x^2)^-1)+-1//2*im*b*(a+b*arctan(c*x^2))*Polylogarithms.polylog(2,1+-2*(1+im*c*x^2)^-1))
@test integrate(x^-3*(a+b*arctan(c*x^2))^2, x) == :(-1//2*x^-2*(a+b*arctan(c*x^2))^2+-1//2*im*c*(a+b*arctan(c*x^2))^2+b*c*(a+b*arctan(c*x^2))*log(2+-2*(1+-1*im*c*x^2)^-1)+-1//2*im*c*b^2*Polylogarithms.polylog(2,-1+2*(1+-1*im*c*x^2)^-1))
@test integrate(x^-5*(a+b*arctan(c*x^2))^2, x) == :(-1//4*c^2*(a+b*arctan(c*x^2))^2+-1//4*x^-4*(a+b*arctan(c*x^2))^2+b^2*c^2*log(x)+-1//4*b^2*c^2*log(1+c^2*x^4)+-1//2*b*c*x^-2*(a+b*arctan(c*x^2)))
@test integrate(x^2*(a+b*arctan(c*x^2))^2, x) == :(1//12*x^3*(2a+im*b*log(1+-1*im*c*x^2))^2+-1//9*b^2*x^3*log(1+-1*im*c*x^2)+-1//12*b^2*x^3*log(1+im*c*x^2)^2+-4//3*a*b*x*c^-1+-4//3*(-1)^3//4*b^2*c^-3//2*arctanh(x*(-1)^3//4*c^(1/2))+-1//3*(-1)^3//4*b^2*c^-3//2*arctanh(x*(-1)^3//4*c^(1/2))^2+-1//6*(-1)^1//4*b^2*c^-3//2*Polylogarithms.polylog(2,1+-1*2^(1/2)*(1+x*(-1)^1//4*c^(1/2))^-1*((-1)^1//4+x*c^(1/2)))+-1//6*(-1)^1//4*b^2*c^-3//2*Polylogarithms.polylog(2,1+-1*(1+x*(-1)^1//4*c^(1/2))^-1*(1+-1im)*(1+x*(-1)^3//4*c^(1/2)))+-1//6*(-1)^3//4*b^2*c^-3//2*Polylogarithms.polylog(2,1+2^(1/2)*(1+x*(-1)^3//4*c^(1/2))^-1*((-1)^3//4+x*c^(1/2)))+-1//6*(-1)^3//4*b^2*c^-3//2*Polylogarithms.polylog(2,1+-1*(1+x*(-1)^3//4*c^(1/2))^-1*(1+im)*(1+x*(-1)^1//4*c^(1/2)))+-1//9*im*b*x^3*(2a+im*b*log(1+-1*im*c*x^2))+1//3*(-1)^1//4*b^2*c^-3//2*arctan(x*(-1)^3//4*c^(1/2))^2+1//3*(-1)^1//4*b^2*c^-3//2*Polylogarithms.polylog(2,1+-2*(1+x*(-1)^1//4*c^(1/2))^-1)+1//3*(-1)^1//4*b^2*c^-3//2*Polylogarithms.polylog(2,1+-2*(1+-1*x*(-1)^1//4*c^(1/2))^-1)+1//3*(-1)^3//4*b^2*c^-3//2*Polylogarithms.polylog(2,1+-2*(1+x*(-1)^3//4*c^(1/2))^-1)+1//3*(-1)^3//4*b^2*c^-3//2*Polylogarithms.polylog(2,1+-2*(1+-1*x*(-1)^3//4*c^(1/2))^-1)+1//6*b^2*x^3*log(1+im*c*x^2)*log(1+-1*im*c*x^2)+2//9*im*a*b*x^3+4//3*(-1)^3//4*b^2*c^-3//2*arctan(x*(-1)^3//4*c^(1/2))+-2//3*im*x*b^2*c^-1*log(1+-1*im*c*x^2)+-2//3*a*b*(-1)^1//4*c^-3//2*arctanh(x*(-1)^3//4*c^(1/2))+-2//3*(-1)^3//4*b^2*c^-3//2*arctan(x*(-1)^3//4*c^(1/2))*log(2*(1+-1*x*(-1)^1//4*c^(1/2))^-1)+-2//3*(-1)^3//4*b^2*c^-3//2*arctanh(x*(-1)^3//4*c^(1/2))*log(2*(1+x*(-1)^3//4*c^(1/2))^-1)+-1//3*im*a*b*x^3*log(1+im*c*x^2)+-1//3*b*(-1)^1//4*c^-3//2*(2a+im*b*log(1+-1*im*c*x^2))*arctan(x*(-1)^3//4*c^(1/2))+-1//3*(-1)^3//4*b^2*c^-3//2*arctan(x*(-1)^3//4*c^(1/2))*log(2^(1/2)*(1+x*(-1)^1//4*c^(1/2))^-1*((-1)^1//4+x*c^(1/2)))+-1//3*(-1)^3//4*b^2*c^-3//2*arctan(x*(-1)^3//4*c^(1/2))*log((1+x*(-1)^1//4*c^(1/2))^-1*(1+-1im)*(1+x*(-1)^3//4*c^(1/2)))+-1//3*(-1)^3//4*b^2*c^-3//2*arctanh(x*(-1)^3//4*c^(1/2))*log(1+-1*im*c*x^2)+1//3*(-1)^3//4*b^2*c^-3//2*arctan(x*(-1)^3//4*c^(1/2))*log(1+im*c*x^2)+1//3*(-1)^3//4*b^2*c^-3//2*arctanh(x*(-1)^3//4*c^(1/2))*log((1+x*(-1)^3//4*c^(1/2))^-1*(1+im)*(1+x*(-1)^1//4*c^(1/2)))+1//3*(-1)^3//4*b^2*c^-3//2*arctanh(x*(-1)^3//4*c^(1/2))*log(-1*2^(1/2)*(1+x*(-1)^3//4*c^(1/2))^-1*((-1)^3//4+x*c^(1/2)))+1//3*(-1)^3//4*b^2*c^-3//2*arctanh(x*(-1)^3//4*c^(1/2))*log(1+im*c*x^2)+2//3*im*x*b^2*c^-1*log(1+im*c*x^2)+2//3*(-1)^3//4*b^2*c^-3//2*arctan(x*(-1)^3//4*c^(1/2))*log(2*(1+x*(-1)^1//4*c^(1/2))^-1)+2//3*(-1)^3//4*b^2*c^-3//2*arctanh(x*(-1)^3//4*c^(1/2))*log(2*(1+-1*x*(-1)^3//4*c^(1/2))^-1))
@test integrate((a+b*arctan(c*x^2))^2, x) == :(x*a^2+-1//4*x*b^2*log(1+im*c*x^2)^2+-1//4*x*b^2*log(1+-1*im*c*x^2)^2+(-1)^1//4*b^2*c^-1//2*Polylogarithms.polylog(2,1+-2*(1+x*(-1)^3//4*c^(1/2))^-1)+(-1)^1//4*b^2*c^-1//2*Polylogarithms.polylog(2,1+-2*(1+-1*x*(-1)^3//4*c^(1/2))^-1)+(-1)^3//4*b^2*c^-1//2*arctan(x*(-1)^3//4*c^(1/2))^2+(-1)^3//4*b^2*c^-1//2*Polylogarithms.polylog(2,1+-2*(1+x*(-1)^1//4*c^(1/2))^-1)+(-1)^3//4*b^2*c^-1//2*Polylogarithms.polylog(2,1+-2*(1+-1*x*(-1)^1//4*c^(1/2))^-1)+(1/2)*x*b^2*log(1+im*c*x^2)*log(1+-1*im*c*x^2)+-1*(-1)^1//4*b^2*c^-1//2*arctanh(x*(-1)^3//4*c^(1/2))^2+-1//2*(-1)^1//4*b^2*c^-1//2*Polylogarithms.polylog(2,1+2^(1/2)*(1+x*(-1)^3//4*c^(1/2))^-1*((-1)^3//4+x*c^(1/2)))+-1//2*(-1)^1//4*b^2*c^-1//2*Polylogarithms.polylog(2,1+-1*(1+x*(-1)^3//4*c^(1/2))^-1*(1+im)*(1+x*(-1)^1//4*c^(1/2)))+-1//2*(-1)^3//4*b^2*c^-1//2*Polylogarithms.polylog(2,1+-1*2^(1/2)*(1+x*(-1)^1//4*c^(1/2))^-1*((-1)^1//4+x*c^(1/2)))+-1//2*(-1)^3//4*b^2*c^-1//2*Polylogarithms.polylog(2,1+-1*(1+x*(-1)^1//4*c^(1/2))^-1*(1+-1im)*(1+x*(-1)^3//4*c^(1/2)))+im*a*b*x*log(1+-1*im*c*x^2)+(-1)^1//4*b^2*c^-1//2*arctan(x*(-1)^3//4*c^(1/2))*log(2^(1/2)*(1+x*(-1)^1//4*c^(1/2))^-1*((-1)^1//4+x*c^(1/2)))+(-1)^1//4*b^2*c^-1//2*arctan(x*(-1)^3//4*c^(1/2))*log((1+x*(-1)^1//4*c^(1/2))^-1*(1+-1im)*(1+x*(-1)^3//4*c^(1/2)))+(-1)^1//4*b^2*c^-1//2*arctan(x*(-1)^3//4*c^(1/2))*log(1+-1*im*c*x^2)+(-1)^1//4*b^2*c^-1//2*arctanh(x*(-1)^3//4*c^(1/2))*log((1+x*(-1)^3//4*c^(1/2))^-1*(1+im)*(1+x*(-1)^1//4*c^(1/2)))+(-1)^1//4*b^2*c^-1//2*arctanh(x*(-1)^3//4*c^(1/2))*log(-1*2^(1/2)*(1+x*(-1)^3//4*c^(1/2))^-1*((-1)^3//4+x*c^(1/2)))+(-1)^1//4*b^2*c^-1//2*arctanh(x*(-1)^3//4*c^(1/2))*log(1+im*c*x^2)+-1*im*a*b*x*log(1+im*c*x^2)+-1*(-1)^1//4*b^2*c^-1//2*arctan(x*(-1)^3//4*c^(1/2))*log(1+im*c*x^2)+-1*(-1)^1//4*b^2*c^-1//2*arctanh(x*(-1)^3//4*c^(1/2))*log(1+-1*im*c*x^2)+-2*a*b*(-1)^3//4*c^-1//2*arctan(x*(-1)^3//4*c^(1/2))+-2*(-1)^1//4*b^2*c^-1//2*arctan(x*(-1)^3//4*c^(1/2))*log(2*(1+x*(-1)^1//4*c^(1/2))^-1)+-2*(-1)^1//4*b^2*c^-1//2*arctanh(x*(-1)^3//4*c^(1/2))*log(2*(1+x*(-1)^3//4*c^(1/2))^-1)+2*a*b*(-1)^3//4*c^-1//2*arctanh(x*(-1)^3//4*c^(1/2))+2*(-1)^1//4*b^2*c^-1//2*arctan(x*(-1)^3//4*c^(1/2))*log(2*(1+-1*x*(-1)^1//4*c^(1/2))^-1)+2*(-1)^1//4*b^2*c^-1//2*arctanh(x*(-1)^3//4*c^(1/2))*log(2*(1+-1*x*(-1)^3//4*c^(1/2))^-1))
@test integrate(x^-2*(a+b*arctan(c*x^2))^2, x) == :(-1//4*x^-1*(2a+im*b*log(1+-1*im*c*x^2))^2+1//4*b^2*x^-1*log(1+im*c*x^2)^2+(-1)^1//4*b^2*c^(1/2)*arctan(x*(-1)^3//4*c^(1/2))^2+(-1)^1//4*b^2*c^(1/2)*Polylogarithms.polylog(2,1+-2*(1+x*(-1)^1//4*c^(1/2))^-1)+(-1)^1//4*b^2*c^(1/2)*Polylogarithms.polylog(2,1+-2*(1+-1*x*(-1)^1//4*c^(1/2))^-1)+(-1)^3//4*b^2*c^(1/2)*Polylogarithms.polylog(2,1+-2*(1+x*(-1)^3//4*c^(1/2))^-1)+(-1)^3//4*b^2*c^(1/2)*Polylogarithms.polylog(2,1+-2*(1+-1*x*(-1)^3//4*c^(1/2))^-1)+-1*(-1)^3//4*b^2*c^(1/2)*arctanh(x*(-1)^3//4*c^(1/2))^2+-1//2*(-1)^1//4*b^2*c^(1/2)*Polylogarithms.polylog(2,1+-1*2^(1/2)*(1+x*(-1)^1//4*c^(1/2))^-1*((-1)^1//4+x*c^(1/2)))+-1//2*(-1)^1//4*b^2*c^(1/2)*Polylogarithms.polylog(2,1+-1*(1+x*(-1)^1//4*c^(1/2))^-1*(1+-1im)*(1+x*(-1)^3//4*c^(1/2)))+-1//2*(-1)^3//4*b^2*c^(1/2)*Polylogarithms.polylog(2,1+2^(1/2)*(1+x*(-1)^3//4*c^(1/2))^-1*((-1)^3//4+x*c^(1/2)))+-1//2*(-1)^3//4*b^2*c^(1/2)*Polylogarithms.polylog(2,1+-1*(1+x*(-1)^3//4*c^(1/2))^-1*(1+im)*(1+x*(-1)^1//4*c^(1/2)))+-1//2*b^2*x^-1*log(1+im*c*x^2)*log(1+-1*im*c*x^2)+im*a*b*x^-1*log(1+im*c*x^2)+(-1)^3//4*b^2*c^(1/2)*arctan(x*(-1)^3//4*c^(1/2))*log(1+im*c*x^2)+(-1)^3//4*b^2*c^(1/2)*arctanh(x*(-1)^3//4*c^(1/2))*log((1+x*(-1)^3//4*c^(1/2))^-1*(1+im)*(1+x*(-1)^1//4*c^(1/2)))+(-1)^3//4*b^2*c^(1/2)*arctanh(x*(-1)^3//4*c^(1/2))*log(-1*2^(1/2)*(1+x*(-1)^3//4*c^(1/2))^-1*((-1)^3//4+x*c^(1/2)))+(-1)^3//4*b^2*c^(1/2)*arctanh(x*(-1)^3//4*c^(1/2))*log(1+im*c*x^2)+-1*b*(-1)^1//4*c^(1/2)*(2a+im*b*log(1+-1*im*c*x^2))*arctan(x*(-1)^3//4*c^(1/2))+-1*(-1)^3//4*b^2*c^(1/2)*arctan(x*(-1)^3//4*c^(1/2))*log(2^(1/2)*(1+x*(-1)^1//4*c^(1/2))^-1*((-1)^1//4+x*c^(1/2)))+-1*(-1)^3//4*b^2*c^(1/2)*arctan(x*(-1)^3//4*c^(1/2))*log((1+x*(-1)^1//4*c^(1/2))^-1*(1+-1im)*(1+x*(-1)^3//4*c^(1/2)))+-1*(-1)^3//4*b^2*c^(1/2)*arctanh(x*(-1)^3//4*c^(1/2))*log(1+-1*im*c*x^2)+-2*a*b*(-1)^1//4*c^(1/2)*arctanh(x*(-1)^3//4*c^(1/2))+-2*(-1)^3//4*b^2*c^(1/2)*arctan(x*(-1)^3//4*c^(1/2))*log(2*(1+-1*x*(-1)^1//4*c^(1/2))^-1)+-2*(-1)^3//4*b^2*c^(1/2)*arctanh(x*(-1)^3//4*c^(1/2))*log(2*(1+x*(-1)^3//4*c^(1/2))^-1)+2*(-1)^3//4*b^2*c^(1/2)*arctan(x*(-1)^3//4*c^(1/2))*log(2*(1+x*(-1)^1//4*c^(1/2))^-1)+2*(-1)^3//4*b^2*c^(1/2)*arctanh(x*(-1)^3//4*c^(1/2))*log(2*(1+-1*x*(-1)^3//4*c^(1/2))^-1))
@test integrate(x^-4*(a+b*arctan(c*x^2))^2, x) == :(-1//12*x^-3*(2a+im*b*log(1+-1*im*c*x^2))^2+1//12*b^2*x^-3*log(1+im*c*x^2)^2+-4//3*(-1)^1//4*b^2*c^3//2*arctan(x*(-1)^3//4*c^(1/2))+-4//3*(-1)^1//4*b^2*c^3//2*arctanh(x*(-1)^3//4*c^(1/2))+-2//3*a*b*c*x^-1+-1//3*b*c*x^-1*(2a+im*b*log(1+-1*im*c*x^2))+-1//3*(-1)^1//4*b^2*c^3//2*arctanh(x*(-1)^3//4*c^(1/2))^2+-1//6*(-1)^1//4*b^2*c^3//2*Polylogarithms.polylog(2,1+2^(1/2)*(1+x*(-1)^3//4*c^(1/2))^-1*((-1)^3//4+x*c^(1/2)))+-1//6*(-1)^1//4*b^2*c^3//2*Polylogarithms.polylog(2,1+-1*(1+x*(-1)^3//4*c^(1/2))^-1*(1+im)*(1+x*(-1)^1//4*c^(1/2)))+-1//6*(-1)^3//4*b^2*c^3//2*Polylogarithms.polylog(2,1+-1*2^(1/2)*(1+x*(-1)^1//4*c^(1/2))^-1*((-1)^1//4+x*c^(1/2)))+-1//6*(-1)^3//4*b^2*c^3//2*Polylogarithms.polylog(2,1+-1*(1+x*(-1)^1//4*c^(1/2))^-1*(1+-1im)*(1+x*(-1)^3//4*c^(1/2)))+-1//6*b^2*x^-3*log(1+im*c*x^2)*log(1+-1*im*c*x^2)+1//3*(-1)^1//4*b^2*c^3//2*Polylogarithms.polylog(2,1+-2*(1+x*(-1)^3//4*c^(1/2))^-1)+1//3*(-1)^1//4*b^2*c^3//2*Polylogarithms.polylog(2,1+-2*(1+-1*x*(-1)^3//4*c^(1/2))^-1)+1//3*(-1)^3//4*b^2*c^3//2*arctan(x*(-1)^3//4*c^(1/2))^2+1//3*(-1)^3//4*b^2*c^3//2*Polylogarithms.polylog(2,1+-2*(1+x*(-1)^1//4*c^(1/2))^-1)+1//3*(-1)^3//4*b^2*c^3//2*Polylogarithms.polylog(2,1+-2*(1+-1*x*(-1)^1//4*c^(1/2))^-1)+-2//3*(-1)^1//4*b^2*c^3//2*arctan(x*(-1)^3//4*c^(1/2))*log(2*(1+x*(-1)^1//4*c^(1/2))^-1)+-2//3*(-1)^1//4*b^2*c^3//2*arctanh(x*(-1)^3//4*c^(1/2))*log(2*(1+x*(-1)^3//4*c^(1/2))^-1)+-1//3*im*c*b^2*x^-1*log(1+-1*im*c*x^2)+-1//3*b*(-1)^3//4*c^3//2*(2a+im*b*log(1+-1*im*c*x^2))*arctan(x*(-1)^3//4*c^(1/2))+-1//3*(-1)^1//4*b^2*c^3//2*arctan(x*(-1)^3//4*c^(1/2))*log(1+im*c*x^2)+-1//3*(-1)^1//4*b^2*c^3//2*arctanh(x*(-1)^3//4*c^(1/2))*log(1+-1*im*c*x^2)+1//3*im*a*b*x^-3*log(1+im*c*x^2)+1//3*(-1)^1//4*b^2*c^3//2*arctan(x*(-1)^3//4*c^(1/2))*log(2^(1/2)*(1+x*(-1)^1//4*c^(1/2))^-1*((-1)^1//4+x*c^(1/2)))+1//3*(-1)^1//4*b^2*c^3//2*arctan(x*(-1)^3//4*c^(1/2))*log((1+x*(-1)^1//4*c^(1/2))^-1*(1+-1im)*(1+x*(-1)^3//4*c^(1/2)))+1//3*(-1)^1//4*b^2*c^3//2*arctanh(x*(-1)^3//4*c^(1/2))*log((1+x*(-1)^3//4*c^(1/2))^-1*(1+im)*(1+x*(-1)^1//4*c^(1/2)))+1//3*(-1)^1//4*b^2*c^3//2*arctanh(x*(-1)^3//4*c^(1/2))*log(-1*2^(1/2)*(1+x*(-1)^3//4*c^(1/2))^-1*((-1)^3//4+x*c^(1/2)))+1//3*(-1)^1//4*b^2*c^3//2*arctanh(x*(-1)^3//4*c^(1/2))*log(1+im*c*x^2)+2//3*im*c*b^2*x^-1*log(1+im*c*x^2)+2//3*a*b*(-1)^3//4*c^3//2*arctanh(x*(-1)^3//4*c^(1/2))+2//3*(-1)^1//4*b^2*c^3//2*arctan(x*(-1)^3//4*c^(1/2))*log(2*(1+-1*x*(-1)^1//4*c^(1/2))^-1)+2//3*(-1)^1//4*b^2*c^3//2*arctanh(x*(-1)^3//4*c^(1/2))*log(2*(1+-1*x*(-1)^3//4*c^(1/2))^-1))
@test integrate(x^-6*(a+b*arctan(c*x^2))^2, x) == :(-1//20*x^-5*(2a+im*b*log(1+-1*im*c*x^2))^2+-8//15*b^2*c^2*x^-1+1//20*b^2*x^-5*log(1+im*c*x^2)^2+-4//15*(-1)^3//4*b^2*c^5//2*arctan(x*(-1)^3//4*c^(1/2))+-2//15*a*b*c*x^-3+-1//5*(-1)^1//4*b^2*c^5//2*arctan(x*(-1)^3//4*c^(1/2))^2+-1//5*(-1)^1//4*b^2*c^5//2*Polylogarithms.polylog(2,1+-2*(1+x*(-1)^1//4*c^(1/2))^-1)+-1//5*(-1)^1//4*b^2*c^5//2*Polylogarithms.polylog(2,1+-2*(1+-1*x*(-1)^1//4*c^(1/2))^-1)+-1//5*(-1)^3//4*b^2*c^5//2*Polylogarithms.polylog(2,1+-2*(1+x*(-1)^3//4*c^(1/2))^-1)+-1//5*(-1)^3//4*b^2*c^5//2*Polylogarithms.polylog(2,1+-2*(1+-1*x*(-1)^3//4*c^(1/2))^-1)+-1//5*b^2*c^2*x^-1*log(1+-1*im*c*x^2)+-1//10*b^2*x^-5*log(1+im*c*x^2)*log(1+-1*im*c*x^2)+-1//15*b*c*x^-3*(2a+im*b*log(1+-1*im*c*x^2))+1//5*(-1)^3//4*b^2*c^5//2*arctanh(x*(-1)^3//4*c^(1/2))^2+1//10*(-1)^1//4*b^2*c^5//2*Polylogarithms.polylog(2,1+-1*2^(1/2)*(1+x*(-1)^1//4*c^(1/2))^-1*((-1)^1//4+x*c^(1/2)))+1//10*(-1)^1//4*b^2*c^5//2*Polylogarithms.polylog(2,1+-1*(1+x*(-1)^1//4*c^(1/2))^-1*(1+-1im)*(1+x*(-1)^3//4*c^(1/2)))+1//10*(-1)^3//4*b^2*c^5//2*Polylogarithms.polylog(2,1+2^(1/2)*(1+x*(-1)^3//4*c^(1/2))^-1*((-1)^3//4+x*c^(1/2)))+1//10*(-1)^3//4*b^2*c^5//2*Polylogarithms.polylog(2,1+-1*(1+x*(-1)^3//4*c^(1/2))^-1*(1+im)*(1+x*(-1)^1//4*c^(1/2)))+4//15*(-1)^3//4*b^2*c^5//2*arctanh(x*(-1)^3//4*c^(1/2))+-2//5*(-1)^3//4*b^2*c^5//2*arctan(x*(-1)^3//4*c^(1/2))*log(2*(1+x*(-1)^1//4*c^(1/2))^-1)+-2//5*(-1)^3//4*b^2*c^5//2*arctanh(x*(-1)^3//4*c^(1/2))*log(2*(1+-1*x*(-1)^3//4*c^(1/2))^-1)+-1//5*im*b*c^2*x^-1*(2a+im*b*log(1+-1*im*c*x^2))+-1//5*(-1)^3//4*b^2*c^5//2*arctan(x*(-1)^3//4*c^(1/2))*log(1+im*c*x^2)+-1//5*(-1)^3//4*b^2*c^5//2*arctanh(x*(-1)^3//4*c^(1/2))*log((1+x*(-1)^3//4*c^(1/2))^-1*(1+im)*(1+x*(-1)^1//4*c^(1/2)))+-1//5*(-1)^3//4*b^2*c^5//2*arctanh(x*(-1)^3//4*c^(1/2))*log(-1*2^(1/2)*(1+x*(-1)^3//4*c^(1/2))^-1*((-1)^3//4+x*c^(1/2)))+-1//5*(-1)^3//4*b^2*c^5//2*arctanh(x*(-1)^3//4*c^(1/2))*log(1+im*c*x^2)+-1//15*im*c*b^2*x^-3*log(1+-1*im*c*x^2)+1//5*im*a*b*x^-5*log(1+im*c*x^2)+1//5*b*(-1)^1//4*c^5//2*(2a+im*b*log(1+-1*im*c*x^2))*arctan(x*(-1)^3//4*c^(1/2))+1//5*(-1)^3//4*b^2*c^5//2*arctan(x*(-1)^3//4*c^(1/2))*log(2^(1/2)*(1+x*(-1)^1//4*c^(1/2))^-1*((-1)^1//4+x*c^(1/2)))+1//5*(-1)^3//4*b^2*c^5//2*arctan(x*(-1)^3//4*c^(1/2))*log((1+x*(-1)^1//4*c^(1/2))^-1*(1+-1im)*(1+x*(-1)^3//4*c^(1/2)))+1//5*(-1)^3//4*b^2*c^5//2*arctanh(x*(-1)^3//4*c^(1/2))*log(1+-1*im*c*x^2)+2//5*im*a*b*c^2*x^-1+2//5*a*b*(-1)^1//4*c^5//2*arctanh(x*(-1)^3//4*c^(1/2))+2//5*(-1)^3//4*b^2*c^5//2*arctan(x*(-1)^3//4*c^(1/2))*log(2*(1+-1*x*(-1)^1//4*c^(1/2))^-1)+2//5*(-1)^3//4*b^2*c^5//2*arctanh(x*(-1)^3//4*c^(1/2))*log(2*(1+x*(-1)^3//4*c^(1/2))^-1)+2//15*im*c*b^2*x^-3*log(1+im*c*x^2))
@test integrate(x^3*(a+b*arctan(c*x^2))^3, x) == :(1//4*c^-2*(a+b*arctan(c*x^2))^3+1//4*x^4*(a+b*arctan(c*x^2))^3+-3//2*b^2*c^-2*(a+b*arctan(c*x^2))*log(2*(1+im*c*x^2)^-1)+-3//4*im*b*c^-2*(a+b*arctan(c*x^2))^2+-3//4*im*b^3*c^-2*Polylogarithms.polylog(2,1+-2*(1+im*c*x^2)^-1)+-3//4*b*c^-1*x^2*(a+b*arctan(c*x^2))^2)
@test integrate(x*(a+b*arctan(c*x^2))^3, x) == :((1/2)*x^2*(a+b*arctan(c*x^2))^3+(1/2)*im*c^-1*(a+b*arctan(c*x^2))^3+3//4*b^3*c^-1*Polylogarithms.polylog(3,1+-2*(1+im*c*x^2)^-1)+3//2*b*c^-1*(a+b*arctan(c*x^2))^2*log(2*(1+im*c*x^2)^-1)+3//2*im*b^2*c^-1*(a+b*arctan(c*x^2))*Polylogarithms.polylog(2,1+-2*(1+im*c*x^2)^-1))
@test integrate(x^-1*(a+b*arctan(c*x^2))^3, x) == :((a+b*arctan(c*x^2))^3*arctanh(1+-2*(1+im*c*x^2)^-1)+-3//4*b^2*(a+b*arctan(c*x^2))*Polylogarithms.polylog(3,1+-2*(1+im*c*x^2)^-1)+-3//8*im*b^3*Polylogarithms.polylog(4,-1+2*(1+im*c*x^2)^-1)+3//4*b^2*(a+b*arctan(c*x^2))*Polylogarithms.polylog(3,-1+2*(1+im*c*x^2)^-1)+3//8*im*b^3*Polylogarithms.polylog(4,1+-2*(1+im*c*x^2)^-1)+-3//4*im*b*(a+b*arctan(c*x^2))^2*Polylogarithms.polylog(2,1+-2*(1+im*c*x^2)^-1)+3//4*im*b*(a+b*arctan(c*x^2))^2*Polylogarithms.polylog(2,-1+2*(1+im*c*x^2)^-1))
@test integrate(x^-3*(a+b*arctan(c*x^2))^3, x) == :(-1//2*x^-2*(a+b*arctan(c*x^2))^3+-1//2*im*c*(a+b*arctan(c*x^2))^3+3//4*c*b^3*Polylogarithms.polylog(3,-1+2*(1+-1*im*c*x^2)^-1)+3//2*b*c*(a+b*arctan(c*x^2))^2*log(2+-2*(1+-1*im*c*x^2)^-1)+-3//2*im*c*b^2*(a+b*arctan(c*x^2))*Polylogarithms.polylog(2,-1+2*(1+-1*im*c*x^2)^-1))
@test integrate(x^-5*(a+b*arctan(c*x^2))^3, x) == :(-1//4*c^2*(a+b*arctan(c*x^2))^3+-1//4*x^-4*(a+b*arctan(c*x^2))^3+-3//4*im*b*c^2*(a+b*arctan(c*x^2))^2+-3//4*im*b^3*c^2*Polylogarithms.polylog(2,-1+2*(1+-1*im*c*x^2)^-1)+-3//4*b*c*x^-2*(a+b*arctan(c*x^2))^2+3//2*b^2*c^2*(a+b*arctan(c*x^2))*log(2+-2*(1+-1*im*c*x^2)^-1))
@test integrate(x^11*(a+b*arctan(c*x^3)), x) == :(1//12*x^12*(a+b*arctan(c*x^3))+-1//12*b*c^-4*arctan(c*x^3)+-1//36*b*c^-1*x^9+1//12*b*c^-3*x^3)
@test integrate(x^8*(a+b*arctan(c*x^3)), x) == :(1//9*x^9*(a+b*arctan(c*x^3))+-1//18*b*c^-1*x^6+1//18*b*c^-3*log(1+c^2*x^6))
@test integrate(x^5*(a+b*arctan(c*x^3)), x) == :(1//6*x^6*(a+b*arctan(c*x^3))+-1//6*b*c^-1*x^3+1//6*b*c^-2*arctan(c*x^3))
@test integrate(x^2*(a+b*arctan(c*x^3)), x) == :(1//3*x^3*(a+b*arctan(c*x^3))+-1//6*b*c^-1*log(1+c^2*x^6))
@test integrate(x^-1*(a+b*arctan(c*x^3)), x) == :(a*log(x)+-1//6*im*b*Polylogarithms.polylog(2,im*c*x^3)+1//6*im*b*Polylogarithms.polylog(2,-1*im*c*x^3))
@test integrate(x^-4*(a+b*arctan(c*x^3)), x) == :(-1//3*x^-3*(a+b*arctan(c*x^3))+b*c*log(x)+-1//6*b*c*log(1+c^2*x^6))
@test integrate(x^-7*(a+b*arctan(c*x^3)), x) == :(-1//6*x^-6*(a+b*arctan(c*x^3))+-1//6*b*c*x^-3+-1//6*b*c^2*arctan(c*x^3))
@test integrate(x^-10*(a+b*arctan(c*x^3)), x) == :(-1//9*x^-9*(a+b*arctan(c*x^3))+-1//3*b*c^3*log(x)+-1//18*b*c*x^-6+1//18*b*c^3*log(1+c^2*x^6))
@test integrate(x^3*(a+b*arctan(c*x^3)), x) == :(1//4*x^4*(a+b*arctan(c*x^3))+-3//4*b*x*c^-1+-1//8*b*c^-4//3*arctan(3^(1/2)+-2*x*c^1//3)+1//4*b*c^-4//3*arctan(x*c^1//3)+1//8*b*c^-4//3*arctan(3^(1/2)+2*x*c^1//3)+-1//16*b*3^(1/2)*c^-4//3*log(1+c^2//3*x^2+-1*x*3^(1/2)*c^1//3)+1//16*b*3^(1/2)*c^-4//3*log(1+c^2//3*x^2+x*3^(1/2)*c^1//3))
@test integrate(a+b*arctan(c*x^3), x) == :(a*x+b*x*arctan(c*x^3)+(1/2)*b*c^-1//3*log(1+c^2//3*x^2)+-1//4*b*c^-1//3*log(1+c^4//3*x^4+-1*c^2//3*x^2)+(1/2)*b*3^(1/2)*c^-1//3*arctan(1//3*3^(1/2)*(1+-2*c^2//3*x^2)))
@test integrate(x^-3*(a+b*arctan(c*x^3)), x) == :(-1//2*x^-2*(a+b*arctan(c*x^3))+(1/2)*b*c^2//3*arctan(x*c^1//3)+-1//4*b*c^2//3*arctan(3^(1/2)+-2*x*c^1//3)+1//4*b*c^2//3*arctan(3^(1/2)+2*x*c^1//3)+-1//8*b*3^(1/2)*c^2//3*log(1+c^2//3*x^2+-1*x*3^(1/2)*c^1//3)+1//8*b*3^(1/2)*c^2//3*log(1+c^2//3*x^2+x*3^(1/2)*c^1//3))
@test integrate(x^-6*(a+b*arctan(c*x^3)), x) == :(-1//5*x^-5*(a+b*arctan(c*x^3))+-3//10*b*c*x^-2+-1//20*b*c^5//3*log(1+c^4//3*x^4+-1*c^2//3*x^2)+1//10*b*c^5//3*log(1+c^2//3*x^2)+1//10*b*3^(1/2)*c^5//3*arctan(1//3*3^(1/2)*(1+-2*c^2//3*x^2)))
@test integrate(x^7*(a+b*arctan(c*x^3)), x) == :(1//8*x^8*(a+b*arctan(c*x^3))+-3//40*b*c^-1*x^5+-1//16*b*c^-8//3*arctan(3^(1/2)+-2*x*c^1//3)+1//8*b*c^-8//3*arctan(x*c^1//3)+1//16*b*c^-8//3*arctan(3^(1/2)+2*x*c^1//3)+-1//32*b*3^(1/2)*c^-8//3*log(1+c^2//3*x^2+x*3^(1/2)*c^1//3)+1//32*b*3^(1/2)*c^-8//3*log(1+c^2//3*x^2+-1*x*3^(1/2)*c^1//3))
@test integrate(x^4*(a+b*arctan(c*x^3)), x) == :(1//5*x^5*(a+b*arctan(c*x^3))+-3//10*b*c^-1*x^2+-1//20*b*c^-5//3*log(1+c^4//3*x^4+-1*c^2//3*x^2)+1//10*b*c^-5//3*log(1+c^2//3*x^2)+-1//10*b*3^(1/2)*c^-5//3*arctan(1//3*3^(1/2)*(1+-2*c^2//3*x^2)))
@test integrate(x*(a+b*arctan(c*x^3)), x) == :((1/2)*x^2*(a+b*arctan(c*x^3))+-1//2*b*c^-2//3*arctan(x*c^1//3)+-1//4*b*c^-2//3*arctan(3^(1/2)+2*x*c^1//3)+1//4*b*c^-2//3*arctan(3^(1/2)+-2*x*c^1//3)+-1//8*b*3^(1/2)*c^-2//3*log(1+c^2//3*x^2+-1*x*3^(1/2)*c^1//3)+1//8*b*3^(1/2)*c^-2//3*log(1+c^2//3*x^2+x*3^(1/2)*c^1//3))
@test integrate(x^-2*(a+b*arctan(c*x^3)), x) == :(-1*x^-1*(a+b*arctan(c*x^3))+(1/2)*b*c^1//3*log(1+c^2//3*x^2)+-1//4*b*c^1//3*log(1+c^4//3*x^4+-1*c^2//3*x^2)+-1//2*b*3^(1/2)*c^1//3*arctan(1//3*3^(1/2)*(1+-2*c^2//3*x^2)))
@test integrate(x^-5*(a+b*arctan(c*x^3)), x) == :(-1//4*x^-4*(a+b*arctan(c*x^3))+-3//4*b*c*x^-1+-1//4*b*c^4//3*arctan(x*c^1//3)+-1//8*b*c^4//3*arctan(3^(1/2)+2*x*c^1//3)+1//8*b*c^4//3*arctan(3^(1/2)+-2*x*c^1//3)+-1//16*b*3^(1/2)*c^4//3*log(1+c^2//3*x^2+-1*x*3^(1/2)*c^1//3)+1//16*b*3^(1/2)*c^4//3*log(1+c^2//3*x^2+x*3^(1/2)*c^1//3))
@test integrate(x^11*(a+b*arctan(c*x^3))^2, x) == :(-1//12*c^-4*(a+b*arctan(c*x^3))^2+1//12*x^12*(a+b*arctan(c*x^3))^2+-1//9*b^2*c^-4*log(1+c^2*x^6)+1//36*b^2*c^-2*x^6+-1//18*b*c^-1*x^9*(a+b*arctan(c*x^3))+1//6*a*b*c^-3*x^3+1//6*b^2*c^-3*x^3*arctan(c*x^3))
@test integrate(x^8*(a+b*arctan(c*x^3))^2, x) == :(1//9*x^9*(a+b*arctan(c*x^3))^2+-1//9*im*c^-3*(a+b*arctan(c*x^3))^2+-1//9*b^2*c^-3*arctan(c*x^3)+1//9*b^2*c^-2*x^3+-2//9*b*c^-3*(a+b*arctan(c*x^3))*log(2*(1+im*c*x^3)^-1)+-1//9*im*b^2*c^-3*Polylogarithms.polylog(2,1+-2*(1+im*c*x^3)^-1)+-1//9*b*c^-1*x^6*(a+b*arctan(c*x^3)))
@test integrate(x^5*(a+b*arctan(c*x^3))^2, x) == :(1//6*c^-2*(a+b*arctan(c*x^3))^2+1//6*x^6*(a+b*arctan(c*x^3))^2+1//6*b^2*c^-2*log(1+c^2*x^6)+-1//3*a*b*c^-1*x^3+-1//3*b^2*c^-1*x^3*arctan(c*x^3))
@test integrate(x^2*(a+b*arctan(c*x^3))^2, x) == :(1//3*x^3*(a+b*arctan(c*x^3))^2+1//3*im*c^-1*(a+b*arctan(c*x^3))^2+1//3*im*b^2*c^-1*Polylogarithms.polylog(2,1+-2*(1+im*c*x^3)^-1)+2//3*b*c^-1*(a+b*arctan(c*x^3))*log(2*(1+im*c*x^3)^-1))
@test integrate(x^-1*(a+b*arctan(c*x^3))^2, x) == :(-1//6*b^2*Polylogarithms.polylog(3,1+-2*(1+im*c*x^3)^-1)+1//6*b^2*Polylogarithms.polylog(3,-1+2*(1+im*c*x^3)^-1)+2//3*(a+b*arctan(c*x^3))^2*arctanh(1+-2*(1+im*c*x^3)^-1)+-1//3*im*b*(a+b*arctan(c*x^3))*Polylogarithms.polylog(2,1+-2*(1+im*c*x^3)^-1)+1//3*im*b*(a+b*arctan(c*x^3))*Polylogarithms.polylog(2,-1+2*(1+im*c*x^3)^-1))
@test integrate(x^-4*(a+b*arctan(c*x^3))^2, x) == :(-1//3*x^-3*(a+b*arctan(c*x^3))^2+-1//3*im*c*(a+b*arctan(c*x^3))^2+-1//3*im*c*b^2*Polylogarithms.polylog(2,-1+2*(1+-1*im*c*x^3)^-1)+2//3*b*c*(a+b*arctan(c*x^3))*log(2+-2*(1+-1*im*c*x^3)^-1))
@test integrate(x^-7*(a+b*arctan(c*x^3))^2, x) == :(-1//6*c^2*(a+b*arctan(c*x^3))^2+-1//6*x^-6*(a+b*arctan(c*x^3))^2+b^2*c^2*log(x)+-1//6*b^2*c^2*log(1+c^2*x^6)+-1//3*b*c*x^-3*(a+b*arctan(c*x^3)))
@test integrate(x^-10*(a+b*arctan(c*x^3))^2, x) == :(-1//9*x^-9*(a+b*arctan(c*x^3))^2+-1//9*b^2*c^2*x^-3+-1//9*b^2*c^3*arctan(c*x^3)+1//9*im*c^3*(a+b*arctan(c*x^3))^2+-2//9*b*c^3*(a+b*arctan(c*x^3))*log(2+-2*(1+-1*im*c*x^3)^-1)+-1//9*b*c*x^-6*(a+b*arctan(c*x^3))+1//9*im*b^2*c^3*Polylogarithms.polylog(2,-1+2*(1+-1*im*c*x^3)^-1))
@test integrate(x^8*(a+b*arctan(c*x^3))^3, x) == :(1//9*x^9*(a+b*arctan(c*x^3))^3+-1//6*b*c^-3*(a+b*arctan(c*x^3))^2+-1//6*b^3*c^-3*Polylogarithms.polylog(3,1+-2*(1+im*c*x^3)^-1)+-1//6*b^3*c^-3*log(1+c^2*x^6)+-1//9*im*c^-3*(a+b*arctan(c*x^3))^3+-1//3*b*c^-3*(a+b*arctan(c*x^3))^2*log(2*(1+im*c*x^3)^-1)+-1//6*b*c^-1*x^6*(a+b*arctan(c*x^3))^2+1//3*a*b^2*c^-2*x^3+1//3*b^3*c^-2*x^3*arctan(c*x^3)+-1//3*im*b^2*c^-3*(a+b*arctan(c*x^3))*Polylogarithms.polylog(2,1+-2*(1+im*c*x^3)^-1))
@test integrate(x^5*(a+b*arctan(c*x^3))^3, x) == :(1//6*c^-2*(a+b*arctan(c*x^3))^3+1//6*x^6*(a+b*arctan(c*x^3))^3+-1*b^2*c^-2*(a+b*arctan(c*x^3))*log(2*(1+im*c*x^3)^-1)+-1//2*im*b*c^-2*(a+b*arctan(c*x^3))^2+-1//2*im*b^3*c^-2*Polylogarithms.polylog(2,1+-2*(1+im*c*x^3)^-1)+-1//2*b*c^-1*x^3*(a+b*arctan(c*x^3))^2)
@test integrate(x^2*(a+b*arctan(c*x^3))^3, x) == :(1//3*x^3*(a+b*arctan(c*x^3))^3+(1/2)*b^3*c^-1*Polylogarithms.polylog(3,1+-2*(1+im*c*x^3)^-1)+1//3*im*c^-1*(a+b*arctan(c*x^3))^3+b*c^-1*(a+b*arctan(c*x^3))^2*log(2*(1+im*c*x^3)^-1)+im*b^2*c^-1*(a+b*arctan(c*x^3))*Polylogarithms.polylog(2,1+-2*(1+im*c*x^3)^-1))
@test integrate(x^-1*(a+b*arctan(c*x^3))^3, x) == :(2//3*(a+b*arctan(c*x^3))^3*arctanh(1+-2*(1+im*c*x^3)^-1)+(1/2)*b^2*(a+b*arctan(c*x^3))*Polylogarithms.polylog(3,-1+2*(1+im*c*x^3)^-1)+-1//2*b^2*(a+b*arctan(c*x^3))*Polylogarithms.polylog(3,1+-2*(1+im*c*x^3)^-1)+-1//4*im*b^3*Polylogarithms.polylog(4,-1+2*(1+im*c*x^3)^-1)+1//4*im*b^3*Polylogarithms.polylog(4,1+-2*(1+im*c*x^3)^-1)+(1/2)*im*b*(a+b*arctan(c*x^3))^2*Polylogarithms.polylog(2,-1+2*(1+im*c*x^3)^-1)+-1//2*im*b*(a+b*arctan(c*x^3))^2*Polylogarithms.polylog(2,1+-2*(1+im*c*x^3)^-1))
@test integrate(x^-4*(a+b*arctan(c*x^3))^3, x) == :(-1//3*x^-3*(a+b*arctan(c*x^3))^3+(1/2)*c*b^3*Polylogarithms.polylog(3,-1+2*(1+-1*im*c*x^3)^-1)+-1//3*im*c*(a+b*arctan(c*x^3))^3+b*c*(a+b*arctan(c*x^3))^2*log(2+-2*(1+-1*im*c*x^3)^-1)+-1*im*c*b^2*(a+b*arctan(c*x^3))*Polylogarithms.polylog(2,-1+2*(1+-1*im*c*x^3)^-1))
@test integrate(x^-7*(a+b*arctan(c*x^3))^3, x) == :(-1//6*c^2*(a+b*arctan(c*x^3))^3+-1//6*x^-6*(a+b*arctan(c*x^3))^3+b^2*c^2*(a+b*arctan(c*x^3))*log(2+-2*(1+-1*im*c*x^3)^-1)+-1//2*im*b*c^2*(a+b*arctan(c*x^3))^2+-1//2*im*b^3*c^2*Polylogarithms.polylog(2,-1+2*(1+-1*im*c*x^3)^-1)+-1//2*b*c*x^-3*(a+b*arctan(c*x^3))^2)
@test integrate(x^3*(a+b*arctan(c*x^-1)), x) == :(1//4*x^4*(a+b*arctan(c*x^-1))+-1//4*b*x*c^3+1//4*b*c^4*arctan(x*c^-1)+1//12*b*c*x^3)
@test integrate(x^2*(a+b*arctan(c*x^-1)), x) == :(1//3*x^3*(a+b*arctan(c*x^-1))+-1//6*b*c^3*log(c^2+x^2)+1//6*b*c*x^2)
@test integrate(x*(a+b*arctan(c*x^-1)), x) == :((1/2)*x^2*(a+b*arctan(c*x^-1))+(1/2)*b*c*x+-1//2*b*c^2*arctan(x*c^-1))
@test integrate(a+b*arctan(c*x^-1), x) == :(a*x+b*x*arctan(c*x^-1)+(1/2)*b*c*log(c^2+x^2))
@test integrate(x^-1*(a+b*arctan(c*x^-1)), x) == :(a*log(x)+(1/2)*im*b*Polylogarithms.polylog(2,im*c*x^-1)+-1//2*im*b*Polylogarithms.polylog(2,-1*im*c*x^-1))
@test integrate(x^-2*(a+b*arctan(c*x^-1)), x) == :(-1*x^-1*(a+b*arctan(c*x^-1))+(1/2)*b*c^-1*log(1+c^2*x^-2))
@test integrate(x^-3*(a+b*arctan(c*x^-1)), x) == :(-1//2*x^-2*(a+b*arctan(c*x^-1))+(1/2)*b*c^-1*x^-1+(1/2)*b*c^-2*arctan(x*c^-1))
@test integrate(x^-4*(a+b*arctan(c*x^-1)), x) == :(-1//3*x^-3*(a+b*arctan(c*x^-1))+-1//6*b*c^-3*log(c^2+x^2)+1//3*b*c^-3*log(x)+1//6*b*c^-1*x^-2)
@test integrate(x^3*(a+b*arctan(c*x^-1))^2, x) == :(-1//4*c^4*(a+b*arccot(x*c^-1))^2+1//4*x^4*(a+b*arccot(x*c^-1))^2+-2//3*b^2*c^4*log(x)+-1//3*b^2*c^4*log(1+c^2*x^-2)+1//12*b^2*c^2*x^2+-1//2*b*x*c^3*(a+b*arccot(x*c^-1))+1//6*b*c*x^3*(a+b*arccot(x*c^-1)))
@test integrate(x^2*(a+b*arctan(c*x^-1))^2, x) == :(1//3*x^3*(a+b*arccot(x*c^-1))^2+-1//3*im*c^3*(a+b*arccot(x*c^-1))^2+1//3*x*b^2*c^2+1//3*b^2*c^3*arccot(x*c^-1)+-1//3*im*b^2*c^3*Polylogarithms.polylog(2,-1+2*(1+-1*im*c*x^-1)^-1)+1//3*b*c*x^2*(a+b*arccot(x*c^-1))+2//3*b*c^3*(a+b*arccot(x*c^-1))*log(2+-2*(1+-1*im*c*x^-1)^-1))
@test integrate(x*(a+b*arctan(c*x^-1))^2, x) == :((1/2)*c^2*(a+b*arccot(x*c^-1))^2+(1/2)*x^2*(a+b*arccot(x*c^-1))^2+b^2*c^2*log(x)+(1/2)*b^2*c^2*log(1+c^2*x^-2)+b*c*x*(a+b*arccot(x*c^-1)))
@test integrate((a+b*arctan(c*x^-1))^2, x) == :(x*(a+b*arccot(x*c^-1))^2+im*c*(a+b*arccot(x*c^-1))^2+im*c*b^2*Polylogarithms.polylog(2,1+-2*c*(c+im*x)^-1)+-2*b*c*(a+b*arccot(x*c^-1))*log(2*c*(c+im*x)^-1))
@test integrate(x^-1*(a+b*arctan(c*x^-1))^2, x) == :((1/2)*b^2*Polylogarithms.polylog(3,1+-2*(1+im*c*x^-1)^-1)+-2*(a+b*arccot(x*c^-1))^2*arctanh(1+-2*(1+im*c*x^-1)^-1)+-1//2*b^2*Polylogarithms.polylog(3,-1+2*(1+im*c*x^-1)^-1)+im*b*(a+b*arccot(x*c^-1))*Polylogarithms.polylog(2,1+-2*(1+im*c*x^-1)^-1)+-1*im*b*(a+b*arccot(x*c^-1))*Polylogarithms.polylog(2,-1+2*(1+im*c*x^-1)^-1))
@test integrate(x^-2*(a+b*arctan(c*x^-1))^2, x) == :(-1*x^-1*(a+b*arccot(x*c^-1))^2+-1*im*c^-1*(a+b*arccot(x*c^-1))^2+-1*im*b^2*c^-1*Polylogarithms.polylog(2,1+-2*(1+im*c*x^-1)^-1)+-2*b*c^-1*(a+b*arccot(x*c^-1))*log(2*(1+im*c*x^-1)^-1))
@test integrate(x^-3*(a+b*arctan(c*x^-1))^2, x) == :(-1//2*c^-2*(a+b*arccot(x*c^-1))^2+-1//2*x^-2*(a+b*arccot(x*c^-1))^2+-1//2*b^2*c^-2*log(1+c^2*x^-2)+a*b*c^-1*x^-1+b^2*c^-1*x^-1*arccot(x*c^-1))
@test integrate(x^3*(a+b*arctan(c*x^-1))^3, x) == :(-1//4*c^4*(a+b*arccot(x*c^-1))^3+1//4*x^4*(a+b*arccot(x*c^-1))^3+1//4*x*b^3*c^3+1//4*b^3*c^4*arccot(x*c^-1)+-1*im*b*c^4*(a+b*arccot(x*c^-1))^2+-1*im*b^3*c^4*Polylogarithms.polylog(2,-1+2*(1+-1*im*c*x^-1)^-1)+2*b^2*c^4*(a+b*arccot(x*c^-1))*log(2+-2*(1+-1*im*c*x^-1)^-1)+-3//4*b*x*c^3*(a+b*arccot(x*c^-1))^2+1//4*b*c*x^3*(a+b*arccot(x*c^-1))^2+1//4*b^2*c^2*x^2*(a+b*arccot(x*c^-1)))
@test integrate(x^2*(a+b*arctan(c*x^-1))^3, x) == :(1//3*x^3*(a+b*arccot(x*c^-1))^3+b^3*c^3*log(x)+(1/2)*b*c^3*(a+b*arccot(x*c^-1))^2+(1/2)*b^3*c^3*Polylogarithms.polylog(3,-1+2*(1+-1*im*c*x^-1)^-1)+(1/2)*b^3*c^3*log(1+c^2*x^-2)+-1//3*im*c^3*(a+b*arccot(x*c^-1))^3+b*c^3*(a+b*arccot(x*c^-1))^2*log(2+-2*(1+-1*im*c*x^-1)^-1)+x*b^2*c^2*(a+b*arccot(x*c^-1))+(1/2)*b*c*x^2*(a+b*arccot(x*c^-1))^2+-1*im*b^2*c^3*(a+b*arccot(x*c^-1))*Polylogarithms.polylog(2,-1+2*(1+-1*im*c*x^-1)^-1))
@test integrate(x*(a+b*arctan(c*x^-1))^3, x) == :((1/2)*c^2*(a+b*arccot(x*c^-1))^3+(1/2)*x^2*(a+b*arccot(x*c^-1))^3+-3*b^2*c^2*(a+b*arccot(x*c^-1))*log(2+-2*(1+-1*im*c*x^-1)^-1)+3//2*im*b*c^2*(a+b*arccot(x*c^-1))^2+3//2*im*b^3*c^2*Polylogarithms.polylog(2,-1+2*(1+-1*im*c*x^-1)^-1)+3//2*b*c*x*(a+b*arccot(x*c^-1))^2)
@test integrate((a+b*arctan(c*x^-1))^3, x) == :(x*(a+b*arccot(x*c^-1))^3+im*c*(a+b*arccot(x*c^-1))^3+-3//2*c*b^3*Polylogarithms.polylog(3,1+-2*c*(c+im*x)^-1)+-3*b*c*(a+b*arccot(x*c^-1))^2*log(2*c*(c+im*x)^-1)+3*im*c*b^2*(a+b*arccot(x*c^-1))*Polylogarithms.polylog(2,1+-2*c*(c+im*x)^-1))
@test integrate(x^-1*(a+b*arctan(c*x^-1))^3, x) == :(-2*(a+b*arccot(x*c^-1))^3*arctanh(1+-2*(1+im*c*x^-1)^-1)+-3//2*b^2*(a+b*arccot(x*c^-1))*Polylogarithms.polylog(3,-1+2*(1+im*c*x^-1)^-1)+-3//4*im*b^3*Polylogarithms.polylog(4,1+-2*(1+im*c*x^-1)^-1)+3//2*b^2*(a+b*arccot(x*c^-1))*Polylogarithms.polylog(3,1+-2*(1+im*c*x^-1)^-1)+3//4*im*b^3*Polylogarithms.polylog(4,-1+2*(1+im*c*x^-1)^-1)+-3//2*im*b*(a+b*arccot(x*c^-1))^2*Polylogarithms.polylog(2,-1+2*(1+im*c*x^-1)^-1)+3//2*im*b*(a+b*arccot(x*c^-1))^2*Polylogarithms.polylog(2,1+-2*(1+im*c*x^-1)^-1))
@test integrate(x^-2*(a+b*arctan(c*x^-1))^3, x) == :(-1*x^-1*(a+b*arccot(x*c^-1))^3+-1*im*c^-1*(a+b*arccot(x*c^-1))^3+-3//2*b^3*c^-1*Polylogarithms.polylog(3,1+-2*(1+im*c*x^-1)^-1)+-3*b*c^-1*(a+b*arccot(x*c^-1))^2*log(2*(1+im*c*x^-1)^-1)+-3*im*b^2*c^-1*(a+b*arccot(x*c^-1))*Polylogarithms.polylog(2,1+-2*(1+im*c*x^-1)^-1))
@test integrate(x^-3*(a+b*arctan(c*x^-1))^3, x) == :(-1//2*c^-2*(a+b*arccot(x*c^-1))^3+-1//2*x^-2*(a+b*arccot(x*c^-1))^3+3*b^2*c^-2*(a+b*arccot(x*c^-1))*log(2*(1+im*c*x^-1)^-1)+3//2*im*b*c^-2*(a+b*arccot(x*c^-1))^2+3//2*im*b^3*c^-2*Polylogarithms.polylog(2,1+-2*(1+im*c*x^-1)^-1)+3//2*b*c^-1*x^-1*(a+b*arccot(x*c^-1))^2)
@test integrate(x^2*arctan(x^(1/2)), x) == :(-1//3*x^(1/2)+-1//15*x^5//2+1//3*arctan(x^(1/2))+1//9*x^3//2+1//3*x^3*arctan(x^(1/2)))
@test integrate(x*arctan(x^(1/2)), x) == :((1/2)*x^(1/2)+-1//2*arctan(x^(1/2))+-1//6*x^3//2+(1/2)*x^2*arctan(x^(1/2)))
@test integrate(arctan(x^(1/2)), x) == :(-1*x^(1/2)+x*arctan(x^(1/2))+arctan(x^(1/2)))
@test integrate(x^-1*arctan(x^(1/2)), x) == :(im*Polylogarithms.polylog(2,-1*im*x^(1/2))+-1*im*Polylogarithms.polylog(2,im*x^(1/2)))
@test integrate(x^-2*arctan(x^(1/2)), x) == :(-1*x^-1//2+-1*arctan(x^(1/2))+-1*x^-1*arctan(x^(1/2)))
@test integrate(x^-3*arctan(x^(1/2)), x) == :((1/2)*x^-1//2+(1/2)*arctan(x^(1/2))+-1//6*x^-3//2+-1//2*x^-2*arctan(x^(1/2)))
@test integrate(x^3//2*arctan(x^(1/2)), x) == :(-1//5*log(1+x)+-1//10*x^2+1//5*x+2//5*x^5//2*arctan(x^(1/2)))
@test integrate(x^(1/2)*arctan(x^(1/2)), x) == :(-1//3*x+1//3*log(1+x)+2//3*x^3//2*arctan(x^(1/2)))
@test integrate(x^-1//2*arctan(x^(1/2)), x) == :(-1*log(1+x)+2*x^(1/2)*arctan(x^(1/2)))
@test integrate(x^-3//2*arctan(x^(1/2)), x) == :(-1*log(1+x)+-2*x^-1//2*arctan(x^(1/2))+log(x))
@test integrate(x^-5//2*arctan(x^(1/2)), x) == :(-1//3*x^-1+-1//3*log(x)+1//3*log(1+x)+-2//3*x^-3//2*arctan(x^(1/2)))
@test integrate(x^-1*arctan(a*x^5), x) == :(-1//10*im*Polylogarithms.polylog(2,im*a*x^5)+1//10*im*Polylogarithms.polylog(2,-1*im*a*x^5))
@test integrate(x^-1*arctan(a*x^n), x) == :((1/2)*im*n^-1*Polylogarithms.polylog(2,-1*im*a*x^n)+-1//2*im*n^-1*Polylogarithms.polylog(2,im*a*x^n))

max_stars_repo_path: test/5 Inverse trig functions/5.3 Inverse tangent/5.3.2 (d x)^m (a+b arctan(c x^n))^p.jl


text:
# another tree test

using IncrementalInference
using Test

@testset "Variable ordering Bayes tree member check." begin

  fg = generateCanonicalFG_Kaess()
  # Choose specific variable ordering and perform check.
  vo = [:l1, :l2, :x1, :x2, :x3]
  tree = buildTreeReset!(fg, vo)
  @test vo == tree.eliminationOrder
  @test vo == getEliminationOrder(tree)

end

@testset "Test Caesar Ring 1D symbolic tree construction" begin

fg = generateCanonicalFG_CaesarRing1D()
# drawGraph(fg, show=true)

eo = [:x0;:x2;:x4;:x6;:x1;:l1;:x5;:x3;]

tree = buildTreeReset!(fg,eo)
# drawTree(tree, show=true)

@test length(tree.cliques) == 6

C0 = getClique(tree, :x3)
@test intersect( C0 |> getFrontals, [:x3; :x5; :l1]) |> length == 3
@test C0 |> getCliqSeparatorVarIds |> length == 0

cC0 = getChildren(tree, C0)
@test cC0 |> length == 3

C1 = getClique(tree, :x1)
@test C1 in cC0
@test C1 |> getFrontals == [:x1;]
@test intersect(C1 |> getCliqSeparatorVarIds, [:x3;:l1]) |> length == 2

cC1 = getChildren(tree, C1)
@test cC1 |> length == 2

C4 = getClique(tree, :x2)
@test C4 in cC1
@test C4 |> getFrontals == [:x2;]
@test intersect(C4 |> getCliqSeparatorVarIds, [:x3;:x1]) |> length == 2

C5 = getClique(tree, :x0)
@test C5 in cC1
@test C5 |> getFrontals == [:x0;]
@test intersect(C5 |> getCliqSeparatorVarIds, [:l1;:x1]) |> length == 2

C2 = getClique(tree, :x6)
@test C2 in cC0
@test C2 |> getFrontals == [:x6;]
@test intersect(C2 |> getCliqSeparatorVarIds, [:l1;:x5]) |> length == 2

C3 = getClique(tree, :x4)
@test C3 in cC0
@test C3 |> getFrontals == [:x4;]
@test intersect(C3 |> getCliqSeparatorVarIds, [:x3;:x5]) |> length == 2

end

@testset "Test tree formation and consistent APIs" begin

fg = generateCanonicalFG_TestSymbolic()

#writeGraphPdf(fg, show=true)

eo = [:x1; :l3; :l1; :x5; :x2; :l2; :x4; :x3]

tree = buildTreeReset!(fg,eo)
# drawTree(tree, show=true)

@warn "TODO, complete further testing on tree formation"

## test variable order APIs consistent, see issue 499

vo = getEliminationOrder(fg)
tree1 = buildTreeReset!(fg, vo)
# drawTree(tree1, show=true)

tree2 = buildTreeReset!(fg)
# drawTree(tree2, show=true)

@test getEliminationOrder(tree1) == getEliminationOrder(tree1)

end

#


text:
using Documenter
using EvoTrees

push!(LOAD_PATH,"../src/")

pages = ["Home" => "index.md",
    "Examples" => "examples.md"]

makedocs(
    sitename="EvoTrees.jl",
    authors = "Jeremie Desgagne-Bouchard and contributors.",
    format=Documenter.HTML(),
    pages = pages,
    modules = [EvoTrees],)

deploydocs(repo="github.com/Evovest/EvoTrees.jl.git",
    target="build",
    push_preview=false)


text:
import ProximalCore
import RecursiveArrayTools

@inline function ProximalCore.prox(h, x::RecursiveArrayTools.ArrayPartition, gamma...)
    # unwrap
    y, fy = ProximalCore.prox(h, x.x, gamma...)
    # wrap
    return RecursiveArrayTools.ArrayPartition(y), fy
end

@inline function ProximalCore.gradient(h, x::RecursiveArrayTools.ArrayPartition)
    # unwrap
    grad, fx = ProximalCore.gradient(h, x.x)
    # wrap
    return RecursiveArrayTools.ArrayPartition(grad), fx
end

@inline ProximalCore.prox!(y::RecursiveArrayTools.ArrayPartition, h, x::RecursiveArrayTools.ArrayPartition, gamma...) = ProximalCore.prox!(y.x, h, x.x, gamma...)

@inline ProximalCore.gradient!(y::RecursiveArrayTools.ArrayPartition, h, x::RecursiveArrayTools.ArrayPartition) = ProximalCore.gradient!(y.x, h, x.x)


text:
using TrajectoryOptimization
using BenchmarkTools
using ForwardDiff
using LinearAlgebra
using StaticArrays
using Test
using TrajOptPlots
using MeshCat
const TO = TrajectoryOptimization

model = Dynamics.DubinsCar()
if !isdefined(Main,:vis)
    vis = Visualizer(); open(vis);
    set_mesh!(vis, model)
end

models = Dynamics.CopyModel(model, 2)
size(models) == (6,4)
n,m = size(models)
zeros(models)

x,u = rand(models)
x1,u1 = x[1:3], u[1:2]
x2,u2 = x[4:6], u[3:4]
@test !(states(models, x, 1) ≈ states(models, x, 2))
@test !(states(models, x, 1) ≈ states(models, x, 2))

@test dynamics(models, x, u) == [dynamics(model, x1,u1); dynamics(model,x2,u2)]
z = KnotPoint(x,u,0.1)
z1 = KnotPoint(x1,u1,0.1)
z2 = KnotPoint(x2,u2,0.1)
@test discrete_dynamics(RK3, models, z) ==
    [discrete_dynamics(model, z1);discrete_dynamics(model, z2)]

function my_jacobian(z::SVector)
    n,m = size(models)
    ix,iu = 1:n, n .+ (1:m)
    f_aug(z) = dynamics(models, view(z,ix), view(z,iu))
    ForwardDiff.jacobian(f_aug, z)
end
@test jacobian(models, z) ≈ my_jacobian(z.z)

function my_discrete_jacobian(s::SVector{NM1}) where NM1
    t = 0.0
    ix = @SVector [1,2,3,4,5,6]
    iu = @SVector [7,8,9,10]
    idt = NM1
    fd_aug(s) = discrete_dynamics(RK3, models, s[ix], s[iu], t, s[idt])
    ForwardDiff.jacobian(fd_aug, s)
end
s = [z.z; @SVector [z.dt]]
@test my_discrete_jacobian(s) ≈ discrete_jacobian(RK3, models, z)

xf = rand(models)[1]
@test TO.state_diff(models, x, xf) == x - xf
@test TO.state_diff_jacobian(models, x) == I
@test TO.∇²differential(models, x, @SVector rand(6)) == I*0

# Test solve
N = 101
tf = 5.0

x0_1 = @SVector [ 1,0,pi/2]
x0_2 = @SVector [-1,0,pi/2]
x0 = Dynamics.build_state(models, x0_1, x0_2)

xf_1 = @SVector [-1,2,pi/2]
xf_2 = @SVector [ 1,2,pi/2]
xf = Dynamics.build_state(models, xf_1, xf_2)

Q = Diagonal(@SVector fill(1e-2,n))
Rd = @SVector [0.1, 1e-2]
R = Diagonal( [Rd; Rd] )
obj = LQRObjective(Q,R,Q*100,xf,N)

conSet = ConstraintSet(n,m,N)
con = TO.CollisionConstraint(n, (@SVector [1,2]), (@SVector [4,5]), 0.5)
add_constraint!(conSet, con, 1:N-1)

prob = Problem(models, obj, xf, tf, x0=x0, constraints=conSet)
# solver = iLQRSolver(prob)
solver = AugmentedLagrangianSolver(prob)
@test size(solver) == (6,4,N)
solve!(solver)
@test iterations(solver) < 20  # should be 16
@test max_violation(solver) < 1e-10 # should be 0

prob = Problem(models, obj, xf, tf, x0=x0, constraints=conSet)
solver = iLQRSolver2(prob)
solve!(solver)
visualize!(vis, solver)

# Test with rigid bodies
Rot = UnitQuaternion{Float64,CayleyMap}
model = Dynamics.FreeBody{Rot,Float64}()
models = Dynamics.CopyModel(model, 2)

x,u = zeros(models)
x = repeat(zeros(model)[1],2)
u = repeat(zeros(model)[2],2)
x,u = rand(models)
x1,u1 = SVector{13}(x[1:13]), SVector{6}(u[1:6])
x2,u2 = SVector{13}(x[14:26]), SVector{6}(u[7:12])
@test !(states(models, x, 1) ≈ states(models, x, 2))
@test !(states(models, x, 1) ≈ states(models, x, 2))

@test dynamics(models, x, u) == [dynamics(model, x1,u1); dynamics(model,x2,u2)]
z = KnotPoint(x,u,0.1)
z1 = KnotPoint(x1,u1,0.1)
z2 = KnotPoint(x2,u2,0.1)
@test discrete_dynamics(RK3, models, z) ==
    [discrete_dynamics(model, z1);discrete_dynamics(model, z2)]

function my_jacobian(z::SVector)
    n,m = size(models)
    ix,iu = 1:n, n .+ (1:m)
    f_aug(z) = dynamics(models, view(z,ix), view(z,iu))
    ForwardDiff.jacobian(f_aug, z)
end
@test jacobian(models, z) ≈ my_jacobian(z.z)

function my_discrete_jacobian(s::SVector{NM1}) where NM1
    t = 0.0
    ix = z._x
    iu = z._u
    idt = NM1
    fd_aug(s) = discrete_dynamics(RK3, models, s[ix], s[iu], t, s[idt])
    ForwardDiff.jacobian(fd_aug, s)
end
s = [z.z; @SVector [z.dt]]
@test my_discrete_jacobian(s) ≈ discrete_jacobian(RK3, models, z)

xf = rand(models)[1]
x1f = SVector{13}(xf[1:13])
x2f = SVector{13}(xf[14:26])
@test TO.state_diff(models, x, xf) ==
    [TO.state_diff(model, x1, x1f); TO.state_diff(model, x2, x2f)]
G1 = TO.state_diff_jacobian(model, x1)
G2 = TO.state_diff_jacobian(model, x2)
@test cat(G1,G2,dims=[1,2]) ≈ TO.state_diff_jacobian(models, x)
b1 = @SVector rand(12)
b2 = @SVector rand(12)
∇G1 = TO.∇²differential(model, x1, b1)
∇G2 = TO.∇²differential(model, x2, b2)
@test cat(∇G1, ∇G2, dims=[1,2]) ≈ TO.∇²differential(models, x, [b1;b2])

# Test solve
N = 101
tf = 5.0

x01 = Dynamics.build_state(model, [0, 1,0], UnitQuaternion(I), zeros(3), zeros(3))
x02 = Dynamics.build_state(model, [0,-1,0], UnitQuaternion(I), zeros(3), zeros(3))
x0 = [x01; x02]

xf1 = Dynamics.build_state(model, [ 1,0,0], expm(deg2rad( 45)*@SVector [1,0,0.]), zeros(3), zeros(3))
xf2 = Dynamics.build_state(model, [-1,0,0], expm(deg2rad(-45)*@SVector [1,0,0.]), zeros(3), zeros(3))
xf = [xf1; xf2]

Qd = Dynamics.fill_state(model, 1e-1, 1e-1, 1e-2, 1e-2)
Q = Diagonal([Qd; Qd])
Rd = @SVector fill(1e-2, 6)
R = Diagonal([Rd; Rd])
obj = LQRObjective(Q,R,Q*10,xf,N)

prob = Problem(models, obj, xf, tf, x0=x0)
solver = iLQRSolver(prob)
solver.opts.verbose = true
solve!(solver)

if !isdefined(Main,:vis)
    vis = Visualizer(); open(vis);
    set_mesh!(vis, model)
end
visualize!(vis, models, get_trajectory(solver))


text:
using StackOverflow

mutable struct Questions
    link::Vector
    view_count::Vector
    creation_date::Vector
    is_answered::Vector
    owner::Vector
    last_activity_date::Vector
    score::Vector
    accepted_answer_id::Vector
    question_id::Vector
    tags::Vector
    title::Vector
    answer_count::Vector
end

"""
    getrecentquestionsfortag(;tag::String = "Julia", site::String = "stackoverflow", order::String = "desc",
                                    sort::String = "activity", todate::String = "", fromdate::String = "")

Will return the 30 StackOverflow questions given a certain tag and sort parameter.

Possible order values: "desc" or "asc"
Possible sort values: "activity", "votes", "creation", "hot", "week", or "month"
"""
function getrecentquestionsfortag(;tag::String = "Julia", site::String = "stackoverflow", order::String = "desc",
                                    sort::String = "activity", todate::String = "", fromdate::String = "")

    questionholder = Vector()

    if fromdate != "" && todate != ""
        datequery = "fromdate=$(fromdate)&todate=$(todate)&"  #?fromdate=1573776000&todate=1574467200&order
        r = HTTP.request("GET", "https://api.stackexchange.com/2.2/questions?$(datequery)order=$(order)&sort=$(sort)&tagged=$(tag)&site=$(site)")
    else
        r = HTTP.request("GET", "https://api.stackexchange.com/2.2/questions?order=$(order)&sort=$(sort)&tagged=$(tag)&site=$(site)")
    end

    json = convert_HTTP_Response_To_JSON(r)

    for (k, v) in json
        if occursin("items", k)
            for item in v

                link                = Vector()
                view_count          = Vector()
                creation_date       = Vector()
                is_answered         = Vector()
                owner               = Vector()
                last_activity_date  = Vector()
                score               = Vector()
                accepted_answer_id  = Vector()
                question_id         = Vector()
                tags                = Vector()
                title               = Vector()
                answer_count        = Vector()

                for (key, value) in item
                    if occursin("link", key)
                        push!(link, value)
                    elseif occursin("view_count", key)
                        push!(view_count, value)
                    elseif occursin("creation_date", key)
                        push!(creation_date, value)
                    elseif occursin("is_answered", key)
                        push!(is_answered, value)
                    elseif occursin("owner", key)
                        push!(owner, value)
                    elseif occursin("last_activity_date", key)
                        push!(last_activity_date, value)
                    elseif occursin("score", key)
                        push!(score, value)
                    elseif occursin("accepted_answer_id", key)
                        push!(accepted_answer_id, value)
                    elseif occursin("question_id", key)
                        push!(question_id, value)
                    elseif occursin("tags", key)
                        push!(tags, value)
                    elseif occursin("title", key)
                        push!(title, value)
                    elseif occursin("answer_count", key)
                        push!(answer_count, value)
                    end
                    # println(key)
                end
                question = Questions(link, view_count, creation_date, is_answered, owner,
                    last_activity_date, score, accepted_answer_id, question_id, tags, title,
                    answer_count)

                push!(questionholder, question)
            end

            return questionholder
        end
    end

end


text:
#!/usr/bin/env julia
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 4., Zadanie 2. — podstawowy test (przykład z wykładu)
#

include("interpolation.jl")
using .Interpolation

FL = Interpolation.FL

x = Vector{FL}([3, 1, 5, 6])
f = Vector{FL}([1, -3, 2, 4])

fx = Interpolation.ilorazyRoznicowe(x, f)

ft = Interpolation.warNewton(x, fx, 1.0)

println(ft)

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
# this user supplied file creates a dictionary of arguments
# if a key is repeated, the last use of the key is used
# it is a little bit dangerous letting the user run arbitrary code
# as part of the solver
# now that this file is read inside a function, it is better encapsulated

arg_dict = Dict{Any, Any} (
"run_type" => 10,
"order" => 1,
"real_time" => true,
"IC_name" => "ICsinwave",
"numBC" => 1,
"BC1" => [0],
"BC1_name" => "sinwaveBC",
#"BC2" => [4, 10],
#"BC2_name" => "noPenetrationBC",
"delta_t" => 2*pi/5000,
"t_max" => 100*2*pi/5000,
#"smb_name" => "src/mesh_files/channel.smb",
"smb_name" => "src/mesh_files/tri8l.smb",
#"smb_name" => "src/mesh_files/channel_verylarge.smb",
#"dmg_name" => "../../mesh_files/vortex.dmg",
"dmg_name" => ".null",
"res_tol" => 1e-10,
# "step_tol" => 1e-6,
"itermax" => 20,
"res_abstol" => 1e-10,
"res_reltol" => 1e-9,
"step_tol" => 1e-10,
"itermax" => 30,
"writeq" => false,
"itermax" => 30,
"writeq" => false,
"write_edge_vertnums" => false,
"write_face_vertnums" => false,
"write_qic" => false,
"writeboundary" => false,
"write_res" => false,
"print_cond" => false,
"write_counts" => false,
"write_vis" => true,
"output_freq" => 1,
"solve" => true,
"do_postproc" => true,
"exact_soln_func" => "ICsinwave",
"write_face_vertnums" => true
)


text:
@testset "$f" for f in (density_reference,density)
	@testset "twopoints" begin
		x = Float64[0,1]
		xeval = Float64[0]
		C = 1/(√(2π)*length(x))
		@test f(x,1/√2,xeval) ≈ [1+exp(-1)]*C*√2
		@test f(x,1,xeval) ≈ [1+exp(-1/2)]*C
		@test f(x,√2,xeval) ≈ [1+exp(-1/4)]*C/√2

		@test f(x.+1,1,xeval.+1) ≈ [1+exp(-1/2)]*C # translation
		@test f(x.*5,5,xeval) ≈ [1+exp(-1/2)]*C/5 # scaling
		@test f(-x,1,xeval) ≈ [1+exp(-1/2)]*C # mirror
	end

	@testset "fewpoints" begin
		x = [-4, 3.5, 10, 12]
		xeval = [0.0]
		bandwidths = 10.0.^(-3:3)
		ground_truth = [1.9824062450251952342e-2660052, 2.8962330179366389711e-26600, 9.8509906928407312942e-267, 0.00025162823020266138054, 0.029491681668065441022, 0.0039758850413688216497, 0.00039892870428519092503]

		@test all(isapprox.(f(x,bandwidths,xeval), ground_truth; atol=1e-20, rtol=1e-3))
		@test all(isapprox.(f(x.*-3.1.-7.8,bandwidths.*3.1,xeval.-7.8), ground_truth/3.1; atol=1e-20, rtol=1e-3))
	end
end

@testset "$f" for f in (smooth_reference,smooth)
	@testset "twopoints" begin
		x = Float64[0,1]
		y = Float64[0,1]
		xeval = Float64[0]
		@test f(x,y,1/√2,xeval) ≈ [1/(ℯ+1)]
		@test f(x,y,1,xeval) ≈ [1/(√ℯ+1)]
		@test f(x,y,√2,xeval) ≈ [1/(exp(1/4)+1)]

		@test f(x.+1,y,1,xeval.+1) ≈ [1/(√ℯ+1)] # translation
		@test f(x.*5,y,5,xeval) ≈ [1/(√ℯ+1)] # scaling
		@test f(-x,y,1,xeval) ≈ [1/(√ℯ+1)] # mirror

		@test f(x,y.+1,1,xeval) ≈ [1+1/(√ℯ+1)] # translation
		@test f(x,3y,1,xeval) ≈ [3/(√ℯ+1)] # scaling
		@test f(x,-y,1,xeval) ≈ -[1/(√ℯ+1)] # mirror

		# extreme bandwidths
		@test f(x,y,1e-5,[0.5]) ≈ [0.5]
		@test f(x,y,1e5,[0.5]) ≈ [0.5]
		@test f(x,y,1e-100,[0.5]) ≈ [0.5]
		@test f(x,y,1e100,[0.5]) ≈ [0.5]

		@test f(x,y,1e-5,[0.6]) ≈ [1.0]
		@test f(x,y,1e5,[0.6]) ≈ [0.5]
		@test f(x,y,1e-100,[0.6]) ≈ [1.0]
		@test f(x,y,1e100,[0.6]) ≈ [0.5]

		@test f(x,y,1e-5,[0.4]) ≈ [0.0] atol=1e-9
		@test f(x,y,1e5,[0.4]) ≈ [0.5]
		@test f(x,y,1e-100,[0.4]) ≈ [0.0] atol=1e-9
		@test f(x,y,1e100,[0.4]) ≈ [0.5]
	end

	@testset "fewpoints" begin
		x = [-4, 3.5, 10, 12]
		y = [1.0, -2.8, 8, 13]
		xeval = [0.0]
		bandwidths = 10.0.^(-3:3)
		ground_truth = [-2.79999999999999982236, -2.79999999999999982236, -2.79999999999999982236, -2.29473588724824868121, 3.20241290611236881292, 4.78316896715214792770, 4.79983163801801130898]

		@test all(isapprox.(f(x,y,bandwidths,xeval), ground_truth; rtol=1e-3))
		@test all(isapprox.(f(x.*-3.1.-7.8,y,bandwidths.*3.1,xeval.-7.8), ground_truth; rtol=1e-3))
	end
end

@testset "npoints=$N" for N in [99,100,101,995,1000,1005,9997,10000,10003]
	rng = StableRNG(42)
	x = randn(rng, N)
	y = randn(rng, N)
	xeval = randn(rng,10)

	bandwidths = 10.0.^ [0,-1,1,-2,2,-3,3,-5,5,-10,10,-100,100]
	@testset "bandwidth=$bw" for bw in bandwidths
		gtd = density_reference(x,bw,xeval)

		@test all(isapprox.(density(x,bw,xeval,rtol=1e-3), gtd; atol=1e-12, rtol=1e-3))
		@test all(isapprox.(density(x,bw,xeval,rtol=1e-6), gtd; atol=1e-12, rtol=1e-6))
		@test all(isapprox.(density(x,bw,xeval,rtol=1e-9), gtd; atol=1e-12, rtol=1e-9))

		gt = smooth_reference(x,y,bw,xeval)
		@test all(isapprox.(smooth(x,y,bw,xeval;rtol=1e-3), gt; rtol=1e-3))
		@test all(isapprox.(smooth(x,y,bw,xeval;rtol=1e-6), gt; rtol=1e-6))
		@test all(isapprox.(smooth(x,y,bw,xeval;rtol=1e-9), gt; rtol=1e-9))
		@test all(isapprox.(smooth(x,y,bw,xeval;atol=1e-1), gt; atol=1e-1))
		@test all(isapprox.(smooth(x,y,bw,xeval;atol=1e-3), gt; atol=1e-3))
		@test all(isapprox.(smooth(x,y,bw,xeval;atol=1e-6), gt; atol=1e-6))

		gks = GaussianKernelSmoother(x,y)
		@test all(isapprox.(density(gks,bw,xeval;rtol=1e-3), gtd; atol=1e-12, rtol=1e-3))
		@test all(isapprox.(density(gks,bw,xeval;atol=1e-3), gtd; atol=1e-3))

		@test all(isapprox.(smooth.(gks,bw,xeval;rtol=1e-3), gt; rtol=1e-3))
		@test all(isapprox.(smooth.(gks,bw,xeval;atol=1e-3), gt; atol=1e-3))
	end
end


text:
using Pkg
Pkg.activate(@__DIR__)

using JGraph
using Documenter

makedocs(;
    modules=[JGraph],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    sitename="JGraph.jl",
    authors="Vincent Gauthier",
    doctest  = true,
)

deploydocs(;
    target = "build",
    branch = "gh-pages",
    repo="github.com/vgauthier/JGraph.jl.git",
    devbranch = "master",
    devurl = "dev",
    versions = ["stable" => "v^", "v#.#", "dev" => "dev"]
)


text:
module Utils
# check that a package is installed
function checkpkgs(pkglist::Vararg{AbstractString, N}) where N
    missing = []
    for p in pkglist
        pkg = Symbol(p)
        try
            @eval using $pkg
        catch e
            push!(missing, p)
        end
    end
    if length(missing) > 0
        m = join(missing, ", ")
        error("""These packages are missing: $m
                 You can install with e.g. :
                    import Pkg; Pkg.add("$(missing[1])")
            """
        )
    end
end

checkpkgs("Downloads", "Plots", "Distributions", "DataFrames", "Loess", "Gumbo", "HTTP", "Cascadia")

import Downloads, DataFrames, Loess, Plots, Distributions
import Gumbo, Cascadia, HTTP, CSV

function getfile(url::String)
    filename = basename(url)
    if !isfile(filename)
        _local = Downloads.download(url, filename)
        println("Downloaded ", _local)
    end
end

default(linewidth=2, legend=false, size = (480, 320));

function Base.transpose(df::DataFrame; col_names=nothing)
    tdf = DataFrames.DataFrame(collect.(eachrow(df)), :auto)
    isnothing(col_names) || rename!(tdf, col_names)
    return tdf
end

function plotloess!(xs, ys; kwargs...)
    model = Loess.loess(xs, ys)
    us = range(extrema(xs)...; length = 101)
    vs = Loess.predict(model, us)
    Plots.plot!(us, vs; kwargs...)
end

# this is a hack
# unfortunately, there's no Julia equivalent for Pandas' read_html
function retrievetables(url::String; kwargs...)
    res = HTTP.request("GET",url)
    g = Gumbo.parsehtml(String(res.body))

    delim="∘"
    tbls = eachmatch(Cascadia.Selector("tbody"), g.root)
    v = []
    for tbl in tbls
        rows = children(tbl)
        io = IOBuffer()
        # this is an ugly hack: convert to CSV and let CSV deal with the mess
        # (hopefully, parsing types correctly)
        for row in rows
            entries = [strip(nodeText(c)) for c in children(row)]
            println(io, join(entries, delim))
        end
        seekstart(io)
        df = CSV.File(io, delim=delim; kwargs...) |> DataFrame
        push!(v, df)
    end
    return v
end

export checkpkgs, retrievetables, transpose, plotloess!, getfile
end;


text:
export extend, NfToNfMor, automorphisms

mutable struct NfToNfMor <: Map{AnticNumberField, AnticNumberField, HeckeMap, NfToNfMor}
  header::MapHeader{AnticNumberField, AnticNumberField}
  prim_img::nf_elem
  prim_preimg::nf_elem

  function NfToNfMor()
    z = new()
    z.header = MapHeader()
    return r
  end

  function NfToNfMor(K::AnticNumberField, L::AnticNumberField, y::nf_elem, isomorphism::Bool = false)
    z = new()
    z.prim_img = y

    function _image(x::nf_elem)
      g = parent(K.pol)(x)
      return evaluate(g, y)
    end

    if !isomorphism
      z.header = MapHeader(K, L, _image)
      return z
    end

    M = zero_matrix(FlintQQ, degree(L), degree(L))
    b = basis(K)
    for i = 1:degree(L)
      c = _image(b[i])
      for j = 1:degree(L)
        M[j, i] = coeff(c, j - 1)
      end
    end
    t = zero_matrix(FlintQQ, degree(L), 1)
    if degree(L) == 1
      t[1, 1] = coeff(gen(L), 0)
    else
      t[2, 1] = fmpq(1) # coefficient vector of gen(L)
    end

    s = solve(M, t)
    z.prim_preimg = K(parent(K.pol)([ s[i, 1] for i = 1:degree(K) ]))

    function _preimage(x::nf_elem)
      g = parent(L.pol)(x)
      return evaluate(g, z.prim_preimg)
    end

    z.header = MapHeader(K, L, _image, _preimage)
    return z
  end

  function NfToNfMor(K::AnticNumberField, L::AnticNumberField, y::nf_elem, y_inv::nf_elem)
    z = new()
    z.prim_img = y
    z.prim_preimg = y_inv

    function _image(x::nf_elem)
      g = parent(K.pol)(x)
      return evaluate(g, y)
    end

    function _preimage(x::nf_elem)
      g = parent(L.pol)(x)
      return evaluate(g, y_inv)
    end

    z.header = MapHeader(K, L, _image, _preimage)
    return z
  end
end

function inv(f::NfToNfMor)
  if degree(domain(f)) != degree(codomain(f))
    error("The map is not invertible")
  end
  if isdefined(f, :prim_preimg)
    return hom(codomain(f), domain(f), f.prim_preimg, check = false)
  end
  L = codomain(f)
  K = domain(f)
  M = zero_matrix(FlintQQ, degree(L), degree(L))
  a = f.prim_img
  el = one(L)
  M[1, 1] = 1
  for i = 2:degree(L)
    mul!(el, el, a)
    for j = 1:degree(L)
      M[j, i] = coeff(el, j - 1)
    end
  end
  t = zero_matrix(FlintQQ, degree(L), 1)
  if degree(L) == 1
    t[1, 1] = coeff(gen(L), 0)
  else
    t[2, 1] = fmpq(1) # coefficient vector of gen(L)
  end

  s = solve(M, t)
  img = K(parent(K.pol)([ s[i, 1] for i = 1:degree(K) ]))
  return hom(L, K, img, check = false)
end

function _compute_preimg(m::NfToNfMor)
  # build the matrix for the basis change
  K = domain(m)
  L = codomain(m)
  M = zero_matrix(FlintQQ, degree(L), degree(L))
  b = basis(K)
  for i = 1:degree(L)
    c = m(b[i])
    for j = 1:degree(L)
      M[j, i] = coeff(c, j - 1)
    end
  end
  t = zero_matrix(FlintQQ, degree(L), 1)
  t[2, 1] = fmpq(1) # coefficient vector of gen(L)
  s = solve(M, t)
  m.prim_preimg = K(parent(K.pol)([ s[i, 1] for i = 1:degree(K) ]))
  return m.prim_preimg
end

function Base.:(==)(f::NfToNfMor, g::NfToNfMor)
  if (domain(f) != domain(g)) || (codomain(f) != codomain(g))
    return false
  end

  return f.prim_img == g.prim_img
end

#_D = Dict()

function evaluate(f::fmpq_poly, a::nf_elem)
  R = parent(a)
  if iszero(f)
    return zero(R)
  end
  l = length(f) - 1
  s = R(coeff(f, l))
  for i in l-1:-1:0
    #s = s*a + R(coeff(f, i))
    mul!(s, s, a)
    # TODO (easy): Once fmpq_poly_add_fmpq is improved in flint, remove the R(..)
    add!(s, s, R(coeff(f, i)))
  end
  return s
end

function *(f::NfToNfMor, g::NfToNfMor)
#  global _D
#  _s = Base.stacktrace()[2:3]
#  if !(_s in keys(_D))
#    _D[_s] = true
#    println("Called ...")
#    Base.show_backtrace(stdout, Base.stacktrace()[2:3])
#  end
  codomain(f) == domain(g) || throw("Maps not compatible")

  a = gen(domain(f))
  y = g(f(a))

  return hom(domain(f), codomain(g), y, check = false)
end

function ^(f::NfToNfMor, b::Int)
  K = domain(f)
  @assert K == codomain(f)
  d = degree(K)
  b = mod(b, d)
  if b == 0
    return NfToNfMor(K, K, gen(K))
  elseif b == 1
    return f
  else
    bit = ~((~UInt(0)) >> 1)
    while (UInt(bit) & b) == 0
      bit >>= 1
    end
    z = f
    bit >>= 1
    while bit != 0
      z = z * z
      if (UInt(bit) & b) != 0
        z = z * f
      end
      bit >>= 1
    end
    return z
  end
end

Base.copy(f::NfToNfMor) = f

Base.hash(f::NfToNfMor, h::UInt) = Base.hash(f.prim_img, h)

function show(io::IO, h::NfToNfMor)
  if domain(h) == codomain(h)
    println(io, "Automorphism of ", domain(h))
  else
    println(io, "Injection of ", domain(h), " into ", codomain(h))
  end
  println(io, "defined by ", gen(domain(h)), " -> ", h.prim_img)
end

function automorphisms(K::AnticNumberField; copy::Bool = true)
  try
    Aut = _get_automorphisms_nf(K)::Vector{NfToNfMor}
    if copy
      return Base.copy(Aut)
    else
      return Aut
    end
  catch e
    if !isa(e, AccessorNotSetError)
      rethrow(e)
    end
  end
  if degree(K) == 1
    auts = NfToNfMor[hom(K, K, one(K))]
  else
    f = K.pol
    Kt, t = PolynomialRing(K, "t", cached = false)
    f1 = change_ring(f, Kt)
    divpol = Kt(nf_elem[-gen(K), K(1)])
    f1 = divexact(f1, divpol)
    lr = roots(f1, max_roots = div(degree(K), 2))
    Aut1 = Vector{NfToNfMor}(undef, length(lr)+1)
    for i = 1:length(lr)
      Aut1[i] = hom(K, K, lr[i], check = false)
    end
    Aut1[end] = id_hom(K)
    auts = closure(Aut1, degree(K))
  end
  _set_automorphisms_nf(K, auts)
  if copy
    return Base.copy(auts)
  else
    return auts
  end
end

function hom(K::AnticNumberField, L::AnticNumberField, a::nf_elem; check::Bool = true, compute_inverse::Bool = false)
 if check
   if !iszero(evaluate(K.pol, a))
     error("The data does not define a homomorphism")
   end
 end
 return NfToNfMor(K, L, a, compute_inverse)
end

function hom(K::AnticNumberField, L::AnticNumberField, a::nf_elem, a_inv::nf_elem; check::Bool = true)
 if check
   if !iszero(evaluate(K.pol, a))
     error("The data does not define a homomorphism")
   end
   if !iszero(evaluate(L.pol, a_inv))
     error("The data does not define a homomorphism")
   end
 end
 return NfToNfMor(K, L, a, a_inv)
end

id_hom(K::AnticNumberField) = hom(K, K, gen(K), check = false)

morphism_type(::Type{AnticNumberField}) = NfToNfMor

isinjective(m::NfToNfMor) = true

issurjective(m::NfToNfMor) = (degree(domain(m)) == degree(codomain(m)))

isbijective(m::NfToNfMor) = issurjective(m)


text:
using Gurobi, MathProgBase, Test

@testset "Env" begin
    env = Gurobi.Env()

    # Test supplying environment manually

    solver = GurobiSolver(env)
    # Check that the env for each model is the same
    m1 = MathProgBase.LinearQuadraticModel(solver)
    @test m1.inner.env === env
    m2 = MathProgBase.LinearQuadraticModel(solver)
    @test m2.inner.env === env
    # Check that finalizer doesn't touch env when manually provided
    finalize(m1.inner)
    @test Gurobi.is_valid(env)

    # Test creating environment automatically

    solver = GurobiSolver()
    # Check that the env for each model is different
    m3 = MathProgBase.LinearQuadraticModel(solver)
    @test m3.inner.env !== env
    m4 = MathProgBase.LinearQuadraticModel(solver)
    @test m4.inner.env !== env
    @test m3.inner.env !== m4.inner.env
    # Check that env is finalized with model when not supplied manually
    finalize(m3.inner)
    @test !Gurobi.is_valid(m3.inner.env)
end