text:
# here we keep code related to the lattice itself
# note that for the time being I will mostly specialize to 2 dimensions, later I will add code
# for 3 and 4 dimensions. This is being done to avoid splatting.

# TODO should probably implement specialized spacetime vectors using StaticArrays
# TODO fix spacetime indexing scheme

Base.getindex(lat::AbstractLattice, idx::Tuple) = idx
Base.getindex(lat::AbstractLattice, idx...) = getindex(lat, idx)

"""
    fieldsites(::Type{T}, lat::AbstractLattice)

Creates an uninitialized array with elements of type T for every poin on the lattice.
"""
function fieldsites(::Type{T}, lat::AbstractLattice{d}) where {T, d}
    Array{T, d}((lat.L for i ∈ 1:(d-1))..., lat.T)
end
export fieldsites

nsites(l::AbstractLattice{d}) where d = l.T*(l.L^(d-1))
export nsites

# this is an implementation that doesn't involve splatting
function Base.sub2ind(M::Tuple, m::Tuple)
    1 + sum((m[n]-1)*prod(M[1:(n-1)]) for n ∈ 1:length(M))
end

size(lat::AbstractLattice{d}) where d = tuple((lat.L for i ∈ 1:(d-1))..., lat.T)

sizenocheck(lat::AbstractLattice{d}, idx::Integer) where d = (idx == d) ? lat.T : lat.L

"""
    eachspacetimeindex(lat::AbstractLattice)

Return an iterator over every spacetime index in the lattice `lat`.  Spacetime indices will be
in the form of `SVector`s.
"""
function eachspacetimeindex(lat::AbstractLattice{d}) where d
    dims = size(lat)
    (SVector{d}(ind2sub(dims , x)) for x ∈ 1:nsites(lat))
end
export eachspacetimeindex

function size(lat::AbstractLattice{d}, idx::Integer) where d
    if idx == d
        return lat.T
    elseif 1 ≤ idx < d
        return lat.D
    else
        ArgumentError("Lattice only has $d spacetime dimensions.")
    end
end

"""
    LatticeToroidal <: AbstractLattice

A type for storing parameters of a simple Euclidean rectilinear lattice with D spatial dimensions
and toroidal topology in all directions.
"""
struct LatticeToroidal{d} <: AbstractLattice{d}
    a::Float64  # lattice spacing in GeV-1
    L::Int  # lattice length in units of a
    T::Int  # lattice duration in units of a
end
export LatticeToroidal

LatticeToroidal(d::Int, a::Float64, L::Integer, T::Integer) = LatticeToroidal{d}(a, L, T)
LatticeToroidal(d::Int, a::Float64, L::Integer) = LatticeToroidal(d, a, L, L)

# helper function for index computers
@inline function _loopedidx(i::Integer, L::Integer)
    # pretty sure this has to be different for negative and positive?
    # Note to self: think in terms of 0-based indexing, otherwise confusing as fuck
    if i < 1
        L + (i % L)
    else
        ((i-1) % L) + 1
    end
end

# these duplications are to avoid splatting
# WARNING right now these don't check lattice dimensions
# for rank-2
function Base.getindex(lat::LatticeToroidal, i::Integer, j::Integer)
    _loopedidx(i, lat.L), _loopedidx(j, lat.T)
end
# for rank-3
function Base.getindex(lat::LatticeToroidal, i::Integer, j::Integer, k::Integer)
    _loopedidx(i, lat.L), _loopedidx(j, lat.L), _loopedidx(k, lat.T)
end
# for rank-4
function Base.getindex(lat::LatticeToroidal, i::Integer, j::Integer, k::Integer, l::Integer)
    _loopedidx(i, lat.L), _loopedidx(j, lat.L), _loopedidx(k, lat.L), _loopedidx(l, lat.L)
end
Base.getindex(lat::LatticeToroidal, idx...) = getindex(lat, idx)
Base.getindex(lat::LatticeToroidal, idx::Tuple) = getindex(lat, idx...)
Base.getindex(lat::LatticeToroidal, idx::SVector) = getindex(lat, Tuple(idx))

# TODO compiler will not know size of these!
spacetimebasisvec(μ::Integer, N::Integer) = SVector{N}([i ≠ μ ? 0 : 1 for i ∈ 1:N])
spacetimebasisvec(μ::Integer, lat::AbstractLattice{d}) where d = spacetimebasisvec(μ, d)
export spacetimebasisvec

# TODO this is probably expensive and shitty
inc(x::SVector{N}, μ::Integer) where N = SVector{N}([i ≠ μ ? x[i] : x[i]+1 for i ∈ 1:N])
Base.dec(x::SVector{N}, μ::Integer) where N = SVector{N}([i ≠ μ ? x[i] : x[i]-1 for i ∈ 1:N])
export inc


text:
using Test
import ParallelStencil
using ParallelStencil.ParallelKernel
import ParallelStencil.ParallelKernel: @reset_parallel_kernel, @is_initialized, SUPPORTED_PACKAGES, PKG_CUDA, PKG_THREADS
import ParallelStencil.ParallelKernel: @require, @prettystring
import ParallelStencil.ParallelKernel: checknoargs, checkargs_sharedMem, Dim3
using ParallelStencil.ParallelKernel.Exceptions
TEST_PACKAGES = SUPPORTED_PACKAGES
@static if PKG_CUDA in TEST_PACKAGES
    import CUDA
    if !CUDA.functional() TEST_PACKAGES = filter!(x->x≠PKG_CUDA, TEST_PACKAGES) end
end

@static for package in TEST_PACKAGES  eval(:(
    @testset "$(basename(@__FILE__)) (package: $(nameof($package)))" begin
        @testset "1. kernel language macros" begin
            @testset "mapping to package" begin
                @require !@is_initialized()
                @init_parallel_kernel($package, Float64)
                @require @is_initialized()
                if $package == $PKG_CUDA
                    @test @prettystring(1, @gridDim()) == "CUDA.gridDim()"
                    @test @prettystring(1, @blockIdx()) == "CUDA.blockIdx()"
                    @test @prettystring(1, @blockDim()) == "CUDA.blockDim()"
                    @test @prettystring(1, @threadIdx()) == "CUDA.threadIdx()"
                    @test @prettystring(1, @sync_threads()) == "CUDA.sync_threads()"
                    @test @prettystring(1, @sharedMem(Float32, (2,3))) == "CUDA.@cuDynamicSharedMem Float32 (2, 3)"
                    @test @prettystring(1, @pk_show()) == "CUDA.@cushow"
                    @test @prettystring(1, @pk_println()) == "CUDA.@cuprintln"
                elseif $package == $PKG_THREADS
                    @test @prettystring(1, @gridDim()) == "ParallelStencil.ParallelKernel.@gridDim_cpu"
                    @test @prettystring(1, @blockIdx()) == "ParallelStencil.ParallelKernel.@blockIdx_cpu"
                    @test @prettystring(1, @blockDim()) == "ParallelStencil.ParallelKernel.@blockDim_cpu"
                    @test @prettystring(1, @threadIdx()) == "ParallelStencil.ParallelKernel.@threadIdx_cpu"
                    @test @prettystring(1, @sync_threads()) == "ParallelStencil.ParallelKernel.@sync_threads_cpu"
                    @test @prettystring(1, @sharedMem(Float32, (2,3))) == "ParallelStencil.ParallelKernel.@sharedMem_cpu Float32 (2, 3)"
                    @test @prettystring(1, @pk_show()) == "Base.@show"
                    @test @prettystring(1, @pk_println()) == "Base.@println"
                end;
            end;
            @testset "@gridDim, @blockIdx, @blockDim, @threadIdx (1D)" begin
                @static if $package == $PKG_THREADS
                    A  = @zeros(4)
                    @parallel_indices (ix) function test_macros!(A)
                        @test @gridDim() == Dim3(2, 1, 1)
                        @test @blockIdx() == Dim3(ix-2, 1, 1)
                        @test @blockDim() == Dim3(1, 1, 1)
                        @test @threadIdx() == Dim3(1, 1, 1)
                        return
                    end
                    @parallel (3:4) test_macros!(A);
                    @test true # @gridDim test succeeded if this line is reached (the above tests within the kernel are not captured if they succeed, only if they fail...; alternative test implementations would be of course possible, but would be more complex).
                    @test true # @blockIdx ...
                    @test true # @blockDim ...
                    @test true # @threadIdx ...
                end
            end;
            @testset "@gridDim, @blockIdx, @blockDim, @threadIdx (2D)" begin
                @static if $package == $PKG_THREADS
                    A  = @zeros(4, 5)
                    @parallel_indices (ix,iy) function test_macros!(A)
                        @test @gridDim() == Dim3(2, 3, 1)
                        @test @blockIdx() == Dim3(ix-2, iy-1, 1)
                        @test @blockDim() == Dim3(1, 1, 1)
                        @test @threadIdx() == Dim3(1, 1, 1)
                        return
                    end
                    @parallel (3:4, 2:4) test_macros!(A);
                    @test true # @gridDim test succeeded if this line is reached (the above tests within the kernel are not captured if they succeed, only if they fail...; alternative test implementations would be of course possible, but would be more complex).
                    @test true # @blockIdx ...
                    @test true # @blockDim ...
                    @test true # @threadIdx ...
                end
            end;
            @testset "@gridDim, @blockIdx, @blockDim, @threadIdx (3D)" begin
                @static if $package == $PKG_THREADS
                    A  = @zeros(4, 5, 6)
                    @parallel_indices (ix,iy,iz) function test_macros!(A)
                        @test @gridDim() == Dim3(2, 3, 6)
                        @test @blockIdx() == Dim3(ix-2, iy-1, iz)
                        @test @blockDim() == Dim3(1, 1, 1)
                        @test @threadIdx() == Dim3(1, 1, 1)
                        return
                    end
                    @parallel (3:4, 2:4, 1:6) test_macros!(A);
                    @test true # @gridDim test succeeded if this line is reached (the above tests within the kernel are not captured if they succeed, only if they fail...; alternative test implementations would be of course possible, but would be more complex).
                    @test true # @blockIdx ...
                    @test true # @blockDim ...
                    @test true # @threadIdx ...
                end
            end;
            @testset "sync_threads" begin
                @static if $package == $PKG_THREADS
                    @test @prettystring(ParallelStencil.ParallelKernel.@sync_threads_cpu()) == "begin\nend"
                end;
            end;
            @testset "shared memory (allocation)" begin
                @static if $package == $PKG_THREADS
                    @test typeof(@sharedMem(Float32,(2,3))) == typeof(ParallelStencil.ParallelKernel.MArray{Tuple{2,3},   Float32, length((2,3)),   prod((2,3))}(undef))
                    @test typeof(@sharedMem(Bool,(2,3,4)))  == typeof(ParallelStencil.ParallelKernel.MArray{Tuple{2,3,4}, Bool,    length((2,3,4)), prod((2,3,4))}(undef))
                end;
            end;
            @testset "@sharedMem (1D)" begin
                @static if $package == $PKG_THREADS
                    A  = @rand(4)
                    B  = @zeros(4)
                    @parallel_indices (ix) function memcopy!(B, A)
                        tx  = @threadIdx().x
                        A_l = @sharedMem(eltype(A), (@blockDim().x))
                        A_l[tx] = A[ix]
                        @sync_threads()
                        B[ix] = A_l[tx]
                        return
                    end
                    @parallel memcopy!(B, A);
                    @test B == A
                end
            end;
            @testset "@sharedMem (2D)" begin
                @static if $package == $PKG_THREADS
                    A  = @rand(4,5)
                    B  = @zeros(4,5)
                    @parallel_indices (ix,iy) function memcopy!(B, A)
                        tx  = @threadIdx().x
                        ty  = @threadIdx().y
                        A_l = @sharedMem(eltype(A), (@blockDim().x, @blockDim().y))
                        A_l[tx,ty] = A[ix,iy]
                        @sync_threads()
                        B[ix,iy] = A_l[tx,ty]
                        return
                    end
                    @parallel memcopy!(B, A);
                    @test B == A
                end
            end;
            @testset "@sharedMem (3D)" begin
                @static if $package == $PKG_THREADS
                    A  = @rand(4,5,6)
                    B  = @zeros(4,5,6)
                    @parallel_indices (ix,iy,iz) function memcopy!(B, A)
                        tx  = @threadIdx().x
                        ty  = @threadIdx().y
                        tz  = @threadIdx().z
                        A_l = @sharedMem(eltype(A), (@blockDim().x, @blockDim().y, @blockDim().z))
                        A_l[tx,ty,tz] = A[ix,iy,iz]
                        @sync_threads()
                        B[ix,iy,iz] = A_l[tx,ty,tz]
                        return
                    end
                    @parallel memcopy!(B, A);
                    @test B == A
                end
            end;
            @reset_parallel_kernel()
        end;
        @testset "2. Exceptions" begin
            @init_parallel_kernel($package, Float64)
            @require @is_initialized
            @testset "no arguments" begin
                @test_throws ArgumentError checknoargs(:(something));                                      # Error: length(args) != 0
            end;
            @testset "arguments @sharedMem" begin
                @test_throws ArgumentError checkargs_sharedMem();                                          # Error: isempty(args)
                @test_throws ArgumentError checkargs_sharedMem(:(something));                              # Error: length(args) != 2
                @test_throws ArgumentError checkargs_sharedMem(:(something), :(something), :(something));  # Error: length(args) != 2
            end;
            @reset_parallel_kernel()
        end;
    end;
)) end == nothing || true;

ext: jl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
include("../common.jl")

name = "nq"
upstream_version = v"2.5.5" # when you increment this, reset offset to v"0.0.0"
offset = v"0.0.1" # increment this when rebuilding with unchanged upstream_version
version = offset_version(upstream_version, offset)

# nq just produces a binary and does not need GAP for this at all.
# So let's *not* use common.jl

# Collection of sources required to build libsingular-julia
sources = [
    ArchiveSource("https://github.com/gap-packages/$(name)/releases/download/v$(upstream_version)/$(name)-$(upstream_version).tar.bz2",
                  "b7a8fa12c02db78dc05ce047b451e7855345608bbfb0d11eb105056527694e42"),
    DirectorySource("../bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd nq*

atomic_patch -p1 ${WORKSPACE}/srcdir/patches/nq.patch
autoreconf -fiv
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --with-gmp=${prefix}
make -j${nproc}

# copy just the nq executable
mkdir -p ${prefix}/bin/
cp bin/nq ${prefix}/bin/

install_license LICENSE
"""

name = gap_pkg_name(name)

platforms = supported_platforms()
filter!(p -> nbits(p) == 64, platforms) # we only care about 64bit builds
filter!(!Sys.iswindows, platforms)      # Windows is not supported

dependencies = [
    Dependency("GMP_jll", v"6.2.0"),
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("nq", :nq),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat="1.6", preferred_gcc_version=v"7")


text:
export SpeakerListenerEnv

# TODO: review the implementation of this environment

mutable struct SpeakerListenerEnv{T<:Vector{Float64}} <: AbstractEnv
    target::T
    content::T
    player_vel::T
    player_pos::T
    landmarks_pos::Vector{T}
    landmarks_num::Int
    ϵ::Any
    damping::Any
    max_accel::Any
    space_dim::Int
    init_step::Int
    play_step::Int
    max_steps::Int
    continuous::Bool
end

"""
    SpeakerListenerEnv(;kwargs...)

`SpeakerListenerEnv` is a simple cooperative environment of two agents, a `Speaker` and a `Listener`, who are placed in an environment
with `N` landmarks. At each episode, the `Listener` must navigate to a particular landmark(`env.target`) and obtains **reward** based
on its distance to the target. However, while the `Listener` can observe the relative position of the landmarks, it doesn't know which
is the target landmark. Conversely, the `Speaker` can observe the target's landmark, and it can produce a communication output(`env.content`)
at each time step which is observed by the `Listener`.

For more concrete description, you can refer to:
* [Multi-Agent Actor-Critic for Mixed Cooperative-Competitive Environments](https://arxiv.org/abs/1706.02275)
* [multiagent-particle-envs](https://github.com/openai/multiagent-particle-envs)

# Keyword arguments
- `N::Int = 3`, the number of landmarks in the environment.
- `stop = 0.01`, when the distance between the `Listener` and the target is smaller than the `stop`, the game will be terminated.
- `damping = 0.25`, for simulation of the physical space, `Listener`'s action will meet the damping in each step.
- `max_accel = 0.02`, the maximum acceleration of the `Listener` in each step.
- `space_dim::Int = 2`, the dimension of the environment's space.
- `max_steps::Int = 25`, the maximum playing steps in one episode.
- `continuous::Bool = true`, set to `false` if you want the action_space of the players to be discrete. Otherwise, the action_space will be continuous.
"""
function SpeakerListenerEnv(;
    N::Int = 3,
    stop = 0.01,
    damping = 0.25,
    max_accel = 0.5,
    space_dim::Int = 2,
    max_steps::Int = 50,
    continuous::Bool = true,
)
    SpeakerListenerEnv(
        zeros(N),
        zeros(N),
        zeros(space_dim),
        zeros(space_dim),
        [zeros(space_dim) for _ in Base.OneTo(N)],
        N,
        stop,
        damping,
        max_accel,
        space_dim,
        0,
        0,
        max_steps,
        continuous,
    )
end

function RLBase.reset!(env::SpeakerListenerEnv)
    env.init_step = 0
    env.play_step = 0
    env.target = zeros(env.landmarks_num)
    env.content = zeros(env.landmarks_num)
    env.player_pos = zeros(env.space_dim)
    env.landmarks_pos = [zeros(env.space_dim) for _ in Base.OneTo(env.landmarks_num)]
end

RLBase.is_terminated(env::SpeakerListenerEnv) =
    (reward(env) > -env.ϵ) || (env.play_step > env.max_steps)
RLBase.players(::SpeakerListenerEnv) = (:Speaker, :Listener, CHANCE_PLAYER)

RLBase.state(env::SpeakerListenerEnv, ::Observation{Any}, players::Tuple) =
    Dict(p => state(env, p) for p in players)

RLBase.state(env::SpeakerListenerEnv, ::Observation{Any}, player::Symbol) =
# for speaker, it can observe the target and help listener to arrive it.
    if player == :Speaker
        env.target
        # for listener, it can observe current velocity, relative positions of landmarks, and speaker's conveyed information.
    elseif player == :Listener
        vcat(
            env.player_vel...,
            (
                vcat((landmark_pos .- env.player_pos)...) for
                landmark_pos in env.landmarks_pos
            )...,
            env.content...,
        )
    else
        @error "No player $player."
    end

RLBase.state(env::SpeakerListenerEnv, ::Observation{Any}, ::ChancePlayer) =
    vcat(env.landmarks_pos, [env.player_pos])

RLBase.state_space(env::SpeakerListenerEnv, ::Observation{Any}, players::Tuple) =
    Dict(player => state_space(env, player) for player in players)

RLBase.state_space(env::SpeakerListenerEnv, ::Observation{Any}, player::Symbol) =
    if player == :Speaker
        # env.target
        Space(SVector(ntuple(_ -> (0.0, 1.0), env.landmarks_num)))
    elseif player == :Listener
        Space(
            Float64,
            length(env.player_vel) +
            length(env.landmarks_pos) * length(env.player_pos) +
            length(env.content),
        )
    else
        @error "No player $player."
    end

RLBase.state_space(env::SpeakerListenerEnv, ::Observation{Any}, ::ChancePlayer) = [
    # landmarks' positions
    (
        Space((SVector(ntuple(_ -> -Inf .. Inf, env.space_dim)))) for
        _ in Base.OneTo(env.landmarks_num)
    )...,
    # player's position, no bounds.
    Space(SVector(ntuple(_ -> -Inf .. Inf, env.space_dim))),
]

RLBase.action_space(env::SpeakerListenerEnv, players::Tuple) =
    Space(Dict(p => action_space(env, p) for p in players))

RLBase.action_space(env::SpeakerListenerEnv, player::Symbol) =
    if player == :Speaker
        env.continuous ?
        Space([ClosedInterval(0, 1) for _ in Base.OneTo(env.landmarks_num)]) :
        Space([ZeroTo(1) for _ in Base.OneTo(env.landmarks_num)])
    elseif player == :Listener
        # there has two directions in each dimension.
        env.continuous ?
        Space([ClosedInterval(0, 1) for _ in Base.OneTo(2 * env.space_dim)]) :
        Space([ZeroTo(1) for _ in Base.OneTo(2 * env.space_dim)])
    else
        @error "No player $player."
    end

function RLBase.action_space(env::SpeakerListenerEnv, ::ChancePlayer)
    if env.init_step < env.landmarks_num + 1
        Space(SVector(ntuple(_ -> -1.0 .. 1.0, env.space_dim)))
    else
        Space(OneTo(env.landmarks_num))
    end
end

function (env::SpeakerListenerEnv)(action, ::ChancePlayer)
    env.init_step += 1
    if env.init_step <= env.landmarks_num
        env.landmarks_pos[env.init_step] = action
    elseif env.init_step == env.landmarks_num + 1
        env.player_pos = action
    else
        @assert action in Base.OneTo(env.landmarks_num) "The target should be assigned to one of the landmarks."
        env.target[action] = 1.0
    end
end

function (env::SpeakerListenerEnv)(actions::Dict, players::Tuple)
    @assert length(actions) == length(players)
    for p in players
        env(actions[p], p)
    end
    env.play_step += 1
end

function (env::SpeakerListenerEnv)(action::Vector, player::Symbol)
    if player == :Speaker
        # update conveyed information.
        env.content = round.(action)
    elseif player == :Listener
        # update velocity, here env.damping is for simulation physical rule.
        action = round.(action)
        acceleration = [action[2*i] - action[2*i-1] for i in Base.OneTo(env.space_dim)]
        env.player_vel .*= (1 - env.damping)
        env.player_vel .+= (acceleration * env.max_accel)
        # update position
        env.player_pos .+= env.player_vel * 0.1 # velocity * time
    else
        @error "No player $player."
    end
end

RLBase.reward(::SpeakerListenerEnv, ::ChancePlayer) = -Inf

function RLBase.reward(env::SpeakerListenerEnv, p)
    if sum(env.target) == 1
        goal = findfirst(env.target .== 1.0)
        -sum((env.landmarks_pos[goal] .- env.player_pos) .^ 2)
    else
        -Inf
    end
end

RLBase.current_player(env::SpeakerListenerEnv) =
    if env.init_step < env.landmarks_num + 2
        CHANCE_PLAYER
    else
        (:Speaker, :Listener)
    end

RLBase.NumAgentStyle(::SpeakerListenerEnv) = MultiAgent(2)
RLBase.DynamicStyle(::SpeakerListenerEnv) = SIMULTANEOUS
RLBase.ActionStyle(::SpeakerListenerEnv) = MINIMAL_ACTION_SET
RLBase.InformationStyle(::SpeakerListenerEnv) = IMPERFECT_INFORMATION
RLBase.ChanceStyle(::SpeakerListenerEnv) = EXPLICIT_STOCHASTIC


text:
using Asap
# Nodes
a = Node([0., 0., 0.], :frame, :fixed)
b = Node([0., 1., 0.], :frame, :free)
c = Node([1., 1., 0.] / sqrt(2), :frame, :free)
nodes = [a, b]
nodes2 = [a, c]
# elements
begin
    E = 200e6
    A = 9290/1e6
    G = 77e6
    I1 = 113e-6
    I2 = 38.8e-6
    J = 575e-3
end

# Iz is strong axis, Iy is weak axis
e = Element(nodes, [1, 2], E, A, G, I1, I2, J)
elements = [e]

# load
l = Load(nodes, b.position, [0., 0., -10., 0., 0., 0.])
loads = [l]

#structure
s = Structure(nodes, elements, loads)
analyze(s)

u1 = s.U
r1 = s.reactions

# Iz is weak axis, Iy is strong axis
e2 = Element(nodes, [1, 2], E, A, G, I2, I1, J)
e2.Ψ = 0.
elements2 = [e2]

s2 = Structure(nodes, elements2, loads)
analyze(s2)

u2 = s2.U
r2 = s2.reactions

max_stars_repo_name: keithjlee/Asap

max_issues_repo_name: keithjlee/Asap

max_forks_repo_name: keithjlee/Asap


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

@test integrate((im+tan(x))^-1*sin(x)^4, x) == :(-5//32*(im+tan(x))^-2+-1//32*(im+-1*tan(x))^-2+-1*im*(-8*tan(x)+8im)^-1+-3*im*(16im+16*tan(x))^-1+-1//16*im*x+1//24*im*(im+tan(x))^-3)
@test integrate((im+tan(x))^-1*sin(x)^3, x) == :(1//5*sin(x)^5+-1//5*im*cos(x)^5+1//3*im*cos(x)^3)
@test integrate((im+tan(x))^-1*sin(x)^2, x) == :(-1//8*(im+tan(x))^-2+-1*im*(-8*tan(x)+8im)^-1+-1*im*(4im+4*tan(x))^-1+-1//8*im*x)
@test integrate((im+tan(x))^-1*sin(x), x) == :(1//3*sin(x)^3+1//3*im*cos(x)^3)
@test integrate((im+tan(x))^-1*csc(x), x) == :(im*arctanh(cos(x))+-1*im*cos(x)+sin(x))
@test integrate((im+tan(x))^-1*csc(x)^2, x) == :(im*x+im*cot(x)+log(cos(x))+log(tan(x)))
@test integrate((im+tan(x))^-1*csc(x)^3, x) == :(-1*csc(x)+-1//2*im*arctanh(cos(x))+(1/2)*im*cot(x)*csc(x))
@test integrate((im+tan(x))^-1*csc(x)^4, x) == :(-1//2*cot(x)^2+1//3*im*cot(x)^3)
@test integrate((im+tan(x))^-1*csc(x)^5, x) == :(-1//3*csc(x)^3+-1//8*im*arctanh(cos(x))+-1//8*im*cot(x)*csc(x)+1//4*im*csc(x)^3*cot(x))
@test integrate((im+tan(x))^-1*csc(x)^6, x) == :(-1//2*cot(x)^2+-1//4*cot(x)^4+1//3*im*cot(x)^3+1//5*im*cot(x)^5)
@test integrate(sin(c+d*x)^5*(a+b*tan(c+d*x)), x) == :(b*d^-1*arctanh(sin(c+d*x))+-1*a*d^-1*cos(c+d*x)+-1*b*d^-1*sin(c+d*x)+-1//3*b*d^-1*sin(c+d*x)^3+-1//5*a*d^-1*cos(c+d*x)^5+-1//5*b*d^-1*sin(c+d*x)^5+2//3*a*d^-1*cos(c+d*x)^3)
@test integrate(sin(c+d*x)^4*(a+b*tan(c+d*x)), x) == :(3//8*a*x+-1*b*d^-1*log(cos(c+d*x))+-1//4*d^-1*sin(c+d*x)^3*(a+b*tan(c+d*x))*cos(c+d*x)+-1//8*d^-1*(3a+4*b*tan(c+d*x))*cos(c+d*x)*sin(c+d*x))
@test integrate(sin(c+d*x)^3*(a+b*tan(c+d*x)), x) == :(b*d^-1*arctanh(sin(c+d*x))+-1*a*d^-1*cos(c+d*x)+-1*b*d^-1*sin(c+d*x)+-1//3*b*d^-1*sin(c+d*x)^3+1//3*a*d^-1*cos(c+d*x)^3)
@test integrate(sin(c+d*x)^2*(a+b*tan(c+d*x)), x) == :((1/2)*a*x+-1*b*d^-1*log(cos(c+d*x))+-1//2*d^-1*(a+b*tan(c+d*x))*cos(c+d*x)*sin(c+d*x))
@test integrate((a+b*tan(c+d*x))*sin(c+d*x), x) == :(b*d^-1*arctanh(sin(c+d*x))+-1*a*d^-1*cos(c+d*x)+-1*b*d^-1*sin(c+d*x))
@test integrate((a+b*tan(c+d*x))*csc(c+d*x), x) == :(b*d^-1*arctanh(sin(c+d*x))+-1*a*d^-1*arctanh(cos(c+d*x)))
@test integrate(csc(c+d*x)^2*(a+b*tan(c+d*x)), x) == :(b*d^-1*log(tan(c+d*x))+-1*a*d^-1*cot(c+d*x))
@test integrate(csc(c+d*x)^3*(a+b*tan(c+d*x)), x) == :(b*d^-1*arctanh(sin(c+d*x))+-1*b*d^-1*csc(c+d*x)+-1//2*a*d^-1*arctanh(cos(c+d*x))+-1//2*a*d^-1*cot(c+d*x)*csc(c+d*x))
@test integrate(csc(c+d*x)^4*(a+b*tan(c+d*x)), x) == :(b*d^-1*log(tan(c+d*x))+-1*a*d^-1*cot(c+d*x)+-1//2*b*d^-1*cot(c+d*x)^2+-1//3*a*d^-1*cot(c+d*x)^3)
@test integrate(csc(c+d*x)^5*(a+b*tan(c+d*x)), x) == :(b*d^-1*arctanh(sin(c+d*x))+-1*b*d^-1*csc(c+d*x)+-3//8*a*d^-1*arctanh(cos(c+d*x))+-1//3*b*d^-1*csc(c+d*x)^3+-3//8*a*d^-1*cot(c+d*x)*csc(c+d*x)+-1//4*a*d^-1*csc(c+d*x)^3*cot(c+d*x))
@test integrate(csc(c+d*x)^6*(a+b*tan(c+d*x)), x) == :(b*d^-1*log(tan(c+d*x))+-1*a*d^-1*cot(c+d*x)+-1*b*d^-1*cot(c+d*x)^2+-2//3*a*d^-1*cot(c+d*x)^3+-1//4*b*d^-1*cot(c+d*x)^4+-1//5*a*d^-1*cot(c+d*x)^5)
@test integrate((a+b*tan(c+d*x))^2*sin(c+d*x)^4, x) == :(x*(-15//8*b^2+3//8*a^2)+b^2*d^-1*tan(c+d*x)+-2*a*b*d^-1*log(cos(c+d*x))+1//4*d^-1*(a+b*tan(c+d*x))^2*cos(c+d*x)^3*sin(c+d*x)+1//8*d^-1*cos(c+d*x)^2*(a+b*tan(c+d*x))*(7b+-5*a*tan(c+d*x)))
@test integrate((a+b*tan(c+d*x))^2*sin(c+d*x)^3, x) == :(b^2*d^-1*sec(c+d*x)+-1*a^2*d^-1*cos(c+d*x)+2*b^2*d^-1*cos(c+d*x)+-1//3*b^2*d^-1*cos(c+d*x)^3+1//3*a^2*d^-1*cos(c+d*x)^3+-2*a*b*d^-1*sin(c+d*x)+2*a*b*d^-1*arctanh(sin(c+d*x))+-2//3*a*b*d^-1*sin(c+d*x)^3)
@test integrate((a+b*tan(c+d*x))^2*sin(c+d*x)^2, x) == :(x*((1/2)*a^2+-3//2*b^2)+3//2*b^2*d^-1*tan(c+d*x)+-2*a*b*d^-1*log(cos(c+d*x))+-1//2*d^-1*(a+b*tan(c+d*x))^2*cos(c+d*x)*sin(c+d*x))
@test integrate((a+b*tan(c+d*x))^2*sin(c+d*x), x) == :(b^2*d^-1*cos(c+d*x)+b^2*d^-1*sec(c+d*x)+-1*a^2*d^-1*cos(c+d*x)+-2*a*b*d^-1*sin(c+d*x)+2*a*b*d^-1*arctanh(sin(c+d*x)))
@test integrate((a+b*tan(c+d*x))^2*csc(c+d*x), x) == :(b^2*d^-1*sec(c+d*x)+-1*a^2*d^-1*arctanh(cos(c+d*x))+2*a*b*d^-1*arctanh(sin(c+d*x)))
@test integrate((a+b*tan(c+d*x))^2*csc(c+d*x)^2, x) == :(b^2*d^-1*tan(c+d*x)+-1*a^2*d^-1*cot(c+d*x)+2*a*b*d^-1*log(tan(c+d*x)))
@test integrate((a+b*tan(c+d*x))^2*csc(c+d*x)^3, x) == :(b^2*d^-1*sec(c+d*x)+-1*b^2*d^-1*arctanh(cos(c+d*x))+-1//2*a^2*d^-1*arctanh(cos(c+d*x))+-2*a*b*d^-1*csc(c+d*x)+2*a*b*d^-1*arctanh(sin(c+d*x))+-1//2*a^2*d^-1*cot(c+d*x)*csc(c+d*x))
@test integrate((a+b*tan(c+d*x))^2*csc(c+d*x)^4, x) == :(b^2*d^-1*tan(c+d*x)+-1*d^-1*(a^2+b^2)*cot(c+d*x)+-1//3*a^2*d^-1*cot(c+d*x)^3+-1*a*b*d^-1*cot(c+d*x)^2+2*a*b*d^-1*log(tan(c+d*x)))
@test integrate((a+b*tan(c+d*x))^2*csc(c+d*x)^5, x) == :(-3//2*b^2*d^-1*arctanh(cos(c+d*x))+-3//8*a^2*d^-1*arctanh(cos(c+d*x))+3//2*b^2*d^-1*sec(c+d*x)+-2*a*b*d^-1*csc(c+d*x)+2*a*b*d^-1*arctanh(sin(c+d*x))+-3//8*a^2*d^-1*cot(c+d*x)*csc(c+d*x)+-2//3*a*b*d^-1*csc(c+d*x)^3+-1//2*b^2*d^-1*csc(c+d*x)^2*sec(c+d*x)+-1//4*a^2*d^-1*csc(c+d*x)^3*cot(c+d*x))
@test integrate((a+b*tan(c+d*x))^2*csc(c+d*x)^6, x) == :(b^2*d^-1*tan(c+d*x)+-1*d^-1*(a^2+2*b^2)*cot(c+d*x)+-1//3*d^-1*cot(c+d*x)^3*(b^2+2*a^2)+-1//5*a^2*d^-1*cot(c+d*x)^5+-2*a*b*d^-1*cot(c+d*x)^2+2*a*b*d^-1*log(tan(c+d*x))+-1//2*a*b*d^-1*cot(c+d*x)^4)
@test integrate((a+b*tan(c+d*x))^3*sin(c+d*x)^3, x) == :(-1*a^3*d^-1*cos(c+d*x)+-5//2*b^3*d^-1*arctanh(sin(c+d*x))+1//3*a^3*d^-1*cos(c+d*x)^3+5//2*b^3*d^-1*sin(c+d*x)+5//6*b^3*d^-1*sin(c+d*x)^3+(1/2)*b^3*d^-1*sin(c+d*x)^3*tan(c+d*x)^2+-1*a*b^2*d^-1*cos(c+d*x)^3+-1*b*a^2*d^-1*sin(c+d*x)^3+-3*b*a^2*d^-1*sin(c+d*x)+3*a*b^2*d^-1*sec(c+d*x)+3*b*a^2*d^-1*arctanh(sin(c+d*x))+6*a*b^2*d^-1*cos(c+d*x))
@test integrate((a+b*tan(c+d*x))^3*sin(c+d*x)^2, x) == :(b^3*d^-1*tan(c+d*x)^2+(1/2)*a*x*(a^2+-9*b^2)+-1*b*d^-1*(-2*b^2+3*a^2)*log(cos(c+d*x))+-1//2*d^-1*(a+b*tan(c+d*x))^3*cos(c+d*x)*sin(c+d*x)+9//2*a*b^2*d^-1*tan(c+d*x))
@test integrate((a+b*tan(c+d*x))^3*sin(c+d*x), x) == :(-1*a^3*d^-1*cos(c+d*x)+-3//2*b^3*d^-1*arctanh(sin(c+d*x))+3//2*b^3*d^-1*sin(c+d*x)+(1/2)*b^3*d^-1*tan(c+d*x)^2*sin(c+d*x)+-3*b*a^2*d^-1*sin(c+d*x)+3*a*b^2*d^-1*cos(c+d*x)+3*a*b^2*d^-1*sec(c+d*x)+3*b*a^2*d^-1*arctanh(sin(c+d*x)))
@test integrate((a+b*tan(c+d*x))^3*csc(c+d*x), x) == :(-1*a^3*d^-1*arctanh(cos(c+d*x))+-1//2*b^3*d^-1*arctanh(sin(c+d*x))+(1/2)*b^3*d^-1*sec(c+d*x)*tan(c+d*x)+3*a*b^2*d^-1*sec(c+d*x)+3*b*a^2*d^-1*arctanh(sin(c+d*x)))
@test integrate((a+b*tan(c+d*x))^3*csc(c+d*x)^2, x) == :((1/2)*b^3*d^-1*tan(c+d*x)^2+-1*a^3*d^-1*cot(c+d*x)+3*a*b^2*d^-1*tan(c+d*x)+3*b*a^2*d^-1*log(tan(c+d*x)))
@test integrate((a+b*tan(c+d*x))^3*csc(c+d*x)^3, x) == :((1/2)*b^3*d^-1*arctanh(sin(c+d*x))+-1//2*a^3*d^-1*arctanh(cos(c+d*x))+(1/2)*b^3*d^-1*sec(c+d*x)*tan(c+d*x)+-3*a*b^2*d^-1*arctanh(cos(c+d*x))+-3*b*a^2*d^-1*csc(c+d*x)+3*a*b^2*d^-1*sec(c+d*x)+3*b*a^2*d^-1*arctanh(sin(c+d*x))+-1//2*a^3*d^-1*cot(c+d*x)*csc(c+d*x))
@test integrate((a+b*tan(c+d*x))^3*csc(c+d*x)^4, x) == :((1/2)*b^3*d^-1*tan(c+d*x)^2+-1//3*a^3*d^-1*cot(c+d*x)^3+b*d^-1*(b^2+3*a^2)*log(tan(c+d*x))+-1*a*d^-1*(a^2+3*b^2)*cot(c+d*x)+3*a*b^2*d^-1*tan(c+d*x)+-3//2*b*a^2*d^-1*cot(c+d*x)^2)
@test integrate((a+b*tan(c+d*x))^3*csc(c+d*x)^5, x) == :(-3//2*b^3*d^-1*csc(c+d*x)+-3//8*a^3*d^-1*arctanh(cos(c+d*x))+3//2*b^3*d^-1*arctanh(sin(c+d*x))+(1/2)*b^3*d^-1*sec(c+d*x)^2*csc(c+d*x)+-1*b*a^2*d^-1*csc(c+d*x)^3+-3*b*a^2*d^-1*csc(c+d*x)+3*b*a^2*d^-1*arctanh(sin(c+d*x))+-9//2*a*b^2*d^-1*arctanh(cos(c+d*x))+-3//8*a^3*d^-1*cot(c+d*x)*csc(c+d*x)+-1//4*a^3*d^-1*csc(c+d*x)^3*cot(c+d*x)+9//2*a*b^2*d^-1*sec(c+d*x)+-3//2*a*b^2*d^-1*csc(c+d*x)^2*sec(c+d*x))
@test integrate((a+b*tan(c+d*x))^3*csc(c+d*x)^6, x) == :((1/2)*b^3*d^-1*tan(c+d*x)^2+-1//5*a^3*d^-1*cot(c+d*x)^5+b*d^-1*(2*b^2+3*a^2)*log(tan(c+d*x))+-1*a*d^-1*(a^2+6*b^2)*cot(c+d*x)+3*a*b^2*d^-1*tan(c+d*x)+-3//4*b*a^2*d^-1*cot(c+d*x)^4+-1//2*b*d^-1*cot(c+d*x)^2*(b^2+6*a^2)+-1//3*a*d^-1*cot(c+d*x)^3*(2*a^2+3*b^2))
@test integrate((a+b*tan(c+d*x))^4*sin(c+d*x)^3, x) == :(-1*a^4*d^-1*cos(c+d*x)+-3*b^4*d^-1*cos(c+d*x)+-3*b^4*d^-1*sec(c+d*x)+1//3*a^4*d^-1*cos(c+d*x)^3+1//3*b^4*d^-1*cos(c+d*x)^3+1//3*b^4*d^-1*sec(c+d*x)^3+-10*a*b^3*d^-1*arctanh(sin(c+d*x))+-4*b*a^3*d^-1*sin(c+d*x)+-2*a^2*b^2*d^-1*cos(c+d*x)^3+4*b*a^3*d^-1*arctanh(sin(c+d*x))+6*a^2*b^2*d^-1*sec(c+d*x)+10*a*b^3*d^-1*sin(c+d*x)+12*a^2*b^2*d^-1*cos(c+d*x)+-4//3*b*a^3*d^-1*sin(c+d*x)^3+10//3*a*b^3*d^-1*sin(c+d*x)^3+2*a*b^3*d^-1*sin(c+d*x)^3*tan(c+d*x)^2)
@test integrate((a+b*tan(c+d*x))^4*sin(c+d*x)^2, x) == :(x*((1/2)*a^4+5//2*b^4+-9*a^2*b^2)+5//6*b^4*d^-1*tan(c+d*x)^3+(1/2)*b^2*d^-1*(-5*b^2+18*a^2)*tan(c+d*x)+4*a*b^3*d^-1*tan(c+d*x)^2+-1//2*d^-1*(a+b*tan(c+d*x))^4*cos(c+d*x)*sin(c+d*x)+-4*a*b*d^-1*(a^2+-2*b^2)*log(cos(c+d*x)))
@test integrate((a+b*tan(c+d*x))^4*sin(c+d*x), x) == :(-1*a^4*d^-1*cos(c+d*x)+-1*b^4*d^-1*cos(c+d*x)+-2*b^4*d^-1*sec(c+d*x)+1//3*b^4*d^-1*sec(c+d*x)^3+-6*a*b^3*d^-1*arctanh(sin(c+d*x))+-4*b*a^3*d^-1*sin(c+d*x)+4*b*a^3*d^-1*arctanh(sin(c+d*x))+6*a*b^3*d^-1*sin(c+d*x)+6*a^2*b^2*d^-1*cos(c+d*x)+6*a^2*b^2*d^-1*sec(c+d*x)+2*a*b^3*d^-1*tan(c+d*x)^2*sin(c+d*x))
@test integrate((a+b*tan(c+d*x))^4*csc(c+d*x), x) == :(-1*a^4*d^-1*arctanh(cos(c+d*x))+-1*b^4*d^-1*sec(c+d*x)+1//3*b^4*d^-1*sec(c+d*x)^3+-2*a*b^3*d^-1*arctanh(sin(c+d*x))+4*b*a^3*d^-1*arctanh(sin(c+d*x))+6*a^2*b^2*d^-1*sec(c+d*x)+2*a*b^3*d^-1*sec(c+d*x)*tan(c+d*x))
@test integrate((a+b*tan(c+d*x))^4*csc(c+d*x)^2, x) == :(-1*a^4*d^-1*cot(c+d*x)+1//3*b^4*d^-1*tan(c+d*x)^3+2*a*b^3*d^-1*tan(c+d*x)^2+4*b*a^3*d^-1*log(tan(c+d*x))+6*a^2*b^2*d^-1*tan(c+d*x))
@test integrate((a+b*tan(c+d*x))^4*csc(c+d*x)^3, x) == :(-1//2*a^4*d^-1*arctanh(cos(c+d*x))+1//3*b^4*d^-1*sec(c+d*x)^3+-6*a^2*b^2*d^-1*arctanh(cos(c+d*x))+-4*b*a^3*d^-1*csc(c+d*x)+2*a*b^3*d^-1*arctanh(sin(c+d*x))+4*b*a^3*d^-1*arctanh(sin(c+d*x))+6*a^2*b^2*d^-1*sec(c+d*x)+-1//2*a^4*d^-1*cot(c+d*x)*csc(c+d*x)+2*a*b^3*d^-1*sec(c+d*x)*tan(c+d*x))
@test integrate((a+b*tan(c+d*x))^4*csc(c+d*x)^4, x) == :(-1//3*a^4*d^-1*cot(c+d*x)^3+1//3*b^4*d^-1*tan(c+d*x)^3+b^2*d^-1*(b^2+6*a^2)*tan(c+d*x)+-1*a^2*d^-1*(a^2+6*b^2)*cot(c+d*x)+-2*b*a^3*d^-1*cot(c+d*x)^2+2*a*b^3*d^-1*tan(c+d*x)^2+4*a*b*d^-1*(a^2+b^2)*log(tan(c+d*x)))
@test integrate((a+b*tan(c+d*x))^4*csc(c+d*x)^5, x) == :(b^4*d^-1*sec(c+d*x)+-1*b^4*d^-1*arctanh(cos(c+d*x))+-3//8*a^4*d^-1*arctanh(cos(c+d*x))+1//3*b^4*d^-1*sec(c+d*x)^3+-9*a^2*b^2*d^-1*arctanh(cos(c+d*x))+-6*a*b^3*d^-1*csc(c+d*x)+-4*b*a^3*d^-1*csc(c+d*x)+4*b*a^3*d^-1*arctanh(sin(c+d*x))+6*a*b^3*d^-1*arctanh(sin(c+d*x))+9*a^2*b^2*d^-1*sec(c+d*x)+-4//3*b*a^3*d^-1*csc(c+d*x)^3+-3//8*a^4*d^-1*cot(c+d*x)*csc(c+d*x)+-1//4*a^4*d^-1*csc(c+d*x)^3*cot(c+d*x)+-3*a^2*b^2*d^-1*csc(c+d*x)^2*sec(c+d*x)+2*a*b^3*d^-1*sec(c+d*x)^2*csc(c+d*x))
@test integrate((a+b*tan(c+d*x))^4*csc(c+d*x)^6, x) == :(-1*d^-1*(a^4+b^4+12*a^2*b^2)*cot(c+d*x)+-1//5*a^4*d^-1*cot(c+d*x)^5+1//3*b^4*d^-1*tan(c+d*x)^3+-1*b*a^3*d^-1*cot(c+d*x)^4+2*a*b^3*d^-1*tan(c+d*x)^2+2*b^2*d^-1*(b^2+3*a^2)*tan(c+d*x)+-2//3*a^2*d^-1*cot(c+d*x)^3*(a^2+3*b^2)+-2*a*b*d^-1*cot(c+d*x)^2*(b^2+2*a^2)+4*a*b*d^-1*(a^2+2*b^2)*log(tan(c+d*x)))
@test integrate((a+b*tan(c+d*x))^4*csc(c+d*x)^7, x) == :(-5//2*b^4*d^-1*arctanh(cos(c+d*x))+-5//16*a^4*d^-1*arctanh(cos(c+d*x))+5//2*b^4*d^-1*sec(c+d*x)+5//6*b^4*d^-1*sec(c+d*x)^3+-10*a*b^3*d^-1*csc(c+d*x)+-4*b*a^3*d^-1*csc(c+d*x)+4*b*a^3*d^-1*arctanh(sin(c+d*x))+10*a*b^3*d^-1*arctanh(sin(c+d*x))+-45//4*a^2*b^2*d^-1*arctanh(cos(c+d*x))+-10//3*a*b^3*d^-1*csc(c+d*x)^3+-5//16*a^4*d^-1*cot(c+d*x)*csc(c+d*x)+-5//24*a^4*d^-1*csc(c+d*x)^3*cot(c+d*x)+-4//3*b*a^3*d^-1*csc(c+d*x)^3+-4//5*b*a^3*d^-1*csc(c+d*x)^5+-1//2*b^4*d^-1*csc(c+d*x)^2*sec(c+d*x)^3+-1//6*a^4*d^-1*csc(c+d*x)^5*cot(c+d*x)+45//4*a^2*b^2*d^-1*sec(c+d*x)+2*a*b^3*d^-1*csc(c+d*x)^3*sec(c+d*x)^2+-15//4*a^2*b^2*d^-1*csc(c+d*x)^2*sec(c+d*x)+-3//2*a^2*b^2*d^-1*csc(c+d*x)^4*sec(c+d*x))
@test integrate((a+b*tan(c+d*x))^-1*sin(c+d*x)^5, x) == :(b*d^-1*(5*a^2+5*b^2)^-1*sin(c+d*x)^5+-1*a*d^-1*(a^2+b^2)^-1*cos(c+d*x)+-1*a*d^-1*(5*a^2+5*b^2)^-1*cos(c+d*x)^5+2*a*d^-1*(3*a^2+3*b^2)^-1*cos(c+d*x)^3+a*b^2*d^-1*(a^2+b^2)^-2*cos(c+d*x)+b*a^4*d^-1*(a^2+b^2)^-3*sin(c+d*x)+b*a^5*d^-1*(a^2+b^2)^-7//2*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x)))+a^3*b^2*d^-1*(a^2+b^2)^-3*cos(c+d*x)+-1//3*a*b^2*d^-1*(a^2+b^2)^-2*cos(c+d*x)^3+1//3*b*a^2*d^-1*(a^2+b^2)^-2*sin(c+d*x)^3)
@test integrate((a+b*tan(c+d*x))^-1*sin(c+d*x)^4, x) == :(d^-1*(4*a^2+4*b^2)^-1*cos(c+d*x)^4*(b+a*tan(c+d*x))+-1//8*d^-1*(a^2+b^2)^-2*cos(c+d*x)^2*(4*b*(b^2+2*a^2)+a*(b^2+5*a^2)*tan(c+d*x))+1//8*a*x*(a^2+b^2)^-3*(-1*b^4+3*a^4+-6*a^2*b^2)+b*a^4*d^-1*(a^2+b^2)^-3*log(a*cos(c+d*x)+b*sin(c+d*x)))
@test integrate((a+b*tan(c+d*x))^-1*sin(c+d*x)^3, x) == :(a*d^-1*(3*a^2+3*b^2)^-1*cos(c+d*x)^3+b*d^-1*(3*a^2+3*b^2)^-1*sin(c+d*x)^3+-1*a*d^-1*(a^2+b^2)^-1*cos(c+d*x)+a*b^2*d^-1*(a^2+b^2)^-2*cos(c+d*x)+b*a^2*d^-1*(a^2+b^2)^-2*sin(c+d*x)+b*a^3*d^-1*(a^2+b^2)^-5//2*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x))))
@test integrate((a+b*tan(c+d*x))^-1*sin(c+d*x)^2, x) == :((1/2)*a*x*(a^2+b^2)^-2*(a^2+-1*b^2)+-1*d^-1*(2*a^2+2*b^2)^-1*cos(c+d*x)^2*(b+a*tan(c+d*x))+b*a^2*d^-1*(a^2+b^2)^-2*log(a*cos(c+d*x)+b*sin(c+d*x)))
@test integrate((a+b*tan(c+d*x))^-1*sin(c+d*x), x) == :(b*d^-1*(a^2+b^2)^-1*sin(c+d*x)+-1*a*d^-1*(a^2+b^2)^-1*cos(c+d*x)+a*b*d^-1*(a^2+b^2)^-3//2*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x))))
@test integrate((a+b*tan(c+d*x))^-1*csc(c+d*x), x) == :(-1*a^-1*d^-1*arctanh(cos(c+d*x))+b*a^-1*d^-1*(a^2+b^2)^-1//2*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x))))
@test integrate((a+b*tan(c+d*x))^-1*csc(c+d*x)^2, x) == :(-1*a^-1*d^-1*cot(c+d*x)+b*a^-2*d^-1*log(a+b*tan(c+d*x))+-1*b*a^-2*d^-1*log(tan(c+d*x)))
@test integrate((a+b*tan(c+d*x))^-1*csc(c+d*x)^3, x) == :(-1//2*a^-1*d^-1*arctanh(cos(c+d*x))+b*a^-2*d^-1*csc(c+d*x)+-1*a^-3*b^2*d^-1*arctanh(cos(c+d*x))+-1//2*a^-1*d^-1*cot(c+d*x)*csc(c+d*x)+b*a^-3*d^-1*(a^2+b^2)^(1/2)*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x))))
@test integrate((a+b*tan(c+d*x))^-1*csc(c+d*x)^4, x) == :(-1//3*a^-1*d^-1*cot(c+d*x)^3+(1/2)*b*a^-2*d^-1*cot(c+d*x)^2+-1*a^-3*d^-1*(a^2+b^2)*cot(c+d*x)+b*a^-4*d^-1*(a^2+b^2)*log(a+b*tan(c+d*x))+-1*b*a^-4*d^-1*(a^2+b^2)*log(tan(c+d*x)))
@test integrate((a+b*tan(c+d*x))^-1*csc(c+d*x)^6, x) == :(-1//5*a^-1*d^-1*cot(c+d*x)^5+-1*a^-5*d^-1*(a^2+b^2)^2*cot(c+d*x)+-1//3*a^-3*d^-1*cot(c+d*x)^3*(b^2+2*a^2)+1//4*b*a^-2*d^-1*cot(c+d*x)^4+b*a^-6*d^-1*(a^2+b^2)^2*log(a+b*tan(c+d*x))+(1/2)*b*a^-4*d^-1*cot(c+d*x)^2*(b^2+2*a^2)+-1*b*a^-6*d^-1*(a^2+b^2)^2*log(tan(c+d*x)))
@test integrate((a+b*tan(c+d*x))^-2*sin(c+d*x)^6, x) == :(1//16*x*(a^2+b^2)^-5*(b^8+5*a^8+-80*a^6*b^2+8*a^2*b^6+50*a^4*b^4)+-1//6*d^-1*(a^2+b^2)^-2*cos(c+d*x)^6*((a^2+-1*b^2)*tan(c+d*x)+2*a*b)+-1//16*d^-1*(a^2+b^2)^-4*cos(c+d*x)^2*((-1*b^6+11*a^6+-43*a^4*b^2+-7*a^2*b^4)*tan(c+d*x)+48*b*a^5)+1//24*d^-1*(a^2+b^2)^-3*cos(c+d*x)^4*((-7*b^4+13*a^4+-18*a^2*b^2)*tan(c+d*x)+12*a*b*(b^2+3*a^2))+-1*b*a^6*d^-1*(a+b*tan(c+d*x))^-1*(a^2+b^2)^-4+2*b*a^5*d^-1*(a^2+b^2)^-5*(a^2+-3*b^2)*log(a*cos(c+d*x)+b*sin(c+d*x)))
@test integrate((a+b*tan(c+d*x))^-2*sin(c+d*x)^4, x) == :(1//8*x*(a^2+b^2)^-4*(b^6+3*a^6+-33*a^4*b^2+13*a^2*b^4)+-1//8*d^-1*(a^2+b^2)^-3*cos(c+d*x)^2*((-1*b^4+5*a^4+-12*a^2*b^2)*tan(c+d*x)+16*b*a^3)+1//4*d^-1*(a^2+b^2)^-2*cos(c+d*x)^4*((a^2+-1*b^2)*tan(c+d*x)+2*a*b)+-1*b*a^4*d^-1*(a+b*tan(c+d*x))^-1*(a^2+b^2)^-3+2*b*a^3*d^-1*(a^2+b^2)^-4*(a^2+-2*b^2)*log(a*cos(c+d*x)+b*sin(c+d*x)))
@test integrate((a+b*tan(c+d*x))^-2*sin(c+d*x)^2, x) == :((1/2)*x*(a^2+b^2)^-3*(a^4+b^4+-6*a^2*b^2)+-1//2*d^-1*(a^2+b^2)^-2*cos(c+d*x)^2*((a^2+-1*b^2)*tan(c+d*x)+2*a*b)+-1*b*a^2*d^-1*(a+b*tan(c+d*x))^-1*(a^2+b^2)^-2+2*a*b*d^-1*(a^2+b^2)^-3*(a^2+-1*b^2)*log(a*cos(c+d*x)+b*sin(c+d*x)))
@test integrate((a+b*tan(c+d*x))^-2*csc(c+d*x)^2, x) == :(-1*a^-2*d^-1*cot(c+d*x)+-1*b*a^-2*d^-1*(a+b*tan(c+d*x))^-1+-2*b*a^-3*d^-1*log(tan(c+d*x))+2*b*a^-3*d^-1*log(a+b*tan(c+d*x)))
@test integrate((a+b*tan(c+d*x))^-2*csc(c+d*x)^4, x) == :(-1//3*a^-2*d^-1*cot(c+d*x)^3+b*a^-3*d^-1*cot(c+d*x)^2+-1*a^-4*d^-1*(a^2+3*b^2)*cot(c+d*x)+-1*b*a^-4*d^-1*(a+b*tan(c+d*x))^-1*(a^2+b^2)+-2*b*a^-5*d^-1*(a^2+2*b^2)*log(tan(c+d*x))+2*b*a^-5*d^-1*(a^2+2*b^2)*log(a+b*tan(c+d*x)))
@test integrate((a+b*tan(c+d*x))^-2*csc(c+d*x)^6, x) == :(-1//5*a^-2*d^-1*cot(c+d*x)^5+(1/2)*b*a^-3*d^-1*cot(c+d*x)^4+-1//3*a^-4*d^-1*cot(c+d*x)^3*(2*a^2+3*b^2)+-1*b*a^-6*d^-1*(a+b*tan(c+d*x))^-1*(a^2+b^2)^2+-1*a^-6*d^-1*(a^2+b^2)*(a^2+5*b^2)*cot(c+d*x)+2*b*a^-5*d^-1*cot(c+d*x)^2*(a^2+b^2)+-2*b*a^-7*d^-1*(a^2+b^2)*(a^2+3*b^2)*log(tan(c+d*x))+2*b*a^-7*d^-1*(a^2+b^2)*(a^2+3*b^2)*log(a+b*tan(c+d*x)))
@test integrate((a+b*tan(c+d*x))^-2*sin(c+d*x)^5, x) == :(-3//64*d^-1*(a+b*tan(c+d*x))^-2*(a*cos(c+d*x)+b*sin(c+d*x))^2*sec(c+d*x)^2*((a^2+b^2)^-3*(8*a^4+8*b^4+-48*a^2*b^2)*cos(c+d*x)+-1//3*(a^2+b^2)^-2*(-4*b^2+4*a^2)*cos(3c+3*d*x)+a*(a^2+b^2)^-3*(a*cos(c+d*x)+b*sin(c+d*x))^-1*(a^4+5*b^4+-10*a^2*b^2)+10*b*(a^2+b^2)^-7//2*(b^4+5*a^4+-10*a^2*b^2)*arctanh((a^2+b^2)^-1//2*(-1b+a*tan((1/2)*c+(1/2)*d*x)))+8//3*a*b*(a^2+b^2)^-2*sin(3c+3*d*x)+-32*a*b*(a^2+b^2)^-3*(a^2+-1*b^2)*sin(c+d*x))+1//64*d^-1*(a+b*tan(c+d*x))^-2*(a*cos(c+d*x)+b*sin(c+d*x))^2*sec(c+d*x)^2*(-1*(a^2+b^2)^-4*(-12*b^6+12*a^6+-180*a^4*b^2+180*a^2*b^4)*cos(c+d*x)+-1//5*(a^2+b^2)^-2*(-4*b^2+4*a^2)*cos(5c+5*d*x)+1//3*(a^2+b^2)^-3*(8*a^4+8*b^4+-48*a^2*b^2)*cos(3c+3*d*x)+-1*a*(a^2+b^2)^-4*(a*cos(c+d*x)+b*sin(c+d*x))^-1*(a^6+-7*b^6+-21*a^4*b^2+35*a^2*b^4)+14*b*(a^2+b^2)^-9//2*(b^6+-7*a^6+-21*a^2*b^4+35*a^4*b^2)*arctanh((a^2+b^2)^-1//2*(-1b+a*tan((1/2)*c+(1/2)*d*x)))+8//5*a*b*(a^2+b^2)^-2*sin(5c+5*d*x)+24*a*b*(a^2+b^2)^-4*(3*a^4+3*b^4+-10*a^2*b^2)*sin(c+d*x)+-32//3*a*b*(a^2+b^2)^-3*(a^2+-1*b^2)*sin(3c+3*d*x))+1//64*d^-1*(a+b*tan(c+d*x))^-2*(a^2+b^2)^-3*sec(c+d*x)^2*(a*cos(c+d*x)+b*sin(c+d*x))*(-1*(a^2+b^2)*(3*a*(a^2+-3*b^2)+-2*b*(a^2+b^2)*sin(2c+2*d*x)+2*a*(a^2+b^2)*cos(2c+2*d*x))+6*b*(a^2+b^2)^(1/2)*(b^2+-3*a^2)*(a*cos(c+d*x)+b*sin(c+d*x))*arctanh((a^2+b^2)^-1//2*(-1b+a*tan((1/2)*c+(1/2)*d*x))))+5//64*d^-1*(a+b*tan(c+d*x))^-2*(a^2+b^2)^-2*sec(c+d*x)^2*(a*(a^2+b^2)+2*b*(a^2+b^2)^(1/2)*(a*cos(c+d*x)+b*sin(c+d*x))*arctanh((a^2+b^2)^-1//2*(-1b+a*tan((1/2)*c+(1/2)*d*x))))*(a*cos(c+d*x)+b*sin(c+d*x)))
@test integrate((a+b*tan(c+d*x))^-2*sin(c+d*x)^3, x) == :(-1*a^2*d^-1*(a^2+b^2)^-2*cos(c+d*x)+-1//3*b^2*d^-1*(a^2+b^2)^-2*cos(c+d*x)^3+1//3*a^2*d^-1*(a^2+b^2)^-2*cos(c+d*x)^3+a^3*b^2*d^-1*(a^2+b^2)^-3*(a*cos(c+d*x)+b*sin(c+d*x))^-1+-3*a^2*b^3*d^-1*(a^2+b^2)^-7//2*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x)))+-2*a*b^3*d^-1*(a^2+b^2)^-3*sin(c+d*x)+2*b*a^3*d^-1*(a^2+b^2)^-3*sin(c+d*x)+2*b*a^4*d^-1*(a^2+b^2)^-7//2*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x)))+4*a^2*b^2*d^-1*(a^2+b^2)^-3*cos(c+d*x)+2//3*a*b*d^-1*(a^2+b^2)^-2*sin(c+d*x)^3)
@test integrate((a+b*tan(c+d*x))^-2*sin(c+d*x), x) == :(b^2*d^-1*(a^2+b^2)^-2*cos(c+d*x)+-1*a^2*d^-1*(a^2+b^2)^-2*cos(c+d*x)+-1*b^3*d^-1*(a^2+b^2)^-5//2*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x)))+a*b^2*d^-1*(a^2+b^2)^-2*(a*cos(c+d*x)+b*sin(c+d*x))^-1+2*a*b*d^-1*(a^2+b^2)^-2*sin(c+d*x)+2*b*a^2*d^-1*(a^2+b^2)^-5//2*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x))))
@test integrate((a+b*tan(c+d*x))^-2*csc(c+d*x), x) == :(-1*a^-2*d^-1*arctanh(cos(c+d*x))+b*d^-1*(a^2+b^2)^-3//2*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x)))+b*a^-2*d^-1*(a^2+b^2)^-1//2*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x)))+a^-1*b^2*d^-1*(a^2+b^2)^-1*(a*cos(c+d*x)+b*sin(c+d*x))^-1)
@test integrate((a+b*tan(c+d*x))^-2*csc(c+d*x)^3, x) == :(-1//2*a^-2*d^-1*arctanh(cos(c+d*x))+a^-3*b^2*d^-1*(a*cos(c+d*x)+b*sin(c+d*x))^-1+-3*a^-4*b^2*d^-1*arctanh(cos(c+d*x))+2*b*a^-3*d^-1*csc(c+d*x)+-1//2*a^-2*d^-1*cot(c+d*x)*csc(c+d*x)+-1*b*a^-2*d^-1*(a^2+b^2)^-1//2*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x)))+3*b*a^-4*d^-1*(a^2+b^2)^(1/2)*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x))))
@test integrate((a+b*tan(c+d*x))^-3*sin(c+d*x)^6, x) == :(-1//6*d^-1*(a^2+b^2)^-3*cos(c+d*x)^6*(b*(-1*b^2+3*a^2)+a*(a^2+-3*b^2)*tan(c+d*x))+1//16*a*x*(a^2+b^2)^-6*(-3*b^8+5*a^8+-180*a^6*b^2+-68*a^2*b^6+390*a^4*b^4)+1//24*d^-1*(a^2+b^2)^-4*cos(c+d*x)^4*(6*b*(-1*b^4+9*a^4+-4*a^2*b^2)+a*(-3*b^4+13*a^4+-62*a^2*b^2)*tan(c+d*x))+-1//2*b*a^6*d^-1*(a+b*tan(c+d*x))^-2*(a^2+b^2)^-4+-1//16*a*d^-1*(a^2+b^2)^-5*cos(c+d*x)^2*((3*b^6+11*a^6+-119*a^4*b^2+65*a^2*b^4)*tan(c+d*x)+24*b*a^3*(-5*b^2+3*a^2))+b*a^4*d^-1*(a^2+b^2)^-6*(3*a^4+15*b^4+-22*a^2*b^2)*log(a*cos(c+d*x)+b*sin(c+d*x))+-2*b*a^5*d^-1*(a+b*tan(c+d*x))^-1*(a^2+b^2)^-5*(a^2+-3*b^2))
@test integrate((a+b*tan(c+d*x))^-3*sin(c+d*x)^4, x) == :(1//4*d^-1*(a^2+b^2)^-3*cos(c+d*x)^4*(b*(-1*b^2+3*a^2)+a*(a^2+-3*b^2)*tan(c+d*x))+3//8*a*x*(a^2+b^2)^-5*(a^6+-3*b^6+-25*a^4*b^2+35*a^2*b^4)+-1//2*b*a^4*d^-1*(a+b*tan(c+d*x))^-2*(a^2+b^2)^-3+-1//8*a*d^-1*(a^2+b^2)^-4*cos(c+d*x)^2*((5*a^4+9*b^4+-34*a^2*b^2)*tan(c+d*x)+24*a*b*(a^2+-1*b^2))+-2*b*a^3*d^-1*(a+b*tan(c+d*x))^-1*(a^2+b^2)^-4*(a^2+-2*b^2)+3*b*a^2*d^-1*(a^2+b^2)^-5*(a^4+2*b^4+-5*a^2*b^2)*log(a*cos(c+d*x)+b*sin(c+d*x)))
@test integrate((a+b*tan(c+d*x))^-3*sin(c+d*x)^2, x) == :((1/2)*a*x*(a^2+b^2)^-4*(a^4+9*b^4+-14*a^2*b^2)+-1//2*d^-1*(a^2+b^2)^-3*cos(c+d*x)^2*(b*(-1*b^2+3*a^2)+a*(a^2+-3*b^2)*tan(c+d*x))+b*d^-1*(a^2+b^2)^-4*(b^4+3*a^4+-8*a^2*b^2)*log(a*cos(c+d*x)+b*sin(c+d*x))+-1//2*b*a^2*d^-1*(a+b*tan(c+d*x))^-2*(a^2+b^2)^-2+-2*a*b*d^-1*(a+b*tan(c+d*x))^-1*(a^2+b^2)^-3*(a^2+-1*b^2))
@test integrate((a+b*tan(c+d*x))^-3*csc(c+d*x)^2, x) == :(-1*a^-3*d^-1*cot(c+d*x)+-3*b*a^-4*d^-1*log(tan(c+d*x))+-2*b*a^-3*d^-1*(a+b*tan(c+d*x))^-1+3*b*a^-4*d^-1*log(a+b*tan(c+d*x))+-1//2*b*a^-2*d^-1*(a+b*tan(c+d*x))^-2)
@test integrate((a+b*tan(c+d*x))^-3*csc(c+d*x)^4, x) == :(-1//3*a^-3*d^-1*cot(c+d*x)^3+-1*a^-5*d^-1*(a^2+6*b^2)*cot(c+d*x)+3//2*b*a^-4*d^-1*cot(c+d*x)^2+b*a^-6*d^-1*(3*a^2+10*b^2)*log(a+b*tan(c+d*x))+-1*b*a^-6*d^-1*(3*a^2+10*b^2)*log(tan(c+d*x))+-2*b*a^-5*d^-1*(a+b*tan(c+d*x))^-1*(a^2+2*b^2)+-1//2*b*a^-4*d^-1*(a+b*tan(c+d*x))^-2*(a^2+b^2))
@test integrate((a+b*tan(c+d*x))^-3*csc(c+d*x)^6, x) == :(-1//5*a^-3*d^-1*cot(c+d*x)^5+-1*a^-7*d^-1*(a^4+15*b^4+12*a^2*b^2)*cot(c+d*x)+-1//3*a^-5*d^-1*cot(c+d*x)^3*(2*a^2+6*b^2)+3//4*b*a^-4*d^-1*cot(c+d*x)^4+b*a^-8*d^-1*(3*a^4+21*b^4+20*a^2*b^2)*log(a+b*tan(c+d*x))+b*a^-6*d^-1*cot(c+d*x)^2*(3*a^2+5*b^2)+-1*b*a^-8*d^-1*(3*a^4+21*b^4+20*a^2*b^2)*log(tan(c+d*x))+-1//2*b*a^-6*d^-1*(a+b*tan(c+d*x))^-2*(a^2+b^2)^2+-2*b*a^-7*d^-1*(a+b*tan(c+d*x))^-1*(a^2+b^2)*(a^2+3*b^2))
@test integrate((a+b*tan(c+d*x))^-3*sin(c+d*x)^5, x) == :(1//64*d^-1*(a+b*tan(c+d*x))^-3*(a*cos(c+d*x)+b*sin(c+d*x))^3*sec(c+d*x)^3*((a^2+b^2)^-3*(a*cos(c+d*x)+b*sin(c+d*x))^-2*(-2*b*((a^2+b^2)^2*sin(3c+3*d*x)+-5*b^2*(b^2+-5*a^2)*sin(c+d*x))+2*a*(a^2+b^2)^2*cos(3c+3*d*x)+10*a^3*(a^2+-5*b^2)*cos(c+d*x))+60*a*b*(a^2+b^2)^-7//2*(a^2+-1*b^2)*arctanh((a^2+b^2)^-1//2*(-1b+a*tan((1/2)*c+(1/2)*d*x))))+1//64*d^-1*(a+b*tan(c+d*x))^-3*(a*cos(c+d*x)+b*sin(c+d*x))^3*sec(c+d*x)^3*(-1*(a^2+b^2)^-4*(a*cos(c+d*x)+b*sin(c+d*x))^-1*(-3*b^6+6*a^6+-87*a^4*b^2+80*a^2*b^4)+-24*a*(a^2+b^2)^-4*(a^4+5*b^4+-10*a^2*b^2)*cos(c+d*x)+24*b*(a^2+b^2)^-4*(b^4+5*a^4+-10*a^2*b^2)*sin(c+d*x)+8//3*a*(a^2+b^2)^-3*(a^2+-3*b^2)*cos(3c+3*d*x)+8//3*b*(a^2+b^2)^-3*(b^2+-3*a^2)*sin(3c+3*d*x)+-1*b*(a^2+b^2)^-3*(a*cos(c+d*x)+b*sin(c+d*x))^-2*(-3*a^4+-3*b^4+10*a^2*b^2)*sin(c+d*x)+-70*a*b*(a^2+b^2)^-9//2*(3*a^4+3*b^4+-10*a^2*b^2)*arctanh((a^2+b^2)^-1//2*(-1b+a*tan((1/2)*c+(1/2)*d*x))))+-1//320*d^-1*(a+b*tan(c+d*x))^-3*(a^2+b^2)^-6*sec(c+d*x)^3*(a*cos(c+d*x)+b*sin(c+d*x))*((a^2+b^2)*(a^9*cos(7c+7*d*x)+-1*b^9*sin(7c+7*d*x)+-21*b^9*sin(3c+3*d*x)+-3*a^9*cos(5c+5*d*x)+-3*b^9*sin(5c+5*d*x)+105*b^9*sin(c+d*x)+a*b^8*cos(7c+7*d*x)+-1*b*a^8*sin(7c+7*d*x)+-2730*a^2*b^7*sin(c+d*x)+-1470*a^6*b^3*sin(c+d*x)+-105*b*a^8*sin(c+d*x)+-105*b*a^8*sin(3c+3*d*x)+-6*a^4*b^5*sin(7c+7*d*x)+-4*a^2*b^7*sin(7c+7*d*x)+-4*a^6*b^3*sin(7c+7*d*x)+4*a^3*b^6*cos(7c+7*d*x)+4*a^7*b^2*cos(7c+7*d*x)+6*a^5*b^4*cos(7c+7*d*x)+9*a*b^8*cos(5c+5*d*x)+9*b*a^8*sin(5c+5*d*x)+18*a^4*b^5*sin(5c+5*d*x)+18*a^5*b^4*cos(5c+5*d*x)+24*a^3*b^6*cos(5c+5*d*x)+24*a^6*b^3*sin(5c+5*d*x)+105*a*(a^8+-1*b^8+-26*a^6*b^2+-14*a^2*b^6+56*a^4*b^4)*cos(c+d*x)+168*a^2*b^7*sin(3c+3*d*x)+294*a^4*b^5*sin(3c+3*d*x)+5880*a^4*b^5*sin(c+d*x)+21*a*(a^2+b^2)^2*(a^4+5*b^4+-10*a^2*b^2)*cos(3c+3*d*x))+1260*a*b*(a^2+b^2)^(1/2)*(a*cos(c+d*x)+b*sin(c+d*x))^2*(a^6+-1*b^6+-7*a^4*b^2+7*a^2*b^4)*arctanh((a^2+b^2)^-1//2*(-1b+a*tan((1/2)*c+(1/2)*d*x))))+3//64*d^-1*(a+b*tan(c+d*x))^-3*(a^2+b^2)^-3*sec(c+d*x)^3*(a*cos(c+d*x)+b*sin(c+d*x))*((a^2+b^2)*(a*(-1*b^2+2*a^2)*cos(c+d*x)+b*(a^2+-2*b^2)*sin(c+d*x))+6*a*b*(a^2+b^2)^(1/2)*(a*cos(c+d*x)+b*sin(c+d*x))^2*arctanh((a^2+b^2)^-1//2*(-1b+a*tan((1/2)*c+(1/2)*d*x)))))
@test integrate((a+b*tan(c+d*x))^-3*sin(c+d*x)^3, x) == :(b^3*d^-1*(a^2+b^2)^-3*sin(c+d*x)+-1*a^3*d^-1*(a^2+b^2)^-3*cos(c+d*x)+-1//3*b^3*d^-1*(a^2+b^2)^-3*sin(c+d*x)^3+1//3*a^3*d^-1*(a^2+b^2)^-3*cos(c+d*x)^3+b*a^2*d^-1*(a^2+b^2)^-3*sin(c+d*x)^3+-1*a*b^2*d^-1*(a^2+b^2)^-3*cos(c+d*x)^3+-9*a^2*b^3*d^-1*(a^2+b^2)^-4*sin(c+d*x)+-3*a*b^4*d^-1*(a^2+b^2)^-4*cos(c+d*x)+-3*a^2*b^4*d^-1*(a^2+b^2)^-4*(a*cos(c+d*x)+b*sin(c+d*x))^-1+3*a*b^5*d^-1*(a^2+b^2)^-9//2*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x)))+3*b*a^4*d^-1*(a^2+b^2)^-4*sin(c+d*x)+3*b*a^5*d^-1*(a^2+b^2)^-9//2*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x)))+3*a^4*b^2*d^-1*(a^2+b^2)^-4*(a*cos(c+d*x)+b*sin(c+d*x))^-1+9*a^3*b^2*d^-1*(a^2+b^2)^-4*cos(c+d*x)+-23//2*a^3*b^3*d^-1*(a^2+b^2)^-9//2*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x)))+(1/2)*a^3*b^3*d^-1*(a^2+b^2)^-4*(a*cos(c+d*x)+b*sin(c+d*x))^-2*(b*cos(c+d*x)+-1*a*sin(c+d*x)))
@test integrate((a+b*tan(c+d*x))^-3*sin(c+d*x), x) == :(-1//4*d^-1*(a+b*tan(c+d*x))^-3*(a^2+b^2)^-7//2*sec(c+d*x)^3*(a*cos(c+d*x)+b*sin(c+d*x))*((a^2+b^2)^(1/2)*(a*(a^2+b^2)^2*cos(3c+3*d*x)+a*(b^4+3*a^4+-26*a^2*b^2)*cos(c+d*x)+-2*b*(a^4+-3*b^4+(a^2+b^2)^2*cos(2c+2*d*x)+13*a^2*b^2)*sin(c+d*x))+12*a*b*(a*cos(c+d*x)+b*sin(c+d*x))^2*(-3*b^2+2*a^2)*arctanh((a^2+b^2)^-1//2*(-1b+a*tan((1/2)*c+(1/2)*d*x)))))
@test integrate((a+b*tan(c+d*x))^-3*csc(c+d*x), x) == :((1/2)*a^-3*d^-1*(a+b*tan(c+d*x))^-3*sec(c+d*x)^2*(a*cos(c+d*x)+b*sin(c+d*x))*(-2*(a*cos(c+d*x)+b*sin(c+d*x))^2*log(cos((1/2)*c+(1/2)*d*x))*sec(c+d*x)+2*(a*cos(c+d*x)+b*sin(c+d*x))^2*log(sin((1/2)*c+(1/2)*d*x))*sec(c+d*x)+-1*a*b^3*(a^2+b^2)^-1*tan(c+d*x)+-2*b*(a^2+b^2)^-5//2*(a*cos(c+d*x)+b*sin(c+d*x))^2*(2*b^4+6*a^4+5*a^2*b^2)*arctanh((a^2+b^2)^-1//2*(-1b+a*tan((1/2)*c+(1/2)*d*x)))*sec(c+d*x)+3*a*b^2*(a^2+b^2)^-2*(b^2+2*a^2)*(a*cos(c+d*x)+b*sin(c+d*x))*sec(c+d*x)))
@test integrate((a+b*tan(c+d*x))^-3*csc(c+d*x)^3, x) == :(-1//2*a^-3*d^-1*arctanh(cos(c+d*x))+-6*a^-5*b^2*d^-1*arctanh(cos(c+d*x))+3*b*a^-4*d^-1*csc(c+d*x)+3*a^-4*b^2*d^-1*(a*cos(c+d*x)+b*sin(c+d*x))^-1+-1//2*a^-3*d^-1*cot(c+d*x)*csc(c+d*x)+(1/2)*a^-3*b^3*d^-1*(a^2+b^2)^-3//2*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x)))+-3*b*a^-3*d^-1*(a^2+b^2)^-1//2*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x)))+6*b*a^-5*d^-1*(a^2+b^2)^(1/2)*arctanh((a^2+b^2)^-1//2*(b*cos(c+d*x)+-1*a*sin(c+d*x)))+(1/2)*a^-3*b^3*d^-1*(a^2+b^2)^-1*(a*cos(c+d*x)+b*sin(c+d*x))^-2*(b*cos(c+d*x)+-1*a*sin(c+d*x)))
@test integrate((a+b*tan(c+d*x))^-4*sin(c+d*x)^4, x) == :(1//8*x*(a^2+b^2)^-6*(3*a^8+3*b^8+-132*a^2*b^6+-132*a^6*b^2+370*a^4*b^4)+-1//8*d^-1*(a^2+b^2)^-5*cos(c+d*x)^2*((-3*b^6+5*a^6+-65*a^4*b^2+55*a^2*b^4)*tan(c+d*x)+16*a*b*(b^4+2*a^4+-5*a^2*b^2))+1//4*d^-1*(a^2+b^2)^-4*cos(c+d*x)^4*((a^4+b^4+-6*a^2*b^2)*tan(c+d*x)+4*a*b*(a^2+-1*b^2))+-1//3*b*a^4*d^-1*(a+b*tan(c+d*x))^-3*(a^2+b^2)^-3+-1*b*a^3*d^-1*(a+b*tan(c+d*x))^-2*(a^2+b^2)^-4*(a^2+-2*b^2)+-3*b*a^2*d^-1*(a+b*tan(c+d*x))^-1*(a^2+b^2)^-5*(a^4+2*b^4+-5*a^2*b^2)+4*a*b*d^-1*(a^2+b^2)^-6*(a^2+-1*b^2)*(a^4+b^4+-8*a^2*b^2)*log(a*cos(c+d*x)+b*sin(c+d*x)))
@test integrate((a+b*tan(c+d*x))^-4*sin(c+d*x)^2, x) == :((1/2)*x*(a^2+b^2)^-5*(a^6+-3*b^6+-25*a^4*b^2+35*a^2*b^4)+-1//2*d^-1*(a^2+b^2)^-4*cos(c+d*x)^2*((a^4+b^4+-6*a^2*b^2)*tan(c+d*x)+4*a*b*(a^2+-1*b^2))+-1*b*d^-1*(a+b*tan(c+d*x))^-1*(a^2+b^2)^-4*(b^4+3*a^4+-8*a^2*b^2)+-1//3*b*a^2*d^-1*(a+b*tan(c+d*x))^-3*(a^2+b^2)^-2+-1*a*b*d^-1*(a+b*tan(c+d*x))^-2*(a^2+b^2)^-3*(a^2+-1*b^2)+4*a*b*d^-1*(a^2+b^2)^-5*(a^4+2*b^4+-5*a^2*b^2)*log(a*cos(c+d*x)+b*sin(c+d*x)))
@test integrate((a+b*tan(c+d*x))^-4*csc(c+d*x)^2, x) == :(-1*a^-4*d^-1*cot(c+d*x)+-1*b*a^-3*d^-1*(a+b*tan(c+d*x))^-2+-4*b*a^-5*d^-1*log(tan(c+d*x))+-3*b*a^-4*d^-1*(a+b*tan(c+d*x))^-1+4*b*a^-5*d^-1*log(a+b*tan(c+d*x))+-1//3*b*a^-2*d^-1*(a+b*tan(c+d*x))^-3)
@test integrate((a+b*tan(c+d*x))^-4*csc(c+d*x)^4, x) == :(-1//3*a^-4*d^-1*cot(c+d*x)^3+-1*a^-6*d^-1*(a^2+10*b^2)*cot(c+d*x)+2*b*a^-5*d^-1*cot(c+d*x)^2+-1*b*a^-6*d^-1*(a+b*tan(c+d*x))^-1*(3*a^2+10*b^2)+-1*b*a^-5*d^-1*(a+b*tan(c+d*x))^-2*(a^2+2*b^2)+-4*b*a^-7*d^-1*(a^2+5*b^2)*log(tan(c+d*x))+4*b*a^-7*d^-1*(a^2+5*b^2)*log(a+b*tan(c+d*x))+-1//3*b*a^-4*d^-1*(a+b*tan(c+d*x))^-3*(a^2+b^2))
@test integrate((a+b*tan(c+d*x))^-4*csc(c+d*x)^6, x) == :(-1//5*a^-4*d^-1*cot(c+d*x)^5+b*a^-5*d^-1*cot(c+d*x)^4+-1*a^-8*d^-1*(a^4+35*b^4+20*a^2*b^2)*cot(c+d*x)+-1//3*a^-6*d^-1*cot(c+d*x)^3*(2*a^2+10*b^2)+-1*b*a^-8*d^-1*(a+b*tan(c+d*x))^-1*(3*a^4+21*b^4+20*a^2*b^2)+-4*b*a^-9*d^-1*(a^4+14*b^4+10*a^2*b^2)*log(tan(c+d*x))+2*b*a^-7*d^-1*cot(c+d*x)^2*(2*a^2+5*b^2)+4*b*a^-9*d^-1*(a^4+14*b^4+10*a^2*b^2)*log(a+b*tan(c+d*x))+-1//3*b*a^-6*d^-1*(a+b*tan(c+d*x))^-3*(a^2+b^2)^2+-1*b*a^-7*d^-1*(a+b*tan(c+d*x))^-2*(a^2+b^2)*(a^2+3*b^2))
@test integrate((1+tan(x))^-1*csc(x), x) == :(-1*arctanh(cos(x))+(1/2)*2^(1/2)*arctanh((1/2)*2^(1/2)*(-1*sin(x)+cos(x))))

max_stars_repo_path: test/4 Trig functions/4.3 Tangent/4.3.1.3 (d sin)^m (a+b tan)^n.jl


text:
const me = 0.511 # MeV
const ec = 0.302811
const αem = ec^2/(4pi)

#################################################
# VECTOR SELF-ENERGY IN DILUTE NON-RELATIVISTIC GAS

# Plasma dispersion relation Z(xi) in a form valid for real arguments
function Zi(x)
	sqrt(pi)*exp(-x^2)
end

function Zr(x)
	-2*dawson(x)
end

# Real part of longitudinal photon polarization tensor in dilute, non-relativistic gas.
function ΠLr_dnr(q0,q,m,Z,n,T)
	s = sqrt(T/m)
	xi = q0/(sqrt(2)*q*s)
	qp = (q0^2 - q^2) / (2*q0*m)
	(
	 ((ec^2 * Z^2 * n) / (sqrt(2)*q*s)) *
	 (Zr(xi*(1 + qp)) - Zr(xi*(1 - qp)))
	 )
end

# Real part of transverse photon polarization tensor in dilute, non-relativistic gas.
function ΠTr_dnr(q0,q,m,Z,n,T)
	wp2 = ec^2 * Z^2 * n / m
	s = sqrt(T/m)
	xi = q0/(sqrt(2)*q*s)
	Q2 = q0^2 - q^2
	qp = Q2 / (2*q0*m)
	(
	 wp2*(1 - (s^2 - Q2/(4*m^2))*(m/q0)*xi*
		  (Zr(xi*(1 + qp)) - Zr(xi*(1 - qp))))
	 )
end

# Imaginary part of longitudinal photon polarization tensor in dilute, non-relativistic gas
function ΠLi_dnr(q0,q,m,Z,n,T)
	s = sqrt(T/m)
	xi = q0/(sqrt(2)*q*s)
	qp = (q0^2 - q^2) / (2*q0*m)
	(
	 ((ec^2 * Z^2 * n) / (sqrt(2)*q*s)) *
	 (Zi(xi*(1 + qp)) - Zi(xi*(1 - qp)))
	 )
end

# Imaginary part of transverse photon polarization tensor in dilute, non-relativistic gas
function ΠTi_dnr(q0,q,m,Z,n,T)
	wp2 = ec^2 * Z^2 * n / m
	s = sqrt(T/m)
	xi = q0/(sqrt(2)*q*s)
	Q2 = q0^2 - q^2
	qp = Q2 / (2*q0*m)
	(
	 -wp2*((s^2 - Q2/(4*m^2))*(m/q0)*xi*
		  (Zi(xi*(1 + qp)) - Zi(xi*(1 - qp))))
	 )
end

#################################################
# VECTOR SELF-ENERGY IN DEGENERATE FERMION GAS

# Imaginary part of longitudinal photon polarization tensor for degenerate electron gas.
function ΠLi_e_degen(q0,q,mue;mf=me)
	Q2 = q0^2 - q^2
	pre = -(ec^2/(4pi)) * (Q2/q^3)
	ep = max(mue,mf + q0,0.5*(q0 + q*sqrt(1 - 4*mf^2/Q2)))
	if ep > mue + q0
		0.
	else
		pre * (-((2*ep^3)/3) + (2*mue^3)/3 + ep^2*q0 + mue^2*q0 - q0^3/3
		 - (ep*Q2)/2 + (mue*Q2)/2 + (q0*Q2)/2)
	end
end

# Imaginary part of transverse photon polarization tensor for degenerate electron gas.
function ΠTi_e_degen(q0,q,mue;mf=me)
	Q2 = q0^2 - q^2
	pre = -(ec^2/(4pi)) * (1/q)
	ep = max(mue,mf + q0,0.5*(q0 + q*sqrt(1 - 4*mf^2/Q2)))
	if ep > mue + q0
		0.
	else
		(pre/12)*(-12*mf^2*(-ep + mue + q0) + 4*(-ep^3 + (mue + q0)^3) -
 (4*q0^2*(-ep^3 + (mue + q0)^3))/q^2 - 6*(-ep + mue + q0)*Q2 +
 (6*q0*(-ep^2 + (mue + q0)^2)*Q2)/q^2 - (3*(-ep + mue + q0)*Q2^2)/q^2)
	end
end

# Real part of longitudinal photon polarization tensor for degenerate electron gas.
# For internal integral definitions, see below.
function ΠLr_e_degen(q0,q,mue;mf=me)
	Q2 = q0^2 - q^2

	pf = sqrt(mue^2 - mf^2)

	ik = q*(-(mue*pf) + mf^2*log((mue + pf)/mf))

	i0m = integ0m(q0,q,mue,pf,mf,Q2)
	i1m = integ1m(q0,q,mue,pf,mf,Q2)
	i2m = integ2m(q0,q,mue,pf,mf,Q2)

	i0p = integ0p(q0,q,mue,pf,mf,Q2)
	i1p = integ1p(q0,q,mue,pf,mf,Q2)
	i2p = integ2p(q0,q,mue,pf,mf,Q2)

	int = (ik - (Q2/4)*i0m - i2m + q0*i1m
		   + i2p + (Q2/4)*i0p + q0*i1p)
	(Q2/q^3)*(ec^2/(2*pi^2))*int
end

# Real part of transverse photon polarization tensor for degenerate electron gas.
# For internal integral definitions, see below.
function ΠTr_e_degen(q0,q,mue;mf=me)
	Q2 = q0^2 - q^2

	pf = sqrt(mue^2 - mf^2)

	ik = -(1 + Q2/(2*q^2))*q*(-(mue*pf) + mf^2*log((mue + pf)/mf))

	i0m = integ0m(q0,q,mue,pf,mf,Q2)
	i1m = integ1m(q0,q,mue,pf,mf,Q2)
	i2m = integ2m(q0,q,mue,pf,mf,Q2)

	i0p = integ0p(q0,q,mue,pf,mf,Q2)
	i1p = integ1p(q0,q,mue,pf,mf,Q2)
	i2p = integ2p(q0,q,mue,pf,mf,Q2)

	int = (ik + (Q2/4 + Q2^2/(8*q^2) + mf^2/2)*(i0m - i0p)
		   + (Q2/(2*q^2))*(i2m-i2p) - q0*Q2*(i1m + i1p)/(2*q^2))
	(ec^2/(2*pi^2))*int/q
end

# Helper functions and internal integrals for real parts of polarization tensors in degenerate electron gas.

function logabs(x)
	log(abs(x))
end

function arcTanh(x)
	0.5*log(abs((1 + x)/(1 - x)))
end

function integ0m(q0,q,mue,pf,mf,Q2)
	QQ2 = q^2 + q0^2
	C1 = sqrt(Q2*(Q2 - 4*mf^2))
	((2*q*logabs((mue + pf)/mf) +
	 2*mue*logabs((-2*pf*q - 2*mue*q0 + Q2)/(2*pf*q - 2*mue*q0 + Q2)) +
	 ((4*mf^2*q*q0 - 2*q*q0*Q2 + C1*QQ2)*(logabs(mf) + logabs(-(C1*q) - 2*mue*Q2 + q0*Q2)))/
	 sqrt(Q2*(-2*q0*(C1*q + 2*mf^2*q0) + Q2*QQ2)) +
	 ((2*q*q0*(-2*mf^2 + Q2) + C1*QQ2)*(-logabs(mf) - logabs(C1*q - 2*mue*Q2 + q0*Q2)))/
	 sqrt(Q2*(2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2)) -
	 ((4*mf^2*q*q0 - 2*q*q0*Q2 + C1*QQ2)*logabs(C1*mue*q + 2*mf^2*Q2 - mue*q0*Q2 +
											 pf*sqrt(Q2*(-2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2))))/
	 sqrt(Q2*(-2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2)) +
	 ((-4*mf^2*q*q0 + 2*q*q0*Q2 + C1*QQ2)*logabs(C1*mue*q - 2*mf^2*Q2 + mue*q0*Q2 -
											  pf*sqrt(Q2*(2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2))))/
	 sqrt(Q2*(2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2)))/2)
end

function integ0p(q0,q,mue,pf,mf,Q2)
	QQ2 = q^2 + q0^2
	C1 = sqrt(Q2*(Q2 - 4*mf^2))

	((-2*q*logabs((mue + pf)/mf) +
  2*mue*logabs((2*pf*q + 2*mue*q0 + Q2)/(-2*pf*q + 2*mue*q0 + Q2)) +
  ((4*mf^2*q*q0 - 2*q*q0*Q2 - C1*QQ2)*(logabs(mf) + logabs(-(C1*q) - (2*mue + q0)*Q2)))/
   sqrt(Q2*(2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2)) -
  ((4*mf^2*q*q0 - 2*q*q0*Q2 + C1*QQ2)*(-logabs(mf) - logabs(C1*q - (2*mue + q0)*Q2)))/
   sqrt(Q2*(-2*q0*(C1*q + 2*mf^2*q0) + Q2*QQ2)) -
  ((4*mf^2*q*q0 - 2*q*q0*Q2 + C1*QQ2)*logabs(C1*mue*q - 2*mf^2*Q2 - mue*q0*Q2 -
      pf*sqrt(Q2*(-2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2))))/
   sqrt(Q2*(-2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2)) +
  ((-4*mf^2*q*q0 + 2*q*q0*Q2 + C1*QQ2)*logabs(C1*mue*q + 2*mf^2*Q2 + mue*q0*Q2 +
      pf*sqrt(Q2*(2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2))))/
  sqrt(Q2*(2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2)))/2)
end

function integ1m(q0,q,mue,pf,mf,Q2)
	QQ2 = q^2 + q0^2
	C1 = sqrt(Q2*(Q2 - 4*mf^2))
	((pf*q + (q*q0*(-2*mf^2 + Q2)*arcTanh(mue/pf))/Q2 -
		((8*q0*Q2*(-2*mf^2 + Q2)*(4*mf^2*q^2 + Q2^2) + (4*C1*q - 4*q0*Q2)*
		  (4*mf^2*(q^4 + q^2*q0^2 - 2*q0^4) + (q^2 + 3*q0^2)*Q2^2))*
		 sqrt(Q2*(-2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2))*
		 arcTanh((-(C1*mue*q) - 2*mf^2*Q2 + mue*q0*Q2)/(sqrt(-mf^2 + mue^2)*
														sqrt(Q2*(-2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2)))))/
		(C1*Q2*(-64*mf^2*Q2^2 + (4*C1*q - 4*q0*Q2)^2)) -
		((16*mf^4*q0*(q^3 - q*q0^2)^2 + Q2^3*(C1*(q^3 + 3*q*q0^2) + q0*(3*q^2 + q0^2)*Q2) -
		  4*mf^2*Q2^2*(C1*(q^3 + 2*q*q0^2) + q0*(4*q^2 + q0^2)*Q2))*
		 arcTanh((C1*mue*q - 2*mf^2*Q2 + mue*q0*Q2)/(sqrt(-mf^2 + mue^2)*
													 sqrt(Q2*(2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2)))))/
		(4*C1*Q2^2*sqrt(Q2*(2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2))) +
		mue^2*logabs((-2*pf*q - 2*mue*q0 + Q2)/(2*pf*q - 2*mue*q0 + Q2)))/2)
end

function integ1p(q0,q,mue,pf,mf,Q2)
	QQ2 = q^2 + q0^2
	C1 = sqrt(Q2*(Q2 - 4*mf^2))
	((-(pf*q) + (q*q0*(-2*mf^2 + Q2)*arcTanh(mue/pf))/Q2 +
  ((8*q0*Q2*(-2*mf^2 + Q2)*(4*mf^2*q^2 + Q2^2) - (-4*C1*q + 4*q0*Q2)*
      (4*mf^2*(q^4 + q^2*q0^2 - 2*q0^4) + (q^2 + 3*q0^2)*Q2^2))*
    sqrt(Q2*(-2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2))*
    arcTanh((C1*mue*q - (2*mf^2 + mue*q0)*Q2)/
      (pf*sqrt(Q2*(-2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2)))))/
   (C1*Q2*(-64*mf^2*Q2^2 + (4*C1*q - 4*q0*Q2)^2)) +
  ((16*mf^4*q0*(q^3 - q*q0^2)^2 + Q2^3*(C1*(q^3 + 3*q*q0^2) + q0*(3*q^2 + q0^2)*Q2) -
     4*mf^2*Q2^2*(C1*(q^3 + 2*q*q0^2) + q0*(4*q^2 + q0^2)*Q2))*
    arcTanh((-2*mf^2*Q2 - mue*(C1*q + q0*Q2))/
      (pf*sqrt(Q2*(2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2)))))/
   (4*C1*Q2^2*sqrt(Q2*(2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2))) +
  mue^2*logabs((2*pf*q + 2*mue*q0 + Q2)/(-2*pf*q + 2*mue*q0 + Q2)))/2)
end

function integ2m(q0,q,mue,pf,mf,Q2)
	QQ2 = q^2 + q0^2
	C1 = sqrt(Q2*(Q2 - 4*mf^2))
	((2*mue*pf*q + (4*pf*q*q0*(-2*mf^2 + Q2))/Q2 +
  2*mf^2*q*arcTanh(mue/pf) +
  (q*(4*mf^2*(q^4 + q^2*q0^2 - 2*q0^4) + (q^2 + 3*q0^2)*Q2^2)*
    arcTanh(mue/pf))/Q2^2 -
  (8*sqrt(Q2*(-2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2))*
    (Q2^2*(4*mf^2*q^2 + Q2^2)*(4*mf^2*(q^4 + q^2*q0^2 - 2*q0^4) +
       (q^2 + 3*q0^2)*Q2^2) + q0*(4*C1*q - 4*q0*Q2)*(4*mf^4*(q^3 - q*q0^2)^2 +
       mf^2*(5*q^4 - 2*q^2*q0^2 - 3*q0^4)*Q2^2 + Q2^4*QQ2))*
    arcTanh((-(C1*mue*q) - 2*mf^2*Q2 + mue*q0*Q2)/(sqrt(-mf^2 + mue^2)*
       sqrt(Q2*(-2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2)))))/
   (C1*Q2^3*(-64*mf^2*Q2^2 + (4*C1*q - 4*q0*Q2)^2)) +
  (sqrt(Q2*(2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2))*
    (Q2^2*(4*mf^2*q^2 + Q2^2)*(4*mf^2*(q^4 + q^2*q0^2 - 2*q0^4) +
       (q^2 + 3*q0^2)*Q2^2) - 4*q0*(C1*q + q0*Q2)*(4*mf^4*(q^3 - q*q0^2)^2 +
       mf^2*(5*q^4 - 2*q^2*q0^2 - 3*q0^4)*Q2^2 + Q2^4*QQ2))*
    arcTanh((C1*mue*q - 2*mf^2*Q2 + mue*q0*Q2)/(sqrt(-mf^2 + mue^2)*
       sqrt(Q2*(2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2)))))/
   (2*C1*Q2^3*(-4*mf^2*Q2^2 + (C1*q + q0*Q2)^2)) +
   4*mue^3*logabs((-2*pf*q - 2*mue*q0 + Q2)/(2*pf*q - 2*mue*q0 + Q2)))/12)
end

function integ2p(q0,q,mue,pf,mf,Q2)
	QQ2 = q^2 + q0^2
	C1 = sqrt(Q2*(Q2 - 4*mf^2))
	((-2*mue*pf*q + (4*pf*q*q0*(-2*mf^2 + Q2))/Q2 - 2*mf^2*q*arcTanh(mue/pf) -
  (q*(4*mf^2*(q^4 + q^2*q0^2 - 2*q0^4) + (q^2 + 3*q0^2)*Q2^2)*arcTanh(mue/pf))/Q2^2 -
  (sqrt(Q2*(-2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2))*
    (8*Q2^2*(4*mf^2*q^2 + Q2^2)*(4*mf^2*(q^4 + q^2*q0^2 - 2*q0^4) +
       (q^2 + 3*q0^2)*Q2^2) - 8*q0*(-4*C1*q + 4*q0*Q2)*(4*mf^4*(q^3 - q*q0^2)^2 +
       mf^2*(5*q^4 - 2*q^2*q0^2 - 3*q0^4)*Q2^2 + Q2^4*QQ2))*
    arcTanh((C1*mue*q - (2*mf^2 + mue*q0)*Q2)/
      (pf*sqrt(Q2*(-2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2)))))/
   (C1*Q2^3*(-64*mf^2*Q2^2 + (4*C1*q - 4*q0*Q2)^2)) +
  (sqrt(Q2*(2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2))*
    (8*Q2^2*(4*mf^2*q^2 + Q2^2)*(4*mf^2*(q^4 + q^2*q0^2 - 2*q0^4) +
       (q^2 + 3*q0^2)*Q2^2) - 32*q0*(C1*q + q0*Q2)*(4*mf^4*(q^3 - q*q0^2)^2 +
       mf^2*(5*q^4 - 2*q^2*q0^2 - 3*q0^4)*Q2^2 + Q2^4*QQ2))*
    arcTanh((-2*mf^2*Q2 - mue*(C1*q + q0*Q2))/
      (pf*sqrt(Q2*(2*C1*q*q0 - 4*mf^2*q0^2 + Q2*QQ2)))))/
   (16*C1*Q2^3*(-4*mf^2*Q2^2 + (C1*q + q0*Q2)^2)) +
  4*mue^3*logabs((2*pf*q + 2*mue*q0 + Q2)/(-2*pf*q + 2*mue*q0 + Q2)))/12)
end

max_stars_repo_name: wderocco/DarkScatter

max_issues_repo_name: wderocco/DarkScatter

max_forks_repo_name: wderocco/DarkScatter


text:
include("utils.jl")

@testset "fixed Rank" begin
    M = FixedRankMatrices(3,2,2)
    Mc = FixedRankMatrices(3,2,2,ℂ)
    x = SVDMPoint([1.0 0.0; 0.0 1.0; 0.0 0.0])
    v = UMVTVector([0. 0.; 0. 0.; 1. 1.], [1. 0.; 0. 1.],zeros(2,2))
    @test repr(M) == "FixedRankMatrices(3, 2, 2, ℝ)"
    @test repr(Mc) == "FixedRankMatrices(3, 2, 2, ℂ)"
    @test sprint(show, "text/plain", x) == """
    SVDMPoint{Array{Float64,2},Array{Float64,1},Array{Float64,2}}
    U factor:
     3×2 Array{Float64,2}:
      1.0  0.0
      0.0  1.0
      0.0  0.0
    singular values:
     2-element Array{Float64,1}:
      1.0
      1.0
    Vt factor:
     2×2 Array{Float64,2}:
      1.0  0.0
      0.0  1.0"""
    @test sprint(show, "text/plain", v) == """
    UMVTVector{Array{Float64,2},Array{Float64,2},Array{Float64,2}}
    U factor:
     3×2 Array{Float64,2}:
      0.0  0.0
      0.0  0.0
      1.0  1.0
    M factor:
     2×2 Array{Float64,2}:
      1.0  0.0
      0.0  1.0
    Vt factor:
     2×2 Array{Float64,2}:
      0.0  0.0
      0.0  0.0"""
    @test inner(M,x,v,v) == norm(M,x,v)^2
    @test x == SVDMPoint(x.U,x.S,x.Vt)
    @test v == UMVTVector(v.U, v.M, v.Vt)
    @testset "Fixed Rank Matrices – Basics" begin
        @test representation_size(M) == (3,2)
        @test representation_size(Mc) == (3,2)
        @test manifold_dimension(M) == 6
        @test manifold_dimension(Mc) == 12
        @test !is_manifold_point(M, SVDMPoint([1. 0.; 0. 0.],2))
        @test_throws DomainError is_manifold_point(M, SVDMPoint([1. 0.; 0. 0.],2), true)

        @test !is_tangent_vector(M, SVDMPoint( [1. 0.; 0. 1.; 0. 0.] ), UMVTVector( zeros(2,1), zeros(1,2), zeros(2,2) )  )
        @test !is_tangent_vector(M, SVDMPoint([1. 0.; 0. 0.],2),v)
        @test_throws DomainError is_tangent_vector(M, SVDMPoint([1. 0.; 0. 0.],2), v, true)
        @test !is_tangent_vector(M, x, UMVTVector(x.U, v.M, x.Vt,2))
        @test_throws DomainError is_tangent_vector(M, x, UMVTVector(x.U, v.M, x.Vt,2), true)
        @test !is_tangent_vector(M, x, UMVTVector(v.U, v.M, x.Vt,2))
        @test_throws DomainError is_tangent_vector(M, x, UMVTVector(v.U, v.M, x.Vt,2), true)

        @test is_manifold_point(M,x)
        @test is_tangent_vector(M,x,v)
    end
    types = [
       [ Matrix{Float64}, Vector{Float64}, Matrix{Float64} ],
    ]
    TEST_FLOAT32 && push!(types, [ Matrix{Float32}, Vector{Float32}, Matrix{Float32} ] )

    for T in types
        @testset "Type $T" begin
            y = retract(M, x, v, PolarRetraction())
            z = SVDMPoint( [1/sqrt(2) 1/sqrt(2); 1/sqrt(2) -1/sqrt(2); 0. 0.] )
            pts = []
            for p in [x,y,z]
                push!(pts, SVDMPoint( convert.( T , [ p.U, p.S, p.Vt ] )... )   )
            end
            for p in pts
                @test is_manifold_point(M,p)
            end
            @testset "SVD MPoint Basics" begin
                s = svd(x.U*Diagonal(x.S)*x.Vt)
                x2 = SVDMPoint(s)
                x3 = SVDMPoint(s.U, s.S, s.Vt)
                @test SVDMPoint(x.U, x.S, x.Vt) == x
                @test x.S == x2.S
                @test x.U == x2.U
                @test x.Vt == x2.Vt
                @test x == x2
                @test x2.U == x3.U
                @test x2.S == x3.S
                @test x2.Vt == x3.Vt
                @test x2 == x3
                y = SVDMPoint([1. 0.; 0. 0.; 0. 0.], 1)
                s2 = svd([1. 0.; 0. 0.; 0. 0.])
                y2 = SVDMPoint(s2,1)
                y3 = SVDMPoint(s2.U, s2.S, s2.Vt,1)
                @test y.S == y2.S
                @test y.U == y2.U
                @test y.Vt == y2.Vt
                @test y == y2
                @test y2.U == y3.U
                @test y2.S == y3.S
                @test y2.Vt == y3.Vt
                @test y2==y3

                @test is_manifold_point(M,x)
                xM = x.U*Diagonal(x.S)*x.Vt
                @test is_manifold_point(M,xM)
                @test !is_manifold_point(M,xM[1:2,:])
                @test_throws DomainError is_manifold_point(M,xM[1:2,:], true)
                @test_throws DomainError is_manifold_point(FixedRankMatrices(3,2,1), x, true)
                @test_throws DomainError is_manifold_point(FixedRankMatrices(3,2,1), xM, true)
                xF1 = SVDMPoint(2*x.U, x.S, x.Vt)
                @test !is_manifold_point(M,xF1)
                @test_throws DomainError is_manifold_point(M,xF1,true)
                xF2 = SVDMPoint(x.U, x.S, 2*x.Vt)
                @test !is_manifold_point(M,xF2)
                @test_throws DomainError is_manifold_point(M,xF2,true)
            end
            @testset "UMV TVector Basics" begin
                w = UMVTVector(v.U, 2*v.M, v.Vt)
                @test v+w == UMVTVector(2*v.U, 3*v.M, 2*v.Vt)
                @test v-w == UMVTVector(0*v.U, -v.M, 0*v.Vt)
                @test 2*v == UMVTVector(2*v.U, 2*v.M, 2*v.Vt)
                @test v*2 == UMVTVector(v.U*2, v.M*2, v.Vt*2)
                @test 2\v == UMVTVector(2\v.U, 2\v.M, 2\v.Vt)
                @test v/2 == UMVTVector(v.U/2, v.M/2, v.Vt/2)
                @test +v == v
                @test -v == UMVTVector(-v.U, -v.M, -v.Vt)
                w = UMVTVector(v.U, v.M, v.Vt)
                @test v == w
                w = allocate(v, number_eltype(v))
                zero_tangent_vector!(M,w,x)
                oneP = SVDMPoint(one(zeros(3, 3)), ones(2), one(zeros(2, 2)), 2)
                @test oneP == one(x)
                oneV = UMVTVector(one(zeros(3,3)),one(zeros(2,2)),one(zeros(2,2)),2)
                @test oneV == one(v)
            end
            test_manifold(
                M,
                pts,
                test_exp_log = false,
                default_inverse_retraction_method = nothing,
                test_injectivity_radius = false,
                default_retraction_method = PolarRetraction(),
                test_is_tangent = false,
                test_project_tangent = true,
                test_vector_transport = false,
                test_forward_diff = false,
                test_reverse_diff = false,
                test_vector_spaces = false,
                test_vee_hat = false,
                test_tangent_vector_broadcasting = false, #broadcast not so easy for 3 matrix type
                projection_atol_multiplier = 15,
                retraction_methods = [PolarRetraction()]
            )
        end
    end
end


text:
function advance!(f!::Function, x::Array, p, t::Number, tn::Number, h::Number)
  n = length(x)

  mem = fill(zero(eltype(x)), n * 5)
  k1 = view(mem, 1:n)
  k2 = view(mem, (n + 1):(2 * n))
  k3 = view(mem, (2 * n + 1):(3 * n))
  k4 = view(mem, (3 * n + 1):(4 * n))
  work = view(mem, (4 * n + 1):(5 * n))

  end_loop = false
  while !end_loop && t < tn
    if tn - t < h
      h = tn - t
      end_loop = true
    end

    for i = 1:n
      work[i] = x[i]
    end
    f!(k1, work, p, t)

    for i = 1:n
      work[i] = x[i] + 0.5 * h * k1[i]
    end
    f!(k2, work, p, t + 0.5 * h)

    for i = 1:n
      work[i] = x[i] + 0.5 * h * k2[i]
    end
    f!(k3, work, p, t + 0.5 * h)

    for i = 1:n
      work[i] = x[i] + h * k3[i]
    end
    f!(k4, work, p, t + h)

    for i = 1:n
      x[i] += h / 6.0 * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i])
    end

    t += h
  end

  return
end


text:
using DataFrames

df = readtable(STDIN)

io = IOBuffer()
print(io, df)


text:
module DensestCommonSubgraph

# package code goes here
include("bitstr.jl")
include("dcs_graphs.jl")
include("dcs.jl")
include("dcs_parallel.jl")
include("dcs_lp.jl")
include("dcs_common.jl")

end # module


text:
"""

```
random_move!(x_ref::AbstractVector{T},
             irefatom::Int,
             sides::T,
             x_new::AbstractVector{T}, RNG) where T
```

Function that generates a new random position for a molecule.

The new position is returned in `x`, a previously allocated array.

`x_solvent_random` might be a view of the array that contains all the solvent molecules.

"""
function random_move!(x_ref::AbstractVector{T},
                      irefatom::Int,
                      sides::T,
                      x_new::AbstractVector{T}, RNG) where T

  # To avoid boundary problems, the center of coordinates are generated in a
  # much larger region, and wrapped aftwerwards
  scale = 100.

  # Generate random coordiantes for the center of mass
  ncm1 = -scale*sides[1]/2 + random(RNG,Float64)*scale*sides[1]
  ncm2 = -scale*sides[2]/2 + random(RNG,Float64)*scale*sides[2]
  ncm3 = -scale*sides[3]/2 + random(RNG,Float64)*scale*sides[3]
  newcm = T(ncm1,ncm2,ncm3)

  # Generate random rotation angles
  beta = (2*π)*random(RNG,Float64)
  gamma = (2*π)*random(RNG,Float64)
  theta = (2*π)*random(RNG,Float64)

  # Copy the coordinates of the molecule chosen to the random-coordinates vector
  @. x_new = x_ref

  # Take care that this molecule is not split by periodic boundary conditions, by
  # wrapping its coordinates around its reference atom
  wrap!(x_new,sides,x_ref[irefatom])

  # Move molecule to new position
  move!(x_new,newcm,beta,gamma,theta)

  nothing
end


text:
using Bio

function extract_seqs(seq_IDs, filenames, format)
	output_array = []
	for filename in filenames
		print("Checking "*filename)
		for entry in open(filename, format)
			print(record.id)
			if record.id in seq_IDs
				push!(output_array, record)
			end
		end
	end
	return output_array
end

function extract_nt_seqs(models,project,filenames,output;format=Bio.Seq.FASTQ)
	seqs_to_retrieve = []
	for model in models
		append!(seqs_to_retrieve, project["seq_info_table"][project["seq_info_table"][:model] .== model, :seqID])
	end
	output_sequences = extract_seqs(seqs_to_retrieve, filenames, format)
	output =open(output, format)
	write(output, output_sequences)
end

function extract_aa_seqs(models,project,filenames,output;format=Bio.Seq.FASTA)
	nt_seqs_to_retrieve = []
	translations = []
	aa_seqs_to_retrieve = []
	for model in models
		append!(nt_seqs_to_retrieve, project["seq_info_table"][project["seq_info_table"][:model] .== model, :seqID])
		append!(translations, project["seq_info_table"][project["seq_info_table"][:model] .== model, :translation])
		for (i, seqID) in enumerate(nt_seqs_to_retrieve)
			push!(aa_seqs_to_retrieve, join([seqID,translations[i]],'.'))
		end
	end
	output_sequences = extract_seqs(aa_seqs_to_retrieve, filenames, format)
	output =open(output, format)
	write(output, output_sequences)
end

max_stars_repo_name: ianpgm/apfamet

max_issues_repo_name: ianpgm/apfamet

max_forks_repo_name: ianpgm/apfamet


text:
abstract type AbstractGravity <: AbstractTerm end

struct Gravity <: AbstractGravity end
struct Buoyancy{𝒯} <: AbstractGravity end
struct FluctuationGravity <: AbstractGravity end

@inline calc_component!(source, ::Nothing, state, _...) = nothing

@inline function calc_component!(source, ::Gravity, state, aux, physics)
    ρ = state.ρ
    ∇Φ = aux.∇Φ

    source.ρu -= ρ * ∇Φ

    nothing
end

@inline function calc_component!(source, ::Buoyancy{(:ρ, :ρu, :ρθ)}, state, aux, physics)
    ρθ = state.ρθ
    α = physics.parameters.α
    g = physics.parameters.g
    orientation = physics.orientation

    k = vertical_unit_vector(orientation, aux)

    source.ρu -= -α * g * k * ρθ

    nothing
end

# FluctuationGravity Components
@inline calc_fluctuation_component!(source, _...) = nothing
@inline calc_component!(source, ::FluctuationGravity, _...) = nothing

@inline function calc_fluctuation_component!(source, ::FluctuationGravity, state_1, state_2, aux_1, aux_2)
    ρ_1, ρ_2 = state_1.ρ, state_2.ρ
    Φ_1, Φ_2 = aux_1.Φ, aux_2.Φ
    α = ave(ρ_1, ρ_2) * 0.5
    source.ρu -= α * (Φ_1 - Φ_2) * I

    nothing
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
module ProfileSVG

using FlameGraphs, Colors, UUIDs, Profile
using Base.StackTraces: StackFrame
const FlameGraph = FlameGraphs.Node{FlameGraphs.NodeData}

export @profview

"""
    @profview f(args...)

Clear the Profile buffer, profile `f(args...)`, and view the result graphically.
"""
macro profview(ex)
    return quote
        Profile.clear()
        @profile $(esc(ex))
        view()
    end
end

struct FGConfig
    g::FlameGraph
    fcolor
    fontsize::Int
end

include("svgwriter.jl")

"""
    ProfileSVG.view([fcolor], data=Profile.fetch(); fontsize=12, kwargs...)

View profiling results. See [FlameGraphs](https://github.com/timholy/FlameGraphs.jl)
for options for `kwargs` and the default value of `fcolor`.
"""
function view(fcolor, data::Vector{UInt64}=Profile.fetch(); fontsize::Integer=12, kwargs...)
    g = flamegraph(data; kwargs...)
    FGConfig(g, fcolor, fontsize)
end
function view(data::Vector{UInt64}=Profile.fetch(); fontsize::Integer=12, kwargs...)
    view(FlameGraphs.default_colors, data; fontsize=fontsize, kwargs...)
end

"""
    ProfileSVG.save([fcolor], io, g=flamegraph(); kwargs...)

Save profile results as an SVG file. See [FlameGraphs](https://github.com/timholy/FlameGraphs.jl)
for options for `kwargs` and the default value of `fcolor`.
"""
function save(fcolor, io::IO, g::FlameGraph; fontsize::Integer=12)
    show(io, MIME("image/svg+xml"), FGConfig(g, fcolor, fontsize))
end
function save(fcolor, filename::AbstractString, g::FlameGraph; kwargs...)
    open(filename, "w") do file
        save(fcolor, file, g; kwargs...)
    end
    return nothing
end
function save(fcolor, io::IO; fontsize::Integer=12, kwargs...)
    g = flamegraph(; kwargs...)
    save(fcolor, io, g; fontsize=fontsize)
end
function save(fcolor, filename::AbstractString; kwargs...)
    open(filename, "w") do file
        save(fcolor, file; kwargs...)
    end
    return nothing
end
save(io::IO, args...; kwargs...) = save(FlameGraphs.default_colors, io, args...; kwargs...)
save(filename::AbstractString, args...; kwargs...) = save(FlameGraphs.default_colors, filename, args...; kwargs...)

Base.showable(::MIME"image/svg+xml", fg::FGConfig) = true

function Base.show(io::IO, ::MIME"image/svg+xml", fg::FGConfig)
    g, fcolor, fontsize = fg.g, fg.fcolor, fg.fontsize
    ncols, nrows = length(g.data.span), FlameGraphs.depth(g)
    leftmargin = rightmargin = 10
    width = 1000
    topmargin = 30
    botmargin = 40
    rowheight = 15
    height = ceil(rowheight*nrows + botmargin + topmargin)
    xstep = (width - (leftmargin + rightmargin)) / ncols
    ystep = (height - (topmargin + botmargin)) / nrows
    avgcharwidth = 6  # for Verdana 12 pt font

    function printrec(io::IO, samples, xstart, xend, y, sf::StackFrame, rgb, fontsize)
        width = xend - xstart
        info = "$(sf.func) in $(sf.file):$(sf.line)"
        shortinfo = "$(sf.func) in $(basename(string(sf.file))):$(sf.line)"
        info = eschtml(info)
        shortinfo = eschtml(shortinfo)
        #if avgcharwidth*3 > width
        #    shortinfo = ""
        #elseif length(shortinfo) * avgcharwidth > width
        #    nchars = int(width/avgcharwidth)-2
        #    shortinfo = eschtml(info[1:nchars] * "..")
        #end
        r = round(Integer, 255*red(rgb))
        g = round(Integer, 255*green(rgb))
        b = round(Integer, 255*blue(rgb))
        print(io, """<rect vector-effect="non-scaling-stroke" x="$xstart" y="$y" width="$width" height="$ystep" fill="rgb($r,$g,$b)" rx="2" ry="2" data-shortinfo="$shortinfo" data-info="$info"/>\n""")
        #if shortinfo != ""
        println(io, """\n<text text-anchor="" x="$(xstart+4)" y="$(y+11.5)" font-size="$fontsize" font-family="Verdana" fill="rgb(0,0,0)" ></text>""")
        # end
    end

    function eschtml(str)
        s = replace(str, '<' => "<")
        s = replace(s, '>' => ">")
        s = replace(s, '&' => "&")
        s
    end

    function flamerects(fcolor, io::IO, g, j, nextidx)
        ndata = g.data
        thiscolor = fcolor(nextidx, j, ndata)
        y = height - j*ystep - botmargin
        xstart = (first(ndata.span)-1) * xstep + leftmargin
        xend   = ( last(ndata.span)-1) * xstep + leftmargin
        printrec(io, length(ndata.span), xstart, xend, y, ndata.sf, thiscolor, fontsize)

        for c in g
            flamerects(fcolor, io, c, j+1, nextidx)
        end
        return nothing
    end

    fig_id = string("fig-", replace(string(uuid4()), "-" => ""))
    svgheader(io, fig_id, width=width, height=height)

    nextidx = fill(1, nrows)
    flamerects(fcolor, io, g, 1, nextidx)

    svgfinish(io, fig_id)
end

end # module


text:
module RandomSeed

using Test
using TreeParzen
using Random

SEED = 12355

objective(params) = params[:x]^2
space = Dict(:x => HP.Uniform(:x, -10.0, 10.0))
nsamples = 10

# Figure out what the best parameter is with this seed.
#
# Note: this is not hard-coded on purpose since setting a random seed only guarantes you get the random state if
# and only if all other extraneous variables (i.e. bug fixes and speed improvements may change results) did not change.
Random.seed!(SEED)
best = fmin(objective, space, nsamples)
best_x = best[:x]

# Ensure equivalent calls to fmin, returns the same value.
for i in 1:10
    Random.seed!(SEED)
    best = fmin(objective, space, nsamples)
    @test isequal(best[:x], best_x)
end

end


text:
eigvals(A::Operator,n::Integer;tolerance::Float64=100eps()) =
    eigs(A,n;tolerance=tolerance)[1]

function eigs(A::Operator,n::Integer;tolerance::Float64=100eps())
    typ = eltype(A)

    ds=domainspace(A)
    C = Conversion(ds,rangespace(A))

    A1,C1 = zeros(typ,n,n),zeros(typ,n,n)
    A1[1:end,1:n] = A[1:n,1:n]
    C1[1:end,1:n] = C[1:n,1:n]

    λ,V=eigen(A1,C1)

    pruneeigs(λ,V,ds,tolerance)
end

eigvals(Bcs::Operator,A::Operator,n::Integer;tolerance::Float64=100eps()) =
    eigs(Bcs,A,n;tolerance=tolerance)[1]

function eigs(Bcs_in::Operator,A_in::Operator,n::Integer;tolerance::Float64=100eps())
    Bcs, A = promotedomainspace([Bcs_in, A_in])

    nf = size(Bcs,1)
    @assert isfinite(nf)

    typ = promote_type(eltype(Bcs),eltype(A))

    ds=domainspace(A)
    C = Conversion(ds,rangespace(A))

    A1,C1 = zeros(typ,n,n),zeros(typ,n,n)
    A1[1:nf,1:n] = Bcs[1:nf,1:n]
    A1[nf+1:end,1:n] = A[1:n-nf,1:n]
    C1[nf+1:end,1:n] = C[1:n-nf,1:n]

    λ,V = eigen(A1,C1)

    λ, V = pruneeigs(λ,V,ds,tolerance)
    p = sortperm(λ; lt=(x,y) -> isless(abs(x),abs(y)))
    λ[p], V[p]
end

function pruneeigs(λ,V,ds,tolerance)
    retλ=eltype(λ)[]
    retV=VFun{typeof(ds),eltype(V)}[]
    n=length(λ)
    for k=1:n
        if slnorm(V,n-3:n,k)≤tolerance
            push!(retλ,λ[k])
            push!(retV,Fun(ds,V[:,k]))
        end
    end
    retλ,retV
end


text:
module PolyhedraTests

using Test

using Gridap
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Arrays
using STLCutters

using STLCutters: Polyhedron
using STLCutters: STL
using STLCutters: restrict
using STLCutters: clip
using STLCutters: split
using STLCutters: surface, volume
using STLCutters: get_facet_planes
using STLCutters: compute_distances!
using STLCutters: plot
using STLCutters: check_graph
using STLCutters: compute_stl_model

p = Polyhedron(TRI)
@test check_graph(p)
#@test volume(p) ≈ 1/2
#@test surface(p) ≈ 2+√2

p = Polyhedron(QUAD)
@test check_graph(p)
#@test volume(p) ≈ 1
#@test surface(p) ≈ 4

p = Polyhedron(TET)
@test check_graph(p)
@test volume(p) ≈ 1/6
@test surface(p) ≈ 3/2 + (√3)/2

p = Polyhedron(HEX)
@test check_graph(p)
@test volume(p) ≈ 1
@test surface(p) ≈ 6

plot(p)

f, = writevtk(p,"p")
rm.(f)

vertices = [
  Point(0.1,-0.2,0.5),
  Point(1.0,-0.3,0.5),
  Point(-0.2,1.2,0.5)]

facet_to_vertices =
[[1,2,3]]

stlmodel = compute_stl_model( Table(facet_to_vertices), vertices )
stl = STL(stlmodel)

Πf = get_facet_planes(stl)
Π = [get_offset(stl,2)+1]
compute_distances!(p,Πf,Π)

p⁻,p⁺ = split(p,Π[1])

pn = [p⁻,p⁺]

@test volume(pn) ≈ 2*volume(p⁻) ≈ 1
@test surface(pn) ≈ 2*surface(p⁻) ≈ 8
@test surface(p⁻,stl) ≈ 1

end # module


text:
import Base.Iterators.Flatten

function puzzle19(path::String=joinpath(@__DIR__, "..", "data/19.txt"))
    data = readlines(path)
    n = length(data)
    m = Array{Any,2}(nothing, n, n)
    x = split.(data, "")
    for i in 1:n
        m[i, :] = x[i]
    end

    i = 1
    j = findall(m[1, :] .== "|")[1]
    path = []
    found = []
    direction = "down"
    while true
        # Keep track of the path
        append!(path, [(i, j)])
        # Continue along the direction
        if ismatch(r"-|\||[A-Z]", m[i, j])
            # Keep track of the letters we find
            if ismatch(r"[A-Z]", m[i, j])
                append!(found, m[i, j])
            end
            if direction=="down"
                i += 1
            elseif direction=="up"
                i -= 1
            elseif direction=="right"
                j += 1
            elseif direction=="left"
                j -= 1
            else
                error("invalid direction")
            end
        # Change direction
        elseif m[i, j] == "+"
            # There are no + on the outside edge
            if in(direction, ["down", "up"])
                # check other directions
                if ismatch(r"-|\|", m[i, j + 1])
                    j += 1
                    direction = "right"
                elseif ismatch(r"-|\|", m[i, j - 1])
                    j -= 1
                    direction = "left"
                else
                    print("wtf")
                end
            elseif in(direction, ["right", "left"])
                # check other directions
                if ismatch(r"-|\|", m[i + 1, j])
                    i += 1
                    direction = "down"
                elseif ismatch(r"-|\|", m[i - 1, j])
                    i -= 1
                    direction = "up"
                else
                    print("wtf")
                end
            else
                error("couldn't change direction!")
            end
        else
            # Must be at the end
            break
        end
    end
    partone = reduce(*, string.(found))
    parttwo = length(path) - 1

    return [partone, parttwo]
end


text:
### Data types

import Base.+, Base.-, Base.*, Base./, Base.==, Base.size, LinearAlgebra.norm, LinearAlgebra.dot, Statistics.mean, Statistics.var
export ModelGeom, Model, DataParams, SrcParams, Data, ABC_Geom, FunctionMod2Dat, LinFunctionMod2Mod

## Model geometry

mutable struct ModelGeom
	nz::Int64
	dz::Float32
end

function ==(geom1::ModelGeom, geom2::ModelGeom)
	return (geom1.nz == geom2.nz) && (geom1.dz == geom2.dz)
end

## Model

mutable struct Model
	geom::ModelGeom
	m::Array{Float32, 2}
end

function +(m1::Model, m2::Model)
	if m1.geom == m2.geom
		return Model(m1.geom, m1.m.+m2.m)
	end
end

function +(m1::Model, m2::Array{Float32})
	return Model(m1.geom, m1.m.+m2)
end

function +(m1::Model, c::Float32)
	return Model(m1.geom, m1.m.+c)
end

function -(m1::Model, m2::Model)
	if m1.geom == m2.geom
		return Model(m1.geom, m1.m.-m2.m)
	end
end

function /(m::Model, c::Float32)
	return Model(m.geom, m.m/c)
end

function /(m1::Model, m2::Array{Float32})
	return Model(m1.geom, m1.m./m2)
end

function dot(m1::Model, m2::Model)
	if m1.geom == m2.geom
		return vec(sum(m1.m.*m2.m, dims = 1))
	end
end

function dot(m1::Model, m2::Array{Float32, 2})
	return vec(sum(m1.m.*m2, dims = 1))
end

function dot(m1::Array{Float32, 2}, m2::Model)
	return vec(sum(m1.*m2.m, dims = 1))
end

function size(m::Model)
	return size(m.m)
end

function norm(m::Model)
	return vec(sqrt.(sum(m.m.^2f0, dims = 1)))
end

function mean(m::Model)
	return Model(m.geom, mean(m.m, dims = 2))
end

function var(m::Model)
	return Model(m.geom, var(m.m, dims = 2))
end

## Data parameters

mutable struct DataParams
	nt::Int64
	dt::Float32
	p::Array{Float32, 1}
end

function ==(dat_pars1::DataParams, dat_pars2::DataParams)
	return (dat_pars1.nt == dat_pars2.nt) && (dat_pars1.dt == dat_pars2.dt) && (dat_pars1.p == dat_pars2.p)
end

function size(dat_pars::DataParams)
	return dat_pars.nt, length(dat_pars.p)
end

## Source parameters

mutable struct SrcParams
	wav::Array{Float32, 1}
	dt::Float32
	zsrc::Float32
end

## Data

mutable struct Data
	d::Array{Float32, 3}
	pars::DataParams
	dttu::Union{Nothing, Array{Float32, 4}}
end

function +(d1::Float32, d2::Data)
	return Data(d1.+d2.d, d2.pars, nothing)
end

function +(d1::Data, d2::Data)
	d1.pars == d2.pars && (return Data(d1.d.+d2.d, d1.pars, nothing))
end

function +(d1::Data, d2::Array{Float32, 3})
	return Data(d1.d.+d2, d1.pars, nothing)
end

function -(d1::Array{Float32, 3}, d2::Data)
	return Data(d1.-d2.d, d2.pars, nothing)
end

function -(d1::Data, d2::Array{Float32, 3})
	return Data(d1.d.-d2, d1.pars, nothing)
end

function -(d1::Data, d2::Data)
	return Data(d1.d.-d2.d, d1.pars, nothing)
end

function *(d1::Float32, d2::Data)
	return Data(d1*d2.d, d2.pars, d2.dttu)
end

function norm(d::Data)
	return vec(sqrt.(sum(d.d.^2f0, dims = (1, 2))))
end

## Absorbing boundary condition options

mutable struct ABC_Geom
	nz::Union{Nothing, Int64}
	size::NTuple{2, Int64}
	fact::Float32
end

function ABC_Geom(; nz::Union{Nothing, Int64} = nothing, size::NTuple{2, Int64} = (40, 40), fact::Float32 = 1f0)
	return ABC_Geom(nz, size, fact)
end

## Model-to-data functions

mutable struct FunctionMod2Dat
	eval::Function
	jac_adj::Function
	back::Function
end

## Model-to-model linear functions

mutable struct LinFunctionMod2Mod
	eval::Function
	eval_adj::Function
end


text:
module FileIO

using DelimitedFiles, MatrixMarket, MAT, CSV, DataFrames, HDF5
export readtxt,   # Text file
	   readmtx,   # MatrixMarket file
	   readmat,   # MATLAB mat file
	   readhdf,   # HDF5 file
       readgenelist

function readmtx(filename)
    # read MatrixMarket file
    X=mmread(filename);
end

function readcsv(filename)
    # read CSV file with header and row name
    df = CSV.File(filename; datarow=2) |> DataFrame!
    X=convert(Matrix, df[:,2:end])
    genelist=df[:,1]
    return X,genelist
end

function readmat(filename)
    # read Matlab Mat file
    file=matopen(filename)
    X=read(file,"X")
    genelist=read(file,"genelist")
    close(file)
    X = convert(Array{Float64,2}, X)
    return X,genelist
end

function readtxt(filename)
    # read DLM text file
    X=readdlm(filename,',',Int16)
end

function readhdf(filename)
    # https://anndata.readthedocs.io/en/latest/anndata.AnnData.html
end

function readgenelist(filename::String,colidx::Integer=1)
    # read genelist
    genelist=readdlm(filename,'\t',String)
    genelist=vec(genelist[:,colidx])
end

# using UnicodePlots
# function showx(X)
#    spy(X)
# end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "FFMPEG"
version = v"4.1.0"

# Collection of sources required to build FFMPEG
sources = [
    "https://ffmpeg.org/releases/ffmpeg-4.1.tar.bz2" =>
    "b684fb43244a5c4caae652af9022ed5d85ce15210835bce054a33fb26033a1a5",
]

# Bash recipe for building across all platforms
# TODO: Theora and Opus once their releases are available
script = raw"""
cd $WORKSPACE/srcdir
cd ffmpeg-4.1/
sed -i 's/-lflite"/-lflite -lasound"/' configure
apk add coreutils yasm

if [[ "${target}" == *-linux-* ]]; then
    export ccOS="linux"
elif [[ "${target}" == *-apple-* ]]; then
    export ccOS="darwin"
elif [[ "${target}" == *-w32-* ]]; then
    export ccOS="mingw32"
elif [[ "${target}" == *-w64-* ]]; then
    export ccOS="mingw64"
elif [[ "${target}" == *-unknown-freebsd* ]]; then
    export ccOS="freebsd"
else
    export ccOS="linux"
fi

if [[ "${target}" == x86_64-* ]]; then
    export ccARCH="x86_64"
elif [[ "${target}" == i686-* ]]; then
    export ccARCH="i686"
elif [[ "${target}" == arm-* ]]; then
    export ccARCH="arm"
elif [[ "${target}" == aarch64-* ]]; then
    export ccARCH="aarch64"
elif [[ "${target}" == powerpc64le-* ]]; then
    export ccARCH="powerpc64le"
else
    export ccARCH="x86_64"
fi

export PKG_CONFIG_PATH="${prefix}/lib/pkgconfig"
pkg-config --list-all

./configure            \
  --enable-cross-compile \
  --cross-prefix=/opt/${target}/bin/${target}- \
  --arch=${ccARCH}     \
  --target-os=${ccOS}  \
  --cc="${CC}"         \
  --cxx="${CXX}"       \
  --dep-cc="${CC}"     \
  --ar="${AR}"         \
  --nm="${NM}"         \
  --objcc="${CC} -ObjC" \
  --sysinclude=${prefix}/include \
  --pkg-config=$(which pkg-config) \
  --pkg-config-flags=--static \
  --prefix=$prefix     \
  --sysroot=/opt/${target}/${target}/sys-root \
  --extra-libs=-lpthread \
  --enable-gpl         \
  --enable-version3    \
  --enable-nonfree     \
  --disable-static     \
  --enable-shared      \
  --enable-pic         \
  --disable-debug      \
  --disable-doc        \
  --enable-avresample  \
  --enable-libass      \
  --enable-libfdk-aac  \
  --enable-libfreetype \
  --enable-libmp3lame  \
  --enable-libvorbis   \
  --enable-libx264     \
  --enable-libx265     \
  --enable-libvpx      \
  --enable-encoders    \
  --enable-decoders    \
  --enable-muxers      \
  --enable-demuxers    \
  --enable-parsers     \
  --enable-openssl     \
  --disable-schannel   \
  --extra-cflags="-I${prefix}/include" \
  --extra-ldflags="-L${prefix}/lib"
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
        # glibc Linuces
        Linux(:i686),
        Linux(:x86_64),
        Linux(:aarch64),
        Linux(:armv7l),
        Linux(:powerpc64le),

        # musl Linuces
        Linux(:i686, :musl),
        Linux(:x86_64, :musl),
        Linux(:aarch64, :musl),
        Linux(:armv7l, :musl),

        # BSDs
        MacOS(:x86_64),
        FreeBSD(:x86_64),

        # Windows
        Windows(:i686),
        Windows(:x86_64),
]

# The products that we will ensure are always built
products(prefix) = [
    ExecutableProduct(prefix, "ffmpeg", :ffmpeg),
    ExecutableProduct(prefix, "ffprobe", :ffprobe),
    ExecutableProduct(prefix, "x264", :x264),
    ExecutableProduct(prefix, "x265", :x265),
    LibraryProduct(prefix, "libavformat", :libavformat),
    LibraryProduct(prefix, "libavcodec", :libavcodec),
    LibraryProduct(prefix, "libavutil", :libavutil),
    LibraryProduct(prefix, "libpostproc", :libpostproc),
    LibraryProduct(prefix, "libswresample", :libswresample),
    LibraryProduct(prefix, "libavdevice", :libavdevice),
    LibraryProduct(prefix, "libavresample", :libavresample),
    LibraryProduct(prefix, "libavfilter", :libavfilter),
    LibraryProduct(prefix, "libswscale", :libswscale)
]

# Dependencies that must be installed before this package can be built
# TODO: Theora and Opus once their releases are available
dependencies = [
    "https://github.com/JuliaIO/LibassBuilder/releases/download/v0.14.0-2/build_libass.v0.14.0.jl",
    "https://github.com/SimonDanisch/FDKBuilder/releases/download/0.1.6/build_libfdk.v0.1.6.jl",
    "https://github.com/SimonDanisch/FribidiBuilder/releases/download/0.14.0/build_fribidi.v0.14.0.jl",
    "https://github.com/JuliaGraphics/FreeTypeBuilder/releases/download/v2.9.1-4/build_FreeType2.v2.10.0.jl",
    "https://github.com/JuliaIO/LAMEBuilder/releases/download/v3.100.0-2/build_liblame.v3.100.0.jl",
    "https://github.com/JuliaIO/LibVorbisBuilder/releases/download/v1.3.6-2/build_libvorbis.v1.3.6.jl",
    "https://github.com/JuliaIO/OggBuilder/releases/download/v1.3.3-7/build_Ogg.v1.3.3.jl",
    "https://github.com/JuliaIO/LibVPXBuilder/releases/download/v1.8.0/build_LibVPX.v1.8.0.jl",
    "https://github.com/JuliaIO/x264Builder/releases/download/v2019.5.25-static/build_x264Builder.v2019.5.25.jl",
    "https://github.com/JuliaIO/x265Builder/releases/download/v3.0.0-static/build_x265Builder.v3.0.0.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/Bzip2-v1.0.6-0/build_Bzip2.v1.0.6.jl",
    "https://github.com/ianshmean/ZlibBuilder/releases/download/v1.2.11/build_Zlib.v1.2.11.jl",
    "https://github.com/giordano/Yggdrasil/releases/download/OpenSSL-v1.1.1%2B0/build_OpenSSL.v1.1.1+c.jl",
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Connectivity{PL,N}

A connectivity list of `N` indices representing a [`Polytope`](@ref)
of type `PL`. Indices are taken from a global vector of [`Point`](@ref).

Connectivity objects are constructed with the [`connect`](@ref) function.
"""
struct Connectivity{PL<:Polytope,N}
  indices::NTuple{N,Int}

  function Connectivity{PL,N}(indices) where {PL,N}
    @assert nvertices(PL) == N "invalid connectivity list"
    new(indices)
  end
end

"""
    paramdim(connectivity)

Return the parametric dimension of the `connectivity`.
"""
paramdim(::Type{Connectivity{PL,N}}) where {PL,N} = paramdim(PL)
paramdim(c::Connectivity) = paramdim(typeof(c))

"""
    issimplex(connectivity)

Tells whether or not the `connectivity` is simplex.
"""
issimplex(::Type{Connectivity{PL,N}}) where {PL,N} = issimplex(PL)
issimplex(c::Connectivity) = issimplex(typeof(c))

"""
    indices(connectivity)

Return the list of indices of the `connectivity`.
"""
indices(c::Connectivity) = c.indices

"""
    connect(indices, [PL])

Connect a list of `indices` from a global vector of [`Point`](@ref)
into a [`Polytope`](@ref) of type `PL`.

The type `PL` can be a [`Ngon`](@ref) in which case the length of
the indices is used to identify the actual N-gon type (e.g. Triangle).

Finally, the type `PL` can be ommitted. In this case, the indices are
assumed to be connected as a [`Ngon`](@ref) or as a [`Segment`](@ref).

## Examples

Connect indices into a Triangle:

```julia
connect((1,2,3), Triangle)
```

Connect indices into N-gons, a `Triangle` and a `Quadrangle`:

```julia
connect.([(1,2,3), (2,3,4,5)], Ngon)
```

Connect indices into N-gon or segment:

```julia
connect((1,2)) # Segment
connect((1,2,3)) # Triangle
connect((1,2,3,4)) # Quadrangle
```
"""
connect(indices::Tuple, PL::Type{<:Polytope}) =
  Connectivity{PL,length(indices)}(indices)

function connect(indices::Tuple, ::Type{Ngon})
  N = length(indices)
  Connectivity{Ngon{N},N}(indices)
end

function connect(indices::Tuple)
  N = length(indices)
  N > 2 ? connect(indices, Ngon) : connect(indices, Segment)
end

"""
    materialize(connec, points)

Materialize a face using the `connec` list and a global vector of `points`.
"""
function materialize(connec::Connectivity{PL},
                     points::AbstractVector{P}) where {PL<:Polytope,P<:Point}
  PL(view(points, SVector(connec.indices...)))
end

function Base.show(io::IO, c::Connectivity{PL}) where {PL}
  name = prettyname(PL)
  inds = c.indices
  print(io, "$name$inds")
end


text:
using Revise
using AdFem
using PyPlot
using LinearAlgebra
using Statistics

if pdata == 1
    mmesh = Mesh(50, 50, 1/50)
elseif pdata == 2
    mmesh = Mesh(joinpath(PDATA, "pipe.msh"))
    mmesh = Mesh(mmesh.nodes / 5, mmesh.elems)
elseif pdata == 3
    mmesh = Mesh(joinpath(PDATA, "twoholes_large.stl"))
    mmesh = Mesh(mmesh.nodes * 10, mmesh.elems)
end