text:
# black box, uses score function estimator
function single_sample_gradient_estimate!(
        var_model::GenerativeFunction, var_model_args::Tuple,
        model::GenerativeFunction, model_args::Tuple, observations::ChoiceMap,
        scale_factor=1.)

    # sample from variational approximation
    var_trace = simulate(var_model, var_model_args)

    # compute learning signal
    constraints = merge(observations, get_choices(var_trace)) # TODO characterize what it means when not all var_model choices are in the model..
    (model_trace, model_log_weight) = generate(model, model_args, constraints)
    log_weight = model_log_weight - get_score(var_trace)

    # accumulate the weighted gradient
    accumulate_param_gradients!(var_trace, nothing, log_weight * scale_factor)

    # unbiased estimate of objective function, and trace
    (log_weight, var_trace, model_trace)
end

function vimco_geometric_baselines(log_weights)
    num_samples = length(log_weights)
    s = sum(log_weights)
    baselines = Vector{Float64}(undef, num_samples)
    for i=1:num_samples
        temp = log_weights[i]
        log_weights[i] = (s - log_weights[i]) / (num_samples - 1)
        baselines[i] = logsumexp(log_weights) - log(num_samples)
        log_weights[i] = temp
    end
    baselines
end

function logdiffexp(x, y)
    m = max(x, y)
    m + log(exp(x - m) - exp(y - m))
end

function vimco_arithmetic_baselines(log_weights)
    num_samples = length(log_weights)
    log_total_weight = logsumexp(log_weights)
    baselines = Vector{Float64}(undef, num_samples)
    for i=1:num_samples
        log_sum_f_without_i = logdiffexp(log_total_weight, log_weights[i])
        log_f_hat = log_sum_f_without_i - log(num_samples - 1)
        baselines[i] = logsumexp(log_sum_f_without_i, log_f_hat) - log(num_samples)
    end
    baselines
end

# black box, VIMCO gradient estimator
# for use in training models
function multi_sample_gradient_estimate!(
        var_model::GenerativeFunction, var_model_args::Tuple,
        model::GenerativeFunction, model_args::Tuple, observations::ChoiceMap,
        num_samples::Int, scale_factor=1., geometric=true)

    # sample from variational approximation multiple times
    traces = Vector{Any}(undef, num_samples)
    log_weights = Vector{Float64}(undef, num_samples)
    for i=1:num_samples
        traces[i] = simulate(var_model, var_model_args)
        constraints = merge(observations, get_choices(traces[i])) # TODO characterize as above
        model_weight, = assess(model, model_args, constraints)
        log_weights[i] = model_weight - get_score(traces[i])
    end

    # multi-sample log marginal likelihood estimate
    log_total_weight = logsumexp(log_weights)
    L = log_total_weight - log(num_samples)

    # baselines
    if geometric
        baselines = vimco_geometric_baselines(log_weights)
    else
        baselines = vimco_arithmetic_baselines(log_weights)
    end

    weights_normalized = exp.(log_weights .- log_total_weight)
    for i=1:num_samples
        learning_signal = (L - baselines[i]) - weights_normalized[i]
        accumulate_param_gradients!(traces[i], nothing, learning_signal * scale_factor)
    end

    # collection of traces and normalized importance weights, and estimate of
    # objective function
    (L, traces, weights_normalized)
end

function _maybe_accumulate_param_grad!(trace, update::ParamUpdate, scale_factor::Real)
    return accumulate_param_gradients!(trace, nothing, scale_factor)
end

function _maybe_accumulate_param_grad!(trace, update::Nothing, scale_factor::Real)
end

"""
    (elbo_estimate, traces, elbo_history) = black_box_vi!(
        model::GenerativeFunction, model_args::Tuple,
        [model_update::ParamUpdate,]
        observations::ChoiceMap,
        var_model::GenerativeFunction, var_model_args::Tuple,
        var_model_update::ParamUpdate;
        options...)

Fit the parameters of a variational model (`var_model`) to the posterior
distribution implied by the given `model` and `observations` using stochastic
gradient methods. Users may optionally specify a `model_update` to jointly
update the parameters of `model`.

# Additional arguments:
- `iters=1000`: Number of iterations of gradient descent.
- `samples_per_iter=100`: Number of samples from the variational and generative
        model to accumulate gradients over before a single gradient step.
- `verbose=false`: If `true`, print information about the progress of fitting.
- `callback`: Callback function that takes `(iter, traces, elbo_estimate)`
        as input, where `iter` is the iteration number and `traces` are samples
        from `var_model` for that iteration.
"""
function black_box_vi!(
        model::GenerativeFunction, model_args::Tuple,
        model_update::Union{ParamUpdate,Nothing},
        observations::ChoiceMap,
        var_model::GenerativeFunction, var_model_args::Tuple,
        var_model_update::ParamUpdate;
        iters=1000, samples_per_iter=100, verbose=false,
        callback=(iter, traces, elbo_estimate) -> nothing)

    var_traces = Vector{Any}(undef, samples_per_iter)
    model_traces = Vector{Any}(undef, samples_per_iter)
    elbo_history = Vector{Float64}(undef, iters)
    for iter=1:iters

        # compute gradient estimate and objective function estimate
        elbo_estimate = 0.0
        # TODO multithread (note that this would require accumulate_param_gradients! to be threadsafe)
        for sample=1:samples_per_iter

            # accumulate the variational family gradients
            (log_weight, var_trace, model_trace) = single_sample_gradient_estimate!(
                var_model, var_model_args,
                model, model_args, observations, 1/samples_per_iter)
            elbo_estimate += (log_weight / samples_per_iter)

            # accumulate the generative model gradients
            _maybe_accumulate_param_grad!(model_trace, model_update, 1.0 / samples_per_iter)

            # record the traces
            var_traces[sample] = var_trace
            model_traces[sample] = model_trace
        end
        elbo_history[iter] = elbo_estimate

        # print it
        verbose && println("iter $iter; est objective: $elbo_estimate")

        # callback
        callback(iter, var_traces, elbo_estimate)

        # update parameters of variational family
        apply!(var_model_update)

        # update parameters of generative model
        if !isnothing(model_update)
            apply!(model_update)
        end
    end

    (elbo_history[end], var_traces, elbo_history, model_traces)
end

black_box_vi!(model::GenerativeFunction, model_args::Tuple,
              observations::ChoiceMap,
              var_model::GenerativeFunction, var_model_args::Tuple,
              var_model_update::ParamUpdate; options...) =
    black_box_vi!(model, model_args, nothing, observations,
                  var_model, var_model_args, var_model_update; options...)

"""
    (iwelbo_estimate, traces, iwelbo_history) = black_box_vimco!(
        model::GenerativeFunction, model_args::Tuple,
        [model_update::ParamUpdate,]
        observations::ChoiceMap,
        var_model::GenerativeFunction, var_model_args::Tuple,
        var_model_update::ParamUpdate,
        grad_est_samples::Int; options...)

Fit the parameters of a variational model (`var_model`) to the posterior
distribution implied by the given `model` and `observations` using stochastic
gradient methods applied to the [Variational Inference with Monte Carlo
Objectives](https://arxiv.org/abs/1602.06725) (VIMCO) lower bound on the
marginal likelihood. Users may optionally specify a `model_update` to jointly
update the parameters of `model`.

# Additional arguments:
- `grad_est_samples::Int`: Number of samples for the VIMCO gradient estimate.
- `iters=1000`: Number of iterations of gradient descent.
- `samples_per_iter=100`: Number of samples from the variational and generative
        model to accumulate gradients over before a single gradient step.
- `geometric=true`: Whether to use the geometric or arithmetric baselines
        described in [Variational Inference with Monte Carlo
        Objectives](https://arxiv.org/abs/1602.06725)
- `verbose=false`: If `true`, print information about the progress of fitting.
- `callback`: Callback function that takes `(iter, traces, elbo_estimate)`
        as input, where `iter` is the iteration number and `traces` are samples
        from `var_model` for that iteration.
"""
function black_box_vimco!(
        model::GenerativeFunction, model_args::Tuple,
        model_update::Union{ParamUpdate,Nothing}, observations::ChoiceMap,
        var_model::GenerativeFunction, var_model_args::Tuple,
        var_model_update::ParamUpdate, grad_est_samples::Int;
        iters=1000, samples_per_iter=100, geometric=true, verbose=false,
        callback=(iter, traces, elbo_estimate) -> nothing)

    resampled_var_traces = Vector{Any}(undef, samples_per_iter)
    model_traces = Vector{Any}(undef, samples_per_iter)

    iwelbo_history = Vector{Float64}(undef, iters)
    for iter=1:iters

        # compute gradient estimate and objective function estimate
        iwelbo_estimate = 0.
        for sample=1:samples_per_iter

            # accumulate the variational family gradients
            (est, original_var_traces, weights) = multi_sample_gradient_estimate!(
                var_model, var_model_args,
                model, model_args, observations, grad_est_samples,
                1/samples_per_iter, geometric)
            iwelbo_estimate += (est / samples_per_iter)

            # record a variational trace obtained by resampling from the weighted collection
            resampled_var_traces[sample] = original_var_traces[categorical(weights)]

            # accumulate the generative model gradient estimator
            for (var_trace, weight) in zip(original_var_traces, weights)
                constraints = merge(observations, get_choices(var_trace))
                (model_trace, _) = generate(model, model_args, constraints)
                _maybe_accumulate_param_grad!(model_trace, model_update, weight / samples_per_iter)
            end
        end
        iwelbo_history[iter] = iwelbo_estimate

        # print it
        verbose && println("iter $iter; est objective: $iwelbo_estimate")

        # callback
        callback(iter, resampled_var_traces, iwelbo_estimate)

        # update parameters of variational family
        apply!(var_model_update)

        # update parameters of generative model
        if !isnothing(model_update)
            apply!(model_update)
        end

    end

    (iwelbo_history[end], resampled_var_traces, iwelbo_history, model_traces)
end

black_box_vimco!(model::GenerativeFunction, model_args::Tuple,
                 observations::ChoiceMap,
                 var_model::GenerativeFunction, var_model_args::Tuple,
                 var_model_update::ParamUpdate,
                 grad_est_samples::Int; options...) =
    black_box_vimco!(model, model_args, nothing, observations,
                     var_model, var_model_args, var_model_update,
                     grad_est_samples; options...)

export black_box_vi!, black_box_vimco!


text:
# Convert a fasta formatted file of sequences to a csv file.

using CSV
using DataFrames
using BioSequences

df = DataFrame(id=[], seq=[])
reader = FASTA.Reader(open("/home/dillon/data/cUR50/cur50.fasta", "r"))
count = 0
for (i, record) in enumerate(reader)
    try
        if FASTA.hasidentifier(record) && FASTA.hassequence(record)
            push!(df, [FASTA.identifier(record), FASTA.sequence(record)])
        end
        if i % 10^5 == 0
            print("$i read\n")
        end
    catch exc
        print("Exception encountered: ", exc)
        print(record)
    end
end
close(reader)
print("Writing csv.")
CSV.write("/home/dillon/data/cUR50/cur50.csv", df)

ext: jl


text:
module DFTranspose

using DataFrames

export df_transpose

include("transpose.jl")

end


text:
# Autogenerated wrapper script for Antlr4CppRuntime_jll for armv6l-linux-musleabihf-cxx03
export libantlr4_cpp_runtime

using Libuuid_jll
JLLWrappers.@generate_wrapper_header("Antlr4CppRuntime")
JLLWrappers.@declare_library_product(libantlr4_cpp_runtime, "libantlr4-runtime.so.4.9.3")
function __init__()
    JLLWrappers.@generate_init_header(Libuuid_jll)
    JLLWrappers.@init_library_product(
        libantlr4_cpp_runtime,
        "lib/libantlr4-runtime.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "context_implementations.jl" begin
    # https://github.com/TuringLang/DynamicPPL.jl/issues/129
    @testset "#129" begin
        @model function test(x)
            μ ~ MvNormal(fill(0, 2), 2.0)
            z = Vector{Int}(undef, length(x))
            # `z .~ Categorical.(ps)` cannot be parsed by Julia 1.0
            (.~)(z, Categorical.(fill([0.5, 0.5], length(x))))
            for i in 1:length(x)
                x[i] ~ Normal(μ[z[i]], 0.1)
            end
        end

        test([1, 1, -1])(VarInfo(), SampleFromPrior(), LikelihoodContext())
    end

    # https://github.com/TuringLang/DynamicPPL.jl/issues/28#issuecomment-829223577
    @testset "dot tilde: arrays of distributions" begin
        @testset "assume" begin
            @model function test(x, size)
                y = Array{Float64,length(size)}(undef, size...)
                y .~ Normal.(x)
                return y, getlogp(__varinfo__)
            end

            for ysize in ((2,), (2, 3), (2, 3, 4))
                for x in (
                    # scalar
                    randn(),
                    # drop trailing dimensions
                    ntuple(i -> randn(ysize[1:i]), length(ysize))...,
                    # singleton dimensions
                    ntuple(
                        i -> randn(ysize[1:(i - 1)]..., 1, ysize[(i + 1):end]...),
                        length(ysize),
                    )...,
                )
                    model = test(x, ysize)
                    y, lp = model()
                    @test lp ≈ sum(logpdf.(Normal.(x), y))

                    ys = [first(model()) for _ in 1:10_000]
                    @test norm(mean(ys) .- x, Inf) < 0.1
                    @test norm(std(ys) .- 1, Inf) < 0.1
                end
            end
        end

        @testset "observe" begin
            @model function test(x, y)
                return y .~ Normal.(x)
            end

            for ysize in ((2,), (2, 3), (2, 3, 4))
                for x in (
                    # scalar
                    randn(),
                    # drop trailing dimensions
                    ntuple(i -> randn(ysize[1:i]), length(ysize))...,
                    # singleton dimensions
                    ntuple(
                        i -> randn(ysize[1:(i - 1)]..., 1, ysize[(i + 1):end]...),
                        length(ysize),
                    )...,
                )
                    y = randn(ysize)
                    z = logjoint(test(x, y), VarInfo())
                    @test z ≈ sum(logpdf.(Normal.(x), y))
                end
            end
        end
    end
end


text:
# NOTE: Import libraries here ...

# DP GMM model under stick-breaking construction
@model dp_gmm_sb(y, K) = begin
    nobs = length(y)

    mu ~ filldist(Normal(0, 3), K)
    sig ~ filldist(Gamma(1, 1/10), K)  # mean = 0.1

    alpha ~ Gamma(1, 1/10)  # mean = 0.1
    crm = DirichletProcess(alpha)
    v ~ filldist(StickBreakingProcess(crm), K - 1)
    eta = stickbreak(v)

    y .~ UnivariateGMM(mu, sig, Categorical(eta))
end

# NOTE: Read data y here ...
# Here, y (a vector of length 500) is noisy univariate draws from a
# mixture distribution with 4 components.

# Fit DP-SB-GMM with ADVI
advi = ADVI(1, 2000)  # num_elbo_samples, max_iters
q = vi(dp_gmm_sb(y, 10), advi, optimizer=Flux.ADAM(1e-2));

# Fit DP-SB-GMM with HMC
hmc_chain = sample(dp_gmm_sb(y, 500),  # data, number of mixture components
                   HMC(0.01, 100),  # stepsize, number of leapfrog steps
                   1000)  # iterations

# Fit DP-SB-GMM with NUTS
@time nuts_chain = begin
    n_samples = 500  # number of MCMC samples
    nadapt = 500  # number of iterations to adapt tuning parameters in NUTS
    iterations = n_samples + nadapt
    target_accept_ratio = 0.8

    sample(dp_gmm_sb(y, 10),  # data, number of mixture components.
           NUTS(nadapt, target_accept_ratio, max_depth=10),
           iterations);
end

max_stars_repo_name: luiarthur/TuringBnpBenchmarks

max_issues_repo_name: luiarthur/TuringBnpBenchmarks

max_forks_repo_name: luiarthur/TuringBnpBenchmarks


text:
"""
    Proof of concept for Caesar-ROS integration
    (check Caesar wiki for details/instructions)

    To do:
    - re-enable JSON replies
    - switch to a pattern that does not require use of global variables
    - periodic export of factor graph object
"""

using RobotOS

# standard types
@rosimport geometry_msgs.msg: Point, Pose2D
# bespoke types
@rosimport caesar_ros.srv: AddVariable, AddPrior3D, AddFactorZPR, AddFactorXYH

rostypegen()
using .geometry_msgs.msg
using .caesar_ros.srv

### Load Caesar mmisam stuff
using RoME
using DistributedFactorGraphs

using JSON2
using DocStringExtensions

global fg

"""
    $SIGNATURES

Service callback for adding a variable to the factor graph.

Notes
- See https://github.com/JuliaRobotics/DistributedFactorGraphs.jl/issues/134
"""
function add_variable(req::AddVariableRequest)
    @info("received add_variable request", req)
    global fg
    reply = AddVariableResponse("failed to add variable")
    type_symbol = Symbol(req.type)
    if in(type_symbol, [:Pose2, :Pose3, :Point3])
        global fg
        vars = ls(fg)
        id = Symbol(req.id)
        if(!in(id, vars))
            retv = addVariable!(fg, id, getfield(Main, type_symbol))
            reply.status = JSON2.write(packVariable(fg, retv))
        else
            @warn("Failed to add variable - repeat variable")
        end
    else
        @warn("Failed to add variable - unsupported type", type_symbol)
    end
    return(reply)
end

"""
    $SIGNATURES

Service callback for adding a binary XYH partial factor.
"""
function add_factor_xyh(req::AddFactorXYHRequest )
    @info("received add_factor_xyh request", req)
    id0=Symbol(req.id0)
    id1=Symbol(req.id1)
    global fg
    vars = ls(fg)
    reply = AddFactorXYHResponse()
    if (in(id0, vars) & in(id1, vars))
        mu_xy = [req.x, req.y]
        Sigma = reshape(req.covariance,3,3)
        Sigma_xy = Sigma[1:2,1:2]
        z_xyh = Pose3Pose3XYYaw(MvNormal(mu_xy, Sigma[1:2,1:2]), Normal(req.yaw, Sigma[3,3]))
        retv = addFactor!(fg, [id0, id1], z_xyh)
        # @show reply.status = JSON2.write(packVariable(fg, retv))
    else
        @warn("Failed to add FactorXYH - missing pose")
        reply.status="failed"
    end

    return(reply)
end

"""
    $SIGNATURES

Service callback for adding a unary  partial factor.
"""
function add_factor_zpr(req::AddFactorZPRRequest )
    @info("received add_factor_zpr request", req)
    id = Symbol(req.id)
    global fg
    vars = ls(fg)
    reply = AddFactorZPRResponse()
    if (in(id,vars) )
        mu_z = req.z
        mu_rp = [req.roll, req.pitch]
        Sigma_zpr = reshape(req.covariance, 3, 3)
        sigma_z = Sigma_zpr[1,1]
        Sigma_rp = Sigma_zpr[3:-1:2,3:-1:2]
        z_zpr = PriorPose3ZRP(Normal(mu_z, sigma_z), MvNormal(mu_rp, Sigma_rp))
        retv = addFactor!(fg, [id], z_zpr)
        # reply.status = JSON2.write(packVariable(fg, retv)) # breaks
    else
        @warn("Failed to add FactorZPR - missing pose")
        reply.status="failed"
    end

    return(reply)
end

function main()
    init_node("caesar_srv_endpoint")

    global fg
    fg = initfg()

    # services - variable
    add_variable_srv = Service("AddVariable", caesar_ros.srv.AddVariable, add_variable)
    # services - factors
    add_factor_zpr_srv = Service("AddFactorZPR",caesar_ros.srv.AddFactorZPR, add_factor_zpr)
    add_factor_xyh_srv = Service("AddFactorXYH",caesar_ros.srv.AddFactorXYH, add_factor_xyh )

    # main loop
    @info "services have been set up; entering main loop"
    loop_rate = Rate(5.0)
    while ! is_shutdown()
        rossleep(loop_rate)
    end
end

main()


text:
#=
Run test locally using
$ julia test/runtest.jl /dev/ttyXYZ

/dev/ttyXYZ can be:
/dev/ttyS0  (for Travis)
/dev/ttyS4
/dev/ttyUSB0
=#

using LibSerialPort
using Test

if haskey(ENV, "CI")
    @testset "LibSerialPort CI port listing" begin
        list_ports()
    end
else
    if length(ARGS) == 0
        port = "/dev/ttyS0"  # /dev/ttyS4 /dev/ttyUSB0
    else
        port = ARGS[1]
    end

    baudrate = length(ARGS) >= 2 ? ARGS[2] : "115200"

    @testset "LibSerialPort" begin
        @testset "Low level API" begin
            include("test-low-level-api.jl")
            @test test_low_level_api() == nothing
            @test test_low_level_api(port, baudrate) == nothing
        end

        @testset "High level API" begin
            include("test-high-level-api.jl")
            @test test_high_level_api() == nothing
            @test test_high_level_api(port, baudrate) == nothing
        end

        @testset "Reading with timeouts" begin
            LibSerialPort.open(port, 115200) do s
                #Tests assume serial port being tested won't output anything after being flushed
                #TODO: Find a better way to test this
                flush(s)
                @test readline(s, 1.0) == "" #readline with a 1 second timeout
                @test readuntil(s, 'a', 1.0) == "" #readuntil 'a' with a 1 second timeout
            end
        end

        # console.jl runs forever, thus isn't amenable to unit testing
        # @testset "Examples" begin
        #     include("../examples/console.jl")
        #     @test console() == nothing
        #     @test console(port, baudrate) == nothing
        # end
    end
end


text:
@testset "Table" begin
    t = @inferred(Table(a = [1, 2, 3], b = [2.0, 4.0, 6.0]))::Table

    @test Table(t) == t
    @test Table(t; c = [true,false,true]) == Table(a = [1,2,3], b = [2.0,4.0,6.0], c = [true,false,true])
    @test Table(t, Table(c = [true,false,true])) == Table(a = [1,2,3], b = [2.0,4.0,6.0], c = [true,false,true])
    @test Table(t; b = nothing) == Table(a = [1,2,3])

    @test columnnames(t) == (:a, :b)
    @test propertynames(t) == (:a, :b)
    @test columns(t) == (a = [1,2,3], b = [2.0, 4.0, 6.0])
    @test rows(t) === t
    @test Tables.istable(t) === true
    @test Tables.rowaccess(t) === true
    @test Tables.columnaccess(t) === true
    @test Tables.schema(t) == Tables.Schema(NamedTuple{(:a,:b),Tuple{Int,Float64}})
    @test @inferred(Tables.materializer(t)(map(x -> 2*x, Tables.columns(t)))) isa Table
    @test Tables.materializer(t)(map(x -> 2*x, Tables.columns(t))) == Table(a = [2,4,6], b = [4.0,8.0,12.0])

    @test t.a == [1,2,3]
    @test axes(t) === (Base.OneTo(3),)
    @test size(t) === (3,)
    @test length(t) === 3
    @test @inferred(t[2]) === (a = 2, b = 4.0)
    @test @inferred(t[1:2])::Table == Table(a = [1,2], b = [2.0,4.0])
    @test @inferred(t[:])::Table == t
    @test @inferred(view(t, 1:2))::Table == Table(a = [1,2], b = [2.0,4.0])
    @test @inferred(view(t, :))::Table == t
    @test_throws BoundsError t[4]

    @test similar(t) isa typeof(t)
    @test axes(similar(t)) === axes(t)

    @test @inferred(vcat(t))::Table == t
    @test @inferred(vcat(t, t))::Table == Table(a = [1,2,3,1,2,3], b = [2.0, 4.0, 6.0, 2.0, 4.0, 6.0])
    @test @inferred(vcat(t, t, t))::Table == Table(a = [1,2,3,1,2,3,1,2,3], b = [2.0, 4.0, 6.0, 2.0, 4.0, 6.0, 2.0, 4.0, 6.0])

    @test @inferred(hcat(t))::Table == Table(a = hcat([1; 2; 3]), b = hcat([2.0; 4.0; 6.0]))
    @test @inferred(hcat(t, t))::Table == Table(a = [1 1;2 2;3 3], b = [2.0 2.0; 4.0 4.0; 6.0 6.0])
    @test @inferred(hcat(t, t, t))::Table == Table(a = [1 1 1;2 2 2;3 3 3], b = [2.0 2.0 2.0; 4.0 4.0 4.0; 6.0 6.0 6.0])

    @test [t t; t t]::Table == Table(a = [1 1;2 2;3 3;1 1;2 2;3 3], b = [2.0 2.0; 4.0 4.0; 6.0 6.0; 2.0 2.0; 4.0 4.0; 6.0 6.0])
    @test @inferred(vec(t))::Table == t

    io = IOBuffer()
    show(io, MIME"text/plain"(), t)
    str = String(take!(io))
    @test str == """
        Table with 2 columns and 3 rows:
             a  b
           ┌───────
         1 │ 1  2.0
         2 │ 2  4.0
         3 │ 3  6.0"""

    @test @inferred(getproperties((:b, :a))(t))::Table == Table(b = [2.0, 4.0, 6.0], a = [1, 2, 3])

    @test t == t
    @test isequal(t, t)
    @test !isless(t, t)

    t2 = empty(t)
    @test t2 isa typeof(t)
    @test isempty(t2)

    t3 = copy(t)
    @test t3 isa typeof(t)
    @test t3 == t
    @test @inferred(t3[1:2])::Table == Table(a = [1,2], b = [2.0,4.0])
    t3[3] = (a = 4, b = 8.0)
    @test t3 == Table(a = [1,2,4], b = [2.0, 4.0, 8.0])
    t3[1:2] .= Ref((a = 0, b = 0.0))
    @test t3 == Table(a = [0,0,4], b = [0.0, 0.0, 8.0])
    t3[1:2] = Table(a = [2,3], b = [4.0, 6.0])
    @test t3 == Table(a = [2,3,4], b = [4.0, 6.0, 8.0])
    v3 = view(t3, :)
    v3[3] = (a = 5, b = 10.0)
    @test t3 == Table(a = [2,3,5], b = [4.0, 6.0, 10.0])

    empty!(t3)
    @test isempty(t3)

    push!(t3, (a=10, b=0.0))
    @test t3 == Table(a = [10], b = [0.0])

    pushfirst!(t3, (a=-10, b=-10.0))
    @test t3 == Table(a = [-10, 10], b = [-10.0, 0.0])

    insert!(t3, 2, (a=5,b=5.0))
    @test t3 == Table(a = [-10, 5, 10], b = [-10.0, 5.0, 0.0])

    @test splice!(t3, 2:1) == empty(t)
    @test t3 == Table(a = [-10, 5, 10], b = [-10.0, 5.0, 0.0])
    @test splice!(t3, 2:1, (a=1, b=1.0)) == empty(t)
    @test t3 == Table(a = [-10, 1, 5, 10], b = [-10.0, 1.0, 5.0, 0.0])
    @test splice!(t3, 3:2, Table(a=[2], b=[2.0])) == empty(t)
    @test t3 == Table(a = [-10, 1, 2, 5, 10], b = [-10.0, 1.0, 2.0, 5.0, 0.0])
    @test splice!(t3, 2) === (a = 1, b = 1.0)
    @test t3 == Table(a = [-10, 2, 5, 10], b = [-10.0, 2.0, 5.0, 0.0])
    @test splice!(t3, 2, (a = 3, b = 3.0)) === (a = 2, b = 2.0)
    @test t3 == Table(a = [-10, 3, 5, 10], b = [-10.0, 3.0, 5.0, 0.0])
    @test splice!(t3, 2, Table(a = [4], b = [4.0])) === (a = 3, b = 3.0)
    @test t3 == Table(a = [-10, 4, 5, 10], b = [-10.0, 4.0, 5.0, 0.0])

    @test pop!(t3) === (a = 10, b = 0.0)
    @test t3 == Table(a = [-10, 4, 5], b = [-10.0, 4.0, 5.0])

    @test popfirst!(t3) === (a = -10, b = -10.0)
    @test t3 == Table(a = [4, 5], b = [4.0, 5.0])

    append!(t3, Table(a = [6], b = [6.0]))
    @test t3 == Table(a = [4,5,6], b = [4.0, 5.0, 6.0])

    prepend!(t3, Table(a = [3], b = [3.0]))
    @test t3 == Table(a = [3,4,5,6], b = [3.0, 4.0, 5.0, 6.0])

    deleteat!(t3, 4)
    @test t3 == Table(a = [3,4,5], b = [3.0, 4.0, 5.0])

    @test length(resize!(Table(a=[]), 100)) == 100

    @test sort(Table(a=collect(100:-1:1))) == Table(a=1:100)

    @testset "Merging Tables" begin
        t1 = Table(a = [1,2,3],)
        t2 = Table(b = [2.0, 4.0, 6.0],)
        t3 = Table(a = [1,2,3], b = [2.0, 4.0, 6.0])

        @test @inferred(map(merge, t1, t2))::Table == t3
        @test @inferred(mapview(merge, t1, t2))::Table == t3
        @test @inferred(broadcast(merge, t1, t2))::Table == t3
    end

    @testset "GetProperty on Tables" begin
        t = Table(a = [1,2,3], b = [2.0, 4.0, 6.0], c = [true, false, true])

        @test @inferred(map(getproperty(:a), t))::Vector == [1,2,3]
        @test @inferred(mapview(getproperty(:a), t))::Vector == [1,2,3]
        @test @inferred(broadcast(getproperty(:a), t))::Vector == [1,2,3]
        @test @inferred(mapreduce(getproperty(:a), +, t)) === 6
        @test @inferred(filter(getproperty(:c), t))::Table == Table(a = [1,3], b = [2.0, 6.0], c = [true, true])
        @test @inferred(findall(getproperty(:c), t))::Vector{Int} == [1, 3]

        @test @inferred(map(getproperties((:a,)), t))::Table == Table(a = [1,2,3])
        @test @inferred(mapview(getproperties((:a,)), t))::Table == Table(a = [1,2,3])
        @test @inferred(broadcast(getproperties((:a,)), t))::Table == Table(a = [1,2,3])
        @test @inferred(mapreduce(getproperties((:a,)), (acc, row) -> acc + row.a, t; init = 0)) === 6
    end

    @testset "@Select and @Compute on Tables" begin
        t = Table(a = [1,2,3], b = [2.0, 4.0, 6.0], c = [true, false, true])

        c = @Compute(2*$a)
        @test @inferred(c(t))::Vector == [2, 4, 6] # (Works because 2 * vector works)
        @test @inferred(map(c, t))::Vector == [2, 4, 6]
        @test @inferred(mapview(c, t))::MappedArray == [2, 4, 6]
        @test @inferred(broadcast(c, t))::Vector == [2, 4, 6]
        @test @inferred(mapreduce(c, +, t)) === 12

        c2 = @Compute($b > 3.0)
        @test @inferred(filter(c2, t))::Table == Table(a = [2,3], b = [4.0,6.0], c = [false,true])
        @test @inferred(findall(c2, t))::Vector{Int} == [2, 3]

        s = @Select(sum = $a + $b)
        @test @inferred(s(t))::Table == Table(sum = [3.0, 6.0, 9.0]) # (Works because vector + vector works)
        @test @inferred(map(s, t))::Table == Table(sum = [3.0, 6.0, 9.0])
        @test @inferred(mapview(s, t))::Table == Table(sum = [3.0, 6.0, 9.0])
        @test @inferred(broadcast(s, t))::Table == Table(sum = [3.0, 6.0, 9.0])
        @test @inferred(mapreduce(s, (acc, row) -> acc + row.sum, t; init = 0.0)) === 18.0
    end

    @testset "missing in Tables" begin
        t = Table(a = [1, 2, 3], b = [2.0, 4.0, missing])

        @test t[1]::eltype(t) == (a = 1, b = 2.0)
        @test isequal(t[3]::eltype(t), (a = 3, b = missing))

        @test (t == t) === missing
        @test isequal(t, t)
        @test !isless(t, t)
    end

    @testset "Tables.jl" begin
        t = Table(a = [1, 2, 3], b = [2.0, 4.0, missing])
        @test isequal(t |> columntable, getfield(t, :data))
        r = t |> rowtable
        @test length(r) == 3
        for (a, b) in zip(t, r)
            @test isequal(a, b)
        end
        @test isequal(Table(t |> columntable), t)
        @test isequal(Table(t |> rowtable), t)
    end

    @testset "group" begin
        t = Table(a = [1, 2, 1], b = [2.0, 4.0, 6.0])
        out = group(getproperty(:a), t)
        @test typeof(out) <: SplitApplyCombine.AbstractDictionary{Int, <:Table}
        @test out == SplitApplyCombine.dictionary([
            1 => Table(a=[1, 1], b=[2.0, 6.0]),
            2 => Table(a=[2],    b=[4.0])
        ])
    end

    @testset "innerjoin" begin
        customers = Table(id = 1:3, name = ["Alice", "Bob", "Charlie"], address = ["12 Beach Street", "163 Moon Road", "6 George Street"])
        orders = Table(customer_id = [2, 2, 3, 3], items = ["Socks", "Tie", "Shirt", "Underwear"])

        t = innerjoin(getproperty(:id), getproperty(:customer_id), customers, orders)

        @test t isa Table
        @test t == Table(id = [2, 2, 3, 3],
                         name = ["Bob", "Bob", "Charlie", "Charlie"],
                         address = ["163 Moon Road", "163 Moon Road", "6 George Street", "6 George Street"],
                         customer_id = [2, 2, 3, 3],
                         items = ["Socks", "Tie", "Shirt", "Underwear"])

        # Issue #34
        orders2 = Table(customer_id = [2, 2, 3, 3], items = ["Socks", "Tie", "Shirt", missing])
        t2 = innerjoin(getproperty(:id), getproperty(:customer_id), customers, orders2)
        @test t2 isa Table
        @test isequal(t2, Table(id = [2, 2, 3, 3],
                          name = ["Bob", "Bob", "Charlie", "Charlie"],
                          address = ["163 Moon Road", "163 Moon Road", "6 George Street", "6 George Street"],
                          customer_id = [2, 2, 3, 3],
                          items = ["Socks", "Tie", "Shirt", missing]))
    end

    @testset "adapt" begin
        tbl = Table(a = randn(Float32, 10^2), b = rand(Float64, 10^2), c = rand(1:100, 10^2))
        @test typeof(@inferred adapt(TestArrayConverter(), tbl)) == typeof(tbl)
        adapted_tbl = adapt(TestArrayConverter(), tbl)
        @test propertynames(adapted_tbl) == propertynames(tbl)
        @test adapted_tbl == tbl
    end
end


text:
using TensorOperations, IndexNotation
using Base.Test

if VERSION.minor < 3
    vecnorm(x) = norm(vec(x))
    include("tests2.jl")
else
    include("tests3.jl")
end


text:
type FastCodingData{T}
  Dtx_fft::Array{Complex{T}, 3};
  DtD_fft::Array{Complex{T}, 3};
  convH::Int;
  convW::Int;
end
type FastCodingOpts{T}
  z::Array{T, 2};
  z_fft::Array{Complex{T}, 3};
  t::Array{T, 2};
  t_fft::Array{Complex{T}, 3};
  lambda_t::Array{T, 2};
  lambda_t_fft::Array{Complex{T}, 3};
  beta::T;       # 1
  mu_t::T;       # 1
  max_iter::Int; # 1000
  tol::T;        # 1e-10
  DEBUG::Bool;   # false
  debugData::Dict;
end
FastCodingOpts(T) = FastCodingOpts{T}(zeros(0, 0), zeros(0, 0, 0), zeros(0, 0), zeros(0, 0, 0), zeros(0, 0), zeros(0, 0, 0),
                      1, 1, 1000, 1e-10, false, Dict());

type LearningData{T}
  Ztx_fft::Array{Complex{T}, 2};
  ZtZ_fft::Array{Complex{T}, 3};
  convH::Int;
  convW::Int;
  filterSize::Int;
end
type LearningOpts{T}
  Dinit::Array{T, 3};
  s::Array{T, 3};
  lambda_s::Array{T, 3};
  mu_s::T;       # 1
  max_iter::Int; # 1000
  tol::T;        # 1e-10
  DEBUG::Bool;   # false
  debugData::Dict;
end
type LearningDebug{T}
  rerr::Array{T, 1};
  obj_d::Array{T, 1};
  obj_s::Array{T, 1};
  obj::Array{T, 1};
end
LearningOpts(T) = LearningOpts{T}(zeros(0, 0, 0), zeros(0, 0, 0), zeros(0, 0, 0), 1., 1000, 1e-10, false, Dict());

type FCSCOpts{T}
  Dinit::Array{T, 3};
  beta::T;       # 1
  mu_s::T;       # 1
  mu_t::T;       # 1
  max_iter::Int; # 1000
  tol::T;        # 1e-10
  DEBUG::Bool;   # false
  FIXED_D::Bool; # false
end
type FCSCDebug{T}
  norm_z::Array{T, 1};
  rerr::Array{T, 1};
  obj::Array{T, 1};
  dif_ds::Array{T, 1};
  dif_zt::Array{T, 1};
end
FCSCOpts(T) = FCSCOpts{T}(zeros(0, 0, 0), 1., 1., 1., 1000, 1e-10, false, false);

max_stars_repo_name: bkong/FCSC


text:
mutable struct Universe
    min::Float64
    max::Float64
    width::Float64
    halfwidth::Float64
    bins::Int
    spectra::Array{Float64,1}
end

"""
    Universe(mini, maxi; width = nothing, bins = nothing)

Creates a 1-D discretized segment that starts at mini and ends at maxi. The width of the bins for the discretization
can either be provided or inferred from the number of bins. Returns a Universe object.
"""
function Universe(mini, maxi; width = nothing, bins = nothing)
    nowidth = isa( width, Nothing )
    nobins =  isa( bins, Nothing )
    @assert ( nowidth + nobins ) < 2
    if nowidth
        width = (maxi - mini ) / bins
    end
    if nobins
        bins = round((maxi - mini ) / width) |> Int
    end
    return Universe( mini, maxi, width, width / 2.0, bins, zeros( bins ) )
end

"""
    SpectralArray(Universes::Array{Universe,1})

Takes an array of Universe types and returns a 2-Array of the spectra.
"""
function SpectralArray(Universes::Array{Universe,1})
    RetArray = zeros( length(Universes), length(Universes[1].spectra) )
    for (spectra, uni) in enumerate(Universes)
        RetArray[ spectra, : ] = uni.spectra
    end
    return RetArray
end

"""
    GaussianBand(sigma,amplitude,center)

Constructs a Gaussian kernel generator.
"""
struct GaussianBand
    sigma::Float64
    amplitude::Float64
    center::Float64
end

"""
    (B::GaussianBand)(X::Float64)

Returns the scalar probability associated with a GaussianBand object (kernel) at a location in space(`X`).
"""
(B::GaussianBand)(X::Float64) = B.amplitude * exp( -0.5 * ( ( X - B.center ) / B.sigma) ^ 2 ) / sqrt( 2.0 * pi * B.sigma )

"""
    LorentzianBand(gamma,amplitude,center)

Constructs a Lorentzian kernel generator.
"""
struct LorentzianBand
    gamma::Float64
    amplitude::Float64
    center::Float64
end

"""
    (B::LorentzianBand)(X::Float64)

Returns the probability associated with a LorentzianBand object (kernel) at a location in space(`X`).
"""
(B::LorentzianBand)(X::Float64) = B.amplitude / ( pi * B.gamma * (1 + ( ( X - B.center ) / B.gamma) ^ 2 ) )

"""
    (U::Universe)(Band::Union{ GaussianBand, LorentzianBand})

A Universe objects internal "spectra" can be updated to include the additive contribution of any Band-like object.
"""
function (U::Universe)(Band::Union{ GaussianBand, LorentzianBand})
    for b in 1 : U.bins
        U.spectra[b] += Band( U.min + ( b * U.width ) + U.halfwidth )
    end
end

"""
    (U::Universe)(Band...)

A Universe objects internal "spectra" can be updated to include the additive contribution of many Band-like objects.
"""
function (U::Universe)(Band...)
    for bin in 1 : U.bins, band in Band
        U.spectra[bin] += band( U.min + ( bin * U.width ) + U.halfwidth )
    end
end


text:
# set up

include("set_up.jl")

using JLD
using HDF5

################################################################################
##      parameters for training data                                          ##
################################################################################

# set nbr parameters
set_nbr_params = 2

# nbr particles
nbr_particels = 25

# set burn-in
burn_in = 2000

# set nbr of cores
nbr_of_cores = 1

# set data to use
sim_data = true
mcmc_alg = "MCWM"  # set MCWM or PMCMC
data_set = "old"

# scale for prior dist
log_scale_prior = false

# set dt's
dt = 0.035 # new = 0.5 old = 0.03
dt_U = 1. # new = 1 old = 1

# set length for training and test data
length_training_data = 1000
length_test_data = 1000
length_training_data = length_training_data + length_test_data

# set nbr of iterations
nbr_iterations = burn_in+length_training_data

################################################################################
##                         set model parameters  for training                 ##
################################################################################

problem_training = set_up_problem(nbr_of_unknown_parameters=set_nbr_params, use_sim_data = sim_data, data_set = data_set)
problem_training.alg_param.alg = mcmc_alg
problem_training.alg_param.R = nbr_iterations
problem_training.alg_param.N = nbr_particels
problem_training.alg_param.burn_in = burn_in
problem_training.alg_param.nbr_of_cores = nbr_of_cores
problem_training.alg_param.dt = dt
problem_training.alg_param.dt_U = dt_U
problem_training.alg_param.alg = mcmc_alg
problem_training.adaptive_update =  AMUpdate_gen(eye(set_nbr_params), 1/sqrt(set_nbr_params), 0.15, 1, 0.8, 25)

################################################################################
##                generate training data                                     ###
################################################################################

if !log_scale_prior
	tic()
	res_training, theta_training, loglik_training, cov_matrix = MCMC(problem_training, true, true)
	time_pre_er = toc()
	#export_parameters(res_problem_normal_prior_est_AM_gen[2],jobname)
else
	tic()
	res_training, theta_training, loglik_training, cov_matrix  = @time MCMC(problem_training_nonlog, true, true)
	time_pre_er = toc()
end

################################################################################
##                store generated data                                       ###
################################################################################

# export generated data
export_data(problem_training, res_training[1],"gp_training_$(set_nbr_params)_par")

# split training and test data

theta_test = theta_training[1:length_training_data+1:end]
loglik_test = loglik_training[1:length_training_data+1:end]

theta_training = theta_training[1:length_training_data]
loglik_training = loglik_training[1:length_training_data]

# save training data, test data, and covaraince matrix to a Julia workspace file
save("gp_training_$(set_nbr_params)_par_training_and_test_data.jld", "res_training", res_training, "theta_training", theta_training, "loglik_training", loglik_training, "theta_test", theta_test, "loglik_test", loglik_test,"cov_matrix",cov_matrix)

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
#=
Given an unordered list of flights taken by someone, each represented as (origin, destination) pairs, and a starting airport, compute the person's itinerary. If no such itinerary exists, return null. If there are multiple possible itineraries, return the lexicographically smallest one. All flights must be used in the itinerary.

For example, given the list of flights [('SFO', 'HKO'), ('YYZ', 'SFO'), ('YUL', 'YYZ'), ('HKO', 'ORD')] and starting airport 'YUL', you should return the list ['YUL', 'YYZ', 'SFO', 'HKO', 'ORD'].

Given the list of flights [('SFO', 'COM'), ('COM', 'YYZ')] and starting airport 'COM', you should return null.

Given the list of flights [('A', 'B'), ('A', 'C'), ('B', 'C'), ('C', 'A')] and starting airport 'A', you should return the list ['A', 'B', 'C', 'A', 'C'] even though ['A', 'C', 'A', 'B', 'C'] is also a valid itinerary. However, the first one is lexicographically smaller.
=#
using Test
include("Solutions/problem41_find_flight_itinerary.jl")

@test finditinerary([("SFO", "HKO"), ("YYZ", "SFO"), ("YUL", "YYZ"), ("HKO", "ORD")], "YUL") == ["YUL", "YYZ", "SFO", "HKO", "ORD"]
@test isnothing(finditinerary([("SFO", "COM"), ("COM", "YYZ")], "COM"))
@test finditinerary([("A", "B"), ("A", "C"), ("B", "C"), ("C", "A")], "A") == ["A", "B", "C", "A", "C"]

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


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

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

"""

Get the specified network interface ip configuration in a virtual machine scale set.
Param: resourceGroupName::String (required)
Param: virtualMachineScaleSetName::String (required)
Param: virtualmachineIndex::String (required)
Param: networkInterfaceName::String (required)
Param: ipConfigurationName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: expand::String
Return: NetworkInterfaceIPConfiguration
"""
function _swaggerinternal_networkInterfacesGetVirtualMachineScaleSetIpConfiguration(_api::NetworkInterfacesApi, resourceGroupName::String, virtualMachineScaleSetName::String, virtualmachineIndex::String, networkInterfaceName::String, ipConfigurationName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkInterfaceIPConfiguration, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/microsoft.Compute/virtualMachineScaleSets/{virtualMachineScaleSetName}/virtualMachines/{virtualmachineIndex}/networkInterfaces/{networkInterfaceName}/ipConfigurations/{ipConfigurationName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualMachineScaleSetName", virtualMachineScaleSetName)  # type String
    Swagger.set_param(_ctx.path, "virtualmachineIndex", virtualmachineIndex)  # type String
    Swagger.set_param(_ctx.path, "networkInterfaceName", networkInterfaceName)  # type String
    Swagger.set_param(_ctx.path, "ipConfigurationName", ipConfigurationName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfacesGetVirtualMachineScaleSetIpConfiguration(_api::NetworkInterfacesApi, resourceGroupName::String, virtualMachineScaleSetName::String, virtualmachineIndex::String, networkInterfaceName::String, ipConfigurationName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesGetVirtualMachineScaleSetIpConfiguration(_api, resourceGroupName, virtualMachineScaleSetName, virtualmachineIndex, networkInterfaceName, ipConfigurationName, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfacesGetVirtualMachineScaleSetIpConfiguration(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, virtualMachineScaleSetName::String, virtualmachineIndex::String, networkInterfaceName::String, ipConfigurationName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesGetVirtualMachineScaleSetIpConfiguration(_api, resourceGroupName, virtualMachineScaleSetName, virtualmachineIndex, networkInterfaceName, ipConfigurationName, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Get the specified network interface in a virtual machine scale set.
Param: resourceGroupName::String (required)
Param: virtualMachineScaleSetName::String (required)
Param: virtualmachineIndex::String (required)
Param: networkInterfaceName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: expand::String
Return: NetworkInterface
"""
function _swaggerinternal_networkInterfacesGetVirtualMachineScaleSetNetworkInterface(_api::NetworkInterfacesApi, resourceGroupName::String, virtualMachineScaleSetName::String, virtualmachineIndex::String, networkInterfaceName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkInterface, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/microsoft.Compute/virtualMachineScaleSets/{virtualMachineScaleSetName}/virtualMachines/{virtualmachineIndex}/networkInterfaces/{networkInterfaceName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualMachineScaleSetName", virtualMachineScaleSetName)  # type String
    Swagger.set_param(_ctx.path, "virtualmachineIndex", virtualmachineIndex)  # type String
    Swagger.set_param(_ctx.path, "networkInterfaceName", networkInterfaceName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfacesGetVirtualMachineScaleSetNetworkInterface(_api::NetworkInterfacesApi, resourceGroupName::String, virtualMachineScaleSetName::String, virtualmachineIndex::String, networkInterfaceName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesGetVirtualMachineScaleSetNetworkInterface(_api, resourceGroupName, virtualMachineScaleSetName, virtualmachineIndex, networkInterfaceName, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfacesGetVirtualMachineScaleSetNetworkInterface(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, virtualMachineScaleSetName::String, virtualmachineIndex::String, networkInterfaceName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesGetVirtualMachineScaleSetNetworkInterface(_api, resourceGroupName, virtualMachineScaleSetName, virtualmachineIndex, networkInterfaceName, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Get the specified network interface ip configuration in a virtual machine scale set.
Param: resourceGroupName::String (required)
Param: virtualMachineScaleSetName::String (required)
Param: virtualmachineIndex::String (required)
Param: networkInterfaceName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: expand::String
Return: NetworkInterfaceIPConfigurationListResult
"""
function _swaggerinternal_networkInterfacesListVirtualMachineScaleSetIpConfigurations(_api::NetworkInterfacesApi, resourceGroupName::String, virtualMachineScaleSetName::String, virtualmachineIndex::String, networkInterfaceName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkInterfaceIPConfigurationListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/microsoft.Compute/virtualMachineScaleSets/{virtualMachineScaleSetName}/virtualMachines/{virtualmachineIndex}/networkInterfaces/{networkInterfaceName}/ipConfigurations", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualMachineScaleSetName", virtualMachineScaleSetName)  # type String
    Swagger.set_param(_ctx.path, "virtualmachineIndex", virtualmachineIndex)  # type String
    Swagger.set_param(_ctx.path, "networkInterfaceName", networkInterfaceName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfacesListVirtualMachineScaleSetIpConfigurations(_api::NetworkInterfacesApi, resourceGroupName::String, virtualMachineScaleSetName::String, virtualmachineIndex::String, networkInterfaceName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesListVirtualMachineScaleSetIpConfigurations(_api, resourceGroupName, virtualMachineScaleSetName, virtualmachineIndex, networkInterfaceName, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfacesListVirtualMachineScaleSetIpConfigurations(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, virtualMachineScaleSetName::String, virtualmachineIndex::String, networkInterfaceName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesListVirtualMachineScaleSetIpConfigurations(_api, resourceGroupName, virtualMachineScaleSetName, virtualmachineIndex, networkInterfaceName, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets all network interfaces in a virtual machine scale set.
Param: resourceGroupName::String (required)
Param: virtualMachineScaleSetName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkInterfaceListResult
"""
function _swaggerinternal_networkInterfacesListVirtualMachineScaleSetNetworkInterfaces(_api::NetworkInterfacesApi, resourceGroupName::String, virtualMachineScaleSetName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkInterfaceListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/microsoft.Compute/virtualMachineScaleSets/{virtualMachineScaleSetName}/networkInterfaces", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualMachineScaleSetName", virtualMachineScaleSetName)  # 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"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfacesListVirtualMachineScaleSetNetworkInterfaces(_api::NetworkInterfacesApi, resourceGroupName::String, virtualMachineScaleSetName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesListVirtualMachineScaleSetNetworkInterfaces(_api, resourceGroupName, virtualMachineScaleSetName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfacesListVirtualMachineScaleSetNetworkInterfaces(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, virtualMachineScaleSetName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesListVirtualMachineScaleSetNetworkInterfaces(_api, resourceGroupName, virtualMachineScaleSetName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets information about all network interfaces in a virtual machine in a virtual machine scale set.
Param: resourceGroupName::String (required)
Param: virtualMachineScaleSetName::String (required)
Param: virtualmachineIndex::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkInterfaceListResult
"""
function _swaggerinternal_networkInterfacesListVirtualMachineScaleSetVMNetworkInterfaces(_api::NetworkInterfacesApi, resourceGroupName::String, virtualMachineScaleSetName::String, virtualmachineIndex::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkInterfaceListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/microsoft.Compute/virtualMachineScaleSets/{virtualMachineScaleSetName}/virtualMachines/{virtualmachineIndex}/networkInterfaces", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualMachineScaleSetName", virtualMachineScaleSetName)  # type String
    Swagger.set_param(_ctx.path, "virtualmachineIndex", virtualmachineIndex)  # 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"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfacesListVirtualMachineScaleSetVMNetworkInterfaces(_api::NetworkInterfacesApi, resourceGroupName::String, virtualMachineScaleSetName::String, virtualmachineIndex::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesListVirtualMachineScaleSetVMNetworkInterfaces(_api, resourceGroupName, virtualMachineScaleSetName, virtualmachineIndex, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfacesListVirtualMachineScaleSetVMNetworkInterfaces(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, virtualMachineScaleSetName::String, virtualmachineIndex::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesListVirtualMachineScaleSetVMNetworkInterfaces(_api, resourceGroupName, virtualMachineScaleSetName, virtualmachineIndex, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export networkInterfacesGetVirtualMachineScaleSetIpConfiguration, networkInterfacesGetVirtualMachineScaleSetNetworkInterface, networkInterfacesListVirtualMachineScaleSetIpConfigurations, networkInterfacesListVirtualMachineScaleSetNetworkInterfaces, networkInterfacesListVirtualMachineScaleSetVMNetworkInterfaces


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

@test integrate(x^3*cot(a+b*x), x) == :(-1//4*im*x^4+b^-1*x^3*log(1+-1*exp(2*im*(a+b*x)))+3//2*x*b^-3*Polylogarithms.polylog(3,exp(2*im*(a+b*x)))+3//4*im*b^-4*Polylogarithms.polylog(4,exp(2*im*(a+b*x)))+-3//2*im*b^-2*x^2*Polylogarithms.polylog(2,exp(2*im*(a+b*x))))
@test integrate(x^2*cot(a+b*x), x) == :((1/2)*b^-3*Polylogarithms.polylog(3,exp(2*im*(a+b*x)))+-1//3*im*x^3+b^-1*x^2*log(1+-1*exp(2*im*(a+b*x)))+-1*im*x*b^-2*Polylogarithms.polylog(2,exp(2*im*(a+b*x))))
@test integrate(x*cot(a+b*x), x) == :(-1//2*im*x^2+x*b^-1*log(1+-1*exp(2*im*(a+b*x)))+-1//2*im*b^-2*Polylogarithms.polylog(2,exp(2*im*(a+b*x))))
@test integrate(x^3*cot(a+b*x)^2, x) == :(-1//4*x^4+3//2*b^-4*Polylogarithms.polylog(3,exp(2*im*(a+b*x)))+-1*im*b^-1*x^3+-1*b^-1*x^3*cot(a+b*x)+3*b^-2*x^2*log(1+-1*exp(2*im*(a+b*x)))+-3*im*x*b^-3*Polylogarithms.polylog(2,exp(2*im*(a+b*x))))
@test integrate(x^2*cot(a+b*x)^2, x) == :(-1//3*x^3+-1*im*b^-1*x^2+-1*im*b^-3*Polylogarithms.polylog(2,exp(2*im*(a+b*x)))+-1*b^-1*x^2*cot(a+b*x)+2*x*b^-2*log(1+-1*exp(2*im*(a+b*x))))
@test integrate(x*cot(a+b*x)^2, x) == :(-1//2*x^2+b^-2*log(sin(a+b*x))+-1*x*b^-1*cot(a+b*x))
@test integrate(x^3*cot(a+b*x)^3, x) == :(-1//2*b^-1*x^3+1//4*im*x^4+-1*b^-1*x^3*log(1+-1*exp(2*im*(a+b*x)))+3*x*b^-3*log(1+-1*exp(2*im*(a+b*x)))+-3//2*im*b^-4*Polylogarithms.polylog(2,exp(2*im*(a+b*x)))+-3//2*im*b^-2*x^2+-3//2*x*b^-3*Polylogarithms.polylog(3,exp(2*im*(a+b*x)))+-3//2*b^-2*x^2*cot(a+b*x)+-3//4*im*b^-4*Polylogarithms.polylog(4,exp(2*im*(a+b*x)))+-1//2*b^-1*x^3*cot(a+b*x)^2+3//2*im*b^-2*x^2*Polylogarithms.polylog(2,exp(2*im*(a+b*x))))
@test integrate(x^2*cot(a+b*x)^3, x) == :(b^-3*log(sin(a+b*x))+-1//2*b^-1*x^2+-1//2*b^-3*Polylogarithms.polylog(3,exp(2*im*(a+b*x)))+1//3*im*x^3+-1*x*b^-2*cot(a+b*x)+-1*b^-1*x^2*log(1+-1*exp(2*im*(a+b*x)))+-1//2*b^-1*x^2*cot(a+b*x)^2+im*x*b^-2*Polylogarithms.polylog(2,exp(2*im*(a+b*x))))
@test integrate(x*cot(a+b*x)^3, x) == :((1/2)*im*x^2+-1//2*x*b^-1+-1//2*b^-2*cot(a+b*x)+(1/2)*im*b^-2*Polylogarithms.polylog(2,exp(2*im*(a+b*x)))+-1*x*b^-1*log(1+-1*exp(2*im*(a+b*x)))+-1//2*x*b^-1*cot(a+b*x)^2)
@test integrate((a+im*a*cot(e+f*x))^-1*(c+d*x)^3, x) == :(-3//8*d^3*f^-4*(a+im*a*cot(e+f*x))^-1+1//8*a^-1*d^-1*(c+d*x)^4+-3//8*d*a^-1*f^-2*(c+d*x)^2+-1//2*im*f^-1*(a+im*a*cot(e+f*x))^-1*(c+d*x)^3+1//4*im*a^-1*f^-1*(c+d*x)^3+3//4*d*f^-2*(a+im*a*cot(e+f*x))^-1*(c+d*x)^2+-3//8*im*x*a^-1*d^3*f^-3+3//4*im*d^2*f^-3*(a+im*a*cot(e+f*x))^-1*(c+d*x))
@test integrate((a+im*a*cot(e+f*x))^-1*(c+d*x)^2, x) == :(1//6*a^-1*d^-1*(c+d*x)^3+(1/2)*d*f^-2*(a+im*a*cot(e+f*x))^-1*(c+d*x)+-1//2*im*f^-1*(a+im*a*cot(e+f*x))^-1*(c+d*x)^2+-1//4*x*a^-1*d^2*f^-2+1//4*im*a^-1*f^-1*(c+d*x)^2+1//4*im*d^2*f^-3*(a+im*a*cot(e+f*x))^-1)
@test integrate((a+im*a*cot(e+f*x))^-1*(c+d*x), x) == :(1//4*d*f^-2*(a+im*a*cot(e+f*x))^-1+1//4*a^-1*d^-1*(c+d*x)^2+-1//2*im*f^-1*(a+im*a*cot(e+f*x))^-1*(c+d*x)+1//4*im*d*x*a^-1*f^-1)
@test integrate((a+im*a*cot(e+f*x))^-1*(c+d*x)^-1, x) == :((1/2)*a^-1*d^-1*log(c+d*x)+-1//2*a^-1*d^-1*cos(-2*e+2*c*f*d^-1)*cosintegral(2*f*x+2*c*f*d^-1)+-1//2*a^-1*d^-1*sin(-2*e+2*c*f*d^-1)*sinintegral(2*f*x+2*c*f*d^-1)+(1/2)*im*a^-1*d^-1*cosintegral(2*f*x+2*c*f*d^-1)*sin(-2*e+2*c*f*d^-1)+-1//2*im*a^-1*d^-1*cos(-2*e+2*c*f*d^-1)*sinintegral(2*f*x+2*c*f*d^-1))
@test integrate((a+im*a*cot(e+f*x))^-1*(c+d*x)^-2, x) == :(-1*d^-1*(a+im*a*cot(e+f*x))^-1*(c+d*x)^-1+f*a^-1*d^-2*cos(-2*e+2*c*f*d^-1)*sinintegral(2*f*x+2*c*f*d^-1)+-1*f*a^-1*d^-2*cosintegral(2*f*x+2*c*f*d^-1)*sin(-2*e+2*c*f*d^-1)+-1*im*f*a^-1*d^-2*cos(-2*e+2*c*f*d^-1)*cosintegral(2*f*x+2*c*f*d^-1)+-1*im*f*a^-1*d^-2*sin(-2*e+2*c*f*d^-1)*sinintegral(2*f*x+2*c*f*d^-1))
@test integrate((a+im*a*cot(e+f*x))^-1*(c+d*x)^-3, x) == :(-1//2*d^-1*(a+im*a*cot(e+f*x))^-1*(c+d*x)^-2+a^-1*d^-3*f^2*cos(-2*e+2*c*f*d^-1)*cosintegral(2*f*x+2*c*f*d^-1)+a^-1*d^-3*f^2*sin(-2*e+2*c*f*d^-1)*sinintegral(2*f*x+2*c*f*d^-1)+(1/2)*im*f*a^-1*d^-2*(c+d*x)^-1+-1*im*f*d^-2*(a+im*a*cot(e+f*x))^-1*(c+d*x)^-1+im*a^-1*d^-3*f^2*cos(-2*e+2*c*f*d^-1)*sinintegral(2*f*x+2*c*f*d^-1)+-1*im*a^-1*d^-3*f^2*cosintegral(2*f*x+2*c*f*d^-1)*sin(-2*e+2*c*f*d^-1))
@test integrate((a+im*a*cot(e+f*x))^-2*(c+d*x)^3, x) == :(1//16*a^-2*d^-1*(c+d*x)^4+-3//512*a^-2*d^3*f^-4*exp(4*im*e+4*im*f*x)+3//16*a^-2*d^3*f^-4*exp(2*im*e+2*im*f*x)+-3//8*d*a^-2*f^-2*(c+d*x)^2*exp(2*im*e+2*im*f*x)+-1//16*im*a^-2*f^-1*(c+d*x)^3*exp(4*im*e+4*im*f*x)+1//4*im*a^-2*f^-1*(c+d*x)^3*exp(2*im*e+2*im*f*x)+3//64*d*a^-2*f^-2*(c+d*x)^2*exp(4*im*e+4*im*f*x)+-3//8*im*a^-2*d^2*f^-3*(c+d*x)*exp(2*im*e+2*im*f*x)+3//128*im*a^-2*d^2*f^-3*(c+d*x)*exp(4*im*e+4*im*f*x))
@test integrate((a+im*a*cot(e+f*x))^-2*(c+d*x)^2, x) == :(1//12*a^-2*d^-1*(c+d*x)^3+-1//4*d*a^-2*f^-2*(c+d*x)*exp(2*im*e+2*im*f*x)+-1//8*im*a^-2*d^2*f^-3*exp(2*im*e+2*im*f*x)+-1//16*im*a^-2*f^-1*(c+d*x)^2*exp(4*im*e+4*im*f*x)+1//4*im*a^-2*f^-1*(c+d*x)^2*exp(2*im*e+2*im*f*x)+1//32*d*a^-2*f^-2*(c+d*x)*exp(4*im*e+4*im*f*x)+1//128*im*a^-2*d^2*f^-3*exp(4*im*e+4*im*f*x))
@test integrate((a+im*a*cot(e+f*x))^-2*(c+d*x), x) == :(-1//8*d*a^-2*x^2+1//4*x*a^-2*(c+d*x)+1//16*d*f^-2*(a+im*a*cot(e+f*x))^-2+3//16*d*f^-2*(a^2+im*a^2*cot(e+f*x))^-1+-1//4*im*f^-1*(a+im*a*cot(e+f*x))^-2*(c+d*x)+-1//4*im*f^-1*(a^2+im*a^2*cot(e+f*x))^-1*(c+d*x)+3//16*im*d*x*a^-2*f^-1)
@test integrate((a+im*a*cot(e+f*x))^-2*(c+d*x)^-1, x) == :(1//4*a^-2*d^-1*log(c+d*x)+-1//2*a^-2*d^-1*cos(-2*e+2*c*f*d^-1)*cosintegral(2*f*x+2*c*f*d^-1)+-1//2*a^-2*d^-1*sin(-2*e+2*c*f*d^-1)*sinintegral(2*f*x+2*c*f*d^-1)+1//4*a^-2*d^-1*cos(-4*e+4*c*f*d^-1)*cosintegral(4*f*x+4*c*f*d^-1)+1//4*a^-2*d^-1*sin(-4*e+4*c*f*d^-1)*sinintegral(4*f*x+4*c*f*d^-1)+(1/2)*im*a^-2*d^-1*cosintegral(2*f*x+2*c*f*d^-1)*sin(-2*e+2*c*f*d^-1)+-1//2*im*a^-2*d^-1*cos(-2*e+2*c*f*d^-1)*sinintegral(2*f*x+2*c*f*d^-1)+-1//4*im*a^-2*d^-1*cosintegral(4*f*x+4*c*f*d^-1)*sin(-4*e+4*c*f*d^-1)+1//4*im*a^-2*d^-1*cos(-4*e+4*c*f*d^-1)*sinintegral(4*f*x+4*c*f*d^-1))
@test integrate((a+im*a*cot(e+f*x))^-2*(c+d*x)^-2, x) == :(-1//4*a^-2*d^-1*(c+d*x)^-1+(1/2)*a^-2*d^-1*(c+d*x)^-1*cos(2*e+2*f*x)+-1//4*a^-2*d^-1*(c+d*x)^-1*cos(2*e+2*f*x)^2+1//4*a^-2*d^-1*(c+d*x)^-1*sin(2*e+2*f*x)^2+f*a^-2*d^-2*cos(-2*e+2*c*f*d^-1)*sinintegral(2*f*x+2*c*f*d^-1)+f*a^-2*d^-2*cosintegral(4*f*x+4*c*f*d^-1)*sin(-4*e+4*c*f*d^-1)+(1/2)*im*a^-2*d^-1*(c+d*x)^-1*sin(2*e+2*f*x)+-1*f*a^-2*d^-2*cos(-4*e+4*c*f*d^-1)*sinintegral(4*f*x+4*c*f*d^-1)+-1*f*a^-2*d^-2*cosintegral(2*f*x+2*c*f*d^-1)*sin(-2*e+2*c*f*d^-1)+-1//4*im*a^-2*d^-1*(c+d*x)^-1*sin(4*e+4*f*x)+im*f*a^-2*d^-2*cos(-4*e+4*c*f*d^-1)*cosintegral(4*f*x+4*c*f*d^-1)+im*f*a^-2*d^-2*sin(-4*e+4*c*f*d^-1)*sinintegral(4*f*x+4*c*f*d^-1)+-1*im*f*a^-2*d^-2*cos(-2*e+2*c*f*d^-1)*cosintegral(2*f*x+2*c*f*d^-1)+-1*im*f*a^-2*d^-2*sin(-2*e+2*c*f*d^-1)*sinintegral(2*f*x+2*c*f*d^-1))
@test integrate((a+im*a*cot(e+f*x))^-3*(c+d*x)^3, x) == :(1//32*a^-3*d^-1*(c+d*x)^4+-9//1024*a^-3*d^3*f^-4*exp(4*im*e+4*im*f*x)+1//1728*a^-3*d^3*f^-4*exp(6*im*e+6*im*f*x)+9//64*a^-3*d^3*f^-4*exp(2*im*e+2*im*f*x)+-9//32*d*a^-3*f^-2*(c+d*x)^2*exp(2*im*e+2*im*f*x)+-3//32*im*a^-3*f^-1*(c+d*x)^3*exp(4*im*e+4*im*f*x)+-1//96*d*a^-3*f^-2*(c+d*x)^2*exp(6*im*e+6*im*f*x)+1//48*im*a^-3*f^-1*(c+d*x)^3*exp(6*im*e+6*im*f*x)+3//16*im*a^-3*f^-1*(c+d*x)^3*exp(2*im*e+2*im*f*x)+9//128*d*a^-3*f^-2*(c+d*x)^2*exp(4*im*e+4*im*f*x)+-9//32*im*a^-3*d^2*f^-3*(c+d*x)*exp(2*im*e+2*im*f*x)+-1//288*im*a^-3*d^2*f^-3*(c+d*x)*exp(6*im*e+6*im*f*x)+9//256*im*a^-3*d^2*f^-3*(c+d*x)*exp(4*im*e+4*im*f*x))
@test integrate((a+im*a*cot(e+f*x))^-3*(c+d*x)^2, x) == :(1//24*a^-3*d^-1*(c+d*x)^3+-3//16*d*a^-3*f^-2*(c+d*x)*exp(2*im*e+2*im*f*x)+-3//32*im*a^-3*d^2*f^-3*exp(2*im*e+2*im*f*x)+-3//32*im*a^-3*f^-1*(c+d*x)^2*exp(4*im*e+4*im*f*x)+-1//144*d*a^-3*f^-2*(c+d*x)*exp(6*im*e+6*im*f*x)+-1//864*im*a^-3*d^2*f^-3*exp(6*im*e+6*im*f*x)+1//48*im*a^-3*f^-1*(c+d*x)^2*exp(6*im*e+6*im*f*x)+3//16*im*a^-3*f^-1*(c+d*x)^2*exp(2*im*e+2*im*f*x)+3//64*d*a^-3*f^-2*(c+d*x)*exp(4*im*e+4*im*f*x)+3//256*im*a^-3*d^2*f^-3*exp(4*im*e+4*im*f*x))
@test integrate((a+im*a*cot(e+f*x))^-3*(c+d*x), x) == :(-1//16*d*a^-3*x^2+1//8*x*a^-3*(c+d*x)+1//36*d*f^-2*(a+im*a*cot(e+f*x))^-3+11//96*d*f^-2*(a^3+im*a^3*cot(e+f*x))^-1+-1//6*im*f^-1*(a+im*a*cot(e+f*x))^-3*(c+d*x)+-1//8*im*f^-1*(a^3+im*a^3*cot(e+f*x))^-1*(c+d*x)+5//96*d*a^-1*f^-2*(a+im*a*cot(e+f*x))^-2+-1//8*im*a^-1*f^-1*(a+im*a*cot(e+f*x))^-2*(c+d*x)+11//96*im*d*x*a^-3*f^-1)
@test integrate((a+im*a*cot(e+f*x))^-3*(c+d*x)^-1, x) == :(1//8*a^-3*d^-1*log(c+d*x)+-3//8*a^-3*d^-1*cos(-2*e+2*c*f*d^-1)*cosintegral(2*f*x+2*c*f*d^-1)+-3//8*a^-3*d^-1*sin(-2*e+2*c*f*d^-1)*sinintegral(2*f*x+2*c*f*d^-1)+-1//8*a^-3*d^-1*cos(-6*e+6*c*f*d^-1)*cosintegral(6*f*x+6*c*f*d^-1)+-1//8*a^-3*d^-1*sin(-6*e+6*c*f*d^-1)*sinintegral(6*f*x+6*c*f*d^-1)+3//8*a^-3*d^-1*cos(-4*e+4*c*f*d^-1)*cosintegral(4*f*x+4*c*f*d^-1)+3//8*a^-3*d^-1*sin(-4*e+4*c*f*d^-1)*sinintegral(4*f*x+4*c*f*d^-1)+-3//8*im*a^-3*d^-1*cos(-2*e+2*c*f*d^-1)*sinintegral(2*f*x+2*c*f*d^-1)+-3//8*im*a^-3*d^-1*cosintegral(4*f*x+4*c*f*d^-1)*sin(-4*e+4*c*f*d^-1)+-1//8*im*a^-3*d^-1*cos(-6*e+6*c*f*d^-1)*sinintegral(6*f*x+6*c*f*d^-1)+1//8*im*a^-3*d^-1*cosintegral(6*f*x+6*c*f*d^-1)*sin(-6*e+6*c*f*d^-1)+3//8*im*a^-3*d^-1*cos(-4*e+4*c*f*d^-1)*sinintegral(4*f*x+4*c*f*d^-1)+3//8*im*a^-3*d^-1*cosintegral(2*f*x+2*c*f*d^-1)*sin(-2*e+2*c*f*d^-1))
@test integrate((a+im*a*cot(e+f*x))^-3*(c+d*x)^-2, x) == :(-1//8*a^-3*d^-1*(c+d*x)^-1+-3//8*a^-3*d^-1*(c+d*x)^-1*cos(2*e+2*f*x)^2+1//8*a^-3*d^-1*(c+d*x)^-1*cos(2*e+2*f*x)^3+3//8*a^-3*d^-1*(c+d*x)^-1*sin(2*e+2*f*x)^2+3//32*a^-3*d^-1*(c+d*x)^-1*cos(6*e+6*f*x)+9//32*a^-3*d^-1*(c+d*x)^-1*cos(2*e+2*f*x)+-3//2*f*a^-3*d^-2*cos(-4*e+4*c*f*d^-1)*sinintegral(4*f*x+4*c*f*d^-1)+-3//4*f*a^-3*d^-2*cosintegral(2*f*x+2*c*f*d^-1)*sin(-2*e+2*c*f*d^-1)+-3//4*f*a^-3*d^-2*cosintegral(6*f*x+6*c*f*d^-1)*sin(-6*e+6*c*f*d^-1)+-3//8*im*a^-3*d^-1*(c+d*x)^-1*sin(4*e+4*f*x)+-1//8*im*a^-3*d^-1*(c+d*x)^-1*sin(2*e+2*f*x)^3+3//2*f*a^-3*d^-2*cosintegral(4*f*x+4*c*f*d^-1)*sin(-4*e+4*c*f*d^-1)+3//4*f*a^-3*d^-2*cos(-6*e+6*c*f*d^-1)*sinintegral(6*f*x+6*c*f*d^-1)+3//4*f*a^-3*d^-2*cos(-2*e+2*c*f*d^-1)*sinintegral(2*f*x+2*c*f*d^-1)+3//32*im*a^-3*d^-1*(c+d*x)^-1*sin(6*e+6*f*x)+15//32*im*a^-3*d^-1*(c+d*x)^-1*sin(2*e+2*f*x)+-3//4*im*f*a^-3*d^-2*cos(-6*e+6*c*f*d^-1)*cosintegral(6*f*x+6*c*f*d^-1)+-3//4*im*f*a^-3*d^-2*cos(-2*e+2*c*f*d^-1)*cosintegral(2*f*x+2*c*f*d^-1)+-3//4*im*f*a^-3*d^-2*sin(-6*e+6*c*f*d^-1)*sinintegral(6*f*x+6*c*f*d^-1)+-3//4*im*f*a^-3*d^-2*sin(-2*e+2*c*f*d^-1)*sinintegral(2*f*x+2*c*f*d^-1)+3//2*im*f*a^-3*d^-2*cos(-4*e+4*c*f*d^-1)*cosintegral(4*f*x+4*c*f*d^-1)+3//2*im*f*a^-3*d^-2*sin(-4*e+4*c*f*d^-1)*sinintegral(4*f*x+4*c*f*d^-1))
@test integrate((a+im*a*cot(e+f*x))^-1*(c+d*x)^m, x) == :((1/2)*a^-1*d^-1*(1+m)^-1*(c+d*x)^(1+m)+im*2^(-2+-1m)*a^-1*f^-1*(-1*im*f*d^-1*(c+d*x))^(-1m)*(c+d*x)^m*SpecialFunctions.gamma(1+m,-2*im*f*d^-1*(c+d*x))*exp(2*im*(e+-1*c*f*d^-1)))
@test integrate((a+im*a*cot(e+f*x))^-2*(c+d*x)^m, x) == :(1//4*a^-2*d^-1*(1+m)^-1*(c+d*x)^(1+m)+im*2^(-2+-1m)*a^-2*f^-1*(-1*im*f*d^-1*(c+d*x))^(-1m)*(c+d*x)^m*SpecialFunctions.gamma(1+m,-2*im*f*d^-1*(c+d*x))*exp(2*im*(e+-1*c*f*d^-1))+-1*im*4^(-2+-1m)*a^-2*f^-1*(-1*im*f*d^-1*(c+d*x))^(-1m)*(c+d*x)^m*SpecialFunctions.gamma(1+m,-4*im*f*d^-1*(c+d*x))*exp(4*im*(e+-1*c*f*d^-1)))
@test integrate((a+im*a*cot(e+f*x))^-3*(c+d*x)^m, x) == :(1//8*a^-3*d^-1*(1+m)^-1*(c+d*x)^(1+m)+-3*im*2^(-5+-2m)*a^-3*f^-1*(-1*im*f*d^-1*(c+d*x))^(-1m)*(c+d*x)^m*SpecialFunctions.gamma(1+m,-4*im*f*d^-1*(c+d*x))*exp(4*im*(e+-1*c*f*d^-1))+3*im*2^(-4+-1m)*a^-3*f^-1*(-1*im*f*d^-1*(c+d*x))^(-1m)*(c+d*x)^m*SpecialFunctions.gamma(1+m,-2*im*f*d^-1*(c+d*x))*exp(2*im*(e+-1*c*f*d^-1))+im*2^(-4+-1m)*3^(-1+-1m)*a^-3*f^-1*(-1*im*f*d^-1*(c+d*x))^(-1m)*(c+d*x)^m*SpecialFunctions.gamma(1+m,-6*im*f*d^-1*(c+d*x))*exp(6*im*(e+-1*c*f*d^-1)))
@test integrate((c+d*x)^3*(a+b*cot(e+f*x)), x) == :(1//4*a*d^-1*(c+d*x)^4+b*f^-1*(c+d*x)^3*log(1+-1*exp(2*im*(e+f*x)))+-1//4*im*b*d^-1*(c+d*x)^4+3//2*b*d^2*f^-3*(c+d*x)*Polylogarithms.polylog(3,exp(2*im*(e+f*x)))+3//4*im*b*d^3*f^-4*Polylogarithms.polylog(4,exp(2*im*(e+f*x)))+-3//2*im*b*d*f^-2*(c+d*x)^2*Polylogarithms.polylog(2,exp(2*im*(e+f*x))))
@test integrate((c+d*x)^2*(a+b*cot(e+f*x)), x) == :(1//3*a*d^-1*(c+d*x)^3+b*f^-1*(c+d*x)^2*log(1+-1*exp(2*im*(e+f*x)))+(1/2)*b*d^2*f^-3*Polylogarithms.polylog(3,exp(2*im*(e+f*x)))+-1//3*im*b*d^-1*(c+d*x)^3+-1*im*b*d*f^-2*(c+d*x)*Polylogarithms.polylog(2,exp(2*im*(e+f*x))))
@test integrate((a+b*cot(e+f*x))*(c+d*x), x) == :((1/2)*a*d^-1*(c+d*x)^2+b*f^-1*(c+d*x)*log(1+-1*exp(2*im*(e+f*x)))+-1//2*im*b*d^-1*(c+d*x)^2+-1//2*im*b*d*f^-2*Polylogarithms.polylog(2,exp(2*im*(e+f*x))))
@test integrate((a+b*cot(e+f*x))^2*(c+d*x)^3, x) == :(-1//4*b^2*d^-1*(c+d*x)^4+1//4*a^2*d^-1*(c+d*x)^4+-1*im*b^2*f^-1*(c+d*x)^3+-1*b^2*f^-1*(c+d*x)^3*cot(e+f*x)+3//2*b^2*d^3*f^-4*Polylogarithms.polylog(3,exp(2*im*(e+f*x)))+2*a*b*f^-1*(c+d*x)^3*log(1+-1*exp(2*im*(e+f*x)))+3*d*b^2*f^-2*(c+d*x)^2*log(1+-1*exp(2*im*(e+f*x)))+-1//2*im*a*b*d^-1*(c+d*x)^4+-3*im*b^2*d^2*f^-3*(c+d*x)*Polylogarithms.polylog(2,exp(2*im*(e+f*x)))+3*a*b*d^2*f^-3*(c+d*x)*Polylogarithms.polylog(3,exp(2*im*(e+f*x)))+3//2*im*a*b*d^3*f^-4*Polylogarithms.polylog(4,exp(2*im*(e+f*x)))+-3*im*a*b*d*f^-2*(c+d*x)^2*Polylogarithms.polylog(2,exp(2*im*(e+f*x))))
@test integrate((a+b*cot(e+f*x))^2*(c+d*x)^2, x) == :(-1//3*b^2*d^-1*(c+d*x)^3+1//3*a^2*d^-1*(c+d*x)^3+-1*im*b^2*f^-1*(c+d*x)^2+-1*b^2*f^-1*(c+d*x)^2*cot(e+f*x)+a*b*d^2*f^-3*Polylogarithms.polylog(3,exp(2*im*(e+f*x)))+-1*im*b^2*d^2*f^-3*Polylogarithms.polylog(2,exp(2*im*(e+f*x)))+2*a*b*f^-1*(c+d*x)^2*log(1+-1*exp(2*im*(e+f*x)))+2*d*b^2*f^-2*(c+d*x)*log(1+-1*exp(2*im*(e+f*x)))+-2//3*im*a*b*d^-1*(c+d*x)^3+-2*im*a*b*d*f^-2*(c+d*x)*Polylogarithms.polylog(2,exp(2*im*(e+f*x))))
@test integrate((a+b*cot(e+f*x))^2*(c+d*x), x) == :((1/2)*a^2*d^-1*(c+d*x)^2+-1*c*x*b^2+-1//2*d*b^2*x^2+d*b^2*f^-2*log(sin(e+f*x))+-1*b^2*f^-1*(c+d*x)*cot(e+f*x)+-1*im*a*b*d^-1*(c+d*x)^2+2*a*b*f^-1*(c+d*x)*log(1+-1*exp(2*im*(e+f*x)))+-1*im*a*b*d*f^-2*Polylogarithms.polylog(2,exp(2*im*(e+f*x))))
@test integrate((a+b*cot(e+f*x))^3*(c+d*x)^3, x) == :(-1//2*b^3*f^-1*(c+d*x)^3+1//4*a^3*d^-1*(c+d*x)^4+-1*b^3*f^-1*(c+d*x)^3*log(1+-1*exp(2*im*(e+f*x)))+-3//4*a*b^2*d^-1*(c+d*x)^4+-1//2*b^3*f^-1*(c+d*x)^3*cot(e+f*x)^2+1//4*im*b^3*d^-1*(c+d*x)^4+-3*im*a*b^2*f^-1*(c+d*x)^3+-3*a*b^2*f^-1*(c+d*x)^3*cot(e+f*x)+3*b*a^2*f^-1*(c+d*x)^3*log(1+-1*exp(2*im*(e+f*x)))+3*b^3*d^2*f^-3*(c+d*x)*log(1+-1*exp(2*im*(e+f*x)))+-3//2*im*d*b^3*f^-2*(c+d*x)^2+-3//2*im*b^3*d^3*f^-4*Polylogarithms.polylog(2,exp(2*im*(e+f*x)))+-3//2*d*b^3*f^-2*(c+d*x)^2*cot(e+f*x)+-3//2*b^3*d^2*f^-3*(c+d*x)*Polylogarithms.polylog(3,exp(2*im*(e+f*x)))+-3//4*im*b*a^2*d^-1*(c+d*x)^4+-3//4*im*b^3*d^3*f^-4*Polylogarithms.polylog(4,exp(2*im*(e+f*x)))+9//2*a*b^2*d^3*f^-4*Polylogarithms.polylog(3,exp(2*im*(e+f*x)))+9*a*d*b^2*f^-2*(c+d*x)^2*log(1+-1*exp(2*im*(e+f*x)))+3//2*im*d*b^3*f^-2*(c+d*x)^2*Polylogarithms.polylog(2,exp(2*im*(e+f*x)))+9//2*b*a^2*d^2*f^-3*(c+d*x)*Polylogarithms.polylog(3,exp(2*im*(e+f*x)))+9//4*im*b*a^2*d^3*f^-4*Polylogarithms.polylog(4,exp(2*im*(e+f*x)))+-9*im*a*b^2*d^2*f^-3*(c+d*x)*Polylogarithms.polylog(2,exp(2*im*(e+f*x)))+-9//2*im*b*d*a^2*f^-2*(c+d*x)^2*Polylogarithms.polylog(2,exp(2*im*(e+f*x))))
@test integrate((a+b*cot(e+f*x))^3*(c+d*x)^2, x) == :(1//3*a^3*d^-1*(c+d*x)^3+b^3*d^2*f^-3*log(sin(e+f*x))+-1*a*b^2*d^-1*(c+d*x)^3+-1*b^3*f^-1*(c+d*x)^2*log(1+-1*exp(2*im*(e+f*x)))+-1//2*b^3*d^2*f^-1*x^2+-1//2*b^3*d^2*f^-3*Polylogarithms.polylog(3,exp(2*im*(e+f*x)))+-1//2*b^3*f^-1*(c+d*x)^2*cot(e+f*x)^2+1//3*im*b^3*d^-1*(c+d*x)^3+-1*im*b*a^2*d^-1*(c+d*x)^3+-1*c*d*x*b^3*f^-1+-1*d*b^3*f^-2*(c+d*x)*cot(e+f*x)+-3*im*a*b^2*f^-1*(c+d*x)^2+-3*a*b^2*f^-1*(c+d*x)^2*cot(e+f*x)+3*b*a^2*f^-1*(c+d*x)^2*log(1+-1*exp(2*im*(e+f*x)))+3//2*b*a^2*d^2*f^-3*Polylogarithms.polylog(3,exp(2*im*(e+f*x)))+im*d*b^3*f^-2*(c+d*x)*Polylogarithms.polylog(2,exp(2*im*(e+f*x)))+-3*im*a*b^2*d^2*f^-3*Polylogarithms.polylog(2,exp(2*im*(e+f*x)))+6*a*d*b^2*f^-2*(c+d*x)*log(1+-1*exp(2*im*(e+f*x)))+-3*im*b*d*a^2*f^-2*(c+d*x)*Polylogarithms.polylog(2,exp(2*im*(e+f*x))))
@test integrate((a+b*cot(e+f*x))^3*(c+d*x), x) == :((1/2)*a^3*d^-1*(c+d*x)^2+(1/2)*im*b^3*d^-1*(c+d*x)^2+-1*b^3*f^-1*(c+d*x)*log(1+-1*exp(2*im*(e+f*x)))+-3*a*c*x*b^2+-3//2*a*d*b^2*x^2+-1//2*d*x*b^3*f^-1+-1//2*d*b^3*f^-2*cot(e+f*x)+-1//2*b^3*f^-1*cot(e+f*x)^2*(c+d*x)+(1/2)*im*d*b^3*f^-2*Polylogarithms.polylog(2,exp(2*im*(e+f*x)))+-3*a*b^2*f^-1*(c+d*x)*cot(e+f*x)+3*a*d*b^2*f^-2*log(sin(e+f*x))+3*b*a^2*f^-1*(c+d*x)*log(1+-1*exp(2*im*(e+f*x)))+-3//2*im*b*a^2*d^-1*(c+d*x)^2+-3//2*im*b*d*a^2*f^-2*Polylogarithms.polylog(2,exp(2*im*(e+f*x))))
@test integrate((a+b*cot(e+f*x))^-1*(c+d*x)^3, x) == :(d^-1*(c+d*x)^4*(4a+-4*im*b)^-1+-1*b*f^-1*(c+d*x)^3*(a^2+b^2)^-1*log(1+-1*(a+-1*im*b)^-1*(a+im*b)*exp(2*im*(e+f*x)))+-3*im*b*d^3*f^-4*(4*a^2+4*b^2)^-1*Polylogarithms.polylog(4,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*(e+f*x)))+-3*b*d^2*f^-3*(2*a^2+2*b^2)^-1*(c+d*x)*Polylogarithms.polylog(3,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*(e+f*x)))+3*im*b*d*f^-2*(c+d*x)^2*(2*a^2+2*b^2)^-1*Polylogarithms.polylog(2,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*(e+f*x))))
@test integrate((a+b*cot(e+f*x))^-1*(c+d*x)^2, x) == :(d^-1*(c+d*x)^3*(3a+-3*im*b)^-1+-1*b*d^2*f^-3*(2*a^2+2*b^2)^-1*Polylogarithms.polylog(3,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*(e+f*x)))+-1*b*f^-1*(c+d*x)^2*(a^2+b^2)^-1*log(1+-1*(a+-1*im*b)^-1*(a+im*b)*exp(2*im*(e+f*x)))+im*b*d*f^-2*(a^2+b^2)^-1*(c+d*x)*Polylogarithms.polylog(2,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*(e+f*x))))
@test integrate((a+b*cot(e+f*x))^-1*(c+d*x), x) == :(d^-1*(c+d*x)^2*(2a+-2*im*b)^-1+-1*b*f^-1*(a^2+b^2)^-1*(c+d*x)*log(1+-1*(a+-1*im*b)^-1*(a+im*b)*exp(2*im*(e+f*x)))+im*b*d*f^-2*(2*a^2+2*b^2)^-1*Polylogarithms.polylog(2,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*(e+f*x))))
@test integrate((a+b*cot(e+f*x))^-2*(c+d*x)^3, x) == :(1//4*d^-1*(a+im*b)^-2*(c+d*x)^4+-1*b^2*d^-1*(c+d*x)^4*(a^2+b^2)^-2+-1*b*d^-1*(a+im*b)^-2*(b+im*a)^-1*(c+d*x)^4+-2*im*b^2*f^-1*(c+d*x)^3*(a^2+b^2)^-2+3//2*b^2*d^3*f^-4*(a^2+b^2)^-2*Polylogarithms.polylog(3,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x))+3//2*b^2*d^3*f^-4*(a^2+b^2)^-2*Polylogarithms.polylog(4,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x))+-3*d*b^2*f^-2*(c+d*x)^2*(a^2+b^2)^-2*Polylogarithms.polylog(2,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x))+-2*im*b^2*f^-1*(c+d*x)^3*(a^2+b^2)^-2*log(1+-1*(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x))+-2*b*f^-1*(a+im*b)^-2*(a+-1*im*b)^-1*(c+d*x)^3*log(1+-1*(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x))+-2*b^2*f^-1*(a+im*b)^-2*(a+-1*im*b)^-1*(c+d*x)^3*(b+im*a+-1*(-1b+im*a)*exp(2*im*e+2*im*f*x))^-1+3*d*b^2*f^-2*(c+d*x)^2*(a^2+b^2)^-2*log(1+-1*(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x))+3//2*b*d^3*f^-4*(a+im*b)^-2*(b+im*a)^-1*Polylogarithms.polylog(4,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x))+-3*im*b^2*d^2*f^-3*(a^2+b^2)^-2*(c+d*x)*Polylogarithms.polylog(2,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x))+-3*im*b^2*d^2*f^-3*(a^2+b^2)^-2*(c+d*x)*Polylogarithms.polylog(3,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x))+-3*b*d*f^-2*(a+im*b)^-2*(b+im*a)^-1*(c+d*x)^2*Polylogarithms.polylog(2,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x))+-3*b*d^2*f^-3*(a+im*b)^-2*(a+-1*im*b)^-1*(c+d*x)*Polylogarithms.polylog(3,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x)))
@test integrate((a+b*cot(e+f*x))^-2*(c+d*x)^2, x) == :(1//3*d^-1*(a+im*b)^-2*(c+d*x)^3+-4//3*b^2*d^-1*(c+d*x)^3*(a^2+b^2)^-2+-2*im*b^2*f^-1*(c+d*x)^2*(a^2+b^2)^-2+-4//3*b*d^-1*(a+im*b)^-2*(b+im*a)^-1*(c+d*x)^3+-1*im*b^2*d^2*f^-3*(a^2+b^2)^-2*Polylogarithms.polylog(2,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x))+-1*im*b^2*d^2*f^-3*(a^2+b^2)^-2*Polylogarithms.polylog(3,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x))+-1*b*d^2*f^-3*(a+im*b)^-2*(a+-1*im*b)^-1*Polylogarithms.polylog(3,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x))+-2*im*b^2*f^-1*(c+d*x)^2*(a^2+b^2)^-2*log(1+-1*(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x))+-2*b*f^-1*(a+im*b)^-2*(a+-1*im*b)^-1*(c+d*x)^2*log(1+-1*(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x))+-2*d*b^2*f^-2*(a^2+b^2)^-2*(c+d*x)*Polylogarithms.polylog(2,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x))+-2*b^2*f^-1*(a+im*b)^-2*(a+-1*im*b)^-1*(c+d*x)^2*(b+im*a+-1*(-1b+im*a)*exp(2*im*e+2*im*f*x))^-1+2*d*b^2*f^-2*(a^2+b^2)^-2*(c+d*x)*log(1+-1*(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x))+-2*b*d*f^-2*(a+im*b)^-2*(b+im*a)^-1*(c+d*x)*Polylogarithms.polylog(2,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*e+2*im*f*x)))
@test integrate((a+b*cot(e+f*x))^-2*(c+d*x), x) == :(-1*d^-1*(c+d*x)^2*(2*a^2+2*b^2)^-1+b*f^-1*(a+b*cot(e+f*x))^-1*(a^2+b^2)^-1*(c+d*x)+b*f^-2*(a^2+b^2)^-2*(b*d+-2*a*c*f+-2*a*d*f*x)*log(1+-1*(a+-1*im*b)^-1*(a+im*b)*exp(2*im*(e+f*x)))+1//4*a^-1*d^-1*f^-2*(a+im*b)^-1*(a+-1*im*b)^-2*(b*d+-2*a*c*f+-2*a*d*f*x)^2+im*a*b*d*f^-2*(a^2+b^2)^-2*Polylogarithms.polylog(2,(a+-1*im*b)^-1*(a+im*b)*exp(2*im*(e+f*x))))

max_stars_repo_path: test/4 Trig functions/4.4 Cotangent/4.4.10 (c+d x)^m (a+b cot)^n.jl


text:
import Base.write, Base.read

struct Layer{T<:Number}
    W::Matrix{T}
    b::Vector{T}
    f::Function
end

Layer(W, b) = Layer(W, b, identity)

Layer(in::Integer, out::Integer, f::Function) = Layer(randn(out, in), randn(out), f)

function (l::Layer)(x::AbstractArray)
    return l.W*x .+ l.b
end

struct Chain
    layers::Vector{<:Layer}
    function Chain(xs)
        for i in 1:length(xs)-1
            size(xs[i].W, 1) != size(xs[i+1].W, 2) && throw(ArgumentError("The former layers's out must equal to the latter's in."))
        end
        new(xs)
    end
end

Chain(xs...) = Chain([xs...])

Chain(ns::Integer...; f::Function=σ) = Chain([Layer(ns[i], ns[i+1], f) for i in 1:length(ns)-1])

Chain(c::Chain, xs...) = Chain(c.layers..., xs...)

Chain(c1::Chain, c2::Chain) = Chain(c1.layers..., c2.layers...)

function (c::Chain)(x::AbstractArray)
    for l in c.layers
        x = l.f.(l(x))
    end
    return x
end

function write(io::IO, c::Chain)
    write(io, UInt8(length(c.layers)))
    write(io, UInt16(size(c.layers[1].W, 2)))
    for l in c.layers
        write(io, UInt16(length(l.b)))
    end
    for l in c.layers
        write(io, vec(l.W))
        write(io, l.b)
    end
end

getlayersize(io::IO, n::Integer) = [read(io, UInt16) for _ in 1:n+1]

function getlayer(io::IO, lin::Integer, lout::Integer, f::Function=σ)
    w = zeros(lout, lin)
    b = zeros(lout)
    for i in 1:lin, j in 1:lout
        w[j, i] = read(io, Float64)
    end
    for i in 1:lout
       b[i] = read(io, Float64)
    end
    Layer(w, b, f)
end

function read(io::IO, ::Type{Chain}; f::Function=σ)
    n = read(io, UInt8)
    lsize = getlayersize(io, n)
    Chain([getlayer(io, lsize[i], lsize[i+1], f) for i in 1:n])
end

σ(x) = 1.0/(1.0+exp(-x))
dσ(x) = σ(x) * (1 - σ(x))


text:
using Test, ModelConstructors, Distributions, Dates, Nullables, Random

@testset "Ensure transformations to the real line/model space are valid" begin
    for T in subtypes(Transform)
        global u = parameter(:σ_pist, 2.5230, (1e-8, 5.), (1e-8, 5.), T(), fixed=false)
        @test ( transform_to_real_line(u) |> x -> transform_to_model_space(u,x) ) == u.value

        if !isa(T,Type{ModelConstructors.Untransformed})
            # check transform_to_real_line and transform_to_model_space to different things if T is not ModelConstructors.Untransformed
            @test transform_to_real_line(u,u.value) != transform_to_model_space(u,u.value)
        end
    end
end

tomodel_answers = zeros(3)
toreal_answers  = zeros(3)
a = 1e-8; b = 5.; c = 1.
x = 2.5230
cx = 2 * (x - (a+b)/2)/(b-a)
toreal_answers[1] = 1. / (c * (x - a))
toreal_answers[2] = (1/c) * (1. / (1. - cx^2)^(-3/2)) * (2/(b-a))
toreal_answers[3] = 1.
x = transform_to_real_line(parameter(:σ_pist, 2.5230, (1e-8, 5.), (1e-8, 5.),
                                     ModelConstructors.Exponential(), fixed=false))
tomodel_answers[1] = c * exp(c * (x - b))
x = transform_to_real_line(parameter(:σ_pist, 2.5230, (1e-8, 5.), (1e-8, 5.),
                                     ModelConstructors.SquareRoot(), fixed=false))
tomodel_answers[2] = (b - a) / 2 * c / (1 + c^2 * x^2)^(3/2)
tomodel_answers[3] = 1.
@testset "Ensure derivatives of transformations to the real line/model space are valid" begin
    for (i,T) in enumerate(subtypes(Transform))
        global u = parameter_ad(:σ_pist, 2.5230, (1e-8, 5.), (1e-8, 5.), T(), fixed=false)
        @test differentiate_transform_to_real_line(u, u.value) == toreal_answers[i]
        global x = transform_to_real_line(u)
        @test differentiate_transform_to_model_space(u,x) == tomodel_answers[i]

        if !isa(T,Type{ModelConstructors.Untransformed})
            # check transform_to_real_line and transform_to_model_space to different things if T is not ModelConstructors.Untransformed
            @test differentiate_transform_to_real_line(u,u.value) != differentiate_transform_to_model_space(u,u.value)
        end
    end
end

@testset "Check type conversion for `value`, `valuebounds`, and `transform_parameterization`" begin
    @info "The following warning is expected"
    u1 = parameter(:σ_pist, 2.5230, (1, 5), (Float32(1), Float32(5)), Untransformed())
    u2 = parameter(:σ_pist, 2.5230, (1., 5.), (1., 5.), Untransformed())

    @test u1.value == u2.value
    @test u1.valuebounds == u2.valuebounds
    @test u1.transform_parameterization == u2.transform_parameterization
    @test typeof(u1.value) == typeof(u2.value)
    @test eltype(u1.valuebounds) == eltype(u2.valuebounds)
    @test eltype(u1.transform_parameterization) == eltype(u2.transform_parameterization)
end

# probability
N = 10^2
u = parameter(:bloop, 2.5230, (1e-8, 5.), (1e-8, 5.), ModelConstructors.SquareRoot(); fixed = true)
v = parameter(:cat, 2.5230, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), Gamma(2.00, 0.1))

pvec =  ParameterVector{Float64}(undef, N)
for i in 1:length(pvec)
	global pvec[i] = (i%2 == 0) ? u : v
end
@testset "Check logpdf/pdf function approximations" begin
    @test logpdf(pvec) ≈ 50*logpdf(v)
    @test pdf(pvec) ≈ exp(50*logpdf(v))
end

updated = ModelConstructors.update(pvec, ones(length(pvec)))
ModelConstructors.update!(pvec, ones(length(pvec)))

@testset "Check if update! preserves dimensions and values" begin
    @test all(updated .== pvec)
    @test logpdf(pvec) == logpdf(updated)
end

# test we only update non-fixed parameters
@testset "Ensure only non-fixed parameters are updated" begin
    for p in pvec
        if p.fixed
            @test p.value == 2.5230
        elseif isa(p, Parameter)
            @test p.value == one(Float64)
        end
    end
end

# vector of new values must be the same length
@testset "Ensure update! enforces the same length of the parameter vector being updated" begin
    @test_throws BoundsError ModelConstructors.update!(pvec, ones(length(pvec)-1))
end

@testset "Ensure parameters being updated are of the same type." begin
    for w in [parameter_ad(:moop, 3.0, fixed=false), parameter_ad(:moop, 3.0; scaling = log, fixed=false)]
        # new values must be of the same type
        @test_throws ErrorException parameter_ad(w, one(Int))

        # new value is out of bounds
        @test_throws ParamBoundsError parameter_ad(w, -1.)
    end
end

@testset "Ensure parameter value types can be changed when forced by keyword." begin
    for w in [parameter_ad(:moop, 1.0, fixed=false), parameter_ad(:moop, 1.0; scaling = log, fixed=false)]
        # new values must be of the same type
        @test typeof(parameter_ad(w, one(Int); change_value_type = true).value) == Int
    end
end

# Tests moved from DSGE.jl core.jl
# UnscaledParameter, fixed=false
α =  parameter(:α, 0.1596, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Normal(0.30, 0.05), fixed=false)
@testset "Test non-fixed UnscaledParameter" begin
    @test isa(α, UnscaledParameter)
    @test α.key == :α
    @test isa(α.prior.value, Normal)
    @test α.prior.value.μ == 0.3
    @test α.description == "No description available."
    @test α.tex_label == ""
    @test isa(α.transform, ModelConstructors.SquareRoot)
    @test get_untransformed_values(α) == α.value
end

# UnscaledParameter, fixed = true
α_fixed =  parameter(:α_fixed, 0.1596, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.Untransformed(), Normal(0.30, 0.05), fixed=true)
@testset "Test fixed UnscaledParameter" begin
    @test α_fixed.transform_parameterization == (0.1596,0.1596)
    @test isa(α_fixed.transform, ModelConstructors.Untransformed)
    @test get_untransformed_values(α_fixed) == α_fixed.value
end

# UnscaledParameter, fixed = true, transform should be overwritten given fixed
α_fixed =  parameter(:α_fixed, 0.1596, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Normal(0.30, 0.05), fixed=true)
@testset "Test fixed UnscaledParameter, ensuring transform is overwritten" begin
    @test isa(α_fixed.transform, ModelConstructors.Untransformed)
end

# Fixed UnscaledParameter, minimal constructor
δ = parameter(:δ, 0.025)
@testset "Test fixed UnscaledParameter minimal constructor" begin
    @test δ.fixed
    @test δ.transform_parameterization == (0.025, 0.025)
    @test δ.valuebounds == (0.025, 0.025)
end

# Scaled parameter
β = parameter(:β, 0.1402, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), GammaAlt(0.25, 0.1), fixed=false,  scaling = x -> (1 + x/100)\1, description="β: Discount rate.", tex_label="\\beta ")
@testset "Test ScaledParameter constructor" begin
    @test isa(β, ScaledParameter)
    @test isa(β.prior.value, Gamma)
    @test isa(β.transform, ModelConstructors.Exponential)
    @test get_untransformed_values(β) == β.scaledvalue
end

# Invalid transform
@testset "Ensure error thrown on invalid transform" begin
    @test_throws UndefVarError α_bad = parameter(:α, 0.1596, (1e-5, 0.999), (1e-5, 0.999),
                                                 InvalidTransform(), Normal(0.30, 0.05), fixed=false)
end

# Arithmetic with parameters
@testset "Check arithmetic with parameters" begin
    @test promote_type(AbstractParameter{Float64}, Float16) == Float64
    @test promote_type(AbstractParameter{Float64}, Int8) == Float64
    ## @test promote_rule(AbstractParameter{Float64}, Float16) == Float64
    ## @test promote_rule(AbstractParameter{Float64}, Int8) == Float64
    @test δ + δ == 0.05
    @test δ^2 == 0.025^2
    @test -δ == -0.025
    @test log(δ) == log(0.025)
end

# transform_to_real_line and transform_to_model_space
cx = 2 * (α - 1/2)
@testset "Check parameter transformations for optimization" begin
    @test abs(transform_to_real_line(α) - cx / sqrt(1 - cx^2)) < .001
    @test transform_to_real_line(δ) == 0.025
end

# Check parameters can be declared as vectors
@testset "Test vector-valued parameters" begin
    p = parameter(:test, ones(5), (0.0, 5.), (0., 5.0), Untransformed(), MvNormal(ones(5), ones(5)), fixed = false, scaling = x -> x/2)
    @test all(ones(5) ./ 2 .== parameter(p, ones(5)).scaledvalue)
end

# Check we can update only specific parameter values in a ParameterVector
@testset "Test updating specific values of a ParameterVector" begin
    # Check you can overwrite values which are unfixed
    global pvec = ParameterVector{Float64}(undef, 3)
    pvec[1] = parameter(:a, 1., (0., 3.), (0., 3.), fixed = false)
    pvec[2] = parameter(:b, 1.)
    pvec[3] = parameter(:c, 1., (0., 3.), (0., 3.), fixed = false)
    global vals = [2., 2.]
    update!(pvec, vals, BitArray([true, false, true]))
    @test map(x -> x.value, pvec) == [2., 1., 2.]

    # Check that overwriting fixed parameter values does not work, even if you pass
    # the BitArray telling update! to change the values.
    pvec[1] = parameter(:a, 1.)
    pvec[2] = parameter(:b, 1.)
    pvec[3] = parameter(:c, 1.)
    global vals = [2., 2.]
    update!(pvec, vals, BitArray([true, false, true]))
    @test map(x -> x.value, pvec) == [1., 1., 1.]
end

@testset "Broadcasting with Parameter types" begin
    a = parameter(:a, 1., (0., 1.), (0., 1.), Untransformed(), fixed = false)
    b = parameter(:b, 1., (0., 1.), (0., 1.), Untransformed(), fixed = false, scaling = x -> x / 100.)

    rand_vals = rand(3)
    p=SteadyStateParameterGrid(:sspg, rand(3))
    p2=SteadyStateParameterGrid(:sspg, rand(3,2))
    rand1 = rand()

    @test a .* ones(3) == a.value * ones(3)
    @test a .* rand_vals == a.value * rand_vals
    @test b .* ones(3) == b.scaledvalue * ones(3)
    @test b .* rand_vals == b.scaledvalue * rand_vals
    @test p .* 1 == p.value
    @test p .* rand1 == p.value .* rand1
    @test p .* rand_vals == p.value .* rand_vals
    @test a .* ones(3) == a.value * ones(3)
    @test a .+ rand_vals == a.value .+ rand_vals
    @test b .+ ones(3) == b.scaledvalue .+ ones(3)
    @test b .+ rand_vals == b.scaledvalue .+ rand_vals
    @test p .+ 1 == p.value .+ 1.
    @test p .+ rand1 == p.value .+ rand1
    @test p .+ rand_vals == p.value + rand_vals
    @test p[1] == p.value[1]
    @test first(p[1:1]) == p.value[1]
    @test p[2:3] == p.value[2:3]
    @test -p == -p.value
end

nothing

max_stars_repo_name: FRBNY-DSGE/ModelConstructors

max_issues_repo_name: FRBNY-DSGE/ModelConstructors

max_forks_repo_name: FRBNY-DSGE/ModelConstructors


text:
using ReachabilityAnalysis, SparseArrays

LazySets.set_ztol(Float64, 1e-14)

function gearbox(; X0 = Hyperrectangle([0., 0., -0.0167, 0.003,  0., 0.],
                                       [0., 0.,  0.0001, 0.0001, 0., 0.]),
                   init = [(1, X0)])

    # variables
    vx = 1  # x velocity
    vy = 2  # y velocity
    px = 3  # x position
    py = 4  # y position
    I = 5   # accumulated impulse
    t = 6   # time

    # number of variables
    n = 5 + 1

    # constants
    Fs = 70.0  # shifting force on the sleeve [N]
    Tf = 1.0  # resisting moments on the second gear [N⋅m]
    # TODO where do the following values come from?
    ms = 3.2  # mass of the sleeve
    Rs = 0.08  # radius of the sleeve
    Jg₂ = 0.7  # inertia of the second gear
    Δp = -0.003  # horizontal (px) distance
    θ = 0.628318530717959  # included angle of the second gear [°]
    mg₂ = 18.1  # mass of the second gear
    ζ = 0.9  # coefficient of restitution

    # discrete structure (graph)
    automaton = LightAutomaton(2)
    add_transition!(automaton, 1, 1, 1)
    add_transition!(automaton, 1, 1, 2)
    add_transition!(automaton, 1, 2, 3)

    # mode 1 ("free")
    A = zeros(n, n)
    b = zeros(n)
    A[px, vx] = 1.
    A[py, vy] = 1.
    b[vx] = Fs / ms
    b[vy] = - (Rs * Tf) / Jg₂
    b[t] = 1.
    invariant = HPolyhedron([
        HalfSpace(sparsevec([px], [1.], n), Δp),  # px <= Δp
        HalfSpace(sparsevec([px, py], [tan(θ), 1.], n), 0.),    # py <= -px * tan(θ)
        HalfSpace(sparsevec([px, py], [tan(θ), -1.], n), 0.)])  # py >= px * tan(θ)
    m_1 = @system(x' = Ax + b, x ∈ invariant)

    # mode 2 ("meshed")
    A = zeros(n, n)
    b = zeros(n)
    b[t] = 1.
    m_2 = @system(x' = Ax + b, x ∈ Universe(n))

    # common assignment matrix (requires individual modifications)
    A_template = zeros(n, n)
    for i in 3:6
        A_template[i, i] = 1.
    end
    denominator = ms * cos(θ)^2 + mg₂ * sin(θ)^2
    A_template[vx, vx] = (ms * cos(θ)^2 - mg₂ * ζ * sin(θ)^2) / denominator
    A_template[vx, vy] = (-(ζ + 1.) * mg₂ * sin(θ) * cos(θ)) / denominator
    A_template[vy, vx] = (-(ζ + 1.) * ms * sin(θ) * cos(θ)) / denominator
    A_template[vy, vy] = (mg₂ * sin(θ)^2 - ms * ζ * cos(θ)^2) / denominator
    A_template[I, vx] = ((ζ + 1.) * ms * mg₂ * sin(θ)) / denominator
    A_template[I, vy] = ((ζ + 1.) * ms * mg₂ * cos(θ)) / denominator

    # transition l1 -> l1
    # TODO what happened to the term '2nb' and the whole second constraint in the paper?
    guard = HPolyhedron([
        HalfSpace(sparsevec([px, py], [-tan(θ), -1.], n), 0.),     # py >= -px * tan(θ)
        HalfSpace(sparsevec([vx, vy], [-sin(θ), -cos(θ)], n), 0.)  # vx * sin(θ) + vy * cos(θ) >= 0
        ])
    A = copy(A_template)
    t1 = ConstrainedLinearMap(A, guard)

    # transition l1 -> l1
    # TODO same remark as with the other guard
    guard = HPolyhedron([
        HalfSpace(sparsevec([px, py], [-tan(θ), 1.], n), 0.),     # py <= px * tan(θ)
        HalfSpace(sparsevec([vx, vy], [-sin(θ), cos(θ)], n), 0.)  # vx * sin(θ) - vy * cos(θ) >= 0
        ])
    A = copy(A_template)
    A[vx, vy] *= -1.
    A[vy, vx] *= -1.
    A[I, vy] *= -1.
    t2 = ConstrainedLinearMap(A, guard)

    # transition l1 -> l2
    guard = HalfSpace(sparsevec([px], [-1.], n), -Δp)  # px >= Δp
    A = copy(A_template)
    A[vx, vx] = 0.
    A[vx, vy] = 0.
    A[vy, vx] = 0.
    A[vy, vy] = 0.
    A[I, vx] = A[I, vy] = ms
    t3 = ConstrainedLinearMap(A, guard)

    H = HybridSystem(automaton=automaton,
                     modes=[m_1, m_2],
                     resetmaps=[t1, t2, t3])

    return InitialValueProblem(H, init)
end

# Requirements
function req1(sol)
    fp = sol.F
    for idx in 1:length(fp)
        if location(fp[idx]) == 1
            verif = tend(fp[idx]) < 0.2
            !verif && return false
        end
    end
    return true
end

function req2(sol)
    return ρ(sparsevec([5], [1.0], 6), sol) < 20.0
end

# ======= Homgeneized version ==========

using ReachabilityAnalysis: add_dimension

function gearbox_homog(; X0 = Hyperrectangle([0., 0., -0.0167, 0.003,  0., 1.],
                                             [0., 0.,  0.0001, 0.0001, 0., 0.]),
                         init = [(1, X0)])

    # variables
    vx = 1  # x velocity
    vy = 2  # y velocity
    px = 3  # x position
    py = 4  # y position
    I = 5   # accumulated impulse
    η = 6   # accounts for the homgeneized input term

    # number of variables
    n = 5 + 1

    # constants
    Fs = 70.0  # shifting force on the sleeve [N]
    Tf = 1.0  # resisting moments on the second gear [N⋅m]
    # TODO where do the following values come from?
    ms = 3.2  # mass of the sleeve
    Rs = 0.08  # radius of the sleeve
    Jg₂ = 0.7  # inertia of the second gear
    Δp = -0.003  # horizontal (px) distance
    θ = 0.628318530717959  # included angle of the second gear [°]
    mg₂ = 18.1  # mass of the second gear
    ζ = 0.9  # coefficient of restitution

    # discrete structure (graph)
    automaton = LightAutomaton(2)
    add_transition!(automaton, 1, 1, 1)
    add_transition!(automaton, 1, 1, 2)
    add_transition!(automaton, 1, 2, 3)

    # mode 1 ("free")
    A = zeros(n-1, n-1)
    b = zeros(n-1)
    A[px, vx] = 1.
    A[py, vy] = 1.
    b[vx] = Fs / ms
    b[vy] = - (Rs * Tf) / Jg₂
    invariant = HPolyhedron([
        HalfSpace(sparsevec([px], [1.], n), Δp),  # px <= Δp
        HalfSpace(sparsevec([px, py], [tan(θ), 1.], n), 0.),    # py <= -px * tan(θ)
        HalfSpace(sparsevec([px, py], [tan(θ), -1.], n), 0.)])  # py >= px * tan(θ)
    Aext = add_dimension(A)
    Aext[1:n-1, n] = b
    m_1 = @system(x' = Aext * x, x ∈ invariant)

    # mode 2 ("meshed")
    A0 = zeros(n, n)
    m_2 = @system(x' = A0 * x, x ∈ Universe(n))

    # common assignment matrix (requires individual modifications)
    A_template = zeros(n, n)
    for i in 3:6
        A_template[i, i] = 1.
    end
    denominator = ms * cos(θ)^2 + mg₂ * sin(θ)^2
    A_template[vx, vx] = (ms * cos(θ)^2 - mg₂ * ζ * sin(θ)^2) / denominator
    A_template[vx, vy] = (-(ζ + 1.) * mg₂ * sin(θ) * cos(θ)) / denominator
    A_template[vy, vx] = (-(ζ + 1.) * ms * sin(θ) * cos(θ)) / denominator
    A_template[vy, vy] = (mg₂ * sin(θ)^2 - ms * ζ * cos(θ)^2) / denominator
    A_template[I, vx] = ((ζ + 1.) * ms * mg₂ * sin(θ)) / denominator
    A_template[I, vy] = ((ζ + 1.) * ms * mg₂ * cos(θ)) / denominator

    # transition l1 -> l1
    # TODO what happened to the term '2nb' and the whole second constraint in the paper?
    guard = HPolyhedron([
        HalfSpace(sparsevec([px, py], [-tan(θ), -1.], n), 0.),     # py >= -px * tan(θ)
        HalfSpace(sparsevec([vx, vy], [-sin(θ), -cos(θ)], n), 0.)  # vx * sin(θ) + vy * cos(θ) >= 0
        ])
    A = copy(A_template)
    t1 = ConstrainedLinearMap(A, guard)

    # transition l1 -> l1
    # TODO same remark as with the other guard
    guard = HPolyhedron([
        HalfSpace(sparsevec([px, py], [-tan(θ), 1.], n), 0.),     # py <= px * tan(θ)
        HalfSpace(sparsevec([vx, vy], [-sin(θ), cos(θ)], n), 0.)  # vx * sin(θ) - vy * cos(θ) >= 0
        ])
    A = copy(A_template)
    A[vx, vy] *= -1.
    A[vy, vx] *= -1.
    A[I, vy] *= -1.
    t2 = ConstrainedLinearMap(A, guard)

    # transition l1 -> l2
    guard = HalfSpace(sparsevec([px], [-1.], n), -Δp)  # px >= Δp
    A = copy(A_template)
    A[vx, vx] = 0.
    A[vx, vy] = 0.
    A[vy, vx] = 0.
    A[vy, vy] = 0.
    A[I, vx] = A[I, vy] = ms
    t3 = ConstrainedLinearMap(A, guard)

    H = HybridSystem(automaton=automaton,
                     modes=[m_1, m_2],
                     resetmaps=[t1, t2, t3])

    return InitialValueProblem(H, init)
end

function _project_dirs_34(extdirs)
    vec = Vector{Vector{Float64}}()
    for d in extdirs.directions
        d34 = d[3:4]
        if !iszero(d34)
            push!(vec, d34)
        end
    end
    out = CustomDirections(vec)
end


text:
using PyPlot, SparseArrays
rows = [2,3,4,4,5,5,5,6,6,7,7,8,8,9,9,10,11]
cols = [3,2,1,2,2,4,6,2,5,2,5,2,5,2,5,5,5]
A = sparse(rows, cols, 1, 11, 11)
A[sum(A,dims=2)[:] .== 0,:] .= 1  # If a vertex has no outbound edges, all edges to all other vertices
spy(A, marker=".");

function pagerank(A, d=0.85, tol=1e-6)
    N = size(A,2)
    L = sum(A, dims=2)
    M = (A ./ L)'
    R = ones(N) / N
    while true
        newR = d * (M * R) .+ (1 - d) / N
        if maximum(abs.(newR - R)) < tol
            return newR
        end
        R = newR
    end
end

R = pagerank(A, 0.85)
round.(100R, digits=1) # Show percentages


text:
const UNREGISTERED_PACKAGES = MirrorUpdater.Types.SrcDestPair[
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url =
            "https://github.com/DilumAluthge/DelayedErrors.jl",
        destination_repo_name =
            "DelayedErrors.jl-b0cb2517-35d6-4513-a0ca-b0b388a189e4",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url =
            "https://github.com/JuliaMath/Libm.jl",
        destination_repo_name =
            "JuliaMath-Libm.jl",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url =
            "https://github.com/alan-turing-institute/MLJRegistry.jl",
        destination_repo_name =
            "MLJRegistry.jl-a41df0fa-2d8a-11e9-1bfb-8110be68cd3e",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url =
            "https://github.com/UnofficialJuliaMirror/MirrorUpdater.jl",
        destination_repo_name =
            "MirrorUpdater.jl-7f0fc4bf-4ffe-4e21-ab91-7bf1358d5ab3",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url =
            "https://github.com/UnofficialJuliaMirrorSnapshots/RemoveLFS.jl",
        destination_repo_name =
            "RemoveLFS.jl-a262596c-dd27-4d54-b019-3619024addbd",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url =
            "https://github.com/UnofficialJuliaMirrorSnapshots/Snapshots.jl",
        destination_repo_name =
            "Snapshots.jl-44eb87bc-f37b-45e8-9f53-3bcb453a652d",
        ),
    ]