text:
using LinearAlgebra
using JuMP
using Ipopt
using Parameters

## Setting parameters

@with_kw struct Par
    burn::Int64=0 ;
    nsim::Int64=500000 ;
    nfast::Int64=20000 ;
    nboost::Int64=20000;
end

Param=Par()

## Road map

# 1) Warm start
# 2) Chain filling
# 3) Optimization

## Warm start
# min ||g(x,e)|| st. Inequalities

function warm_start(p,q)
    T,N=size(q)
    wmst=Model(Ipopt.Optimizer)
    set_optimizer_attribute(wmst, "print_level", 0)

    # Declaring variables
    @variable(wmst, d[1:T,1:N] >= 0) # δ_it Convex cost function c'(q)
    @variable(wmst, 1.0e-16 <= qs[1:T,1:N] <= 1.0) # unobserved true q_it

    # Non-negative markups
    # @constraint(wmst,mup[t=1:T,j=1:N],d[t,j]<=p[t])
    for j=1:N
        for t=1:T
            set_upper_bound(d[t,j],p[t])
        end
    end

    # Common ratio constraings (p_t-δ_it)/q_it=(p_t-δ_jt)/q_jt
    @NLconstraint(wmst,eq[t=1:T,j=2:N],(p[t]-d[t,1])/qs[t,1]-(p[t]-d[t,j])/qs[t,j]==0)

    # Co-monotone constraint (δ_it'-δ_it)(q_it'-q_it)>=0
    for j=1:N
        @constraint(wmst,[t=1:T,s=1:T],(d[t,j]-d[s,j])*(qs[t,j]-qs[s,j])>=0)
    end

    # Setting objective function
    @NLexpression(wmst, obj, sum((p[t]*(qs[t,j] - q[t,j]))^2 for t=1:T, j=1:N))
    @NLobjective(wmst,Min,obj)

    ## Solving the Model
    JuMP.optimize!(wmst)

    # Checking convergence
    # JuMP.primal_status(wmst)==MOI.FEASIBLE_POINT || error("Unfeasible point in warm start")
    status=JuMP.primal_status(wmst)==MOI.FEASIBLE_POINT

    # Saving variable values
    dsim=value.(d)
    qsim=value.(qs)

    # Memory cleaning
    wmst=nothing
    println("Warm star done.")
    return dsim, qsim, status
end

## Build chain

# Getting random direction vector from uniform distribution on the unit sphere

# Compute interval Α to draw uniformly to get α≥0
# Interval Α has to satisfy all constraints
# 1) Sign constraints
# 2) Markup constraints
# 3) Common ratio property
# 4) Co-monotone property

## jump function

function jump(p, q, eδ, eq; n::Int64=0)
    n < 5 || return q.-eq, eδ, eq
    T, N = size(eδ)
    # Getting random direction vector from uniform distribution on the unit sphere
    v=rand(2,T,N)
    v=v.*2 .-1 # comes from unif*(1-(-1))+(-1)=unif*(ub-lb)+lb
    ξ=v./norm(v) # ξ( δ ,t, i) ξ(q*,t,i)
    unif=rand(N)

    # Initializing matrices
    αmax=ones(N).*1.0e6
    αmin=ones(N).*-1.0e6
    α=zeros(N)

    # Sign constraints for δit
    for i=1:N, t=1:T
        if ξ[1,t,i]<0
            αmax[i]=αmax[i] < -(eδ[t,i]/ξ[1,t,i]) ? αmax[i] : -(eδ[t,i]/ξ[1,t,i])
            αmin[i]=αmin[i] > (p[t]-eδ[t,i]/ξ[1,t,i]) ? αmin[i] : (p[t]-eδ[t,i]/ξ[1,t,i])
        else
            αmax[i]=αmax[i] < (p[t]-eδ[t,i]/ξ[1,t,i]) ? αmax[i] : (p[t]-eδ[t,i]/ξ[1,t,i])
            αmin[i]=αmin[i] > -(eδ[t,i]/ξ[1,t,i]) ? αmin[i] : -(eδ[t,i]/ξ[1,t,i])
        end
    end

    # Sign constraints for q^*it
    for i=1:N, t=1:T
        if ξ[2,t,i]<0
            αmax[i]=αmax[i] < -(eq[t,i]/ξ[2,t,i]) ? αmax[i] : -(eq[t,i]/ξ[2,t,i])
        else
            αmin[i]=αmin[i] > -(eq[t,i]/ξ[2,t,i]) ? αmin[i] : -(eq[t,i]/ξ[2,t,i])
        end
    end

    # Co-monotone constraints
    for i=1:N, t=1:T, s=1:T
        t==s && continue
        num=eq[s,i]-eq[t,i]
        den=ξ[2,s,i]-ξ[2,t,i]

        if  den > 0
            # αmin[i]=αmin[i] > -(num/den) ? αmin[i] : -(num/den)
            if num > 0
                αmin[i]=αmin[i] > -(num/den) ? αmin[i] : -(num/den)
            else
                αmax[i]=αmax[i] < -(num/den) ? αmax[i] : -(num/den) # If negative
            end
        else
            # αmax[i]=αmax[i] < -(num/den) ? αmax[i] : -(num/den)
            if num < 0
                αmin[i]=αmin[i] > -(num/den) ? αmin[i] : -(num/den) # if negative
            else
                αmax[i]=αmax[i] < -(num/den) ? αmax[i] : -(num/den)
            end
        end

        num=eδ[s,i]-eδ[t,i]
        den=ξ[1,s,i]-ξ[1,t,i]

        if  den > 0
            # αmin[i]=αmin[i] > -(num/den) ? αmin[i] : -(num/den)
            if num > 0
                αmin[i]=αmin[i] > -(num/den) ? αmin[i] : -(num/den)
            else
                αmax[i]=αmax[i] < -(num/den) ? αmax[i] : -(num/den) # if negative
            end
        else
            # αmax[i]=αmax[i] < -(num/den) ? αmax[i] : -(num/den)
            if num < 0
                αmin[i]=αmin[i] > -(num/den) ? αmin[i] : -(num/den) #if negative
            else
                αmax[i]=αmax[i] < -(num/den) ? αmax[i] : -(num/den)
            end
        end
    end

    # Draw uniformily from Α
    for i=1:N
        α[i]=αmax[i]>αmin[i] ? unif[i]*(αmax[i]-αmin[i])+αmin[i] : 0.0
    end

    # Generate random δt ∀i in t=0
    δsim=eδ.+α'.*ξ[1,:,:]

    # Common ratio constraint
    qtemp=eq.+α'.*ξ[2,:,:]
    # Pick a country i at random for every t
    pick=rand(collect(1:N),T)
    # estimate ratio for every t
    r=[(p[t]-δsim[t,i])/qtemp[t,i] for (t,i) ∈ zip(collect(1:T), pick)]
    # Use Common ratio property to estimate q^* for every i and t
    qsim=(p.-δsim)./r

    # Estimate W
    w=q.-qsim

    ## Check inequalities
    check=[(δsim[t,i]-δsim[s,i])*(qsim[t,i]-qsim[s,i]) for t=1:T, s=1:T, i=1:N]
    # minimum(check)>=0 || error("Co-monotone inequalities are not satisfied!")
    minimum(check)>=0 || jump(p,q,eδ,eq,n=n+1)

    return w, δsim, qsim
end

function jump_algo(p, q, eδ, eq; n::Int64=0)
    n < 5 || return q.-eq, eδ, eq
    T, N = size(eδ)
    # Getting random direction vector from uniform distribution on the unit sphere
    v=rand(2,T,N)
    v=v.*2 .-1 # comes from unif*(1-(-1))+(-1)=unif*(ub-lb)+lb
    ξ=v./norm(v) # ξ( δ ,t, i) ξ(q*,t,i)
    unif=rand(N)

    # Initializing matrices
    αmax=ones(N).*1.0e6
    αmin=ones(N).*-1.0e6
    βmax=ones(N).*1.0e6
    βmin=ones(N).*-1.0e6
    α=zeros(N)
    β=zeros(N)

    # Sign constraints for δit
    for i=1:N, t=1:T
        if ξ[1,t,i]<0
            αmax[i]=αmax[i] < -(eδ[t,i]/ξ[1,t,i]) ? αmax[i] : -(eδ[t,i]/ξ[1,t,i])
            αmin[i]=αmin[i] > (p[t]-eδ[t,i]/ξ[1,t,i]) ? αmin[i] : (p[t]-eδ[t,i]/ξ[1,t,i])
        else
            αmax[i]=αmax[i] < (p[t]-eδ[t,i]/ξ[1,t,i]) ? αmax[i] : (p[t]-eδ[t,i]/ξ[1,t,i])
            αmin[i]=αmin[i] > -(eδ[t,i]/ξ[1,t,i]) ? αmin[i] : -(eδ[t,i]/ξ[1,t,i])
        end
    end

    # Sign constraints for q^*it
    for i=1:N, t=1:T
        if ξ[2,t,i]<0
            βmax[i]=βmax[i] < -(eq[t,i]/ξ[2,t,i]) ? βmax[i] : -(eq[t,i]/ξ[2,t,i])
        else
            βmin[i]=βmin[i] > -(eq[t,i]/ξ[2,t,i]) ? βmin[i] : -(eq[t,i]/ξ[2,t,i])
        end
    end

    # Co-monotone constraints
    for i=1:N, t=1:T, s=1:T
        t==s && continue
        # For q
        num=eq[s,i]-eq[t,i]
        den=ξ[2,s,i]-ξ[2,t,i]

        if  den > 0
            # αmin[i]=αmin[i] > -(num/den) ? αmin[i] : -(num/den)
            if num > 0
                βmin[i]=βmin[i] > -(num/den) ? βmin[i] : -(num/den)
            else
                βmax[i]=βmax[i] < -(num/den) ? βmax[i] : -(num/den) # If negative
            end
        else
            # βmax[i]=βmax[i] < -(num/den) ? βmax[i] : -(num/den)
            if num < 0
                βmin[i]=βmin[i] > -(num/den) ? βmin[i] : -(num/den) # if negative
            else
                βmax[i]=βmax[i] < -(num/den) ? βmax[i] : -(num/den)
            end
        end

        num=eδ[s,i]-eδ[t,i]
        den=ξ[1,s,i]-ξ[1,t,i]

        if  den > 0
            # αmin[i]=αmin[i] > -(num/den) ? αmin[i] : -(num/den)
            if num > 0
                αmin[i]=αmin[i] > -(num/den) ? αmin[i] : -(num/den)
            else
                αmax[i]=αmax[i] < -(num/den) ? αmax[i] : -(num/den) # if negative
            end
        else
            # αmax[i]=αmax[i] < -(num/den) ? αmax[i] : -(num/den)
            if num < 0
                αmin[i]=αmin[i] > -(num/den) ? αmin[i] : -(num/den) #if negative
            else
                αmax[i]=αmax[i] < -(num/den) ? αmax[i] : -(num/den)
            end
        end
    end

    # Draw uniformily from Α
    for i=1:N
        α[i]=αmax[i]>αmin[i] ? unif[i]*(αmax[i]-αmin[i])+αmin[i] : 0.0
    end

    # Generate random δt ∀i in t=0
    δsim=eδ.+α'.*ξ[1,:,:]

    # Common ratio constraint
    # qsim=zeros(size(eq))
    # for t=1:T, s=1:T
    #     t==s && continue
    #     pick = rand(1:N) # pick a country at random
    #
    #     if δsim[s,pick]-δsim[t,pick] >= 0
    #         # t>=2 & (minimum(qsim[s,:]-qsim[t,:]) >= 0) && continue
    #         diff = 1.0-eq[t,pick]
    #         qs= diff > 0 ? rand()*(diff)+eq[t,pick] : eq[t,pick]
    #         r=(p[s]-δsim[s,pick])/qs
    #         qsim[s,:]=(p[s].-δsim[s,:])./r
    #     else
    #         # t>=2 & (maximum(qsim[s,:]-qsim[t,:]) <= 0) && continue
    #         diff = eq[t,pick]-1e-8
    #         qs = diff > 0 ? rand()*diff+0.0 : eq[t,pick]
    #         r=(p[s]-δsim[s,pick])/qs
    #         qsim[s,:]=(p[s].-δsim[s,:])./r
    #     end
    # end

    # Draw uniformily from B
    for i=1:N
        β[i]=βmax[i]>βmin[i] ? rand()*(βmax[i]-βmin[i])+βmin[i] : 0.0
    end

    # Generate random q^*it ∀i in t=0
    qtemp=eq.+β'.*ξ[2,:,:]
    # Pick a country i at random for every t
    pick=rand(collect(1:N),T)

    # estimate ratio for every t
    r=[(p[t]-δsim[t,i])/qtemp[t,i] for (t,i) ∈ zip(collect(1:T), pick)]
    # Use Common ratio property to estimate q^* for every i and t
    qsim=(p.-δsim)./r

    # Estimate W
    w=q.-qsim

    ## Check inequalities
    check=[(δsim[t,i]-δsim[s,i])*(qsim[t,i]-qsim[s,i]) for t=1:T, s=1:T, i=1:N]
    # minimum(check)>=0 || error("Co-monotone inequalities are not satisfied!")
    minimum(check)>=0 || jump(p,q,eδ,eq,n=n+1)

    return w, δsim, qsim
end

## Fill in MCMC Markov Chain Monte Carlo integration

function gchain(p, q; Par::Par=Param)
    @unpack burn, nsim = Par
    T, N = size(q)
    eδ, eq, status = warm_start(p,q)
    ew = eq-q
    mchain=zeros(N,(nsim-burn))
    status || return mchain, status # If Unfeasible solution in warmstart don't do loop
    for r=burn+1:nsim
        W, δsim, qsim = jump(p,q,eδ,eq)
        logtrydens = norm(p.*W)-norm(p.*ew)

        if logtrydens.>log.(rand())
            ew=W
            eδ=δsim
            eq=qsim
        end

        mchain[:,r-burn]=gM(p,ew)
    end
    println("Chain is ready.")
    return mchain, status
end

function gchain_algo(p, q; Par::Par=Param)
    @unpack burn, nsim = Par
    T, N = size(q)
    eδ, eq, status = warm_start(p,q)
    ew = eq-q
    mchain=zeros(N,(nsim-burn))
    status || return mchain, status # If Unfeasible solution in warmstart don't do loop
    for r=burn+1:nsim
        W, δsim, qsim = jump_algo(p,q,eδ,eq)
        logtrydens = norm(p.*W)-norm(p.*ew)

        if logtrydens.>log.(rand())
            ew=W
            eδ=δsim
            eq=qsim
        end

        mchain[:,r-burn]=gM(p,ew)
    end
    println("Chain is ready.")
    return mchain, status
end

function gchain_alt(p, q; Par::Par=Param)
    @unpack burn, nsim = Par
    T, N = size(q)
    eδ, eq, status = warm_start(p,q)
    ew = eq-q
    mchain=zeros(N,(nsim-burn))
    status || return mchain, status # If Unfeasible solution in warmstart don't do loop
    for r=burn+1:nsim
        W, δsim, qsim = jump(p,q,eδ,eq)
        logtrydens = norm(p.*W)-norm(p.*ew)

        if rand() < 1.4^(logtrydens[1])
            ew=W
            eδ=δsim
            eq=qsim
        end

        mchain[:,r-burn]=gM(p,ew)
    end
    println("Chain is ready.")
    return mchain, status
end

max_stars_repo_name: hans-mtz/SCR

max_issues_repo_name: hans-mtz/SCR

max_forks_repo_name: hans-mtz/SCR


text:
function namedtp_to_vec( x :: NamedTuple )
    l = length(x)
    t = [ as( ( typeof(xx) <: Array) ? Array : typeof(xx), size(xx) ) for xx in x]
    trans = as(NamedTuple{keys(x)}(t))
    return trans, inverse(trans, x)
end

function distribution_wrapper(pdist :: ArbitraryDistribution, trans_x :: TransformVariables.TransformTuple)

    f(x, params) = pdist.distribution( transform(trans_x, x), params )

    return f
end

function distribution_wrapper(pdist :: ArbitraryDistribution, trans_x :: TransformVariables.TransformTuple, trans_params :: TransformVariables.TransformTuple)

    f(x, params) = pdist.distribution( transform(trans_x, x), transform(trans_params, params) )

    return f
end

ext: jl


text:
using InteractiveDynamics
using GLMakie
using Agents
using Random, Statistics

Random.seed!(165) # hide
model, agent_step!, model_step! = Models.daisyworld(; solar_luminosity = 1.0, solar_change = 0.0, scenario = :change)
Daisy, Land = Agents.Models.Daisy, Agents.Models.Land

using AbstractPlotting: to_color
daisycolor(a::Daisy) = RGBAf0(to_color(a.breed))
landcolor = cgrad(:thermal)
daisycolor(a::Land) = to_color(landcolor[(a.temperature+50)/150])

daisyshape(a::Daisy) = '♣'
daisysize(a::Daisy) = 15
daisyshape(a::Land) = :rect
daisysize(a::Land) = 20
landfirst = by_type((Land, Daisy), false) # scheduler

# static daisyworld plot:
fig, abmstepper = abm_plot(
    model;
    ac = daisycolor, am = daisyshape, as = daisysize,
    scheduler = landfirst # crucial to change model scheduler!
)

display(fig)

# %% interactive evolution
fig, abmstepper = abm_play(
    model, agent_step!, model_step!;
    ac = daisycolor, am = daisyshape, as = daisysize,
    scheduler = landfirst # crucial to change model scheduler!
)

# %% daisyworld video:
model, agent_step!, model_step! = Models.daisyworld(; solar_luminosity = 1.0, solar_change = 0.0, scenario = :change)
abm_video(
    "daisyworld.mp4", model, agent_step!, model_step!;
    ac = daisycolor, am = daisyshape, as = daisysize,
    scheduler = landfirst, # crucial to change model scheduler!
    title = "Daisyworld"
)

# %% Parameter exploration and data collection:
params = Dict(
    :solar_change => -0.1:0.01:0.1,
    :surface_albedo => 0:0.01:1,
)

black(a) = a.breed == :black
white(a) = a.breed == :white
daisies(a) = a isa Daisy
land(a) = a isa Land
adata = [(black, count, daisies), (white, count, daisies), (:temperature, mean, land)]
mdata = [:solar_luminosity]

alabels = ["black", "white", "T"]
mlabels = ["L"]

model, agent_step!, model_step! = Models.daisyworld(; solar_luminosity = 1.0, solar_change = 0.0, scenario = :change)

fig, adf, mdf = abm_data_exploration(
    model, agent_step!, model_step!, params;
    ac = daisycolor, am = daisyshape, as = daisysize,
    mdata, adata, alabels, mlabels,
    scheduler = landfirst # crucial to change model scheduler!
)


text:
module Stiefel

# isometric W
# tangent vectors Δ = W*A + Z where Z = W⟂*B = Q*R = U*S*V
# with thus A' = -A, W'*Z = 0, W'*U = 0

using TensorKit
using TensorKit: similarstoragetype, SectorDict
using ..TensorKitManifolds: projecthermitian!, projectantihermitian!,
                            projectisometric!, projectcomplement!, PolarNewton,
                            _stiefelexp, _stiefellog, eleps
import ..TensorKitManifolds: base, checkbase,
                                inner, retract, transport, transport!

# special type to store tangent vectors using A and Z = Q*R,
struct StiefelTangent{T<:AbstractTensorMap, TA<:AbstractTensorMap}
    W::T
    A::TA
    Z::T
end
function Base.copy(Δ::StiefelTangent)
    Δ′ = StiefelTangent(Δ.W, copy(Δ.A), copy(Δ.Z))
    return Δ′
end

Base.getindex(Δ::StiefelTangent) = Δ.W * Δ.A + Δ.Z
base(Δ::StiefelTangent) = Δ.W
checkbase(Δ₁::StiefelTangent, Δ₂::StiefelTangent) = Δ₁.W == Δ₂.W ? Δ₁.W :
    throw(ArgumentError("tangent vectors with different base points"))

# Basic vector space behaviour
Base.:+(Δ₁::StiefelTangent, Δ₂::StiefelTangent) =
    StiefelTangent(checkbase(Δ₁,Δ₂), Δ₁.A + Δ₂.A, Δ₁.Z + Δ₂.Z)
Base.:-(Δ₁::StiefelTangent, Δ₂::StiefelTangent) =
    StiefelTangent(checkbase(Δ₁,Δ₂), Δ₁.A - Δ₂.A, Δ₁.Z - Δ₂.Z)
Base.:-(Δ::StiefelTangent) = (-1)*Δ

Base.:*(Δ::StiefelTangent, α::Real) = rmul!(copy(Δ), α)
Base.:*(α::Real, Δ::StiefelTangent) = lmul!(α, copy(Δ))
Base.:/(Δ::StiefelTangent, α::Real) = rmul!(copy(Δ), inv(α))
Base.:\(α::Real, Δ::StiefelTangent) = lmul!(inv(α), copy(Δ))

Base.zero(Δ::StiefelTangent) = StiefelTangent(Δ.W, zero(Δ.A), zero(Δ.Z))

function TensorKit.rmul!(Δ::StiefelTangent, α::Real)
    rmul!(Δ.A, α)
    rmul!(Δ.Z, α)
    return Δ
end
function TensorKit.lmul!(α::Real, Δ::StiefelTangent)
    lmul!(α, Δ.A)
    lmul!(α, Δ.Z)
    return Δ
end
function TensorKit.axpy!(α::Real, Δx::StiefelTangent, Δy::StiefelTangent)
    checkbase(Δx, Δy)
    axpy!(α, Δx.A, Δy.A)
    axpy!(α, Δx.Z, Δy.Z)
    return Δy
end
function TensorKit.axpby!(α::Real, Δx::StiefelTangent, β::Real, Δy::StiefelTangent)
    checkbase(Δx, Δy)
    axpby!(α, Δx.A, β, Δy.A)
    axpby!(α, Δx.Z, β, Δy.Z)
    return Δy
end

function TensorKit.dot(Δ₁::StiefelTangent, Δ₂::StiefelTangent)
    checkbase(Δ₁, Δ₂)
    dot(Δ₁.A, Δ₂.A) + dot(Δ₁.Z, Δ₂.Z)
end
TensorKit.norm(Δ::StiefelTangent, p::Real = 2) =
    norm((norm(Δ.A, p), norm(Δ.Z, p)), p)

# tangent space methods
function inner(W::AbstractTensorMap, Δ₁::StiefelTangent, Δ₂::StiefelTangent;
                metric = :euclidean)
    if metric == :euclidean
        return inner_euclidean(W, Δ₁, Δ₂)
    elseif metric == :canonical
        return inner_canonical(W, Δ₁, Δ₂)
    else
        throw(ArgumentError("unknown metric: $metric"))
    end
end
function project!(X::AbstractTensorMap, W::AbstractTensorMap; metric = :euclidean)
    if metric == :euclidean
        return project_euclidean!(X, W)
    elseif metric == :canonical
        return project_canonical!(X, W)
    else
        throw(ArgumentError("unknown metric: `metric = $metric`"))
    end
end
project(X::AbstractTensorMap, W::AbstractTensorMap; metric = :euclidean) =
    project!(copy(X), W; metric = metric)

function retract(W::AbstractTensorMap, Δ::StiefelTangent, α::Real; alg = :exp)
    if alg == :exp
        return retract_exp(W, Δ, α)
    elseif alg == :cayley
        return retract_cayley(W, Δ, α)
    else
        throw(ArgumentError("unknown algorithm: `alg = $alg`"))
    end
end
function invretract(Wold::AbstractTensorMap, Wnew::AbstractTensorMap; alg = :exp)
    if alg == :exp
        return invretract_exp(Wold, Wnew)
    elseif alg == :cayley
        return invretract_cayley(Wold, Wnew)
    else
        throw(ArgumentError("unknown algorithm: `alg = $alg`"))
    end
end
function transport!(Θ::StiefelTangent, W::AbstractTensorMap, Δ::StiefelTangent, α::Real, W′;
                    alg = :exp)
    if alg == :exp
        return transport_exp!(Θ, W, Δ, α, W′)
    elseif alg == :cayley
        return transport_cayley!(Θ, W, Δ, α, W′)
    else
        throw(ArgumentError("unknown algorithm: `alg = $alg`"))
    end
end
function transport(Θ::StiefelTangent, W::AbstractTensorMap, Δ::StiefelTangent, α::Real, W′;
                    alg = :exp)
    return transport!(copy(Θ), W, Δ, α, W′; alg = alg)
end

# euclidean metric
function inner_euclidean(W::AbstractTensorMap, Δ₁::StiefelTangent, Δ₂::StiefelTangent)
    Δ₁ === Δ₂ ? norm(Δ₁)^2 : real(dot(Δ₁,Δ₂))
end
function project_euclidean!(X::AbstractTensorMap, W::AbstractTensorMap)
    P = W'*X
    Z = mul!(X, W, P, -1, 1)
    A = projectantihermitian!(P)
    Z = projectcomplement!(Z, W)
    return StiefelTangent(W, A, Z)
end
project_euclidean(X, W) = project_euclidean!(copy(X), W)

# canonical metric
function inner_canonical(W::AbstractTensorMap, Δ₁::StiefelTangent, Δ₂::StiefelTangent)
    if Δ₁ === Δ₂
        return (norm(Δ₁.A)^2)/2 + norm(Δ₁.Z)^2
    else
        return real(dot(Δ₁.A, Δ₂.A)/2 + dot(Δ₁.Z, Δ₂.Z))
    end
end
function project_canonical!(X::AbstractTensorMap, W::AbstractTensorMap)
    P = W'*X
    Z = mul!(X, W, P, -1, 1)
    A = rmul!(projectantihermitian!(P), 2)
    Z = projectcomplement!(Z, W)
    return StiefelTangent(W, A, Z)
end
project_canonical(X, W) = project_canonical!(copy(X), W)

# geodesic retraction for canonical metric using exponential
function stiefelexp(W::AbstractTensorMap,
                    A::AbstractTensorMap,
                    Z::AbstractTensorMap,
                    α::Real)
    S = spacetype(W)
    G = sectortype(W)
    Wdata′ = TensorKit.SectorDict{G, storagetype(W)}()
    Qdata = TensorKit.SectorDict{G, storagetype(W)}()
    Qdata′ = TensorKit.SectorDict{G, storagetype(W)}()
    Rdata′ = TensorKit.SectorDict{G, storagetype(W)}()
    dims = TensorKit.SectorDict{G, Int}()
    for c in blocksectors(W)
        w′, q, q′, r′ = _stiefelexp(block(W, c), block(A, c), block(Z, c), α)
        Wdata′[c] = w′
        Qdata[c] = q
        Qdata′[c] = q′
        Rdata′[c] = r′
        dims[c] = size(q, 2)
    end
    V = S(dims)
    W′ = TensorMap(Wdata′, space(W))
    Q = TensorMap(Qdata, codomain(W)←V)
    Q′ = TensorMap(Qdata′, codomain(W)←V)
    R′ = TensorMap(Rdata′, V←domain(W))
    return W′, Q, Q′, R′
end

# can be computed efficiently: O(np^2) + O(p^3)
function retract_exp(W::AbstractTensorMap, Δ::StiefelTangent, α::Real)
    W == base(Δ) || throw(ArgumentError("not a valid tangent vector at base point"))
    W′, Q, Q′, R′ = stiefelexp(W, Δ.A, Δ.Z, α)
    A′ = Δ.A
    Z′ = projectcomplement!(Q′*R′, W′) # to ensure orthogonality
    return W′, StiefelTangent(W′, A′, Z′)
end
function invretract_exp(Wold::AbstractTensorMap, Wnew::AbstractTensorMap;
                            tol = eleps(Wold)^(2/3))
    space(Wold) == space(Wnew) || throw(SectorMismatch())

    S = spacetype(Wold)
    G = sectortype(Wold)
    Adata = TensorKit.SectorDict{G, storagetype(Wold)}()
    Zdata = TensorKit.SectorDict{G, storagetype(Wold)}()
    for c in blocksectors(Wold)
        a, q, r = _stiefellog(block(Wold, c), block(Wnew, c); tol = tol)
        Adata[c] = a
        Zdata[c] = q*r
    end
    A = TensorMap(Adata, domain(Wold)←domain(Wold))
    Z = TensorMap(Zdata, space(Wold))
    return StiefelTangent(Wold, A, Z)
end

# vector transport compatible with above `retract`: also differentiated retraction
# isometric for both euclidean and canonical metric
# not parallel transport for either metric as the corresponding connection has torsion
# can be computed efficiently: O(np^2) + O(p^3)
function transport_exp!(Θ::StiefelTangent, W::AbstractTensorMap,
                        Δ::StiefelTangent, α::Real, W′::AbstractTensorMap)
    W == checkbase(Δ,Θ) || throw(ArgumentError("not a valid tangent vector at base point"))
    # TODO: stiefelexp call does not depend on Θ
    # cache result or find some other way not to recompute this information
    _W′, Q, Q′, R′ = stiefelexp(W, Δ.A, Δ.Z, α)
    W′ ≈ _W′ || throw(ArgumentError("not a valid tangent vector at end point"))
    A = Θ.A
    Z = Θ.Z
    QZ = Q'*Θ.Z
    A′ = Θ.A
    Z′ = projectcomplement!(mul!(Z, (Q′-Q), QZ, 1, 1), W′)
    return StiefelTangent(W′, A′, Z′)
end
transport_exp(Θ::StiefelTangent, W::AbstractTensorMap, Δ::StiefelTangent, α::Real, W′) =
    transport_exp!(copy(Θ), W, Δ, α, W′)

# Cayley retraction, slightly more efficient than above?
# can be computed efficiently: O(np^2) + O(p^3)
function retract_cayley(W::AbstractTensorMap, Δ::StiefelTangent, α::Real)
    W == base(Δ) || throw(ArgumentError("not a valid tangent vector at base point"))
    A, Z = Δ.A, Δ.Z
    ZdZ = Z'*Z
    X = axpy!(α^2/4, ZdZ, axpy!(-α/2, A, one(A)))
    iX = inv(X)
    W′ = projectisometric!((2*W+α*Z)*iX - W; alg = PolarNewton())
    A′ = projectantihermitian!((A - (α/2)*ZdZ)*iX)
    Z′ = (Z-α*(W+α/2*Z)*(iX*ZdZ))
    Z′ = projectcomplement!(Z′*projecthermitian!(iX), W′)
    return W′, StiefelTangent(W′, A′, Z′)
end
function invretract_cayley(Wold::AbstractTensorMap, Wnew::AbstractTensorMap)
    space(Wnew) == space(Wold) || throw(SpaceMismatch())

    P = Wold'*Wnew
    iX = rmul!(axpy!(1,P,one(P)), 1/2)
    X = inv(iX)
    Z = projectcomplement!(Wnew - Wold*P, Wold)*X
    A = projectantihermitian!(rmul!(axpy!(-1, X, mul!(one(X), Z', Z, 1/4, 1)), 2))
    return StiefelTangent(Wold, A, Z)
end

# vector transport compatible with above `retract_caley`, but not differentiated retraction
# isometric for both euclidean and canonical metric
# can be computed efficiently: O(np^2) + O(p^3)
function transport_cayley!(Θ::StiefelTangent, W::AbstractTensorMap, Δ::StiefelTangent,
                            α::Real, W′)
    W == checkbase(Δ,Θ) || throw(ArgumentError("not a valid tangent vector at base point"))
    A, Z = Δ.A, Δ.Z
    X = axpy!(α^2/4, Z'*Z, axpy!(-α/2, A, one(A)))
    A′ = Θ.A
    ZdZ = Z'*Θ.Z
    Z′ = projectcomplement!(axpy!(-α, (W+(α/2)*Z)*(X\ZdZ), Θ.Z), W′)
    return StiefelTangent(W′, A′, Z′)
end
transport_cayley(Θ::StiefelTangent, W::AbstractTensorMap, Δ::StiefelTangent, α::Real, W′) =
    transport_cayley!(copy(Θ), W, Δ, α, W′)

end


text:
function goalstate(interpreter::ConcreteInterpreter,
                   domain::Domain, problem::GenericProblem)
    types = Set{Term}([Compound(ty, Term[o]) for (o, ty) in problem.objtypes])
    facts = Set{Term}(flatten_conjs(problem.goal))
    state = GenericState(types, facts, Dict{Symbol,Any}())
    return state
end

function goalstate(interpreter::ConcreteInterpreter,
                   domain::Domain, objtypes::AbstractDict, terms)
    types = Set{Term}([Compound(ty, Term[o]) for (o, ty) in objtypes])
    goal = GenericState(types, Set{Term}(), Dict{Symbol,Any}())
    for t in flatten_conjs(terms)
        if t.name == :(==) # Function equality
            @assert length(t.args) == 2 "Assignments must have two arguments."
            term, val = t.args[1], t.args[2]
            @assert(is_func(term, domain) && !is_attached_func(term, domain),
                    "Unrecognized function $(term.name).")
            @assert(is_ground(term), "Assigned terms must be ground.")
            @assert(isa(val, Const), "Terms must be equal to constants.")
            goal[term] = val.name
        elseif is_pred(t, domain) # Predicates
            push!(goal.facts, t)
        else
            error("Term $t in $terms cannot be handled.")
        end
    end
    return goal
end


text:
module ExpectationMaximization

using ArgCheck
using Distributions
using StatsFuns: logsumexp

# Extended functions
import Distributions: fit_mle

export fit_mle

include("fit_em.jl")
include("fit_mle_product_distributions.jl")
end


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

@test integrate(tan(c+d*x)^5*(a+a*sin(c+d*x)), x) == :(-1*a*d^-1*sin(c+d*x)+-1*a^2*d^-1*(a+-1*a*sin(c+d*x))^-1+-23//16*a*d^-1*log(1+-1*sin(c+d*x))+1//8*a^2*d^-1*(a+a*sin(c+d*x))^-1+1//8*a^3*d^-1*(a+-1*a*sin(c+d*x))^-2+7//16*a*d^-1*log(1+sin(c+d*x)))
@test integrate(tan(c+d*x)^3*(a+a*sin(c+d*x)), x) == :(a*d^-1*sin(c+d*x)+(1/2)*a^2*d^-1*(a+-1*a*sin(c+d*x))^-1+-1//4*a*d^-1*log(1+sin(c+d*x))+5//4*a*d^-1*log(1+-1*sin(c+d*x)))
@test integrate((a+a*sin(c+d*x))*tan(c+d*x), x) == :(-1*a*d^-1*log(1+-1*sin(c+d*x))+-1*a*d^-1*sin(c+d*x))
@test integrate((a+a*sin(c+d*x))*cot(c+d*x), x) == :(a*d^-1*log(sin(c+d*x))+a*d^-1*sin(c+d*x))
@test integrate(cot(c+d*x)^3*(a+a*sin(c+d*x)), x) == :(-1*a*d^-1*csc(c+d*x)+-1*a*d^-1*log(sin(c+d*x))+-1*a*d^-1*sin(c+d*x)+-1//2*a*d^-1*csc(c+d*x)^2)
@test integrate(cot(c+d*x)^5*(a+a*sin(c+d*x)), x) == :(a*d^-1*csc(c+d*x)^2+a*d^-1*log(sin(c+d*x))+a*d^-1*sin(c+d*x)+2*a*d^-1*csc(c+d*x)+-1//3*a*d^-1*csc(c+d*x)^3+-1//4*a*d^-1*csc(c+d*x)^4)
@test integrate(cot(c+d*x)^7*(a+a*sin(c+d*x)), x) == :(a*d^-1*csc(c+d*x)^3+-1*a*d^-1*log(sin(c+d*x))+-1*a*d^-1*sin(c+d*x)+-3*a*d^-1*csc(c+d*x)+-3//2*a*d^-1*csc(c+d*x)^2+-1//5*a*d^-1*csc(c+d*x)^5+-1//6*a*d^-1*csc(c+d*x)^6+3//4*a*d^-1*csc(c+d*x)^4)
@test integrate(tan(c+d*x)^6*(a+a*sin(c+d*x)), x) == :(-1*a*x+a*d^-1*cos(c+d*x)+a*d^-1*tan(c+d*x)+-1*a*d^-1*sec(c+d*x)^3+3*a*d^-1*sec(c+d*x)+-1//3*a*d^-1*tan(c+d*x)^3+1//5*a*d^-1*sec(c+d*x)^5+1//5*a*d^-1*tan(c+d*x)^5)
@test integrate(tan(c+d*x)^4*(a+a*sin(c+d*x)), x) == :(a*x+-1*a*d^-1*cos(c+d*x)+-1*a*d^-1*tan(c+d*x)+-2*a*d^-1*sec(c+d*x)+1//3*a*d^-1*sec(c+d*x)^3+1//3*a*d^-1*tan(c+d*x)^3)
@test integrate(tan(c+d*x)^2*(a+a*sin(c+d*x)), x) == :(-1*a*x+a*d^-1*cos(c+d*x)+a*d^-1*(1+-1*sin(c+d*x))^-1*cos(c+d*x))
@test integrate(cot(c+d*x)^2*(a+a*sin(c+d*x)), x) == :(-1*a*x+a*d^-1*cos(c+d*x)+-1*a*d^-1*arctanh(cos(c+d*x))+-1*a*d^-1*cot(c+d*x))
@test integrate(cot(c+d*x)^4*(a+a*sin(c+d*x)), x) == :(a*x+a*d^-1*cot(c+d*x)+-3//2*a*d^-1*cos(c+d*x)+-1//3*a*d^-1*cot(c+d*x)^3+3//2*a*d^-1*arctanh(cos(c+d*x))+-1//2*a*d^-1*cot(c+d*x)^2*cos(c+d*x))
@test integrate(cot(c+d*x)^6*(a+a*sin(c+d*x)), x) == :(-1*a*x+-1*a*d^-1*cot(c+d*x)+-15//8*a*d^-1*arctanh(cos(c+d*x))+-1//5*a*d^-1*cot(c+d*x)^5+1//3*a*d^-1*cot(c+d*x)^3+15//8*a*d^-1*cos(c+d*x)+-1//4*a*d^-1*cot(c+d*x)^4*cos(c+d*x)+5//8*a*d^-1*cot(c+d*x)^2*cos(c+d*x))
@test integrate((a+a*sin(c+d*x))^2*tan(c+d*x)^5, x) == :(-2*a^2*d^-1*sin(c+d*x)+-31//8*a^2*d^-1*log(1+-1*sin(c+d*x))+-9//4*a^3*d^-1*(a+-1*a*sin(c+d*x))^-1+-1//2*a^2*d^-1*sin(c+d*x)^2+-1//8*a^2*d^-1*log(1+sin(c+d*x))+1//4*a^4*d^-1*(a+-1*a*sin(c+d*x))^-2)
@test integrate((a+a*sin(c+d*x))^2*tan(c+d*x)^3, x) == :(a^3*d^-1*(a+-1*a*sin(c+d*x))^-1+(1/2)*a^2*d^-1*sin(c+d*x)^2+2*a^2*d^-1*sin(c+d*x)+3*a^2*d^-1*log(1+-1*sin(c+d*x)))
@test integrate((a+a*sin(c+d*x))^2*tan(c+d*x), x) == :(-2*a^2*d^-1*log(1+-1*sin(c+d*x))+-2*a^2*d^-1*sin(c+d*x)+-1//2*a^2*d^-1*sin(c+d*x)^2)
@test integrate((a+a*sin(c+d*x))^2*cot(c+d*x)^3, x) == :(-1//2*a^-2*d^-1*(a+a*sin(c+d*x))^4*csc(c+d*x)^2)
@test integrate((a+a*sin(c+d*x))^2*cot(c+d*x)^7, x) == :((1/2)*a^2*d^-1*csc(c+d*x)^4+-6*a^2*d^-1*csc(c+d*x)+-2*a^2*d^-1*sin(c+d*x)+2*a^2*d^-1*csc(c+d*x)^3+2*a^2*d^-1*log(sin(c+d*x))+-2//5*a^2*d^-1*csc(c+d*x)^5+-1//2*a^2*d^-1*sin(c+d*x)^2+-1//6*a^2*d^-1*csc(c+d*x)^6)
@test integrate((a+a*sin(c+d*x))^2*tan(c+d*x)^6, x) == :(-9//2*x*a^2+-2*a^2*d^-1*sec(c+d*x)^3+2*a^2*d^-1*cos(c+d*x)+6*a^2*d^-1*sec(c+d*x)+-3//2*a^2*d^-1*tan(c+d*x)^3+2//5*a^2*d^-1*sec(c+d*x)^5+9//2*a^2*d^-1*tan(c+d*x)+9//10*a^2*d^-1*tan(c+d*x)^5+-1//2*a^2*d^-1*sin(c+d*x)^2*tan(c+d*x)^5)
@test integrate((a+a*sin(c+d*x))^2*tan(c+d*x)^4, x) == :(7//2*x*a^2+-16//3*a^2*d^-1*cos(c+d*x)+-7//2*a^2*d^-1*cos(c+d*x)*sin(c+d*x)+-8//3*a^2*d^-1*(1+-1*sin(c+d*x))^-1*sin(c+d*x)^2*cos(c+d*x)+1//3*a^4*d^-1*(a+-1*a*sin(c+d*x))^-2*sin(c+d*x)^3*cos(c+d*x))
@test integrate((a+a*sin(c+d*x))^2*tan(c+d*x)^2, x) == :(-5//2*x*a^2+2*a^2*d^-1*cos(c+d*x)+(1/2)*a^2*d^-1*cos(c+d*x)*sin(c+d*x)+2*a^2*d^-1*(1+-1*sin(c+d*x))^-1*cos(c+d*x))
@test integrate((a+a*sin(c+d*x))^2, x) == :(3//2*x*a^2+-2*a^2*d^-1*cos(c+d*x)+-1//2*a^2*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate((a+a*sin(c+d*x))^2*cot(c+d*x)^2, x) == :(-1//2*x*a^2+-1*a^2*d^-1*cot(c+d*x)+-2*a^2*d^-1*arctanh(cos(c+d*x))+2*a^2*d^-1*cos(c+d*x)+(1/2)*a^2*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate((a+a*sin(c+d*x))^2*cot(c+d*x)^4, x) == :(-1//2*x*a^2+-2*a^2*d^-1*cos(c+d*x)+3*a^2*d^-1*arctanh(cos(c+d*x))+-1//3*a^2*d^-1*cot(c+d*x)^3+-1*a^2*d^-1*cot(c+d*x)*csc(c+d*x)+-1//2*a^2*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate((a+a*sin(c+d*x))^3*tan(c+d*x)^7, x) == :(7*a^3*d^-1*sin(c+d*x)+-13//8*a^5*d^-1*(a+-1*a*sin(c+d*x))^-2+-1//16*a^3*d^-1*log(1+sin(c+d*x))+1//3*a^3*d^-1*sin(c+d*x)^3+1//6*a^6*d^-1*(a+-1*a*sin(c+d*x))^-3+3//2*a^3*d^-1*sin(c+d*x)^2+71//8*a^4*d^-1*(a+-1*a*sin(c+d*x))^-1+209//16*a^3*d^-1*log(1+-1*sin(c+d*x)))
@test integrate((a+a*sin(c+d*x))^3*tan(c+d*x)^3, x) == :(2*a^4*d^-1*(a+-1*a*sin(c+d*x))^-1+5*a^3*d^-1*sin(c+d*x)+7*a^3*d^-1*log(1+-1*sin(c+d*x))+1//3*a^3*d^-1*sin(c+d*x)^3+3//2*a^3*d^-1*sin(c+d*x)^2)
@test integrate((a+a*sin(c+d*x))^3*tan(c+d*x), x) == :(-4*a^3*d^-1*log(1+-1*sin(c+d*x))+-4*a^3*d^-1*sin(c+d*x)+-3//2*a^3*d^-1*sin(c+d*x)^2+-1//3*a^3*d^-1*sin(c+d*x)^3)
@test integrate((a+a*sin(c+d*x))^3*cot(c+d*x)^3, x) == :(-3*a^3*d^-1*csc(c+d*x)+-2*a^3*d^-1*sin(c+d*x)+2*a^3*d^-1*log(sin(c+d*x))+-3//2*a^3*d^-1*sin(c+d*x)^2+-1//2*a^3*d^-1*csc(c+d*x)^2+-1//3*a^3*d^-1*sin(c+d*x)^3)
@test integrate((a+a*sin(c+d*x))^3*tan(c+d*x)^6, x) == :(-23//2*x*a^3+-136//15*a^3*d^-1*cos(c+d*x)^3+136//5*a^3*d^-1*cos(c+d*x)+23//2*a^3*d^-1*cos(c+d*x)*sin(c+d*x)+-13//15*a^5*d^-1*(a+-1*a*sin(c+d*x))^-2*sin(c+d*x)^4*cos(c+d*x)+1//5*a^6*d^-1*(a+-1*a*sin(c+d*x))^-3*sin(c+d*x)^5*cos(c+d*x)+23//3*a^6*d^-1*(a^3+-1*a^3*sin(c+d*x))^-1*sin(c+d*x)^3*cos(c+d*x))
@test integrate((a+a*sin(c+d*x))^3*tan(c+d*x)^4, x) == :(17//2*x*a^3+-6*a^3*d^-1*cos(c+d*x)+1//3*a^3*d^-1*cos(c+d*x)^3+-25//3*a^3*d^-1*(1+-1*sin(c+d*x))^-1*cos(c+d*x)+-3//2*a^3*d^-1*cos(c+d*x)*sin(c+d*x)+2//3*a^3*d^-1*(1+-1*sin(c+d*x))^-2*cos(c+d*x))
@test integrate((a+a*sin(c+d*x))^3*tan(c+d*x)^2, x) == :(-11//2*x*a^3+5*a^3*d^-1*cos(c+d*x)+-1//3*a^3*d^-1*cos(c+d*x)^3+4*a^3*d^-1*(1+-1*sin(c+d*x))^-1*cos(c+d*x)+3//2*a^3*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate((a+a*sin(c+d*x))^3, x) == :(5//2*x*a^3+-4*a^3*d^-1*cos(c+d*x)+1//3*a^3*d^-1*cos(c+d*x)^3+-3//2*a^3*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate((a+a*sin(c+d*x))^3*cot(c+d*x)^2, x) == :((1/2)*x*a^3+-1*a^3*d^-1*cot(c+d*x)+-3*a^3*d^-1*arctanh(cos(c+d*x))+3*a^3*d^-1*cos(c+d*x)+-1//3*a^3*d^-1*cos(c+d*x)^3+3//2*a^3*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate((a+a*sin(c+d*x))^4*tan(c+d*x)^5, x) == :(a^6*d^-1*(a+-1*a*sin(c+d*x))^-2+-25*a^4*d^-1*log(1+-1*sin(c+d*x))+-16*a^4*d^-1*sin(c+d*x)+-11*a^5*d^-1*(a+-1*a*sin(c+d*x))^-1+-9//2*a^4*d^-1*sin(c+d*x)^2+-4//3*a^4*d^-1*sin(c+d*x)^3+-1//4*a^4*d^-1*sin(c+d*x)^4)
@test integrate((a+a*sin(c+d*x))^4*tan(c+d*x)^3, x) == :(4*a^4*d^-1*sin(c+d*x)^2+4*a^5*d^-1*(a+-1*a*sin(c+d*x))^-1+12*a^4*d^-1*sin(c+d*x)+16*a^4*d^-1*log(1+-1*sin(c+d*x))+1//4*a^4*d^-1*sin(c+d*x)^4+4//3*a^4*d^-1*sin(c+d*x)^3)
@test integrate((a+a*sin(c+d*x))^4*tan(c+d*x), x) == :(-8*a^4*d^-1*log(1+-1*sin(c+d*x))+-8*a^4*d^-1*sin(c+d*x)+-7//2*a^4*d^-1*sin(c+d*x)^2+-4//3*a^4*d^-1*sin(c+d*x)^3+-1//4*a^4*d^-1*sin(c+d*x)^4)
@test integrate((a+a*sin(c+d*x))^4*cot(c+d*x)^3, x) == :(-4*a^4*d^-1*csc(c+d*x)+5*a^4*d^-1*log(sin(c+d*x))+-5//2*a^4*d^-1*sin(c+d*x)^2+-4//3*a^4*d^-1*sin(c+d*x)^3+-1//2*a^4*d^-1*csc(c+d*x)^2+-1//4*a^4*d^-1*sin(c+d*x)^4)
@test integrate((a+a*sin(c+d*x))^4*tan(c+d*x)^4, x) == :(163//8*x*a^4+-16*a^4*d^-1*cos(c+d*x)+4//3*a^4*d^-1*cos(c+d*x)^3+-56//3*a^4*d^-1*(1+-1*sin(c+d*x))^-1*cos(c+d*x)+-35//8*a^4*d^-1*cos(c+d*x)*sin(c+d*x)+-1//4*a^4*d^-1*sin(c+d*x)^3*cos(c+d*x)+4//3*a^4*d^-1*(1+-1*sin(c+d*x))^-2*cos(c+d*x))
@test integrate((a+a*sin(c+d*x))^4*tan(c+d*x)^2, x) == :(-95//8*x*a^4+12*a^4*d^-1*cos(c+d*x)+-4//3*a^4*d^-1*cos(c+d*x)^3+8*a^4*d^-1*(1+-1*sin(c+d*x))^-1*cos(c+d*x)+1//4*a^4*d^-1*sin(c+d*x)^3*cos(c+d*x)+31//8*a^4*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate((a+a*sin(c+d*x))^4, x) == :(35//8*x*a^4+-8*a^4*d^-1*cos(c+d*x)+4//3*a^4*d^-1*cos(c+d*x)^3+-27//8*a^4*d^-1*cos(c+d*x)*sin(c+d*x)+-1//4*a^4*d^-1*sin(c+d*x)^3*cos(c+d*x))
@test integrate((a+a*sin(c+d*x))^4*cot(c+d*x)^2, x) == :(17//8*x*a^4+-1*a^4*d^-1*cot(c+d*x)+-4*a^4*d^-1*arctanh(cos(c+d*x))+4*a^4*d^-1*cos(c+d*x)+-4//3*a^4*d^-1*cos(c+d*x)^3+1//4*a^4*d^-1*sin(c+d*x)^3*cos(c+d*x)+23//8*a^4*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate((a+a*sin(c+d*x))^4*cot(c+d*x)^4, x) == :(-61//8*x*a^4+-5*a^4*d^-1*cot(c+d*x)+2*a^4*d^-1*arctanh(cos(c+d*x))+-1//3*a^4*d^-1*cot(c+d*x)^3+4//3*a^4*d^-1*cos(c+d*x)^3+-2*a^4*d^-1*cot(c+d*x)*csc(c+d*x)+-19//8*a^4*d^-1*cos(c+d*x)*sin(c+d*x)+-1//4*a^4*d^-1*sin(c+d*x)^3*cos(c+d*x))
@test integrate((a+a*sin(c+d*x))^4*cot(c+d*x)^6, x) == :(97//8*x*a^4+-4*a^4*d^-1*cos(c+d*x)+10*a^4*d^-1*cot(c+d*x)+-5//3*a^4*d^-1*cot(c+d*x)^3+-4//3*a^4*d^-1*cos(c+d*x)^3+-1//5*a^4*d^-1*cot(c+d*x)^5+5//2*a^4*d^-1*arctanh(cos(c+d*x))+-1*a^4*d^-1*csc(c+d*x)^3*cot(c+d*x)+1//4*a^4*d^-1*sin(c+d*x)^3*cos(c+d*x)+5//2*a^4*d^-1*cot(c+d*x)*csc(c+d*x)+15//8*a^4*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate((a+a*sin(c+d*x))^-1*tan(c+d*x)^7, x) == :(-35//128*a^-1*d^-1*arctanh(sin(c+d*x))+1//8*a^-1*d^-1*tan(c+d*x)^8+-35//192*a^-1*d^-1*tan(c+d*x)^3*sec(c+d*x)+-1//8*a^-1*d^-1*tan(c+d*x)^7*sec(c+d*x)+7//48*a^-1*d^-1*tan(c+d*x)^5*sec(c+d*x)+35//128*a^-1*d^-1*sec(c+d*x)*tan(c+d*x))
@test integrate((a+a*sin(c+d*x))^-1*tan(c+d*x)^5, x) == :(1//6*a^-1*d^-1*tan(c+d*x)^6+5//16*a^-1*d^-1*arctanh(sin(c+d*x))+-5//16*a^-1*d^-1*sec(c+d*x)*tan(c+d*x)+-1//6*a^-1*d^-1*tan(c+d*x)^5*sec(c+d*x)+5//24*a^-1*d^-1*tan(c+d*x)^3*sec(c+d*x))
@test integrate((a+a*sin(c+d*x))^-1*tan(c+d*x)^3, x) == :(-3//8*a^-1*d^-1*arctanh(sin(c+d*x))+1//4*a^-1*d^-1*tan(c+d*x)^4+-1//4*a^-1*d^-1*tan(c+d*x)^3*sec(c+d*x)+3//8*a^-1*d^-1*sec(c+d*x)*tan(c+d*x))
@test integrate((a+a*sin(c+d*x))^-1*cot(c+d*x), x) == :(a^-1*d^-1*log(sin(c+d*x))+-1*a^-1*d^-1*log(1+sin(c+d*x)))
@test integrate((a+a*sin(c+d*x))^-1*cot(c+d*x)^3, x) == :(a^-1*d^-1*csc(c+d*x)+-1//2*a^-1*d^-1*csc(c+d*x)^2)
@test integrate((a+a*sin(c+d*x))^-1*cot(c+d*x)^5, x) == :(-1*a^-1*d^-1*csc(c+d*x)+-1//4*a^-1*d^-1*cot(c+d*x)^4+1//3*a^-1*d^-1*csc(c+d*x)^3)
@test integrate((a+a*sin(c+d*x))^-1*cot(c+d*x)^7, x) == :(a^-1*d^-1*csc(c+d*x)+-2//3*a^-1*d^-1*csc(c+d*x)^3+-1//6*a^-1*d^-1*cot(c+d*x)^6+1//5*a^-1*d^-1*csc(c+d*x)^5)
@test integrate((a+a*sin(c+d*x))^-1*cot(c+d*x)^9, x) == :(a^-1*d^-1*csc(c+d*x)^3+-1*a^-1*d^-1*csc(c+d*x)+-3//5*a^-1*d^-1*csc(c+d*x)^5+-1//8*a^-1*d^-1*cot(c+d*x)^8+1//7*a^-1*d^-1*csc(c+d*x)^7)
@test integrate((a+a*sin(c+d*x))^-1*tan(c+d*x)^6, x) == :(a^-1*d^-1*sec(c+d*x)+-1*a^-1*d^-1*sec(c+d*x)^3+-1//7*a^-1*d^-1*sec(c+d*x)^7+1//7*a^-1*d^-1*tan(c+d*x)^7+3//5*a^-1*d^-1*sec(c+d*x)^5)
@test integrate((a+a*sin(c+d*x))^-1*tan(c+d*x)^4, x) == :(-1*a^-1*d^-1*sec(c+d*x)+-1//5*a^-1*d^-1*sec(c+d*x)^5+1//5*a^-1*d^-1*tan(c+d*x)^5+2//3*a^-1*d^-1*sec(c+d*x)^3)
@test integrate((a+a*sin(c+d*x))^-1*tan(c+d*x)^2, x) == :(a^-1*d^-1*sec(c+d*x)+-1//3*a^-1*d^-1*sec(c+d*x)^3+1//3*a^-1*d^-1*tan(c+d*x)^3)
@test integrate((a+a*sin(c+d*x))^-1, x) == :(-1*d^-1*(a+a*sin(c+d*x))^-1*cos(c+d*x))
@test integrate((a+a*sin(c+d*x))^-1*cot(c+d*x)^2, x) == :(a^-1*d^-1*arctanh(cos(c+d*x))+-1*a^-1*d^-1*cot(c+d*x))
@test integrate((a+a*sin(c+d*x))^-1*cot(c+d*x)^4, x) == :(-1//2*a^-1*d^-1*arctanh(cos(c+d*x))+-1//3*a^-1*d^-1*cot(c+d*x)^3+(1/2)*a^-1*d^-1*cot(c+d*x)*csc(c+d*x))
@test integrate((a+a*sin(c+d*x))^-1*cot(c+d*x)^6, x) == :(-1//5*a^-1*d^-1*cot(c+d*x)^5+3//8*a^-1*d^-1*arctanh(cos(c+d*x))+-3//8*a^-1*d^-1*cot(c+d*x)*csc(c+d*x)+1//4*a^-1*d^-1*cot(c+d*x)^3*csc(c+d*x))
@test integrate((a+a*sin(c+d*x))^-1*cot(c+d*x)^8, x) == :(-5//16*a^-1*d^-1*arctanh(cos(c+d*x))+-1//7*a^-1*d^-1*cot(c+d*x)^7+-5//24*a^-1*d^-1*cot(c+d*x)^3*csc(c+d*x)+1//6*a^-1*d^-1*cot(c+d*x)^5*csc(c+d*x)+5//16*a^-1*d^-1*cot(c+d*x)*csc(c+d*x))
@test integrate((a+a*sin(c+d*x))^-2*tan(c+d*x)^7, x) == :(-1//4*d^-1*(a+a*sin(c+d*x))^-2+-1//32*d^-1*(a+-1*a*sin(c+d*x))^-2+21//256*d^-1*(a^2+-1*a^2*sin(c+d*x))^-1+35//256*d^-1*(a^2+a^2*sin(c+d*x))^-1+-7//128*a^-2*d^-1*arctanh(sin(c+d*x))+-5//64*a^2*d^-1*(a+a*sin(c+d*x))^-4+1//80*a^3*d^-1*(a+a*sin(c+d*x))^-5+1//192*a*d^-1*(a+-1*a*sin(c+d*x))^-3+19//96*a*d^-1*(a+a*sin(c+d*x))^-3)
@test integrate((a+a*sin(c+d*x))^-2*tan(c+d*x)^5, x) == :(-5//32*d^-1*(a^2+a^2*sin(c+d*x))^-1+-5//64*d^-1*(a^2+-1*a^2*sin(c+d*x))^-1+1//4*d^-1*(a+a*sin(c+d*x))^-2+1//64*d^-1*(a+-1*a*sin(c+d*x))^-2+-7//48*a*d^-1*(a+a*sin(c+d*x))^-3+1//32*a^2*d^-1*(a+a*sin(c+d*x))^-4+5//64*a^-2*d^-1*arctanh(sin(c+d*x)))
@test integrate((a+a*sin(c+d*x))^-2*tan(c+d*x)^3, x) == :(-1//4*d^-1*(a+a*sin(c+d*x))^-2+1//16*d^-1*(a^2+-1*a^2*sin(c+d*x))^-1+3//16*d^-1*(a^2+a^2*sin(c+d*x))^-1+-1//8*a^-2*d^-1*arctanh(sin(c+d*x))+1//12*a*d^-1*(a+a*sin(c+d*x))^-3)
@test integrate((a+a*sin(c+d*x))^-2*tan(c+d*x), x) == :(-1//4*d^-1*(a^2+a^2*sin(c+d*x))^-1+1//4*d^-1*(a+a*sin(c+d*x))^-2+1//4*a^-2*d^-1*arctanh(sin(c+d*x)))
@test integrate((a+a*sin(c+d*x))^-2*cot(c+d*x), x) == :(d^-1*(a^2+a^2*sin(c+d*x))^-1+a^-2*d^-1*log(sin(c+d*x))+-1*a^-2*d^-1*log(1+sin(c+d*x)))
@test integrate((a+a*sin(c+d*x))^-2*cot(c+d*x)^3, x) == :(-2*a^-2*d^-1*log(1+sin(c+d*x))+2*a^-2*d^-1*csc(c+d*x)+2*a^-2*d^-1*log(sin(c+d*x))+-1//2*a^-2*d^-1*csc(c+d*x)^2)
@test integrate((a+a*sin(c+d*x))^-2*cot(c+d*x)^5, x) == :(-1//2*a^-2*d^-1*csc(c+d*x)^2+-1//4*a^-2*d^-1*csc(c+d*x)^4+2//3*a^-2*d^-1*csc(c+d*x)^3)
@test integrate((a+a*sin(c+d*x))^-2*cot(c+d*x)^7, x) == :((1/2)*a^-2*d^-1*csc(c+d*x)^2+-2//3*a^-2*d^-1*csc(c+d*x)^3+-1//6*a^-2*d^-1*csc(c+d*x)^6+2//5*a^-2*d^-1*csc(c+d*x)^5)
@test integrate((a+a*sin(c+d*x))^-2*cot(c+d*x)^9, x) == :(-4//5*a^-2*d^-1*csc(c+d*x)^5+-1//2*a^-2*d^-1*csc(c+d*x)^2+-1//8*a^-2*d^-1*csc(c+d*x)^8+1//4*a^-2*d^-1*csc(c+d*x)^4+1//6*a^-2*d^-1*csc(c+d*x)^6+2//3*a^-2*d^-1*csc(c+d*x)^3+2//7*a^-2*d^-1*csc(c+d*x)^7)
@test integrate((a+a*sin(c+d*x))^-2*cot(c+d*x)^11, x) == :((1/2)*a^-2*d^-1*csc(c+d*x)^2+-6//7*a^-2*d^-1*csc(c+d*x)^7+-2//3*a^-2*d^-1*csc(c+d*x)^3+-1//2*a^-2*d^-1*csc(c+d*x)^4+-1//10*a^-2*d^-1*csc(c+d*x)^10+1//4*a^-2*d^-1*csc(c+d*x)^8+2//9*a^-2*d^-1*csc(c+d*x)^9+6//5*a^-2*d^-1*csc(c+d*x)^5)
@test integrate((a+a*sin(c+d*x))^-2*cot(c+d*x)^13, x) == :(-8//5*a^-2*d^-1*csc(c+d*x)^5+-8//9*a^-2*d^-1*csc(c+d*x)^9+-1//2*a^-2*d^-1*csc(c+d*x)^2+-1//3*a^-2*d^-1*csc(c+d*x)^6+-1//4*a^-2*d^-1*csc(c+d*x)^8+-1//12*a^-2*d^-1*csc(c+d*x)^12+2//3*a^-2*d^-1*csc(c+d*x)^3+2//11*a^-2*d^-1*csc(c+d*x)^11+3//4*a^-2*d^-1*csc(c+d*x)^4+3//10*a^-2*d^-1*csc(c+d*x)^10+12//7*a^-2*d^-1*csc(c+d*x)^7)
@test integrate((a+a*sin(c+d*x))^-3*tan(c+d*x)^5, x) == :(-5//128*d^-1*(a^3+a^3*sin(c+d*x))^-1+-1//32*d^-1*(a^3+-1*a^3*sin(c+d*x))^-1+1//6*d^-1*(a+a*sin(c+d*x))^-3+-7//64*a*d^-1*(a+a*sin(c+d*x))^-4+-5//64*a^-1*d^-1*(a+a*sin(c+d*x))^-2+1//40*a^2*d^-1*(a+a*sin(c+d*x))^-5+1//128*a^-1*d^-1*(a+-1*a*sin(c+d*x))^-2+1//128*a^-3*d^-1*arctanh(sin(c+d*x)))
@test integrate((a+a*sin(c+d*x))^-3*tan(c+d*x)^3, x) == :(-1//6*d^-1*(a+a*sin(c+d*x))^-3+1//16*d^-1*(a^3+a^3*sin(c+d*x))^-1+1//32*d^-1*(a^3+-1*a^3*sin(c+d*x))^-1+-1//32*a^-3*d^-1*arctanh(sin(c+d*x))+1//16*a*d^-1*(a+a*sin(c+d*x))^-4+3//32*a^-1*d^-1*(a+a*sin(c+d*x))^-2)
@test integrate((a+a*sin(c+d*x))^-3*tan(c+d*x), x) == :(-1//8*d^-1*(a^3+a^3*sin(c+d*x))^-1+1//6*d^-1*(a+a*sin(c+d*x))^-3+-1//8*a^-1*d^-1*(a+a*sin(c+d*x))^-2+1//8*a^-3*d^-1*arctanh(sin(c+d*x)))
@test integrate((a+a*sin(c+d*x))^-3*cot(c+d*x), x) == :(d^-1*(a^3+a^3*sin(c+d*x))^-1+a^-3*d^-1*log(sin(c+d*x))+(1/2)*a^-1*d^-1*(a+a*sin(c+d*x))^-2+-1*a^-3*d^-1*log(1+sin(c+d*x)))
@test integrate((a+a*sin(c+d*x))^-3*cot(c+d*x)^3, x) == :(2*d^-1*(a^3+a^3*sin(c+d*x))^-1+-5*a^-3*d^-1*log(1+sin(c+d*x))+3*a^-3*d^-1*csc(c+d*x)+5*a^-3*d^-1*log(sin(c+d*x))+-1//2*a^-3*d^-1*csc(c+d*x)^2)
@test integrate((a+a*sin(c+d*x))^-3*cot(c+d*x)^5, x) == :(a^-3*d^-1*csc(c+d*x)^3+-4*a^-3*d^-1*log(1+sin(c+d*x))+-2*a^-3*d^-1*csc(c+d*x)^2+4*a^-3*d^-1*csc(c+d*x)+4*a^-3*d^-1*log(sin(c+d*x))+-1//4*a^-3*d^-1*csc(c+d*x)^4)
@test integrate((a+a*sin(c+d*x))^-3*cot(c+d*x)^7, x) == :(-3//4*a^-3*d^-1*csc(c+d*x)^4+-1//6*a^-3*d^-1*csc(c+d*x)^6+1//3*a^-3*d^-1*csc(c+d*x)^3+3//5*a^-3*d^-1*csc(c+d*x)^5)
@test integrate((a+a*sin(c+d*x))^-3*cot(c+d*x)^9, x) == :(-2//5*a^-3*d^-1*csc(c+d*x)^5+-1//3*a^-3*d^-1*csc(c+d*x)^3+-1//3*a^-3*d^-1*csc(c+d*x)^6+-1//8*a^-3*d^-1*csc(c+d*x)^8+3//4*a^-3*d^-1*csc(c+d*x)^4+3//7*a^-3*d^-1*csc(c+d*x)^7)
@test integrate((a+a*sin(c+d*x))^-3*cot(c+d*x)^11, x) == :(-5//7*a^-3*d^-1*csc(c+d*x)^7+-3//4*a^-3*d^-1*csc(c+d*x)^4+-1//8*a^-3*d^-1*csc(c+d*x)^8+-1//10*a^-3*d^-1*csc(c+d*x)^10+1//3*a^-3*d^-1*csc(c+d*x)^3+1//3*a^-3*d^-1*csc(c+d*x)^9+1//5*a^-3*d^-1*csc(c+d*x)^5+5//6*a^-3*d^-1*csc(c+d*x)^6)
@test integrate((a+a*sin(c+d*x))^-3*cot(c+d*x)^13, x) == :(-8//9*a^-3*d^-1*csc(c+d*x)^9+-4//3*a^-3*d^-1*csc(c+d*x)^6+-1//3*a^-3*d^-1*csc(c+d*x)^3+-1//12*a^-3*d^-1*csc(c+d*x)^12+3//4*a^-3*d^-1*csc(c+d*x)^4+3//4*a^-3*d^-1*csc(c+d*x)^8+3//11*a^-3*d^-1*csc(c+d*x)^11+6//7*a^-3*d^-1*csc(c+d*x)^7)
@test integrate((a+a*sin(c+d*x))^-4*tan(c+d*x)^5, x) == :(-5//256*d^-1*(a^2+a^2*sin(c+d*x))^-2+-3//256*d^-1*(a^4+-1*a^4*sin(c+d*x))^-1+-1//256*d^-1*(a^4+a^4*sin(c+d*x))^-1+1//8*d^-1*(a+a*sin(c+d*x))^-4+1//256*d^-1*(a^2+-1*a^2*sin(c+d*x))^-2+-7//80*a*d^-1*(a+a*sin(c+d*x))^-5+-5//96*a^-1*d^-1*(a+a*sin(c+d*x))^-3+-1//128*a^-4*d^-1*arctanh(sin(c+d*x))+1//48*a^2*d^-1*(a+a*sin(c+d*x))^-6)
@test integrate((a+a*sin(c+d*x))^-4*tan(c+d*x)^3, x) == :(-1//8*d^-1*(a+a*sin(c+d*x))^-4+1//32*d^-1*(a^2+a^2*sin(c+d*x))^-2+1//64*d^-1*(a^4+a^4*sin(c+d*x))^-1+1//64*d^-1*(a^4+-1*a^4*sin(c+d*x))^-1+1//16*a^-1*d^-1*(a+a*sin(c+d*x))^-3+1//20*a*d^-1*(a+a*sin(c+d*x))^-5)
@test integrate((a+a*sin(c+d*x))^-4*tan(c+d*x), x) == :(-1//16*d^-1*(a^2+a^2*sin(c+d*x))^-2+-1//16*d^-1*(a^4+a^4*sin(c+d*x))^-1+1//8*d^-1*(a+a*sin(c+d*x))^-4+-1//12*a^-1*d^-1*(a+a*sin(c+d*x))^-3+1//16*a^-4*d^-1*arctanh(sin(c+d*x)))
@test integrate((a+a*sin(c+d*x))^-4*cot(c+d*x)^3, x) == :(d^-1*(a^2+a^2*sin(c+d*x))^-2+5*d^-1*(a^4+a^4*sin(c+d*x))^-1+-9*a^-4*d^-1*log(1+sin(c+d*x))+4*a^-4*d^-1*csc(c+d*x)+9*a^-4*d^-1*log(sin(c+d*x))+-1//2*a^-4*d^-1*csc(c+d*x)^2)
@test integrate((a+a*sin(c+d*x))^-4*cot(c+d*x)^7, x) == :(-8*a^-4*d^-1*log(1+sin(c+d*x))+-4*a^-4*d^-1*csc(c+d*x)^2+8*a^-4*d^-1*csc(c+d*x)+8*a^-4*d^-1*log(sin(c+d*x))+-7//4*a^-4*d^-1*csc(c+d*x)^4+-1//6*a^-4*d^-1*csc(c+d*x)^6+4//5*a^-4*d^-1*csc(c+d*x)^5+8//3*a^-4*d^-1*csc(c+d*x)^3)
@test integrate((a+a*sin(c+d*x))^-4*tan(c+d*x)^2, x) == :(-8//9*a^-4*d^-1*sec(c+d*x)^9+-4//5*a^-4*d^-1*sec(c+d*x)^5+1//3*a^-4*d^-1*tan(c+d*x)^3+8//9*a^-4*d^-1*tan(c+d*x)^9+9//5*a^-4*d^-1*tan(c+d*x)^5+12//7*a^-4*d^-1*sec(c+d*x)^7+16//7*a^-4*d^-1*tan(c+d*x)^7)
@test integrate((a+a*sin(e+f*x))^(1/2)*tan(e+f*x)^2, x) == :(5*f^-1*(a+a*sin(e+f*x))^(1/2)*sec(e+f*x)+-2*a^-1*f^-1*(a+a*sin(e+f*x))^3//2*sec(e+f*x)+-1//2*2^(1/2)*a^(1/2)*f^-1*arctanh((1/2)*2^(1/2)*a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x)))
@test integrate((a+a*sin(e+f*x))^(1/2)*cot(e+f*x)^2, x) == :(-1*a^(1/2)*f^-1*arctanh(a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))+-1*f^-1*(a+a*sin(e+f*x))^(1/2)*cot(e+f*x)+3*a*f^-1*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))
@test integrate((a+a*sin(e+f*x))^(1/2)*cot(e+f*x)^4, x) == :(11//8*a^(1/2)*f^-1*arctanh(a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))+-2*a*f^-1*(a+a*sin(e+f*x))^-1//2*cos(e+f*x)+-1//3*f^-1*(a+a*sin(e+f*x))^(1/2)*csc(e+f*x)^2*cot(e+f*x)+11//8*a*f^-1*(a+a*sin(e+f*x))^-1//2*cot(e+f*x)+-1//12*a*f^-1*(a+a*sin(e+f*x))^-1//2*cot(e+f*x)*csc(e+f*x))
@test integrate((a+a*sin(e+f*x))^3//2*tan(e+f*x)^2, x) == :(7//3*f^-1*(a+a*sin(e+f*x))^3//2*sec(e+f*x)+-2//3*a^-1*f^-1*(a+a*sin(e+f*x))^5//2*sec(e+f*x)+11//3*a^2*f^-1*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))
@test integrate((a+a*sin(e+f*x))^3//2*cot(e+f*x)^2, x) == :(-1*f^-1*(a+a*sin(e+f*x))^3//2*cot(e+f*x)+-3*a^3//2*f^-1*arctanh(a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))+5//3*a*f^-1*(a+a*sin(e+f*x))^(1/2)*cos(e+f*x)+11//3*a^2*f^-1*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))
@test integrate((a+a*sin(e+f*x))^3//2*cot(e+f*x)^4, x) == :(37//8*a^3//2*f^-1*arctanh(a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))+-8//3*a^2*f^-1*(a+a*sin(e+f*x))^-1//2*cos(e+f*x)+-2//3*a*f^-1*(a+a*sin(e+f*x))^(1/2)*cos(e+f*x)+-1//3*f^-1*(a+a*sin(e+f*x))^3//2*csc(e+f*x)^2*cot(e+f*x)+29//24*a^2*f^-1*(a+a*sin(e+f*x))^-1//2*cot(e+f*x)+-1//4*a*f^-1*(a+a*sin(e+f*x))^(1/2)*cot(e+f*x)*csc(e+f*x))
@test integrate((a+a*sin(e+f*x))^5//2*tan(e+f*x)^2, x) == :(9//5*f^-1*(a+a*sin(e+f*x))^5//2*sec(e+f*x)+-2//5*a^-1*f^-1*(a+a*sin(e+f*x))^7//2*sec(e+f*x)+31//15*a^2*f^-1*(a+a*sin(e+f*x))^(1/2)*cos(e+f*x)+124//15*a^3*f^-1*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))
@test integrate((a+a*sin(e+f*x))^5//2*cot(e+f*x)^2, x) == :(-1*f^-1*(a+a*sin(e+f*x))^5//2*cot(e+f*x)+-5*a^5//2*f^-1*arctanh(a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))+7//5*a*f^-1*(a+a*sin(e+f*x))^3//2*cos(e+f*x)+31//15*a^2*f^-1*(a+a*sin(e+f*x))^(1/2)*cos(e+f*x)+49//15*a^3*f^-1*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))
@test integrate((a+a*sin(e+f*x))^5//2*cot(e+f*x)^4, x) == :(55//8*a^5//2*f^-1*arctanh(a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))+-16//15*a^2*f^-1*(a+a*sin(e+f*x))^(1/2)*cos(e+f*x)+-9//40*a^3*f^-1*(a+a*sin(e+f*x))^-1//2*cos(e+f*x)+-2//5*a*f^-1*(a+a*sin(e+f*x))^3//2*cos(e+f*x)+-1//3*f^-1*(a+a*sin(e+f*x))^5//2*csc(e+f*x)^2*cot(e+f*x)+17//24*a^2*f^-1*(a+a*sin(e+f*x))^(1/2)*cot(e+f*x)+-5//12*a*f^-1*(a+a*sin(e+f*x))^3//2*cot(e+f*x)*csc(e+f*x))
@test integrate((a+a*sin(e+f*x))^-1//2*tan(e+f*x)^2, x) == :(-1//2*f^-1*(a+a*sin(e+f*x))^-1//2*sec(e+f*x)+3//4*a^-1*f^-1*(a+a*sin(e+f*x))^(1/2)*sec(e+f*x)+5//8*2^(1/2)*a^-1//2*f^-1*arctanh((1/2)*2^(1/2)*a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x)))
@test integrate((a+a*sin(e+f*x))^-1//2*cot(e+f*x)^2, x) == :(a^-1//2*f^-1*arctanh(a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))+-1*f^-1*(a+a*sin(e+f*x))^-1//2*cot(e+f*x))
@test integrate((a+a*sin(e+f*x))^-1//2*cot(e+f*x)^4, x) == :(-7//8*a^-1//2*f^-1*arctanh(a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))+9//8*f^-1*(a+a*sin(e+f*x))^-1//2*cot(e+f*x)+-1//3*f^-1*(a+a*sin(e+f*x))^-1//2*csc(e+f*x)^2*cot(e+f*x)+1//12*f^-1*(a+a*sin(e+f*x))^-1//2*cot(e+f*x)*csc(e+f*x))
@test integrate((a+a*sin(e+f*x))^-3//2*tan(e+f*x)^2, x) == :(-1//4*f^-1*(a+a*sin(e+f*x))^-3//2*sec(e+f*x)+1//32*f^-1*(a+a*sin(e+f*x))^-3//2*cos(e+f*x)+1//64*2^(1/2)*a^-3//2*f^-1*arctanh((1/2)*2^(1/2)*a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))+5//8*a^-1*f^-1*(a+a*sin(e+f*x))^-1//2*sec(e+f*x))
@test integrate((a+a*sin(e+f*x))^-3//2*cot(e+f*x)^2, x) == :(3*a^-3//2*f^-1*arctanh(a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))+-1*a^-1*f^-1*(a+a*sin(e+f*x))^-1//2*cot(e+f*x)+-2*2^(1/2)*a^-3//2*f^-1*arctanh((1/2)*2^(1/2)*a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x)))
@test integrate((a+a*sin(e+f*x))^-3//2*cot(e+f*x)^4, x) == :(-1//8*a^-3//2*f^-1*arctanh(a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))+-1//8*a^-1*f^-1*(a+a*sin(e+f*x))^-1//2*cot(e+f*x)+-1//3*a^-2*f^-1*(a+a*sin(e+f*x))^(1/2)*csc(e+f*x)^2*cot(e+f*x)+11//12*a^-1*f^-1*(a+a*sin(e+f*x))^-1//2*cot(e+f*x)*csc(e+f*x))
@test integrate((a+a*sin(e+f*x))^-5//2*tan(e+f*x)^2, x) == :(-1//6*f^-1*(a+a*sin(e+f*x))^-5//2*sec(e+f*x)+-11//128*a^-1*f^-1*(a+a*sin(e+f*x))^-3//2*cos(e+f*x)+-11//256*2^(1/2)*a^-5//2*f^-1*arctanh((1/2)*2^(1/2)*a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))+11//96*a^-2*f^-1*(a+a*sin(e+f*x))^-1//2*sec(e+f*x)+17//48*a^-1*f^-1*(a+a*sin(e+f*x))^-3//2*sec(e+f*x))
@test integrate((a+a*sin(e+f*x))^-5//2*cot(e+f*x)^2, x) == :(5*a^-5//2*f^-1*arctanh(a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))+-1*a^-1*f^-1*(a+a*sin(e+f*x))^-3//2*cot(e+f*x)+-2*a^-1*f^-1*(a+a*sin(e+f*x))^-3//2*cos(e+f*x)+-7//2*2^(1/2)*a^-5//2*f^-1*arctanh((1/2)*2^(1/2)*a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x)))
@test integrate((a+a*sin(e+f*x))^-5//2*cot(e+f*x)^4, x) == :(45//8*a^-5//2*f^-1*arctanh(a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))+-4*2^(1/2)*a^-5//2*f^-1*arctanh((1/2)*2^(1/2)*a^(1/2)*(a+a*sin(e+f*x))^-1//2*cos(e+f*x))+-19//8*a^-2*f^-1*(a+a*sin(e+f*x))^-1//2*cot(e+f*x)+-1//3*a^-2*f^-1*(a+a*sin(e+f*x))^-1//2*csc(e+f*x)^2*cot(e+f*x)+13//12*a^-2*f^-1*(a+a*sin(e+f*x))^-1//2*cot(e+f*x)*csc(e+f*x))
@test integrate((a+a*sin(e+f*x))^1//3*tan(e+f*x)^4, x) == :(-361//126*f^-1*(a+a*sin(e+f*x))^1//3*sec(e+f*x)+361//63*f^-1*(a+a*sin(e+f*x))^1//3*(1+-1*sin(e+f*x))*sec(e+f*x)+-1//42*f^-1*(a+a*sin(e+f*x))^-2//3*(a+-1*a*sin(e+f*x))^-1*(65*a^2+-142*a^2*sin(e+f*x))*sec(e+f*x)+-3*a^2*f^-1*(a+a*sin(e+f*x))^-2//3*(a+-1*a*sin(e+f*x))^-1*sin(e+f*x)^2*tan(e+f*x)+1//63*f^-1*(a+a*sin(e+f*x))^2//3*(2^1//3*a^1//3+-1*(a+a*sin(e+f*x))^1//3*(1+3^(1/2)))^-1*(1+-1*sin(e+f*x))*(361+361*3^(1/2))*sec(e+f*x)+3//2*a^2*f^-1*(a+a*sin(e+f*x))^-2//3*(a+-1*a*sin(e+f*x))^-1*sin(e+f*x)*tan(e+f*x)+-361//63*2^1//3*3^1//4*a^-2//3*f^-1*((2^1//3*a^1//3+-1*(a+a*sin(e+f*x))^1//3*(1+3^(1/2)))^-2*((a+a*sin(e+f*x))^2//3+2^2//3*a^2//3+2^1//3*a^1//3*(a+a*sin(e+f*x))^1//3))^(1/2)*(-1*(a+a*sin(e+f*x))^1//3*(2^1//3*a^1//3+-1*(a+a*sin(e+f*x))^1//3*(1+3^(1/2)))^-2*(-1*(a+a*sin(e+f*x))^1//3+2^1//3*a^1//3))^-1//2*(a+a*sin(e+f*x))^2//3*(-1*(a+a*sin(e+f*x))^1//3+2^1//3*a^1//3)*Elliptic.E(arccos((2^1//3*a^1//3+-1*(a+a*sin(e+f*x))^1//3*(1+3^(1/2)))^-1*(2^1//3*a^1//3+-1*(a+a*sin(e+f*x))^1//3*(1+-1*3^(1/2)))),1/2+1//4*3^(1/2))*sec(e+f*x)+-1//378*2^1//3*3^3//4*a^-2//3*f^-1*((2^1//3*a^1//3+-1*(a+a*sin(e+f*x))^1//3*(1+3^(1/2)))^-2*((a+a*sin(e+f*x))^2//3+2^2//3*a^2//3+2^1//3*a^1//3*(a+a*sin(e+f*x))^1//3))^(1/2)*(-1*(a+a*sin(e+f*x))^1//3*(2^1//3*a^1//3+-1*(a+a*sin(e+f*x))^1//3*(1+3^(1/2)))^-2*(-1*(a+a*sin(e+f*x))^1//3+2^1//3*a^1//3))^-1//2*(a+a*sin(e+f*x))^2//3*(361+-361*3^(1/2))*(-1*(a+a*sin(e+f*x))^1//3+2^1//3*a^1//3)*Elliptic.F(arccos((2^1//3*a^1//3+-1*(a+a*sin(e+f*x))^1//3*(1+3^(1/2)))^-1*(2^1//3*a^1//3+-1*(a+a*sin(e+f*x))^1//3*(1+-1*3^(1/2)))),1/2+1//4*3^(1/2))*sec(e+f*x))
@test integrate((a+a*sin(e+f*x))^-1//3*tan(e+f*x)^4, x) == :(973//396*f^-1*(a+a*sin(e+f*x))^-1//3*sec(e+f*x)+-973//495*f^-1*(a+a*sin(e+f*x))^-1//3*(1+-1*sin(e+f*x))*sec(e+f*x)+-1//132*f^-1*(1+-1*sin(e+f*x))^-1*(a+a*sin(e+f*x))^-4//3*(95a+356*a*sin(e+f*x))*sec(e+f*x)+3*a^2*f^-1*(a+a*sin(e+f*x))^-4//3*(a+-1*a*sin(e+f*x))^-1*sin(e+f*x)^2*tan(e+f*x)+3//4*a^2*f^-1*(a+a*sin(e+f*x))^-4//3*(a+-1*a*sin(e+f*x))^-1*sin(e+f*x)*tan(e+f*x)+973//2970*2^2//3*3^3//4*a^-4//3*f^-1*((2^1//3*a^1//3+-1*(a+a*sin(e+f*x))^1//3*(1+3^(1/2)))^-2*((a+a*sin(e+f*x))^2//3+2^2//3*a^2//3+2^1//3*a^1//3*(a+a*sin(e+f*x))^1//3))^(1/2)*(-1*(a+a*sin(e+f*x))^1//3*(2^1//3*a^1//3+-1*(a+a*sin(e+f*x))^1//3*(1+3^(1/2)))^-2*(-1*(a+a*sin(e+f*x))^1//3+2^1//3*a^1//3))^-1//2*(a+a*sin(e+f*x))^2//3*(-1*(a+a*sin(e+f*x))^1//3+2^1//3*a^1//3)*Elliptic.F(arccos((2^1//3*a^1//3+-1*(a+a*sin(e+f*x))^1//3*(1+3^(1/2)))^-1*(2^1//3*a^1//3+-1*(a+a*sin(e+f*x))^1//3*(1+-1*3^(1/2)))),1/2+1//4*3^(1/2))*sec(e+f*x))
@test integrate(tan(c+d*x)^3*(a+b*sin(c+d*x)), x) == :((1/2)*d^-1*tan(c+d*x)^2*(a+b*sin(c+d*x))+1//4*d^-1*(-3b+2a)*log(1+sin(c+d*x))+1//4*d^-1*(2a+3b)*log(1+-1*sin(c+d*x))+3//2*b*d^-1*sin(c+d*x))
@test integrate((a+b*sin(c+d*x))*tan(c+d*x), x) == :(-1*b*d^-1*sin(c+d*x)+-1//2*d^-1*(a+b)*log(1+-1*sin(c+d*x))+-1//2*d^-1*(a+-1b)*log(1+sin(c+d*x)))
@test integrate((a+b*sin(c+d*x))*cot(c+d*x), x) == :(a*d^-1*log(sin(c+d*x))+b*d^-1*sin(c+d*x))
@test integrate(cot(c+d*x)^3*(a+b*sin(c+d*x)), x) == :(-1*a*d^-1*log(sin(c+d*x))+-1*b*d^-1*csc(c+d*x)+-1*b*d^-1*sin(c+d*x)+-1//2*a*d^-1*csc(c+d*x)^2)
@test integrate(cot(c+d*x)^5*(a+b*sin(c+d*x)), x) == :(a*d^-1*csc(c+d*x)^2+a*d^-1*log(sin(c+d*x))+b*d^-1*sin(c+d*x)+2*b*d^-1*csc(c+d*x)+-1//3*b*d^-1*csc(c+d*x)^3+-1//4*a*d^-1*csc(c+d*x)^4)
@test integrate(tan(c+d*x)^4*(a+b*sin(c+d*x)), x) == :(a*x+-1*a*d^-1*tan(c+d*x)+-1*b*d^-1*cos(c+d*x)+-2*b*d^-1*sec(c+d*x)+1//3*a*d^-1*tan(c+d*x)^3+1//3*b*d^-1*sec(c+d*x)^3)
@test integrate(tan(c+d*x)^2*(a+b*sin(c+d*x)), x) == :(-1*a*x+a*d^-1*tan(c+d*x)+b*d^-1*cos(c+d*x)+b*d^-1*sec(c+d*x))
@test integrate(cot(c+d*x)^2*(a+b*sin(c+d*x)), x) == :(-1*a*x+b*d^-1*cos(c+d*x)+-1*a*d^-1*cot(c+d*x)+-1*b*d^-1*arctanh(cos(c+d*x)))
@test integrate(cot(c+d*x)^4*(a+b*sin(c+d*x)), x) == :(a*x+a*d^-1*cot(c+d*x)+-3//2*b*d^-1*cos(c+d*x)+-1//3*a*d^-1*cot(c+d*x)^3+3//2*b*d^-1*arctanh(cos(c+d*x))+-1//2*b*d^-1*cot(c+d*x)^2*cos(c+d*x))
@test integrate(cot(c+d*x)^6*(a+b*sin(c+d*x)), x) == :(-1*a*x+-1*a*d^-1*cot(c+d*x)+-15//8*b*d^-1*arctanh(cos(c+d*x))+-1//5*a*d^-1*cot(c+d*x)^5+1//3*a*d^-1*cot(c+d*x)^3+15//8*b*d^-1*cos(c+d*x)+-1//4*b*d^-1*cot(c+d*x)^4*cos(c+d*x)+5//8*b*d^-1*cot(c+d*x)^2*cos(c+d*x))
@test integrate((a+b*sin(c+d*x))^2*tan(c+d*x)^3, x) == :((1/2)*b^2*d^-1*sin(c+d*x)^2+(1/2)*d^-1*(a+b*sin(c+d*x))^2*sec(c+d*x)^2+(1/2)*d^-1*(a+b)*(a+2b)*log(1+-1*sin(c+d*x))+(1/2)*d^-1*(a+-1b)*(a+-2b)*log(1+sin(c+d*x))+2*a*b*d^-1*sin(c+d*x))
@test integrate((a+b*sin(c+d*x))^2*tan(c+d*x), x) == :(-1//2*b^2*d^-1*sin(c+d*x)^2+-1//2*d^-1*(a+b)^2*log(1+-1*sin(c+d*x))+-1//2*d^-1*(a+-1b)^2*log(1+sin(c+d*x))+-2*a*b*d^-1*sin(c+d*x))
@test integrate((a+b*sin(c+d*x))^2*cot(c+d*x), x) == :(a^2*d^-1*log(sin(c+d*x))+(1/2)*b^2*d^-1*sin(c+d*x)^2+2*a*b*d^-1*sin(c+d*x))
@test integrate((a+b*sin(c+d*x))^2*cot(c+d*x)^3, x) == :(-1*d^-1*(a^2+-1*b^2)*log(sin(c+d*x))+-1//2*a^2*d^-1*csc(c+d*x)^2+-1//2*b^2*d^-1*sin(c+d*x)^2+-2*a*b*d^-1*csc(c+d*x)+-2*a*b*d^-1*sin(c+d*x))
@test integrate((a+b*sin(c+d*x))^2*cot(c+d*x)^5, x) == :(d^-1*(a^2+-2*b^2)*log(sin(c+d*x))+(1/2)*b^2*d^-1*sin(c+d*x)^2+(1/2)*d^-1*csc(c+d*x)^2*(-1*b^2+2*a^2)+-1//4*a^2*d^-1*csc(c+d*x)^4+2*a*b*d^-1*sin(c+d*x)+4*a*b*d^-1*csc(c+d*x)+-2//3*a*b*d^-1*csc(c+d*x)^3)
@test integrate((a+b*sin(c+d*x))^2*tan(c+d*x)^4, x) == :(x*a^2+5//2*x*b^2+-1*a^2*d^-1*tan(c+d*x)+-5//2*b^2*d^-1*tan(c+d*x)+1//3*a^2*d^-1*tan(c+d*x)^3+5//6*b^2*d^-1*tan(c+d*x)^3+-4*a*b*d^-1*sec(c+d*x)+-2*a*b*d^-1*cos(c+d*x)+-1//2*b^2*d^-1*sin(c+d*x)^2*tan(c+d*x)^3+2//3*a*b*d^-1*sec(c+d*x)^3)
@test integrate((a+b*sin(c+d*x))^2*tan(c+d*x)^2, x) == :(-1*x*a^2+-3//2*x*b^2+a^2*d^-1*tan(c+d*x)+3//2*b^2*d^-1*tan(c+d*x)+2*a*b*d^-1*cos(c+d*x)+2*a*b*d^-1*sec(c+d*x)+-1//2*b^2*d^-1*sin(c+d*x)^2*tan(c+d*x))
@test integrate((a+b*sin(c+d*x))^2*cot(c+d*x)^2, x) == :((1/2)*x*b^2+-1*x*a^2+-1*a^2*d^-1*cot(c+d*x)+(1/2)*b^2*d^-1*cos(c+d*x)*sin(c+d*x)+-2*a*b*d^-1*arctanh(cos(c+d*x))+2*a*b*d^-1*cos(c+d*x))
@test integrate((a+b*sin(c+d*x))^2*cot(c+d*x)^4, x) == :(x*a^2+-3//2*x*b^2+a^2*d^-1*cot(c+d*x)+-3//2*b^2*d^-1*cot(c+d*x)+-1//3*a^2*d^-1*cot(c+d*x)^3+(1/2)*b^2*d^-1*cos(c+d*x)^2*cot(c+d*x)+-3*a*b*d^-1*cos(c+d*x)+3*a*b*d^-1*arctanh(cos(c+d*x))+-1*a*b*d^-1*cot(c+d*x)^2*cos(c+d*x))
@test integrate((a+b*sin(c+d*x))^2*cot(c+d*x)^6, x) == :(-1*x*a^2+5//2*x*b^2+-1*a^2*d^-1*cot(c+d*x)+-5//6*b^2*d^-1*cot(c+d*x)^3+-1//5*a^2*d^-1*cot(c+d*x)^5+1//3*a^2*d^-1*cot(c+d*x)^3+5//2*b^2*d^-1*cot(c+d*x)+(1/2)*b^2*d^-1*cos(c+d*x)^2*cot(c+d*x)^3+-15//4*a*b*d^-1*arctanh(cos(c+d*x))+15//4*a*b*d^-1*cos(c+d*x)+-1//2*a*b*d^-1*cot(c+d*x)^4*cos(c+d*x)+5//4*a*b*d^-1*cot(c+d*x)^2*cos(c+d*x))
@test integrate((a+b*sin(c+d*x))^3*tan(c+d*x)^3, x) == :((1/2)*d^-1*(a+b*sin(c+d*x))^3*sec(c+d*x)^2+1//3*b^3*d^-1*sin(c+d*x)^3+(1/2)*b*d^-1*(5*b^2+6*a^2)*sin(c+d*x)+1//4*d^-1*(a+b)^2*(2a+5b)*log(1+-1*sin(c+d*x))+1//4*d^-1*(a+-1b)^2*(-5b+2a)*log(1+sin(c+d*x))+3//2*a*b^2*d^-1*sin(c+d*x)^2)
@test integrate((a+b*sin(c+d*x))^3*tan(c+d*x), x) == :(-1//2*d^-1*(a+b)^3*log(1+-1*sin(c+d*x))+-1//2*d^-1*(a+-1b)^3*log(1+sin(c+d*x))+-1//3*b^3*d^-1*sin(c+d*x)^3+-1*b*d^-1*(b^2+3*a^2)*sin(c+d*x)+-3//2*a*b^2*d^-1*sin(c+d*x)^2)
@test integrate((a+b*sin(c+d*x))^3*cot(c+d*x), x) == :(a^3*d^-1*log(sin(c+d*x))+1//3*b^3*d^-1*sin(c+d*x)^3+3*b*a^2*d^-1*sin(c+d*x)+3//2*a*b^2*d^-1*sin(c+d*x)^2)
@test integrate((a+b*sin(c+d*x))^3*cot(c+d*x)^3, x) == :(-1//2*a^3*d^-1*csc(c+d*x)^2+-1//3*b^3*d^-1*sin(c+d*x)^3+-1*a*d^-1*(a^2+-3*b^2)*log(sin(c+d*x))+-1*b*d^-1*(-1*b^2+3*a^2)*sin(c+d*x)+-3*b*a^2*d^-1*csc(c+d*x)+-3//2*a*b^2*d^-1*sin(c+d*x)^2)
@test integrate((a+b*sin(c+d*x))^3*cot(c+d*x)^5, x) == :(-1//4*a^3*d^-1*csc(c+d*x)^4+1//3*b^3*d^-1*sin(c+d*x)^3+a*d^-1*(a^2+-6*b^2)*log(sin(c+d*x))+b*d^-1*(-1*b^2+6*a^2)*csc(c+d*x)+b*d^-1*(-2*b^2+3*a^2)*sin(c+d*x)+(1/2)*a*d^-1*csc(c+d*x)^2*(-3*b^2+2*a^2)+-1*b*a^2*d^-1*csc(c+d*x)^3+3//2*a*b^2*d^-1*sin(c+d*x)^2)
@test integrate((a+b*sin(c+d*x))^3*tan(c+d*x)^4, x) == :(x*a^3+-1*a^3*d^-1*tan(c+d*x)+-3*b^3*d^-1*cos(c+d*x)+-3*b^3*d^-1*sec(c+d*x)+1//3*a^3*d^-1*tan(c+d*x)^3+1//3*b^3*d^-1*cos(c+d*x)^3+1//3*b^3*d^-1*sec(c+d*x)^3+15//2*a*x*b^2+b*a^2*d^-1*sec(c+d*x)^3+-6*b*a^2*d^-1*sec(c+d*x)+-3*b*a^2*d^-1*cos(c+d*x)+-15//2*a*b^2*d^-1*tan(c+d*x)+5//2*a*b^2*d^-1*tan(c+d*x)^3+-3//2*a*b^2*d^-1*sin(c+d*x)^2*tan(c+d*x)^3)
@test integrate((a+b*sin(c+d*x))^3*tan(c+d*x)^2, x) == :(-1*x*a^3+a^3*d^-1*tan(c+d*x)+b^3*d^-1*sec(c+d*x)+2*b^3*d^-1*cos(c+d*x)+-9//2*a*x*b^2+-1//3*b^3*d^-1*cos(c+d*x)^3+3*b*a^2*d^-1*cos(c+d*x)+3*b*a^2*d^-1*sec(c+d*x)+9//2*a*b^2*d^-1*tan(c+d*x)+-3//2*a*b^2*d^-1*sin(c+d*x)^2*tan(c+d*x))
@test integrate((a+b*sin(c+d*x))^3*cot(c+d*x)^2, x) == :(-1*x*a^3+-1*a^3*d^-1*cot(c+d*x)+-1//3*b^3*d^-1*cos(c+d*x)^3+3//2*a*x*b^2+-3*b*a^2*d^-1*arctanh(cos(c+d*x))+3*b*a^2*d^-1*cos(c+d*x)+3//2*a*b^2*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate((a+b*sin(c+d*x))^3*cot(c+d*x)^4, x) == :(x*a^3+a^3*d^-1*cot(c+d*x)+b^3*d^-1*cos(c+d*x)+-1*b^3*d^-1*arctanh(cos(c+d*x))+-9//2*a*x*b^2+-1//3*a^3*d^-1*cot(c+d*x)^3+1//3*b^3*d^-1*cos(c+d*x)^3+-9//2*a*b^2*d^-1*cot(c+d*x)+-9//2*b*a^2*d^-1*cos(c+d*x)+9//2*b*a^2*d^-1*arctanh(cos(c+d*x))+-3//2*b*a^2*d^-1*cot(c+d*x)^2*cos(c+d*x)+3//2*a*b^2*d^-1*cos(c+d*x)^2*cot(c+d*x))
@test integrate((a+b*sin(c+d*x))^3*cot(c+d*x)^6, x) == :(-1*x*a^3+-1*a^3*d^-1*cot(c+d*x)+-5//2*b^3*d^-1*cos(c+d*x)+-5//6*b^3*d^-1*cos(c+d*x)^3+-1//5*a^3*d^-1*cot(c+d*x)^5+1//3*a^3*d^-1*cot(c+d*x)^3+5//2*b^3*d^-1*arctanh(cos(c+d*x))+15//2*a*x*b^2+-45//8*b*a^2*d^-1*arctanh(cos(c+d*x))+-5//2*a*b^2*d^-1*cot(c+d*x)^3+-1//2*b^3*d^-1*cos(c+d*x)^3*cot(c+d*x)^2+15//2*a*b^2*d^-1*cot(c+d*x)+45//8*b*a^2*d^-1*cos(c+d*x)+-3//4*b*a^2*d^-1*cot(c+d*x)^4*cos(c+d*x)+3//2*a*b^2*d^-1*cos(c+d*x)^2*cot(c+d*x)^3+15//8*b*a^2*d^-1*cot(c+d*x)^2*cos(c+d*x))
@test integrate((a+b*sin(c+d*x))^-1*tan(c+d*x)^5, x) == :(a^5*d^-1*(a^2+-1*b^2)^-3*log(a+b*sin(c+d*x))+d^-1*(-4*b^2+4*a^2)^-1*sec(c+d*x)^4*(a+-1*b*sin(c+d*x))+-1//8*d^-1*(a^2+-1*b^2)^-2*sec(c+d*x)^2*(4*a*(-1*b^2+2*a^2)+-1*b*(-5*b^2+9*a^2)*sin(c+d*x))+-1//16*d^-1*(a+b)^-3*(3*b^2+8*a^2+9*a*b)*log(1+-1*sin(c+d*x))+-1//16*d^-1*(a+-1b)^-3*(3*b^2+8*a^2+-9*a*b)*log(1+sin(c+d*x)))
@test integrate((a+b*sin(c+d*x))^-1*tan(c+d*x)^3, x) == :(d^-1*(-2*b^2+2*a^2)^-1*sec(c+d*x)^2*(a+-1*b*sin(c+d*x))+-1*a^3*d^-1*(a^2+-1*b^2)^-2*log(a+b*sin(c+d*x))+1//4*d^-1*(a+b)^-2*(b+2a)*log(1+-1*sin(c+d*x))+1//4*d^-1*(a+-1b)^-2*(-1b+2a)*log(1+sin(c+d*x)))
@test integrate((a+b*sin(c+d*x))^-1*tan(c+d*x), x) == :(-1*d^-1*(-2b+2a)^-1*log(1+sin(c+d*x))+-1*d^-1*(2a+2b)^-1*log(1+-1*sin(c+d*x))+a*d^-1*(a^2+-1*b^2)^-1*log(a+b*sin(c+d*x)))
@test integrate((a+b*sin(c+d*x))^-1*cot(c+d*x), x) == :(a^-1*d^-1*log(sin(c+d*x))+-1*a^-1*d^-1*log(a+b*sin(c+d*x)))
@test integrate((a+b*sin(c+d*x))^-1*cot(c+d*x)^3, x) == :(-1//2*a^-1*d^-1*csc(c+d*x)^2+b*a^-2*d^-1*csc(c+d*x)+a^-3*d^-1*(a^2+-1*b^2)*log(a+b*sin(c+d*x))+-1*a^-3*d^-1*(a^2+-1*b^2)*log(sin(c+d*x)))
@test integrate((a+b*sin(c+d*x))^-1*cot(c+d*x)^5, x) == :(-1//4*a^-1*d^-1*csc(c+d*x)^4+a^-5*d^-1*(a^2+-1*b^2)^2*log(sin(c+d*x))+(1/2)*a^-3*d^-1*csc(c+d*x)^2*(-1*b^2+2*a^2)+-1*a^-5*d^-1*(a^2+-1*b^2)^2*log(a+b*sin(c+d*x))+1//3*b*a^-2*d^-1*csc(c+d*x)^3+-1*b*a^-4*d^-1*(-1*b^2+2*a^2)*csc(c+d*x))
@test integrate((a+b*sin(c+d*x))^-1*tan(c+d*x)^4, x) == :(a*d^-1*(-3*b^2+3*a^2)^-1*tan(c+d*x)^3+b*d^-1*(a^2+-1*b^2)^-1*sec(c+d*x)+-1*b*d^-1*(-3*b^2+3*a^2)^-1*sec(c+d*x)^3+-1*a^3*d^-1*(a^2+-1*b^2)^-2*tan(c+d*x)+2*a^4*d^-1*(a^2+-1*b^2)^-5//2*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x)))+b*a^2*d^-1*(a^2+-1*b^2)^-2*sec(c+d*x))
@test integrate((a+b*sin(c+d*x))^-1*tan(c+d*x)^2, x) == :(a*d^-1*(a^2+-1*b^2)^-1*tan(c+d*x)+-1*b*d^-1*(a^2+-1*b^2)^-1*sec(c+d*x)+-2*a^2*d^-1*(a^2+-1*b^2)^-3//2*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x))))
@test integrate((a+b*sin(c+d*x))^-1*cot(c+d*x)^2, x) == :(-1*a^-1*d^-1*cot(c+d*x)+b*a^-2*d^-1*arctanh(cos(c+d*x))+-2*a^-2*d^-1*(a^2+-1*b^2)^(1/2)*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x))))
@test integrate((a+b*sin(c+d*x))^-1*cot(c+d*x)^4, x) == :(2*a^-4*d^-1*(a^2+-1*b^2)^3//2*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x)))+-1//3*a^-1*d^-1*csc(c+d*x)^2*cot(c+d*x)+1//3*a^-3*d^-1*(-3*b^2+4*a^2)*cot(c+d*x)+(1/2)*b*a^-2*d^-1*cot(c+d*x)*csc(c+d*x)+-1//2*b*a^-4*d^-1*(-2*b^2+3*a^2)*arctanh(cos(c+d*x)))
@test integrate((a+b*sin(c+d*x))^-1*cot(c+d*x)^6, x) == :(-1*b^-1*d^-1*cot(c+d*x)*csc(c+d*x)+-2*a^-6*d^-1*(a^2+-1*b^2)^5//2*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x)))+-1//5*a^-1*d^-1*csc(c+d*x)^4*cot(c+d*x)+-1//15*a^-5*d^-1*(15*b^4+23*a^4+-35*a^2*b^2)*cot(c+d*x)+(1/2)*a*b^-2*d^-1*csc(c+d*x)^2*cot(c+d*x)+1//4*b*a^-2*d^-1*csc(c+d*x)^3*cot(c+d*x)+1//8*b*a^-6*d^-1*(8*b^4+15*a^4+-20*a^2*b^2)*arctanh(cos(c+d*x))+-1//30*a^-3*b^-2*d^-1*csc(c+d*x)^2*(10*b^4+15*a^4+-22*a^2*b^2)*cot(c+d*x)+1//8*a^-4*b^-1*d^-1*(4*b^4+8*a^4+-9*a^2*b^2)*cot(c+d*x)*csc(c+d*x))
@test integrate((a+b*sin(c+d*x))^-2*tan(c+d*x)^5, x) == :(-1*a^5*d^-1*(a+b*sin(c+d*x))^-1*(a^2+-1*b^2)^-3+-1//4*d^-1*(a^2+-1*b^2)^-3*sec(c+d*x)^2*(-2*b^4+4*a^4+6*a^2*b^2+-1*a*b*(-1*b^2+9*a^2)*sin(c+d*x))+1//4*d^-1*(a^2+-1*b^2)^-2*sec(c+d*x)^4*(a^2+b^2+-2*a*b*sin(c+d*x))+a^4*d^-1*(a^2+-1*b^2)^-4*(a^2+5*b^2)*log(a+b*sin(c+d*x))+-1//8*a*d^-1*(a+b)^-4*(b+4a)*log(1+-1*sin(c+d*x))+-1//8*a*d^-1*(a+-1b)^-4*(-1b+4a)*log(1+sin(c+d*x)))
@test integrate((a+b*sin(c+d*x))^-2*tan(c+d*x)^3, x) == :(a^3*d^-1*(a+b*sin(c+d*x))^-1*(a^2+-1*b^2)^-2+(1/2)*a*d^-1*(a+b)^-3*log(1+-1*sin(c+d*x))+(1/2)*a*d^-1*(a+-1b)^-3*log(1+sin(c+d*x))+(1/2)*d^-1*(a^2+-1*b^2)^-2*sec(c+d*x)^2*(a^2+b^2+-2*a*b*sin(c+d*x))+-1*a^2*d^-1*(a^2+-1*b^2)^-3*(a^2+3*b^2)*log(a+b*sin(c+d*x)))
@test integrate((a+b*sin(c+d*x))^-2*tan(c+d*x), x) == :(-1//2*d^-1*(a+b)^-2*log(1+-1*sin(c+d*x))+-1//2*d^-1*(a+-1b)^-2*log(1+sin(c+d*x))+d^-1*(a^2+-1*b^2)^-2*(a^2+b^2)*log(a+b*sin(c+d*x))+-1*a*d^-1*(a+b*sin(c+d*x))^-1*(a^2+-1*b^2)^-1)
@test integrate((a+b*sin(c+d*x))^-2*cot(c+d*x), x) == :(a^-1*d^-1*(a+b*sin(c+d*x))^-1+a^-2*d^-1*log(sin(c+d*x))+-1*a^-2*d^-1*log(a+b*sin(c+d*x)))
@test integrate((a+b*sin(c+d*x))^-2*cot(c+d*x)^3, x) == :(-1//2*a^-2*d^-1*csc(c+d*x)^2+a^-4*d^-1*(a^2+-3*b^2)*log(a+b*sin(c+d*x))+-1*a^-4*d^-1*(a^2+-3*b^2)*log(sin(c+d*x))+-1*a^-3*d^-1*(a+b*sin(c+d*x))^-1*(a^2+-1*b^2)+2*b*a^-3*d^-1*csc(c+d*x))
@test integrate((a+b*sin(c+d*x))^-2*cot(c+d*x)^5, x) == :(-1//4*a^-2*d^-1*csc(c+d*x)^4+a^-6*d^-1*(a^4+5*b^4+-6*a^2*b^2)*log(sin(c+d*x))+a^-5*d^-1*(a+b*sin(c+d*x))^-1*(a^2+-1*b^2)^2+(1/2)*a^-4*d^-1*csc(c+d*x)^2*(-3*b^2+2*a^2)+-1*a^-6*d^-1*(a^4+5*b^4+-6*a^2*b^2)*log(a+b*sin(c+d*x))+2//3*b*a^-3*d^-1*csc(c+d*x)^3+-4*b*a^-5*d^-1*(a^2+-1*b^2)*csc(c+d*x))
@test integrate((a+b*sin(c+d*x))^-2*tan(c+d*x)^4, x) == :(2*a^5*d^-1*(a^2+-1*b^2)^-7//2*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x)))+-1//12*d^-1*(1+sin(c+d*x))^-1*(a+-1b)^-2*cos(c+d*x)+-1//12*d^-1*(1+sin(c+d*x))^-2*(a+-1b)^-2*cos(c+d*x)+1//12*d^-1*(1+-1*sin(c+d*x))^-1*(a+b)^-2*cos(c+d*x)+1//12*d^-1*(1+-1*sin(c+d*x))^-2*(a+b)^-2*cos(c+d*x)+8*a^3*b^2*d^-1*(a^2+-1*b^2)^-7//2*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x)))+-1//4*d^-1*(1+-1*sin(c+d*x))^-1*(a+b)^-3*(b+3a)*cos(c+d*x)+1//4*d^-1*(1+sin(c+d*x))^-1*(a+-1b)^-3*(-1b+3a)*cos(c+d*x)+b*a^4*d^-1*(a+b*sin(c+d*x))^-1*(a^2+-1*b^2)^-3*cos(c+d*x))
@test integrate((a+b*sin(c+d*x))^-2*tan(c+d*x)^2, x) == :((1/2)*d^-1*(1+-1*sin(c+d*x))^-1*(a+b)^-2*cos(c+d*x)+-2*a^3*d^-1*(a^2+-1*b^2)^-5//2*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x)))+-1//2*d^-1*(1+sin(c+d*x))^-1*(a+-1b)^-2*cos(c+d*x)+-4*a*b^2*d^-1*(a^2+-1*b^2)^-5//2*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x)))+-1*b*a^2*d^-1*(a+b*sin(c+d*x))^-1*(a^2+-1*b^2)^-2*cos(c+d*x))
@test integrate((a+b*sin(c+d*x))^-2*cot(c+d*x)^2, x) == :(-2*a^-2*d^-1*cot(c+d*x)+a^-1*d^-1*(a+b*sin(c+d*x))^-1*cot(c+d*x)+2*b*a^-3*d^-1*arctanh(cos(c+d*x))+-1*a^-3*d^-1*(a^2+-1*b^2)^-1//2*(-4*b^2+2*a^2)*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x))))
@test integrate((a+b*sin(c+d*x))^-2*cot(c+d*x)^4, x) == :(1//3*a^-4*d^-1*(-12*b^2+7*a^2)*cot(c+d*x)+a^-5*d^-1*(a^2+-1*b^2)^-1//2*(2*a^4+8*b^4+-10*a^2*b^2)*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x)))+-1*b*a^-5*d^-1*(-4*b^2+3*a^2)*arctanh(cos(c+d*x))+-1//3*a^-1*d^-1*(a+b*sin(c+d*x))^-1*csc(c+d*x)^2*cot(c+d*x)+-1*a^-3*b^-1*d^-1*(a^2+-2*b^2)*cot(c+d*x)*csc(c+d*x)+1//3*a^-2*b^-1*d^-1*(a+b*sin(c+d*x))^-1*(-4*b^2+3*a^2)*cot(c+d*x)*csc(c+d*x))
@test integrate((a+b*sin(c+d*x))^-2*cot(c+d*x)^6, x) == :(-1//15*a^-6*d^-1*(38*a^4+90*b^4+-135*a^2*b^2)*cot(c+d*x)+-1*a^-7*d^-1*(a^2+-1*b^2)^3//2*(-12*b^2+2*a^2)*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x)))+-1//2*b^-1*d^-1*(a+b*sin(c+d*x))^-1*cot(c+d*x)*csc(c+d*x)+-1//5*a^-1*d^-1*(a+b*sin(c+d*x))^-1*csc(c+d*x)^4*cot(c+d*x)+1//4*b*a^-7*d^-1*(15*a^4+24*b^4+-40*a^2*b^2)*arctanh(cos(c+d*x))+-1//30*a^-4*b^-2*d^-1*csc(c+d*x)^2*(15*a^4+60*b^4+-82*a^2*b^2)*cot(c+d*x)+1//4*a^-5*b^-1*d^-1*(4*a^4+12*b^4+-17*a^2*b^2)*cot(c+d*x)*csc(c+d*x)+1//6*a*b^-2*d^-1*(a+b*sin(c+d*x))^-1*csc(c+d*x)^2*cot(c+d*x)+3//10*b*a^-2*d^-1*(a+b*sin(c+d*x))^-1*csc(c+d*x)^3*cot(c+d*x)+1//6*a^-3*b^-2*d^-1*(a+b*sin(c+d*x))^-1*csc(c+d*x)^2*(2*a^4+9*b^4+-12*a^2*b^2)*cot(c+d*x))
@test integrate((a+b*sin(c+d*x))^-3*tan(c+d*x)^5, x) == :(-1//2*a^5*d^-1*(a+b*sin(c+d*x))^-2*(a^2+-1*b^2)^-3+-1//8*d^-1*(a^2+-1*b^2)^-4*sec(c+d*x)^2*(8*a^3*(a^2+5*b^2)+-1*b*(-1*b^4+27*a^4+22*a^2*b^2)*sin(c+d*x))+-1//16*d^-1*(a+b)^-5*(-1*b^2+8*a^2+-5*a*b)*log(1+-1*sin(c+d*x))+-1//16*d^-1*(a+-1b)^-5*(-1*b^2+8*a^2+5*a*b)*log(1+sin(c+d*x))+1//4*d^-1*(a^2+-1*b^2)^-3*sec(c+d*x)^4*(a*(a^2+3*b^2)+-1*b*(b^2+3*a^2)*sin(c+d*x))+a^3*d^-1*(a^2+-1*b^2)^-5*(a^4+10*b^4+13*a^2*b^2)*log(a+b*sin(c+d*x))+-1*a^4*d^-1*(a+b*sin(c+d*x))^-1*(a^2+-1*b^2)^-4*(a^2+5*b^2))
@test integrate((a+b*sin(c+d*x))^-3*tan(c+d*x)^3, x) == :((1/2)*a^3*d^-1*(a+b*sin(c+d*x))^-2*(a^2+-1*b^2)^-2+(1/2)*d^-1*(a^2+-1*b^2)^-3*sec(c+d*x)^2*(a*(a^2+3*b^2)+-1*b*(b^2+3*a^2)*sin(c+d*x))+1//4*d^-1*(a+b)^-4*(-1b+2a)*log(1+-1*sin(c+d*x))+1//4*d^-1*(a+-1b)^-4*(b+2a)*log(1+sin(c+d*x))+a^2*d^-1*(a+b*sin(c+d*x))^-1*(a^2+-1*b^2)^-3*(a^2+3*b^2)+-1*a*d^-1*(a^2+-1*b^2)^-4*(a^4+3*b^4+8*a^2*b^2)*log(a+b*sin(c+d*x)))
@test integrate((a+b*sin(c+d*x))^-3*tan(c+d*x), x) == :(-1//2*d^-1*(a+b)^-3*log(1+-1*sin(c+d*x))+-1//2*d^-1*(a+-1b)^-3*log(1+sin(c+d*x))+-1*a*d^-1*(a+b*sin(c+d*x))^-2*(-2*b^2+2*a^2)^-1+-1*d^-1*(a+b*sin(c+d*x))^-1*(a^2+-1*b^2)^-2*(a^2+b^2)+a*d^-1*(a^2+-1*b^2)^-3*(a^2+3*b^2)*log(a+b*sin(c+d*x)))
@test integrate((a+b*sin(c+d*x))^-3*cot(c+d*x), x) == :(a^-3*d^-1*log(sin(c+d*x))+a^-2*d^-1*(a+b*sin(c+d*x))^-1+(1/2)*a^-1*d^-1*(a+b*sin(c+d*x))^-2+-1*a^-3*d^-1*log(a+b*sin(c+d*x)))
@test integrate((a+b*sin(c+d*x))^-3*cot(c+d*x)^3, x) == :(-1//2*a^-3*d^-1*csc(c+d*x)^2+a^-5*d^-1*(a^2+-6*b^2)*log(a+b*sin(c+d*x))+-1*a^-5*d^-1*(a^2+-6*b^2)*log(sin(c+d*x))+-1*a^-4*d^-1*(a+b*sin(c+d*x))^-1*(a^2+-3*b^2)+3*b*a^-4*d^-1*csc(c+d*x)+-1//2*a^-3*d^-1*(a+b*sin(c+d*x))^-2*(a^2+-1*b^2))
@test integrate((a+b*sin(c+d*x))^-3*cot(c+d*x)^5, x) == :(-1//4*a^-3*d^-1*csc(c+d*x)^4+b*a^-4*d^-1*csc(c+d*x)^3+a^-7*d^-1*(a^4+15*b^4+-12*a^2*b^2)*log(sin(c+d*x))+a^-6*d^-1*(a+b*sin(c+d*x))^-1*(a^4+5*b^4+-6*a^2*b^2)+a^-5*d^-1*csc(c+d*x)^2*(a^2+-3*b^2)+(1/2)*a^-5*d^-1*(a+b*sin(c+d*x))^-2*(a^2+-1*b^2)^2+-1*a^-7*d^-1*(a^4+15*b^4+-12*a^2*b^2)*log(a+b*sin(c+d*x))+-2*b*a^-6*d^-1*(-5*b^2+3*a^2)*csc(c+d*x))
@test integrate((a+b*sin(c+d*x))^-3*tan(c+d*x)^4, x) == :(-1//12*d^-1*(1+sin(c+d*x))^-1*(a+-1b)^-3*cos(c+d*x)+-1//12*d^-1*(1+sin(c+d*x))^-2*(a+-1b)^-3*cos(c+d*x)+1//12*d^-1*(1+-1*sin(c+d*x))^-1*(a+b)^-3*cos(c+d*x)+1//12*d^-1*(1+-1*sin(c+d*x))^-2*(a+b)^-3*cos(c+d*x)+a^4*d^-1*(a^2+-1*b^2)^-9//2*(b^2+2*a^2)*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x)))+8*a^4*b^2*d^-1*(a^2+-1*b^2)^-9//2*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x)))+-3//4*a*d^-1*(1+-1*sin(c+d*x))^-1*(a+b)^-4*cos(c+d*x)+3//4*a*d^-1*(1+sin(c+d*x))^-1*(a+-1b)^-4*cos(c+d*x)+(1/2)*b*a^4*d^-1*(a+b*sin(c+d*x))^-2*(a^2+-1*b^2)^-3*cos(c+d*x)+4*a^3*b^3*d^-1*(a+b*sin(c+d*x))^-1*(a^2+-1*b^2)^-4*cos(c+d*x)+12*a^2*b^2*d^-1*(a^2+-1*b^2)^-9//2*(a^2+b^2)*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x)))+3//2*b*a^5*d^-1*(a+b*sin(c+d*x))^-1*(a^2+-1*b^2)^-4*cos(c+d*x))
@test integrate((a+b*sin(c+d*x))^-3*tan(c+d*x)^2, x) == :((1/2)*d^-1*(1+-1*sin(c+d*x))^-1*(a+b)^-3*cos(c+d*x)+-1//2*d^-1*(1+sin(c+d*x))^-1*(a+-1b)^-3*cos(c+d*x)+-1*a^2*d^-1*(a^2+-1*b^2)^-7//2*(b^2+2*a^2)*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x)))+-4*a^2*b^2*d^-1*(a^2+-1*b^2)^-7//2*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x)))+-2*b^2*d^-1*(a^2+-1*b^2)^-7//2*(b^2+3*a^2)*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x)))+-2*a*b^3*d^-1*(a+b*sin(c+d*x))^-1*(a^2+-1*b^2)^-3*cos(c+d*x)+-3//2*b*a^3*d^-1*(a+b*sin(c+d*x))^-1*(a^2+-1*b^2)^-3*cos(c+d*x)+-1//2*b*a^2*d^-1*(a+b*sin(c+d*x))^-2*(a^2+-1*b^2)^-2*cos(c+d*x))
@test integrate((a+b*sin(c+d*x))^-3*cot(c+d*x)^2, x) == :((1/2)*a^-1*d^-1*(a+b*sin(c+d*x))^-2*cot(c+d*x)+3*b*a^-4*d^-1*arctanh(cos(c+d*x))+-1*a^-4*d^-1*(a^2+-1*b^2)^-3//2*(2*a^4+6*b^4+-9*a^2*b^2)*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x)))+-1//2*a^-3*d^-1*(a^2+-1*b^2)^-1*(-6*b^2+5*a^2)*cot(c+d*x)+(1/2)*a^-2*d^-1*(a+b*sin(c+d*x))^-1*(a^2+-1*b^2)^-1*(-3*b^2+2*a^2)*cot(c+d*x))
@test integrate((a+b*sin(c+d*x))^-3*cot(c+d*x)^4, x) == :(1//6*a^-5*d^-1*(-60*b^2+17*a^2)*cot(c+d*x)+a^-6*d^-1*(a^2+-1*b^2)^-1//2*(2*a^4+20*b^4+-19*a^2*b^2)*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x)))+-1//2*b*a^-6*d^-1*(-20*b^2+9*a^2)*arctanh(cos(c+d*x))+-1//3*a^-1*d^-1*(a+b*sin(c+d*x))^-2*csc(c+d*x)^2*cot(c+d*x)+-1*a^-4*b^-1*d^-1*(a^2+-5*b^2)*cot(c+d*x)*csc(c+d*x)+1//6*a^-3*b^-1*d^-1*(a+b*sin(c+d*x))^-1*(-20*b^2+3*a^2)*cot(c+d*x)*csc(c+d*x)+1//6*a^-2*b^-1*d^-1*(a+b*sin(c+d*x))^-2*(-5*b^2+3*a^2)*cot(c+d*x)*csc(c+d*x))
@test integrate((a+b*sin(c+d*x))^-3*cot(c+d*x)^6, x) == :(-1//30*a^-7*d^-1*(91*a^4+630*b^4+-645*a^2*b^2)*cot(c+d*x)+-1*a^-8*d^-1*(a^2+-1*b^2)^(1/2)*(2*a^4+42*b^4+-29*a^2*b^2)*arctan((a^2+-1*b^2)^-1//2*(b+a*tan((1/2)*c+(1/2)*d*x)))+-1//3*b^-1*d^-1*(a+b*sin(c+d*x))^-2*cot(c+d*x)*csc(c+d*x)+-1//5*a^-1*d^-1*(a+b*sin(c+d*x))^-2*csc(c+d*x)^4*cot(c+d*x)+1//8*b*a^-8*d^-1*(45*a^4+168*b^4+-200*a^2*b^2)*arctanh(cos(c+d*x))+-1//30*a^-5*b^-2*d^-1*csc(c+d*x)^2*(15*a^4+210*b^4+-187*a^2*b^2)*cot(c+d*x)+1//8*a^-6*b^-1*d^-1*(8*a^4+84*b^4+-79*a^2*b^2)*cot(c+d*x)*csc(c+d*x)+1//12*a*b^-2*d^-1*(a+b*sin(c+d*x))^-2*csc(c+d*x)^2*cot(c+d*x)+7//20*b*a^-2*d^-1*(a+b*sin(c+d*x))^-2*csc(c+d*x)^3*cot(c+d*x)+1//12*a^-4*b^-2*d^-1*(a+b*sin(c+d*x))^-1*csc(c+d*x)^2*(4*a^4+63*b^4+-54*a^2*b^2)*cot(c+d*x)+1//60*a^-3*b^-2*d^-1*(a+b*sin(c+d*x))^-2*csc(c+d*x)^2*(5*a^4+63*b^4+-60*a^2*b^2)*cot(c+d*x))

max_stars_repo_path: test/4 Trig functions/4.1 Sine/4.1.1.3 (g tan)^p (a+b sin)^m.jl


text:
using CorpusLoaders
using Test
using Base.Iterators
using MultiResolutionIterators
using InternedStrings

if haskey(ENV, "CI") && ENV["CI"] == "true"
    @warn "WikiCorpus tests disabled on CI"
else
    # Do the tests
    @testset "basic use" begin
        wk_gen = load(WikiCorpus())
        docs = collect(take(wk_gen, 5));

        @test typeof(docs) ==
        Vector{CorpusLoaders.Document{@NestedVector(String, 4), String}}
        #    @test all(isa.(words, AbstractString))
        #    @test "a" ∈ words
        #    @test "the" ∈ words
        words = collect(flatten_levels(docs, (!lvls)(WikiCorpus, :word)))
        @test typeof(words) == Vector{String}
        @test length(words) >= length(docs)

        vocab = Set(words)
        @test "the" ∈ vocab
        @test "a" ∈ vocab

        #Make sure we're not catching anything we shouldn't
        @test "ENDOFARTICLE" ∉ vocab
        @test "ENDOFARTICLE." ∉ vocab
        @test "<doc" ∉ Set([x[1:min(end, 4)] for x in vocab])

        docs_of_words = full_consolidate(flatten_levels(docs, (!lvls)(WikiCorpus, :doc, :word)))
        @test typeof(docs_of_words) == Vector{CorpusLoaders.Document{Vector{String}, String}}
        @test length(docs_of_words) == length(docs)
        @test sum(length.(docs_of_words)) == length(words)
    end
end


text:
mutable struct Friction{T} <: Bound{T}
    Nx::Adjoint{T,SVector{6,T}}
    D::SMatrix{2,6,T,12}
    cf::T
    b::SVector{2,T}
    offset::SVector{6,T}

    function Friction(body::Body{T}, normal::AbstractVector{T}, cf::T;offset::AbstractVector{T} = zeros(3)) where T
        @assert cf>0
        normal = normal / norm(normal)

        # Derived from plane equation a*v1 + b*v2 + distance*v3 = p - offset
        A = Array(svd(skew(normal)).V) # gives two plane vectors
        A[:,3] = normal # to ensure correct sign
        Ainv = inv(A)
        ainv3 = Ainv[3,:]
        Nx = [ainv3;0;0;0]'
        D = [A[:,1:2];zeros(3,2)]'
        offset = [offset;0;0;0]

        new{T}(Nx, D, cf, zeros(2),offset), body.id
    end
end

@inline function g(friction::Friction, body::Body, Δt, No)
    friction.Nx[SVector(1, 2, 3)]' * (getx3(body, Δt) - friction.offset[SVector(1, 2, 3)])
end

@inline ∂g∂pos(friction::Friction, No) = friction.Nx
@inline ∂g∂vel(friction::Friction, Δt, No) = friction.Nx * Δt

@inline function schurf(ineqc, friction::Friction, i, body::Body, μ, Δt, No)
    φ = g(friction, body, Δt, No)

    γ1 = ineqc.γ1[i]
    s1 = ineqc.s1[i]

    return friction.Nx' * (γ1 / s1 * φ - μ / s1)
end

@inline function schurD(ineqc, friction::Friction, i, body::Body, Δt)
    Nx = friction.Nx
    Nv = Δt * Nx

    γ1 = ineqc.γ1[i]
    s1 = ineqc.s1[i]

    return Nx' * γ1 / s1 * Nv
end

# Direct stuff
@inline function setFrictionForce!(mechanism, ineqc, friction::Friction, i, body::Body)
    No = mechanism.No
    cf = friction.cf
    γ1 = ineqc.γ1[i]
    D = friction.D

    B = D'*friction.b
    F = body.F[No] - B[SVector(1,2,3)]
    τ = body.τ[No] - B[SVector(4,5,6)]
    setForce!(body,F,τ,No)

    f = body.f
    v = body.s1
    body.s1 = @SVector zeros(6)
    dyn = dynamics(mechanism, body)
    body.s1 = v
    body.f = f

    b0 = D*dyn

    if norm(b0) > cf*γ1
        friction.b = b0/norm(b0)*cf*γ1
    else
        friction.b = b0
    end

    B = D'*friction.b
    F += B[SVector(1,2,3)]
    τ += B[SVector(4,5,6)]
    setForce!(body,F,τ,No)

    return
end

# # Smooth stuff
# @inline function setFrictionForce!(mechanism, ineqc, friction::Friction, i, body::Body)
#     Δt = mechanism.Δt
#     No = mechanism.No
#     M = getM(body)
#     v = body.s1
#     cf = friction.cf
#     γ1 = ineqc.γ1[i]
#     D = friction.D

#     B = D'*friction.b
#     F = body.F[No] - B[SVector(1,2,3)]
#     τ = body.τ[No] - B[SVector(4,5,6)]
#     setForce!(body,F,τ,No)

#     ψ = Δt*norm(D*v)

#     f = body.f
#     body.s1 = @SVector zeros(6)
#     dyn = D/M*dynamics(body,mechanism)*Δt^2
#     body.s1 = v
#     body.f = f

#     X = D/M*D' * Δt^2 + I*(ψ/(cf*γ1))

#     friction.b = X\dyn
#     B = D'*friction.b
#     F += B[SVector(1,2,3)]
#     τ += B[SVector(4,5,6)]
#     setForce!(body,F,τ,No)
#     return
# end


text:
using DelimitedFiles

function ReadData(fname::AbstractString)
    return readdlm(fname, ',', Int)
end

function UpdateData!(data::Array)
    for i = 1:length(data)
        if data[i] == 0
            data[i] = 6
            push!(data, 8)
        else
            data[i] = data[i] - 1
        end
    end
end

function Part1(data::Array)
    for i = 1:80
        UpdateData!(data)
    end
    return length(data)
end

function UpdateDDict!(ddict::Dict)
    oldDict = Dict([(i, ddict[i]) for i = 0:8])

    ddict[8] = oldDict[0]
    ddict[7] = oldDict[8]
    ddict[6] = oldDict[7] + oldDict[0]
    for i = 0:5
        ddict[i] = oldDict[i+1]
    end
end

function Part2(data::Array)
    # store data in a dictionary
    ddict = Dict([(i, count(x -> x==i, data)) for i = 0:8])

    for i = 1:256
        UpdateDDict!(ddict)
    end
    return sum(values(ddict))
end

data = [3, 4, 3, 1, 2]

data = vec(ReadData("../data/day6.txt"))

#println(data)

println(Part2(data))


text:
abstract type PSRKRuleModel <: MixingRule end

struct PSRKRule{γ} <: PSRKRuleModel
    components::Array{String,1}
    activity::γ
end

@registermodel PSRKRule
export PSRKRule
function PSRKRule(components::Vector{String}; activity = Wilson, userlocations::Vector{String}=String[],activity_userlocations::Vector{String}=String[], verbose::Bool=false)
    init_activity = activity(components;userlocations = activity_userlocations,verbose)

    model = PSRKRule(components, init_activity)
    return model
end

function mixing_rule(model::RKModel,V,T,z,mixing_model::PSRKRuleModel,α,a,b,c)
    n = sum(z)
    invn = (one(n)/n)
    invn2 = invn^2
    g_E = excess_gibbs_free_energy(mixing_model.activity,1e5,T,z) / n
    b̄ = dot(z,Symmetric(b),z) * invn2
    c̄ = dot(z,c)/n
    Σlogb = sum(z[i]*log(b̄/b[i,i]) for i ∈ @comps)*invn
    Σab = sum(z[i]*a[i,i]*α[i]/b[i,i]/(R̄*T) for i ∈ @comps)*invn
    ā = b̄*R̄*T*(Σab-1/0.64663*(g_E/(R̄*T)+Σlogb))
    return ā,b̄,c̄
end


text:
Cassette.@context LensCtx

lens(nm, arg) = arg
lens(nm, args...) = args
lens(nm, args::NamedTuple) = args
lens(nm; kwargs...) = lens(nm, kwargs)

apl(f, args...) = f(args...)

function apl(fs::Tuple, args...)
  foreach(f -> f(args...), fs)
end

function Cassette.prehook(ctx::LensCtx, ::typeof(lens), nm::Symbol, args...)
  lmap = ctx.metadata
  if haskey(lmap, nm)
    fs = getfield(ctx.metadata, nm)
    apl(fs, args...)
  else
    args
  end
end

"Cassette based Lenscall f(args...)"
function lenscall(lmap::NamedTuple, f, args...; kwargs...)
  f_(args...) = f(args...; kwargs...)
  Cassette.overdub(LensCtx(metadata = lmap), f_, args...)
end

"Cassette based Lenscall f(args...)"
function lenscallnk(lmap::NamedTuple, f, arg)
  Cassette.overdub(LensCtx(metadata = lmap), f, arg)
end

# "Lens call with keyword args"
# function lenscall(lmap::NamedTuple, f, args...; kwargs...)
#   f_(args...) = f(args...; kwargs...)
#   lenscall(lmap, f_, args...)
# end


text:
using Documenter
using Weave
using NeXLMatrixCorrection

include("../weave/buildweave.jl")

pages = [
            "Example" => "example.md",
            "Coating Thickness" => "coatingthickness.md",
            "Compare To PAP" => "testagainstpap.md",
            "Compare to Heinrich" => "testagainstheinrich.md"
        ]

makedocs(modules = [NeXLMatrixCorrection], sitename = "NeXLMatrixCorrection.jl", pages = pages)

function addNISTHeaders(htmlfile::String)
    # read HTML
    html = transcode(String,read(htmlfile))
    # Find </head>
    i = findfirst(r"</[Hh][Ee][Aa][Dd]>", html)
    # Already added???
    j = findfirst("nist-header-footer", html)
    if isnothing(j) && (!isnothing(i))
        # Insert nist-pages links right before </head>
        res = html[1:i.start-1]*
            "<link rel=\"stylesheet\" href=\"https://pages.nist.gov/nist-header-footer/css/nist-combined.css\">\n"*
            "<script src=\"https://pages.nist.gov/nist-header-footer/js/jquery-1.9.0.min.js\" type=\"text/javascript\" defer=\"defer\"></script>\n"*
            "<script src=\"https://pages.nist.gov/nist-header-footer/js/nist-header-footer.js\" type=\"text/javascript\" defer=\"defer\"></script>\n"*
            html[i.start:end]
        write(htmlfile, res)
        println("Inserting NIST header/footer into $htmlfile")
    end
    return htmlfile
end

addNISTHeaders(joinpath(@__DIR__, "build","index.html"))
addNISTHeaders.(map(name->joinpath(@__DIR__, "build", splitext(name)[1], "index.html"), map(p->p.second, pages)))

# deploydocs(repo = "github.com/NeXLMatrixCorrection.jl.git")

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
include("testKNN.jl")
include("testSpMatBuilder.jl");
include("testRotation.jl");
include("testParamLevelSet.jl");
include("testParamLevelSet_extended.jl");
include("testMeshFreeParamLevelSet.jl");
include("testMeshFreeParamLevelSet_extended.jl");
include("testRegularization.jl");
include("testRotationParamLevelSet.jl");
include("testRotationParamLevelSet_extended.jl");


text:
# This file is a part of MGVI.jl, licensed under the MIT License (MIT).

function fisher_information(dist::Normal)
    _, σval = params(dist)
    inv_σ = inv(σval)
    inv_σ_2 = inv_σ * inv_σ
    res = _svector((inv_σ_2, 2*inv_σ_2))
    PDLinMapWithChol(Diagonal(res))
end

function fisher_information(dist::MvNormal)
    μval, Σ = params(dist)
    σval = without_chol(Σ)
    invσ = inv(σval)

    μdof = length(μval)
    σdof = μdof*(μdof+1) ÷ 2

    T = promote_type(eltype(μval), eltype(σval))
    covpart = UpperTriangular(fill(T(0), σdof, σdof))
    for my in 1:μdof
        m_flat_base = my*(my-1)÷2
        for mx in 1:my
            m_flat = m_flat_base + mx
            for ny in 1:my
                n_flat_base = ny*(ny-1)÷2
                for nx in 1:min(ny, (m_flat - n_flat_base))
                    n_flat = n_flat_base + nx
                    covpart[n_flat, m_flat] = (invσ[ny, mx] * invσ[my, nx] + invσ[mx, nx] * invσ[my, ny])
                end
            end
        end
    end

    for x in 1:μdof
        xflat = x*(x+1)÷2
        covpart[1:xflat, xflat] ./= 2
        covpart[xflat, xflat:end] ./= 2
    end

    sqrt_meanpart = cholesky(PositiveFactorizations.Positive, invσ).L
    meanpart_map = PDLinMapWithChol(invσ, sqrt_meanpart)

    sym_covpart = Symmetric(covpart)
    sqrt_covpart = cholesky(PositiveFactorizations.Positive, sym_covpart).L
    covpart_map = PDLinMapWithChol(sym_covpart, sqrt_covpart)

    blockdiag(meanpart_map, covpart_map)
end

function fisher_information(dist::MvNormal{<:Real,<:PDiagMat})
    Σ⁻¹ = inv(Diagonal(dist.Σ))
    mean_fisher_map = PDLinMapWithChol(Σ⁻¹)
    cov_fisher_map = PDLinMapWithChol(2*Σ⁻¹)
    blockdiag(mean_fisher_map, cov_fisher_map)
end

function fisher_information(dist::Exponential)
    λ = params(dist)[1]
    inv_l = inv(λ)
    res = _svector((inv_l * inv_l,))
    PDLinMapWithChol(Diagonal(res))
end

function fisher_information(dist::Poisson)
    λ = params(dist)[1]
    res = _svector((inv(λ),))
    PDLinMapWithChol(Diagonal(res))
end

function fisher_information(dist::Product)
    dists = dist.v
    λinformations = fisher_information.(dists)
    _blockdiag(λinformations)
end

function fisher_information(d::NamedTupleDist)
    dists = values(d)
    λinformations = map(fisher_information, dists)
    _blockdiag(λinformations)
end


text:
module InteractiveControlPlots

using ControlSystems, GtkInteract, Immerse
import Colors, Reactive, Measures, Compose

export bode_immerse, nyquist_immerse, pzmap_immerse, timeplot_immerse
export labelf, null_theme
export example_minimal, example_fo_sys

# package code goes here
include("framework.jl")
include("plots.jl")

examples = ["minimal", "fo_sys"]
for ex in examples
    funcname = Symbol("example_$(ex)")
    @eval begin
        function $funcname()
          pdir = joinpath(Pkg.dir("InteractiveControlPlots"),"src")
          include(joinpath(Pkg.dir("InteractiveControlPlots"),"src/examples/"*$ex*".jl"))
        end
    end
end
# function example_minimal()
#     include(joinpath(pdir,"examples/minimal.jl"))
# end
# function example_fo_sys()
#     include(joinpath(pdir,"examples/fo_sys.jl"))
# end

end # module


text:
"""
Stiffness problem imported from a .inp file.
"""
struct InpStiffness{dim, N, TF, M, TI, TBool, GO, TInds <: AbstractVector{TI}, TMeta<:Metadata} <: StiffnessTopOptProblem{dim, TF}
    inp_content::InpContent{dim, TF, N, TI}
    geom_order::Type{Val{GO}}
    ch::ConstraintHandler{DofHandler{dim, N, TF, M}, TF}
    black::TBool
    white::TBool
    varind::TInds
    metadata::TMeta
end

"""
Imports stiffness problem from a .inp file.
"""
function InpStiffness(filepath_with_ext::AbstractString)
    problem = Parser.extract_inp(filepath_with_ext)
    return InpStiffness(problem)
end
function InpStiffness(problem::Parser.InpContent)
    ch = Parser.inp_to_juafem(problem)
    black, white = find_black_and_white(ch.dh)
    varind = find_varind(black, white)
    metadata = Metadata(ch.dh)
    geom_order = JuAFEM.getorder(ch.dh.field_interpolations[1])
    return InpStiffness(problem, Val{geom_order}, ch, black, white, varind, metadata)
end

getE(p::InpStiffness) = p.inp_content.E
getν(p::InpStiffness) = p..inp_content.ν
nnodespercell(::InpStiffness{dim, N}) where {dim, N} = N
getgeomorder(p::InpStiffness{dim, N, TF, M, TI, GO}) where {dim, N, TF, M, TI, GO} = GO
getdensity(p::InpStiffness) = p.inp_content.density
getpressuredict(p::InpStiffness) = p.inp_content.dloads
getcloaddict(p::InpStiffness) = p.inp_content.cloads
getfacesets(p::InpStiffness) = p.inp_content.facesets


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct ExpressRouteCircuitPeeringsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Creates or updates a peering in the specified express route circuits.
Param: resourceGroupName::String (required)
Param: circuitName::String (required)
Param: peeringName::String (required)
Param: peeringParameters::ExpressRouteCircuitPeering (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ExpressRouteCircuitPeering
"""
function expressRouteCircuitPeeringsCreateOrUpdate(_api::ExpressRouteCircuitPeeringsApi, resourceGroupName::String, circuitName::String, peeringName::String, peeringParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", ExpressRouteCircuitPeering, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/expressRouteCircuits/{circuitName}/peerings/{peeringName}", ["azure_auth"], peeringParameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "circuitName", circuitName)  # type String
    Swagger.set_param(_ctx.path, "peeringName", peeringName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Deletes the specified peering from the specified express route circuit.
Param: resourceGroupName::String (required)
Param: circuitName::String (required)
Param: peeringName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function expressRouteCircuitPeeringsDelete(_api::ExpressRouteCircuitPeeringsApi, resourceGroupName::String, circuitName::String, peeringName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/expressRouteCircuits/{circuitName}/peerings/{peeringName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "circuitName", circuitName)  # type String
    Swagger.set_param(_ctx.path, "peeringName", peeringName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets the specified authorization from the specified express route circuit.
Param: resourceGroupName::String (required)
Param: circuitName::String (required)
Param: peeringName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ExpressRouteCircuitPeering
"""
function expressRouteCircuitPeeringsGet(_api::ExpressRouteCircuitPeeringsApi, resourceGroupName::String, circuitName::String, peeringName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ExpressRouteCircuitPeering, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/expressRouteCircuits/{circuitName}/peerings/{peeringName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "circuitName", circuitName)  # type String
    Swagger.set_param(_ctx.path, "peeringName", peeringName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets all peerings in a specified express route circuit.
Param: resourceGroupName::String (required)
Param: circuitName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ExpressRouteCircuitPeeringListResult
"""
function expressRouteCircuitPeeringsList(_api::ExpressRouteCircuitPeeringsApi, resourceGroupName::String, circuitName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ExpressRouteCircuitPeeringListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/expressRouteCircuits/{circuitName}/peerings", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "circuitName", circuitName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export expressRouteCircuitPeeringsCreateOrUpdate, expressRouteCircuitPeeringsDelete, expressRouteCircuitPeeringsGet, expressRouteCircuitPeeringsList