text:
using ContinuousPOMDPTreeSearchExperiments
using ParticleFilters
using ARDESPOT
using BasicPOMCP
using POMCPOW
using POMDPs
using DiscreteValueIteration
using QMDP
using MCTS
using VDPTag2
using POMDPToolbox
using DataFrames

file_contents = readstring(@__FILE__())

pomdp = VDPTagPOMDP()
dpomdp = AODiscreteVDPTagPOMDP(pomdp, 30, 0.5)

@show max_time = 1.0
@show max_depth = 10

solvers = Dict{String, Union{Solver,Policy}}(
    "to_next" => ToNextML(mdp(pomdp)),
    "manage_uncertainty" => ManageUncertainty(pomdp, 0.01),

    "pomcpow" => begin
        rng = MersenneTwister(13)
        ro = ToNextMLSolver(rng)
        solver = POMCPOWSolver(tree_queries=10_000_000,
                               criterion=MaxUCB(65.0),
                               final_criterion=MaxTries(),
                               max_depth=max_depth,
                               max_time=max_time,
                               k_action=12.0,
                               alpha_action=1/8,
                               k_observation=1.0,
                               alpha_observation=1/30,
                               estimate_value=FORollout(ro),
                               next_action=RootToNextMLFirst(rng),
                               check_repeat_obs=false,
                               check_repeat_act=false,
                               default_action=ReportWhenUsed(TagAction(false, 0.0)),
                               rng=rng
                              )
    end,

    "pft" => begin
        rng = MersenneTwister(13)
        m = 10
        node_updater = ObsAdaptiveParticleFilter(deepcopy(pomdp),
                                           LowVarianceResampler(m),
                                           0.05, rng)
        ro = ToNextML(mdp(pomdp), rng)
        solver = DPWSolver(n_iterations=typemax(Int),
                           exploration_constant=65.0,
                           depth=max_depth,
                           max_time=max_time,
                           k_action = 12.0,
                           alpha_action = 1/8,
                           k_state = 4.0,
                           alpha_state = 1/20,
                           check_repeat_state=false,
                           check_repeat_action=false,
                           estimate_value=RolloutEstimator(ro),
                           next_action=RootToNextMLFirst(rng),
                           rng=rng
                          )
        belief_mdp = GenerativeBeliefMDP(deepcopy(pomdp), node_updater)
        solve(solver, belief_mdp)
    end,

    "d_despot" => begin
        rng = MersenneTwister(13)
        ro = ToNextMLSolver(rng)
        b = IndependentBounds(DefaultPolicyLB(ro), VDPUpper())
        sol = DESPOTSolver(lambda=0.01,
                     K=500,
                     D=10,
                     max_trials=1_000_000,
                     T_max=0.1,
                     bounds=b,
                     random_source=MemorizingSource(500, 10, rng, min_reserve=8),
                     rng=rng)
        planner = solve(sol, dpomdp)
        translate_policy(planner, dpomdp, pomdp, dpomdp)
    end

)

@show N=500

alldata = DataFrame()
# for (k, solver) in solvers
test = ["pft"]
for (k, solver) in [(s, solvers[s]) for s in test]
    @show k
    if isa(solver, Solver)
        planner = solve(solver, pomdp)
    else
        planner = solver
    end
    sims = []
    for i in 1:N
        srand(planner, i+50_000)
        filter = SIRParticleFilter(deepcopy(pomdp), 100_000, rng=MersenneTwister(i+90_000))

        md = Dict(:solver=>k, :i=>i)
        sim = Sim(deepcopy(pomdp),
                  planner,
                  filter,
                  rng=MersenneTwister(i+70_000),
                  max_steps=100,
                  metadata=md
                 )

        push!(sims, sim)
    end

    data = run_parallel(sims)
    # data = run(sims)

    rs = data[:reward]
    println(@sprintf("reward: %6.3f ± %6.3f", mean(rs), std(rs)/sqrt(length(rs))))
end

#=
datestring = Dates.format(now(), "E_d_u_HH_MM")
copyname = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "icaps_2018", "data", "subhunt_table_$(datestring).jl")
write(copyname, file_contents)
filename = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "icaps_2018", "data", "subhunt_$(datestring).csv")
println("saving to $filename...")
writetable(filename, alldata)
println("done.")
=#


text:
using FeedbackParticleFilters, Distributions, Random, PDMats

println("Testing SemigroupMethod.jl:")

@testset "SemigroupMethod.jl" begin
    print("  solver")
    x_pf=[-2.422480820086937 -0.592332203303167 -2.017301296096984 -1.5151245392598531 0.02565906919199346 0.15161614796874012;]
    testens=UnweightedParticleEnsemble(x_pf)
    eq = PoissonEquation(x->[x[1],exp(x[1])], testens)
    update!(eq, testens)
    solve!(eq, SemigroupMethod(0.1,0.01,100));
    print("-")
    @test abs(eq.potential[1,1] - -2.479601997482972) < 1E-9
    print("-")
    @test abs(eq.potential[1,2] - 1.7209004619411277) < 1E-9
    print("-")
    @test abs(eq.potential[1,3] - -2.275929639206152) < 1E-9
    print("-")
    @test abs(eq.potential[1,4] - -1.7270260296018782) < 1E-9
    print("-")
    @test abs(eq.potential[1,5] - 2.354258246950214) < 1E-9
    print("-")
    @test abs(eq.potential[1,6] - 2.4073989573996615) < 1E-9
    print("-")
    @test abs(eq.gain[1,1,1] - 0.2693220783443806) < 1E-9
    print("-")
    @test abs(eq.gain[1,2,1] - 1.976582040902698) < 1E-9
    print("-")
    @test abs(eq.gain[1,3,1] - 0.5635782390559818) < 1E-9
    print("-")
    @test abs(eq.gain[1,4,1] - 1.8205324118211719) < 1E-9
    print("-")
    @test abs(eq.gain[1,5,1] - 0.37877152214068527) < 1E-9
    print("-")
    @test abs(eq.gain[1,6,1] - 0.2722780175889505) < 1E-9
    print("-")
    @test abs(eq.potential[2,1] - -1.0121206801162257) < 1E-9
    print("-")
    @test abs(eq.potential[2,2] - 0.6859801423517747) < 1E-9
    print("-")
    @test abs(eq.potential[2,3] - -0.9514708484563714) < 1E-9
    print("-")
    @test abs(eq.potential[2,4] - -0.75223951474933) < 1E-9
    print("-")
    @test abs(eq.potential[2,5] - 0.9979730732996036) < 1E-9
    print("-")
    @test abs(eq.potential[2,6] - 1.0318778276705487) < 1E-9
    print("-")
    @test abs(eq.gain[1,1,2] - 0.0899649944564862) < 1E-9
    print("-")
    @test abs(eq.gain[1,2,2] - 0.8691936023467642) < 1E-9
    print("-")
    @test abs(eq.gain[1,3,2] - 0.20001331130430455) < 1E-9
    print("-")
    @test abs(eq.gain[1,4,2] - 0.7299601286786728) < 1E-9
    print("-")
    @test abs(eq.gain[1,5,2] - 0.18712062788471728) < 1E-9
    print("-")
    @test abs(eq.gain[1,6,2] - 0.135741019165256) < 1E-9
    println("DONE")
end#SemigroupMethod.jl

ext: jl

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
### A Pluto.jl notebook ###
# v0.16.0

using Markdown
using InteractiveUtils

# ╔═╡ a7fc3207-40eb-4ca0-af05-5f6c09f80fd4
begin
    let
        using Pkg
        Pkg.activate(mktempdir())
        Pkg.Registry.update()
        Pkg.add("PlutoUI")
        Pkg.add("Tar")
        Pkg.add("MLDataPattern")
        Pkg.add("Glob")
        Pkg.add("NIfTI")
        Pkg.add("DataAugmentation")
        Pkg.add("CairoMakie")
        Pkg.add("ImageCore")
        Pkg.add("DataLoaders")
        # Pkg.add("CUDA")
        Pkg.add("FastAI")
        Pkg.add("StaticArrays")
		# Pkg.add(url="https://github.com/Dale-Black/DistanceTransforms.jl")
    end

    using PlutoUI
    using Tar
    using MLDataPattern
    using Glob
    using NIfTI
    using DataAugmentation
    using DataAugmentation: OneHot, Image
    using CairoMakie
    using ImageCore
    using DataLoaders
    # using CUDA
    using FastAI
    using StaticArrays
	# using DistanceTransforms
end

# ╔═╡ 55b595e1-b289-4e9e-a631-ccdb4fb790f0
TableOfContents()

# ╔═╡ 352addad-5193-4769-ad30-be73119e927e
md"""
## Load data
Part of the [Medical Decathlon Dataset](http://medicaldecathlon.com/)
"""

# ╔═╡ a1506bb2-7b2f-4c24-9ef9-6a6b09ac2360
data_dir = raw"/Users/daleblack/Google Drive/Datasets/Task02_Heart"

# ╔═╡ 957ca346-209d-4152-a5b1-48f234d1b390
function loadfn_label(p)
    a = NIfTI.niread(string(p)).raw
    convert_a = convert(Array{UInt8}, a)
    convert_a = convert_a .+ 1
    return convert_a
end

# ╔═╡ 7e2ee61d-4c7c-4237-97e2-4606f5c52192
function loadfn_image(p)
    a = NIfTI.niread(string(p)).raw

    convert_a = convert(Array{Float32}, a)
    convert_a = convert_a/max(convert_a...)
    return convert_a
end

# ╔═╡ ab7fb11a-e94e-48ab-8d25-e76591dd52a7
begin
    niftidata_image(dir) = mapobs(loadfn_image, Glob.glob("*.nii*", dir))
    niftidata_label(dir) =  mapobs(loadfn_label, Glob.glob("*.nii*", dir))
    data = (
        niftidata_image(joinpath(data_dir, "imagesTr")),
        niftidata_label(joinpath(data_dir, "labelsTr")),
    )
end

# ╔═╡ ebf8306c-6c80-47c6-bac9-5a773981282a
md"""
## BlockMethod
Mimic [link](https://fluxml.ai/FastAI.jl/dev/notebooks/imagesegmentation.ipynb.html)
"""

# ╔═╡ 9c02bdc5-d46e-49f0-8848-cdcc3f77bff7
testmethod = BlockMethod(
    (FastAI.Image{3}(), Mask{3}(1:2)),
    (
        ProjectiveTransforms((112, 112, 96)),
        ImagePreprocessing(; C=Gray{N0f8}, means=SVector(1), stds=SVector(1)),
        FastAI.OneHot()
    )
)

# ╔═╡ c70b66c7-2f01-4bbd-9a20-3723f7d8168f
md"""
## Check Method
"""

# ╔═╡ 492f96d5-ffd5-427e-bae6-ea09626b450b
begin
    image, mask = sample = getobs(data, 1);
    checkblock(testmethod.blocks, sample)
end

# ╔═╡ 70170b0b-6320-428b-8ef0-64a09249a54c
describemethod(testmethod)

# ╔═╡ 23742f4f-bace-4575-8ce2-906c0cef3898
md"""
## Create model
"""

# ╔═╡ 2cec3af0-dbb6-4cb7-8050-c56425605a00
begin
    # 3D layer utilities
    conv = (stride, in, out) -> Flux.Conv((3, 3, 3), in=>out, stride=stride, pad=Flux.SamePad())
    tran = (stride, in, out) -> Flux.ConvTranspose((3, 3, 3), in=>out, stride=stride, pad=Flux.SamePad())

    conv1 = (in, out) -> Flux.Chain(conv(1, in, out), Flux.BatchNorm(out), x -> leakyrelu.(x))
    conv2 = (in, out) -> Flux.Chain(conv(2, in, out), Flux.BatchNorm(out), x -> leakyrelu.(x))
    tran2 = (in, out) -> Flux.Chain(tran(2, in, out), Flux.BatchNorm(out), x -> leakyrelu.(x))
end

# ╔═╡ 1986f393-36ed-499b-b73a-eed65f36f876
begin
    function unet3D(in_chs, lbl_chs)
        # Contracting layers
        l1 = Flux.Chain(conv1(in_chs, 4))
        l2 = Flux.Chain(l1, conv1(4, 4), conv2(4, 16))
        l3 = Flux.Chain(l2, conv1(16, 16), conv2(16, 32))
        l4 = Flux.Chain(l3, conv1(32, 32), conv2(32, 64))
        l5 = Flux.Chain(l4, conv1(64, 64), conv2(64, 128))

        # Expanding layers
        l6 = Flux.Chain(l5, tran2(128, 64), conv1(64, 64))
        l7 = Flux.Chain(Flux.Parallel(+, l6, l4), tran2(64, 32), conv1(32, 32))
        l8 = Flux.Chain(Flux.Parallel(+, l7, l3), tran2(32, 16), conv1(16, 16))
        l9 = Flux.Chain(Flux.Parallel(+, l8, l2), tran2(16, 4), conv1(4, 4))
        l10 = Flux.Chain(l9, conv1(4, lbl_chs))
    end
end

# ╔═╡ a5b998de-a93a-4f87-8569-b3da7f07642d
model = unet3D(1, 2)

# ╔═╡ e2b6a66f-aaf3-4673-99ca-7af37abd677c
md"""
## Helper functions
"""

# ╔═╡ d0cb59f8-26f6-4aab-b799-55d81d3e57b2
function dice_metric(ŷ, y)
    dice = 2 * sum(ŷ .& y) / (sum(ŷ) + sum(y))
    return dice
end

# ╔═╡ 47f8d85c-7ab3-4d54-ad42-6280cacc76e0
function as_discrete(array, logit_threshold)
    array = array .>= logit_threshold
    return array
end

# ╔═╡ 431015e9-3269-492c-9d7c-a95bc3414152
md"""
## Loss functions
"""

# ╔═╡ bdaef8df-8a85-489b-8a1e-1e527d7aede5
function dice_loss(ŷ, y)
    ϵ = 1e-5
    return loss = 1 - ((2 * sum(ŷ .* y) + ϵ) / (sum(ŷ .* ŷ) + sum(y .* y) + ϵ))
end

# ╔═╡ 0eb81afb-a055-42f4-97c3-eddd8fa1258e
function hd_loss(ŷ, y, ŷ_dtm, y_dtm)
    M = (ŷ .- y) .^ 2 .* (ŷ_dtm .^ 2 .+ y_dtm .^ 2)
    return loss = mean(M)
end

# ╔═╡ 321387d9-59b8-452d-a44c-8b4b9211d79b
md"""
## Training
"""

# ╔═╡ 610c5d25-8ad8-47ee-a530-b76640261e9d
ps = Flux.params(model);

# ╔═╡ 97799e4a-fd24-4615-a763-72b7384a8c25
loss_function = Flux.Losses.dice_coeff_loss

# ╔═╡ 40be8808-40f2-4e87-b040-a554f0e2e10f
optimizer = Flux.ADAM(0.01)

# ╔═╡ 55e0f48e-90f6-4007-837b-d93ee138fa66
traindl, validdl = methoddataloaders(data, testmethod, 2)

# ╔═╡ 4997a1eb-6c50-463e-bf3e-e359634db323
learner = Learner(model, (traindl, validdl), optimizer, loss_function)

# ╔═╡ 61f55c98-c10f-41e1-9a73-ef4b31187a55
# fitonecycle!(learner, 1, 0.033)

# ╔═╡ d5e633c5-8913-4795-827a-2e6ffaf3549f

# ╔═╡ 99f2c227-e681-41bf-b165-b92ae9664de6
# begin
#   for (xs, ys) in traindl
#       @assert size(xs) == (112, 112, 96, 1, 2)
#       @assert size(ys) == (112, 112, 96, 2, 2)
#   end
# end

# ╔═╡ acc87036-0e03-4a2a-9bc7-241f3c8491e1
# begin
#   for (xs, ys) in validdl
#       @assert size(xs) == (112, 112, 96, 1, 4)
#       @assert size(ys) == (112, 112, 96, 2, 4)
#   end
# end

# ╔═╡ bf3b5453-4a86-4849-9937-1f2e6dadab65
# with_terminal() do
#   for (xs, ys) in validdl
#   @show size(ys)
#   end
# end

# ╔═╡ 3ba4ba5a-894b-4ab3-a43f-60ec8fb9823b
# begin
#   max_epochs = 2
#   val_interval = 1
#   epoch_loss_values = []
#   val_epoch_loss_values = []
#   dice_metric_values = []
# end

# ╔═╡ 5345afc0-a149-4daa-bda9-7c36d49b2efc
# begin
#   for epoch in 1:max_epochs
#       step = 0
#       @show epoch

#       # Loop through training data
#       for (xs, ys) in traindl
#           step += 1
#           @show step

#           gs = Flux.gradient(ps) do
#               ŷs = model(xs)
#               loss = loss_function(ŷs[:, :, :, 2, :], ys[:, :, :, 2, :])
#               return loss
#           end
#           Flux.update!(optimizer, ps, gs)
#       end

#       # Loop through validation data
#       if (epoch + 1) % val_interval == 0
#           val_step = 0
#           for (val_xs, val_ys) in validdl
#               val_step += 1
#               @show val_step

#               local val_ŷs = model(val_xs)
#               local val_loss = loss_function(val_ŷs[:, :, :, 2, :], val_ys[:, :, :, 2, :])

#               val_ŷs, val_ys = as_discrete(val_ŷs, 0.5), as_discrete(val_ys, 0.5)
#           end
#       end
#   end
# end

# ╔═╡ c999e14b-cc35-4b00-9e91-ce7d5644f118
# for epoch in 1:max_epochs
#     epoch_loss = 0
#     step = 0
#     @show epoch

#     # Loop through training data
#     for (xs, ys) in train_loader
#         step += 1
#         @show step

#         xs, ys = xs |> gpu, ys |> gpu

#         ŷs = model(xs)

#         outputs_soft = softmax(ŷs; dims = 2)
#         loss_seg_dice = dice_loss(outputs_soft[:, 2, :, :, :], ys== 1)

#         gt_dtm = compute_dtm(ys)
#         seg_dtm = compute_dtm(outputs_soft[:, 2, :, :, :] .> 0.5)

#         gs = Flux.gradient(ps) do
#             loss_hd = hd_loss(outputs_soft, seg_dtm, gt_dtm)
#             loss = alpha*(loss_seg_dice) + (1 - alpha) * loss_hd
#             return loss
#         end
#         Flux.update!(optimizer, ps, gs)
#         epoch_loss += loss_function(ŷs[:, :, :, 2, :], ys[:, :, :, 2, :])

#         println("$(step/(len(train_ds) // train_loader.batchsize))")
#         println("train_loss: $(loss:.4f)")

#     end

#     epoch_loss = (epoch_loss / step)
#     push!(epoch_loss_values, epoch_loss)

#     println("epoch $(epoch_num + 1) average loss: (epoch_loss:.4f)")
#     alpha -= 0.001
#     if alpha <= 0.001:
#         alpha = 0.001
#     end

#     # Loop through validation data
#     if (epoch + 1) % val_interval == 0
#         val_step = 0
#         val_epoch_loss = 0
#         metric_step = 0
#         dice = 0
#         for (val_xs, val_ys) in val_loader
#             val_step += 1
#             @show val_step

#             val_xs, val_ys = val_xs |> gpu, val_ys |> gpu
#             val_ŷs = model(val_xs)

#             outputs_soft = softmax(val_ŷs; dims = 2)
#             loss_seg_dice = dice_loss(outputs_soft[:, 2, :, :, :], ys== 1)

#             gt_dtm = compute_dtm(val_ys)
#             seg_dtm = compute_dtm(outputs_soft[:, 2, :, :, :] .> 0.5)

#             loss_hd = hd_loss(outputs_soft, seg_dtm, gt_dtm)
#             val_loss = alpha*(loss_seg_dice) + (1 - alpha) * loss_hd

#             val_epoch_loss += val_loss

#             val_ŷs, val_ys = val_ŷs |> cpu, val_ys |> cpu
#             metric_step += 1
#             metric = dice_metric(val_ŷs[:, :, :, 2, :], val_ys[:, :, :, 2, :])
#             dice += metric
#         end

#         val_epoch_loss = (val_epoch_loss / val_step)
#         push!(val_epoch_loss_values, val_epoch_loss)

#         dice = dice / metric_step
#         push!(metric_values, dice)

#     end
# end

# ╔═╡ Cell order:
# ╠═a7fc3207-40eb-4ca0-af05-5f6c09f80fd4
# ╠═55b595e1-b289-4e9e-a631-ccdb4fb790f0
# ╟─352addad-5193-4769-ad30-be73119e927e
# ╠═a1506bb2-7b2f-4c24-9ef9-6a6b09ac2360
# ╠═957ca346-209d-4152-a5b1-48f234d1b390
# ╠═7e2ee61d-4c7c-4237-97e2-4606f5c52192
# ╠═ab7fb11a-e94e-48ab-8d25-e76591dd52a7
# ╟─ebf8306c-6c80-47c6-bac9-5a773981282a
# ╠═9c02bdc5-d46e-49f0-8848-cdcc3f77bff7
# ╟─c70b66c7-2f01-4bbd-9a20-3723f7d8168f
# ╠═492f96d5-ffd5-427e-bae6-ea09626b450b
# ╠═70170b0b-6320-428b-8ef0-64a09249a54c
# ╟─23742f4f-bace-4575-8ce2-906c0cef3898
# ╠═2cec3af0-dbb6-4cb7-8050-c56425605a00
# ╠═1986f393-36ed-499b-b73a-eed65f36f876
# ╠═a5b998de-a93a-4f87-8569-b3da7f07642d
# ╟─e2b6a66f-aaf3-4673-99ca-7af37abd677c
# ╠═d0cb59f8-26f6-4aab-b799-55d81d3e57b2
# ╠═47f8d85c-7ab3-4d54-ad42-6280cacc76e0
# ╟─431015e9-3269-492c-9d7c-a95bc3414152
# ╠═bdaef8df-8a85-489b-8a1e-1e527d7aede5
# ╠═0eb81afb-a055-42f4-97c3-eddd8fa1258e
# ╟─321387d9-59b8-452d-a44c-8b4b9211d79b
# ╠═610c5d25-8ad8-47ee-a530-b76640261e9d
# ╠═97799e4a-fd24-4615-a763-72b7384a8c25
# ╠═40be8808-40f2-4e87-b040-a554f0e2e10f
# ╠═55e0f48e-90f6-4007-837b-d93ee138fa66
# ╠═4997a1eb-6c50-463e-bf3e-e359634db323
# ╠═61f55c98-c10f-41e1-9a73-ef4b31187a55
# ╠═d5e633c5-8913-4795-827a-2e6ffaf3549f
# ╠═99f2c227-e681-41bf-b165-b92ae9664de6
# ╠═acc87036-0e03-4a2a-9bc7-241f3c8491e1
# ╠═bf3b5453-4a86-4849-9937-1f2e6dadab65
# ╠═3ba4ba5a-894b-4ab3-a43f-60ec8fb9823b
# ╠═5345afc0-a149-4daa-bda9-7c36d49b2efc
# ╠═c999e14b-cc35-4b00-9e91-ce7d5644f118


text:
export PoissonSparseAMG,solve!
using AlgebraicMultigrid
using IterativeSolvers

struct PoissonSparseAMG
    spLxy::SparseMatrixCSC{Float64,Int64}
    rgspLxy::Any

    function PoissonSparseAMG(sp)
        tmg=@elapsed rs=ruge_stuben(sp)
        println("AMG init time:",tmg)
        new(sp,rs)
    end
end

function solve!(Xxy,Bxy,sp::PoissonSparseAMG)
    ml = sp.rgspLxy
    p = aspreconditioner(ml)
    b=view(Bxy,1:length(Bxy))
    x=view(Xxy,1:length(Xxy))
    a,cgamglog=IterativeSolvers.cg!(x,sp.spLxy,b,Pl = p;log=true,reltol = 1e-16,maxiter=1000)
    @show cgamglog
    cgamglog
end


text:
#Defines the mission system that contains the reference trajectory or other mission related
#properties.

#Define Abstract Mission System for customisability
abstract type AbstractMissionSystem end

########### Template function ###########################
mutable struct MissionSystem{T <: Number} <: AbstractMissionSystem
    Active::Bool
    Ts::T
    MissionFunction::Function
    Tprev::T
end
MissionSystem() = MissionSystem(false, 0.0, x -> nothing, 0.0)
MissionSystem(Ts::Number) = MissionSystem(false, Float64(Ts), x->nothing, 0.0)
MissionSystem(Ts::Number, missfunc::Function) = MissionSystem(true, Float64(Ts), missfunc, 0.0)
MissionSystem(Active::Bool, Ts::Number, missfunc::Function) = MissionSystem(Active, Float64(Ts), missfunc, 0.0)

#Set the callback affect condition of the ControlSystem
#placeholder!
function affect!(int,dts::Type{<:AbstractMissionSystem})
    for cu in user_cache(int)
        cu.u = [dts.MissionFunction(int.u)...]
    end
  dts.Tprev = int.t
end


text:
function icn_benchmark_unit(params)
    @info "Running a benchmark unit with" params

    search_space_size = params[:domains_size]^params[:domains_size]
    if params[:search] == :complete
        if search_space_size > params[:complete_search_limit]
            @warn "Unit benchmark aborted (complete) search space is too large" search_space_size params[:complete_search_limit]
            return nothing
        end
        if search_space_size > params[:loss_sampling_threshold]
            if isnothing(params[:loss_sampler])
                @warn "Unit benchmark aborted (complete) search space is too large, and loss function is deterministically evaluated" search_space_size params[:loss_sampling_threshold] params[:loss_sampler]
                return nothing
            end
        end
    end
    if params[:search] == :partial
        if search_space_size < params[:partial_search_limit]
            @warn "Unit benchmark aborted (partial) search space is too small" search_space_size params[:partial_search_limit]
            return nothing
        end
        if search_space_size ≤ params[:loss_sampling_threshold]
            if !isnothing(params[:loss_sampler])
                @warn "Unit benchmark aborted (complete) search space is too small, and loss function is stochastically evaluated" search_space_size params[:loss_sampling_threshold] params[:loss_sampler]
                return nothing
            end
        end
    end
    json_name = savename(
        (
            con=string(params[:concept][1]),
            par=params[:concept][2],
            csl=params[:complete_search_limit],
            dom=params[:domains_size],
            gen=params[:generations],
            iter=params[:icn_iterations],
            lang=string(params[:language]),
            lst=params[:loss_sampling_threshold],
            ls=string(params[:loss_sampler]),
            metric=string(params[:metric]),
            psl=params[:partial_search_limit],
            pop=params[:population],
            sampling=params[:sampling],
            search=string(params[:search]),
            mem=params[:memoize],
        ),
        "json",
    )
    save_results = joinpath(datadir("compositions"), json_name)
    if isfile(save_results)
        @warn "The result file already exist" save_results
    else

        # Generate an appropriate parameter for the concept if relevant
        param = if isnothing(params[:concept][2])
            nothing
        elseif params[:concept][2] == 1
            rand(1:params[:domains_size])
        else
            rand(1:params[:domains_size], params[:concept][2])
        end

        # assign parameters
        constraint_concept = concept(BENCHED_CONSTRAINTS[params[:concept][1]])
        metric = params[:metric]
        domain_size = params[:domains_size]
        domains = fill(domain(1:domain_size), domain_size)

        # Time the data retrieval/generation
        t = @timed search_space(
            domain_size,
            constraint_concept,
            param;
            search=params[:search],
            complete_search_limit=params[:complete_search_limit],
            solutions_limit=params[:sampling],
        )
        solutions, non_sltns, has_data = t.value

        bench = @timed explore_learn_compose(
            domains,
            constraint_concept,
            param;
            global_iter=params[:icn_iterations],
            local_iter=params[:generations],
            metric,
            pop_size=params[:population],
            configurations=(solutions, non_sltns),
            sampler=params[:loss_sampler],
            memoize=params[:memoize],
        )
        _, _, all_compos = bench.value

        results = Dict{Any,Any}()

        # Global results
        push!(results, :data => has_data ? :loaded : :explored)
        push!(results, :data_time => t.time)
        push!(results, :data_bytes => t.bytes)
        push!(results, :data_gctime => t.gctime)
        push!(results, :data_gcstats => t.gcstats)
        push!(results, :icn_time => bench.time)
        push!(results, :icn_bytes => bench.bytes)
        push!(results, :icn_gctime => bench.gctime)
        push!(results, :icn_gcstats => bench.gcstats)
        push!(results, :total_time => t.time + bench.time)
        push!(results, :nthreads => Threads.nthreads())

        for (id, (compo, occurence)) in enumerate(pairs(all_compos))
            local_results = Dict{Symbol,Any}()

            # selection rate
            push!(local_results, :selection_rate => occurence / params[:icn_iterations])

            # Code composition
            for lang in (params[:language], :maths)
                push!(
                    local_results,
                    lang => CompositionalNetworks.code(
                        compo, lang; name=string(params[:concept][1])
                    ),
                )
            end
            push!(local_results, :symbols => CompositionalNetworks.symbols(compo))

            push!(results, :params => params)
            push!(results, id => local_results)
        end
        write(save_results, json(results, 2))
        @info "Temp results" results json_name
    end
    return nothing
end

function icn_benchmark(params=ALL_PARAMETERS; clear_results=false)
    # Ensure the folders for data output exist
    clear_results && rm(datadir(); recursive=true, force=true)
    mkpath(datadir("compositions"))

    # Run all the benchmarks for all the unit configuration from params
    configs = dict_list(params)
    @warn "Number of benchmark units is $(length(configs))"
    for (u, c) in enumerate(configs)
        @info "Starting the $u/$(length(configs)) benchmark unit"
        GC.gc()
        icn_benchmark_unit(c)
    end
    return nothing
end


text:
abstract type AbstractOptimizer end

function print_header(method::AbstractOptimizer)
    println("Iter     Function value")
end
population_size(method::AbstractOptimizer) =
    error("`population_size` is not implemented for $(summary(method)).")
metrics(method::AbstractOptimizer) = method.metrics


text:
## constructors.jl
## (c) 2014 David A. van Leeuwen

## Constructors related to the types in namedarraytypes.jl

## This code is licensed under the MIT license
## See the file LICENSE.md in this distribution

@compat letter(i) = string(Char((64+i) % 256))

## helpers for constructing names dictionaries
defaultnamesdict(names::Vector) = OrderedDict(zip(names, 1:length(names)))
defaultnamesdict(dim::Integer) = defaultnamesdict([string(i) for i in 1:dim])
defaultnamesdict(dims::NTuple) = map(defaultnamesdict, dims)

## disambiguation (Argh...)
if VERSION ≥ v"0.4-dev"
    NamedArray{T,N}(a::AbstractArray{T,N}, names::Tuple{}, dimnames::NTuple{N}) = NamedArray{T,N,typeof(a),Tuple{}}(a, (), ())
    NamedArray{T,N}(a::AbstractArray{T,N}, names::Tuple{}) = NamedArray{T,N,typeof(a),Tuple{}}(a, (), ())
end

## Basic constructor: array, tuple of dicts, tuple
## This calls the inner constructor with the appropriate types
function NamedArray{T,N}(a::AbstractArray{T,N}, names::NTuple{N,OrderedDict}, dimnames::NTuple{N})
    NamedArray{T, N, typeof(a), typeof(names)}(a, names, dimnames) ## inner constructor
end

## dimnames created as default, then inner constructor called
function NamedArray{T,N}(array::AbstractArray{T,N}, names::NTuple{N,OrderedDict})
    dimnames = [Symbol(letter(i)) for i=1:ndims(array)]
    NamedArray{T, N, typeof(array), typeof(names)}(array, names, tuple(dimnames...)) ## inner constructor
end

## constructor with array, names and dimnames (dict is created from names)
function NamedArray{T,N}(array::AbstractArray{T,N}, names::NTuple{N,Vector}, dimnames::NTuple{N})
    dicts = defaultnamesdict(names)
    NamedArray(array, dicts, dimnames)
end

## constructor with array, names (dict is created from names), dimnames created as default
function NamedArray{T,N}(array::AbstractArray{T,N}, names::NTuple{N,Vector})
    dicts = defaultnamesdict(names)
    dimnames = [Symbol(letter(i)) for i=1:length(names)]
    NamedArray(array, dicts, tuple(dimnames...))
end

## vectors instead of tuples, with defaults
function NamedArray{T,N,VT}(a::AbstractArray{T,N},
                            names::Vector{VT}=[String[string(i) for i=1:d] for d in size(a)],
                            dimnames::Vector = Symbol[Symbol(letter(i)) for i=1:N])
    length(names) == length(dimnames) == N || error("Dimension mismatch")
    if VT <: OrderedDict
        dicts = tuple(names...)
    else
        dicts = defaultnamesdict(tuple(names...))
    end
    NamedArray(a, dicts, tuple(dimnames...))
end

## Type and dimensions
"""
`NamedArray(T::Type, dims::Int...)` creates an uninitialized array with default names
for the dimensions (`:A`, `:B`, ...) and indices (`"1"`, `"2"`, ...).
"""
function NamedArray(T::DataType, dims::Int...)
    ld = length(dims)
    names = [[string(j) for j=1:i] for i=dims]
    dimnames = [@compat Symbol(letter(i)) for i=1:ld]
    a = Array(T, dims...)
    NamedArray(a, tuple(names...), tuple(dimnames...))
end


text:
using Thebes, Luxor

include(dirname(pathof(Thebes)) * "/../data/moreobjects.jl")

Drawing(800, 800, "/tmp/test1.png")
origin()
background("grey95")
setopacity(0.95)

eyepoint(Point3D(800 ,600, 500))

sethue("grey90")
carpet(600)
drawcube(500)

setline(.5)
setopacity(1)
for i in -20:10:20, j in -20:10:20, z in -20:10:20
    object = make(sphere2)
    scaleby!(object, 50, 50, 50)
    moveby!(object, i * 15, j * 15, z * 15)
    d = distance(object.vertices[1], Point3D(0, 0, 0))
    if d < 300
        pin(object)
    end
end

finish()

preview()


text:
################################################################################
#
#  Small groups
#
################################################################################

export number_of_small_groups, small_group, small_groups_limit

const small_groups_limit = 63

include("small_groups_1_63")

function number_of_small_groups(i)
  return length(small_groups_1_63[i])
end

# Data in the DB:
# <gens, rels, nontrivrels, orderdis, ordersubdis, IsAbelian(G), IsCyclic(G),
# IsSolvable(G), IsNilpotent(G), #AutomorphismGroup(G);

function small_group(i, j)
  i < 1 || i > 63 && error("Group order ($i) must be between 1 and $small_groups_limit")
  j < 1 || j > number_of_small_groups(i) && error("Index ($j) must be between 1 and $(number_of_small_groups(i))")
  P = PermGroup(i)
  Gens = [P()]
  for p in small_groups_1_63[i][j][1]
    push!(Gens, P(p))
  end
  G, _, _ = generic_group(closure(Gens, *, P()), *)
  G.issolvable = small_groups_1_63[i][j][8] ? 1 : 2
  G.isnilpotent = small_groups_1_63[i][j][9] ? 1 : 2
  return G
end


text:
using H5Sparse
using Test
using HDF5, SparseArrays

@testset "H5Sparse.jl" begin
    B = sprand(10, 10, 0.5)
    filename = tempname()
    name = "A"

    # constructors
    A = H5SparseMatrixCSC(filename, name, B)
    @test size(A) == size(B)
    @test size(A, 1) == size(B, 1)
    @test size(A, 2) == size(B, 2)
    @test eltype(A) == eltype(B)
    @test A == B

    A = H5SparseMatrixCSC(filename, name)
    @test A == B

    fid = h5open(filename, "cw")
    A = H5SparseMatrixCSC(fid, name)
    @test A == B

    A = H5SparseMatrixCSC(fid, name, 1:size(B, 1), 1:size(B, 2))
    @test sparse(A) ≈ B

    A = H5SparseMatrixCSC(fid, name, :, :)
    @test sparse(A) ≈ B

    A = H5SparseMatrixCSC(fid, name, :, 3:6)
    @test sparse(A) ≈ B[:, 3:6]

    # getitem
    A = H5SparseMatrixCSC(filename, name)
    for col in 1:size(B, 2)
        @test A[:, col] == B[:, col]
    end
    for row in 1:size(B, 1)
        @test A[row, :] == B[row, :]
    end
    for row in 1:size(B, 1)
        for col in 1:size(B, 2)
            @test A[row, col] ≈ B[row, col]
        end
    end
    for col in 1:size(B, 2)
        @test A[:, 1:col] ≈ B[:, 1:col]
    end

    # conversion
    @test sparse(A) == B
    @test sparse(A[:, 1:size(A, 2)]) == B
    @test sparse(A[:, 2:size(A, 2)-1]) == B[:, 2:size(B, 2)-1]
    @test Matrix(A) == Matrix(B)

    # concatenation
    C = sprand(10, 5, 0.5)
    A = hcat(A, C)
    @test sparse(A) ≈ hcat(B, C)
    D = sprand(10, 11, 0.1)
    A = hcat(A, D)
    @test sparse(A) ≈ hcat(B, C, D)
    A = hcat(A, C, D)
    @test sparse(A) ≈ hcat(B, C, D, C, D)

    # h5isvalidcsc
    @test H5Sparse.h5isvalidcsc(fid, name)
    @test H5Sparse.h5isvalidcsc(filename, name)
    delete_object(fid[name], "m")
    @test !H5Sparse.h5isvalidcsc(fid, name)
    fid["foo"] = 10
    @test !H5Sparse.h5isvalidcsc(fid, "foo")
end


text:
import BioMASS:isinstalled

@testset "Parameter Estimation" begin
    model_ode = Model("../examples/fos_model")
    output = []
    @testset "optimization" begin
        optimize(
            model_ode, 1, max_generation=3, popsize=3,
            local_search_method="mutation", n_children=15
        )
        lines = open(joinpath(model_ode.path, "logs", "1.log"), "r") do f
            readlines(f)
        end
        @test lines[end][1:13] == "Generation3: "
        push!(output, "logs")
        push!(output, "fitparam")
    end
    @testset "optimization_continue" begin
        optimize_continue(
            model_ode, 1, max_generation=6, popsize=3,
            local_search_method="CMAES", maxiter=30
        )
        lines = open(joinpath(model_ode.path, "logs", "1.log"), "r") do f
            readlines(f)
        end
        @test lines[end][1:13] == "Generation6: "
        # test differential_evolution
        if isinstalled("scipy.optimize")
            @testset "Differential evolution" begin
                optimize_continue(
                    model_ode, 1, max_generation=9, popsize=3,
                    local_search_method="DE", maxiter=10
                )
                lines = open(joinpath(model_ode.path, "logs", "1.log"), "r") do f
                    readlines(f)
                end
                @test lines[end][1:13] == "Generation9: "
            end
            @testset "Modified Powell's method" begin
                optimize_continue(
                    model_ode, 1, max_generation=12, popsize=3,
                    local_search_method="Powell", maxiter=5
                )
                lines = open(joinpath(model_ode.path, "logs", "1.log"), "r") do f
                    readlines(f)
                end
                @test lines[end][1:14] == "Generation12: "
            end
        end
    end
    if isinstalled("matplotlib")
        @testset "visualization" begin
            @test run_simulation(model_ode, viz_type="best") === nothing
            files = readdir(joinpath(model_ode.path, "figure", "simulation", "best"))
            n_pdf = 0
            for file in files
                if occursin(".pdf", file)
                    n_pdf += 1
                end
            end
            @test n_pdf == 8  # length(observables)
            push!(output, "figure")
        end
    end
    if isinstalled("numpy")
        @testset "conversion" begin
            @test param2biomass(model_ode.path) === nothing
            @test isdir(joinpath(model_ode.path, "dat2npy", "out", "1"))
            push!(output, "dat2npy")
        end
    end
    for dir in output
        rm(joinpath(model_ode.path, "$dir"), recursive=true, force=true)
    end
end


text:
module TestGraphviz

using Copulas
using Base.Test
using JFinM_Charts

############################
## graphviz general tests ##
############################

## get data
tr = Copulas.CTreeParRef([0, 1, 2, 2])

fname = tempname()
Copulas.render(tr, Copulas.CTreeChart(), "png", fname)

## default charts
##---------------

tr = Copulas.CTreeParRef([0, 1, 2, 2])
Copulas.render(tr)

tr2 = convert(Copulas.CTreePaths, tr)
Copulas.render(tr2)

vn = Copulas.testvine(2)
Copulas.render(vn)

################
## CTreeChart ##
################

## get data
tr = Copulas.CTreeParRef([0, 1, 2, 2])

## default chart
gviz = Copulas.GViz(tr, Copulas.CTreeChart())
Copulas.render(gviz)

@test_throws Exception gviz = Copulas.GViz(tr,
                                           JFinM_Charts.chart(Copulas.CTreeChart,
                                                              emph1 = [2, 4],
                                                              emph2 = [2, 5]))

## visualize directly
gviz = Copulas.render(tr, JFinM_Charts.chart(Copulas.CTreeChart))

## two emphasized sets
gviz = Copulas.GViz(tr, JFinM_Charts.chart(Copulas.CTreeChart,
                                           emph1 = [1, 3],
                                           emph2 = [2, 4],
                                           emphFillColor2 = "red"
                                           ))
Copulas.render(gviz)

## visualize CTreePaths instance
tr2 = convert(Copulas.CTreePaths, tr)
gviz = Copulas.GViz(tr2, JFinM_Charts.chart(Copulas.CTreeChart))
Copulas.viz(gviz)

## render data directly
Copulas.render(tr)

## visualize with chart options
gviz = Copulas.render(tr,
                   JFinM_Charts.chart(Copulas.CTreeChart,
                                      emph1 = [2],
                                      emphFillColor1 = "orange",
                                      emphShape1 = "ellipse"))

## visualize gviz object
##----------------------

## first create gviz object
gviz = Copulas.GViz(tr, JFinM_Charts.chart(Copulas.CTreeChart))
# Copulas.viz(gviz)
Copulas.render(gviz)

####################
## VineTreesChart ##
####################

vn = Copulas.testvine(3)

gviz = Copulas.GViz(vn, JFinM_Charts.chart(Copulas.VineTreesChart,
                                   emph1 = [1, 3]))
Copulas.render(gviz)

## default chart
Copulas.render(vn)

## with emphasized variables
gviz = Copulas.GViz(vn, JFinM_Charts.chart(Copulas.VineTreesChart,
                                           emph1 = [1, 3],
                                           emph2 = [4],
                                           emphFillColor2 = "red",
                                           rootEmph = [2],
                                           rootEmphFillColor1 = "orange"
                                           ))

Copulas.render(gviz)

####################
## VineGraphChart ##
####################

vn = Copulas.testvine(3)
gviz = Copulas.GViz(vn, JFinM_Charts.chart(Copulas.VineGraphChart))
Copulas.render(gviz)

chrt = JFinM_Charts.chart(Copulas.VineGraphChart,
                          maxLayer=3,
                          styles=["bold", "filled", "dotted"],
                          colors=["red", "black", "black"],
                          lens=["1.0", "3.0", "5.0"],
                          weights=["1", "1", "1"])

Copulas.render(Copulas.GViz(vn, chrt))

chrt = JFinM_Charts.chart(Copulas.VineGraphChart,
                          maxLayer=2,
                          styles=["bold", "dotted"],
                          colors=["black", "black"],
                          lens=["1.0", "1.8"],
                          weights=["1", "0"])
Copulas.render(vn, chrt)

vn = Copulas.testvine(2)
Copulas.render(vn, Copulas.VineGraphChart())

## Copulas.viz(tr)

## expGraphvizString = """
## digraph {
## node [shape=circle];
## edge [arrowhead=none];
## 1 -> 2 ;
## 2 -> 3 ;
## 2 -> 4 ;
## }"""

## @test expGraphvizString == Copulas.toConsole(tr)

## tr = Copulas.CTreePaths(6, [4,2], [5,3,1], [5,7], [5,8], [5,9],
##                         [5,10], [11], [12], [13], [14])

## Copulas.viz(tr)

## Copulas.render(tr)
## ## Copulas.render(tr, "png", "/home/chris/treePic")

end


text:
using Polylogarithms
using SpecialFunctions
using Test
using DataFrames, CSV
import Base.MathConstants: π, pi, ℯ, e, γ, eulergamma, catalan, φ, golden
include("test_defs.jl")

@testset "Harmonic numbers" begin

    @testset "    throws errors" begin
        @test_throws DomainError harmonic(-1)
        @test_throws MethodError harmonic(Float32(1.0))
    end

    @testset "    types" begin
        @test typeof(harmonic(1)) == Float64
        @test typeof(harmonic(1.0)) == Float64
        @test typeof(harmonic(complex(1.0))) == Complex{Float64}
    end

    @testset "    basics" begin
        @test harmonic(1) ≅ 1.0
        @test harmonic(2) ≅ 1.5
        @test harmonic(3) ≅ 11.0/6.0
        @test harmonic(4) ≅ 25.0/12.0
        @test harmonic(5) ≅ 137.0/60.0
    end

    @testset "    identities" begin
        for n=4:10
            @test harmonic(n) ≅ 1.0/n + harmonic(n-1)
        end

        for n=20:30
            @test harmonic(n) ≅ sum( 1.0./(1:n) )
        end
    end

    @testset "    dataset 1 (integer arguments)" begin
        data1 = CSV.read("../data/harmonic_test_data_1.csv", DataFrame)
        m = size(data1,1)
        H = Symbol("real(H_n)")
        for i=1:m
            @test harmonic( data1[i,:n] ) ≅  data1[i,H]
        end
    end

    @testset "    dataset 2 (real arguments)" begin
        data2 = CSV.read("../data/harmonic_test_data_2.csv", DataFrame)
        m = size(data2,1)
        H = Symbol("real(H_n)")
        for i=1:m
            @test harmonic( data2[i,:x] ) ≈  data2[i,H] # note the lower accuracy requirement here
        end
    end

    @testset "    dataset 3 (complex arguments)" begin
        data3 = CSV.read("../data/harmonic_test_data_3.csv", DataFrame)
        m = size(data3,1)
        zr = Symbol("Re(z)")
        zi = Symbol("Im(z)")
        Hr = Symbol("Re(H(z))")
        Hi = Symbol("Im(H(z))")
        for i=1:m
            z = data3[i,zr] + im*data3[i,zi]
            Hz = data3[i,Hr] + im*data3[i,Hi]
            @test harmonic( z ) ≈  Hz # note the lower accuracy requirement here
        end
    end
end

@testset "Generalised Harmonic numbers" begin

    @testset "    throws errors" begin
        @test_throws DomainError harmonic(-1, 1.2)
        @test_throws MethodError harmonic(1.0, 1)
        @test_throws MethodError harmonic(1.0, 1.0)
        @test_throws MethodError harmonic(complex(1.0), 1)
    end

    @testset "    output types" begin
        @test typeof(harmonic(1, 1)) == Float64
        @test typeof(harmonic(1, 1.0)) == Float64
    end

    @testset "     identities" begin
        for n=1:10
            @test harmonic(n,1.0) ≅ harmonic(n)
            @test harmonic(n,0.0) ≅ n

            for r = 3:2:9
                @test harmonic(n,r) ≅ Float64(n)^(-r) + polygamma(r-1,n)/gamma(r) + zeta(r)
            end
        end
    end

    @testset "    dataset 4 (integer arguments)" begin
        H = Symbol("real(H_{n,r})")
        data4 = CSV.read("../data/harmonic_test_data_4.csv", DataFrame)
        m = size(data4,1)
        for i=1:m
            @test harmonic( data4[i,:n], data4[i,:r] ) ≅ data4[i,H]
        end
    end

    @testset "    dataset 5 (integer n, real r)" begin
        H = Symbol("real(H_{n,r})")
        data5 = CSV.read("../data/harmonic_test_data_5.csv", DataFrame)
        m = size(data5,1)
        for i=1:m
            @test harmonic( data5[i,:n], data5[i,:r] )  ≅ data5[i,H]
        end
    end
end


text:
#!/usr/bin/env julia

module MirrorUpdater

### Blacklisted Packages. ###

const BLACKLIST = [

]

### Constants. ###

const METADATA   = Pkg.dir("METADATA")
const ORG_NAME   = "JuliaPackageMirrors"
const PACKAGES   = readdir(METADATA)
const TOTAL      = length(PACKAGES)

### Functions. ###

function rename_url(url)
    url = replace(url, "git://", "https://")
    url = replace(url, ".git", "")
end

function print_header(name, number, width)
    println("[$(lpad(number, width)) / $(lpad(TOTAL, width))] $name")
    println("="^Base.tty_size()[2])
end

function update(name)
    dir = joinpath(METADATA, name)
    if isdir(dir)
        file = joinpath(dir, "url")
        if isfile(file)
            url = strip(readall(file))
            isdir(name) || run(`git clone --mirror $(url) $(name)`)
            cd(name) do
                update_local(url)
                update_mirror(name, rename_url(url))
            end
        end
    end
end

function update_local(url)
    try
        info("fetching from origin: '$url'")
        run(`git fetch -p origin`)
    catch
        warn("could not fetch from '$(url)'")
    end
end

function update_mirror(name, url)
    # Check if mirror exists already.
    if !success(`git ls-remote git@github.com:$(ORG_NAME)/$(name).jl.git`)
        info("creating new repo: '$name'")
        run(`hub create $(ORG_NAME)/$(name).jl -d 'Julia package mirror.' -h $url`)
    end
    success(`git remote add $(ORG_NAME) git@github.com:$(ORG_NAME)/$(name).jl.git`)
    # Update the mirrored repo.
    info("updating mirror: '$name'")
    success(`git push $(ORG_NAME) --mirror`)
end

### Main. ###

function main()
    # Update METADATA repo first to find newly added packages.
    Pkg.update()

    info("updating Julia package mirrors...")
    cd(joinpath(dirname(@__FILE__), "..")) do
        width = length(digits(length(PACKAGES)))
        for (number, name) in enumerate(PACKAGES)
            name in BLACKLIST && (warn("skipping $name"); continue)
            print_header(name, number, width)
            update(name)
            println()
        end
    end
end

end

max_stars_repo_name: JuliaPackageMirrors/mirror-updater

max_issues_repo_name: JuliaPackageMirrors/mirror-updater

max_forks_repo_name: JuliaPackageMirrors/mirror-updater


text:
# Tests of Discrete samplers

using Distributions

function est_p(sampler, K::Int, n::Int)
	cnts = zeros(Int, K)
	for k = 1:n
		cnts[rand(sampler)] += 1
	end
	cnts * inv(n)
end

function est_dev(p::Vector{Float64}, sampler, n::Int)
	max(abs(est_p(sampler, length(p), n) - p))
end

const m = 10
const K = 20
const n = 10^6

dev_naive = 0.
dev_dtable = 0.
dev_alias = 0.

for i = 1:m
	p = rand(K)
	p = p / sum(p)

	dev_naive  += est_dev(p, Categorical(p), n)
	dev_dtable += est_dev(p, Distributions.DiscreteDistributionTable(p), n)
	dev_alias  += est_dev(p, Distributions.AliasTable(p), n)

	print('.')
end
println()

dev_naive /= m
dev_dtable /= m
dev_alias /= m

function print_result(name, avg_dev)
	@printf("%-15s : avg.deviation = %.6e\n", name, avg_dev)
end

print_result("Naive", dev_naive)
print_result("DiscreteTable", dev_dtable)
print_result("AliasTable", dev_alias)


text:
using Test

using InvariantSets
using LinearAlgebra: I
using Polyhedra
using CDDLib

@testset "Scaling with `*ᶜ` operator" begin
    P = Hyperrectangle(low=[-.6, -2.0],
                       high=[1.0, 2.0])
    P1 = 2.0*I *ᶜ P
    P2 = 2.0 *ᶜ P
    @test P1 ⊆ P2 && P2 ⊆ P1
    @test any(2*canonical_length(P) .== canonical_length(P2))
    P3 = convert(HPolytope, P)
    P4 = 2.0 *ᶜ P3
    @test any(2*canonical_length(P3) .== canonical_length(P4))
end

@testset "Concrete Addition operator +ᶜ" begin
    P = Hyperrectangle(low=[-.6, -2.0],
                       high=[1.0, 2.0])
    P1 = convert(HPolytope, P)
    vec =  [1.0,-10]
    P3 = P1 +ᶜ vec
    P4 = vec +ᶜ P1
    @test any(canonical_length(P3) .≈ vec .+ canonical_length(P1))
    @test any(canonical_length(P3) .== canonical_length(P4))
    P5 = P1 +ᶜ P1
    @test P5 ⊆ 2 *ᶜ P1 && 2 *ᶜ P1 ⊆ P5
end

@testset "Concrete Subtraction operator -ᶜ" begin
    P = Hyperrectangle(low=[-.6, -2.0],
                       high=[1.0, 2.0])
    P1 = convert(HPolytope, P)
    vec =  [1.0, -0.5]
    P3 = P1 -ᶜ vec
    @test any(canonical_length(P3) .≈ canonical_length(P1) .- vec )
    P4 = vec -ᶜ P1
    @test any(canonical_length(P4) .≈ vec .+ canonical_length(reflect(P1)))
    P5 = P1 -ᶜ P1
    @test P5 ⊆ BallInf(zeros(2), 0.0) &&  BallInf(zeros(2), 0.0) ⊆ P5
end

max_stars_repo_name: mforets/InvariantSets

max_issues_repo_name: mforets/InvariantSets


text:
module PolymakeOscar

using Polymake, Hecke

Hecke.nrows(A::Polymake.pm_MatrixAllocated) = Int(size(A)[1])
Hecke.ncols(A::Polymake.pm_MatrixAllocated) = Int(size(A)[2])

#solves Ax <= b
function solve_ineq(A::fmpz_mat, b::fmpz_mat)
  bA = Array{BigInt, 2}(hcat(-b, A))
  p = @pm Polytope.Polytope( :INEQUALITIES => bA)
  inner = p.INTERIOR_LATTICE_POINTS
  out = p.BOUNDARY_LATTICE_POINTS

  res = zero_matrix(FlintZZ, nrows(inner) + nrows(out), ncols(A))
  for i=1:nrows(out)
    @assert out[i,1] == 1
    for j=1:ncols(A)
      res[i,j] = out[i, j+1]
    end
  end
  for i=1:nrows(inner)
    @assert inner[i,1] == 1
    for j=1:ncols(A)
      res[i+nrows(out), j] = inner[i, j+1]
    end
  end
  return res
end

function solve_non_negative(A::fmpz_mat, b::fmpz_mat)
  bA = Array{BigInt, 2}(hcat(-b, A))
  zI = Array{BigInt, 2}(hcat(zero_matrix(FlintZZ, ncols(A), 1), MatrixSpace(FlintZZ, ncols(A), ncols(A))(1)))
  p = @pm Polytope.Polytope(:EQUATIONS => deepcopy(bA), :INEQUALITIES => deepcopy(zI))
  inner = p.INTERIOR_LATTICE_POINTS
  out = p.BOUNDARY_LATTICE_POINTS

  res = zero_matrix(FlintZZ, nrows(inner) + nrows(out), ncols(A))
  for i=1:nrows(out)
    @assert out[i,1] == 1
    for j=1:ncols(A)
      res[i,j] = out[i, j+1]
    end
  end
  for i=1:nrows(inner)
    @assert inner[i,1] == 1
    for j=1:ncols(A)
      res[i+nrows(out), j] = inner[i, j+1]
    end
  end
  return res
end

function solve_mixed(A::fmpz_mat, b::fmpz_mat, C::fmpz_mat)  # Ax == b && Cx >= 0
  bA = Array{BigInt, 2}(hcat(-b, A))
  z = findall(i->!iszero_row(C, i), 1:nrows(C))
  zI = Array{BigInt, 2}(hcat(zero_matrix(FlintZZ, nrows(C), 1), C))[z, :]
  p = @pm Polytope.Polytope(:EQUATIONS => bA, :INEQUALITIES => zI)
  inner = p.INTERIOR_LATTICE_POINTS
  out = p.BOUNDARY_LATTICE_POINTS

  res = zero_matrix(FlintZZ, nrows(inner) + nrows(out), ncols(A))
  for i=1:nrows(out)
    @assert out[i,1] == 1
    for j=1:ncols(A)
      res[i,j] = out[i, j+1]
    end
  end
  for i=1:nrows(inner)
    @assert inner[i,1] == 1
    for j=1:ncols(A)
      res[i+nrows(out), j] = inner[i, j+1]
    end
  end
  return res
end

function Hecke.valuation(a::FacElem{fmpz, FlintIntegerRing}, p::fmpz)
  return sum(k*valuation(b, p) for (b, k) = a.fac)
end

function norm_equation2_fac_elem(R::NfAbsOrd, k::fmpz; abs::Bool = false)
  @assert Hecke.ismaximal(R)
  lp = factor(k*R)
  s, ms = Hecke.sunit_mod_units_group_fac_elem(collect(keys(lp)))
  C = vcat([matrix(FlintZZ, 1, ngens(s), [valuation(ms(s[i]), p) for i=1:ngens(s)]) for p = keys(lp)])

  lp = factor(k)
  A = vcat([matrix(FlintZZ, 1, ngens(s), [valuation(Hecke.factored_norm(ms(s[i])), p) for i=1:ngens(s)]) for p = keys(lp.fac)])
  b = matrix(FlintZZ, length(lp.fac), 1, [valuation(k, p) for p = keys(lp.fac)])

  so = solve_mixed(A, b, C)
  sol = [ms(s(sub(so, i:i, 1:ncols(so)))) for i=1:nrows(so)]

  if !abs
    u, mu = unit_group_fac_elem(R)
    i = findfirst(x -> norm(mu(x)) == -1, gens(u))
    ns = [norm(x) for x = sol]
    if i === nothing
      return [sol[i] for i in 1:length(sol) if ns[i] == k]
    end
    U = mu(u[i])
    return [ ns[i] == k ? sol[i] : sol[i] * U for i = 1:length(sol)]
  end
  return sol
end

function norm_equation_fac_elem(R::NfAbsOrd, k::fmpz; abs::Bool = false)
  @assert Hecke.ismaximal(R)
  lp = factor(k)
  S = []
  for (p, k) = lp.fac
    P = prime_decomposition(R, p)
    s = solve_non_negative(matrix(FlintZZ, 1, length(P), [degree(x[1]) for x = P]), matrix(FlintZZ, 1, 1, [k]))
    push!(S, (P, [view(s, i:i, 1:ncols(s)) for i=1:nrows(s)]))
  end
  sol = []
  for x in Base.Iterators.ProductIterator(Tuple(t[2] for t = S))
    I = ideal(R, 1)
    for i = 1:length(S)
      I *= prod(S[i][1][j][1]^Int(x[i][j]) for j=1:length(S[i][1]))
    end
    fl, g = Hecke.isprincipal_fac_elem(I)
    if fl
      push!(sol, g)
    end
  end
  if !abs
    u, mu = unit_group_fac_elem(R)
    i = findfirst(x -> norm(mu(x)) == -1, gens(u))
    ns = [norm(x) for x = sol]
    if i === nothing
      return [sol[i] for i in 1:length(sol) if ns[i] == k]
    end
    U = mu(u[i])
    return [ ns[i] == k ? sol[i] : sol[i] * U for i = 1:length(sol)]
  end
  return sol
end

norm_equation_fac_elem(R::NfAbsOrd, k::Integer; abs::Bool = false) =
                            norm_equation_fac_elem(R, fmpz(k), abs = abs)

function norm_equation(R::NfAbsOrd, k::fmpz; abs::Bool = false)
  s = norm_equation_fac_elem(R, k, abs = abs)
  return [R(evaluate(x)) for x = s]
end

norm_equation(R::NfAbsOrd, k::Integer; abs::Bool = false) = norm_equation(R, fmpz(k), abs = abs)

function norm_equation_fac_elem(R::Hecke.NfRelOrd{nf_elem,Hecke.NfOrdFracIdl}, a::NfAbsOrdElem{AnticNumberField,nf_elem})

  @assert Hecke.ismaximal(R)
  Ka, mKa, mkK = absolute_field(nf(R))
  Ra = maximal_order(Ka)
  class_group(Ra)
  k = nf(parent(a))
  class_group(parent(a))

  lp = factor(Ra(mkK(k(a)))*Ra)
  la = factor(a*parent(a))
  S, mS = Hecke.sunit_mod_units_group_fac_elem(collect(keys(lp)))
  s, ms = Hecke.sunit_mod_units_group_fac_elem(collect(keys(la)))
  No = hom(S, s, elem_type(s)[preimage(ms, norm(mkK, mS(g))) for g = gens(S)])
  q, mq = quo(S, kernel(No)[1])
  q, mms = snf(q)
  mq = mq*inv(mms)

  C = vcat([matrix(FlintZZ, 1, ngens(q), [valuation(mS(preimage(mq, q[i])), p) for i=1:ngens(q)]) for p = keys(lp)])

  A = vcat([matrix(FlintZZ, 1, ngens(q), [valuation(norm(mkK, mS(preimage(mq, g))), p) for g in gens(q)]) for p = keys(la)])
  b = matrix(FlintZZ, length(la), 1, [valuation(a, p) for p = keys(la)])

  so = solve_mixed(A, b, C)
  u, mu = Hecke.unit_group_fac_elem(parent(a))
  U, mU = Hecke.unit_group_fac_elem(Ra)
  No = hom(U, u, elem_type(u)[preimage(mu, norm(mkK, mU(g))) for g = gens(U)])
  sol = []
  for i = 1:nrows(so)
    aa = mS(preimage(mq, q(sub(so, i:i, 1:ncols(so)))))
    b = norm(mkK, aa)
    c = b*inv(FacElem(k(a)))
    d = preimage(mu, c)
    fl, p = haspreimage(No, d)
    if fl
      push!(sol, FacElem(Dict(mKa(x) => v for (x, v) = (aa*inv(mU(p))).fac)))
    end
  end

  return sol
end

end


text:
# test array functions
N = 10000 # data length

# function to test similarity of vectors
cosinedist(a,b) = sum(a .* b, dims=1) ./ sqrt.(sum(a .^ 2, dims=1)) ./ sqrt.(sum(b .^ 2,dims=1))

## 1D detrending
@testset "Test detrending" begin
A = rand(N)
A .-= SeisNoise.mean(A) # remove mean
x = collect(1:N)
trendA = A .+ x         # add trend

# test 1D array in-place
@test !isa(try detrend!(trendA) catch ex ex end, Exception)    # Passes
@test cosinedist(A,trendA)[1]  > 0.999
@test eltype(trendA) == eltype(A) # check type stability
@test size(trendA) == size(A)     # check size remains same

# test allocating 1D
trendB = detrend(A .+ x)
@test cosinedist(A,trendB)[1]  > 0.999
@test eltype(trendB) == eltype(A) # check type stability
@test size(trendB) == size(A)     # check size remains same

## 2D detrending
A = rand(N,10) # detrend 10 columns
A .-= SeisNoise.mean(A,dims=1) # remove mean
trendA = A .+ x

# test 2D array in-place
@test !isa(try detrend!(trendA) catch ex ex end, Exception)    # Passes
@test all(cosinedist(A,trendA)  .> 0.999)

# test allocating 2D
trendB = detrend(A .+ x)
@test all(cosinedist(A,trendB)  .> 0.999)

## test RawData
# test in-place
R = RawData()
A = rand(N,10) # detrend 10 columns
A .-= SeisNoise.mean(A,dims=1) # remove mean
R.x = A .+ x
@test !isa(try detrend!(R) catch ex ex end, Exception)    # Passes
@test isa(detrend!(R),Nothing)
@test all(cosinedist(A,R.x)  .> 0.999)

# test allocating
R = RawData()
A = rand(N,10) # detrend 10 columns
A .-= SeisNoise.mean(A,dims=1) # remove mean
R.x = A .+ x
@test !isa(try detrend(R) catch ex ex end, Exception)    # Passes
@test isa(detrend(R),RawData)
Rnew = detrend(R)
@test all(cosinedist(A,Rnew.x)  .> 0.999)

## test CorrData
# test in-place
C = CorrData()
A = rand(N,10) # detrend 10 columns
A .-= SeisNoise.mean(A,dims=1) # remove mean
C.corr = A .+ x
@test !isa(try detrend!(C) catch ex ex end, Exception)    # Passes
@test isa(detrend!(C),Nothing)
@test all(cosinedist(A,C.corr)  .> 0.999)

# test allocating
C = CorrData()
A = rand(N,10) # detrend 10 columns
A .-= SeisNoise.mean(A,dims=1) # remove mean
C.corr = A .+ x
@test !isa(try detrend(C) catch ex ex end, Exception)    # Passes
@test isa(detrend(C),CorrData)
Cnew = detrend(C)
@test all(cosinedist(A,Cnew.corr)  .> 0.999)

## test failure on Integer data
A = rand(-2^30 : 2^30, N)
@test_throws MethodError detrend(A)
end

## test demean
@testset "Test demeaning" begin
# 1D demean in-place
A = Float64.(collect(-N:N))
x = rand() * N
meanA = A .+ x
@test !isa(try demean!(meanA) catch ex ex end, Exception)    # Passes
@test isapprox(A,meanA) # test all close
@test eltype(meanA) == eltype(A)
@test size(meanA) == size(A)

# test allocating 1D
meanB = demean(A .+ x)
@test isapprox(A,meanB) # test all close

# 2D demean
A = repeat(collect(-N:N),outer=(1,10))
x = rand(1,10) * N
meanA = A .+ x
@test !isa(try demean!(meanA) catch ex ex end, Exception)    # Passes
@test isapprox(A,meanA) # test all close

# test allocating 2D
meanB = demean(A .+ x)
@test isapprox(A,meanB) # test all close

# test failure on integer data
A = rand(-2^30 : 2^30, N)
@test_throws MethodError demean(A)

## test RawData
# test in-place
R = RawData()
A = repeat(collect(-N:N),outer=(1,10))
x = rand(1,10) * N
R.x = A .+ x
@test !isa(try demean!(R) catch ex ex end, Exception)    # Passes
@test isa(demean!(R),Nothing)
@test isapprox(A,R.x)

# test allocating
R = RawData()
R.x = A .+ x
@test !isa(try detrend(R) catch ex ex end, Exception)    # Passes
@test isa(detrend(R),RawData)
Rnew = demean(R)
@test isapprox(A,Rnew.x)

## test CorrData
# test in-place
C = CorrData()
C.corr = A .+ x
@test !isa(try demean!(C) catch ex ex end, Exception)    # Passes
@test isa(demean!(C),Nothing)
@test isapprox(A,C.corr)

# test allocating
C = CorrData()
C.corr = A .+ x
@test !isa(try demean(C) catch ex ex end, Exception)    # Passes
@test isa(detrend(C),CorrData)
Cnew = demean(C)
@test isapprox(A,Cnew.corr)
end

## test taper
@testset "Test tapering" begin
# test 1D in-place
A = ones(N)
taperA = deepcopy(A)
max_percentage = 0.05
fs = 100.
max_length = 20.

# check default taper lengths
Nperc = Int(max_percentage * N)
Nlength = Int(max_length * fs)

taperA = deepcopy(A)
@test !isa(try taper!(taperA,fs) catch ex ex end, Exception)
@test all(taperA[1:min(Nperc,Nlength)] .!= 1)
@test all(taperA[end-min(Nperc,Nlength):end] .!= 1)
@test all(taperA[min(Nperc,Nlength)+1:end-min(Nperc,Nlength)-1] .== 1)
@test eltype(taperA) == eltype(A) # check type stability
@test size(taperA) == size(A)     # check size remains same

# test 1D allocating
@test !isa(try taper(A,fs) catch ex ex end, Exception)
@test isa(taper(A,fs),typeof(A))
taperB = taper(A,fs)
@test all(taperA .== taperB)
@test eltype(taperB) == eltype(A) # check type stability

# change max_percentage
max_percentage = 0.01
Nperc = Int(max_percentage * N)
taperA = deepcopy(A)
taper!(taperA,fs,max_percentage=max_percentage)
@test all(taperA[1:Nperc] .!= 1)
@test all(taperA[end-Nperc:end] .!= 1)
@test all(taperA[Nperc+1:end-Nperc-1] .== 1)

# change max_percentage
max_length = 1.
Nlength = Int(max_length * fs)
taperA = deepcopy(A)
taper!(taperA,fs,max_length=max_length)
@test all(taperA[1:Nlength] .!= 1)
@test all(taperA[end-Nlength:end] .!= 1)
@test all(taperA[Nlength+1:end-Nlength-1] .== 1)

# test 2D in-place
A = ones(N,10)
taperA = deepcopy(A)
max_percentage = 0.05
fs = 100.
max_length = 20.

# check default taper lengths
Nperc = Int(max_percentage * N)
Nlength = Int(max_length * fs)

taperA = deepcopy(A)
@test !isa(try taper!(taperA,fs) catch ex ex end, Exception)
@test all(taperA[1:min(Nperc,Nlength),:] .!= 1)
@test all(taperA[end-min(Nperc,Nlength):end,:] .!= 1)
@test all(taperA[min(Nperc,Nlength)+1:end-min(Nperc,Nlength)-1,:] .== 1)

# change max_percentage
max_percentage = 0.01
Nperc = Int(max_percentage * N)
taperA = deepcopy(A)
taper!(taperA,fs,max_percentage=max_percentage)
@test all(taperA[1:Nperc,:] .!= 1)
@test all(taperA[end-Nperc:end,:] .!= 1)
@test all(taperA[Nperc+1:end-Nperc-1,:] .== 1)

# change max_percentage
max_length = 1.
Nlength = Int(max_length * fs)
taperA = deepcopy(A)
taper!(taperA,fs,max_length=max_length)
@test all(taperA[1:Nlength,:] .!= 1)
@test all(taperA[end-Nlength:end,:] .!= 1)
@test all(taperA[Nlength+1:end-Nlength-1,:] .== 1)

# test using Real for fs and max_length
taperA = deepcopy(A)
taper!(taperA,Int(fs),max_length=Int(max_length))
@test all(taperA[1:Nlength,:] .!= 1)
@test all(taperA[end-Nlength:end,:] .!= 1)
@test all(taperA[Nlength+1:end-Nlength-1,:] .== 1)

# test RawData
# test in-place
R = RawData()
A = ones(N,10) # detrend 10 columns
max_percentage = 0.05
fs = 100.
max_length = 20.
Nperc = Int(max_percentage * N)
Nlength = Int(max_length * fs)
R.x = A
R.fs = fs
@test !isa(try taper!(R) catch ex ex end, Exception)    # Passes
@test isa(taper!(R),Nothing)
@test all(R.x[1:min(Nperc,Nlength),:] .!= 1)
@test all(R.x[end-min(Nperc,Nlength):end,:] .!= 1)
@test all(R.x[min(Nperc,Nlength)+1:end-min(Nperc,Nlength)-1,:] .== 1)

# test allocating
R = RawData()
R.fs = fs
R.x = A
@test !isa(try taper(R) catch ex ex end, Exception)    # Passes
@test isa(taper(R),RawData)
Rnew = taper(R)
@test all(Rnew.x[1:min(Nperc,Nlength),:] .!= 1)
@test all(Rnew.x[end-min(Nperc,Nlength):end,:] .!= 1)
@test all(Rnew.x[min(Nperc,Nlength)+1:end-min(Nperc,Nlength)-1,:] .== 1)
@test all(Rnew.x[1:min(Nperc,Nlength),:] .!= 1)
@test all(Rnew.x[end-min(Nperc,Nlength):end,:] .!= 1)
@test all(Rnew.x[min(Nperc,Nlength)+1:end-min(Nperc,Nlength)-1,:] .== 1)

# test with integers
Rnew = taper(R,max_length=Int(max_length))

## test CorrData
# test in-place
C = CorrData()
C.fs = fs
C.corr = A
@test !isa(try taper!(C) catch ex ex end, Exception)    # Passes
@test isa(taper!(C),Nothing)
@test all(C.corr[1:min(Nperc,Nlength),:] .!= 1)
@test all(C.corr[end-min(Nperc,Nlength):end,:] .!= 1)
@test all(C.corr[min(Nperc,Nlength)+1:end-min(Nperc,Nlength)-1,:] .== 1)

# test allocating
C = CorrData()
C.fs = fs
C.corr = A
@test !isa(try taper(C) catch ex ex end, Exception)    # Passes
@test isa(taper(C),CorrData)
Cnew = taper(C)
@test all(Cnew.corr[1:min(Nperc,Nlength),:] .!= 1)
@test all(Cnew.corr[end-min(Nperc,Nlength):end,:] .!= 1)
@test all(Cnew.corr[min(Nperc,Nlength)+1:end-min(Nperc,Nlength)-1,:] .== 1)

# test with Integers
Cnew = taper(C,max_length=Int(max_length))
@test all(Cnew.corr[1:min(Nperc,Nlength),:] .!= 1)
@test all(Cnew.corr[end-min(Nperc,Nlength):end,:] .!= 1)
@test all(Cnew.corr[min(Nperc,Nlength)+1:end-min(Nperc,Nlength)-1,:] .== 1)

## test failure on Integer data
A = ones(Int,N)
@test_throws MethodError taper(A,fs)

## test hanningwindow
A = ones(N)
@test !isa(try hanningwindow(A,2*Nlength) catch ex ex end, Exception)
win = hanningwindow(A,2 * Nlength)
@test eltype(win) == eltype(A)     # test type stability
@test length(win) == 2 * Nlength   # test vector length
@test isapprox(win[1:Nlength],win[end:-1:Nlength+1]) # test taper mirroring
end


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

    #= ::Package:: =#

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

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

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

    @test_int [csc(x)^5/(a + a*csc(x)), x, 6, (3*atanh(cos(x)))/(2*a) - (4*cot(x))/a - (4*cot(x)^3)/(3*a) + (3*cot(x)*csc(x))/(2*a) + (cot(x)*csc(x)^3)/(a + a*csc(x))]
    @test_int [csc(x)^4/(a + a*csc(x)), x, 6, -((3*atanh(cos(x)))/(2*a)) + (2*cot(x))/a - (3*cot(x)*csc(x))/(2*a) + (cot(x)*csc(x)^2)/(a + a*csc(x))]
    @test_int [csc(x)^3/(a + a*csc(x)), x, 4, atanh(cos(x))/a - cot(x)/a - cot(x)/(a + a*csc(x))]
    @test_int [csc(x)^2/(a + a*csc(x)), x, 3, -(atanh(cos(x))/a) + cot(x)/(a + a*csc(x))]
    @test_int [csc(x)^1/(a + a*csc(x)), x, 1, -(cot(x)/(a + a*csc(x)))]
    @test_int [csc(c + d*x)^0/(a + a*csc(c + d*x)), x, 2, x/a + cot(c + d*x)/(d*(a + a*csc(c + d*x)))]
    @test_int [sin(x)^1/(a + a*csc(x)), x, 4, -(x/a) - (2*cos(x))/a + cos(x)/(a + a*csc(x))]
    @test_int [sin(x)^2/(a + a*csc(x)), x, 5, (3*x)/(2*a) + (2*cos(x))/a - (3*cos(x)*sin(x))/(2*a) + (cos(x)*sin(x))/(a + a*csc(x))]
    @test_int [sin(x)^3/(a + a*csc(x)), x, 6, -((3*x)/(2*a)) - (4*cos(x))/a + (4*cos(x)^3)/(3*a) + (3*cos(x)*sin(x))/(2*a) + (cos(x)*sin(x)^2)/(a + a*csc(x))]
    @test_int [sin(x)^4/(a + a*csc(x)), x, 7, (15*x)/(8*a) + (4*cos(x))/a - (4*cos(x)^3)/(3*a) - (15*cos(x)*sin(x))/(8*a) - (5*cos(x)*sin(x)^3)/(4*a) + (cos(x)*sin(x)^3)/(a + a*csc(x))]

    @test_int [csc(c + d*x)^0/(a + a*csc(c + d*x))^2, x, 3, x/a^2 + (4*cot(c + d*x))/(3*a^2*d*(1 + csc(c + d*x))) + cot(c + d*x)/(3*d*(a + a*csc(c + d*x))^2)]

    @test_int [csc(c + d*x)^0/(a + a*csc(c + d*x))^3, x, 4, x/a^3 + cot(c + d*x)/(5*d*(a + a*csc(c + d*x))^3) + (7*cot(c + d*x))/(15*a*d*(a + a*csc(c + d*x))^2) + (22*cot(c + d*x))/(15*d*(a^3 + a^3*csc(c + d*x)))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*csc(e+f*x)^n*(a+a*csc(e+f*x))^(m/2)=#

    @test_int [(a + a*csc(x))^(5/2), x, 5, -2*a^(5/2)*atan((sqrt(a)*cot(x))/sqrt(a + a*csc(x))) - (14*a^3*cot(x))/(3*sqrt(a + a*csc(x))) - (2/3)*a^2*cot(x)*sqrt(a + a*csc(x))]
    @test_int [(a + a*csc(x))^(3/2), x, 4, -2*a^(3/2)*atan((sqrt(a)*cot(x))/sqrt(a + a*csc(x))) - (2*a^2*cot(x))/sqrt(a + a*csc(x))]
    @test_int [(a + a*csc(x))^(1/2), x, 2, -2*sqrt(a)*atan((sqrt(a)*cot(x))/sqrt(a + a*csc(x)))]
    @test_int [1/(a + a*csc(x))^(1/2), x, 5, -((2*atan((sqrt(a)*cot(x))/sqrt(a + a*csc(x))))/sqrt(a)) + (sqrt(2)*atan((sqrt(a)*cot(x))/(sqrt(2)*sqrt(a + a*csc(x)))))/sqrt(a)]
    @test_int [1/(a + a*csc(x))^(3/2), x, 6, -((2*atan((sqrt(a)*cot(x))/sqrt(a + a*csc(x))))/a^(3/2)) + (5*atan((sqrt(a)*cot(x))/(sqrt(2)*sqrt(a + a*csc(x)))))/(2*sqrt(2)*a^(3/2)) + cot(x)/(2*(a + a*csc(x))^(3/2))]
    @test_int [1/(a + a*csc(x))^(5/2), x, 7, -((2*atan((sqrt(a)*cot(x))/sqrt(a + a*csc(x))))/a^(5/2)) + (43*atan((sqrt(a)*cot(x))/(sqrt(2)*sqrt(a + a*csc(x)))))/(16*sqrt(2)*a^(5/2)) + cot(x)/(4*(a + a*csc(x))^(5/2)) + (11*cot(x))/(16*a*(a + a*csc(x))^(3/2))]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*(d*csc(e+f*x))^(n/2)*(a+a*csc(e+f*x))^m=#

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

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

    @test_int [sqrt(csc(e + f*x))*sqrt(a + a*csc(e + f*x)), x, 2, -((2*sqrt(a)*asinh((sqrt(a)*cot(e + f*x))/sqrt(a + a*csc(e + f*x))))/f)]
    @test_int [sqrt(-csc(e + f*x))*sqrt(a - a*csc(e + f*x)), x, 2, -((2*sqrt(a)*asinh((sqrt(a)*cot(e + f*x))/sqrt(a - a*csc(e + f*x))))/f)]

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

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(d*csc(e+f*x))^(n/3)*(a+a*csc(e+f*x))^(m/2)=#

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

    @test_int [sqrt(a + a*csc(c + d*x))*csc(c + d*x)^(4/3), x, 4, -((6*a*cos(c + d*x)*csc(c + d*x)^(4/3))/(5*d*sqrt(a + a*csc(c + d*x)))) - (4*3^(3/4)*sqrt(2 + sqrt(3))*a^2*cot(c + d*x)*(1 - csc(c + d*x)^(1/3))*sqrt((1 + csc(c + d*x)^(1/3) + csc(c + d*x)^(2/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*Elliptic.F(asin((1 - sqrt(3) - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))), -7 - 4*sqrt(3)))/(5*d*sqrt((1 - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*(a - a*csc(c + d*x))*sqrt(a + a*csc(c + d*x)))]
    @test_int [sqrt(a + a*csc(c + d*x))*csc(c + d*x)^(1/3), x, 3, -((2*3^(3/4)*sqrt(2 + sqrt(3))*a^2*cot(c + d*x)*(1 - csc(c + d*x)^(1/3))*sqrt((1 + csc(c + d*x)^(1/3) + csc(c + d*x)^(2/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*Elliptic.F(asin((1 - sqrt(3) - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))), -7 - 4*sqrt(3)))/(d*sqrt((1 - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*(a - a*csc(c + d*x))*sqrt(a + a*csc(c + d*x))))]
    @test_int [sqrt(a + a*csc(c + d*x))/csc(c + d*x)^(2/3), x, 4, -((3*a*cos(c + d*x)*csc(c + d*x)^(1/3))/(2*d*sqrt(a + a*csc(c + d*x)))) - (3^(3/4)*sqrt(2 + sqrt(3))*a^2*cot(c + d*x)*(1 - csc(c + d*x)^(1/3))*sqrt((1 + csc(c + d*x)^(1/3) + csc(c + d*x)^(2/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*Elliptic.F(asin((1 - sqrt(3) - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))), -7 - 4*sqrt(3)))/(2*d*sqrt((1 - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*(a - a*csc(c + d*x))*sqrt(a + a*csc(c + d*x)))]

    @test_int [sqrt(a + a*csc(c + d*x))*csc(c + d*x)^(5/3), x, 6, (24*a*cot(c + d*x))/(7*d*(1 + sqrt(3) - csc(c + d*x)^(1/3))*sqrt(a + a*csc(c + d*x))) - (6*a*cos(c + d*x)*csc(c + d*x)^(5/3))/(7*d*sqrt(a + a*csc(c + d*x))) - (12*3^(1/4)*sqrt(2 - sqrt(3))*a^2*cot(c + d*x)*(1 - csc(c + d*x)^(1/3))*sqrt((1 + csc(c + d*x)^(1/3) + csc(c + d*x)^(2/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*Elliptic.E(asin((1 - sqrt(3) - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))), -7 - 4*sqrt(3)))/(7*d*sqrt((1 - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*(a - a*csc(c + d*x))*sqrt(a + a*csc(c + d*x))) + (8*sqrt(2)*3^(3/4)*a^2*cot(c + d*x)*(1 - csc(c + d*x)^(1/3))*sqrt((1 + csc(c + d*x)^(1/3) + csc(c + d*x)^(2/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*Elliptic.F(asin((1 - sqrt(3) - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))), -7 - 4*sqrt(3)))/(7*d*sqrt((1 - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*(a - a*csc(c + d*x))*sqrt(a + a*csc(c + d*x)))]
    @test_int [sqrt(a + a*csc(c + d*x))*csc(c + d*x)^(2/3), x, 5, (6*a*cot(c + d*x))/(d*(1 + sqrt(3) - csc(c + d*x)^(1/3))*sqrt(a + a*csc(c + d*x))) - (3*3^(1/4)*sqrt(2 - sqrt(3))*a^2*cot(c + d*x)*(1 - csc(c + d*x)^(1/3))*sqrt((1 + csc(c + d*x)^(1/3) + csc(c + d*x)^(2/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*Elliptic.E(asin((1 - sqrt(3) - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))), -7 - 4*sqrt(3)))/(d*sqrt((1 - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*(a - a*csc(c + d*x))*sqrt(a + a*csc(c + d*x))) + (2*sqrt(2)*3^(3/4)*a^2*cot(c + d*x)*(1 - csc(c + d*x)^(1/3))*sqrt((1 + csc(c + d*x)^(1/3) + csc(c + d*x)^(2/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*Elliptic.F(asin((1 - sqrt(3) - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))), -7 - 4*sqrt(3)))/(d*sqrt((1 - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*(a - a*csc(c + d*x))*sqrt(a + a*csc(c + d*x)))]
    @test_int [sqrt(a + a*csc(c + d*x))/csc(c + d*x)^(1/3), x, 6, -((3*a*cot(c + d*x))/(d*(1 + sqrt(3) - csc(c + d*x)^(1/3))*sqrt(a + a*csc(c + d*x)))) - (3*a*cos(c + d*x)*csc(c + d*x)^(2/3))/(d*sqrt(a + a*csc(c + d*x))) + (3*3^(1/4)*sqrt(2 - sqrt(3))*a^2*cot(c + d*x)*(1 - csc(c + d*x)^(1/3))*sqrt((1 + csc(c + d*x)^(1/3) + csc(c + d*x)^(2/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*Elliptic.E(asin((1 - sqrt(3) - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))), -7 - 4*sqrt(3)))/(2*d*sqrt((1 - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*(a - a*csc(c + d*x))*sqrt(a + a*csc(c + d*x))) - (sqrt(2)*3^(3/4)*a^2*cot(c + d*x)*(1 - csc(c + d*x)^(1/3))*sqrt((1 + csc(c + d*x)^(1/3) + csc(c + d*x)^(2/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*Elliptic.F(asin((1 - sqrt(3) - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))), -7 - 4*sqrt(3)))/(d*sqrt((1 - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*(a - a*csc(c + d*x))*sqrt(a + a*csc(c + d*x)))]
    @test_int [sqrt(a + a*csc(c + d*x))/csc(c + d*x)^(4/3), x, 7, -((15*a*cot(c + d*x))/(8*d*(1 + sqrt(3) - csc(c + d*x)^(1/3))*sqrt(a + a*csc(c + d*x)))) - (3*a*cos(c + d*x))/(4*d*csc(c + d*x)^(1/3)*sqrt(a + a*csc(c + d*x))) - (15*a*cos(c + d*x)*csc(c + d*x)^(2/3))/(8*d*sqrt(a + a*csc(c + d*x))) + (15*3^(1/4)*sqrt(2 - sqrt(3))*a^2*cot(c + d*x)*(1 - csc(c + d*x)^(1/3))*sqrt((1 + csc(c + d*x)^(1/3) + csc(c + d*x)^(2/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*Elliptic.E(asin((1 - sqrt(3) - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))), -7 - 4*sqrt(3)))/(16*d*sqrt((1 - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*(a - a*csc(c + d*x))*sqrt(a + a*csc(c + d*x))) - (5*3^(3/4)*a^2*cot(c + d*x)*(1 - csc(c + d*x)^(1/3))*sqrt((1 + csc(c + d*x)^(1/3) + csc(c + d*x)^(2/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*Elliptic.F(asin((1 - sqrt(3) - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))), -7 - 4*sqrt(3)))/(4*sqrt(2)*d*sqrt((1 - csc(c + d*x)^(1/3))/(1 + sqrt(3) - csc(c + d*x)^(1/3))^2)*(a - a*csc(c + d*x))*sqrt(a + a*csc(c + d*x)))]

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

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

    @test_int [csc(c + d*x)^n*sqrt(a + a*csc(c + d*x)), x, 2, -((2*a*cot(c + d*x)*HypergeometricFunctions._₂F₁(1/2, 1 - n, 3/2, 1 - csc(c + d*x)))/(d*sqrt(a + a*csc(c + d*x))))]
    @test_int [csc(c + d*x)^n*sqrt(a - a*csc(c + d*x)), x, 3, -((2*a*cos(c + d*x)*csc(c + d*x)^(1 + n)*HypergeometricFunctions._₂F₁(1/2, 1 - n, 3/2, 1 + csc(c + d*x)))/((-csc(c + d*x))^n*(d*sqrt(a - a*csc(c + d*x)))))]

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

    @test_int [csc(e + f*x)^3*(a + a*csc(e + f*x))^m, x, 5, If(13>=8, (cot(e + f*x)*(a + a*csc(e + f*x))^m)/(f*(2 + 3*m + m^2)) - (cot(e + f*x)*(a + a*csc(e + f*x))^(1 + m))/(a*f*(2 + m)) - (2^(1/2 + m)*(1 + m + m^2)*cot(e + f*x)*(1 + csc(e + f*x))^(-(1/2) - m)*(a + a*csc(e + f*x))^m*HypergeometricFunctions._₂F₁(1/2, 1/2 - m, 3/2, (1/2)*(1 - csc(e + f*x))))/(f*(1 + m)*(2 + m)), (cot(e + f*x)*(a + a*csc(e + f*x))^m)/(f*(2 + 3*m + m^2)) - (cot(e + f*x)*(a + a*csc(e + f*x))^(1 + m))/(a*f*(2 + m)) - (2^(1/2 + m)*(1 + m + m^2)*cot(e + f*x)*(1 + csc(e + f*x))^(-(1/2) - m)*(a + a*csc(e + f*x))^m*HypergeometricFunctions._₂F₁(1/2, 1/2 - m, 3/2, (1/2)*(1 - csc(e + f*x))))/(f*(2 + 3*m + m^2)))]
    @test_int [csc(e + f*x)^2*(a + a*csc(e + f*x))^m, x, 4, -((cot(e + f*x)*(a + a*csc(e + f*x))^m)/(f*(1 + m))) - (2^(1/2 + m)*m*cot(e + f*x)*(1 + csc(e + f*x))^(-(1/2) - m)*(a + a*csc(e + f*x))^m*HypergeometricFunctions._₂F₁(1/2, 1/2 - m, 3/2, (1/2)*(1 - csc(e + f*x))))/(f*(1 + m))]
    @test_int [csc(e + f*x)^1*(a + a*csc(e + f*x))^m, x, 3, -((2^(1/2 + m)*cot(e + f*x)*(1 + csc(e + f*x))^(-(1/2) - m)*(a + a*csc(e + f*x))^m*HypergeometricFunctions._₂F₁(1/2, 1/2 - m, 3/2, (1/2)*(1 - csc(e + f*x))))/f)]
    @test_int [csc(e + f*x)^0*(a + a*csc(e + f*x))^m, x, 3, -((sqrt(2)*AppellF1(1/2 + m, 1/2, 1, 3/2 + m, (1/2)*(1 + csc(e + f*x)), 1 + csc(e + f*x))*cot(e + f*x)*(a + a*csc(e + f*x))^m)/(f*(1 + 2*m)*sqrt(1 - csc(e + f*x))))]
    @test_int [sin(e + f*x)^1*(a + a*csc(e + f*x))^m, x, 3, (sqrt(2)*AppellF1(1/2 + m, 1/2, 2, 3/2 + m, (1/2)*(1 + csc(e + f*x)), 1 + csc(e + f*x))*cot(e + f*x)*(a + a*csc(e + f*x))^m)/(f*(1 + 2*m)*sqrt(1 - csc(e + f*x)))]
    @test_int [sin(e + f*x)^2*(a + a*csc(e + f*x))^m, x, 3, -((sqrt(2)*AppellF1(1/2 + m, 1/2, 3, 3/2 + m, (1/2)*(1 + csc(e + f*x)), 1 + csc(e + f*x))*cot(e + f*x)*(a + a*csc(e + f*x))^m)/(f*(1 + 2*m)*sqrt(1 - csc(e + f*x))))]

    #= ::Section:: =#
    #=Integrands*of*the*form*(a+a*csc(e+f*x))^m*(d*sin(e+f*x))^n=#

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

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

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

    @test_int [(a + b*csc(c + d*x))^4, x, 6, a^4*x - (2*a*b*(2*a^2 + b^2)*atanh(cos(c + d*x)))/d - (b^2*(17*a^2 + 2*b^2)*cot(c + d*x))/(3*d) - (4*a*b^3*cot(c + d*x)*csc(c + d*x))/(3*d) - (b^2*cot(c + d*x)*(a + b*csc(c + d*x))^2)/(3*d)]
    @test_int [(a + b*csc(c + d*x))^3, x, 5, a^3*x - (b*(6*a^2 + b^2)*atanh(cos(c + d*x)))/(2*d) - (5*a*b^2*cot(c + d*x))/(2*d) - (b^2*cot(c + d*x)*(a + b*csc(c + d*x)))/(2*d)]
    @test_int [(a + b*csc(c + d*x))^2, x, 4, a^2*x - (2*a*b*atanh(cos(c + d*x)))/d - (b^2*cot(c + d*x))/d]

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

    @test_int [csc(x)^5/(a + b*csc(x)), x, 9, (a*(2*a^2 + b^2)*atanh(cos(x)))/(2*b^4) - (2*a^4*atanh((a + b*tan(x/2))/sqrt(a^2 - b^2)))/(b^4*sqrt(a^2 - b^2)) - ((3*a^2 + 2*b^2)*cot(x))/(3*b^3) + (a*cot(x)*csc(x))/(2*b^2) - (cot(x)*csc(x)^2)/(3*b)]
    @test_int [csc(x)^4/(a + b*csc(x)), x, 8, -(((2*a^2 + b^2)*atanh(cos(x)))/(2*b^3)) + (2*a^3*atanh((a + b*tan(x/2))/sqrt(a^2 - b^2)))/(b^3*sqrt(a^2 - b^2)) + (a*cot(x))/b^2 - (cot(x)*csc(x))/(2*b)]
    @test_int [csc(x)^3/(a + b*csc(x)), x, 7, (a*atanh(cos(x)))/b^2 - (2*a^2*atanh((a + b*tan(x/2))/sqrt(a^2 - b^2)))/(b^2*sqrt(a^2 - b^2)) - cot(x)/b]
    @test_int [csc(x)^2/(a + b*csc(x)), x, 6, -(atanh(cos(x))/b) + (2*a*atanh((a + b*tan(x/2))/sqrt(a^2 - b^2)))/(b*sqrt(a^2 - b^2))]
    @test_int [csc(x)^1/(a + b*csc(x)), x, 4, -((2*atanh((a + b*tan(x/2))/sqrt(a^2 - b^2)))/sqrt(a^2 - b^2))]
    @test_int [csc(c + d*x)^0/(a + b*csc(c + d*x)), x, 4, x/a + (2*b*atanh((a + b*tan((1/2)*(c + d*x)))/sqrt(a^2 - b^2)))/(a*sqrt(a^2 - b^2)*d)]
    @test_int [sin(x)^1/(a + b*csc(x)), x, 6, -((b*x)/a^2) - (2*b^2*atanh((a + b*tan(x/2))/sqrt(a^2 - b^2)))/(a^2*sqrt(a^2 - b^2)) - cos(x)/a]
    @test_int [sin(x)^2/(a + b*csc(x)), x, 7, ((a^2 + 2*b^2)*x)/(2*a^3) + (2*b^3*atanh((a + b*tan(x/2))/sqrt(a^2 - b^2)))/(a^3*sqrt(a^2 - b^2)) + (b*cos(x))/a^2 - (cos(x)*sin(x))/(2*a)]
    @test_int [sin(x)^3/(a + b*csc(x)), x, 8, -((b*(a^2 + 2*b^2)*x)/(2*a^4)) - (2*b^4*atanh((a + b*tan(x/2))/sqrt(a^2 - b^2)))/(a^4*sqrt(a^2 - b^2)) - ((2*a^2 + 3*b^2)*cos(x))/(3*a^3) + (b*cos(x)*sin(x))/(2*a^2) - (cos(x)*sin(x)^2)/(3*a)]
    @test_int [sin(x)^4/(a + b*csc(x)), x, 9, ((3*a^4 + 4*a^2*b^2 + 8*b^4)*x)/(8*a^5) + (2*b^5*atanh((a + b*tan(x/2))/sqrt(a^2 - b^2)))/(a^5*sqrt(a^2 - b^2)) + (b*(2*a^2 + 3*b^2)*cos(x))/(3*a^4) - ((3*a^2 + 4*b^2)*cos(x)*sin(x))/(8*a^3) + (b*cos(x)*sin(x)^2)/(3*a^2) - (cos(x)*sin(x)^3)/(4*a)]

    @test_int [1/(a + b*csc(c + d*x))^2, x, 6, x/a^2 + (2*b*(2*a^2 - b^2)*atanh((a + b*tan((1/2)*(c + d*x)))/sqrt(a^2 - b^2)))/(a^2*(a^2 - b^2)^(3/2)*d) - (b^2*cot(c + d*x))/(a*(a^2 - b^2)*d*(a + b*csc(c + d*x)))]
    @test_int [1/(a + b*csc(c + d*x))^3, x, 7, x/a^3 + (b*(6*a^4 - 5*a^2*b^2 + 2*b^4)*atanh((a + b*tan((1/2)*(c + d*x)))/sqrt(a^2 - b^2)))/(a^3*(a^2 - b^2)^(5/2)*d) - (b^2*cot(c + d*x))/(2*a*(a^2 - b^2)*d*(a + b*csc(c + d*x))^2) - (b^2*(5*a^2 - 2*b^2)*cot(c + d*x))/(2*a^2*(a^2 - b^2)^2*d*(a + b*csc(c + d*x)))]
    @test_int [1/(a + b*csc(c + d*x))^4, x, 8, x/a^4 + (b*(8*a^6 - 8*a^4*b^2 + 7*a^2*b^4 - 2*b^6)*atanh((a + b*tan((1/2)*(c + d*x)))/sqrt(a^2 - b^2)))/(a^4*(a^2 - b^2)^(7/2)*d) - (b^2*cot(c + d*x))/(3*a*(a^2 - b^2)*d*(a + b*csc(c + d*x))^3) - (b^2*(8*a^2 - 3*b^2)*cot(c + d*x))/(6*a^2*(a^2 - b^2)^2*d*(a + b*csc(c + d*x))^2) - (b^2*(26*a^4 - 17*a^2*b^2 + 6*b^4)*cot(c + d*x))/(6*a^3*(a^2 - b^2)^3*d*(a + b*csc(c + d*x)))]

    @test_int [1/(3 + 5*csc(c + d*x)), x, 2, -(x/12) - (5*atan(cos(c + d*x)/(3 + sin(c + d*x))))/(6*d)]
    @test_int [1/(5 + 3*csc(c + d*x)), x, 5, x/5 + (3*log(3*cos((1/2)*(c + d*x)) + sin((1/2)*(c + d*x))))/(20*d) - (3*log(cos((1/2)*(c + d*x)) + 3*sin((1/2)*(c + d*x))))/(20*d)]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*csc(e+f*x)^m*(a+b*csc(e+f*x))^(n/2)=#

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

    @test_int [csc(e + f*x)^3*(a + b*csc(e + f*x))^m, x, 8, -((cot(e + f*x)*(a + b*csc(e + f*x))^(1 + m))/(b*f*(2 + m))) + (sqrt(2)*a*(a + b)*AppellF1(1/2, 1/2, -1 - m, 3/2, (1/2)*(1 - csc(e + f*x)), (b*(1 - csc(e + f*x)))/(a + b))*cot(e + f*x)*(a + b*csc(e + f*x))^m)/(((a + b*csc(e + f*x))/(a + b))^m*(b^2*f*(2 + m)*sqrt(1 + csc(e + f*x)))) - (sqrt(2)*(a^2 + b^2*(1 + m))*AppellF1(1/2, 1/2, -m, 3/2, (1/2)*(1 - csc(e + f*x)), (b*(1 - csc(e + f*x)))/(a + b))*cot(e + f*x)*(a + b*csc(e + f*x))^m)/(((a + b*csc(e + f*x))/(a + b))^m*(b^2*f*(2 + m)*sqrt(1 + csc(e + f*x))))]
    @test_int [csc(e + f*x)^2*(a + b*csc(e + f*x))^m, x, 7, -((sqrt(2)*(a + b)*AppellF1(1/2, 1/2, -1 - m, 3/2, (1/2)*(1 - csc(e + f*x)), (b*(1 - csc(e + f*x)))/(a + b))*cot(e + f*x)*(a + b*csc(e + f*x))^m)/(((a + b*csc(e + f*x))/(a + b))^m*(b*f*sqrt(1 + csc(e + f*x))))) + (sqrt(2)*a*AppellF1(1/2, 1/2, -m, 3/2, (1/2)*(1 - csc(e + f*x)), (b*(1 - csc(e + f*x)))/(a + b))*cot(e + f*x)*(a + b*csc(e + f*x))^m)/(((a + b*csc(e + f*x))/(a + b))^m*(b*f*sqrt(1 + csc(e + f*x))))]
    @test_int [csc(e + f*x)^1*(a + b*csc(e + f*x))^m, x, 3, -((sqrt(2)*AppellF1(1/2, 1/2, -m, 3/2, (1/2)*(1 - csc(e + f*x)), (b*(1 - csc(e + f*x)))/(a + b))*cot(e + f*x)*(a + b*csc(e + f*x))^m)/(((a + b*csc(e + f*x))/(a + b))^m*(f*sqrt(1 + csc(e + f*x)))))]
    @test_int [csc(e + f*x)^0*(a + b*csc(e + f*x))^m, x, 0, Unintegrable((a + b*csc(e + f*x))^m, x)]
    @test_int [sin(e + f*x)^1*(a + b*csc(e + f*x))^m, x, 0, Unintegrable((a + b*csc(e + f*x))^m*sin(e + f*x), x)]
    @test_int [sin(e + f*x)^2*(a + b*csc(e + f*x))^m, x, 0, Unintegrable((a + b*csc(e + f*x))^m*sin(e + f*x)^2, x)]

    #= ::Section:: =#
    #=Integrands*of*the*form*(a+b*csc(e+f*x))^m*(d*sin(e+f*x))^n=#
    end

max_stars_repo_path: test/4 Trig functions/4.6 Cosecant/4.6.1.2 (d csc)^n (a+b csc)^m.jl


text:
"""
    t_nowiggle(k, ωm, fbaryon)

Transfer function for the linear matter power spectrum without Baryon Acoustic Oscillation (the so-called *no-wiggle* transfer function).

*Reference*: Equation (29-31) of Eisenstein and Hu, ApJ, 496, 605 (1998)

# Arguments
- `k::Real`: the comoving wavenumber **in units of 1/Mpc**. I.e., there is **no h** in the wavenumber.
- `ωm::Real`: the physical baryon density parameter, `ωm` = Ωm h^2
- `fbaryon::Real`: the baryon fraction,  `fbaryon` = Ωb/Ωb

This function is based on [Cosmology Routine Library (CRL)](https://wwwmpa.mpa-garching.mpg.de/~komatsu/crl/).
"""
function t_nowiggle(k::Real, ωm::Real, fb::Real) # k is in units of Mpc^{-1} [without h]
    α = 1 - 0.328fb * log(431ωm) + 0.38 * fb^2 * log(22.3ωm)
    s = 44.5 * log(9.83 / ωm) / √(1 + 10 * (fb * ωm)^(3 / 4))
    Γ = ωm * (α + (1 - α) / (1 + (0.43 * k * s)^4))
    q = k * (2.725 / 2.7)^2 / Γ
    log(2 * exp(1) + 1.8q) /
    (log(2 * exp(1) + 1.8q) + (14.2 + 731 / (1 + 62.5q)) * q^2)
end


text:
@testset "SE3PDController" begin
    val = Valkyrie()
    mechanism = val.mechanism
    state = MechanismState(mechanism)
    rng = MersenneTwister(1)
    base = rand(rng, bodies(mechanism))
    body = rand(rng, bodies(mechanism))
    angular = let angle = π / 2, axis = SVector(1.0, 0.0, 0.0), y0 = one(Quat), yf = Quat(AngleAxis(angle, axis...))
        Interpolated(0.0, 1.0, y0, yf)
    end
    linear = Interpolated(0.0, 1.0, SVector(0.0, 1.0, 2.0), SVector(2.0, 3.0, 4.0))
    frame = default_frame(body)
    gains = SE3PDGains(FramePDGains(frame, PDGains(100.0, 20.0)), FramePDGains(frame, PDGains(1000.0, 200.0)))
    trajectory = SE3Trajectory(default_frame(body), default_frame(base), angular, linear)
    weight = Diagonal(vcat(zeros(3), fill(10.0, 3)))
    controller = SE3PDController(BodyID(base), BodyID(body), trajectory, weight, gains)
    controller(0.5, state)

    # FIXME: proper tests
end

@testset "parameterized contacts" begin
    # Construct a mechanism consisting of a single body which can
    # rotate about its origin
    world = RigidBody{Float64}("world")
    mechanism = Mechanism(world, gravity=SVector(0, 0, -9.81))
    frame = CartesianFrame3D("body")
    inertia = SpatialInertia(frame, SDiagonal(1., 1, 1), SVector(0., 0, 0), 10.0)
    body = RigidBody(inertia)
    joint = Joint("rx", Revolute(SVector(1., 0, 0)))
    attach!(mechanism, world, body, joint)

    contactmodel = SoftContactModel(hunt_crossley_hertz(k = 500e3), ViscoelasticCoulombModel(0.8, 20e3, 100.))
    add_contact_point!(body, Contact.ContactPoint(Point3D(default_frame(body), 0., 0, 0), contactmodel))

    # Make the contact normal a parameter, so that it updates its representation
    # in body frame to match a fixed orientation in world frame
    state = MechanismState(mechanism)
    model = Parametron.Model(defaultoptimizer())
    position = Point3D(default_frame(body), 0., 0, 0)
    μ = 1.0
    normal = let state = state, body = body
        Parameter(model) do
            transform(state, FreeVector3D(root_frame(state.mechanism), 0., 0, 1), default_frame(body))
        end
    end
    controller_contact = QPControl.ContactPoint{4}(position, normal, μ, state, model)
    controller_contact.maxnormalforce[] = 1e3
    @test QPControl.isenabled(controller_contact)

    # Constrain the contact force to be non-zero for testing
    @constraint(model, controller_contact.force_local.v == [0.0, 0.0, 1.0])

    for θ in range(-π, stop=π, length=10)
        set_configuration!(state, [θ])
        solve!(model)
        # Sanity check our constraint
        @test value.(Ref(model), controller_contact.force_local.v) ≈ [0.0, 0.0, 1.0]

        # No matter how we rotate the robot, the contact-aligned frame will still
        # be aligned to the contact normal, which is fixed in world frame. So the
        # linear component of the contact wrench in world frame will always be
        # along [0, 0, 1].
        @test value.(Ref(model), linear(controller_contact.wrench_world)) ≈ [0.0, 0.0, 1.0]
    end
end

@testset "fixed base joint space control, constrained = $constrained" for constrained in [true, false]
    Random.seed!(42)
    mechanism = rand_tree_mechanism(Float64, Prismatic{Float64}, Revolute{Float64}, Revolute{Float64})
    N = 4
    controller = MomentumBasedController{N}(mechanism, defaultoptimizer())
    tasks = Dict{Joint{Float64}, JointAccelerationTask}()
    for joint in tree_joints(mechanism)
        task = JointAccelerationTask(joint)
        tasks[joint] = task
        if constrained
            addtask!(controller, task)
        else
            weight = 1.0
            addtask!(controller, task, weight)
        end
        setdesired!(task, rand(num_velocities(joint)))
    end
    state = MechanismState(mechanism)
    rand!(state)
    τ = similar(velocity(state))
    controller(τ, 0.0, state)
    allocs = @allocated controller(τ, 0.0, state)
    @test allocs == 0

    result = DynamicsResult(mechanism)
    dynamics!(result, state, τ)
    for joint in tree_joints(mechanism)
        @test result.v̇[joint] ≈ tasks[joint].desired atol = 1e-10
    end
end

"""
Automatically load contact points from each body in the mechanism and add them
to the controller. If `parametric_contact_surface=true`, then the `normal` and
`μ` associated with each contact will be Parameters set to random, but fixed,
values.
"""
function set_up_valkyrie_contacts!(controller::MomentumBasedController; parametric_contact_surface=false)
    valmechanism = controller.state.mechanism
    for body in bodies(valmechanism)
        for point in RBD.contact_points(body)
            position = RBD.Contact.location(point)
            if parametric_contact_surface
                normal = let frame = position.frame, direction = normalize(randn(SVector{3}))
                    Parameter(controller.qpmodel) do
                        FreeVector3D(frame, direction)
                    end
                end
                μ = let μ = rand()
                    Parameter(() -> μ, controller.qpmodel)
                end
            else
                normal = FreeVector3D(position.frame, 0.0, 0.0, 1.0)
                μ = point.model.friction.μ
            end
            addcontact!(controller, body, position, normal, μ)
        end
    end
end

@testset "zero velocity free fall" begin
    Random.seed!(5354)
    val = Valkyrie()
    mechanism = val.mechanism
    floatingjoint = val.basejoint
    N = 4
    controller = MomentumBasedController{N}(mechanism, defaultoptimizer(), floatingjoint=floatingjoint)
    set_up_valkyrie_contacts!(controller)
    state = MechanismState(mechanism)
    τ = similar(velocity(state))

    zero!(state)
    rand_configuration!(state)
    for joint in tree_joints(mechanism)
        joint == floatingjoint && continue
        regularize!(controller, joint, 1.0)
    end
    controller(τ, 0., state)
    result = DynamicsResult(mechanism)
    accels = result.accelerations
    RBD.spatial_accelerations!(accels, state, controller.result.v̇)
    for joint in tree_joints(mechanism)
        if joint == floatingjoint
            baseaccel = relative_acceleration(accels, val.pelvis, root_body(mechanism))
            baseaccel = transform(state, baseaccel, frame_after(floatingjoint))
            angularaccel = FreeVector3D(baseaccel.frame, baseaccel.angular)
            @test isapprox(angularaccel, FreeVector3D(frame_after(floatingjoint), zeros(SVector{3})), atol = 1e-4)
            linearaccel = FreeVector3D(baseaccel.frame, baseaccel.linear)
            linearaccel = transform_to_root(state, linearaccel.frame) * linearaccel
            @test isapprox(linearaccel, mechanism.gravitational_acceleration; atol = 1e-4)
        else
            v̇joint = controller.result.v̇[velocity_range(state, joint)]
            @test isapprox(v̇joint, zeros(num_velocities(joint)); atol = 1e-4)
        end
    end
    allocs = @allocated controller(τ, 0., state)
    @test allocs == 0
end

const MAX_NORMAL_FORCE_FIXME = 1e9

@testset "achievable momentum rate" begin
    Random.seed!(533454)
    val = Valkyrie()
    mechanism = val.mechanism
    floatingjoint = val.basejoint
    state = MechanismState(mechanism)
    τ = similar(velocity(state))

    N = 4
    controller = MomentumBasedController{N}(mechanism, defaultoptimizer(), floatingjoint=floatingjoint)

    set_up_valkyrie_contacts!(controller; parametric_contact_surface=true)
    ḣtask = MomentumRateTask(mechanism, centroidal_frame(controller))
    addtask!(controller, ḣtask) #, 1.0)

    for joint in tree_joints(mechanism)
        regularize!(controller, joint, 1e-6)
    end

    for p in range(0., stop=1., length=5)
        rand!(state)
        com = center_of_mass(state)
        centroidal_to_world = Transform3D(centroidal_frame(controller), com.frame, com.v)
        world_to_centroidal = inv(centroidal_to_world)

        # set random active contacts and random achievable wrench
        fg = world_to_centroidal * (mass(mechanism) * mechanism.gravitational_acceleration)
        ḣdes = Wrench(zero(fg), fg)
        for body in keys(controller.contacts)
            for contact in controller.contacts[body]
                active = rand() < p
                if active
                    normal = Parametron.evalarg(contact.normal)
                    μ = Parametron.evalarg(contact.μ)
                    contact.weight[] = 1e-6
                    contact.maxnormalforce[] = MAX_NORMAL_FORCE_FIXME
                    fnormal = 50. * rand()
                    μreduced = sqrt(2) / 2 * μ # due to polyhedral inner approximation; assumes 4 basis vectors or more
                    ftangential = μreduced * fnormal * rand() * cross(normal, FreeVector3D(normal.frame, normalize(randn(SVector{3}))))
                    f = fnormal * normal + ftangential
                    @assert isapprox(ftangential ⋅ normal, 0., atol = 1e-12)
                    @assert norm(f - (normal ⋅ f) * normal) ≤ μreduced * (normal ⋅ f)
                    wrench = Wrench(contact.position, f)
                    ḣdes += transform(wrench, world_to_centroidal * transform_to_root(state, wrench.frame))
                else
                    disable!(contact)
                end
            end
        end
        setdesired!(ḣtask, ḣdes)

        controller(τ, 0., state)

        # Ensure that desired momentum rate is achieved.
        ḣ = Wrench(momentum_matrix(state), controller.result.v̇) + momentum_rate_bias(state)
        ḣ = transform(ḣ, world_to_centroidal)
        @test isapprox(ḣdes, ḣ; atol = 1e-3)

        allocs = @allocated controller(τ, 0., state)
        @test allocs == 0
    end
end

@enum TaskMode begin
    CONSTRAINT
    OBJECTIVE_TERM_SCALAR_WEIGHT
    OBJECTIVE_TERM_SCALAR_WEIGHT_PARAMETER
    OBJECTIVE_TERM_MATRIX_WEIGHT
    OBJECTIVE_TERM_MATRIX_WEIGHT_PARAMETER
end

@testset "spatial acceleration, mode = $mode" for mode in instances(TaskMode)
    Random.seed!(533)
    val = Valkyrie()
    mechanism = val.mechanism
    floatingjoint = val.basejoint
    state = MechanismState(mechanism)
    τ = similar(velocity(state))

    rand!(state)
    N = 4
    controller = MomentumBasedController{N}(mechanism, defaultoptimizer(), floatingjoint=floatingjoint)
    body = val.feet[left]
    base = val.palms[right]
    frame = default_frame(base)
    task = SpatialAccelerationTask(mechanism, path(mechanism, base, body), frame=frame)

    result = DynamicsResult(mechanism)
    accels = result.accelerations
    if mode == CONSTRAINT
        addtask!(controller, task)
        regularize!.(Ref(controller), tree_joints(mechanism), 1.0)
    elseif mode == OBJECTIVE_TERM_SCALAR_WEIGHT
        addtask!(controller, task, 1.0)
    elseif mode == OBJECTIVE_TERM_SCALAR_WEIGHT_PARAMETER
        weight = Parameter(controller.qpmodel, val=1.0)
        addtask!(controller, task, weight)
    elseif mode == OBJECTIVE_TERM_MATRIX_WEIGHT
        addtask!(controller, task, Matrix(I, 6, 6))
    elseif mode == OBJECTIVE_TERM_MATRIX_WEIGHT_PARAMETER
        weight = Parameter(controller.qpmodel, val=Matrix(I, 6, 6))
        addtask!(controller, task, weight)
    else
        error()
    end
    desiredaccel = rand(SpatialAcceleration{Float64}, default_frame(body), default_frame(base), frame)
    setdesired!(task, desiredaccel)
    controller(τ, 0., state)
    v̇ = controller.result.v̇
    spatial_accelerations!(accels, state, v̇)
    accel = relative_acceleration(accels, body, base)
    accel = transform(state, accel, frame)
    @test isapprox(accel, desiredaccel, atol = 1e-8)

    allocs = @allocated controller(τ, 0., state)
    @test allocs == 0
end


text:
module BigO

using BenchmarkTools
using RecipesBase
using Statistics
using Zygote: gradient

export RunReport
export bigos, guessbigo, best

include("run_report.jl")
include("order_estimation.jl")

end


text:
#===============================================================================

    FermiFCI.jl - LR, June 2021

===============================================================================#
module FermiFCI

# Those are the exported functions.
# (the ones that can be called without the namespace prefix)
export Orbital, FullState, SpinState, OneBodyCoeffTensor, TwoBodyCoeffTensor, construct_hamiltonian, diagonalize, mb_state_energy

# Holds type definitions, should be loaded first.
include("typedefs.jl")
include("state_reps/state_functions.jl")

# Some utilities for basis preparation.
include("utils/plain_hilbert_construction.jl")

# Density matrices.
include("utils/density_matrix.jl")
include("utils/density_profile.jl")

# Stuff for constructing the actual Hamiltonian.
include("construction.jl")

# Stuff for diagonalization.
include("diagonalization.jl")

# Define a submodule IO.
module Utils
    include("io/io_helper.jl")
end # module IO

end # module FermiFCI