text:
# ROCm-specific runtime libraries

## GPU runtime library

# reset the runtime cache from global scope, so that any change triggers recompilation
GPUCompiler.reset_runtime()

# load or build the runtime for the most likely compilation job given a compute capability
function load_runtime(dev_isa::String)
    target = GCNCompilerTarget(; dev_isa=dev_isa)
    dummy_source = FunctionSpec(()->return, Tuple{})
    params = ROCCompilerParams()
    job = CompilerJob(target, dummy_source, params)
    GPUCompiler.load_runtime(job)
end

#@inline exception_flag() = ccall("extern julia_exception_flag", llvmcall, Ptr{Cvoid}, ())

function signal_exception()
#=
    ptr = exception_flag()
    if ptr !== C_NULL
        unsafe_store!(convert(Ptr{Int}, ptr), 1)
        threadfence_system()
    else
        @rocprintf("""
            WARNING: could not signal exception status to the host, execution will continue.
                     Please file a bug.
            """)
    end
=#
    return
end

function report_exception(ex)
#=
    @rocprintf("""
        ERROR: a %s was thrown during kernel execution.
               Run Julia on debug level 2 for device stack traces.
        """, ex)
=#
    return
end

report_oom(sz) = nothing #@rocprintf("ERROR: Out of dynamic GPU memory (trying to allocate %i bytes)\n", sz)

function report_exception_name(ex)
#=
    @rocprintf("""
        ERROR: a %s was thrown during kernel execution.
        Stacktrace:
        """, ex)
=#
    return
end

function report_exception_frame(idx, func, file, line)
    #@rocprintf(" [%i] %s at %s:%i\n", idx, func, file, line)
    return
end

## ROCm device library

const libcache = Dict{String, LLVM.Module}()

function load_device_libs(dev_isa)
    device_libs_path === nothing && return

    isa_short = replace(dev_isa, "gfx"=>"")
    device_libs = LLVM.Module[]
    bitcode_files = (
        "hc.amdgcn.bc",
        "hip.amdgcn.bc",
        "irif.amdgcn.bc",
        "ockl.amdgcn.bc",
        "oclc_isa_version_$isa_short.amdgcn.bc",
        "opencl.amdgcn.bc",
        "ocml.amdgcn.bc",
    )

    for file in bitcode_files
        ispath(joinpath(device_libs_path, file)) || continue
        name, ext = splitext(file)
        lib = get!(libcache, name) do
            file_path = joinpath(device_libs_path, file)
            open(file_path) do io
                parse(LLVM.Module, read(file_path), JuliaContext())
            end
        end
        push!(device_libs, lib)
    end

    @assert !isempty(device_libs) "No device libs detected!"
    return device_libs
end

function link_device_libs!(mod::LLVM.Module, dev_isa::String, undefined_fns)
    libs::Vector{LLVM.Module} = load_device_libs(dev_isa)

    # TODO: only link if used
    # TODO: make these globally/locally configurable
    link_oclc_defaults!(mod, dev_isa)
    for lib in libs
        # override libdevice's triple and datalayout to avoid warnings
        triple!(lib, triple(mod))
        datalayout!(lib, datalayout(mod))
    end

    GPUCompiler.link_library!(mod, libs)
end

function link_oclc_defaults!(mod::LLVM.Module, dev_isa::String; finite_only=false,
                             unsafe_math=false, correctly_rounded_sqrt=true, daz=false)
    # link in some defaults for OCLC knobs, to prevent undefined variable errors
    lib = LLVM.Module("OCLC")
    triple!(lib, triple(mod))
    datalayout!(lib, datalayout(mod))

    isa_short = replace(dev_isa, "gfx"=>"")
    for (name,value) in (
            "__oclc_ISA_version"=>parse(Int32, isa_short),
            "__oclc_finite_only_opt"=>Int32(finite_only),
            "__oclc_unsafe_math_opt"=>Int32(unsafe_math),
            "__oclc_correctly_rounded_sqrt32"=>Int32(correctly_rounded_sqrt),
            "__oclc_daz_opt"=>Int32(daz))
        gvtype = convert(LLVMType, typeof(value))
        gv = GlobalVariable(lib, gvtype, name, 4)
        init = ConstantInt(Int32(0), JuliaContext())
        initializer!(gv, init)
        unnamed_addr!(gv, true)
        constant!(gv, true)
    end

    link!(mod, lib)
end


text:
using Plots, Printf
pyplot()
viz = true

@views function sia_2D()
	# physics
	lx   = 10.0
	ly   = 10.0
	n    = 3
	nτ   = 10000
	# numerics
	nx   = 127
	ny   = 127
	nout = 100
	γ    = 0.96
	ε    = 1e-8
	∂x   = lx/nx
	∂y   = ly/ny
	xc   = LinRange(∂x/2, lx-∂x/2, nx)
	yc   = LinRange(∂y/2, ly-∂y/2, ny)
	❄    = zeros(nx  ,ny  )
	qx   = zeros(nx+1,ny  )
	qy   = zeros(nx  ,ny+1)
	∂τ   = zeros(nx  ,ny  )
	Res❄ = zeros(nx  ,ny  )
	Err❄ = zeros(nx  ,ny  )
	b    = 0.5*ones(nx  ,ny  )
	❄    = exp.(.-(xc.-lx./2.0).^2 .-(yc.-ly./2.0)'.^2)
	rad  = (xc.-lx./2.0).^2 .+(yc.-ly./2.0)'.^2; b[rad.>lx/4] .= -0.5
	# action
	t0   = Base.time()
	for iτ = 1:nτ
		Err❄ .= ❄
		qx[2:end-1,:] .= .-0.5.*(❄[1:end-1,:].+❄[2:end,:]).^n .*diff(❄, dims=1)./∂x
		qy[:,2:end-1] .= .-0.5.*(❄[:,1:end-1].+❄[:,2:end]).^n .*diff(❄, dims=2)./∂y
		∂τ            .= min(∂x^2,∂y^2)./(1.0.+❄.^n)./4.1./4.0
		Res❄          .= (.-diff(qx, dims=1)./∂x .-diff(qy, dims=2)./∂y .+ b) .+ γ.*Res❄
		❄             .= max.(0.0, ❄ .+ ∂τ.*Res❄)
		Err❄ .-= ❄
		if mod(iτ,nout)==0 && viz
			p1 = heatmap(xc, yc, ❄', xlabel="lx", ylabel="ly", title="shallow ice, iτ=$iτ", aspect_ratio=1, xlims=(xc[1], xc[end]), ylims=(yc[1], yc[end]), c=:viridis)
			p2 = plot(xc, ❄[:,Int(round(ny/2))], xlabel="lx", ylabel="height", xlims=(xc[1], xc[end]), ylims=(0., 1.), legend=false, framestyle=:box, aspect_ratio=2.5)
			l = @layout [ a{0.8h} ; b{0.2h} ]
			display(plot(p1, p2, layout=l ))
		end
		if mod(iτ,nout)==0 maxErr=maximum(abs.(Err❄)); @printf("iτ=%d, max(err)=%1.2e \n", iτ, maxErr); if maxErr<ε global iτg=iτ; break; end end
	end
	@printf("T_eff = %1.2e GB/s \n", (2/1e9*nx*ny*sizeof(lx))/((Base.time()-t0)/iτg))
	return
end

@time sia_2D()

ext: jl


text:
value(x) = x
cuify(x) = error("To use LinSolveGPUFactorize, you must do `using CuArrays`")
promote_u0(u0,p,t0) = u0
promote_tspan(u0,p,tspan,prob,kwargs) = tspan
get_tmp(x) = nothing
isdistribution(u0) = false

function SciMLBase.tmap(args...)
  error("Zygote must be added to differentiate Zygote? If you see this error, report it.")
end

const IS_OPENBLAS = Ref(true)

function __init__()
  @static if VERSION < v"1.7beta"
    blas = BLAS.vendor()
    IS_OPENBLAS[] = blas == :openblas64 || blas == :openblas
  else
    IS_OPENBLAS[] = occursin("openblas", BLAS.get_config().loaded_libs[1].libname)
  end
  @require ApproxFun="28f2ccd6-bb30-5033-b560-165f7b14dc2f" begin
    eval_u0(u0::ApproxFun.Fun) = false
  end

  @require Measurements="eff96d63-e80a-5855-80a2-b1b0885c5ab7" begin

    promote_u0(u0::AbstractArray{<:Measurements.Measurement},p::AbstractArray{<:Measurements.Measurement},t0) = u0
    promote_u0(u0,p::AbstractArray{<:Measurements.Measurement},t0) = eltype(p).(u0)

    value(x::Type{Measurements.Measurement{T}}) where {T} = T
    value(x::Measurements.Measurement) = Measurements.value(x)

    @inline fastpow(x::Measurements.Measurement, y::Measurements.Measurement) = x^y

    # Support adaptive steps should be errorless
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:Measurements.Measurement,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::Array{<:Measurements.Measurement,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline ODE_DEFAULT_NORM(u::Measurements.Measurement,t) = abs(Measurements.value(u))
  end

  @require MonteCarloMeasurements="0987c9cc-fe09-11e8-30f0-b96dd679fdca" begin

    promote_u0(u0::AbstractArray{<:MonteCarloMeasurements.AbstractParticles},p::AbstractArray{<:MonteCarloMeasurements.AbstractParticles},t0) = u0
    promote_u0(u0,p::AbstractArray{<:MonteCarloMeasurements.AbstractParticles},t0) = eltype(p).(u0)

    value(x::Type{MonteCarloMeasurements.AbstractParticles{T,N}}) where {T,N} = T
    value(x::MonteCarloMeasurements.AbstractParticles) = mean(x.particles)

    @inline fastpow(x::MonteCarloMeasurements.AbstractParticles, y::MonteCarloMeasurements.AbstractParticles) = x^y

    # Support adaptive steps should be errorless
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:MonteCarloMeasurements.AbstractParticles,N},t) where {N}
      sqrt(mean(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))))
    end
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:MonteCarloMeasurements.AbstractParticles,N},t::AbstractArray{<:MonteCarloMeasurements.AbstractParticles,N}) where {N}
      sqrt(mean(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(value.(t)))))
    end
    @inline ODE_DEFAULT_NORM(u::MonteCarloMeasurements.AbstractParticles,t) = abs(value(u))
  end

  @require Unitful="1986cc42-f94f-5a68-af5c-568840ba703d" begin
    # Support adaptive errors should be errorless for exponentiation
    value(x::Type{Unitful.AbstractQuantity{T,D,U}}) where {T,D,U} = T
    value(x::Unitful.AbstractQuantity) = x.val
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:Unitful.AbstractQuantity,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::Array{<:Unitful.AbstractQuantity,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline ODE_DEFAULT_NORM(u::Unitful.AbstractQuantity,t) = abs(value(u))
    @inline UNITLESS_ABS2(x::Unitful.AbstractQuantity) = real(abs2(x)/oneunit(x)*oneunit(x))
  end

  @require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin
    include("tracker.jl")
  end

  @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin
    cuify(x::AbstractArray) = CUDA.CuArray(x)
    default_factorize(A::CUDA.CuArray) = qr(A)

    ODE_DEFAULT_NORM(u::CUDA.CuArray,t) = sqrt(real(sum(abs2,u))/length(u))

    @inline function ODE_DEFAULT_NORM(u::CUDA.CuArray{<:ForwardDiff.Dual},t)
      sqrt(sum(abs2,value.(u)) / length(u))
    end

    @inline function ODE_DEFAULT_NORM(u::CUDA.CuArray{<:ForwardDiff.Dual},t::ForwardDiff.Dual)
      sqrt(sum(abs2,u) / length(u))
    end
  end

  @require ReverseDiff="37e2e3b7-166d-5795-8a7a-e32c996b4267" begin
    include("reversediff.jl")
  end

  @require Zygote="e88e6eb3-aa80-5325-afca-941959d7151f" begin
    function ∇tmap(cx, f, args...)
      ys_and_backs = SciMLBase.tmap((args...) -> Zygote._pullback(cx, f, args...), args...)
      if isempty(ys_and_backs)
        ys_and_backs, _ -> (NoTangent(),NoTangent())
      else
        ys, backs = Zygote.unzip(ys_and_backs)
        function ∇tmap_internal(Δ)
          Δf_and_args_zipped = SciMLBase.tmap((f, δ) -> f(δ), backs, Δ)
          Δf_and_args = Zygote.unzip(Δf_and_args_zipped)
          Δf = reduce(Zygote.accum, Δf_and_args[1])
          (Δf, Δf_and_args[2:end]...)
        end
        ys,∇tmap_internal
      end
    end

    function ∇responsible_map(cx, f, args...)
      ys_and_backs = SciMLBase.responsible_map((args...) -> Zygote._pullback(cx, f, args...), args...)
      if isempty(ys_and_backs)
        ys_and_backs, _ -> (NoTangent(),NoTangent())
      else
        ys, backs = Zygote.unzip(ys_and_backs)
        ys, function ∇responsible_map_internal(Δ)
          # Apply pullbacks in reverse order. Needed for correctness if `f` is stateful.
          Δf_and_args_zipped = SciMLBase.responsible_map((f, δ) -> f(δ), Zygote._tryreverse(SciMLBase.responsible_map, backs, Δ)...)
          Δf_and_args = Zygote.unzip(Zygote._tryreverse(SciMLBase.responsible_map, Δf_and_args_zipped))
          Δf = reduce(Zygote.accum, Δf_and_args[1])
          (Δf, Δf_and_args[2:end]...)
        end
      end
    end
  end

  @require GeneralizedGenerated="6b9d7cbe-bcb9-11e9-073f-15a7a543e2eb" begin
    SciMLBase.numargs(::GeneralizedGenerated.RuntimeFn{Args}) where Args = GeneralizedGenerated.from_type(Args) |> length
  end
end


text:
# ---
# title: 133. Clone Graph
# id: problem133
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search, Breadth-first Search, Graph
# link: <https://leetcode.com/problems/clone-graph/description/>
# hidden: true
# ---
#
# Given a reference of a node in a
# **[connected](https://en.wikipedia.org/wiki/Connectivity_\(graph_theory\)#Connected_graph)**
# undirected graph.
#
# Return a [**deep
# copy**](https://en.wikipedia.org/wiki/Object_copying#Deep_copy) (clone) of the
# graph.
#
# Each node in the graph contains a val (`int`) and a list (`List[Node]`) of its
# neighbors.
#
#
#
#     class Node {
#         public int val;
#         public List<Node> neighbors;
#     }
#
#
#
#
# **Test case format:**
#
# For simplicity sake, each node's value is the same as the node's index
# (1-indexed). For example, the first node with `val = 1`, the second node with
# `val = 2`, and so on. The graph is represented in the test case using an
# adjacency list.
#
# **Adjacency list**  is a collection of unordered  **lists**  used to represent
# a finite graph. Each list describes the set of neighbors of a node in the
# graph.
#
# The given node will always be the first node with `val = 1`. You must return
# the **copy of the given node** as a reference to the cloned graph.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/11/04/133_clone_graph_question.png)
#
#
#
#     Input: adjList = [[2,4],[1,3],[2,4],[1,3]]
#     Output: [[2,4],[1,3],[2,4],[1,3]]
#     Explanation: There are 4 nodes in the graph.
#     1st node (val = 1)'s neighbors are 2nd node (val = 2) and 4th node (val = 4).
#     2nd node (val = 2)'s neighbors are 1st node (val = 1) and 3rd node (val = 3).
#     3rd node (val = 3)'s neighbors are 2nd node (val = 2) and 4th node (val = 4).
#     4th node (val = 4)'s neighbors are 1st node (val = 1) and 3rd node (val = 3).
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/01/07/graph.png)
#
#
#
#     Input: adjList = [[]]
#     Output: [[]]
#     Explanation: Note that the input contains one empty list. The graph consists of only one node with val = 1 and it does not have any neighbors.
#
#
# **Example 3:**
#
#
#
#     Input: adjList = []
#     Output: []
#     Explanation: This an empty graph, it does not have any nodes.
#
#
# **Example 4:**
#
# ![](https://assets.leetcode.com/uploads/2020/01/07/graph-1.png)
#
#
#
#     Input: adjList = [[2],[1]]
#     Output: [[2],[1]]
#
#
#
#
# **Constraints:**
#
#   * `1 <= Node.val <= 100`
#   * `Node.val` is unique for each node.
#   * Number of Nodes will not exceed 100.
#   * There is no repeated edges and no self-loops in the graph.
#   * The Graph is connected and all nodes can be visited starting from the given node.
#
#
## @lc code=start
using LeetCode

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


text:
using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils
using ..MainLoopKernel, ..WorkQueueUtils, ..Housdorff
using CUDA, Logging,..CUDAGpuUtils, ..ResultListUtils,..WorkQueueUtils,..ScanForDuplicates, Logging,StaticArrays, ..IterationUtils, ..ReductionUtils, ..CUDAAtomicUtils,..MetaDataUtils

mainArrDims= (67,177,90);
dataBdim = (32,32,32);
robustnessPercent= 0.9;
numberToLooFor= 2;
goldGPU,segmGPU= CUDA.zeros(mainArrDims),CUDA.zeros(mainArrDims);

boolKernelArgs, mainKernelArgs,threadsBoolKern,blocksBoolKern ,threadsMainKern,blocksMainKern ,shmemSizeBool,shmemSizeMain= preparehousedorfKernel(goldGPU,segmGPU,robustnessPercent,numberToLooFor)

function addToWorkQueueKernel(dilatationArrsA,dilatationArrsB, mainArrDims,dataBdim ,metaDataDims,metaData,iterThrougWarNumb,robustnessPercent   ,shmemSumLengthMaxDiv4,globalFpResOffsetCounter,globalFnResOffsetCounter   ,globalIterationNumber,globalCurrentFnCount,globalCurrentFpCount   ,globalIterationNumb,workQueueEEE,workQueueEEEcounter,workQueueEEO,workQueueEEOcounter    ,workQueueEOE,workQueueEOEcounter,workQueueOEE,workQueueOEEcounter   ,workQueueOOE,workQueueOOEcounter,workQueueEOO,workQueueEOOcounter   ,workQueueOEO,workQueueOEOcounter,workQueueOOO,workQueueOOOcounter   ,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed   ,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY   ,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength, fn,fp,resList,inBlockLoopXZIterWithPadding,paddingStore,shmemblockDataLenght,shmemblockDataLoop)

    if (blockIdxX()==1)
      @ifXY 1 1 WorkQueueUtils.@appendToWorkQueue(1,1,1,0)

    end
      return
  end
  @cuda threads=(32,32) blocks=20 addToWorkQueueKernel(dilatationArrsA,dilatationArrsB, mainArrDims,dataBdim ,metaDataDims,metaData,iterThrougWarNumb,robustnessPercent   ,shmemSumLengthMaxDiv4,globalFpResOffsetCounter,globalFnResOffsetCounter   ,globalIterationNumber,globalCurrentFnCount,globalCurrentFpCount   ,globalIterationNumb,workQueueEEE,workQueueEEEcounter,workQueueEEO,workQueueEEOcounter    ,workQueueEOE,workQueueEOEcounter,workQueueOEE,workQueueOEEcounter   ,workQueueOOE,workQueueOOEcounter,workQueueEOO,workQueueEOOcounter   ,workQueueOEO,workQueueOEOcounter,workQueueOOO,workQueueOOOcounter   ,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed   ,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY   ,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength, fn,fp,resList,inBlockLoopXZIterWithPadding,paddingStore,shmemblockDataLenght,shmemblockDataLoop)

  @test workQueueEEEcounter[1]==1
  @test workQueueEEOcounter[1]==1
  @test workQueueEOEcounter[1]==1
  @test workQueueOEEcounter[1]==1
  @test workQueueOOEcounter[1]==1
  @test workQueueEOOcounter[1]==1
  @test workQueueOEOcounter[1]==1
  @test workQueueOOOcounter[1]==1

Int64(workQueueEEOcounter[1])
Int64(workQueueEOOcounter[1])


text:
@testset "BLAS" begin
    @testset "dot" begin
        @testset "all entries" begin
            n = 10
            x, y = randn(n), randn(n)
            ẋ, ẏ = randn(n), randn(n)
            x̄, ȳ = randn(n), randn(n)
            frule_test(BLAS.dot, (x, ẋ), (y, ẏ))
            rrule_test(BLAS.dot, randn(), (x, x̄), (y, ȳ))
        end

        @testset "over strides" begin
            n = 10
            incx = 2
            incy = 3
            x, y = randn(n * incx), randn(n * incy)
            x̄, ȳ = randn(n * incx), randn(n * incy)
            rrule_test(
                BLAS.dot,
                randn(),
                (n, nothing),
                (x, x̄),
                (incx, nothing),
                (y, ȳ),
                (incy, nothing),
            )
        end
    end

    @testset "nrm2" begin
        @testset "all entries" begin
            @testset "$T" for T in (Float64,ComplexF64)
                n = 10
                x, ẋ, x̄ = randn(T, n), randn(T, n), randn(T, n)
                frule_test(BLAS.nrm2, (x, ẋ))
                rrule_test(BLAS.nrm2, randn(), (x, x̄))
            end
        end

        @testset "over strides" begin
            dims = (3, 2, 1)
            incx = 2
            @testset "Array{$T,$N}" for N in 1:length(dims), T in (Float64,ComplexF64)
                s = (dims[1] * incx, dims[2:N]...)
                n = div(prod(s), incx)
                x, x̄ = randn(T, s), randn(T, s)
                rrule_test(
                    BLAS.nrm2,
                    randn(),
                    (n, nothing),
                    (x, x̄),
                    (incx, nothing);
                    atol=0,
                    rtol=1e-5,
                )
            end
        end
    end

    @testset "asum" begin
        @testset "all entries" begin
            @testset "$T" for T in (Float64,ComplexF64)
                n = 6
                x, ẋ, x̄ = randn(T, n), randn(T, n), randn(T, n)
                frule_test(BLAS.asum, (x, ẋ))
                rrule_test(BLAS.asum, randn(), (x, x̄))
            end
        end

        @testset "over strides" begin
            dims = (2, 2, 1)
            incx = 2
            @testset "Array{$T,$N}" for N in 1:length(dims), T in (Float64,ComplexF64)
                s = (dims[1] * incx, dims[2:N]...)
                n = div(prod(s), incx)
                x, x̄ = randn(T, s), randn(T, s)
                rrule_test( BLAS.asum, randn(), (n, nothing), (x, x̄), (incx, nothing))
            end
        end
    end

    @testset "gemm" begin
        dims = 3:5
        for m in dims, n in dims, p in dims, tA in ('N', 'C', 'T'), tB in ('N', 'C', 'T'), T in (Float64, ComplexF64)
            α = randn(T)
            A = randn(T, tA === 'N' ? (m, n) : (n, m))
            B = randn(T, tB === 'N' ? (n, p) : (p, n))
            C = gemm(tA, tB, α, A, B)
            ȳ = randn(T, size(C)...)
            rrule_test(
                gemm,
                ȳ,
                (tA, nothing),
                (tB, nothing),
                (α, randn(T)),
                (A, randn(T, size(A))),
                (B, randn(T, size(B)));
                check_inferred=false,
            )

            rrule_test(
                gemm,
                ȳ,
                (tA, nothing),
                (tB, nothing),
                (A, randn(T, size(A))),
                (B, randn(T, size(B)));
                check_inferred=false,
            )
        end
    end

    @testset "gemv" begin
        for n in 3:5, m in 3:5, t in ('N', 'C', 'T'), T in (Float64, ComplexF64)
            α = randn(T)
            A = randn(T, m, n)
            x = randn(T, t === 'N' ? n : m)
            y = gemv(t, α, A, x)
            ȳ = randn(T, size(y)...)
            rrule_test(
                gemv,
                ȳ,
                (t, nothing),
                (α, randn(T)),
                (A, randn(T, size(A))),
                (x, randn(T, size(x)));
                check_inferred=false,
            )
        end
    end
end


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using ContactImplicitMPC
using LinearAlgebra

# ## Simulation
s = get_simulation("walledcartpole", "flat_2D_lc", "flat");
model = s.model
env = s.env

# ## Reference Trajectory
ref_traj = contact_trajectory(model, env, H, h);
qref = [0.0; 0.0; 0.0; 0.0]
ur = zeros(model.nu)
γr = zeros(model.nc)
br = zeros(model.nc * friction_dim(env))
ψr = zeros(model.nc)
ηr = zeros(model.nc * friction_dim(env))
wr = zeros(model.nw)
h = 0.04
H = 50

# ## Set Reference
for t = 1:H
	ref_traj.z[t] = pack_z(model, env, qref, γr, br, ψr, ηr)
	ref_traj.θ[t] = pack_θ(model, qref, qref, ur, wr, model.μ_world, ref_traj.h)
end

# ## Simulation Time
h_sim = h / N_sample
H_sim = 1000

# ## No Policy
p = open_loop_policy(fill(zeros(nu), H_sim); N_sample = N_sample)

# ## MPC setup
N_sample = 2
H_mpc = 10
κ_mpc = 2.0e-4

obj = TrackingVelocityObjective(model, env, H_mpc,
	q = [[Diagonal([1e-1;1e-3;1e-8;1e-8]) for t = 1:H_mpc-1];
		 [Diagonal([1e+1;1e+0;1e-8;1e-8]) for t = 1:1]],
	v = [[Diagonal([1e-0;3e+1;1e-8;1e-8]) for t = 1:H_mpc-1];
		 [Diagonal([1e+1;1e+1;1e-8;1e-8]) for t = 1:1]],
	u = [Diagonal([3e-2]) for t = 1:H_mpc])

p = ci_mpc_policy(ref_traj, s, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
    mpc_opts =CIMPCOptions(),
	n_opts = NewtonOptions(
		r_tol = 3e-4,
		max_iter = 5,
		));

# ## Disturbances
idx_d1 = 20
idx_d2 = idx_d1 + 200
idx_d3 = idx_d2 + 150
idx_d4 = idx_d3 + 200
idx_d5 = idx_d4 + 150
idx = [idx_d1, idx_d2, idx_d3, idx_d4, idx_d5]
impulses = [[+0.2; 0; 0; 0], [+0.2; 0; 0; 0], [-0.2; 0; 0; 0], [+0.2; 0; 0; 0], [+0.2; 0; 0; 0]]
d = impulse_disturbances(impulses, idx);

# ## Initial Conditions
q1_sim = zeros(4)
v1_sim = zeros(4)

# ## Simulator
sim = simulator(s, H_sim, h=h_sim, policy=p, dist=d)

# ## Simulate
status = simulate!(sim, q1_sim, v1_sim);

# ## Visualizer
vis = ContactImplicitMPC.Visualizer()
ContactImplicitMPC.render(vis)

# ## Visualize
anim = visualize_robot!(vis, model, sim.traj, h=h_sim, sample = 1)

# ## Timing result
# Julia is [JIT-ed](https://en.wikipedia.org/wiki/Just-in-time_compilation) so re-run the MPC setup through Simulate for correct timing results.
process!(sim.stats, N_sample) # Time budget
H_sim * h_sim / sum(sim.stats.policy_time) # Speed ratio


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: codestar_notifications
using AWS.Compat
using AWS.UUIDs

"""
    create_notification_rule(detail_type, event_type_ids, name, resource, targets)
    create_notification_rule(detail_type, event_type_ids, name, resource, targets, params::Dict{String,<:Any})

Creates a notification rule for a resource. The rule specifies the events you want
notifications about and the targets (such as SNS topics) where you want to receive them.

# Arguments
- `detail_type`: The level of detail to include in the notifications for this resource.
  BASIC will include only the contents of the event as it would appear in AWS CloudWatch.
  FULL will include any supplemental information provided by AWS CodeStar Notifications
  and/or the service for the resource for which the notification is created.
- `event_type_ids`: A list of event types associated with this notification rule. For a
  list of allowed events, see EventTypeSummary.
- `name`: The name for the notification rule. Notifictaion rule names must be unique in
  your AWS account.
- `resource`: The Amazon Resource Name (ARN) of the resource to associate with the
  notification rule. Supported resources include pipelines in AWS CodePipeline, repositories
  in AWS CodeCommit, and build projects in AWS CodeBuild.
- `targets`: A list of Amazon Resource Names (ARNs) of SNS topics to associate with the
  notification rule.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ClientRequestToken"`: A unique, client-generated idempotency token that, when provided
  in a request, ensures the request cannot be repeated with a changed parameter. If a request
  with the same parameters is received and a token is included, the request returns
  information about the initial request that used that token.  The AWS SDKs prepopulate
  client request tokens. If you are using an AWS SDK, an idempotency token is created for
  you.
- `"Status"`: The status of the notification rule. The default value is ENABLED. If the
  status is set to DISABLED, notifications aren't sent for the notification rule.
- `"Tags"`: A list of tags to apply to this notification rule. Key names cannot start with
  \"aws\".
"""
function create_notification_rule(
    DetailType,
    EventTypeIds,
    Name,
    Resource,
    Targets;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return codestar_notifications(
        "POST",
        "/createNotificationRule",
        Dict{String,Any}(
            "DetailType" => DetailType,
            "EventTypeIds" => EventTypeIds,
            "Name" => Name,
            "Resource" => Resource,
            "Targets" => Targets,
            "ClientRequestToken" => string(uuid4()),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_notification_rule(
    DetailType,
    EventTypeIds,
    Name,
    Resource,
    Targets,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return codestar_notifications(
        "POST",
        "/createNotificationRule",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "DetailType" => DetailType,
                    "EventTypeIds" => EventTypeIds,
                    "Name" => Name,
                    "Resource" => Resource,
                    "Targets" => Targets,
                    "ClientRequestToken" => string(uuid4()),
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_notification_rule(arn)
    delete_notification_rule(arn, params::Dict{String,<:Any})

Deletes a notification rule for a resource.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the notification rule you want to delete.

"""
function delete_notification_rule(Arn; aws_config::AbstractAWSConfig=global_aws_config())
    return codestar_notifications(
        "POST",
        "/deleteNotificationRule",
        Dict{String,Any}("Arn" => Arn);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_notification_rule(
    Arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return codestar_notifications(
        "POST",
        "/deleteNotificationRule",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("Arn" => Arn), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_target(target_address)
    delete_target(target_address, params::Dict{String,<:Any})

Deletes a specified target for notifications.

# Arguments
- `target_address`: The Amazon Resource Name (ARN) of the SNS topic to delete.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ForceUnsubscribeAll"`: A Boolean value that can be used to delete all associations with
  this SNS topic. The default value is FALSE. If set to TRUE, all associations between that
  target and every notification rule in your AWS account are deleted.
"""
function delete_target(TargetAddress; aws_config::AbstractAWSConfig=global_aws_config())
    return codestar_notifications(
        "POST",
        "/deleteTarget",
        Dict{String,Any}("TargetAddress" => TargetAddress);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_target(
    TargetAddress,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return codestar_notifications(
        "POST",
        "/deleteTarget",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("TargetAddress" => TargetAddress), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_notification_rule(arn)
    describe_notification_rule(arn, params::Dict{String,<:Any})

Returns information about a specified notification rule.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the notification rule.

"""
function describe_notification_rule(Arn; aws_config::AbstractAWSConfig=global_aws_config())
    return codestar_notifications(
        "POST",
        "/describeNotificationRule",
        Dict{String,Any}("Arn" => Arn);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_notification_rule(
    Arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return codestar_notifications(
        "POST",
        "/describeNotificationRule",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("Arn" => Arn), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_event_types()
    list_event_types(params::Dict{String,<:Any})

Returns information about the event types available for configuring notifications.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Filters"`: The filters to use to return information by service or resource type.
- `"MaxResults"`: A non-negative integer used to limit the number of returned results. The
  default number is 50. The maximum number of results that can be returned is 100.
- `"NextToken"`: An enumeration token that, when provided in a request, returns the next
  batch of the results.
"""
function list_event_types(; aws_config::AbstractAWSConfig=global_aws_config())
    return codestar_notifications(
        "POST", "/listEventTypes"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_event_types(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return codestar_notifications(
        "POST",
        "/listEventTypes",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_notification_rules()
    list_notification_rules(params::Dict{String,<:Any})

Returns a list of the notification rules for an AWS account.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Filters"`: The filters to use to return information by service or resource type. For
  valid values, see ListNotificationRulesFilter.  A filter with the same name can appear more
  than once when used with OR statements. Filters with different names should be applied with
  AND statements.
- `"MaxResults"`: A non-negative integer used to limit the number of returned results. The
  maximum number of results that can be returned is 100.
- `"NextToken"`: An enumeration token that, when provided in a request, returns the next
  batch of the results.
"""
function list_notification_rules(; aws_config::AbstractAWSConfig=global_aws_config())
    return codestar_notifications(
        "POST",
        "/listNotificationRules";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_notification_rules(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return codestar_notifications(
        "POST",
        "/listNotificationRules",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_tags_for_resource(arn)
    list_tags_for_resource(arn, params::Dict{String,<:Any})

Returns a list of the tags associated with a notification rule.

# Arguments
- `arn`: The Amazon Resource Name (ARN) for the notification rule.

"""
function list_tags_for_resource(Arn; aws_config::AbstractAWSConfig=global_aws_config())
    return codestar_notifications(
        "POST",
        "/listTagsForResource",
        Dict{String,Any}("Arn" => Arn);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_tags_for_resource(
    Arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return codestar_notifications(
        "POST",
        "/listTagsForResource",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("Arn" => Arn), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_targets()
    list_targets(params::Dict{String,<:Any})

Returns a list of the notification rule targets for an AWS account.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Filters"`: The filters to use to return information by service or resource type. Valid
  filters include target type, target address, and target status.  A filter with the same
  name can appear more than once when used with OR statements. Filters with different names
  should be applied with AND statements.
- `"MaxResults"`: A non-negative integer used to limit the number of returned results. The
  maximum number of results that can be returned is 100.
- `"NextToken"`: An enumeration token that, when provided in a request, returns the next
  batch of the results.
"""
function list_targets(; aws_config::AbstractAWSConfig=global_aws_config())
    return codestar_notifications(
        "POST", "/listTargets"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_targets(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return codestar_notifications(
        "POST",
        "/listTargets",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    subscribe(arn, target)
    subscribe(arn, target, params::Dict{String,<:Any})

Creates an association between a notification rule and an SNS topic so that the associated
target can receive notifications when the events described in the rule are triggered.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the notification rule for which you want to
  create the association.
- `target`:

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ClientRequestToken"`: An enumeration token that, when provided in a request, returns
  the next batch of the results.
"""
function subscribe(Arn, Target; aws_config::AbstractAWSConfig=global_aws_config())
    return codestar_notifications(
        "POST",
        "/subscribe",
        Dict{String,Any}("Arn" => Arn, "Target" => Target);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function subscribe(
    Arn,
    Target,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return codestar_notifications(
        "POST",
        "/subscribe",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("Arn" => Arn, "Target" => Target), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    tag_resource(arn, tags)
    tag_resource(arn, tags, params::Dict{String,<:Any})

Associates a set of provided tags with a notification rule.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the notification rule to tag.
- `tags`: The list of tags to associate with the resource. Tag key names cannot start with
  \"aws\".

"""
function tag_resource(Arn, Tags; aws_config::AbstractAWSConfig=global_aws_config())
    return codestar_notifications(
        "POST",
        "/tagResource",
        Dict{String,Any}("Arn" => Arn, "Tags" => Tags);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function tag_resource(
    Arn,
    Tags,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return codestar_notifications(
        "POST",
        "/tagResource",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("Arn" => Arn, "Tags" => Tags), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    unsubscribe(arn, target_address)
    unsubscribe(arn, target_address, params::Dict{String,<:Any})

Removes an association between a notification rule and an Amazon SNS topic so that
subscribers to that topic stop receiving notifications when the events described in the
rule are triggered.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the notification rule.
- `target_address`: The ARN of the SNS topic to unsubscribe from the notification rule.

"""
function unsubscribe(Arn, TargetAddress; aws_config::AbstractAWSConfig=global_aws_config())
    return codestar_notifications(
        "POST",
        "/unsubscribe",
        Dict{String,Any}("Arn" => Arn, "TargetAddress" => TargetAddress);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function unsubscribe(
    Arn,
    TargetAddress,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return codestar_notifications(
        "POST",
        "/unsubscribe",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("Arn" => Arn, "TargetAddress" => TargetAddress),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    untag_resource(arn, tag_keys)
    untag_resource(arn, tag_keys, params::Dict{String,<:Any})

Removes the association between one or more provided tags and a notification rule.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the notification rule from which to remove the
  tags.
- `tag_keys`: The key names of the tags to remove.

"""
function untag_resource(Arn, TagKeys; aws_config::AbstractAWSConfig=global_aws_config())
    return codestar_notifications(
        "POST",
        "/untagResource",
        Dict{String,Any}("Arn" => Arn, "TagKeys" => TagKeys);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function untag_resource(
    Arn,
    TagKeys,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return codestar_notifications(
        "POST",
        "/untagResource",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("Arn" => Arn, "TagKeys" => TagKeys), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_notification_rule(arn)
    update_notification_rule(arn, params::Dict{String,<:Any})

Updates a notification rule for a resource. You can change the events that trigger the
notification rule, the status of the rule, and the targets that receive the notifications.
To add or remove tags for a notification rule, you must use TagResource and UntagResource.

# Arguments
- `arn`: The Amazon Resource Name (ARN) of the notification rule.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"DetailType"`: The level of detail to include in the notifications for this resource.
  BASIC will include only the contents of the event as it would appear in AWS CloudWatch.
  FULL will include any supplemental information provided by AWS CodeStar Notifications
  and/or the service for the resource for which the notification is created.
- `"EventTypeIds"`: A list of event types associated with this notification rule.
- `"Name"`: The name of the notification rule.
- `"Status"`: The status of the notification rule. Valid statuses include enabled (sending
  notifications) or disabled (not sending notifications).
- `"Targets"`: The address and type of the targets to receive notifications from this
  notification rule.
"""
function update_notification_rule(Arn; aws_config::AbstractAWSConfig=global_aws_config())
    return codestar_notifications(
        "POST",
        "/updateNotificationRule",
        Dict{String,Any}("Arn" => Arn);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_notification_rule(
    Arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return codestar_notifications(
        "POST",
        "/updateNotificationRule",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("Arn" => Arn), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
random_particles(particle_medium::PhysicalMedium{T,Dim}, particle_shape::Shape{T,Dim}; kws...) where {T<:AbstractFloat,Dim} = random_particles(particle_medium, [particle_shape]; kws...)

function random_particles(particle_medium::PhysicalMedium{T,Dim}, particle_shapes::Vector{S};
        num_particles::Int = length(particle_shapes), volume_fraction::T = zero(T),
        region_shape::Shape{T,Dim} = Box(zeros(T,Dim), ones(T,Dim) .* T(10)*sum(outer_radius.(particle_shapes))),
        current_particles::Vector{AbstractParticle{T,Dim}} = AbstractParticle{T,Dim}[],
        kws...
) where {T<:AbstractFloat,Dim, S<:Shape{T,Dim}}

    particles = if volume_fraction == zero(T)
        num_each = Int(round(num_particles / length(particle_shapes)))
        for s in particle_shapes
            add_particles = random_particles(particle_medium, s, region_shape, num_each; current_particles = current_particles, kws...)
            current_particles = [current_particles; add_particles]
        end
        current_particles
    else
        vol_each = volume_fraction / length(particle_shapes)
        for s in particle_shapes
            add_particles = random_particles(particle_medium, s, region_shape, vol_each; current_particles = current_particles, kws...)
            current_particles = [current_particles; add_particles]
        end
        current_particles
    end

    return particles
end

# random_particles(particle_medium, particle_shapes)

#
# a = 1
# function f1(; a = [])
#     if volume_fraction != zero(T)
#         volume_fraction_each = volume_fraction / length(particle_shapes)
#         for s in particle_shapes
#             b = a
#             a = [a; b]
#             println(a)
#         end
#         println(a)
#     end
#     println(a)
# end
# f1(; a = [2])
# a

# function random_particles(particle_medium::PhysicalMedium{T,Dim}, particle_shape::Shape{T,Dim};
#         region_shape::Shape{T,Dim} = Rectangle(zeros(T,2), T(10)*outer_radius(particle_shape), T(10)*outer_radius(particle_shape)),
#         num_particles::Int = 0, volume_fraction::T = zero(T), kws...) where {T<:AbstractFloat,Dim}
#
#     if volume_fraction == zero(T)
#         if num_particles == 0 num_particles = 5 end
#         random_particles(particle_medium, particle_shape,
#             region_shape, num_particles; kws...)
#     else
#         random_particles(particle_medium, particle_shape,
#             region_shape, volume_fraction; kws...)
#     end
#
# end

function random_particles(particle_medium::P, particle_shape::S,
    region_shape::Shape{T,Dim}, volfrac::AbstractFloat; kws...
) where {T,Dim,P<:PhysicalMedium{T,Dim},S<:Shape{T,Dim}}

    N = Int(round(volfrac * volume(region_shape) / volume(particle_shape)))
    return random_particles(particle_medium, particle_shape, region_shape, N; kws...)
end

"""
    random_particles(particle_medium, particle_shapes::Vector{Shape}, region_shape, volume_fraction::Number;
        seed=Random.make_seed())
    random_particles(particle_medium, particle_shape::Shape, region_shape, volume_fraction::Number;
        seed=Random.make_seed())
    random_particles(particle_medium, particle_shape::Shape, region_shape, N::Integer;
        seed=Random.make_seed())

Generate `N` random particles that fit inside `region_shape` (or fill with `volume_fraction`)

Specify seed to make output deterministic. Algorithm places particles unifomly randomly inside
the bounding box of `region_shape` and discards particle if it overlaps (based on outer radius)
or does not lies completely in box.

When passing particle_shapes::Vector{Shape} we assume each element is equally likely to occur. Repeating the same shape will lead to it being placed more often.
"""
function random_particles(particle_medium::P, particle_shape::S, region_shape::Shape{T,Dim}, N::Integer;
        seed=Random.make_seed(),
        verbose::Bool = false,
        separation_ratio::T = T(1.005), # Min distance between particle centres relative to their outer radiuses.
        max_attempts_to_place_particle::Int = 3000, # Maximum number of attempts to place a particle
        current_particles::Vector{AbstractParticle{T,Dim}} = AbstractParticle{T,Dim}[] # Particles already present.
) where {T<:AbstractFloat,Dim,P<:PhysicalMedium{T,Dim},S<:Shape{T,Dim}}

    # Check volume fraction is not impossible
    volfrac = N * volume(particle_shape) / volume(region_shape)
    max_packing = if length(current_particles) > 0
        0.7854 - sum(volume.(current_particles))/volume(region_shape)
    else 0.7854
    end
    if volfrac  > max_packing
        error("Specified volume fraction is larger than optimal packing of circles.")
    end

    # Create pseudorandom device with specified seed
    randgen = MersenneTwister(seed)

    box = bounding_box(region_shape)

    if verbose
        @printf("""\n
            Generating %d randomly positioned %s shaped particles
            Total particle volume: %0.5g
            Inside %s of volume: %0.5g
            Particle volume fraction: %0.5g
            Bounding box volume: %0.5g
            """, N, name(particle_shape), N*volume(particle_shape), name(region_shape),
            volume(region_shape), volfrac, volume(box)
        )
    end

    # Allocate memory for particles
    L = length(current_particles)
    particles = Vector{Particle{T,Dim,P,S}}(undef, N + L)
    particles[1:L] = current_particles

    for n = L+1:L+N

        num_attempts = 0
        overlapping = true
        while overlapping

            outside_box = true
            while outside_box
                x = (box.dimensions ./ 2) .* (1 .- 2 .* rand(randgen,T,Dim)) + origin(box)
                particles[n] = Particle(particle_medium, congruent(particle_shape, x))
                outside_box = !(particles[n] ⊆ region_shape)
            end

            overlapping = false
            for i = 1:(n-1) #compare with previous particles
                if norm(origin(particles[n]) - origin(particles[i])) < separation_ratio * (outer_radius(particles[n]) + outer_radius(particles[i]))
                    overlapping = true
                    break
                end
            end

            num_attempts += 1
            if num_attempts > max_attempts_to_place_particle
                error("Tried to place a scatterer $max_attempts_to_place_particle times unsuccessfully. You could try increasing max_attempts_to_place_particle")
            end

        end
    end

    return particles[L+1:L+N]
end


text:
function energy_level(mol_idx, v, J)
    E = 0.0
    for k = 0:4
        tmpE1 = pow(v+0.5, k)
        for l =0:4
            tmpE2 = tmpE1 * Y[k,l] * pow(J*(J+1),l)
            E = E + tmpE2
        end
    end
    return E; #in cm^{-1}
end

max_stars_repo_name: RadiativeTransfer/InelasticScattering

max_issues_repo_name: RadiativeTransfer/InelasticScattering

max_forks_repo_name: RadiativeTransfer/InelasticScattering


text:
# ---
# title: 24. Swap Nodes in Pairs
# id: problem24
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Linked List
# link: <https://leetcode.com/problems/swap-nodes-in-pairs/description/>
# hidden: true
# ---
#
# Given a linked list, swap every two adjacent nodes and return its head.
#
# You may **not** modify the values in the list's nodes. Only nodes itself may
# be changed.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/03/swap_ex1.jpg)
#
#
#
#     Input: head = [1,2,3,4]
#     Output: [2,1,4,3]
#
#
# **Example 2:**
#
#
#
#     Input: head = []
#     Output: []
#
#
# **Example 3:**
#
#
#
#     Input: head = [1]
#     Output: [1]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the list is in the range `[0, 100]`.
#   * `0 <= Node.val <= 100`
#
#
## @lc code=start
using LeetCode

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


text:
using Documenter, Literate
using OneHot

function clear_md_files(dir::String)
    for file in readdir(dir; join=true)
        if endswith(file, ".md")
            rm(file)
        end
    end
end

const literate_dir = joinpath(@__DIR__, "src/literate")
clear_md_files(literate_dir)

for file in readdir(literate_dir; join=true)
    if endswith(file, ".jl")
        Literate.markdown(file, literate_dir)
    end
end

makedocs(
    sitename = "OneHot",
    modules = [OneHot],
    pages = [
        "Home" => "index.md",
        "Examples" => [
            "Matrix multiplication" => "literate/matmul.md"
        ],
        "Reference" => "reference.md"
    ],
    strict = true
)

deploydocs(
    repo = "github.com/cossio/OneHot.jl.git",
    devbranch = "master"
)


text:
using RedBlackGraph
using Test

function lookup_relationships()
    @test lookup_relationship(0, 0) == "self"
    @test lookup_relationship(0, 1) == "parent"
    @test lookup_relationship(0, 2) == "grandparent"
    @test lookup_relationship(0, 3) == "great grandparent"
    @test lookup_relationship(0, 4) == "2nd great grandparent"
    @test lookup_relationship(1, 1) == "sibling"
    @test lookup_relationship(1, 2) == "aunt/uncle"
    @test lookup_relationship(1, 3) == "great aunt/uncle"
    @test lookup_relationship(1, 4) == "2nd great aunt/uncle"
    @test lookup_relationship(2, 1) == "aunt/uncle"
    @test lookup_relationship(2, 2) == "1st cousin"
    @test lookup_relationship(2, 3) == "1st cousin 1 removed"
    @test lookup_relationship(2, 4) == "1st cousin 2 removed"
end

function test_calc_relationships()
    r = red_one(1)
    o = one(r)
    # my daughter, me
    u = [r, 0]
    v = [ 2, o]
    @test lookup_relationship(calculate_relationship(u, v)) == "parent"

    # my daughter, my father
    u = [r, 0]
    v = [4, o]
    @test lookup_relationship(calculate_relationship(u, v)) == "grandparent"

    # my son, my maternal grandmother
    u = [ o,  0]
    v = [ 11, r]
    @test lookup_relationship(calculate_relationship(u, v)) == "great grandparent"

    # 2nd great grandparent
    u = [ o,  0]
    v = [ 16, r]
    @test lookup_relationship(calculate_relationship(u, v)) == "2nd great grandparent"

    # 3rd great grandparent
    u = [ o,  0]
    v = [ 32, r]
    @test lookup_relationship(calculate_relationship(u, v)) == "3rd great grandparent"

    # 4th great grandparent
    u = [ o,  0]
    v = [ 64, r]
    @test lookup_relationship(calculate_relationship(u, v)) == "4th great grandparent"

    # myself, my brother
    u = [ r, 0, 2]
    v = [ 0, r, 2]
    @test lookup_relationship(calculate_relationship(u, v)) == "sibling"

    # myself, my uncle
    u = [ r, 0, 4]
    v = [ 0, r, 2]
    @test lookup_relationship(calculate_relationship(u, v)) == "aunt/uncle"

    # myself, my cousin
    u = [ r, 0, 4]
    v = [ 0, r, 4]
    @test lookup_relationship(calculate_relationship(u, v)) == "1st cousin"

    # myself, my 2nd cousin
    u = [ r, 0, 8]
    v = [ 0, r, 8]
    @test lookup_relationship(calculate_relationship(u, v)) == "2nd cousin"

    # myself, my 3rd cousin
    u = [ r, 0, 16]
    v = [ 0, r, 16]
    @test lookup_relationship(calculate_relationship(u, v)) == "3rd cousin"

    # myself, my 4th cousin
    u = [ r, 0, 32]
    v = [ 0, r, 32]
    @test lookup_relationship(calculate_relationship(u, v)) == "4th cousin"

    # myself, my cousin once remove
    u = [ r, 0, 4]
    v = [ 0, r, 8]
    @test lookup_relationship(calculate_relationship(u, v)) == "1st cousin 1 removed"

    # myself, my 2nd cousin
    u = [ r, 0, 8]
    v = [ 0, r, 16]
    @test lookup_relationship(calculate_relationship(u, v)) == "2nd cousin 1 removed"

    # myself, my 3rd cousin
    u = [ r, 0, 16]
    v = [ 0, r, 32]
    @test lookup_relationship(calculate_relationship(u, v)) == "3rd cousin 1 removed"

    # myself, my 4th cousin
    u = [ r, 0, 32]
    v = [ 0, r, 64]
    @test lookup_relationship(calculate_relationship(u, v)) == "4th cousin 1 removed"

    # no relationship
    u = [ r, 0]
    v = [ 0, r]
    @test lookup_relationship(calculate_relationship(u, v)) == "No Relationship"
end

@testset "relationship" begin
    lookup_relationships()
    test_calc_relationships()
end


text:
int_rules_4_5_2_2 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.5.2.2*(g*sec)^p*(a+b*sec)^m*(c+d*sec)^n =#
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (~b * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n) / (~a * ~f * (2 * ~m + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (EqQ(~m + ~n + 1, 0) && NeQ(2 * ~m + 1, 0))))
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (~b * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n) / (~a * ~f * (2 * ~m + 1)) + ((~m + ~n + 1) / (~a * (2 * ~m + 1))) * Antiderivative(csc(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (ILtQ(~m + ~n + 1, 0) && (NeQ(2 * ~m + 1, 0) && (Not(LtQ(~n, 0)) && Not(IGtQ(~n + 1 / 2, 0) && LtQ(~n + 1 / 2, -((~m + ~n)))))))))
    @apply_utils Antiderivative((csc(~(e') + ~(f') * ~x) * sqrt(~c + ~(d') * csc(~(e') + ~(f') * ~x))) / sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)), ~x) => (~a * ~c * log(1 + (~b / ~a) * csc(~e + ~f * ~x)) * cot(~e + ~f * ~x)) / (~b * ~f * sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * sqrt(~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => (2 * ~a * ~c * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m) / (~b * ~f * (2 * ~m + 1) * sqrt(~c + ~d * csc(~e + ~f * ~x))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ(~m, -1 / 2)))
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (2 * ~a * ~c * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m * (~c + ~d * csc(~e + ~f * ~x)) ^ (~n - 1)) / (~b * ~f * (2 * ~m + 1)) - ((~d * (2 * ~n - 1)) / (~b * (2 * ~m + 1))) * Antiderivative(csc(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * csc(~e + ~f * ~x)) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IGtQ(~n - 1 / 2, 0) && LtQ(~m, -1 / 2))))
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (-(~d) * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m * (~c + ~d * csc(~e + ~f * ~x)) ^ (~n - 1)) / (~f * (~m + ~n)) + ((~c * (2 * ~n - 1)) / (~m + ~n)) * Antiderivative(csc(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m * (~c + ~d * csc(~e + ~f * ~x)) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IGtQ(~n - 1 / 2, 0) && (Not(LtQ(~m, -1 / 2)) && Not(IGtQ(~m - 1 / 2, 0) && LtQ(~m, ~n))))))
    @apply_utils Antiderivative((csc(~(e') + ~(f') * ~x) * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n')) / sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)), ~x) => (-2 * ~d * cot(~e + ~f * ~x) * (~c + ~d * csc(~e + ~f * ~x)) ^ (~n - 1)) / (~f * (2 * ~n - 1) * sqrt(~a + ~b * csc(~e + ~f * ~x))) + ((2 * ~c * (2 * ~n - 1)) / (2 * ~n - 1)) * Antiderivative((csc(~e + ~f * ~x) * (~c + ~d * csc(~e + ~f * ~x)) ^ (~n - 1)) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && IGtQ(~n, 0)))
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (2 * ~a * ~c * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m * (~c + ~d * csc(~e + ~f * ~x)) ^ (~n - 1)) / (~b * ~f * (2 * ~m + 1)) - ((~d * (2 * ~n - 1)) / (~b * (2 * ~m + 1))) * Antiderivative(csc(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * csc(~e + ~f * ~x)) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IGtQ(~n, 0) && (LtQ(~m, -1 / 2) && IntegerQ(2 * ~m)))))
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (-(~a) * ~c) ^ ~m * Antiderivative(ExpandTrig(csc(~e + ~f * ~x) * cot(~e + ~f * ~x) ^ (2 * ~m), (~c + ~d * csc(~e + ~f * ~x)) ^ (~n - ~m), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IntegersQ(~m, ~n) && (GeQ(~n - ~m, 0) && GtQ(~m * ~n, 0)))))
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (((-(~a) * ~c) ^ (~m + 1 / 2) * cot(~e + ~f * ~x)) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x)))) * Antiderivative(csc(~e + ~f * ~x) * cot(~e + ~f * ~x) ^ (2 * ~m), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && IntegerQ(~m + 1 / 2)))
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (~b * cot(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n) / (~a * ~f * (2 * ~m + 1)) + ((~m + ~n + 1) / (~a * (2 * ~m + 1))) * Antiderivative(csc(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (ILtQ(~m, 0) && ILtQ(~n - 1 / 2, 0) || ILtQ(~m - 1 / 2, 0) && (ILtQ(~n - 1 / 2, 0) && LtQ(~m, ~n)))))
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((~a * ~c * cot(~e + ~f * ~x)) / (~f * sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x)))) * Subst(Antiderivative((~a + ~b * ~x) ^ (~m - 1 / 2) * (~c + ~d * ~x) ^ (~n - 1 / 2), ~x), ~x, csc(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((~(g') * csc(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m') * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~(n'), ~x) => (-(~a) * ~c) ^ ~m * Antiderivative(ExpandTrig((~g * csc(~e + ~f * ~x)) ^ ~p * cot(~e + ~f * ~x) ^ (2 * ~m), (~c + ~d * csc(~e + ~f * ~x)) ^ (~n - ~m), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n, ~p], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IntegersQ(~m, ~n) && (GeQ(~n - ~m, 0) && GtQ(~m * ~n, 0)))))
    @apply_utils Antiderivative((~(g') * csc(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (((-(~a) * ~c) ^ (~m + 1 / 2) * cot(~e + ~f * ~x)) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x)))) * Antiderivative((~g * csc(~e + ~f * ~x)) ^ ~p * cot(~e + ~f * ~x) ^ (2 * ~m), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~p], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && IntegerQ(~m + 1 / 2)))
    @apply_utils Antiderivative((~(g') * csc(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((~a * ~c * ~g * cot(~e + ~f * ~x)) / (~f * sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x)))) * Subst(Antiderivative((~g * ~x) ^ (~p - 1) * (~a + ~b * ~x) ^ (~m - 1 / 2) * (~c + ~d * ~x) ^ (~n - 1 / 2), ~x), ~x, csc(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((sqrt(~(g') * csc(~(e') + ~(f') * ~x)) * sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x))) / (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => ((-2 * ~b * ~g) / ~f) * Subst(Antiderivative(1 / ((~b * ~c + ~a * ~d) - ~c * ~g * (~x) ^ 2), ~x), ~x, (~b * cot(~e + ~f * ~x)) / (sqrt(~g * csc(~e + ~f * ~x)) * sqrt(~a + ~b * csc(~e + ~f * ~x)))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((sqrt(~(g') * csc(~(e') + ~(f') * ~x)) * sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x))) / (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => (~a / ~c) * Antiderivative(sqrt(~g * csc(~e + ~f * ~x)) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) + ((~b * ~c - ~a * ~d) / (~c * ~g)) * Antiderivative((~g * csc(~e + ~f * ~x)) ^ (3 / 2) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * (~c + ~d * csc(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && NeQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((csc(~(e') + ~(f') * ~x) * sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x))) / (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => ((-2 * ~b) / ~f) * Subst(Antiderivative(1 / (~b * ~c + ~a * ~d + ~d * (~x) ^ 2), ~x), ~x, (~b * cot(~e + ~f * ~x)) / sqrt(~a + ~b * csc(~e + ~f * ~x))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((csc(~(e') + ~(f') * ~x) * sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x))) / (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => ((-(sqrt(~a + ~b * csc(~e + ~f * ~x))) * sqrt(~c / (~c + ~d * csc(~e + ~f * ~x)))) / (~d * ~f * sqrt((~c * ~d * (~a + ~b * csc(~e + ~f * ~x))) / ((~b * ~c + ~a * ~d) * (~c + ~d * csc(~e + ~f * ~x)))))) * Elliptic.E(asin((~c * cot(~e + ~f * ~x)) / (~c + ~d * csc(~e + ~f * ~x))), -((~b * ~c - ~a * ~d)) / (~b * ~c + ~a * ~d)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && EqQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative((csc(~(e') + ~(f') * ~x) * sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x))) / (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => (~b / ~d) * Antiderivative(csc(~e + ~f * ~x) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) - ((~b * ~c - ~a * ~d) / ~d) * Antiderivative(csc(~e + ~f * ~x) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * (~c + ~d * csc(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(((~(g') * csc(~(e') + ~(f') * ~x)) ^ (3 / 2) * sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x))) / (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => (~g / ~d) * Antiderivative(sqrt(~g * csc(~e + ~f * ~x)) * sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) - ((~c * ~g) / ~d) * Antiderivative((sqrt(~g * csc(~e + ~f * ~x)) * sqrt(~a + ~b * csc(~e + ~f * ~x))) / (~c + ~d * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative(((~(g') * csc(~(e') + ~(f') * ~x)) ^ (3 / 2) * sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x))) / (~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => (~b / ~d) * Antiderivative((~g * csc(~e + ~f * ~x)) ^ (3 / 2) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) - ((~b * ~c - ~a * ~d) / ~d) * Antiderivative((~g * csc(~e + ~f * ~x)) ^ (3 / 2) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * (~c + ~d * csc(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && NeQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) / (sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~c + ~(d') * csc(~(e') + ~(f') * ~x))), ~x) => (~b / (~b * ~c - ~a * ~d)) * Antiderivative(csc(~e + ~f * ~x) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) - (~d / (~b * ~c - ~a * ~d)) * Antiderivative((csc(~e + ~f * ~x) * sqrt(~a + ~b * csc(~e + ~f * ~x))) / (~c + ~d * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) || EqQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) / (sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~c + ~(d') * csc(~(e') + ~(f') * ~x))), ~x) => ((-2 * cot(~e + ~f * ~x)) / (~f * (~c + ~d) * sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(-(cot(~e + ~f * ~x) ^ 2)))) * sqrt((~a + ~b * csc(~e + ~f * ~x)) / (~a + ~b)) * Elliptic.Pi((2 * ~d) / (~c + ~d), asin(sqrt(1 - csc(~e + ~f * ~x)) / sqrt(2)), (2 * ~b) / (~a + ~b)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative((~(g') * csc(~(e') + ~(f') * ~x)) ^ (3 / 2) / (sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~c + ~(d') * csc(~(e') + ~(f') * ~x))), ~x) => ((-(~a) * ~g) / (~b * ~c - ~a * ~d)) * Antiderivative(sqrt(~g * csc(~e + ~f * ~x)) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) + ((~c * ~g) / (~b * ~c - ~a * ~d)) * Antiderivative((sqrt(~g * csc(~e + ~f * ~x)) * sqrt(~a + ~b * csc(~e + ~f * ~x))) / (~c + ~d * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((~(g') * csc(~(e') + ~(f') * ~x)) ^ (3 / 2) / (sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~c + ~(d') * csc(~(e') + ~(f') * ~x))), ~x) => ((~g * sqrt(~g * csc(~e + ~f * ~x)) * sqrt(~b + ~a * sin(~e + ~f * ~x))) / sqrt(~a + ~b * csc(~e + ~f * ~x))) * Antiderivative(1 / (sqrt(~b + ~a * sin(~e + ~f * ~x)) * (~d + ~c * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && NeQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) ^ 2 / (sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~c + ~(d') * csc(~(e') + ~(f') * ~x))), ~x) => (-(~a) / (~b * ~c - ~a * ~d)) * Antiderivative(csc(~e + ~f * ~x) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) + (~c / (~b * ~c - ~a * ~d)) * Antiderivative((csc(~e + ~f * ~x) * sqrt(~a + ~b * csc(~e + ~f * ~x))) / (~c + ~d * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) || EqQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) ^ 2 / (sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~c + ~(d') * csc(~(e') + ~(f') * ~x))), ~x) => (1 / ~d) * Antiderivative(csc(~e + ~f * ~x) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) - (~c / ~d) * Antiderivative(csc(~e + ~f * ~x) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * (~c + ~d * csc(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative((~(g') * csc(~(e') + ~(f') * ~x)) ^ (5 / 2) / (sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~c + ~(d') * csc(~(e') + ~(f') * ~x))), ~x) => ((-((~c) ^ 2) * (~g) ^ 2) / (~d * (~b * ~c - ~a * ~d))) * Antiderivative((sqrt(~g * csc(~e + ~f * ~x)) * sqrt(~a + ~b * csc(~e + ~f * ~x))) / (~c + ~d * csc(~e + ~f * ~x)), ~x) + ((~g) ^ 2 / (~d * (~b * ~c - ~a * ~d))) * Antiderivative((sqrt(~g * csc(~e + ~f * ~x)) * (~a * ~c + (~b * ~c - ~a * ~d) * csc(~e + ~f * ~x))) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((~(g') * csc(~(e') + ~(f') * ~x)) ^ (5 / 2) / (sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~c + ~(d') * csc(~(e') + ~(f') * ~x))), ~x) => (~g / ~d) * Antiderivative((~g * csc(~e + ~f * ~x)) ^ (3 / 2) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) - ((~c * ~g) / ~d) * Antiderivative((~g * csc(~e + ~f * ~x)) ^ (3 / 2) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * (~c + ~d * csc(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && NeQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((csc(~(e') + ~(f') * ~x) * sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x))) / sqrt(~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => ((-2 * ~b) / ~f) * Subst(Antiderivative(1 / (1 - ~b * ~d * (~x) ^ 2), ~x), ~x, cot(~e + ~f * ~x) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x)))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative((csc(~(e') + ~(f') * ~x) * sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x))) / sqrt(~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => (-((~b * ~c - ~a * ~d)) / ~d) * Antiderivative(csc(~e + ~f * ~x) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x))), ~x) + (~b / ~d) * Antiderivative((csc(~e + ~f * ~x) * sqrt(~c + ~d * csc(~e + ~f * ~x))) / sqrt(~a + ~b * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && EqQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative((csc(~(e') + ~(f') * ~x) * sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x))) / sqrt(~c + ~(d') * csc(~(e') + ~(f') * ~x)), ~x) => ((-2 * (~a + ~b * csc(~e + ~f * ~x))) / (~d * ~f * sqrt((~a + ~b) / (~c + ~d)) * cot(~e + ~f * ~x))) * sqrt((-((~b * ~c - ~a * ~d)) * (1 - csc(~e + ~f * ~x))) / ((~c + ~d) * (~a + ~b * csc(~e + ~f * ~x)))) * sqrt(((~b * ~c - ~a * ~d) * (1 + csc(~e + ~f * ~x))) / ((~c - ~d) * (~a + ~b * csc(~e + ~f * ~x)))) * Elliptic.Pi((~b * (~c + ~d)) / (~d * (~a + ~b)), asin((sqrt((~a + ~b) / (~c + ~d)) * sqrt(~c + ~d * csc(~e + ~f * ~x))) / sqrt(~a + ~b * csc(~e + ~f * ~x))), ((~a - ~b) * (~c + ~d)) / ((~a + ~b) * (~c - ~d))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) / (sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * sqrt(~c + ~(d') * csc(~(e') + ~(f') * ~x))), ~x) => ((-2 * ~a) / (~b * ~f)) * Subst(Antiderivative(1 / (2 + (~a * ~c - ~b * ~d) * (~x) ^ 2), ~x), ~x, cot(~e + ~f * ~x) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x)))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) / (sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * sqrt(~c + ~(d') * csc(~(e') + ~(f') * ~x))), ~x) => ((-2 * (~c + ~d * csc(~e + ~f * ~x))) / (~f * (~b * ~c - ~a * ~d) * Rt((~c + ~d) / (~a + ~b), 2) * cot(~e + ~f * ~x))) * sqrt(((~b * ~c - ~a * ~d) * (1 - csc(~e + ~f * ~x))) / ((~a + ~b) * (~c + ~d * csc(~e + ~f * ~x)))) * sqrt((-((~b * ~c - ~a * ~d)) * (1 + csc(~e + ~f * ~x))) / ((~a - ~b) * (~c + ~d * csc(~e + ~f * ~x)))) * Elliptic.F(asin(Rt((~c + ~d) / (~a + ~b), 2) * (sqrt(~a + ~b * csc(~e + ~f * ~x)) / sqrt(~c + ~d * csc(~e + ~f * ~x)))), ((~a + ~b) * (~c - ~d)) / ((~a - ~b) * (~c + ~d))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) ^ 2 / (sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * sqrt(~c + ~(d') * csc(~(e') + ~(f') * ~x))), ~x) => (-(~a) / ~b) * Antiderivative(csc(~e + ~f * ~x) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x))), ~x) + (1 / ~b) * Antiderivative((csc(~e + ~f * ~x) * sqrt(~a + ~b * csc(~e + ~f * ~x))) / sqrt(~c + ~d * csc(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && NeQ(~b * ~c - ~a * ~d, 0)
    @apply_utils Antiderivative((csc(~(e') + ~(f') * ~x) * sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x))) / (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ (3 / 2), ~x) => ((~a - ~b) / (~c - ~d)) * Antiderivative(csc(~e + ~f * ~x) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x))), ~x) + ((~b * ~c - ~a * ~d) / (~c - ~d)) * Antiderivative((csc(~e + ~f * ~x) * (1 + csc(~e + ~f * ~x))) / (sqrt(~a + ~b * csc(~e + ~f * ~x)) * (~c + ~d * csc(~e + ~f * ~x)) ^ (3 / 2)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ((~c) ^ 2 - (~d) ^ 2, 0)))
    @apply_utils Antiderivative((~(g') * csc(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (((~a) ^ 2 * ~g * cot(~e + ~f * ~x)) / (~f * sqrt(~a + ~b * csc(~e + ~f * ~x)) * sqrt(~a - ~b * csc(~e + ~f * ~x)))) * Subst(Antiderivative(((~g * ~x) ^ (~p - 1) * (~a + ~b * ~x) ^ (~m - 1 / 2) * (~c + ~d * ~x) ^ ~n) / sqrt(~a - ~b * ~x), ~x), ~x, csc(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && (EqQ(~p, 1) || IntegerQ(~m - 1 / 2)))))
    @apply_utils Antiderivative((~(g') * csc(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (1 / (~g) ^ (~m + ~n)) * Antiderivative((~g * csc(~e + ~f * ~x)) ^ (~m + ~n + ~p) * (~b + ~a * sin(~e + ~f * ~x)) ^ ~m * (~d + ~c * sin(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~p], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (IntegerQ(~m) && IntegerQ(~n)))
    @apply_utils Antiderivative((~(g') * csc(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (((~g * csc(~e + ~f * ~x)) ^ (~m + ~p) * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n) / ((~g) ^ ~m * (~d + ~c * sin(~e + ~f * ~x)) ^ ~n)) * Antiderivative((~b + ~a * sin(~e + ~f * ~x)) ^ ~m * (~d + ~c * sin(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n, ~p], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ(~m + ~n + ~p, 0) && IntegerQ(~m)))
    @apply_utils Antiderivative((~(g') * csc(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => (((~g * csc(~e + ~f * ~x)) ^ ~p * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n) / ((~b + ~a * sin(~e + ~f * ~x)) ^ ~m * (~d + ~c * sin(~e + ~f * ~x)) ^ ~n)) * Antiderivative((~b + ~a * sin(~e + ~f * ~x)) ^ ~m * (~d + ~c * sin(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ(~m + ~n + ~p, 0) && Not(IntegerQ(~m))))
    @apply_utils Antiderivative(csc(~(e') + ~(f') * ~x) ^ ~(p') * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => ((sqrt(~d + ~c * sin(~e + ~f * ~x)) * sqrt(~a + ~b * csc(~e + ~f * ~x))) / (sqrt(~b + ~a * sin(~e + ~f * ~x)) * sqrt(~c + ~d * csc(~e + ~f * ~x)))) * Antiderivative(((~b + ~a * sin(~e + ~f * ~x)) ^ ~m * (~d + ~c * sin(~e + ~f * ~x)) ^ ~n) / sin(~e + ~f * ~x) ^ (~m + ~n + ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (IntegerQ(~m - 1 / 2) && (IntegerQ(~n - 1 / 2) && (IntegerQ(~p) && LeQ(-2, ~m + ~n + ~p, -1)))))
    @apply_utils Antiderivative((~(g') * csc(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => Antiderivative(ExpandTrig((~g * csc(~e + ~f * ~x)) ^ ~p * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (IntegersQ(~m, ~n) || (IntegersQ(~m, ~p) || IntegersQ(~n, ~p))))
    @apply_utils Antiderivative((~(g') * csc(~(e') + ~(f') * ~x)) ^ ~(p') * (~(a') + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * csc(~(e') + ~(f') * ~x)) ^ ~n, ~x) => Unintegrable((~g * csc(~e + ~f * ~x)) ^ ~p * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m * (~c + ~d * csc(~e + ~f * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((sec(~(e') + ~(f') * ~x) * (~A + ~(B') * sec(~(e') + ~(f') * ~x))) / (sqrt(~a + ~(b') * sec(~(e') + ~(f') * ~x)) * (~c + ~(d') * sec(~(e') + ~(f') * ~x)) ^ (3 / 2)), ~x) => ((2 * ~A * (1 + sec(~e + ~f * ~x)) * sqrt(((~b * ~c - ~a * ~d) * (1 - sec(~e + ~f * ~x))) / ((~a + ~b) * (~c + ~d * sec(~e + ~f * ~x))))) / (~f * (~b * ~c - ~a * ~d) * Rt((~c + ~d) / (~a + ~b), 2) * tan(~e + ~f * ~x) * sqrt((-((~b * ~c - ~a * ~d)) * (1 + sec(~e + ~f * ~x))) / ((~a - ~b) * (~c + ~d * sec(~e + ~f * ~x)))))) * Elliptic.E(asin((Rt((~c + ~d) / (~a + ~b), 2) * sqrt(~a + ~b * sec(~e + ~f * ~x))) / sqrt(~c + ~d * sec(~e + ~f * ~x))), ((~a + ~b) * (~c - ~d)) / ((~a - ~b) * (~c + ~d))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && EqQ(~A, ~B))))
    @apply_utils Antiderivative((csc(~(e') + ~(f') * ~x) * (~A + ~(B') * csc(~(e') + ~(f') * ~x))) / (sqrt(~a + ~(b') * csc(~(e') + ~(f') * ~x)) * (~c + ~(d') * csc(~(e') + ~(f') * ~x)) ^ (3 / 2)), ~x) => ((-2 * ~A * (1 + csc(~e + ~f * ~x)) * sqrt(((~b * ~c - ~a * ~d) * (1 - csc(~e + ~f * ~x))) / ((~a + ~b) * (~c + ~d * csc(~e + ~f * ~x))))) / (~f * (~b * ~c - ~a * ~d) * Rt((~c + ~d) / (~a + ~b), 2) * cot(~e + ~f * ~x) * sqrt((-((~b * ~c - ~a * ~d)) * (1 + csc(~e + ~f * ~x))) / ((~a - ~b) * (~c + ~d * csc(~e + ~f * ~x)))))) * Elliptic.E(asin((Rt((~c + ~d) / (~a + ~b), 2) * sqrt(~a + ~b * csc(~e + ~f * ~x))) / sqrt(~c + ~d * csc(~e + ~f * ~x))), ((~a + ~b) * (~c - ~d)) / ((~a - ~b) * (~c + ~d))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ((~c) ^ 2 - (~d) ^ 2, 0) && EqQ(~A, ~B))))
end

max_stars_repo_path: src/rules/4 Trig functions/4.5 Secant/4.5.2.2 (g sec)^p (a+b sec)^m (c+d sec)^n.jl


text:
cd(@__DIR__)
using DiscreteEvents, DataStructures, Distributions, Random, BenchmarkTools
using Plots, Printf
const _bench = [false]

Random.seed!(8710) # set random number seed for reproducibility
num_customers = 10 # total number of customers generated
num_servers = 2 # number of servers
mu = 1.0 / 2 # service rate
lam = 0.9 # arrival rate
arrival_dist = Exponential(1 / lam) # interarrival time distriubtion
service_dist = Exponential(1 / mu); # service time distribution

# Define Customer Behavior
function customer(clk::Clock, server::Channel, id::Int, ta::Float64, ds::Distribution)
    delay!(clk, ta)       # customer arrives
    _bench[end] || now!(clk, ()->@printf("%5.3f: customer %d arrived\n", tau(clk), id))
    s = take!(server)     # customer starts service
    _bench[end] || now!(clk, ()->@printf("%5.3f: customer %d entered service\n", tau(clk), id))
    delay!(clk, rand(ds)) # server is busy
    put!(server, s)       # customer exits service
    _bench[end] || now!(clk, ()->@printf("%5.3f: customer %d exited service\n", tau(clk), id))
end

function initialize!(clock, arrival_dist, service_dist, num_customers, server)
    arrival_time = 0.0
    for i = 1:num_customers # initialize customers
        arrival_time += rand(arrival_dist)
        process!(clock, Prc(i, customer, server, i, arrival_time, service_dist), 1)
    end
end

function run_model(arrival_dist, service_dist, num_customers, num_servers, t)
    clock = Clock() # initialize simulation environment
    server = Channel{Int}(num_servers)  # initialize servers
    for i in 1:num_servers
        put!(server, i)
    end
    initialize!(clock, arrival_dist, service_dist, num_customers, server)
    run!(clock, t)
end

function run_model(num_customers)
    result = @benchmark run_model($arrival_dist, $service_dist, $num_customers, $num_servers, 2000000)
    return mean(result).time*1e-9
end

_bench[end] = true
N = [10,1000,2000,4000]

times = run_model.(N)
plot(N, times, xlabel="Customers", ylabel="Time (seconds)", leg=false, grid=false)
savefig("DiscreteEvents_Queue_Channel.png")