text:
using Revise
using CmdStan
using Plots
using StatsPlots
using Random

################################################################################
############################# PREPARATION ######################################

#%% Input Data #################################################################
y = [607, 583, 521, 494, 369, 782, 570, 678, 467, 620, 425, 395, 346, 361, 310,
300, 382, 294, 315, 323, 421, 339, 398, 328, 335, 291, 329, 310, 294, 321, 286,
349, 279, 268, 293, 310, 259, 241, 243, 272, 247, 275, 220, 245, 268, 357, 273,
301, 322, 276, 401, 368, 149, 507, 411, 362, 358, 355, 362, 324, 332, 268, 259,
274, 248, 254, 242, 286, 276, 237, 259, 251, 239, 247, 260, 237, 206, 242, 361,
267, 245, 331, 357, 284, 263, 244, 317, 225, 254, 253, 251, 314, 239, 248, 250,
200, 256, 233, 427, 391, 331, 395, 337, 392, 352, 381, 330, 368, 381, 316, 335,
316, 302, 375, 361, 330, 351, 186, 221, 278, 244, 218, 126, 269, 238, 194, 384,
154, 555, 387, 317, 365, 357, 390, 320, 316, 297, 354, 266, 279, 327, 285, 258,
267, 226, 237, 264, 510, 490, 458, 425, 522, 927, 555, 550, 516, 548, 560, 545,
633, 496, 498, 223, 222, 309, 244, 207, 258, 255, 281, 258, 226, 257, 263, 266,
238, 249, 340, 247, 216, 241, 239, 226, 273, 235, 251, 290, 473, 416, 451, 475,
406, 349, 401, 334, 446, 401, 252, 266, 210, 228, 250, 265, 236, 289, 244, 327,
274, 223, 327, 307, 338, 345, 381, 369, 445, 296, 303, 326, 321, 309, 307, 319,
288, 299, 284, 278, 310, 282, 275, 372, 295, 306, 303, 285, 316, 294, 284, 324,
264, 278, 369, 254, 306, 237, 439, 287, 285, 261, 299, 311, 265, 292, 282, 271,
268, 270, 259, 269, 249, 261, 425, 291, 291, 441, 222, 347, 244, 232, 272, 264,
190, 219, 317, 232, 256, 185, 210, 213, 202, 226, 250, 238, 252, 233, 221, 220,
287, 267, 264, 273, 304, 294, 236, 200, 219, 276, 287, 365, 438, 420, 396, 359,
405, 397, 383, 360, 387, 429, 358, 459, 371, 368, 452, 358, 371];

ind = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8,
 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,
 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11,
  11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
   12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14,
   14, 14, 14, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16,
   16, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 19, 19, 19, 19, 19, 20, 20, 20,
   20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 21, 21,
   21, 21, 21, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22,
   22, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 24, 24, 24,
   24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 25, 25, 25, 25, 25, 25, 25,
   25, 25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 28, 28, 28,
   28, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 30, 30, 30,
   30, 30, 30, 31, 31, 31, 31, 31, 32, 32, 32, 32, 32, 33, 34, 34, 34, 34, 34,
   34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34];

#%% Functions ##################################################################
function mySliceSampler(pdf_log, x_0, w, m, N = 2000, burnIn = 1000)
    # REQUIRES MODULES:
    # Random
    #
    # INPUT:
    # pdf_log --> log-pdf of the target distribution (a function)
    # x_init --> inital vector of x (dimension 1xD), where D dimension of pdf
    # w --> typical slice size, for a window (dimension 1xD)
    # m --> integer limiting the slice size to m*w
    # N --> number of sample points
    # burnIn --> number of burn-in samples (optional, results will be discarded)
    #
    # OUTPUT:
    # x_s --> sampled x values (NxD), only keep values after full permutations along dimensions
    # pdflog_x_s --> log-pdf of sampled values

    D = length(x_0)  # the dimension of the distribution (nr of coordinates)
    x_s = Array{Float64}(undef,N,D)   # Samples (NxD), preallocate array, with undefined values
    pdflog_x_s = Array{Float64}(undef,N)   # (unnormalised) log-prob of sampled points (Nx1)
                                           # for reducing evaluations
    pdflog_x_1 = 0.0;

    # to be able to input scalars:
    #if typeof(x_0) == Float64
    #    x_0 = [x_0]
    #end

    for ii in 1:(N+burnIn)
        # Update the new x progressively for each dimension:
        L   = 1*x_0  # R and L need to be of same dimension, because we apply the pdf upon those!
        R   = 1*x_0
        x_1 = 1*x_0  # "1*" needed, otherwise x_0 suddendly takes mysterious values
        for d in randperm(D)
            ### 1) Make the Slice
            # random "vertical" position
            if (ii == 1)
                #vertical = rand() * pdf_log(x_0) # in normal space
                z = pdf_log(x_0) + log(rand()) # in log-space
            else
                # saving evaluations, pdflog_x_1 was updated for current permutation in the end
                z = pdflog_x_1 + log(rand())
            end

            ### 2) Make the window, "Stepping-out" alogrithm
            # Randomly place a window, containing x_0:
            L[d] = x_0[d] - rand() * w[d]
            R[d] = L[d] + w[d]

            # Randomly share the max window size among left/right:
            J = floor(m*rand());
            K = (m-1) - J;

            # Extend window to the left, until outside slice or allowance seizes:
            while ((J > 0) && (z < pdf_log(L)))
                #println("Lefting")
                L[d] -= w[d]
                J -= 1
            end

            # Extend window to the right, until outside slice or allowance seizes:
            while ((K > 0) && (z < pdf_log(R)))
                #println("Righting")
                R[d] += w[d]
                K -= 1
            end

            ### 3) Sample from window
            # finding an allowable point:
            while true
                x_1[d] = L[d] + rand() * (R[d] - L[d])

                # this + breaking out of loop reduces the amount of log-pdf evaluations:
                pdflog_x_1 = pdf_log(x_1)

                if (pdflog_x_1 >= z)     # new value found
                    x_0[d] = x_1[d]     # update value for this dimension
                    break
                # Value was not within slice, shrink the interval:
                elseif (x_1[d] < x_0[d])
                    L[d] = x_1[d]
                elseif (x_1[d] > x_0[d])
                    R[d] = x_1[d]
                else
                    throw(ErrorException("Error during shrinking the interval"))
                end
            end

        end # end of permutating through dimensions

        ### 4) Update the chain:
        # Just overwrite the burnIn by using modular arithmetic:
        i = 1 + (ii-1)%N    # gives 1,2,3,...N,1,2,3,...,(N+burnIn)
        x_s[i,:] = x_1    # sample point, after one round of permutations
        pdflog_x_s[i] = pdflog_x_1 # now contains log_pdf of all new coordinates of new point, since has been updated at each permutation
        # note that x_0 has been updated during the permutations

    end

    return x_s, pdflog_x_s
end

## Jespers HDI:
function hdi(theta_samp,alpha=0.05)
    cred_mass = 1.0-alpha
    ci = zeros(2)
    if length(size(theta_samp))>1   # multidimensional theta
        K,N = size(theta_samp)
        cis = zeros(2,K)
        for k in 1:K    # for each dimension

            ts = theta_samp[k,:]
            sind = sortperm(ts)     # get series of indices that sorts theta values
            sts = ts[sind]  # sorted theta values

            # Shifting a 95% index block from far left to far right and use the
            # samllest interval that can be found during the shift
            N = length(sind)    # number of sample points
            length_ci = Inf     # start with inf HDI interval and shrink
            for i in 1:Int(floor(N*alpha))  # iterate over the 5% lowest values (low to high)
                i2 = Int(floor(N*cred_mass)+i)  # the 5% highest values (low to high), one 95% block away on the other side
                prop_ci = [sts[i],sts[i2]]      # proposed interval
                length_prop_ci = prop_ci[2]-prop_ci[1]
                if length_prop_ci < length_ci
                    ci = prop_ci
                    length_ci = ci[2]-ci[1]
                end
            end
            cis[:,k] = ci

        end
        return cis
    else
        N = length(theta_samp)

        ts = theta_samp
        sind = sortperm(ts)
        sts = ts[sind]

        N = length(sind)
        length_ci = Inf
        for i in 1:Int(floor(N*alpha))
            i2 = Int(floor(N*cred_mass)+i)
            prop_ci = [sts[i],sts[i2]]
            length_prop_ci = prop_ci[2]-prop_ci[1]
            if length_prop_ci < length_ci
                ci = prop_ci
                length_ci = ci[2]-ci[1]
            end
        end
        return ci
    end
end

# A recursive application of the HDI algorithm to find the mode:
function find_mode(sample_pts;runs=2)
    # re-apply the hdi algorithm, until there is only two points or so left
    # so the ci interval becomes smaller and smaller
    # in the end, return an average over the interval
    ci = zeros(2)
    i_left = 0
    i_right = 0

    # Make 50% credibility intervals
    alpha = 0.5;
    cred_mass = 1.0 - alpha;

    # sort the pts:
    sind = sortperm(sample_pts)
    sts = sample_pts[sind]
    N = length(sind)

    if N > runs
        # HDI algorithm, but return indices too:
        length_ci = Inf
        for i in 1:Int(floor(N*alpha))
            i2 = Int(floor(N*cred_mass)+i)
            prop_ci = [sts[i], sts[i2]]
            length_prop_ci = prop_ci[2] - prop_ci[1]
            if length_prop_ci < length_ci
                i_left  = i
                i_right = i2
                ci = prop_ci
                length_ci = ci[2] - ci[1]
            end
        end

        # Recursion:
        ci = find_mode(sts[i_left:i_right],runs=runs)

    else
        ci = [sts[1],sts[end]]
    end

    return ci
end

function log_prior_flat(X::Float64)
    return 0.0  # or any other constant value, for the posterior it doesn't matter
end

function log_prior_flat_positiveOnly(X::Float64)
    if X > 0
        return 0.0  # or any other constant value, for the posterior it doesn't matter
    else
        return -Inf
    end
    # Alternative of Jesper:
    #return log(tau > 0)
end

function log_normal(θ::Array{Float64},μ::Float64,τ::Float64)
    # Return the joint pdf of an array of thetas, i.e. the sum in log-space

    # We can leave out the 1/sqrt(2π) since it does not change sampling from the posterior
    return sum(-log(τ) .- 0.5 .* ((θ .- μ) ./ τ).^2)
end

function log_lklhd_fct(jointParams::Array{Float64},zlogy::Array{Float64,1},ind::Array{Int64,1})
    # return the joint lklhd of the given input data
    J = length(jointParams) - 3
    θ = jointParams[1:J]
    μ, σ, τ  = jointParams[(J+1):end]

    jointLklhd = 0.0

    if (σ > 0)
        for i in 1:length(zlogy)
            jointLklhd += -log(σ) - 0.5 * ((zlogy[i] - θ[ind[i]]) / σ)^2
        end
    else
        jointLklhd = -Inf
    end

    return jointLklhd
end

function log_priors(jointParams::Array{Float64})
    # The joint log-pdf of all the priors in the model
    J = length(jointParams) - 3
    θ = jointParams[1:J]
    μ, σ, τ  = jointParams[(J+1):end]

    if (σ > 0) && (τ > 0)
        return log_normal(θ,μ,τ) + log_prior_flat(μ) + log_prior_flat_positiveOnly(σ) + log_prior_flat_positiveOnly(τ)
    else
        return -Inf
    end
end

#%% Logarithmise and Transform Data ############################################
J = maximum(ind)    # number of individuals
I = length(y)       # number of observations

# Go to log-space:
logy = log.(y)

logMean = mean(logy)
logStd = std(logy)

# Mean-centre and scale each element:
zlogy = (logy .- logMean) ./ logStd

#%% Concretise the functions by including the data #############################
log_lklhd(jointParams::Array{Float64}) = log_lklhd_fct(jointParams::Array{Float64}, zlogy, ind)
log_posterior(jointParams::Array{Float64}) = log_lklhd(jointParams::Array{Float64}) + log_priors(jointParams::Array{Float64})

################################################################################
############################# MCMC #############################################

#%% Chain Setup
N = 10^6
burnIn = 10^4
x_start = ones(Float64,J+3) * 0.1
w = ones(Float64,J+3) * 0.1
m = 100

#%% MCMC run
@time chain, pdf = mySliceSampler(log_posterior,x_start,w,m,N,burnIn);
# takes around 550 seconds on my machine (4 cores, i7-7500U 2,7 GHz x 2 each, on Linux)

#%% Transform back (undo mean-centering and scaling and go to non-log space)
#
# we assumed mean-centered and scaled log-normal distributions for theta
# by our "tricks" above, i.e. μ_y, σ_y fixed
#                    z_logy   ~  Normal(θ, σ^2)
#       (log(y) - m_y) / σ_y  ~  Normal(θ, σ^2)
#       (log(y) - m_y) / σ_y  =  θ + σ * x
# where y is the observation and x the predictor, then
#       log(y) = (θ + σ * x) * σ_y + m_y
#       log(y) = θ*σ_y + σ*σ_y*x  + m_y
#       log(y) = θ*σ_y + m_y + σ*σ_y*x
#       log(y) = θ_trans + σ_trans * x
#       log(y) ~ Normal(θ_trans, σ_trans^2)
#           y  = exp(θ_trans + σ_trans * x)
#           y  ~ LogNormal(θ_trans, σ_trans^2)
#         E[y] = exp(θ_trans + σ_trans^2 / 2)
#
# The same is true for the prior distribution yielding theta:
# ()μ and τ come from hyperdistributions)
#                         θ ~ Normal(μ,τ^2)
#     (θ_trans - m_y) / σ_y = μ + τ * xi
#                   θ_trans = μ*σ_y + m_y + τ*σ_y*xi
#                   θ_trans = μ_trans + τ_trans*xi
#                   θ_trans ~ Normal(μ_trans, τ_trans^2)
#
#                    log(y) = μ_trans + τ_trans*xi + σ_trans*x
# xi and x are both ~ N(0,1)
# a sum of two distributions, i.e. here N(0,τ_trans^2) + N(0,σ_trans^2), are a
# convolution of the two, thus their sum is: N(0,τ_trans^2 + σ_trans^2)
# assuming independence, and therefore
#                    log(y) = μ_trans + sqrt(τ_trans^2 + σ_trans^2)*x
#                    log(y) ~ Normal(μ_trans, (σ_trans^2 + τ_trans^2))
#                      E[y] = exp(μ_trans + (σ_trans^2 + τ_trans^2) / 2)
#                      E[y] = exp(μ_trans + σ_trans^2 / 2 + τ_trans^2 / 2)

# Originals
θ = chain[:,1:J]
μ = chain[:,(J+1)]
σ = chain[:,(J+2)]
τ = chain[:,(J+3)]

# Un-scale and un-mean-centre:
θ_trans = θ .* logStd .+ logMean
μ_trans = μ .* logStd .+ logMean
σ_trans = σ .* logStd
τ_trans = τ .* logStd

# Get into non-log space:
θ_trans_unLog = exp.(θ_trans .+ 0.5 .* repeat(σ_trans,1,J).^2)
μ_trans_unLog = exp.(μ_trans .+ 0.5 .* σ_trans.^2 .+ 0.5 .* τ_trans.^2)

################################################################################
############################# TASKS ############################################

######################## Task A-1 (the Dude)####################################
#%% Mode, Mean, HDIs
dude = θ_trans_unLog[:,4]
ω = mean(find_mode(dude))
μ_bar = mean(dude);
left,right = hdi(dude);

#%% Plotting
histogram(dude, bins=100, normalize=:pdf, label="MCMC", alpha=0.3, linealpha=0.1)   # Comes from StatsPlots now
dPlot = density!(dude,linewidth=3,label="density estimate")
dCurve = filter(!isnan,dPlot.series_list[1].plotattributes[:y])
dTopPoint = maximum(dCurve)
plot!([(left,0),(left,dTopPoint/2)], linewidth=3, color="green", label="HDI",
       annotations = (left, dTopPoint/2, text("$(Int(round(left)))",:green,:bottom)))
plot!([(right,0),(right,dTopPoint/2)],linewidth=3,color="green",label="HDI",
       annotations = (right, dTopPoint/2, text("$(Int(round(right)))",:green,:bottom)))
plot!([(μ_bar,0),(μ_bar,dTopPoint)], linewidth=3, color="blue", label="mean",
       annotations = (μ_bar, 0, text("$(Int(round(μ_bar)))",:blue,:top)))
plot!([(ω,0),(ω,dTopPoint)], linewidth=3, color="red", label="mode",
       annotations = (ω, dTopPoint, text("$(Int(round(ω)))",:red,:bottom)))

Plots.savefig("/home/johhub/Desktop/ABDA/A5/figs/A1-Dude-Slice.pdf")

######################## Task A-2-a-i (expectation new individual) #############
#%% Mode, Mean, HDIs
samplePts = μ_trans_unLog
ω = mean(find_mode(samplePts))
μ_bar = mean(samplePts);
left,right = hdi(samplePts);

#%% Plotting
histogram(samplePts, bins=100, normalize=:pdf, label="MCMC", alpha=0.3, linealpha=0.1)   # Comes from StatsPlots now
dPlot = density!(samplePts,linewidth=3,label="density estimate")
dCurve = filter(!isnan,dPlot.series_list[1].plotattributes[:y])
dTopPoint = maximum(dCurve)
plot!([(left,0),(left,dTopPoint/2)], linewidth=3, color="green", label="HDI",
       annotations = (left, dTopPoint/2, text("$(Int(round(left)))",:green,:bottom)))
plot!([(right,0),(right,dTopPoint/2)],linewidth=3,color="green",label="HDI",
       annotations = (right, dTopPoint/2, text("$(Int(round(right)))",:green,:bottom)))
plot!([(μ_bar,0),(μ_bar,dTopPoint)], linewidth=3, color="blue", label="mean",
       annotations = (μ_bar, 0, text("$(Int(round(μ_bar)))",:blue,:top)))
plot!([(ω,0),(ω,dTopPoint)], linewidth=3, color="red", label="mode",
       annotations = (ω, dTopPoint, text("$(Int(round(ω)))",:red,:bottom)))

Plots.savefig("/home/johhub/Desktop/ABDA/A5/figs/A2-Group-Slice.pdf")

########################  Task A-2-a-ii (new individual prediction)  ############
#%% Predict a reaction time for a single measurement (y)
zlogy_sim = zeros(N)
for i in 1:N
    # 1) Pick a posterior sample from mu, tau, and sigma
    idx = Int(ceil(rand()*N))   # choose a random index which will pick from the simulated posterior
    μ_pick = μ[idx]
    τ_pick = τ[idx]
    σ_pick = σ[idx]
    # 2) Simulate a new theta given these samples, i.e. theta ~ N(μ,τ)
    θ_sim  = μ_pick + randn()*τ_pick
    # 3) Simulate a reaction time measurement (y) given the picked theta and sigma, i.e. zlogy ~ N(θ,σ)
    zlogy_sim[i] = θ_sim + randn()*σ_pick
end

# 4) Calculate y from zlogy (undo mean-centering and scaling, and go to non-log space)
y_sim = exp.(zlogy_sim .* logStd .+ logMean)
samplePts = y_sim

ω = mean(find_mode(samplePts))
μ_bar = mean(samplePts);
med = median(samplePts);
left,right = hdi(samplePts);

histogram(samplePts, bins=100, normalize=:pdf, label="Simulation", alpha=0.3, linealpha=0.1)   # Comes from StatsPlots now
dPlot = density!(samplePts,linewidth=3,label="density estimate")
dCurve = filter(!isnan,dPlot.series_list[1].plotattributes[:y])
dTopPoint = maximum(dCurve)
plot!([(left,0),(left,dTopPoint/2)], linewidth=3, color="green", label="HDI",
       annotations = (left, dTopPoint/2, text("$(Int(round(left)))",:green,:bottom)))
plot!([(right,0),(right,dTopPoint/2)],linewidth=3,color="green",label="HDI",
       annotations = (right, dTopPoint/2, text("$(Int(round(right)))",:green,:bottom)))
plot!([(μ_bar,0),(μ_bar,dTopPoint)], linewidth=3, color="blue", label="mean",
       annotations = (μ_bar, 0, text("$(Int(round(μ_bar)))",:blue,:top)))
plot!([(ω,0),(ω,dTopPoint)], linewidth=3, color="red", label="mode",
       annotations = (ω, dTopPoint, text("$(Int(round(ω)))",:red,:bottom)))
plot!([(med,0),(med,dTopPoint)], linewidth=3, color="purple", label="median",
       annotations = (med, 0, text("$(Int(round(med)))",:purple,:bottom)))

Plots.savefig("/home/johhub/Desktop/ABDA/A5/figs/A2-Pred-Slice.pdf")

########################  Task A-2-b (compare to website statistics)  ##########
#%% Taken from the homepage:
median_web = 273.0
mean_web = 284.0
println("Δ Median = ",median_web - med)
println("Δ Mean = ",mean_web - μ_bar)

########################  Task A-3 #############################################
# Compare hierarchical theta to individual theta using sample means

#%% Logarithmic sample means
θ_means_log = zeros(J)
#for i in 1:I
#    n = length(findall(x -> x==ind[i],ind)) # finding the number of data points for each individual
#    θ_means_log[ind[i]] += logy[i] / n  # sum up the weighted results
#end
# or use this alternative:
for j in 1:J
    θ_means_log[j] = mean(logy[ind .== j])
end

#%% Plot into one figure
# Get ordered indices:
sortIdx = sortperm(θ_means_log)
# Limits for the figure
myXlims = (minimum(θ_trans),maximum(θ_trans))
# Initialise the subplots
StatsPlots.plot(layout=(J, 1),size = (1000, 1500))
# Plot each theta:
for i in 1:(J-1)
    j = sortIdx[i]
    # Sampled thetas and their mean:
    histogram!(θ_trans[:,j], bins=100, normalize=:pdf, legend=false, alpha=0.3, linealpha=0.0,
            ann=(myXlims[1]+.05,4,"ind $j:"),ticks=nothing, yaxis=false, subplot=i, xlims=myXlims)
    vline!([mean(θ_trans[:,j])],linewidth=3, color="black", subplot=i, legend=false)

    # The (log) sample means of the initial data:
    vline!([θ_means_log[j]],linewidth=3, color="red", subplot=i, legend=false)
end

# The last one separately so I can see it in Hydrogen:
j = sortIdx[J]
histogram!(θ_trans[:,j], bins=100, normalize=:pdf, legend=false, alpha=0.3, linealpha=0.0,
            ann=(myXlims[1]+.05,10,"ind $J:"),ticks=nothing, yaxis=false, subplot=J, xlims=myXlims)
vline!([mean(θ_trans[:,j])],linewidth=3, color="black", subplot=J, legend=false)
vline!([θ_means_log[j]],linewidth=3, color="red", subplot=J, legend=false)

Plots.savefig("/home/johhub/Desktop/ABDA/A5/figs/A2-Comp-MLE-Slice-All.pdf")

# --------------- Old Code (please ignore) -------------------------
# theta ~ N(θ_means_log,sigma)

#histogram(θ_trans[:,1], bins=100, normalize=:pdf,label="theta from MCMC",alpha=0.2, linealpha=0.0)
#plot([(θ_means_log[1],0),(θ_means_log[1],.5)], linewidth=3, color="red", label="sample mean",
#       annotations = (med, 0, text("$(Int(round(θ_means_log[1])))",:red,:bottom)))
#plot()
#for i in 1:12
#    histogram!(θ_trans[:,i], bins=100, normalize=:pdf, legend=false, alpha=0.3, linealpha=0.0,
#            inset = (1, bbox(0.0,0.1*i,1,0.1,:bottom,:left)),  # x,y, width, height, origin
#            ticks=nothing, subplot=1+i, bg_inside=nothing)
#vline!([θ_means_log[i]],linewidth=3, color="red", subplot=i+1, legend=false)
#end

# Get a color map:
#curColor = get_color_palette(:auto, plot_color(:white), J)
#for n = 1:Int(ceil(J/5))
#    plot()
#    for j in ((n-1)*5+1):min((n*5),J)
#        histogram!(θ_trans[:,j], bins=100, normalize=:pdf,legend=false,alpha=0.2, linealpha=0.0, color=curColor[j])
#        vline!([θ_means_log[j]], linewidth=3, color=curColor[j])
#    end

#Plots.savefig("/home/johhub/Desktop/ABDA/A5/figs/A2-Comp-MLE-Slice-$n.pdf")

#end


text:
using Torch
using Base.Test

xs = rand(5,3)

th = THArray(xs)

@test xs == th

ext: jl


text:
@testset "4.median-of-two-sorted-arrays.jl" begin
    @test find_median_sorted_arrays([1, 3], [2]) == 2.0
    @test find_median_sorted_arrays([1, 2], [3, 4]) == 2.5
    for _ in 1:40
        nums1 = rand(1:200, 50) |> sort!
        nums2 = rand(1:200, 49) |> sort!
        tt = [nums1; nums2] |> sort!
        @test find_median_sorted_arrays(nums1, nums2) == tt[50]
    end
end


text:
# This example solution is partially derived from the Julia base implementation of rational numbers (MIT license)
# https://github.com/JuliaLang/julia/blob/52bafeb981bac548afd2264edb518d8d86944dca/base/rational.jl
# https://github.com/JuliaLang/julia/blob/52bafeb981bac548afd2264edb518d8d86944dca/LICENSE.md

import Base: +, *, ^, -, /, <, <=, inv, show, promote_rule

struct RationalNumber{T<:Integer} <: Real
    num::T
    den::T

    # reduce to lowest terms
    function RationalNumber{T}(num::Integer, den::Integer) where T<:Integer
        num == den == zero(T) && throw(ArgumentError("zero($T)//zero($T) is not a valid rational number"))
        num2, den2 = sign(den) < 0 ? Base.divgcd(-num, -den) : Base.divgcd(num, den)
        new(num2, den2)
    end
end

RationalNumber(num::T, den::T) where {T<:Integer} = RationalNumber{T}(num, den)

numerator(x::RationalNumber) = x.num
denominator(x::RationalNumber) = x.den

RationalNumber{T}(x::Integer) where {T<:Integer} = RationalNumber{T}(convert(T, x), convert(T, 1))

promote_rule(::Type{RationalNumber{T}}, ::Type{S}) where {T<:Integer,S<:Integer} = RationalNumber{promote_type(T,S)}

show(io::IO, x::RationalNumber) = (show(io, numerator(x)); print(io, "//"); show(io, denominator(x)))

# Helper method
numden(x::RationalNumber) = numerator(x), denominator(x)

# Arithmetics
function +(x::RationalNumber, y::RationalNumber)
    a, b = numden(x)
    c, d = numden(y)

    RationalNumber(a * d + b * c, b * d)
end

-(x::RationalNumber, y::RationalNumber) = x + (-1) * y
-(x::RationalNumber) = -1 * x

function *(x::RationalNumber, y::RationalNumber)
    a, b = numden(x)
    c, d = numden(y)

    RationalNumber(a * c, b * d)
end

inv(x::RationalNumber) = RationalNumber(denominator(x), numerator(x))
/(x::RationalNumber, y::RationalNumber) = x * inv(y)

^(x::Real, y::RationalNumber) = x^(numerator(y) / denominator(y))

# Ordering
<(x::RationalNumber, y::RationalNumber) = x.den == y.den ? x.num < y.num : x.num * y.den < y.num * x.den
<=(x::RationalNumber, y::RationalNumber) = x.den == y.den ? x.num <= y.num : x.num * y.den <= y.num * x.den


text:
using Pkg
haskey(Pkg.installed(), "Coverage") || Pkg.add("Coverage")

using Coverage
Codecov.submit(Codecov.process_folder())


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

using GZip

function search_replace{T<:AbstractString}(files::Vector{T}, src::AbstractString,
    dst::AbstractString; outdir::AbstractString="./converted")
    mkpath(outdir)
    for file in files
        fileroot, fileext = splitext(file)
        if fileext == ".gz" #Gzip format
            replace_text(file, src, dst, outdir, fopen=GZip.open)
        else #assume ascii-compatible
            replace_text(file, src, dst, outdir, fopen=open)
        end
    end
end

function replace_text(file::AbstractString, src::AbstractString, dst::AbstractString,
    outdir::AbstractString; fopen::Function=open)
    text = fopen(readall, file)
    text = replace(text, src, dst)
    fout = fopen(joinpath(outdir, basename(file)), "w")
    write(fout, text)
    close(fout)
end

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
import Base.Random: srand
import RNG: gen_seed

for (w, T, Td) in ((32, UInt32, UInt64), (64, UInt64, UInt128))
    @eval @inline function philox_mulhilo(a::$T, b::$T)
        product = (a % $Td) * (b % $Td)
        hi = (product >> $w) % $T
        lo = product % $T
        (hi, lo)
    end
end

@inline PHILOX_M2x_0(::Type{UInt64}) = 0xD2B74407B1CE6E93
@inline PHILOX_M4x_0(::Type{UInt64}) = 0xD2E7470EE14C6C93
@inline PHILOX_M4x_1(::Type{UInt64}) = 0xCA5A826395121157

@inline PHILOX_M2x_0(::Type{UInt32}) = 0xd256d193
@inline PHILOX_M4x_0(::Type{UInt32}) = 0xD2511F53
@inline PHILOX_M4x_1(::Type{UInt32}) = 0xCD9E8D57

@inline PHILOX_W_0(::Type{UInt64}) = 0x9E3779B97F4A7C15
@inline PHILOX_W_1(::Type{UInt64}) = 0xBB67AE8584CAA73B

@inline PHILOX_W_0(::Type{UInt32}) = 0x9E3779B9
@inline PHILOX_W_1(::Type{UInt32}) = 0xBB67AE85

"""
```julia
Philox2x{T, R} <: R123Generator2x{T}
Philox2x([seed, R])
Philox2x(T[, seed, R])
```

Philox2x is one kind of Philox Counter-Based RNGs. It generates two numbers at a time.

`T` is `UInt32` or `UInt64`(default).

`seed` is an `Integer` which will be automatically converted to `T`.

`R` denotes to the Rounds which must be at least 1 and no more than 16. With 10 rounds (by default), it has a
considerable safety margin over the minimum number of rounds with no known statistical flaws, but still has
excellent performance.
"""
type Philox2x{T<:Union{UInt32, UInt64}, R} <: R123Generator2x{T}
    x1::T
    x2::T
    key::T
    ctr1::T
    ctr2::T
    p::Int
end

function Philox2x{T<:Union{UInt32, UInt64}}(::Type{T}=UInt64, seed::Integer=gen_seed(T), R::Integer=10)
    @assert 1 <= R <= 16
    r = Philox2x{T, Int(R)}(0, 0, 0, 0, 0, 0)
    srand(r, seed)
end
Philox2x(seed::Integer, R::Integer=10) = Philox2x(UInt64, seed, R)

function srand{T<:Union{UInt32, UInt64}}(r::Philox2x{T}, seed::Integer=gen_seed(T))
    r.x1 = r.x2 = 0
    r.key = seed % T
    r.ctr1 = r.ctr2 = 0
    random123_r(r)
    r
end

@inline function philox2x_round{T<:Union{UInt32, UInt64}}(ctr1::T, ctr2::T, key::T)
    hi, lo = philox_mulhilo(PHILOX_M2x_0(T), ctr1)
    hi $ key $ ctr2, lo
end

@inline function philox2x_bumpkey{T<:Union{UInt32, UInt64}}(key::T)
    key + PHILOX_W_0(T)
end

@inline function random123_r{T<:Union{UInt32, UInt64}, R}(r::Philox2x{T, R})
    ctr1, ctr2, key = r.ctr1, r.ctr2, r.key
    if R > 0                               ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 1  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 2  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 3  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 4  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 5  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 6  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 7  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 8  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 9  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 10 key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 11 key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 12 key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 13 key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 14 key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 15 key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    r.x1, r.x2 = ctr1, ctr2
end

"""
```julia
Philox4x{T, R} <: R123Generator4x{T}
Philox4x([seed, R])
Philox4x(T[, seed, R])
```

Philox4x is one kind of Philox Counter-Based RNGs. It generates four numbers at a time.

`T` is `UInt32` or `UInt64`(default).

`seed` is a `Tuple` of two `Integer`s which will both be automatically converted to `T`.

`R` denotes to the Rounds which must be at least 1 and no more than 16. With 10 rounds (by default), it has a
considerable safety margin over the minimum number of rounds with no known statistical flaws, but still has
excellent performance.
"""
type Philox4x{T<:Union{UInt32, UInt64}, R} <: R123Generator4x{T}
    x1::T
    x2::T
    x3::T
    x4::T
    key1::T
    key2::T
    ctr1::T
    ctr2::T
    ctr3::T
    ctr4::T
    p::Int
end

function Philox4x{T<:Union{UInt32, UInt64}}(::Type{T}=UInt64,
        seed::NTuple{2, Integer}=gen_seed(T, 2), R::Integer=10)
    @assert 1 <= R <= 16
    r = Philox4x{T, Int(R)}(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
    srand(r, seed)
end
Philox4x(seed::NTuple{2, Integer}, R::Integer=10) = Philox4x(UInt64, seed, R)

function srand{T<:Union{UInt32, UInt64}}(r::Philox4x{T}, seed::NTuple{2, Integer}=gen_seed(T, 2))
    r.x1 = r.x2 = r.x3 = r.x4 = 0
    r.key1 = seed[1] % T
    r.key2 = seed[2] % T
    r.ctr1 = r.ctr2 = r.ctr3 = r.ctr4 = 0
    random123_r(r)
    r.p = 0
    r
end

@inline function philox4x_round{T<:Union{UInt32, UInt64}}(
        ctr1::T, ctr2::T, ctr3::T, ctr4::T, key1::T, key2::T)
    hi1, lo1 = philox_mulhilo(PHILOX_M4x_0(T), ctr1)
    hi2, lo2 = philox_mulhilo(PHILOX_M4x_1(T), ctr3)
    hi2 $ ctr2 $ key1, lo2, hi1 $ ctr4 $ key2, lo1
end

@inline function philox4x_bumpkey{T<:Union{UInt32, UInt64}}(key1::T, key2::T)
    key1 + PHILOX_W_0(T), key2 + PHILOX_W_1(T)
end

@inline function random123_r{T<:Union{UInt32, UInt64}, R}(r::Philox4x{T, R})
    ctr1, ctr2, ctr3, ctr4 = r.ctr1, r.ctr2, r.ctr3, r.ctr4
    key1, key2 = r.key1, r.key2
    if R > 0
        ctr1, ctr2, ctr3, ctr4 = philox4x_round(ctr1, ctr2, ctr3, ctr4, key1, key2);
    end
    if R > 1
        key1, key2 = philox4x_bumpkey(key1, key2);
        ctr1, ctr2, ctr3, ctr4 = philox4x_round(ctr1, ctr2, ctr3, ctr4, key1, key2);
    end
    if R > 2
        key1, key2 = philox4x_bumpkey(key1, key2);
        ctr1, ctr2, ctr3, ctr4 = philox4x_round(ctr1, ctr2, ctr3, ctr4, key1, key2);
    end
    if R > 3
        key1, key2 = philox4x_bumpkey(key1, key2);
        ctr1, ctr2, ctr3, ctr4 = philox4x_round(ctr1, ctr2, ctr3, ctr4, key1, key2);
    end
    if R > 4
        key1, key2 = philox4x_bumpkey(key1, key2);
        ctr1, ctr2, ctr3, ctr4 = philox4x_round(ctr1, ctr2, ctr3, ctr4, key1, key2);
    end
    if R > 5
        key1, key2 = philox4x_bumpkey(key1, key2);
        ctr1, ctr2, ctr3, ctr4 = philox4x_round(ctr1, ctr2, ctr3, ctr4, key1, key2);
    end
    if R > 6
        key1, key2 = philox4x_bumpkey(key1, key2);
        ctr1, ctr2, ctr3, ctr4 = philox4x_round(ctr1, ctr2, ctr3, ctr4, key1, key2);
    end
    if R > 7
        key1, key2 = philox4x_bumpkey(key1, key2);
        ctr1, ctr2, ctr3, ctr4 = philox4x_round(ctr1, ctr2, ctr3, ctr4, key1, key2);
    end
    if R > 8
        key1, key2 = philox4x_bumpkey(key1, key2);
        ctr1, ctr2, ctr3, ctr4 = philox4x_round(ctr1, ctr2, ctr3, ctr4, key1, key2);
    end
    if R > 9
        key1, key2 = philox4x_bumpkey(key1, key2);
        ctr1, ctr2, ctr3, ctr4 = philox4x_round(ctr1, ctr2, ctr3, ctr4, key1, key2);
    end
    if R > 10
        key1, key2 = philox4x_bumpkey(key1, key2);
        ctr1, ctr2, ctr3, ctr4 = philox4x_round(ctr1, ctr2, ctr3, ctr4, key1, key2);
    end
    if R > 11
        key1, key2 = philox4x_bumpkey(key1, key2);
        ctr1, ctr2, ctr3, ctr4 = philox4x_round(ctr1, ctr2, ctr3, ctr4, key1, key2);
    end
    if R > 12
        key1, key2 = philox4x_bumpkey(key1, key2);
        ctr1, ctr2, ctr3, ctr4 = philox4x_round(ctr1, ctr2, ctr3, ctr4, key1, key2);
    end
    if R > 13
        key1, key2 = philox4x_bumpkey(key1, key2);
        ctr1, ctr2, ctr3, ctr4 = philox4x_round(ctr1, ctr2, ctr3, ctr4, key1, key2);
    end
    if R > 14
        key1, key2 = philox4x_bumpkey(key1, key2);
        ctr1, ctr2, ctr3, ctr4 = philox4x_round(ctr1, ctr2, ctr3, ctr4, key1, key2);
    end
    if R > 15
        key1, key2 = philox4x_bumpkey(key1, key2);
        ctr1, ctr2, ctr3, ctr4 = philox4x_round(ctr1, ctr2, ctr3, ctr4, key1, key2);
    end
    r.x1, r.x2, r.x3, r.x4 = ctr1, ctr2, ctr3, ctr4
end


text:
using Documenter, Algorithms

makedocs(sitename="Algorithms",
pages=[
    "index.md",
    "Maths.md",
    "Algebra.md"
],
doctest=true
)


text:
#
# @description Unit tests for DES block-cipher.
#
# @author Alejandro Claro (alejandro.claro@gmail.com)
#
# Copyright 2017 All rights reserved.
# Use of this source code is governed by a MIT-style license that can be found in the LICENSE file.
#
@testset "DES cipher tests" begin
  const DES_PLAINTEXT = UInt8[ 0x01, 0x23, 0x45, 0x67, 0x89, 0xAB, 0xCD, 0xEF ]
  const DES_R0        = Cryptography.unpack_bits(UInt8[ 0xF0, 0xAA, 0xF0, 0xAA ])

  cipher = DesCipher([ 0x13, 0x34, 0x57, 0x79, 0x9B, 0xBC, 0xDF, 0xF1 ])

  @test block_size(cipher) == 8
  @test key_size(cipher)   == 8

  @test cipher.subkeys == [
    UInt8[ 0x00, 0x00, 0x00, 0x01, 0x01, 0x00, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x01, 0x00, 0x01, 0x01, 0x01, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x00, 0x01, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x01, 0x01, 0x01, 0x00, 0x00, 0x01, 0x00 ],
    UInt8[ 0x00, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x01, 0x01, 0x00, 0x01, 0x00, 0x01, 0x01, 0x01, 0x00, 0x01, 0x01, 0x00, 0x01, 0x01, 0x00, 0x00, 0x01, 0x01, 0x01, 0x00, 0x01, 0x01, 0x00, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x01, 0x00, 0x00, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x01, 0x00, 0x01 ],
    UInt8[ 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x01, 0x00, 0x01, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x01, 0x00, 0x01, 0x01, 0x00, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x01, 0x01, 0x00, 0x00, 0x01 ],
    UInt8[ 0x00, 0x01, 0x01, 0x01, 0x00, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x01, 0x00, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00, 0x01, 0x01, 0x00, 0x01, 0x01, 0x00, 0x01, 0x01, 0x00, 0x01, 0x01, 0x00, 0x00, 0x01, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x00, 0x00, 0x01, 0x01, 0x01, 0x00, 0x01 ],
    UInt8[ 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x01, 0x01, 0x01, 0x00, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00, 0x01, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x00, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x00, 0x00 ],
    UInt8[ 0x00, 0x01, 0x01, 0x00, 0x00, 0x00, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00, 0x00, 0x01, 0x00, 0x01, 0x00, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x01, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x01, 0x00, 0x00, 0x01, 0x00, 0x01, 0x01, 0x01, 0x01 ],
    UInt8[ 0x01, 0x01, 0x01, 0x00, 0x01, 0x01, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x01, 0x00, 0x00, 0x01, 0x00, 0x01, 0x01, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x01, 0x01, 0x00, 0x00, 0x00, 0x01, 0x00, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00 ],
    UInt8[ 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x00, 0x01, 0x00, 0x01, 0x00, 0x00, 0x00, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x01, 0x00, 0x00, 0x01, 0x01, 0x01, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x01 ],
    UInt8[ 0x01, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x01, 0x01, 0x00, 0x01, 0x01, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x01, 0x00, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x01 ],
    UInt8[ 0x01, 0x00, 0x01, 0x01, 0x00, 0x00, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x01, 0x01, 0x00, 0x01, 0x00, 0x00, 0x00, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x01, 0x01, 0x00, 0x00, 0x01, 0x00, 0x00, 0x01, 0x01, 0x01, 0x01 ],
    UInt8[ 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00, 0x00, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x01, 0x00, 0x01, 0x00, 0x00, 0x01, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x01, 0x01, 0x00 ],
    UInt8[ 0x00, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x01, 0x01, 0x00, 0x00, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00, 0x01, 0x01, 0x00, 0x00, 0x01, 0x00, 0x01, 0x00, 0x00, 0x00, 0x01, 0x01, 0x00, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00, 0x00, 0x01 ],
    UInt8[ 0x01, 0x00, 0x00, 0x01, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x00, 0x01, 0x00, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00, 0x00, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x01 ],
    UInt8[ 0x00, 0x01, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x01, 0x01, 0x01, 0x00, 0x01, 0x01, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x01, 0x00, 0x01, 0x01, 0x01, 0x00, 0x00, 0x01, 0x01, 0x01, 0x00, 0x00, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00 ],
    UInt8[ 0x01, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x01, 0x01, 0x00, 0x00, 0x00, 0x01, 0x01, 0x00, 0x01, 0x00, 0x00, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x01, 0x00, 0x01, 0x00 ],
    UInt8[ 0x01, 0x01, 0x00, 0x00, 0x01, 0x00, 0x01, 0x01, 0x00, 0x00, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x01, 0x00, 0x00, 0x00, 0x01, 0x00, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x01, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x01, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00, 0x01 ]
  ]

  @test_throws ArgumentError DesCipher(UInt8[])
  @test_throws ArgumentError DesCipher([ 0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07 ])

  @test_throws ArgumentError encrypt(cipher, UInt8[])
  @test_throws ArgumentError encrypt(cipher, [ 0x05, 0x06 ])

  @test_throws ArgumentError decrypt(cipher, UInt8[])
  @test_throws ArgumentError decrypt(cipher, [ 0x05, 0x06 ])

  @test Cryptography.permutate(Cryptography.unpack_bits(DES_PLAINTEXT), Cryptography.DES_INITIAL_PERMUTATION) == UInt8[
    0x01, 0x01, 0x00, 0x00, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x01, 0x01, 0x00, 0x00, 0x01, 0x01, 0x00, 0x00, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01,
    0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00,
    0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00 ]

  @test Cryptography.permutate(DES_R0, Cryptography.DES_EXPANSION_FUNCTION) == UInt8[
    0x00, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01,
    0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x01, 0x01, 0x01, 0x00, 0x01, 0x00,
    0x00, 0x00, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01 ]

  @test encrypt(cipher, DES_PLAINTEXT) == UInt8[ 0x85, 0xE8, 0x13, 0x54, 0x0F, 0x0A, 0xB4, 0x05 ]

  @test decrypt(cipher, UInt8[ 0x85, 0xE8, 0x13, 0x54, 0x0F, 0x0A, 0xB4, 0x05 ]) == DES_PLAINTEXT
end


text:
mutable struct Entropy{D} <: InfoDist
    data::Vector{Int}
    bs::NTuple{D, Int}
    N::Int
    function Entropy(b::Int, bs::Int...)
        if b ≤ zero(b) || any(b -> b ≤ zero(b), bs)
            throw(ArgumentError("all bases must be greater than zero"))
        end
        new{length(bs) + 1}(zeros(Int, b * prod(bs)), tuple(b, bs...), 0)
    end
end

function Entropy(xs::AbstractArray{Int,3})
    iszero(length(xs)) && throw(ArgumentError("no observations provided"))
    any(b -> b < 1, xs) && throw(ArgumentError("observations must be 1 or greater"))

    dist = Entropy(maximum(xs; dims=(2,3))...)
    observe!(dist, xs)
end
function Entropy(xs::AbstractArray{Int,2})
    iszero(length(xs)) && throw(ArgumentError("no observations provided"))
    any(b -> b < 1, xs) && throw(ArgumentError("observations must be 1 or greater"))

    dist = Entropy(maximum(xs; dims=2)...)
    observe!(dist, xs)
end
function Entropy(xs::AbstractArray{Int,1})
    iszero(length(xs)) && throw(ArgumentError("no observations provided"))
    any(b -> b < 1, xs) && throw(ArgumentError("observations must be 1 or greater"))

    dist = Entropy(maximum(xs))
    observe!(dist, xs)
end

function observe!(dist::Entropy, xs::AbstractArray{Int,3})
    @views for i in 1:size(xs, 3)
        @inbounds observe!(dist, xs[:,:,i])
    end
    dist
end
function observe!(dist::Entropy, xs::AbstractArray{Int,2})
    dist.N += size(xs, 2)
    @views for t in 1:size(xs, 2)
        @inbounds x = index(xs[:,t], dist.bs)
        dist.data[x] += 1
    end
    dist
end
function observe!(dist::Entropy, xs::AbstractArray{Int,1})
    dist.N += length(xs)
    @views for t in 1:length(xs)
        @inbounds dist.data[xs[t]] += 1
    end
    dist
end

function entropy(xs::AbstractArray{Int}, N::Int)
    h = N * log2(N)
    for i in eachindex(xs)
        @inbounds n = xs[i]
        if !iszero(n)
            h -= n * log2(n)
        end
    end
    h / N
end

estimate(dist::Entropy) = entropy(dist.data, dist.N)

entropy!(dist::Entropy, xs::AbstractArray{Int}) = estimate(observe!(dist, xs))

entropy(xs::AbstractArray{Int}) = estimate(Entropy(xs))

function clear!(dist::Entropy)
    fill!(dist.data, 0)
    dist.N = 0
end

Base.length(dist::Entropy) = prod(size(dist))

Base.size(dist::Entropy) = dist.bs

Base.getindex(dist::Entropy, idx...) = getindex(dist.data, idx...)

function entropy(::Type{Kozachenko}, xs::AbstractMatrix{Float64}; metric::Metric=Euclidean())
    D, T = size(xs)
    δ = minimumdistances(xs; metric=metric)
    D*(mean(log.(δ)) + log(2.0)) + sdterm(D) + eulermascheroni(T)
    D * (mean(log.(δ)) + log(2.0)) + sdterm(D) + eulermascheroni(T)
end

eulermascheroni(x) = (x < zero(x)) ? zero(Float64) : digamma(x) - digamma(1)

function sdterm(dim::Int)
    c = (π / 4.0)^(dim / 2)
    if iseven(dim)
        log(c) - sum(log.(2:(dim ÷ 2)))
    else
        c = (c * 2.0^((dim + 1) / 2.)) / sqrt(π)
        log(c) - sum(log.(3:2:dim))
    end
end

function minimumdistances(xs; metric::Metric=Euclidean())
    b = BallTree(xs)
    _, δ = knn(b, xs, 2)
    first.(δ)
end

entropy(xs::AbstractMatrix{Float64}; kwargs...) = entropy(Kozachenko, xs; kwargs...)

entropy(xs::AbstractVector{Float64}; kwargs...) = entropy(reshape(xs, 1, length(xs)))


text:
export TransitionMixture

"""
Description:

    A mixture of discrete transitions:

    f(out, in1, z, A1, A2, ...) = Cat(out|A1*in1)^z_1 * Cat(out|A2*in1)^z_2 * ...

Interfaces:

    1. out
    2. in1
    3. z (switch)
    4. A1 (transition matrix)
    5. A2 (transition matrix)
    ...

Construction:

    TransitionMixture(out, in1, z, A1, A2, ..., id=:some_id)
"""
mutable struct TransitionMixture <: SoftFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function TransitionMixture(out, in1, z, args::Vararg; id=generateId(TransitionMixture))
        n_factors = length(args)
        @ensureVariables(out, in1, z)
        for i=1:n_factors
            @ensureVariables(args[i])
        end
        self = new(id, Array{Interface}(undef, n_factors + 3), Dict{Symbol,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:in1] = self.interfaces[2] = associate!(Interface(self), in1)
        self.i[:z] = self.interfaces[3] = associate!(Interface(self), z)
        for k = 1:n_factors
            self.i[:a*k] = self.interfaces[k + 3] = associate!(Interface(self), args[k])
        end

        return self
    end
end

slug(::Type{TransitionMixture}) = "TM"

# Average energy functional
function averageEnergy(::Type{TransitionMixture},
                       dist_out_in1_switch::ProbabilityDistribution{Multivariate, Contingency},
                       dist_factors::Vararg{ProbabilityDistribution})

    n_factors = length(dist_factors)
    U = 0.0
    for k = 1:n_factors
        U += -tr(dist_out_in1_switch.params[:p][k]'*unsafeLogMean(dist_factors[k]))
    end

    return U
end


text:
# # Five dim sys

# ## Model

module five_dim_sys #jl

using ReachabilityAnalysis

## system matrix
D = [-1.0 -4.0  0.0  0.0  0.0;
      4.0 -1.0  0.0  0.0  0.0;
      0.0  0.0 -3.0  1.0  0.0;
      0.0  0.0 -1.0 -3.0  0.0;
      0.0  0.0  0.0  0.0 -2.0]
P = [ 0.6 -0.1  0.1  0.7 -0.2;
     -0.5  0.7 -0.1 -0.8  0.0;
      0.9 -0.5  0.3 -0.6  0.1;
      0.5 -0.7  0.5  0.6  0.3;
      0.8  0.7  0.6 -0.3  0.2]
A = P * D * inv(P)

## state domain
X = Universe(5)

## input domain
U = Ball2(zeros(5), 0.01)

function model(X0)
    S = @system(x' = Ax + u, x ∈ X, u ∈ U)
    return IVP(S, X0)
end

end #jl

# ## References

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
include("common.jl")
@testset "Parameterized functions" begin
    # Test parameterized delayed logistic equation

    # delayed logistic equation
    f_inplace(du, u, h, p, t) = (du[1] = p[1] * u[1] * (1 - h(p, t-1; idxs = 1)))
    f_scalar(u, h, p, t) = p[1] * u * (1 - h(p, t-1))

    # simple history function (simplicity here comes at a price, see below)
    h(p, t; idxs=nothing) = 0.1

    @testset for inplace in (true, false)
        # define problem
        prob = DDEProblem(inplace ? f_inplace : f_scalar,
                          inplace ? [0.1] : 0.1,
                          h, (0.0, 50.0), [0.3]; constant_lags = [1])

        # solve problem with initial parameter:
        # since h(p, 0) = u0 holds only in the scalar case we have to specify
        # initial_order=1 to indicate that the discontinuity at t = 0 is of first
        # order and to obtain the same results in both cases
        sol1 = solve(prob, MethodOfSteps(Tsit5()); initial_order = 1)
        @test length(sol1) == 26
        @test first(sol1(12)) ≈ 0.884 atol=1e-4
        @test first(sol1[end]) ≈ 1 atol=1e-5

        # solve problem with updated parameter
        prob.p[1] = 1.4
        sol2 = solve(prob, MethodOfSteps(Tsit5()); initial_order = 1)
        @test length(sol2) == 52
        @test first(sol2(12)) ≈ 1.127 atol=4e-4
        @test first(sol2[end]) ≈ 0.995 atol=3e-4
    end
end


text:
@checked_lib librsync_constructors "/Users/jgoldfar/Documents/misc/julia/Librsync/deps/constructors.o"
function construct_rs_stats_t()
  ccall((:construct_rs_stats_t, librsync_constructors), Ptr{Librsync.rs_stats_t}, ())
end


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

@test integrate(x^3*arctan(a+b*x^4), x) == :(-1//8*b^-1*log(1+(a+b*x^4)^2)+1//4*b^-1*(a+b*x^4)*arctan(a+b*x^4))
@test integrate(x^(-1+n)*arctan(a+b*x^n), x) == :(-1//2*b^-1*n^-1*log(1+(a+b*x^n)^2)+b^-1*n^-1*(a+b*x^n)*arctan(a+b*x^n))
@test integrate(x^5*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(1//6*x^6*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+1//36*x^5*(-1*e)^-1//2*(d+e*x^2)^(1/2)+5//96*x*d^2*(-1*e)^-5//2*(d+e*x^2)^(1/2)+5//96*d^3*e^-7//2*(-1*e)^(1/2)*arctanh(x*e^(1/2)*(d+e*x^2)^-1//2)+5//144*d*x^3*(-1*e)^-3//2*(d+e*x^2)^(1/2))
@test integrate(x^3*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(1//4*x^4*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+1//16*x^3*(-1*e)^-1//2*(d+e*x^2)^(1/2)+-3//32*d^2*e^-5//2*(-1*e)^(1/2)*arctanh(x*e^(1/2)*(d+e*x^2)^-1//2)+3//32*d*x*(-1*e)^-3//2*(d+e*x^2)^(1/2))
@test integrate(x*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :((1/2)*x^2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+1//4*x*(-1*e)^-1//2*(d+e*x^2)^(1/2)+1//4*d*e^-3//2*(-1*e)^(1/2)*arctanh(x*e^(1/2)*(d+e*x^2)^-1//2))
@test integrate(x^-1*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)*log(x)+(1/2)*d^(1/2)*e^-1//2*(-1*e)^(1/2)*(1+e*d^-1*x^2)^(1/2)*(d+e*x^2)^-1//2*Polylogarithms.polylog(2,exp(2*arcsinh(x*d^-1//2*e^(1/2))))+-1//2*d^(1/2)*e^-1//2*(-1*e)^(1/2)*(1+e*d^-1*x^2)^(1/2)*(d+e*x^2)^-1//2*arcsinh(x*d^-1//2*e^(1/2))^2+d^(1/2)*e^-1//2*(-1*e)^(1/2)*(1+e*d^-1*x^2)^(1/2)*(d+e*x^2)^-1//2*arcsinh(x*d^-1//2*e^(1/2))*log(1+-1*exp(2*arcsinh(x*d^-1//2*e^(1/2))))+-1*d^(1/2)*e^-1//2*(-1*e)^(1/2)*(1+e*d^-1*x^2)^(1/2)*(d+e*x^2)^-1//2*arcsinh(x*d^-1//2*e^(1/2))*log(x))
@test integrate(x^-3*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(-1//2*x^-2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+-1//2*d^-1*x^-1*(-1*e)^(1/2)*(d+e*x^2)^(1/2))
@test integrate(x^-5*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(-1//4*x^-4*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+-1//6*d^-2*x^-1*(-1*e)^3//2*(d+e*x^2)^(1/2)+-1//12*d^-1*x^-3*(-1*e)^(1/2)*(d+e*x^2)^(1/2))
@test integrate(x^-7*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(-1//6*x^-6*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+-4//45*d^-3*x^-1*(-1*e)^5//2*(d+e*x^2)^(1/2)+-2//45*d^-2*x^-3*(-1*e)^3//2*(d+e*x^2)^(1/2)+-1//30*d^-1*x^-5*(-1*e)^(1/2)*(d+e*x^2)^(1/2))
@test integrate(x^-9*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(-1//8*x^-8*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+-3//140*d^-2*x^-5*(-1*e)^3//2*(d+e*x^2)^(1/2)+-2//35*d^-4*x^-1*(-1*e)^7//2*(d+e*x^2)^(1/2)+-1//35*d^-3*x^-3*(-1*e)^5//2*(d+e*x^2)^(1/2)+-1//56*d^-1*x^-7*(-1*e)^(1/2)*(d+e*x^2)^(1/2))
@test integrate(x^6*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(-1//49*(-1*e)^-7//2*(d+e*x^2)^7//2+1//7*x^7*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+-1//7*d^2*(-1*e)^-7//2*(d+e*x^2)^3//2+1//7*d^3*(-1*e)^-7//2*(d+e*x^2)^(1/2)+3//35*d*(-1*e)^-7//2*(d+e*x^2)^5//2)
@test integrate(x^4*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(1//5*x^5*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+1//25*(-1*e)^-5//2*(d+e*x^2)^5//2+-2//15*d*(-1*e)^-5//2*(d+e*x^2)^3//2+1//5*d^2*(-1*e)^-5//2*(d+e*x^2)^(1/2))
@test integrate(x^2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(-1//9*(-1*e)^-3//2*(d+e*x^2)^3//2+1//3*x^3*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+1//3*d*(-1*e)^-3//2*(d+e*x^2)^(1/2))
@test integrate(arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(x*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+(-1*e)^-1//2*(d+e*x^2)^(1/2))
@test integrate(x^-2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(-1*x^-1*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+-1*d^-1//2*(-1*e)^(1/2)*arctanh(d^-1//2*(d+e*x^2)^(1/2)))
@test integrate(x^-4*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(-1//3*x^-3*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+-1//6*d^-3//2*(-1*e)^3//2*arctanh(d^-1//2*(d+e*x^2)^(1/2))+-1//6*d^-1*x^-2*(-1*e)^(1/2)*(d+e*x^2)^(1/2))
@test integrate(x^-6*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(-1//5*x^-5*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+-3//40*d^-5//2*(-1*e)^5//2*arctanh(d^-1//2*(d+e*x^2)^(1/2))+-3//40*d^-2*x^-2*(-1*e)^3//2*(d+e*x^2)^(1/2)+-1//20*d^-1*x^-4*(-1*e)^(1/2)*(d+e*x^2)^(1/2))
@test integrate(x^9//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(2//11*x^11//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+4//121*x^9//2*(-1*e)^-1//2*(d+e*x^2)^(1/2)+36//847*d*x^5//2*(-1*e)^-3//2*(d+e*x^2)^(1/2)+60//847*d^2*x^(1/2)*(-1*e)^-5//2*(d+e*x^2)^(1/2)+30//847*d^11//4*e^-13//4*(-1*e)^(1/2)*((d^(1/2)+x*e^(1/2))^-2*(d+e*x^2))^(1/2)*(d+e*x^2)^-1//2*(d^(1/2)+x*e^(1/2))*Elliptic.F(2*arctan(d^-1//4*e^1//4*x^(1/2)),1/2))
@test integrate(x^5//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(2//7*x^7//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+4//49*x^5//2*(-1*e)^-1//2*(d+e*x^2)^(1/2)+20//147*d*x^(1/2)*(-1*e)^-3//2*(d+e*x^2)^(1/2)+-10//147*d^7//4*e^-9//4*(-1*e)^(1/2)*((d^(1/2)+x*e^(1/2))^-2*(d+e*x^2))^(1/2)*(d+e*x^2)^-1//2*(d^(1/2)+x*e^(1/2))*Elliptic.F(2*arctan(d^-1//4*e^1//4*x^(1/2)),1/2))
@test integrate(x^(1/2)*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(2//3*x^3//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+4//9*x^(1/2)*(-1*e)^-1//2*(d+e*x^2)^(1/2)+2//9*d^3//4*e^-5//4*(-1*e)^(1/2)*((d^(1/2)+x*e^(1/2))^-2*(d+e*x^2))^(1/2)*(d+e*x^2)^-1//2*(d^(1/2)+x*e^(1/2))*Elliptic.F(2*arctan(d^-1//4*e^1//4*x^(1/2)),1/2))
@test integrate(x^-3//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(-2*x^-1//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+2*d^-1//4*e^-1//4*(-1*e)^(1/2)*((d^(1/2)+x*e^(1/2))^-2*(d+e*x^2))^(1/2)*(d+e*x^2)^-1//2*(d^(1/2)+x*e^(1/2))*Elliptic.F(2*arctan(d^-1//4*e^1//4*x^(1/2)),1/2))
@test integrate(x^-7//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(-2//5*x^-5//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+-4//15*d^-1*x^-3//2*(-1*e)^(1/2)*(d+e*x^2)^(1/2)+-2//15*d^-5//4*e^3//4*(-1*e)^(1/2)*((d^(1/2)+x*e^(1/2))^-2*(d+e*x^2))^(1/2)*(d+e*x^2)^-1//2*(d^(1/2)+x*e^(1/2))*Elliptic.F(2*arctan(d^-1//4*e^1//4*x^(1/2)),1/2))
@test integrate(x^-11//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(-2//9*x^-9//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+-20//189*d^-2*x^-3//2*(-1*e)^3//2*(d+e*x^2)^(1/2)+-4//63*d^-1*x^-7//2*(-1*e)^(1/2)*(d+e*x^2)^(1/2)+10//189*d^-9//4*e^7//4*(-1*e)^(1/2)*((d^(1/2)+x*e^(1/2))^-2*(d+e*x^2))^(1/2)*(d+e*x^2)^-1//2*(d^(1/2)+x*e^(1/2))*Elliptic.F(2*arctan(d^-1//4*e^1//4*x^(1/2)),1/2))
@test integrate(x^-15//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(-2//13*x^-13//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+-60//1001*d^-3*x^-3//2*(-1*e)^5//2*(d+e*x^2)^(1/2)+-36//1001*d^-2*x^-7//2*(-1*e)^3//2*(d+e*x^2)^(1/2)+-4//143*d^-1*x^-11//2*(-1*e)^(1/2)*(d+e*x^2)^(1/2)+-30//1001*d^-13//4*e^11//4*(-1*e)^(1/2)*((d^(1/2)+x*e^(1/2))^-2*(d+e*x^2))^(1/2)*(d+e*x^2)^-1//2*(d^(1/2)+x*e^(1/2))*Elliptic.F(2*arctan(d^-1//4*e^1//4*x^(1/2)),1/2))
@test integrate(x^7//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(2//9*x^9//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+4//81*x^7//2*(-1*e)^-1//2*(d+e*x^2)^(1/2)+28//405*d*x^3//2*(-1*e)^-3//2*(d+e*x^2)^(1/2)+-28//135*d^2*e^-5//2*x^(1/2)*(-1*e)^(1/2)*(d+e*x^2)^(1/2)*(d^(1/2)+x*e^(1/2))^-1+-14//135*d^9//4*e^-11//4*(-1*e)^(1/2)*((d^(1/2)+x*e^(1/2))^-2*(d+e*x^2))^(1/2)*(d+e*x^2)^-1//2*(d^(1/2)+x*e^(1/2))*Elliptic.F(2*arctan(d^-1//4*e^1//4*x^(1/2)),1/2)+28//135*d^9//4*e^-11//4*(-1*e)^(1/2)*((d^(1/2)+x*e^(1/2))^-2*(d+e*x^2))^(1/2)*(d+e*x^2)^-1//2*(d^(1/2)+x*e^(1/2))*Elliptic.E(2*arctan(d^-1//4*e^1//4*x^(1/2)),1/2))
@test integrate(x^3//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(2//5*x^5//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+4//25*x^3//2*(-1*e)^-1//2*(d+e*x^2)^(1/2)+12//25*d*e^-3//2*x^(1/2)*(-1*e)^(1/2)*(d+e*x^2)^(1/2)*(d^(1/2)+x*e^(1/2))^-1+-12//25*d^5//4*e^-7//4*(-1*e)^(1/2)*((d^(1/2)+x*e^(1/2))^-2*(d+e*x^2))^(1/2)*(d+e*x^2)^-1//2*(d^(1/2)+x*e^(1/2))*Elliptic.E(2*arctan(d^-1//4*e^1//4*x^(1/2)),1/2)+6//25*d^5//4*e^-7//4*(-1*e)^(1/2)*((d^(1/2)+x*e^(1/2))^-2*(d+e*x^2))^(1/2)*(d+e*x^2)^-1//2*(d^(1/2)+x*e^(1/2))*Elliptic.F(2*arctan(d^-1//4*e^1//4*x^(1/2)),1/2))
@test integrate(x^-1//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(2*x^(1/2)*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+-4*e^-1//2*x^(1/2)*(-1*e)^(1/2)*(d+e*x^2)^(1/2)*(d^(1/2)+x*e^(1/2))^-1+-2*d^1//4*e^-3//4*(-1*e)^(1/2)*((d^(1/2)+x*e^(1/2))^-2*(d+e*x^2))^(1/2)*(d+e*x^2)^-1//2*(d^(1/2)+x*e^(1/2))*Elliptic.F(2*arctan(d^-1//4*e^1//4*x^(1/2)),1/2)+4*d^1//4*e^-3//4*(-1*e)^(1/2)*((d^(1/2)+x*e^(1/2))^-2*(d+e*x^2))^(1/2)*(d+e*x^2)^-1//2*(d^(1/2)+x*e^(1/2))*Elliptic.E(2*arctan(d^-1//4*e^1//4*x^(1/2)),1/2))
@test integrate(x^-5//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(-2//3*x^-3//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+-4//3*d^-1*x^-1//2*(-1*e)^(1/2)*(d+e*x^2)^(1/2)+4//3*d^-1*x^(1/2)*(-1*e^2)^(1/2)*(d+e*x^2)^(1/2)*(d^(1/2)+x*e^(1/2))^-1+-4//3*d^-3//4*e^1//4*(-1*e)^(1/2)*((d^(1/2)+x*e^(1/2))^-2*(d+e*x^2))^(1/2)*(d+e*x^2)^-1//2*(d^(1/2)+x*e^(1/2))*Elliptic.E(2*arctan(d^-1//4*e^1//4*x^(1/2)),1/2)+2//3*d^-3//4*e^1//4*(-1*e)^(1/2)*((d^(1/2)+x*e^(1/2))^-2*(d+e*x^2))^(1/2)*(d+e*x^2)^-1//2*(d^(1/2)+x*e^(1/2))*Elliptic.F(2*arctan(d^-1//4*e^1//4*x^(1/2)),1/2))
@test integrate(x^-9//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2), x) == :(-2//7*x^-7//2*arctan(x*(-1*e)^(1/2)*(d+e*x^2)^-1//2)+-12//35*d^-2*x^-1//2*(-1*e)^3//2*(d+e*x^2)^(1/2)+-4//35*d^-1*x^-5//2*(-1*e)^(1/2)*(d+e*x^2)^(1/2)+-12//35*d^-2*e^3//2*x^(1/2)*(-1*e)^(1/2)*(d+e*x^2)^(1/2)*(d^(1/2)+x*e^(1/2))^-1+-6//35*d^-7//4*e^5//4*(-1*e)^(1/2)*((d^(1/2)+x*e^(1/2))^-2*(d+e*x^2))^(1/2)*(d+e*x^2)^-1//2*(d^(1/2)+x*e^(1/2))*Elliptic.F(2*arctan(d^-1//4*e^1//4*x^(1/2)),1/2)+12//35*d^-7//4*e^5//4*(-1*e)^(1/2)*((d^(1/2)+x*e^(1/2))^-2*(d+e*x^2))^(1/2)*(d+e*x^2)^-1//2*(d^(1/2)+x*e^(1/2))*Elliptic.E(2*arctan(d^-1//4*e^1//4*x^(1/2)),1/2))
@test integrate(x^-2*arctan(1+x+x^2), x) == :((1/2)*arctan(1+x)+(1/2)*log(x)+-1//2*log(1+x^2)+1//4*log(2+x^2+2x)+-1*x^-1*arctan(1+x+x^2))
@test integrate((1+-1*c^2*x^2)^-1*(a+b*arctan((1+c*x)^-1//2*(1+-1*c*x)^(1/2)))^3, x) == :(-2*c^-1*(a+b*arctan((1+c*x)^-1//2*(1+-1*c*x)^(1/2)))^3*arctanh(1+-2*(1+im*(1+c*x)^-1//2*(1+-1*c*x)^(1/2))^-1)+-3//2*b^2*c^-1*(a+b*arctan((1+c*x)^-1//2*(1+-1*c*x)^(1/2)))*Polylogarithms.polylog(3,-1+2*(1+im*(1+c*x)^-1//2*(1+-1*c*x)^(1/2))^-1)+-3//4*im*b^3*c^-1*Polylogarithms.polylog(4,1+-2*(1+im*(1+c*x)^-1//2*(1+-1*c*x)^(1/2))^-1)+3//2*b^2*c^-1*(a+b*arctan((1+c*x)^-1//2*(1+-1*c*x)^(1/2)))*Polylogarithms.polylog(3,1+-2*(1+im*(1+c*x)^-1//2*(1+-1*c*x)^(1/2))^-1)+3//4*im*b^3*c^-1*Polylogarithms.polylog(4,-1+2*(1+im*(1+c*x)^-1//2*(1+-1*c*x)^(1/2))^-1)+-3//2*im*b*c^-1*(a+b*arctan((1+c*x)^-1//2*(1+-1*c*x)^(1/2)))^2*Polylogarithms.polylog(2,-1+2*(1+im*(1+c*x)^-1//2*(1+-1*c*x)^(1/2))^-1)+3//2*im*b*c^-1*(a+b*arctan((1+c*x)^-1//2*(1+-1*c*x)^(1/2)))^2*Polylogarithms.polylog(2,1+-2*(1+im*(1+c*x)^-1//2*(1+-1*c*x)^(1/2))^-1))
@test integrate((1+-1*c^2*x^2)^-1*(a+b*arctan((1+c*x)^-1//2*(1+-1*c*x)^(1/2)))^2, x) == :((1/2)*b^2*c^-1*Polylogarithms.polylog(3,1+-2*(1+im*(1+c*x)^-1//2*(1+-1*c*x)^(1/2))^-1)+-2*c^-1*(a+b*arctan((1+c*x)^-1//2*(1+-1*c*x)^(1/2)))^2*arctanh(1+-2*(1+im*(1+c*x)^-1//2*(1+-1*c*x)^(1/2))^-1)+-1//2*b^2*c^-1*Polylogarithms.polylog(3,-1+2*(1+im*(1+c*x)^-1//2*(1+-1*c*x)^(1/2))^-1)+im*b*c^-1*(a+b*arctan((1+c*x)^-1//2*(1+-1*c*x)^(1/2)))*Polylogarithms.polylog(2,1+-2*(1+im*(1+c*x)^-1//2*(1+-1*c*x)^(1/2))^-1)+-1*im*b*c^-1*(a+b*arctan((1+c*x)^-1//2*(1+-1*c*x)^(1/2)))*Polylogarithms.polylog(2,-1+2*(1+im*(1+c*x)^-1//2*(1+-1*c*x)^(1/2))^-1))
@test integrate((1+-1*c^2*x^2)^-1*(a+b*arctan((1+c*x)^-1//2*(1+-1*c*x)^(1/2))), x) == :(-1*a*c^-1*log((1+c*x)^-1//2*(1+-1*c*x)^(1/2))+(1/2)*im*b*c^-1*Polylogarithms.polylog(2,im*(1+c*x)^-1//2*(1+-1*c*x)^(1/2))+-1//2*im*b*c^-1*Polylogarithms.polylog(2,-1*im*(1+c*x)^-1//2*(1+-1*c*x)^(1/2)))
@test integrate(x^m*arctan(tan(a+b*x)), x) == :(x^(1+m)*(1+m)^-1*arctan(tan(a+b*x))+-1*b*x^(2+m)*(2+m^2+3m)^-1)
@test integrate(x^2*arctan(tan(a+b*x)), x) == :(-1//12*b*x^4+1//3*x^3*arctan(tan(a+b*x)))
@test integrate(x*arctan(tan(a+b*x)), x) == :((1/2)*x^2*arctan(tan(a+b*x))+-1//6*b*x^3)
@test integrate(arctan(tan(a+b*x)), x) == :((1/2)*b^-1*arctan(tan(a+b*x))^2)
@test integrate(x^-1*arctan(tan(a+b*x)), x) == :(b*x+-1*(-1*arctan(tan(a+b*x))+b*x)*log(x))
@test integrate(x^m*arctan(cot(a+b*x)), x) == :(b*x^(2+m)*(2+m^2+3m)^-1+x^(1+m)*(1+m)^-1*arctan(cot(a+b*x)))
@test integrate(x^2*arctan(cot(a+b*x)), x) == :(1//3*x^3*arctan(cot(a+b*x))+1//12*b*x^4)
@test integrate(x*arctan(cot(a+b*x)), x) == :((1/2)*x^2*arctan(cot(a+b*x))+1//6*b*x^3)
@test integrate(arctan(cot(a+b*x)), x) == :(-1//2*b^-1*arctan(cot(a+b*x))^2)
@test integrate(x^-1*arctan(cot(a+b*x)), x) == :((b*x+arctan(cot(a+b*x)))*log(x)+-1*b*x)
@test integrate(arctan(tan(a+b*x)), x) == :((1/2)*b^-1*arctan(tan(a+b*x))^2)
@test integrate(x^2*arctan(c+d*tan(a+b*x)), x) == :(-1//8*b^-3*Polylogarithms.polylog(4,-1*(1+-1d+im*c)^-1*(1+d+im*c)*exp(2*im*a+2*im*b*x))+1//3*x^3*arctan(c+d*tan(a+b*x))+1//8*b^-3*Polylogarithms.polylog(4,-1*(c+im*(1+d))^-1*(c+im*(1+-1d))*exp(2*im*a+2*im*b*x))+-1//4*b^-1*x^2*Polylogarithms.polylog(2,-1*(c+im*(1+d))^-1*(c+im*(1+-1d))*exp(2*im*a+2*im*b*x))+-1//6*im*x^3*log(1+(c+im*(1+d))^-1*(c+im*(1+-1d))*exp(2*im*a+2*im*b*x))+1//4*b^-1*x^2*Polylogarithms.polylog(2,-1*(1+-1d+im*c)^-1*(1+d+im*c)*exp(2*im*a+2*im*b*x))+1//6*im*x^3*log(1+(1+-1d+im*c)^-1*(1+d+im*c)*exp(2*im*a+2*im*b*x))+-1//4*im*x*b^-2*Polylogarithms.polylog(3,-1*(c+im*(1+d))^-1*(c+im*(1+-1d))*exp(2*im*a+2*im*b*x))+1//4*im*x*b^-2*Polylogarithms.polylog(3,-1*(1+-1d+im*c)^-1*(1+d+im*c)*exp(2*im*a+2*im*b*x)))
@test integrate(x*arctan(c+d*tan(a+b*x)), x) == :((1/2)*x^2*arctan(c+d*tan(a+b*x))+-1//4*im*x^2*log(1+(c+im*(1+d))^-1*(c+im*(1+-1d))*exp(2*im*a+2*im*b*x))+-1//4*x*b^-1*Polylogarithms.polylog(2,-1*(c+im*(1+d))^-1*(c+im*(1+-1d))*exp(2*im*a+2*im*b*x))+-1//8*im*b^-2*Polylogarithms.polylog(3,-1*(c+im*(1+d))^-1*(c+im*(1+-1d))*exp(2*im*a+2*im*b*x))+1//4*im*x^2*log(1+(1+-1d+im*c)^-1*(1+d+im*c)*exp(2*im*a+2*im*b*x))+1//4*x*b^-1*Polylogarithms.polylog(2,-1*(1+-1d+im*c)^-1*(1+d+im*c)*exp(2*im*a+2*im*b*x))+1//8*im*b^-2*Polylogarithms.polylog(3,-1*(1+-1d+im*c)^-1*(1+d+im*c)*exp(2*im*a+2*im*b*x)))
@test integrate(arctan(c+d*tan(a+b*x)), x) == :(x*arctan(c+d*tan(a+b*x))+-1//4*b^-1*Polylogarithms.polylog(2,-1*(c+im*(1+d))^-1*(c+im*(1+-1d))*exp(2*im*a+2*im*b*x))+1//4*b^-1*Polylogarithms.polylog(2,-1*(1+-1d+im*c)^-1*(1+d+im*c)*exp(2*im*a+2*im*b*x))+(1/2)*im*x*log(1+(1+-1d+im*c)^-1*(1+d+im*c)*exp(2*im*a+2*im*b*x))+-1//2*im*x*log(1+(c+im*(1+d))^-1*(c+im*(1+-1d))*exp(2*im*a+2*im*b*x)))
@test integrate(x^2*arctan(c+(1+im*c)*tan(a+b*x)), x) == :(-1//12*b*x^4+1//3*x^3*arctan(c+(1+im*c)*tan(a+b*x))+1//8*b^-3*Polylogarithms.polylog(4,im*c*exp(2*im*a+2*im*b*x))+-1//4*b^-1*x^2*Polylogarithms.polylog(2,im*c*exp(2*im*a+2*im*b*x))+-1//6*im*x^3*log(1+-1*im*c*exp(2*im*a+2*im*b*x))+-1//4*im*x*b^-2*Polylogarithms.polylog(3,im*c*exp(2*im*a+2*im*b*x)))
@test integrate(x*arctan(c+(1+im*c)*tan(a+b*x)), x) == :((1/2)*x^2*arctan(c+(1+im*c)*tan(a+b*x))+-1//6*b*x^3+-1//4*im*x^2*log(1+-1*im*c*exp(2*im*a+2*im*b*x))+-1//4*x*b^-1*Polylogarithms.polylog(2,im*c*exp(2*im*a+2*im*b*x))+-1//8*im*b^-2*Polylogarithms.polylog(3,im*c*exp(2*im*a+2*im*b*x)))
@test integrate(arctan(c+(1+im*c)*tan(a+b*x)), x) == :(x*arctan(c+(1+im*c)*tan(a+b*x))+-1//2*b*x^2+-1//4*b^-1*Polylogarithms.polylog(2,im*c*exp(2*im*a+2*im*b*x))+-1//2*im*x*log(1+-1*im*c*exp(2*im*a+2*im*b*x)))
@test integrate(x^2*arctan(c+(-1+im*c)*tan(a+b*x)), x) == :(-1//8*b^-3*Polylogarithms.polylog(4,-1*im*c*exp(2*im*a+2*im*b*x))+1//3*x^3*arctan(c+-1*(1+-1*im*c)*tan(a+b*x))+1//12*b*x^4+1//4*b^-1*x^2*Polylogarithms.polylog(2,-1*im*c*exp(2*im*a+2*im*b*x))+1//6*im*x^3*log(1+im*c*exp(2*im*a+2*im*b*x))+1//4*im*x*b^-2*Polylogarithms.polylog(3,-1*im*c*exp(2*im*a+2*im*b*x)))
@test integrate(x*arctan(c+(-1+im*c)*tan(a+b*x)), x) == :((1/2)*x^2*arctan(c+-1*(1+-1*im*c)*tan(a+b*x))+1//6*b*x^3+1//4*im*x^2*log(1+im*c*exp(2*im*a+2*im*b*x))+1//4*x*b^-1*Polylogarithms.polylog(2,-1*im*c*exp(2*im*a+2*im*b*x))+1//8*im*b^-2*Polylogarithms.polylog(3,-1*im*c*exp(2*im*a+2*im*b*x)))
@test integrate(arctan(c+(-1+im*c)*tan(a+b*x)), x) == :(x*arctan(c+-1*(1+-1*im*c)*tan(a+b*x))+(1/2)*b*x^2+1//4*b^-1*Polylogarithms.polylog(2,-1*im*c*exp(2*im*a+2*im*b*x))+(1/2)*im*x*log(1+im*c*exp(2*im*a+2*im*b*x)))
@test integrate(arctan(cot(a+b*x)), x) == :(-1//2*b^-1*arctan(cot(a+b*x))^2)
@test integrate(x^2*arctan(c+d*cot(a+b*x)), x) == :(-1//8*b^-3*Polylogarithms.polylog(4,(1+d+im*c)^-1*(1+-1d+im*c)*exp(2*im*a+2*im*b*x))+1//3*x^3*arctan(c+d*cot(a+b*x))+1//8*b^-3*Polylogarithms.polylog(4,(c+im*(1+-1d))^-1*(c+im*(1+d))*exp(2*im*a+2*im*b*x))+-1//4*b^-1*x^2*Polylogarithms.polylog(2,(c+im*(1+-1d))^-1*(c+im*(1+d))*exp(2*im*a+2*im*b*x))+-1//6*im*x^3*log(1+-1*(c+im*(1+-1d))^-1*(c+im*(1+d))*exp(2*im*a+2*im*b*x))+1//4*b^-1*x^2*Polylogarithms.polylog(2,(1+d+im*c)^-1*(1+-1d+im*c)*exp(2*im*a+2*im*b*x))+1//6*im*x^3*log(1+-1*(1+d+im*c)^-1*(1+-1d+im*c)*exp(2*im*a+2*im*b*x))+-1//4*im*x*b^-2*Polylogarithms.polylog(3,(c+im*(1+-1d))^-1*(c+im*(1+d))*exp(2*im*a+2*im*b*x))+1//4*im*x*b^-2*Polylogarithms.polylog(3,(1+d+im*c)^-1*(1+-1d+im*c)*exp(2*im*a+2*im*b*x)))
@test integrate(x*arctan(c+d*cot(a+b*x)), x) == :((1/2)*x^2*arctan(c+d*cot(a+b*x))+-1//4*im*x^2*log(1+-1*(c+im*(1+-1d))^-1*(c+im*(1+d))*exp(2*im*a+2*im*b*x))+-1//4*x*b^-1*Polylogarithms.polylog(2,(c+im*(1+-1d))^-1*(c+im*(1+d))*exp(2*im*a+2*im*b*x))+-1//8*im*b^-2*Polylogarithms.polylog(3,(c+im*(1+-1d))^-1*(c+im*(1+d))*exp(2*im*a+2*im*b*x))+1//4*im*x^2*log(1+-1*(1+d+im*c)^-1*(1+-1d+im*c)*exp(2*im*a+2*im*b*x))+1//4*x*b^-1*Polylogarithms.polylog(2,(1+d+im*c)^-1*(1+-1d+im*c)*exp(2*im*a+2*im*b*x))+1//8*im*b^-2*Polylogarithms.polylog(3,(1+d+im*c)^-1*(1+-1d+im*c)*exp(2*im*a+2*im*b*x)))
@test integrate(arctan(c+d*cot(a+b*x)), x) == :(x*arctan(c+d*cot(a+b*x))+-1//4*b^-1*Polylogarithms.polylog(2,(c+im*(1+-1d))^-1*(c+im*(1+d))*exp(2*im*a+2*im*b*x))+1//4*b^-1*Polylogarithms.polylog(2,(1+d+im*c)^-1*(1+-1d+im*c)*exp(2*im*a+2*im*b*x))+(1/2)*im*x*log(1+-1*(1+d+im*c)^-1*(1+-1d+im*c)*exp(2*im*a+2*im*b*x))+-1//2*im*x*log(1+-1*(c+im*(1+-1d))^-1*(c+im*(1+d))*exp(2*im*a+2*im*b*x)))
@test integrate(x^2*arctan(c+(1+-1*im*c)*cot(a+b*x)), x) == :(-1//8*b^-3*Polylogarithms.polylog(4,im*c*exp(2*im*a+2*im*b*x))+1//3*x^3*arctan(c+(1+-1*im*c)*cot(a+b*x))+1//12*b*x^4+1//4*b^-1*x^2*Polylogarithms.polylog(2,im*c*exp(2*im*a+2*im*b*x))+1//6*im*x^3*log(1+-1*im*c*exp(2*im*a+2*im*b*x))+1//4*im*x*b^-2*Polylogarithms.polylog(3,im*c*exp(2*im*a+2*im*b*x)))
@test integrate(x*arctan(c+(1+-1*im*c)*cot(a+b*x)), x) == :((1/2)*x^2*arctan(c+(1+-1*im*c)*cot(a+b*x))+1//6*b*x^3+1//4*im*x^2*log(1+-1*im*c*exp(2*im*a+2*im*b*x))+1//4*x*b^-1*Polylogarithms.polylog(2,im*c*exp(2*im*a+2*im*b*x))+1//8*im*b^-2*Polylogarithms.polylog(3,im*c*exp(2*im*a+2*im*b*x)))
@test integrate(arctan(c+(1+-1*im*c)*cot(a+b*x)), x) == :(x*arctan(c+(1+-1*im*c)*cot(a+b*x))+(1/2)*b*x^2+1//4*b^-1*Polylogarithms.polylog(2,im*c*exp(2*im*a+2*im*b*x))+(1/2)*im*x*log(1+-1*im*c*exp(2*im*a+2*im*b*x)))
@test integrate(x^2*arctan(c+(-1+-1*im*c)*cot(a+b*x)), x) == :(-1//12*b*x^4+1//3*x^3*arctan(c+-1*(1+im*c)*cot(a+b*x))+1//8*b^-3*Polylogarithms.polylog(4,-1*im*c*exp(2*im*a+2*im*b*x))+-1//4*b^-1*x^2*Polylogarithms.polylog(2,-1*im*c*exp(2*im*a+2*im*b*x))+-1//6*im*x^3*log(1+im*c*exp(2*im*a+2*im*b*x))+-1//4*im*x*b^-2*Polylogarithms.polylog(3,-1*im*c*exp(2*im*a+2*im*b*x)))
@test integrate(x*arctan(c+(-1+-1*im*c)*cot(a+b*x)), x) == :((1/2)*x^2*arctan(c+-1*(1+im*c)*cot(a+b*x))+-1//6*b*x^3+-1//4*im*x^2*log(1+im*c*exp(2*im*a+2*im*b*x))+-1//4*x*b^-1*Polylogarithms.polylog(2,-1*im*c*exp(2*im*a+2*im*b*x))+-1//8*im*b^-2*Polylogarithms.polylog(3,-1*im*c*exp(2*im*a+2*im*b*x)))
@test integrate(arctan(c+(-1+-1*im*c)*cot(a+b*x)), x) == :(x*arctan(c+-1*(1+im*c)*cot(a+b*x))+-1//2*b*x^2+-1//4*b^-1*Polylogarithms.polylog(2,-1*im*c*exp(2*im*a+2*im*b*x))+-1//2*im*x*log(1+im*c*exp(2*im*a+2*im*b*x)))
@test integrate(arctan(sinh(x)), x) == :(im*Polylogarithms.polylog(2,-1*im*exp(x))+x*arctan(sinh(x))+-1*im*Polylogarithms.polylog(2,im*exp(x))+-2*x*arctan(exp(x)))
@test integrate(x*arctan(sinh(x)), x) == :(im*Polylogarithms.polylog(3,im*exp(x))+(1/2)*x^2*arctan(sinh(x))+-1*im*Polylogarithms.polylog(3,-1*im*exp(x))+-1*x^2*arctan(exp(x))+im*x*Polylogarithms.polylog(2,-1*im*exp(x))+-1*im*x*Polylogarithms.polylog(2,im*exp(x)))
@test integrate(x^2*arctan(sinh(x)), x) == :(-2*im*Polylogarithms.polylog(4,im*exp(x))+2*im*Polylogarithms.polylog(4,-1*im*exp(x))+-2//3*x^3*arctan(exp(x))+1//3*x^3*arctan(sinh(x))+im*x^2*Polylogarithms.polylog(2,-1*im*exp(x))+-1*im*x^2*Polylogarithms.polylog(2,im*exp(x))+-2*im*x*Polylogarithms.polylog(3,-1*im*exp(x))+2*im*x*Polylogarithms.polylog(3,im*exp(x)))
@test integrate((e+f*x)^3*arctan(tanh(a+b*x)), x) == :(-1//4*f^-1*(e+f*x)^4*arctan(exp(2a+2*b*x))+1//4*f^-1*(e+f*x)^4*arctan(tanh(a+b*x))+-3//16*im*b^-4*f^3*Polylogarithms.polylog(5,-1*im*exp(2a+2*b*x))+-1//4*im*b^-1*(e+f*x)^3*Polylogarithms.polylog(2,im*exp(2a+2*b*x))+1//4*im*b^-1*(e+f*x)^3*Polylogarithms.polylog(2,-1*im*exp(2a+2*b*x))+3//16*im*b^-4*f^3*Polylogarithms.polylog(5,im*exp(2a+2*b*x))+-3//8*im*f*b^-2*(e+f*x)^2*Polylogarithms.polylog(3,-1*im*exp(2a+2*b*x))+-3//8*im*b^-3*f^2*(e+f*x)*Polylogarithms.polylog(4,im*exp(2a+2*b*x))+3//8*im*f*b^-2*(e+f*x)^2*Polylogarithms.polylog(3,im*exp(2a+2*b*x))+3//8*im*b^-3*f^2*(e+f*x)*Polylogarithms.polylog(4,-1*im*exp(2a+2*b*x)))
@test integrate((e+f*x)^2*arctan(tanh(a+b*x)), x) == :(-1//3*f^-1*(e+f*x)^3*arctan(exp(2a+2*b*x))+1//3*f^-1*(e+f*x)^3*arctan(tanh(a+b*x))+-1//4*im*b^-1*(e+f*x)^2*Polylogarithms.polylog(2,im*exp(2a+2*b*x))+-1//8*im*b^-3*f^2*Polylogarithms.polylog(4,im*exp(2a+2*b*x))+1//4*im*b^-1*(e+f*x)^2*Polylogarithms.polylog(2,-1*im*exp(2a+2*b*x))+1//8*im*b^-3*f^2*Polylogarithms.polylog(4,-1*im*exp(2a+2*b*x))+-1//4*im*f*b^-2*(e+f*x)*Polylogarithms.polylog(3,-1*im*exp(2a+2*b*x))+1//4*im*f*b^-2*(e+f*x)*Polylogarithms.polylog(3,im*exp(2a+2*b*x)))
@test integrate((e+f*x)*arctan(tanh(a+b*x)), x) == :((1/2)*f^-1*(e+f*x)^2*arctan(tanh(a+b*x))+-1//2*f^-1*(e+f*x)^2*arctan(exp(2a+2*b*x))+-1//4*im*b^-1*(e+f*x)*Polylogarithms.polylog(2,im*exp(2a+2*b*x))+-1//8*im*f*b^-2*Polylogarithms.polylog(3,-1*im*exp(2a+2*b*x))+1//4*im*b^-1*(e+f*x)*Polylogarithms.polylog(2,-1*im*exp(2a+2*b*x))+1//8*im*f*b^-2*Polylogarithms.polylog(3,im*exp(2a+2*b*x)))
@test integrate(arctan(tanh(a+b*x)), x) == :(x*arctan(tanh(a+b*x))+-1*x*arctan(exp(2a+2*b*x))+-1//4*im*b^-1*Polylogarithms.polylog(2,im*exp(2a+2*b*x))+1//4*im*b^-1*Polylogarithms.polylog(2,-1*im*exp(2a+2*b*x)))
@test integrate(x^2*arctan(c+d*tanh(a+b*x)), x) == :(1//3*x^3*arctan(c+d*tanh(a+b*x))+-1//6*im*x^3*log(1+(im+c+-1d)^-1*(im+c+d)*exp(2a+2*b*x))+-1//8*im*b^-3*Polylogarithms.polylog(4,-1*(im+c+-1d)^-1*(im+c+d)*exp(2a+2*b*x))+1//6*im*x^3*log(1+(im+d+-1c)^-1*(im+-1c+-1d)*exp(2a+2*b*x))+1//8*im*b^-3*Polylogarithms.polylog(4,-1*(im+d+-1c)^-1*(im+-1c+-1d)*exp(2a+2*b*x))+-1//4*im*x*b^-2*Polylogarithms.polylog(3,-1*(im+d+-1c)^-1*(im+-1c+-1d)*exp(2a+2*b*x))+-1//4*im*b^-1*x^2*Polylogarithms.polylog(2,-1*(im+c+-1d)^-1*(im+c+d)*exp(2a+2*b*x))+1//4*im*x*b^-2*Polylogarithms.polylog(3,-1*(im+c+-1d)^-1*(im+c+d)*exp(2a+2*b*x))+1//4*im*b^-1*x^2*Polylogarithms.polylog(2,-1*(im+d+-1c)^-1*(im+-1c+-1d)*exp(2a+2*b*x)))
@test integrate(x*arctan(c+d*tanh(a+b*x)), x) == :((1/2)*x^2*arctan(c+d*tanh(a+b*x))+-1//4*im*x^2*log(1+(im+c+-1d)^-1*(im+c+d)*exp(2a+2*b*x))+-1//8*im*b^-2*Polylogarithms.polylog(3,-1*(im+d+-1c)^-1*(im+-1c+-1d)*exp(2a+2*b*x))+1//4*im*x^2*log(1+(im+d+-1c)^-1*(im+-1c+-1d)*exp(2a+2*b*x))+1//8*im*b^-2*Polylogarithms.polylog(3,-1*(im+c+-1d)^-1*(im+c+d)*exp(2a+2*b*x))+-1//4*im*x*b^-1*Polylogarithms.polylog(2,-1*(im+c+-1d)^-1*(im+c+d)*exp(2a+2*b*x))+1//4*im*x*b^-1*Polylogarithms.polylog(2,-1*(im+d+-1c)^-1*(im+-1c+-1d)*exp(2a+2*b*x)))
@test integrate(arctan(c+d*tanh(a+b*x)), x) == :(x*arctan(c+d*tanh(a+b*x))+(1/2)*im*x*log(1+(im+d+-1c)^-1*(im+-1c+-1d)*exp(2a+2*b*x))+-1//2*im*x*log(1+(im+c+-1d)^-1*(im+c+d)*exp(2a+2*b*x))+-1//4*im*b^-1*Polylogarithms.polylog(2,-1*(im+c+-1d)^-1*(im+c+d)*exp(2a+2*b*x))+1//4*im*b^-1*Polylogarithms.polylog(2,-1*(im+d+-1c)^-1*(im+-1c+-1d)*exp(2a+2*b*x)))
@test integrate(x^2*arctan(c+(im+c)*tanh(a+b*x)), x) == :(1//3*x^3*arctan(c+(im+c)*tanh(a+b*x))+-1//12*im*b*x^4+1//6*im*x^3*log(1+im*c*exp(2a+2*b*x))+1//8*im*b^-3*Polylogarithms.polylog(4,-1*im*c*exp(2a+2*b*x))+-1//4*im*x*b^-2*Polylogarithms.polylog(3,-1*im*c*exp(2a+2*b*x))+1//4*im*b^-1*x^2*Polylogarithms.polylog(2,-1*im*c*exp(2a+2*b*x)))
@test integrate(x*arctan(c+(im+c)*tanh(a+b*x)), x) == :((1/2)*x^2*arctan(c+(im+c)*tanh(a+b*x))+-1//6*im*b*x^3+-1//8*im*b^-2*Polylogarithms.polylog(3,-1*im*c*exp(2a+2*b*x))+1//4*im*x^2*log(1+im*c*exp(2a+2*b*x))+1//4*im*x*b^-1*Polylogarithms.polylog(2,-1*im*c*exp(2a+2*b*x)))
@test integrate(arctan(c+(im+c)*tanh(a+b*x)), x) == :(x*arctan(c+(im+c)*tanh(a+b*x))+(1/2)*im*x*log(1+im*c*exp(2a+2*b*x))+-1//2*im*b*x^2+1//4*im*b^-1*Polylogarithms.polylog(2,-1*im*c*exp(2a+2*b*x)))
@test integrate(x^2*arctan(c+-1*(im+-1c)*tanh(a+b*x)), x) == :(1//3*x^3*arctan(c+-1*(im+-1c)*tanh(a+b*x))+-1//6*im*x^3*log(1+-1*im*c*exp(2a+2*b*x))+-1//8*im*b^-3*Polylogarithms.polylog(4,im*c*exp(2a+2*b*x))+1//12*im*b*x^4+-1//4*im*b^-1*x^2*Polylogarithms.polylog(2,im*c*exp(2a+2*b*x))+1//4*im*x*b^-2*Polylogarithms.polylog(3,im*c*exp(2a+2*b*x)))
@test integrate(x*arctan(c+-1*(im+-1c)*tanh(a+b*x)), x) == :((1/2)*x^2*arctan(c+-1*(im+-1c)*tanh(a+b*x))+-1//4*im*x^2*log(1+-1*im*c*exp(2a+2*b*x))+1//6*im*b*x^3+1//8*im*b^-2*Polylogarithms.polylog(3,im*c*exp(2a+2*b*x))+-1//4*im*x*b^-1*Polylogarithms.polylog(2,im*c*exp(2a+2*b*x)))
@test integrate(arctan(c+-1*(im+-1c)*tanh(a+b*x)), x) == :(x*arctan(c+-1*(im+-1c)*tanh(a+b*x))+(1/2)*im*b*x^2+-1//2*im*x*log(1+-1*im*c*exp(2a+2*b*x))+-1//4*im*b^-1*Polylogarithms.polylog(2,im*c*exp(2a+2*b*x)))
@test integrate((e+f*x)^3*arctan(coth(a+b*x)), x) == :(1//4*f^-1*(e+f*x)^4*arctan(coth(a+b*x))+1//4*f^-1*(e+f*x)^4*arctan(exp(2a+2*b*x))+-3//16*im*b^-4*f^3*Polylogarithms.polylog(5,im*exp(2a+2*b*x))+-1//4*im*b^-1*(e+f*x)^3*Polylogarithms.polylog(2,-1*im*exp(2a+2*b*x))+1//4*im*b^-1*(e+f*x)^3*Polylogarithms.polylog(2,im*exp(2a+2*b*x))+3//16*im*b^-4*f^3*Polylogarithms.polylog(5,-1*im*exp(2a+2*b*x))+-3//8*im*f*b^-2*(e+f*x)^2*Polylogarithms.polylog(3,im*exp(2a+2*b*x))+-3//8*im*b^-3*f^2*(e+f*x)*Polylogarithms.polylog(4,-1*im*exp(2a+2*b*x))+3//8*im*f*b^-2*(e+f*x)^2*Polylogarithms.polylog(3,-1*im*exp(2a+2*b*x))+3//8*im*b^-3*f^2*(e+f*x)*Polylogarithms.polylog(4,im*exp(2a+2*b*x)))
@test integrate((e+f*x)^2*arctan(coth(a+b*x)), x) == :(1//3*f^-1*(e+f*x)^3*arctan(coth(a+b*x))+1//3*f^-1*(e+f*x)^3*arctan(exp(2a+2*b*x))+-1//4*im*b^-1*(e+f*x)^2*Polylogarithms.polylog(2,-1*im*exp(2a+2*b*x))+-1//8*im*b^-3*f^2*Polylogarithms.polylog(4,-1*im*exp(2a+2*b*x))+1//4*im*b^-1*(e+f*x)^2*Polylogarithms.polylog(2,im*exp(2a+2*b*x))+1//8*im*b^-3*f^2*Polylogarithms.polylog(4,im*exp(2a+2*b*x))+-1//4*im*f*b^-2*(e+f*x)*Polylogarithms.polylog(3,im*exp(2a+2*b*x))+1//4*im*f*b^-2*(e+f*x)*Polylogarithms.polylog(3,-1*im*exp(2a+2*b*x)))
@test integrate((e+f*x)*arctan(coth(a+b*x)), x) == :((1/2)*f^-1*(e+f*x)^2*arctan(coth(a+b*x))+(1/2)*f^-1*(e+f*x)^2*arctan(exp(2a+2*b*x))+-1//4*im*b^-1*(e+f*x)*Polylogarithms.polylog(2,-1*im*exp(2a+2*b*x))+-1//8*im*f*b^-2*Polylogarithms.polylog(3,im*exp(2a+2*b*x))+1//4*im*b^-1*(e+f*x)*Polylogarithms.polylog(2,im*exp(2a+2*b*x))+1//8*im*f*b^-2*Polylogarithms.polylog(3,-1*im*exp(2a+2*b*x)))
@test integrate(arctan(coth(a+b*x)), x) == :(x*arctan(coth(a+b*x))+x*arctan(exp(2a+2*b*x))+-1//4*im*b^-1*Polylogarithms.polylog(2,-1*im*exp(2a+2*b*x))+1//4*im*b^-1*Polylogarithms.polylog(2,im*exp(2a+2*b*x)))
@test integrate(x^2*arctan(c+d*coth(a+b*x)), x) == :(1//3*x^3*arctan(c+d*coth(a+b*x))+-1//6*im*x^3*log(1+-1*(im+c+-1d)^-1*(im+c+d)*exp(2a+2*b*x))+-1//8*im*b^-3*Polylogarithms.polylog(4,(im+c+-1d)^-1*(im+c+d)*exp(2a+2*b*x))+1//6*im*x^3*log(1+-1*(im+d+-1c)^-1*(im+-1c+-1d)*exp(2a+2*b*x))+1//8*im*b^-3*Polylogarithms.polylog(4,(im+d+-1c)^-1*(im+-1c+-1d)*exp(2a+2*b*x))+-1//4*im*x*b^-2*Polylogarithms.polylog(3,(im+d+-1c)^-1*(im+-1c+-1d)*exp(2a+2*b*x))+-1//4*im*b^-1*x^2*Polylogarithms.polylog(2,(im+c+-1d)^-1*(im+c+d)*exp(2a+2*b*x))+1//4*im*x*b^-2*Polylogarithms.polylog(3,(im+c+-1d)^-1*(im+c+d)*exp(2a+2*b*x))+1//4*im*b^-1*x^2*Polylogarithms.polylog(2,(im+d+-1c)^-1*(im+-1c+-1d)*exp(2a+2*b*x)))
@test integrate(x*arctan(c+d*coth(a+b*x)), x) == :((1/2)*x^2*arctan(c+d*coth(a+b*x))+-1//4*im*x^2*log(1+-1*(im+c+-1d)^-1*(im+c+d)*exp(2a+2*b*x))+-1//8*im*b^-2*Polylogarithms.polylog(3,(im+d+-1c)^-1*(im+-1c+-1d)*exp(2a+2*b*x))+1//4*im*x^2*log(1+-1*(im+d+-1c)^-1*(im+-1c+-1d)*exp(2a+2*b*x))+1//8*im*b^-2*Polylogarithms.polylog(3,(im+c+-1d)^-1*(im+c+d)*exp(2a+2*b*x))+-1//4*im*x*b^-1*Polylogarithms.polylog(2,(im+c+-1d)^-1*(im+c+d)*exp(2a+2*b*x))+1//4*im*x*b^-1*Polylogarithms.polylog(2,(im+d+-1c)^-1*(im+-1c+-1d)*exp(2a+2*b*x)))
@test integrate(arctan(c+d*coth(a+b*x)), x) == :(x*arctan(c+d*coth(a+b*x))+(1/2)*im*x*log(1+-1*(im+d+-1c)^-1*(im+-1c+-1d)*exp(2a+2*b*x))+-1//2*im*x*log(1+-1*(im+c+-1d)^-1*(im+c+d)*exp(2a+2*b*x))+-1//4*im*b^-1*Polylogarithms.polylog(2,(im+c+-1d)^-1*(im+c+d)*exp(2a+2*b*x))+1//4*im*b^-1*Polylogarithms.polylog(2,(im+d+-1c)^-1*(im+-1c+-1d)*exp(2a+2*b*x)))
@test integrate(x^2*arctan(c+(im+c)*coth(a+b*x)), x) == :(1//3*x^3*arctan(c+(im+c)*coth(a+b*x))+-1//12*im*b*x^4+1//6*im*x^3*log(1+-1*im*c*exp(2a+2*b*x))+1//8*im*b^-3*Polylogarithms.polylog(4,im*c*exp(2a+2*b*x))+-1//4*im*x*b^-2*Polylogarithms.polylog(3,im*c*exp(2a+2*b*x))+1//4*im*b^-1*x^2*Polylogarithms.polylog(2,im*c*exp(2a+2*b*x)))
@test integrate(x*arctan(c+(im+c)*coth(a+b*x)), x) == :((1/2)*x^2*arctan(c+(im+c)*coth(a+b*x))+-1//6*im*b*x^3+-1//8*im*b^-2*Polylogarithms.polylog(3,im*c*exp(2a+2*b*x))+1//4*im*x^2*log(1+-1*im*c*exp(2a+2*b*x))+1//4*im*x*b^-1*Polylogarithms.polylog(2,im*c*exp(2a+2*b*x)))
@test integrate(arctan(c+(im+c)*coth(a+b*x)), x) == :(x*arctan(c+(im+c)*coth(a+b*x))+(1/2)*im*x*log(1+-1*im*c*exp(2a+2*b*x))+-1//2*im*b*x^2+1//4*im*b^-1*Polylogarithms.polylog(2,im*c*exp(2a+2*b*x)))
@test integrate(x^2*arctan(c+-1*(im+-1c)*coth(a+b*x)), x) == :(1//3*x^3*arctan(c+-1*(im+-1c)*coth(a+b*x))+-1//6*im*x^3*log(1+im*c*exp(2a+2*b*x))+-1//8*im*b^-3*Polylogarithms.polylog(4,-1*im*c*exp(2a+2*b*x))+1//12*im*b*x^4+-1//4*im*b^-1*x^2*Polylogarithms.polylog(2,-1*im*c*exp(2a+2*b*x))+1//4*im*x*b^-2*Polylogarithms.polylog(3,-1*im*c*exp(2a+2*b*x)))
@test integrate(x*arctan(c+-1*(im+-1c)*coth(a+b*x)), x) == :((1/2)*x^2*arctan(c+-1*(im+-1c)*coth(a+b*x))+-1//4*im*x^2*log(1+im*c*exp(2a+2*b*x))+1//6*im*b*x^3+1//8*im*b^-2*Polylogarithms.polylog(3,-1*im*c*exp(2a+2*b*x))+-1//4*im*x*b^-1*Polylogarithms.polylog(2,-1*im*c*exp(2a+2*b*x)))
@test integrate(arctan(c+-1*(im+-1c)*coth(a+b*x)), x) == :(x*arctan(c+-1*(im+-1c)*coth(a+b*x))+(1/2)*im*b*x^2+-1//2*im*x*log(1+im*c*exp(2a+2*b*x))+-1//4*im*b^-1*Polylogarithms.polylog(2,-1*im*c*exp(2a+2*b*x)))
@test integrate(arctan(exp(x)), x) == :((1/2)*im*Polylogarithms.polylog(2,-1*im*exp(x))+-1//2*im*Polylogarithms.polylog(2,im*exp(x)))
@test integrate(x*arctan(exp(x)), x) == :((1/2)*im*Polylogarithms.polylog(3,im*exp(x))+-1//2*im*Polylogarithms.polylog(3,-1*im*exp(x))+(1/2)*im*x*Polylogarithms.polylog(2,-1*im*exp(x))+-1//2*im*x*Polylogarithms.polylog(2,im*exp(x)))
@test integrate(x^2*arctan(exp(x)), x) == :(im*Polylogarithms.polylog(4,-1*im*exp(x))+-1*im*Polylogarithms.polylog(4,im*exp(x))+im*x*Polylogarithms.polylog(3,im*exp(x))+(1/2)*im*x^2*Polylogarithms.polylog(2,-1*im*exp(x))+-1*im*x*Polylogarithms.polylog(3,-1*im*exp(x))+-1//2*im*x^2*Polylogarithms.polylog(2,im*exp(x)))
@test integrate(arctan(exp(a+b*x)), x) == :((1/2)*im*b^-1*Polylogarithms.polylog(2,-1*im*exp(a+b*x))+-1//2*im*b^-1*Polylogarithms.polylog(2,im*exp(a+b*x)))
@test integrate(x*arctan(exp(a+b*x)), x) == :((1/2)*im*b^-2*Polylogarithms.polylog(3,im*exp(a+b*x))+-1//2*im*b^-2*Polylogarithms.polylog(3,-1*im*exp(a+b*x))+(1/2)*im*x*b^-1*Polylogarithms.polylog(2,-1*im*exp(a+b*x))+-1//2*im*x*b^-1*Polylogarithms.polylog(2,im*exp(a+b*x)))
@test integrate(x^2*arctan(exp(a+b*x)), x) == :(im*b^-3*Polylogarithms.polylog(4,-1*im*exp(a+b*x))+-1*im*b^-3*Polylogarithms.polylog(4,im*exp(a+b*x))+im*x*b^-2*Polylogarithms.polylog(3,im*exp(a+b*x))+(1/2)*im*b^-1*x^2*Polylogarithms.polylog(2,-1*im*exp(a+b*x))+-1*im*x*b^-2*Polylogarithms.polylog(3,-1*im*exp(a+b*x))+-1//2*im*b^-1*x^2*Polylogarithms.polylog(2,im*exp(a+b*x)))
@test integrate(arctan(a+b*f^(c+d*x)), x) == :(d^-1*log(f)^-1*arctan(a+b*f^(c+d*x))*log(2*b*f^(c+d*x)*(1+-1*im*(a+b*f^(c+d*x)))^-1*(im+-1a)^-1)+(1/2)*im*d^-1*log(f)^-1*Polylogarithms.polylog(2,1+-2*(1+-1*im*(a+b*f^(c+d*x)))^-1)+-1*d^-1*log(f)^-1*arctan(a+b*f^(c+d*x))*log(2*(1+-1*im*(a+b*f^(c+d*x)))^-1)+-1//2*im*d^-1*log(f)^-1*Polylogarithms.polylog(2,1+-2*b*f^(c+d*x)*(1+-1*im*(a+b*f^(c+d*x)))^-1*(im+-1a)^-1))
@test integrate(arctan(exp(x))*exp(-1x), x) == :(x+-1//2*log(1+exp(2x))+-1*arctan(exp(x))*exp(-1x))
@test integrate((-1+x)^-3*arctan(x), x) == :((4+-4x)^-1+-1//4*log(1+-1x)+1//8*log(1+x^2)+-1//2*(1+-1x)^-2*arctan(x))
@test integrate((4+3x)^-3*arctan(1+2x), x) == :(-1*(136+102x)^-1+-5//578*log(1+2x+2*x^2)+5//289*log(4+3x)+8//867*arctan(1+2x)+-1//6*(4+3x)^-2*arctan(1+2x))
@test integrate(arctan((1+x)^(1/2)), x) == :(-1*(1+x)^(1/2)+2*arctan((1+x)^(1/2))+x*arctan((1+x)^(1/2)))
@test integrate((1+x^2)^-1*(2+arctan(x))^-1, x) == :(log(2+arctan(x)))
@test integrate((a+a*x^2)^-1*(b+-2*b*arctan(x))^-1, x) == :(-1//2*a^-1*b^-1*log(1+-2*arctan(x)))
@test integrate((1+x)^-2*(1+x^2)^-1*(x+x^3+(1+x)^2*arctan(x)), x) == :((1+x)^-1+(1/2)*arctan(x)^2+log(1+x))
@test integrate(x^3*arctan((1+x)^(1/2)+-1*x^(1/2)), x) == :(-1//8*x^(1/2)+-1//40*x^5//2+1//8*arctan(x^(1/2))+1//24*x^3//2+1//56*x^7//2+-1//8*x^4*arctan(x^(1/2))+1//16*pi*x^4)
@test integrate(x^2*arctan((1+x)^(1/2)+-1*x^(1/2)), x) == :(-1//6*arctan(x^(1/2))+-1//18*x^3//2+1//6*x^(1/2)+1//30*x^5//2+-1//6*x^3*arctan(x^(1/2))+1//12*pi*x^3)
@test integrate(x*arctan((1+x)^(1/2)+-1*x^(1/2)), x) == :(-1//4*x^(1/2)+1//4*arctan(x^(1/2))+1//12*x^3//2+-1//4*x^2*arctan(x^(1/2))+1//8*pi*x^2)
@test integrate(arctan((1+x)^(1/2)+-1*x^(1/2)), x) == :((1/2)*x^(1/2)+-1//2*arctan(x^(1/2))+-1//2*x*arctan(x^(1/2))+1//4*pi*x)
@test integrate(x^-1*arctan((1+x)^(1/2)+-1*x^(1/2)), x) == :((1/2)*im*Polylogarithms.polylog(2,im*x^(1/2))+-1//2*im*Polylogarithms.polylog(2,-1*im*x^(1/2))+1//4*pi*log(x))
@test integrate(x^-2*arctan((1+x)^(1/2)+-1*x^(1/2)), x) == :((1/2)*x^-1//2+(1/2)*arctan(x^(1/2))+(1/2)*x^-1*arctan(x^(1/2))+-1//4*pi*x^-1)
@test integrate(x^-3*arctan((1+x)^(1/2)+-1*x^(1/2)), x) == :(-1//4*x^-1//2+-1//4*arctan(x^(1/2))+1//12*x^-3//2+-1//8*pi*x^-2+1//4*x^-2*arctan(x^(1/2)))
@test integrate(x^-4*arctan((1+x)^(1/2)+-1*x^(1/2)), x) == :(-1//18*x^-3//2+1//6*x^-1//2+1//6*arctan(x^(1/2))+1//30*x^-5//2+-1//12*pi*x^-3+1//6*x^-3*arctan(x^(1/2)))
@test integrate((d+-1*d*a^-1*c^2*x^2)^-1//2*arctan(c*x*(a+-1*c^2*x^2)^-1//2)^m, x) == :(c^-1*(1+m)^-1*(a+-1*c^2*x^2)^(1/2)*(d+-1*d*a^-1*c^2*x^2)^-1//2*arctan(c*x*(a+-1*c^2*x^2)^-1//2)^(1+m))
@test integrate((d+-1*d*a^-1*c^2*x^2)^-1//2*arctan(c*x*(a+-1*c^2*x^2)^-1//2)^2, x) == :(1//3*c^-1*(a+-1*c^2*x^2)^(1/2)*(d+-1*d*a^-1*c^2*x^2)^-1//2*arctan(c*x*(a+-1*c^2*x^2)^-1//2)^3)
@test integrate((d+-1*d*a^-1*c^2*x^2)^-1//2*arctan(c*x*(a+-1*c^2*x^2)^-1//2), x) == :((1/2)*c^-1*(a+-1*c^2*x^2)^(1/2)*(d+-1*d*a^-1*c^2*x^2)^-1//2*arctan(c*x*(a+-1*c^2*x^2)^-1//2)^2)
@test integrate((d+-1*d*a^-1*c^2*x^2)^-1//2*arctan(c*x*(a+-1*c^2*x^2)^-1//2)^-1, x) == :(c^-1*(a+-1*c^2*x^2)^(1/2)*(d+-1*d*a^-1*c^2*x^2)^-1//2*log(arctan(c*x*(a+-1*c^2*x^2)^-1//2)))
@test integrate((d+-1*d*a^-1*c^2*x^2)^-1//2*arctan(c*x*(a+-1*c^2*x^2)^-1//2)^-2, x) == :(-1*c^-1*(a+-1*c^2*x^2)^(1/2)*(d+-1*d*a^-1*c^2*x^2)^-1//2*arctan(c*x*(a+-1*c^2*x^2)^-1//2)^-1)
@test integrate((d+-1*d*a^-1*c^2*x^2)^-1//2*arctan(c*x*(a+-1*c^2*x^2)^-1//2)^-3, x) == :(-1//2*c^-1*(a+-1*c^2*x^2)^(1/2)*(d+-1*d*a^-1*c^2*x^2)^-1//2*arctan(c*x*(a+-1*c^2*x^2)^-1//2)^-2)
@test integrate((a+b*x^2)^-1//2*arctan(e*x*(-1*e^2*x^2+-1*a*b^-1*e^2)^-1//2)^m, x) == :(e^-1*(1+m)^-1*(a+b*x^2)^-1//2*(-1*e^2*x^2+-1*a*b^-1*e^2)^(1/2)*arctan(e*x*(-1*e^2*x^2+-1*a*b^-1*e^2)^-1//2)^(1+m))
@test integrate((a+b*x^2)^-1//2*arctan(e*x*(-1*e^2*x^2+-1*a*b^-1*e^2)^-1//2)^2, x) == :(1//3*e^-1*(a+b*x^2)^-1//2*(-1*e^2*x^2+-1*a*b^-1*e^2)^(1/2)*arctan(e*x*(-1*e^2*x^2+-1*a*b^-1*e^2)^-1//2)^3)
@test integrate((a+b*x^2)^-1//2*arctan(e*x*(-1*e^2*x^2+-1*a*b^-1*e^2)^-1//2), x) == :((1/2)*e^-1*(a+b*x^2)^-1//2*(-1*e^2*x^2+-1*a*b^-1*e^2)^(1/2)*arctan(e*x*(-1*e^2*x^2+-1*a*b^-1*e^2)^-1//2)^2)
@test integrate((a+b*x^2)^-1//2*arctan(e*x*(-1*e^2*x^2+-1*a*b^-1*e^2)^-1//2)^-1, x) == :(e^-1*(a+b*x^2)^-1//2*(-1*e^2*x^2+-1*a*b^-1*e^2)^(1/2)*log(arctan(e*x*(-1*e^2*x^2+-1*a*b^-1*e^2)^-1//2)))
@test integrate((a+b*x^2)^-1//2*arctan(e*x*(-1*e^2*x^2+-1*a*b^-1*e^2)^-1//2)^-2, x) == :(-1*e^-1*(a+b*x^2)^-1//2*(-1*e^2*x^2+-1*a*b^-1*e^2)^(1/2)*arctan(e*x*(-1*e^2*x^2+-1*a*b^-1*e^2)^-1//2)^-1)
@test integrate((a+b*x^2)^-1//2*arctan(e*x*(-1*e^2*x^2+-1*a*b^-1*e^2)^-1//2)^-3, x) == :(-1//2*e^-1*(a+b*x^2)^-1//2*(-1*e^2*x^2+-1*a*b^-1*e^2)^(1/2)*arctan(e*x*(-1*e^2*x^2+-1*a*b^-1*e^2)^-1//2)^-2)
@test integrate((a+b*x)^-1*arctan(c*(a+b*x))*log(d*(a+b*x)), x) == :((1/2)*im*b^-1*Polylogarithms.polylog(3,im*c*(a+b*x))+-1//2*im*b^-1*Polylogarithms.polylog(3,-1*im*c*(a+b*x))+(1/2)*im*b^-1*Polylogarithms.polylog(2,-1*im*c*(a+b*x))*log(d*(a+b*x))+-1//2*im*b^-1*Polylogarithms.polylog(2,im*c*(a+b*x))*log(d*(a+b*x)))
@test integrate(arctan(sinh(a*c+b*c*x))*exp(c*(a+b*x)), x) == :(-1*b^-1*c^-1*log(1+exp(2*c*(a+b*x)))+b^-1*c^-1*arctan(sinh(c*(a+b*x)))*exp(a*c+b*c*x))
@test integrate(arctan(cosh(a*c+b*c*x))*exp(c*(a+b*x)), x) == :(b^-1*c^-1*arctan(cosh(c*(a+b*x)))*exp(a*c+b*c*x)+-1//2*b^-1*c^-1*(1+2^(1/2))*log(3+2*2^(1/2)+exp(2*c*(a+b*x)))+-1//2*b^-1*c^-1*(1+-1*2^(1/2))*log(3+-2*2^(1/2)+exp(2*c*(a+b*x))))
@test integrate(arctan(tanh(a*c+b*c*x))*exp(c*(a+b*x)), x) == :(b^-1*c^-1*arctan(tanh(c*(a+b*x)))*exp(a*c+b*c*x)+(1/2)*2^(1/2)*b^-1*c^-1*arctan(1+-1*2^(1/2)*exp(a*c+b*c*x))+-1//2*2^(1/2)*b^-1*c^-1*arctan(1+2^(1/2)*exp(a*c+b*c*x))+-1//4*2^(1/2)*b^-1*c^-1*log(1+-1*2^(1/2)*exp(a*c+b*c*x)+exp(2*c*(a+b*x)))+1//4*2^(1/2)*b^-1*c^-1*log(1+2^(1/2)*exp(a*c+b*c*x)+exp(2*c*(a+b*x))))
@test integrate(arctan(coth(a*c+b*c*x))*exp(c*(a+b*x)), x) == :(b^-1*c^-1*arctan(coth(c*(a+b*x)))*exp(a*c+b*c*x)+(1/2)*2^(1/2)*b^-1*c^-1*arctan(1+2^(1/2)*exp(a*c+b*c*x))+-1//2*2^(1/2)*b^-1*c^-1*arctan(1+-1*2^(1/2)*exp(a*c+b*c*x))+-1//4*2^(1/2)*b^-1*c^-1*log(1+2^(1/2)*exp(a*c+b*c*x)+exp(2*c*(a+b*x)))+1//4*2^(1/2)*b^-1*c^-1*log(1+-1*2^(1/2)*exp(a*c+b*c*x)+exp(2*c*(a+b*x))))
@test integrate(arctan(sech(a*c+b*c*x))*exp(c*(a+b*x)), x) == :(b^-1*c^-1*arctan(sech(c*(a+b*x)))*exp(a*c+b*c*x)+(1/2)*b^-1*c^-1*(1+2^(1/2))*log(3+2*2^(1/2)+exp(2*c*(a+b*x)))+(1/2)*b^-1*c^-1*(1+-1*2^(1/2))*log(3+-2*2^(1/2)+exp(2*c*(a+b*x))))
@test integrate(arctan(csch(a*c+b*c*x))*exp(c*(a+b*x)), x) == :(b^-1*c^-1*log(1+exp(2*c*(a+b*x)))+b^-1*c^-1*arctan(csch(c*(a+b*x)))*exp(a*c+b*c*x))
@test integrate(x^-1*(a+b*arctan(c*x^n))*(d+e*log(f*x^m)), x) == :(a*d*log(x)+(1/2)*a*e*m^-1*log(f*x^m)^2+(1/2)*im*b*d*n^-1*Polylogarithms.polylog(2,-1*im*c*x^n)+-1//2*im*b*d*n^-1*Polylogarithms.polylog(2,im*c*x^n)+(1/2)*im*b*e*m*n^-2*Polylogarithms.polylog(3,im*c*x^n)+(1/2)*im*b*e*n^-1*Polylogarithms.polylog(2,-1*im*c*x^n)*log(f*x^m)+-1//2*im*b*e*m*n^-2*Polylogarithms.polylog(3,-1*im*c*x^n)+-1//2*im*b*e*n^-1*Polylogarithms.polylog(2,im*c*x^n)*log(f*x^m))

max_stars_repo_path: test/5 Inverse trig functions/5.3 Inverse tangent/5.3.7 Inverse tangent functions.jl


text:
module Darknet
using Libdl, CEnum

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("Darknet not installed properly, run `] build Darknet', restart Julia and try again")
end
include(depsjl_path)

# Module initialization function
function __init__()
    check_deps()
end

include("ctypes.jl")
export Ctm, Ctime_t, Cclock_t

include(joinpath(@__DIR__, "..", "gen", "libdarknet_common.jl"))
include(joinpath(@__DIR__, "..", "gen", "libdarknet_api.jl"))

include("testing.jl")

include("utils.jl")

end # module Darknet


text:
# ---
# title: 1268. Search Suggestions System
# id: problem1268
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/search-suggestions-system/description/>
# hidden: true
# ---
#
# Given an array of strings `products` and a string `searchWord`. We want to
# design a system that suggests at most three product names from `products`
# after each character of `searchWord` is typed. Suggested products should have
# common prefix with the searchWord. If there are more than three products with
# a common prefix return the three lexicographically minimums products.
#
# Return _list of lists_ of the suggested `products` after each character of
# `searchWord` is typed.
#
#
#
# **Example 1:**
#
#
#
#     Input: products = ["mobile","mouse","moneypot","monitor","mousepad"], searchWord = "mouse"
#     Output: [
#     ["mobile","moneypot","monitor"],
#     ["mobile","moneypot","monitor"],
#     ["mouse","mousepad"],
#     ["mouse","mousepad"],
#     ["mouse","mousepad"]
#     ]
#     Explanation: products sorted lexicographically = ["mobile","moneypot","monitor","mouse","mousepad"]
#     After typing m and mo all products match and we show user ["mobile","moneypot","monitor"]
#     After typing mou, mous and mouse the system suggests ["mouse","mousepad"]
#
#
# **Example 2:**
#
#
#
#     Input: products = ["havana"], searchWord = "havana"
#     Output: [["havana"],["havana"],["havana"],["havana"],["havana"],["havana"]]
#
#
# **Example 3:**
#
#
#
#     Input: products = ["bags","baggage","banner","box","cloths"], searchWord = "bags"
#     Output: [["baggage","bags","banner"],["baggage","bags","banner"],["baggage","bags"],["bags"]]
#
#
# **Example 4:**
#
#
#
#     Input: products = ["havana"], searchWord = "tatiana"
#     Output: [[],[],[],[],[],[],[]]
#
#
#
#
# **Constraints:**
#
#   * `1 <= products.length <= 1000`
#   * There are no repeated elements in `products`.
#   * `1 <= Σ products[i].length <= 2 * 10^4`
#   * All characters of `products[i]` are lower-case English letters.
#   * `1 <= searchWord.length <= 1000`
#   * All characters of `searchWord` are lower-case English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
export Activity

"""
An [`Activity`](@ref)'s type. Available values are `AT_GAME`, `AT_STREAMING`,
`AT_LISTENING`, and `AT_WATCHING`.
More details [here](https://discordapp.com/developers/docs/topics/gateway#activity-object-activity-types).
"""
@enum ActivityType AT_GAME AT_STREAMING AT_LISTENING AT_WATCHING
@boilerplate ActivityType :export :lower

"""
Flags which indicate what an [`Activity`](@ref) payload contains.
More details [here](https://discordapp.com/developers/docs/topics/gateway#activity-object-activity-flags).
"""
@enum ActivityFlags begin
    AF_INSTANCE=1<<0
    AF_JOIN=1<<1
    AF_SPECTATE=1<<2
    AF_JOIN_REQUEST=1<<3
    AF_SYNC=1<<4
    AF_PLAY=1<<5
end
@boilerplate ActivityFlags :export :lower

"""
The start and stop times of an [`Activity`](@ref).
More details [here](https://discordapp.com/developers/docs/topics/gateway#activity-object-activity-timestamps).
"""
struct ActivityTimestamps
    start::Optional{DateTime}
    stop::Optional{DateTime}
end
@boilerplate ActivityTimestamps :constructors :docs :lower :merge :mock

"""
The current party of an [`Activity`](@ref)'s player.
More details [here](https://discordapp.com/developers/docs/topics/gateway#activity-object-activity-party).
"""
struct ActivityParty
    id::Optional{String}
    size::Optional{Vector{Int}}
end
@boilerplate ActivityParty :constructors :docs :lower :merge :mock

"""
Images and hover text for an [`Activity`](@ref).
More details [here](https://discordapp.com/developers/docs/topics/gateway#activity-object-activity-assets).
"""
struct ActivityAssets
    large_image::Optional{String}
    large_text::Optional{String}
    small_image::Optional{String}
    small_text::Optional{String}
end
@boilerplate ActivityAssets :constructors :docs :lower :merge :mock

"""
Secrets for Rich Presence joining and spectating of an [`Activity`](@ref).
More details [here](https://discordapp.com/developers/docs/topics/gateway#activity-object-activity-secrets).
"""
struct ActivitySecrets
    join::Optional{String}
    spectate::Optional{String}
    match::Optional{String}
end
@boilerplate ActivitySecrets :constructors :docs :lower :merge :mock

"""
A [`User`](@ref) activity.
More details [here](https://discordapp.com/developers/docs/topics/gateway#activity-object).
"""
struct Activity
    name::String
    type::ActivityType
    url::OptionalNullable{String}
    timestamps::Optional{ActivityTimestamps}
    application_id::Optional{Snowflake}
    details::OptionalNullable{String}
    state::OptionalNullable{String}
    party::Optional{ActivityParty}
    assets::Optional{ActivityAssets}
    secrets::Optional{ActivitySecrets}
    instance::Optional{Bool}
    flags::Optional{Int}
end
@boilerplate Activity :constructors :docs :lower :merge :mock


text:
"""
normalizePixels(img)

Normalize pixels using val-min/max-min to span the interval [0,1].
...
# Arguments
- `img`: a grayscale image.

# Output
- a normalized image with pixel intensity between 0 and 1
...
"""
function normalizePixels(img)

	arr = Float64.(img)
	return Gray.((arr .- minimum(arr)) ./ (maximum(arr) - minimum(arr)))

end

max_stars_repo_name: neuro-myoung/CellSegmenter

max_issues_repo_name: neuro-myoung/CellSegmenter

max_forks_repo_name: neuro-myoung/CellSegmenter


text:
@testset "FBM" begin
    rng = MersenneTwister(42)
    h = 0.3
    k = FBMKernel(h = h)
    v1 = rand(rng, 3)
    v2 = rand(rng, 3)
    @test k(v1,v2) ≈ (sqeuclidean(v1, zero(v1))^h + sqeuclidean(v2, zero(v2))^h - sqeuclidean(v1-v2, zero(v1-v2))^h)/2 atol=1e-5
    @test repr(k) == "Fractional Brownian Motion Kernel (h = $(h))"

    test_interface(k)
    @test repr(k) == "Fractional Brownian Motion Kernel (h = $(h))"
    test_ADs(FBMKernel, ADs = [:ReverseDiff, :Zygote])
    @test_broken "Tests failing for kernelmatrix(k, x) for ForwardDiff"
    test_params(k, ([h],))
end


text:
using Diban

@time df = read_parquet("c:/scratch/nihao.parquet", multithreaded = false);

@time df = read_parquet("c:/scratch/test.parquet", [:V1, :V2, :V31]);

@time df = read_parquet("c:/scratch/test.parquet", ["V1", "V2", "V31"]);

using DataFrames

@time df2=DataFrame(df);

@time df = read_parquet("C:/git/parquet-data-collection/synthetic_data.parquet")

@time df = read_parquet("C:/scratch/nihao_zstd.parquet")

using DataFrames
@time df = DataFrame(df, copycols=false)


text:
using Makie, GeoMakie
using GeoMakie.GeoJSON

# Acquire data
states = download("https://github.com/openpolis/geojson-italy/raw/master/geojson/limits_IT_provinces.geojson")
states_bytes = read(states)
geo = GeoJSONTables.read(states_bytes)

states_str = read(states, String)
using JSON

geo = GeoJSON.dict2geo(JSON.parse(states_str))

mesh(geo, strokecolor = :blue, strokewidth = 1, color = (blue, 0.5), shading = false)

polys = AbstractPlotting.convert_arguments(AbstractPlotting.Mesh, geo)

geoms = geo |> GeoInterface.features .|> GeoInterface.geometry

sc = Scene(; scale_plot = false)
meshes = GeoMakie.toMeshes.(geoms)
poly!.(meshes; shading = :false, color = (:blue, 0.5), strokecolor = :blue, strokewidth = 2)

AbstractPlotting.current_scene()