text:
## packages
using Pkg; Pkg.activate(pwd())
using Revise
using LinearAlgebra, Random, Statistics, Distributions
using FileIO, JLD2, CSV
using Plots; gr()
dir = pwd()

## parallel packages
using Distributed
nworkers() == 1 && addprocs()
@everywhere using Pkg
@everywhere Pkg.activate(pwd())
@everywhere using Revise
@everywhere dir = pwd()
@everywhere include(joinpath(dir, "src\\ops_utility.jl"))
@everywhere include(joinpath(dir, "src\\ops_build.jl"))
@everywhere include(joinpath(dir, "src\\ops_methods.jl"))
@everywhere include(joinpath(dir, "src\\cfrops_solve.jl"))
@everywhere include(joinpath(dir, "src\\dbr_cfrops_solve.jl"))
@everywhere include(joinpath(dir, "src\\ccfrops_solve.jl"))
@everywhere include(joinpath(dir, "src\\optops_solve.jl"))
@everywhere include(joinpath(dir, "src\\ops_results.jl"))

## build g, gs, gos
function build_robust_game()
    println("Building robust game...")
    fn = joinpath(dir, "data\\vars_opt_rewardfunction_expected_15city_flipped.jld2")
    fload = load(fn) # loads g and reward_exp
    g, (reward_exp, _) = fload["g"], fload["reward_exp"]
    gs = GameSet(g)
    gos = GameOptSet(g, gs, reward_exp)
    return g, gs, gos
end

# build opponent strategy and uncertainty set
@everywhere function build_opponent_strat(g, gs, gos, σfix, conf::Float64, sd::Float64)
    noise = Normal(0, sd)
    hw = quantile(noise, 1 - (1 - conf) / 2)
    σlb = [max.(σi .- hw, 0.0) for σi in σfix]
    σub = [min.(σi .+ hw, 1.0) for σi in σfix]
    mean(mean.(σub[i] - σlb[i] for i in 1:length(σlb)))
    rp2fix = real_strat(σfix, gs, gos)[2]
    rp2lb = real_strat(σlb, gs, gos)[2]
    rp2ub = real_strat(σub, gs, gos)[2]
    rp2lb_ccfr = real_strat_struct_combined(σlb, g)
    rp2ub_ccfr = real_strat_struct_combined(σub, g)
    !all(rp2lb .<= rp2fix .<= rp2ub) && @warn "Negative realizaiton probabilities."
    !all(rp2lb_ccfr .<= rp2ub_ccfr) && @warn "Negative realizaiton probabilities."
    return rp2fix, rp2lb, rp2ub, rp2lb_ccfr, rp2ub_ccfr, noise
end

@everywhere function calc_robust_strat(g, gs, gos, rp2fix, rp2lb, rp2ub, rp2lb_ccfr, rp2ub_ccfr, σfix, pconf)
    T_min = 30
    T_sec = T_min * 60
    println("Calculating NE LP strategy...")
    status_ne, u_ne, r1_ne, r2_ne, t_ne = lp_nash(gos, timelimit = T_sec)
    println("Calculating BR LP strategy...")
    status_br, u_br, r1_br = lp_best_response(gos, rp2fix, fixedplayer = 2)
    println("Calculating Robust BR LP strategy...")
    status_r, u_r, r1_r, r2_r, t_r = lp_robust_br(gos, rp2lb, rp2ub, timelimit = T_sec)
    println("Calculating Robust DBR CFR strategy...")
    (u_dbr, _, _, σ_dbr, status_dbr), t_dbr = @timed cfr_dbr(g, gs, pconf, σfix,
            tol = 5e-9, timelimit = T_min, iterlimit = 25_000)
    r1_dbr = real_strat(σ_dbr, gs, gos)[1]
    println("Calculating Robust Constrained CFR strategy...")
    (u_ccfr, _, _, σ_ccfr, _, status_ccfr, _), t_ccfr = @timed ccfr_full(g, gs, rp2lb_ccfr, rp2ub_ccfr,
            timelimit = T_min, tol = 5e-9, iterlimit = 25_000,
            α = 7.4, β = -2.9, γ = 6.3, discounted = true,
            λmax = 2000, λscale = 5)
    r1_ccfr = real_strat(σ_ccfr, gs, gos)[1]
    return (r1_ne, r1_br, r1_r, r1_dbr, r1_ccfr), (t_ne, 0.0, t_r, t_dbr, t_ccfr),
        (status_ne, status_br, status_r, status_dbr, status_ccfr)
end

@everywhere function build_robust(g, gs, gos, conf::Float64, sd::Float64, opplevel::Int)
    pconf = fill(conf, gs.ni)
    _, _, _, σfix = cfr(opplevel, g, gs)
    σfix .= [zeros(length(x)) for x in σfix]
    for i = 1:length(σfix)
        σfix[i][1] = 1.0
    end
    println("Building opponent strategies...")
    rp2fix, rp2lb, rp2ub, rp2lb_ccfr, rp2ub_ccfr, noise = build_opponent_strat(g, gs, gos, σfix, conf, sd)
    r1s, ts, statuses = calc_robust_strat(g, gs, gos, rp2fix, rp2lb, rp2ub, rp2lb_ccfr, rp2ub_ccfr, σfix, pconf)
    return r1s, ts, statuses, σfix, noise
end

@everywhere function simrobust(σfix_rand, gs, gos, r1s, ind_nz; reps = 10)
    r1_ne, r1_br, r1_r, r1_dbr, r1_ccfr = r1s
    r2_rand = real_strat(σfix_rand, gs, gos)[2]
    u1_ne = sum(gos.reward[I] * r1_ne[I[1]] * r2_rand[I[2]] for I in ind_nz)
    u1_br = sum(gos.reward[I] * r1_br[I[1]] * r2_rand[I[2]] for I in ind_nz)
    u1_r = sum(gos.reward[I] * r1_r[I[1]] * r2_rand[I[2]] for I in ind_nz)
    u1_dbr = sum(gos.reward[I] * r1_dbr[I[1]] * r2_rand[I[2]] for I in ind_nz)
    u1_ccfr = sum(gos.reward[I] * r1_ccfr[I[1]] * r2_rand[I[2]] for I in ind_nz)
    return u1_ne, u1_br, u1_r, u1_dbr, u1_ccfr
end

@everywhere function calc_robust(gs, gos, r1s, σfix, noise, nreps)
    println("Building noisy opponent strategies...")
    σfix_rands = [[normalize(clamp.(x .+ rand(noise, length(x)), 0.0, 1.0), 1) for x in σfix] for _ in 1:nreps]
    !all(all(all(y .>= 0) for y in x) for x in σfix_rands) && @warn "Negative probabilities."
    println("Calculating utility of robust strategies...")
    ind_nz = findall(!iszero, gos.reward)
    u1s = pmap(x -> simrobust(x, gs, gos, r1s, ind_nz, reps = nreps), σfix_rands)
    u1_nes = [x[1] for x in u1s]
    u1_brs = [x[2] for x in u1s]
    u1_rs = [x[3] for x in u1s]
    u1_dbrs = [x[4] for x in u1s]
    u1_ccfrs = [x[5] for x in u1s]
    return (u1_nes, u1_brs, u1_rs, u1_dbrs, u1_ccfrs)
end

function plot_robust(u1tuple)
    u1_nes, u1_brs, u1_rs, u1_dbrs, u1_ccfrs = u1tuple
    fig = histogram(u1_nes, bins=:auto, color=:mediumseagreen, alpha=0.2,
                xlabel = "Defender Utility", ylabel = "Count",
                label = "Nash Equilibrium Strategy",
                xlims = (-0.420, -0.390))
    # histogram!(fig, u1_brs, bins=:auto, color=:tomato, alpha=0.2,
    #             label = "Best Response Strategy")
    histogram!(fig, u1_rs, bins=:auto, color=:tomato, alpha=0.95,
                label = "Robust LP Strategy")
    histogram!(fig, u1_dbrs, bins=:auto, color=:steelblue, alpha=0.4,
                label = "Data Biased CFR Strategy")
    histogram!(fig, u1_ccfrs, bins=:auto, color=:gold, alpha=0.6,
                label = "Constrained CFR Strategy")
    fig
end

@everywhere signerr(u::Number, u_true::Number) = (u - u_true) / abs(u_true)

@everywhere function process_robust(u1tuple, r1s, gos, nreps)
    means = mean.(u1tuple)
    stds = std.(u1tuple)
    tstar = 1.98
    confintervals = (means .- tstar .* stds ./ sqrt(nreps), means .+ tstar .* stds ./ sqrt(nreps))
    u1_nes, u1_brs, u1_rs, u1_dbrs, u1_ccfrs = u1tuple
    df_robust = DataFrame(NE = u1_nes, RLP = u1_rs, DBR = u1_dbrs, CCFR = u1_ccfrs)
    r1_ne, r1_br, r1_r, r1_dbr, r1_ccfr = r1s
    br_results = pmap(x -> lp_best_response(gos, x, fixedplayer = 1), [r1_ne, r1_r, r1_dbr, r1_ccfr])
    br_ne, br_r, br_dbr, br_ccfr = collect(x[2] for x in br_results)
    # exploitability
    exp_neg = zeros(3)
    exp_neg[1] = signerr(br_r, br_ne)
    exp_neg[2] = signerr(br_dbr, br_ne)
    exp_neg[3] = signerr(br_ccfr, br_ne)
    # exploitation
    exp_pos = zeros(3)
    exp_pos[1] = signerr(means[3], means[1]) # r to ne
    exp_pos[2] = signerr(means[4], means[1]) # dbr to ne
    exp_pos[3] = signerr(means[5], means[1]) # ccfr to ne
    return means[3], means[4], means[5], exp_neg[1], exp_neg[2], exp_neg[3],
        exp_pos[1], exp_pos[2], exp_pos[3]
end

@everywhere function size_experiment_vector(gvec, gsvec, gosvec)
    nruns = length(gvec)
    df_size = DataFrame(size = Int[],
        status_r = Symbol[], status_dbr = Bool[], status_ccfr = Bool[],
        t_r = Float64[], t_dbr = Float64[], t_ccfr = Float64[],
        u_r = Float64[], u_dbr = Float64[], u_ccfr = Float64[],
        exp_neg_r = Float64[], exp_neg_dbr = Float64[], exp_neg_ccfr = Float64[],
        exp_pos_r = Float64[], exp_pos_dbr = Float64[], exp_pos_ccfr = Float64[])
    conf = 0.99
    sd = 0.02
    opplevel = 1
    nreps = 180
    for i in 1:length(gvec)
        print("Run $i:                                   ")
        g = gvec[i]
        gs = gsvec[i]
        gos = gosvec[i]
        print("\rRun $i: Building strategies...")
        (r1s, ts, statuses, σfix, noise), buildtime = @timed build_robust(g, gs, gos, conf, sd, opplevel)
        print("\rRun $i: Simulating value...")
        u1tuple, calctime = @timed calc_robust(gs, gos, r1s, σfix, noise, nreps)
        summary = vcat(i + 6, statuses[3:5]..., ts[3:5]..., process_robust(u1tuple, r1s, gos, nreps)...)
        push!(df_size, summary)
    end
    return df_size
end

@everywhere function size_experiment(g, gs, gos)
    conf = 0.99
    sd = 0.02
    opplevel = 1
    nreps = 180
    print("\rBuilding strategies...")
    (r1s, ts, statuses, σfix, noise), buildtime = @timed build_robust(g, gs, gos, conf, sd, opplevel)
    print("\rSimulating value...")
    u1tuple, calctime = @timed calc_robust(gs, gos, r1s, σfix, noise, nreps)
    return vcat(i + 6, statuses[3:5]..., ts[3:5]..., process_robust(u1tuple, r1s, gos, nreps)...)
end

## run
dirlocal = "C:\\Users\\AKeith\\JuliaProjectsLocal\\StrategyGamesLocal"
fnlocal = joinpath(dirlocal, "data\\vars_opt_rewardfunction_expected_7to15city_temp.jld2")
@load joinpath(dirlocal, "data\\vars_opt_rewardfunction_expected_7to15city.jld2") results # approx 5 min
ncities = [7, 8, 9, 10, 11, 12, 13, 14, 15] # number of cities (all other game params fixed)
n = length(ncities)
gvec = [ri[1] for ri in results]
gsvec = GameSet.(gvec)
rexps = [ri[4] for ri in results]
gosvec = [GameOptSet(gvec[i], gsvec[i], rexps[i]) for i in 1:n]
gsizes = [gamesize(gi)[4] for gi in gvec]
seqsizes = [size(ri[7][1])[1] * size(ri[7][2])[1] for ri in results]
statuses = [ri[5] for ri in results]
runtimes = [ri[end] for ri in results] ./ (60 * 60) # in hours

# result_size = pmap(x -> size_experiment(gvec[x], gsvec[x], gosvec[x]), 1:2)
df_size, t_exp = @timed size_experiment_vector(gvec, gsvec, gosvec)
println("Robust Experiment Total Time: $t_exp")
fn = "data\\robust_results_temp.csv"
CSV.write(joinpath(dir, fn), df_size)

# dfs = df_size[1:8, :]
dfs = df_size[1:8, :]
function plot_rel_robust(dir)
    plot(dfs[:size], signerr.(dfs[:u_r], dfs[:u_r]) .* 100,
        seriestype = :line, linecolor = :black, linestyle = :dash,
        legend = :bottomright, label = "Robust LP", ylim = (-6, 1),
        xlabel = "Problem Size (Number of Cities)", ylabel = "Relative Utility (Percent)")
    plot!(dfs[:size], signerr.(dfs[:u_dbr], dfs[:u_r]) .* 100,
        label = "Data-Biased CFR", linecolor = :lightblue)
    plot!(dfs[:size], signerr.(dfs[:u_dbr], dfs[:u_r]) .* 100,
        lab = "",
        seriestype = :scatter, markercolor = :lightblue, markerstrokealpha = 0.0)
    plot!(dfs[:size], signerr.(dfs[:u_ccfr], dfs[:u_r]) .* 100,
        label = "Constrained CFR", linecolor = :steelblue)
    plot!(dfs[:size], signerr.(dfs[:u_ccfr], dfs[:u_r]) .* 100,
        lab = "",
        seriestype = :scatter, markercolor = :steelblue, markerstrokealpha = 0.0)
    fn = "data\\plot_rel_robust_temp.pdf"
    savefig(joinpath(dir, fn))
end
plot_rel_robust(dir)

plot(dfr[:u_r] - dfr[:u_ccfr], label = "CCFR")
plot!(dfr[:u_r] - dfr[:u_dbr], label = "DBR")

plot(dfr[:size], dfr[:u_r], label = "Robust LP",
    ylim = (-1.0, 0.0), palette = :viridis)
plot!(dfr[:size], dfr[:u_dbr], label = "Data-Biased CFR")
plot!(dfr[:size], dfr[:u_ccfr], label = "Constrained CFR")

plot(dfr[:exp_neg_r] - dfr[:exp_neg_ccfr], label = "CCFR")
plot!(dfr[:exp_neg_r] - dfr[:exp_neg_dbr], label = "DBR")

plot(dfr[:exp_pos_r] - dfr[:exp_pos_ccfr], label = "CCFR")
plot!(dfr[:exp_pos_r] - dfr[:exp_pos_dbr], label = "DBR")

plot(dfr[:exp_pos_r] ./ dfr[:exp_neg_r], label = "R")
plot!(dfr[:exp_pos_dbr] ./ dfr[:exp_neg_dbr], label = "DBR")
plot!(dfr[:exp_pos_ccfr] ./ dfr[:exp_neg_ccfr], label = "CCFR")

##########################################
# old processing
conf = 0.99
sd = 0.02
opplevel = 1
nreps = 180
# (g, gs, gos), gametime = @timed build_robust_game()
(r1s, ts, statuses, σfix, noise), buildtime = @timed build_robust(g, gs, gos, conf, sd, opplevel)
@show (buildtime / 60)
# fn = "data\\vars_robust_builds_temp.pdf"
# @save joinpath(dir, fn) r1s, σfix, noise, buildtime

u1tuple, calctime = @timed calc_robust(gs, gos, r1s, σfix, noise, nreps)
@show (calctime / 60)
# fn = "data\\results_robust_calc_temp.pdf"
# @save joinpath(dir, fn) u1tuple, calctime

means = mean.(u1tuple)
stds = std.(u1tuple)
tstar = 1.98
confintervals = (means .- tstar .* stds ./ sqrt(nreps), means .+ tstar .* stds ./ sqrt(nreps))
minimum.(u1tuple)
u1_nes, u1_brs, u1_rs, u1_dbrs, u1_ccfrs = u1tuple
df_robust = DataFrame(NE = u1_nes, RLP = u1_rs, DBR = u1_dbrs, CCFR = u1_ccfrs)
fn = "data\\results_robust_data.csv"
using CSV
CSV.write(joinpath(dir, fn), df_robust)
plot_robust(u1tuple)
fn = "data\\plot_robust_hist_v2_temp.pdf"
# savefig(joinpath(dir, fn))

# _, br_ne, _ = lp_best_response(gos, r1_ne, fixedplayer = 1)
# _, br_r, _ = lp_best_response(gos, r1_r, fixedplayer = 1)
# _, br_dbr, _ = lp_best_response(gos, r1_dbr, fixedplayer = 1)
# _, br_ccfr, _ = lp_best_response(gos, r1_ccfr, fixedplayer = 1)
r1_ne, r1_br, r1_r, r1_dbr, r1_ccfr = r1s
br_results = pmap(x -> lp_best_response(gos, x, fixedplayer = 1), [r1_ne, r1_r, r1_dbr, r1_ccfr])
br_ne, br_r, br_dbr, br_ccfr = collect(x[2] for x in br_results)

# exploitability
rerr(br_r, br_ne)
rerr(br_dbr, br_ne)
rerr(br_ccfr, br_ne)

# exploitation
rerr(means[3], means[1]) # r to ne
rerr(means[4], means[1]) # dbr to ne
rerr(means[5], means[1]) # ccfr to ne

##########################################################
# Robust CFR DOE

using Pkg; Pkg.activate(pwd())
using Distributed, SharedArrays, Statistics, Dates
using CSV
# using Plots; gr()
nprocs() == 1 && addprocs()
@everywhere using Pkg
@everywhere Pkg.activate(pwd())
@everywhere(using Revise, BenchmarkTools, Printf, ProgressMeter, FileIO, JLD2)
@everywhere(using Random, LinearAlgebra, Statistics, Distances, IterTools, DataFrames, DataFramesMeta)
@everywhere(using StaticArrays, SparseArrays)
@everywhere(using JuMP, Clp)
@everywhere(using Clp: ClpCInterface)
@everywhere dir = pwd()
@everywhere include(joinpath(dir, "src\\ops_utility.jl"))
@everywhere include(joinpath(dir, "src\\ops_build.jl"))
@everywhere include(joinpath(dir, "src\\ops_methods.jl"))
@everywhere include(joinpath(dir, "src\\cfrops_solve.jl"))
@everywhere include(joinpath(dir, "src\\dcfrops_solve.jl"))
@everywhere include(joinpath(dir, "src\\ccfrops_solve.jl"))
@everywhere include(joinpath(dir, "src\\optops_solve.jl"))
@everywhere include(joinpath(dir, "src\\ops_results.jl"))

# see scratch for screening design and RSM design
@everywhere function build_design(fn::String)
    factor_names = ["alpha", "beta", "gamma", "lambda_max", "lambda_scale"]
    col_types = [Int64, Int64, Int64, Int64, Int64]
    design = CSV.read(joinpath(dir, fn), types = col_types)
    rename!(design, old => new for (old, new) = zip(names(design), Symbol.(factor_names)))
    return design
end

@everywhere function build_robust_doe(simdata)
    g, gs, gos = build_robust_game()
    conf, sd, nreps = simdata
    opplevel = 1
    _, _, _, σfix = cfr(opplevel, g, gs)
    σfix .= [zeros(length(x)) for x in σfix]
    for i = 1:length(σfix)
        σfix[i][1] = 1.0
    end
    _, _, _, rp2lb, rp2ub, noise = build_opponent_strat(g, gs, gos, σfix, conf, sd)
    gamedata = (g, gs, gos, σfix, rp2lb, rp2ub, noise)
    return gamedata
end

@everywhere function strategy_doe(gamedata, simdata, params)
    g, gs, gos, σfix, rp2lb, rp2ub, noise = gamedata
    conf, sd, nreps = simdata
    a, b, c, lmax, lscale = params
    status_ne, u_ne, r1_ne, r2_ne, t_ne = lp_nash(gos, timelimit = 5 * 60)
    (u_ccfr, _, _, σ_ccfr, _, status_ccfr, iter_ccfr), t_cfr = @timed ccfr_full(g, gs, rp2lb, rp2ub,
            timelimit = 20, tol = 5e-9, α = a, β = b, γ = c, discounted = true,
            λmax = lmax, λscale = lscale)
    r1_ccfr = real_strat(σ_ccfr, gs, gos)[1]
    return (u_ne, u_ccfr), (r1_ne, r1_ccfr), (status_ccfr, iter_ccfr, t_cfr)
end

@everywhere function simrobust_doe(r1s, gamedata, simdata, ind_nz)
    r1_ne, r1_ccfr = r1s
    _, gs, gos, σfix, _, _, noise = gamedata
    _, _, nreps = simdata
    u1_ne = zeros(nreps)
    u1_ccfr = zeros(nreps)
    for i = 1:nreps
        σfix_rand = [normalize(clamp.(x .+ rand(noise, length(x)), 0.0, 1.0), 1) for x in σfix]
        r2_rand = real_strat(σfix_rand, gs, gos)[2]
        u1_ne[i] = sum(gos.reward[I] * r1_ne[I[1]] * r2_rand[I[2]] for I in ind_nz)
        u1_ccfr[i] = sum(gos.reward[I] * r1_ccfr[I[1]] * r2_rand[I[2]] for I in ind_nz)
    end
    return u1_ne, u1_ccfr
end

@everywhere function process_doe(u1s, r1s, gamedata, simdata)
    g, gs, gos, σfix, rp2lb, rp2ub, noise = gamedata
    conf, sd, nreps = simdata
    u1_ne, u1_ccfr = u1s
    r1_ne, r1_ccfr = r1s
    means = mean.(u1s)
    stds = std.(u1s)
    tstar = 1.98
    confintervals = (means .- tstar .* stds ./ sqrt(nreps), means .+ tstar .* stds ./ sqrt(nreps))
    df_robust = DataFrame(NE = u1_ne, CCFR = u1_ccfr)
    _, br_ne, _ = lp_best_response(gos, r1_ne, fixedplayer = 1)
    _, br_ccfr, _ = lp_best_response(gos, r1_ccfr, fixedplayer = 1)
    exp_neg = signerr(br_ccfr, br_ne) # exploitability
    exp_pos = signerr(means[2], means[1]) # exploitation
    return means[2], confintervals[2][1], confintervals[2][2], exp_neg, exp_pos
end

@everywhere function calc_run(df_row, gamedata, simdata, runnum)
    g, gs, gos, σfix, rp2lb, rp2ub = gamedata
    conf, sd, nreps = simdata
    a = df_row.alpha[1]
    b = df_row.beta[1]
    c = df_row.gamma[1]
    lmax = df_row.lambda_max[1]
    lscale = df_row.lambda_scale[1]
    params = (a, b, c, lmax, lscale)
    println("----------- Start Run $runnum -----------")
    uraw, r1s, (status_ccfr, iter_ccfr, t_ccfr) = strategy_doe(gamedata, simdata, params)
    ind_nz = findall(!iszero, gos.reward)
    println("----------- Sim Robust $runnum -----------")
    u1s = simrobust_doe(r1s, gamedata, simdata, ind_nz)
    umean, ulb, uub, exp_neg, exp_pos = process_doe(u1s, r1s, gamedata, simdata)
    println("------##### End Run $runnum #####------")
    return umean, ulb, uub, exp_neg, exp_pos, status_ccfr, iter_ccfr, t_ccfr
end

design = build_design("data\\design_robust.csv")
nruns = size(design, 1)
simdata = 0.99, 0.02, 90
gamedata = build_robust_doe(simdata)
results_robust_doe, t_robust_doe = @timed pmap(i -> calc_run(design[i:i, :], gamedata, simdata, i), 1:nruns)
fn = joinpath(dir, "data\\results_robust_doe_temp.jld2")
@save fn results_robust_doe, t_robust_doe, gamedata

dfr = hcat(design, DataFrame(utility = [ri[1] for ri in results_robust_doe],
                ulb = [ri[2] for ri in results_robust_doe],
                uub = [ri[3] for ri in results_robust_doe],
                exploitability = [ri[4] for ri in results_robust_doe],
                exploitation = [ri[5] for ri in results_robust_doe],
                status = [ri[6] for ri in results_robust_doe],
                iter = [ri[7] for ri in results_robust_doe],
                time = [ri[8] for ri in results_robust_doe]))
fn = "data\\results_robust_doe_temp.csv"
CSV.write(joinpath(dir, fn), dfr)

max_stars_repo_name: ajkeith/Cyber-Air-Defense

max_issues_repo_name: ajkeith/Cyber-Air-Defense


text:
## module TestEconDatasets

using EconDatasets
using Base.Test
using DataFrames
using Dates

# write your own tests here
tests = [
         "readYahooFinance.jl",
         "runtests_ijulia.jl",
         "getDataset.jl", # time consuming
         "dataset.jl"
         ]

println("Running EconDatasets tests:")

for t in tests
    include(string(Pkg.dir("EconDatasets"), "/test/", t))
end

## end

ext: jl


text:
export RobotBasedMagneticFieldStaticProtocolParams, RobotBasedMagneticFieldStaticProtocol, measurement, filename

Base.@kwdef struct RobotBasedMagneticFieldStaticProtocolParams <: RobotBasedProtocolParams
  positions::Union{Positions, Missing} = missing
  postMoveWaitTime::typeof(1.0u"s") = 0.5u"s"
  numCooldowns::Integer = 0
  robotVelocity::typeof(1.0u"m/s") = 0.01u"m/s"
  switchBrakes::Bool = false
end
RobotBasedMagneticFieldStaticProtocolParams(dict::Dict) = createRobotBasedProtocolParams(RobotBasedMagneticFieldStaticProtocolParams, dict)

Base.@kwdef mutable struct RobotBasedMagneticFieldStaticProtocol <: RobotBasedProtocol
  name::AbstractString
  description::AbstractString
  scanner::MPIScanner
  params::RobotBasedMagneticFieldStaticProtocolParams

  biChannel::BidirectionalChannel{ProtocolEvent}
  done::Bool = false
  cancelled::Bool = false
  finishAcknowledged::Bool = false

  executeTask::Union{Task, Nothing} = nothing
  measurement::Union{MagneticFieldMeasurement, Missing} = missing
  filename::Union{AbstractString, Missing} = missing
  safetyTask::Union{Task, Nothing} = nothing
  safetyChannel::Union{Channel{ProtocolEvent}, Nothing} = nothing
end

measurement(protocol::RobotBasedMagneticFieldStaticProtocol) = protocol.measurement
measurement(protocol::RobotBasedMagneticFieldStaticProtocol, measurement::Union{MagneticFieldMeasurement, Missing}) = protocol.measurement = measurement
filename(protocol::RobotBasedMagneticFieldStaticProtocol) = protocol.filename
filename(protocol::RobotBasedMagneticFieldStaticProtocol, filename::String) = protocol.filename = filename

function init(protocol::RobotBasedMagneticFieldStaticProtocol)
  measurement_ = MagneticFieldMeasurement()
  MPIFiles.description(measurement_, "Generated by protocol $(name(protocol)) with the following description: $(description(protocol))")
  MPIFiles.positions(measurement_, positions(protocol))
  measurement(protocol, measurement_)
  # For inner protocol communication
  protocol.safetyChannel = Channel{ProtocolEvent}(32)
  # I'd prefer to only start task during execution and also close it in cleanup
  protocol.safetyTask = @tspawnat protocol.scanner.generalParams.protocolThreadID watchTemperature(protocol.scanner, protocol.safetyChannel)
  # For user I/O
  return BidirectionalChannel{ProtocolEvent}(protocol.biChannel)
end

function preMoveAction(protocol::RobotBasedMagneticFieldStaticProtocol, pos::Vector{<:Unitful.Length})
  if isready(protocol.safetyChannel)
    # Shut down system appropiately and then end execution
    close(protocol.safetyChannel)
    throw(IllegalStateException(take!(protocol.safetyChannel).message))
  end
  @info "moving to position" pos
end

function postMoveAction(protocol::RobotBasedMagneticFieldStaticProtocol, pos::Vector{<:Unitful.Length})
  gaussmeter = getGaussMeter(scanner(protocol))
  field_ = getXYZValues(gaussmeter)
  fieldError_ = calculateFieldError(gaussmeter, field_)
  fieldFrequency_ = getFrequency(gaussmeter)
  timestamp_ = now()
  temperature_ = getTemperature(gaussmeter)

  addMeasuredPosition(measurement(protocol), pos, field=field_, fieldError=fieldError_, fieldFrequency=fieldFrequency_, timestamp=timestamp_, temperature=temperature_)
end

function cleanup(protocol::RobotBasedMagneticFieldStaticProtocol)
  saveMagneticFieldAsHDF5(measurement(protocol), filename(protocol))
  close(protocol.scanner)
end

# Here I'd like to also dispatch on protocol and not only scanner
function watchTemperature(scanner::MPIScanner, channel::Channel)
  while isopen(channel)
    temp = getTemperature(scanner)
    if temp > 1.0
      put!(channel, IllegaleStateEvent("Temperature is too high. Aborting protocol"))
    end
    sleep(0.05)
  end
end

function handleEvent(protocol::RobotBasedMagneticFieldStaticProtocol, event::FinishedAckEvent)
  # End SafetyTask
  close(protocol.safetyChannel)
end


text:
## Main function to run Bayesian OPTIMISATION
function settingsBayesOptMCBU(oedmc_def)

    # Manage inputs order
    inpnam = oedmc_def["Model"]["inpName"];

    # Optimisation section
    steps = length(oedmc_def["switchT"])-1;
    induc = oedmc_def["Model"]["nInp"];
    fiii = length(oedmc_def["fixedInp"]);

    doof = (steps*(induc-fiii))+fiii;

    # THE THIRD MOUNTAINS!
    uppe = Array{Any,1}(undef, doof); # Vector of upper bounds. Length set with doof but deppending on the user settings from fixedStep and equalStep this can change.
    lowe = Array{Any,1}(undef, doof); # Vector of lower bounds. Same ^
    coun = 1; # Index for the entry of ustr that will advance only if certain conditions are meet (that that position has been already checked and a variable or number has been introduced)
                # Tecnichally this indicated the first entry of uppe and lowwe for each inducer variables.
    fs = [oedmc_def["fixedStep"][j][1] for j in 1:length(oedmc_def["fixedStep"])]; # Fixed Steps (the actual step number)
    us = [oedmc_def["fixedStep"][j][2] for j in 1:length(oedmc_def["fixedStep"])]; # Input values (or Any entry) for each step specified in fixedStep
    fixin = 0; # Index to account for fixedInp. If the fixedInp is the first one in the inpam list this is the best way (that I have found) to account for it. This is because of the use of i (input index). This is because first I consider the case where there is a fisedInp and then check the rest.

    for i in 1:length(inpnam) # Loop over the number of inducers
        if inpnam[i] in oedmc_def["fixedInp"] # If the current inducer is fixed then add only one bound and move index (coun). Fixin works as described in the other MOUNTAINS
            uppe[coun] = oedmc_def["uUpper"][i];
            lowe[coun] = oedmc_def["uLower"][i];
            coun += 1;
            fixin += 1;
        else
            try # Using coun you take all the indexes corresponding to the current inducer. The you fill it with the specified uper and lower bounds. Added the try just in case someoen runns it in a machien that does not support Float 64? (I am not sure)
                uppe[coun:coun+(steps-1)] = repeat([convert(Float64, oedmc_def["uUpper"][i])], outer = steps);
                lowe[coun:coun+(steps-1)] = repeat([convert(Float64, oedmc_def["uLower"][i])], outer = steps);
            catch
                uppe[coun:coun+(steps-1)] = repeat([convert(Float32, oedmc_def["uUpper"][i])], outer = steps);
                lowe[coun:coun+(steps-1)] = repeat([convert(Float32, oedmc_def["uLower"][i])], outer = steps);
            end

            for j in 1:length(uppe[coun:coun+(steps-1)]) # Now we loop over each entry that we filled to check if we need to delete it according to the users fixedStep and equalStep
                if j in fs # Since this is not a fixed input j also indicates the step index, so check if this is in the list of fixedStep
                    for w in 1:length(fs) # As commentetd in other loops, this for and if is to consider the case the user does not introduce fixedStep entries sorted by step.
                        if j == fs[w]
                            if us[w][i-fixin] != Any # Check if current input contained in a fixedStep is defined as Any (to be optimised) or not (has been fixed. )
                                uppe[coun-1+j] = "";# Since it was fixed (different from Any) ad an empty string there to remove later in both bounds.
                                lowe[coun-1+j] = "";
                            end
                        end
                    end
                end
                if oedmc_def["equalStep"] != [] # If equalStep is empty then we do not need to do anything else, if it has something then there are entries to remove from the bounds
                    for k in 1:length(oedmc_def["equalStep"]) # Loop over each entry of equalStep
                        tmpk = zeros(length(oedmc_def["equalStep"][k])); # Temporal variable that will contain the same values as the equalStep entry except the ones that have been fixed for that step (this will happen if the user set an Any for the current inducer)
                        for h in 1:length(oedmc_def["equalStep"][k]) # Loop over each value in the current entry
                            if oedmc_def["equalStep"][k][h] in fs # Check if the current step considered is also in fixedStep
                                for w in 1:length(fs) # Account for ht euser not introducng steps sorted
                                    if oedmc_def["equalStep"][k][h] == fs[w]
                                        if us[w][i-fixin] == Any # If the current step considered was defined as Any for the current inducer add the index, if it was fixed to a value ignore (this has been taken into account before)
                                            tmpk[h] = oedmc_def["equalStep"][k][h];
                                        end
                                    end
                                end
                            else
                                tmpk[h] = oedmc_def["equalStep"][k][h] # If the current step was not in fixedStep just add it
                            end
                        end
                        tmpk = convert.(Int, filter!(e->e∉[0],tmpk)); # Remove 0 entries (fixed values)
                        if j in tmpk && j != tmpk[1] # Now, check if the current step is in our temporary index. If it is the first entry in the vector then do nothing, if not remove (it will be the same as the first, so no bound for it)
                            uppe[coun-1+j] = "";
                            lowe[coun-1+j] = "";
                        end
                    end
                end
            end
            coun += (steps); # Add the number of steps to the varable coun so we can move to the next inducer
        end
    end

    uppe = [uppe[i] for i in 1:length(uppe) if uppe[i] != ""]; # Remove empty strings so we are left with a vector of numbers with the bounds for the inducers.
    lowe = [lowe[i] for i in 1:length(lowe) if lowe[i] != ""];

    #################################### OPTIMISER SETINGS THAT CAN BE MODIFIED START HERE ####################################

    model = ElasticGPE(convert(Int, length(uppe)),                            #
           mean = MeanConst(0.),
           kernel = SEArd(zeros(length(uppe)), 5.),
           logNoise = 0.,
           capacity = 3000);

    modeloptimizer = MAPGPOptimizer(every = 50, noisebounds = [-4, 3],       # bounds of the logNoise
                   kernbounds = [[-1*ones(length(uppe)); 0], [4*ones(length(uppe)); 10]],  # bounds of the 3 parameters GaussianProcesses.get_param_names(model.kernel)
                   maxeval = 40);

    sut = Symbol(string(oedmc_def["Model"]["NameF"],"Utility"));
    opt = BOpt((@eval $sut),
          model,
          ExpectedImprovement(),                   # type of acquisition
          modeloptimizer,
          lowe,
          uppe,                                     # lowerbounds, upperbounds
          repetitions = 1,                          # evaluate the function for each input 1 times
          maxiterations = oedmc_def["maxiter"],     # evaluate at 50 input positions
          sense = Max,                              # maximise the function
          verbosity = Progress);

    #################################### OPTIMISER SETINGS THAT CAN BE MODIFIED STOP HERE ####################################

    return opt

end

max_stars_repo_name: DavidGomezC/BOMBSbeta

max_issues_repo_name: DavidGomezC/BOMBSbeta

max_forks_repo_name: DavidGomezC/BOMBSbeta


text:
"""
Combine a feature along with the number of scenario successes and failures.
"""
struct FeatureSuccessAndFailure
    feature::Feature
    n_success::UInt
    n_failure::UInt
end

"""
Accumulate results from executed features as they are being executed. Keep track of whether the
total run is a success of a failure.
"""
mutable struct ResultAccumulator
    isaccumsuccess::Bool
    features::Vector{FeatureSuccessAndFailure}

    ResultAccumulator() = new(true, [])
end

"""
    accumulateresult!(acc::ResultAccumulator, result::FeatureResult)

Check for success or failure in this feature result and update the accumulator accordingly.
"""
function accumulateresult!(acc::ResultAccumulator, result::FeatureResult)
    n_success::UInt = 0
    n_failure::UInt = 0
    # Count the number of successes and failures for each step in each scenario. A scenario is
    # successful if all its steps are successful.
    for scenarioresult in result.scenarioresults
        arestepssuccessful = [issuccess(step) for step in scenarioresult.steps]
        isscenariosuccessful = all(arestepssuccessful)
        if isscenariosuccessful
            n_success += 1
        else
            n_failure += 1
        end

        acc.isaccumsuccess &= isscenariosuccessful
    end

    push!(acc.features, FeatureSuccessAndFailure(result.feature, n_success, n_failure))
end

"""
    issuccess(acc::ResultAccumulator)

True if all scenarios in all accumulated features are successful.
"""
issuccess(acc::ResultAccumulator) = acc.isaccumsuccess

"""
    featureresults(accumulator::ResultAccumulator)

A public getter for the results of all features.
"""
function featureresults(accumulator::ResultAccumulator)
    accumulator.features
end


text:
function img_similarity(pathtoimage)
    imagelist = glob("*.jpg", pathtoimage)
    append!(imagelist,  glob("*.JPG", pathtoimage))
    append!(imagelist,  glob("*.png", pathtoimage))
    append!(imagelist,  glob("*.PNG", pathtoimage))

    temp = Array{Array{Float64,2},1}(undef,length(imagelist))
    @showprogress "Loading Images " for i in 1:length(imagelist)
        temp[i] = Float64.(Gray.(imresize(load(imagelist[i]),128,128)))
    end

    similarity_matrix = Matrix{Float64}(undef, length(temp), length(temp))
    for i in 1:length(temp)
        for j in 1:length(temp)
            similarity_matrix[j,i] = assess_ssim(temp[j],temp[i])
        end
    end
    return similarity_matrix
end


text:
#!/usr/local/bin/julia

using ArgParse

function readFile(filepath)
	Sequences = []
	for Sequence in eachline(filepath)
		push!(Sequences, Sequence)
	end

	return Sequences
end

HammingDistance(Seq₁, Seq₂) = sum(collect(Seq₁) .!= collect(Seq₂))

function main(args)
	# Argument parsing
	configs = ArgParseSettings()

	add_arg_group!(configs, "Options:")
	@add_arg_table! configs begin
		"--input", "-i"
			help	 = "Input filename"
			action	 = :store_arg
			arg_type = String
		"--output", "-o"
			help	 = "Output filename"
			action	 = :store_arg
			arg_type = String
		end
	# Store arguments to variables
	parsed_args = parse_args(args, configs)
	infile      = parsed_args["input"]
	outfile     = parsed_args["output"]

	# Read file and compute problem via function
	@show s,t = readFile(infile)
	@show output = HammingDistance(s,t)
	open(f->write(f, "$output"), outfile, "w+")

end

main(ARGS)


text:
"""
    WrappedGP{Tgp<:AbstractGP} <: SthenoAbstractGP

A thin wrapper around an AbstractGP that does some book-keeping.

```julia
f = wrap(GP(SEKernel()), GPC())
```
builds a `WrappedGP` that `Stheno` knows how to work with.
"""
struct WrappedGP{Tgp<:AbstractGP} <: SthenoAbstractGP
    gp::Tgp
    n::Int
    gpc::GPC
    function WrappedGP{Tgp}(gp::Tgp, gpc::GPC) where {Tgp<:GP}
        wgp = new{Tgp}(gp, next_index(gpc), gpc)
        gpc.n += 1
        return wgp
    end
end

wrap(gp::Tgp, gpc::GPC) where {Tgp<:GP} = WrappedGP{Tgp}(gp, gpc)

mean(f::WrappedGP, x::AbstractVector) = mean(f.gp, x)

cov(f::WrappedGP, x::AbstractVector) = cov(f.gp, x)
var(f::WrappedGP, x::AbstractVector) = var(f.gp, x)

cov(f::WrappedGP, x::AbstractVector, x′::AbstractVector) = cov(f.gp, x, x′)
var(f::WrappedGP, x::AbstractVector, x′::AbstractVector) = var(f.gp, x, x′)

function cov(f::WrappedGP, f′::WrappedGP, x::AbstractVector, x′::AbstractVector)
    return f === f′ ? cov(f, x, x′) : zeros(length(x), length(x′))
end
function var(f::WrappedGP, f′::WrappedGP, x::AbstractVector, x′::AbstractVector)
    return f === f′ ? var(f, x, x′) : zeros(length(x))
end


text:
function unsafe_update_v1_to_v2!(v1::Vector{T}, v2::Vector{T}, n::Int) where T<:Real
    @inbounds @simd for i = 1:n
        v1[i] = v2[i]
    end
end

function unsafe_update_vv1_to_vv2!(vv1::Vector{Vector{T}}, vv2::Vector{Vector{T}}, n1::Int , n2::Int) where T<:Real

    @inbounds for i = 1:n1
        @simd for j = 1:n2
            vv1[i][j] = vv2[i][j]
        end
    end
end

function my_copy(vec::Vector{T})::Vector{T} where T<:Real
    n = length(vec)
    new_vec = Vector{T}(undef, n)

    @inbounds @simd for i = 1:n
        new_vec[i] = vec[i]
    end

    return new_vec
end

# test
# a = [[rand(1:5) for j in 1:5] for i in 1:200]
function my_deep_copy_unsafe(vec_vec::Vector{Vector{T}}, n1::Int , n2::Int)::Vector{Vector{T}} where T<:Real
    new_vec = Vector{Vector{T}}(undef, n1)

    @inbounds for i = 1:n1
        new_vec[i] = Vector{T}(undef, n2)

        @simd for j = 1:n2
            new_vec[i][j] = vec_vec[i][j]
        end
    end

    return new_vec
end


text:
using GPUCompiler

if !@isdefined(TestRuntime)
    include("../util.jl")
end

# create a native test compiler, and generate reflection methods for it

NativeCompilerJob = CompilerJob{NativeCompilerTarget,TestCompilerParams}

# local method table for device functions
@static if isdefined(Base.Experimental, Symbol("@overlay"))
Base.Experimental.@MethodTable(method_table)
else
const method_table = nothing
end

GPUCompiler.method_table(@nospecialize(job::NativeCompilerJob)) = method_table

function native_job(@nospecialize(func), @nospecialize(types); kernel::Bool=false, kwargs...)
    source = FunctionSpec(func, Base.to_tuple_type(types), kernel)
    target = NativeCompilerTarget(always_inline=true)
    params = TestCompilerParams()
    CompilerJob(target, source, params), kwargs
end

function native_code_typed(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kwargs...)
    GPUCompiler.code_typed(job; kwargs...)
end

function native_code_warntype(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kwargs...)
    GPUCompiler.code_warntype(io, job; kwargs...)
end

function native_code_llvm(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kwargs...)
    GPUCompiler.code_llvm(io, job; kwargs...)
end

function native_code_native(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kwargs...)
    GPUCompiler.code_native(io, job; kwargs...)
end

# aliases without ::IO argument
for method in (:code_warntype, :code_llvm, :code_native)
    native_method = Symbol("native_$(method)")
    @eval begin
        $native_method(@nospecialize(func), @nospecialize(types); kwargs...) =
            $native_method(stdout, func, types; kwargs...)
    end
end

# simulates codegen for a kernel function: validates by default
function native_code_execution(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = native_job(func, types; kernel=true, kwargs...)
    GPUCompiler.compile(:asm, job; kwargs...)
end

module LazyCodegen
    using LLVM
    using LLVM.Interop
    using GPUCompiler

    import ..native_job

    @static if Base.libllvm_version < v"12"
        # We have one global JIT and TM
        const jit = Ref{LLVM.OrcJIT}()
        const tm  = Ref{LLVM.TargetMachine}()

        function __init__()
            optlevel = LLVM.API.LLVMCodeGenLevelDefault

            tm[] = GPUCompiler.JITTargetMachine(optlevel=optlevel)
            LLVM.asm_verbosity!(tm[], true)

            jit[] = LLVM.OrcJIT(tm[]) # takes ownership of tm
            atexit() do
                dispose(jit[])
            end
        end

        import GPUCompiler: deferred_codegen_jobs, CompilerJob
        mutable struct CallbackContext
            job::CompilerJob
            stub::String
            compiled::Bool
        end

        const outstanding = IdDict{CallbackContext, Nothing}()

        # Setup the lazy callback for creating a module
        function callback(orc_ref::LLVM.API.LLVMOrcJITStackRef, callback_ctx::Ptr{Cvoid})
            orc = LLVM.OrcJIT(orc_ref)
            cc = Base.unsafe_pointer_to_objref(callback_ctx)::CallbackContext

            @assert !cc.compiled
            job = cc.job

            ir, meta = GPUCompiler.codegen(:llvm, job; validate=false)
            entry_name = name(meta.entry)

            jitted_mod = compile!(orc, ir)

            addr = addressin(orc, jitted_mod, entry_name)
            ptr  = pointer(addr)

            cc.compiled = true
            delete!(outstanding, cc)

            # 4. Update the stub pointer to point to the recently compiled module
            set_stub!(orc, cc.stub, ptr)

            # 5. Return the address of the implementation, since we are going to call it now
            return UInt64(reinterpret(UInt, ptr))
        end

        function get_trampoline(job)
            cc = CallbackContext(job, String(gensym(:trampoline)), false)
            outstanding[cc] = nothing

            c_callback = @cfunction(callback, UInt64, (LLVM.API.LLVMOrcJITStackRef, Ptr{Cvoid}))

            orc = jit[]
            initial_addr = callback!(orc, c_callback, pointer_from_objref(cc))
            create_stub!(orc, cc.stub, initial_addr)
            return address(orc, cc.stub)
        end
    else # LLVM >=12

        struct CompilerInstance
            jit::LLVM.LLJIT
            lctm::LLVM.LazyCallThroughManager
            ism::LLVM.IndirectStubsManager
        end
        const jit = Ref{CompilerInstance}()

        function __init__()
            optlevel = LLVM.API.LLVMCodeGenLevelDefault
            tm = GPUCompiler.JITTargetMachine(optlevel=optlevel)
            LLVM.asm_verbosity!(tm, true)

            lljit = LLJIT(;tm)

            jd_main = JITDylib(lljit)

            prefix = LLVM.get_prefix(lljit)
            dg = LLVM.CreateDynamicLibrarySearchGeneratorForProcess(prefix)
            add!(jd_main, dg)

            es = ExecutionSession(lljit)

            lctm = LLVM.LocalLazyCallThroughManager(triple(lljit), es)
            ism = LLVM.LocalIndirectStubsManager(triple(lljit))

            jit[] = CompilerInstance(lljit, lctm, ism)
            atexit() do
                ci = jit[]
                dispose(ci.ism)
                dispose(ci.lctm)
                dispose(ci.jit)
            end
        end

        function get_trampoline(job)
            compiler = jit[]
            lljit = compiler.jit
            lctm  = compiler.lctm
            ism   = compiler.ism

            # We could also use one dylib per job
            jd = JITDylib(lljit)

            entry_sym = String(gensym(:entry))
            target_sym = String(gensym(:target))
            flags = LLVM.API.LLVMJITSymbolFlags(
                LLVM.API.LLVMJITSymbolGenericFlagsCallable |
                LLVM.API.LLVMJITSymbolGenericFlagsExported, 0)
            entry = LLVM.API.LLVMOrcCSymbolAliasMapPair(
                mangle(lljit, entry_sym),
                LLVM.API.LLVMOrcCSymbolAliasMapEntry(
                    mangle(lljit, target_sym), flags))

            mu = LLVM.reexports(lctm, ism, jd, Ref(entry))
            LLVM.define(jd, mu)

            # 2. Lookup address of entry symbol
            addr = lookup(lljit, entry_sym)

            # 3. add MU that will call back into the compiler
            sym = LLVM.API.LLVMOrcCSymbolFlagsMapPair(mangle(lljit, target_sym), flags)

            function materialize(mr)
                ir, meta = GPUCompiler.codegen(:llvm, job; validate=false)

                # Rename entry to match target_sym
                LLVM.name!(meta.entry, target_sym)

                # So 1. serialize the module
                buf = convert(MemoryBuffer, ir)

                # 2. deserialize and wrap by a ThreadSafeModule
                ctx = ThreadSafeContext()
                mod = parse(LLVM.Module, buf; ctx=context(ctx))
                tsm = ThreadSafeModule(mod; ctx)

                il = LLVM.IRTransformLayer(lljit)
                LLVM.emit(il, mr, tsm)

                return nothing
            end

            function discard(jd, sym)
            end

            mu = LLVM.CustomMaterializationUnit(entry_sym, Ref(sym), materialize, discard)
            LLVM.define(jd, mu)
            return addr
        end
    end

    import GPUCompiler: deferred_codegen_jobs
    @generated function deferred_codegen(::Val{f}, ::Val{tt}) where {f,tt}
        job, _ = native_job(f, tt)

        addr = get_trampoline(job)
        trampoline = pointer(addr)
        id = Base.reinterpret(Int, trampoline)

        deferred_codegen_jobs[id] = job

        quote
            ptr = ccall("extern deferred_codegen", llvmcall, Ptr{Cvoid}, (Ptr{Cvoid},), $trampoline)
            assume(ptr != C_NULL)
            return ptr
        end
    end

    @generated function abi_call(f::Ptr{Cvoid}, rt::Type{RT}, tt::Type{T}, args::Vararg{Any, N}) where {T, RT, N}
        argtt    = tt.parameters[1]
        rettype  = rt.parameters[1]
        argtypes = DataType[argtt.parameters...]

        argexprs = Union{Expr, Symbol}[]
        ccall_types = DataType[]

        before = :()
        after = :(ret)

        # Note this follows: emit_call_specfun_other
        Context() do ctx
            T_jlvalue = LLVM.StructType(LLVMType[],;ctx)
            T_prjlvalue = LLVM.PointerType(T_jlvalue, #= AddressSpace::Tracked =# 10)

            for (source_i, source_typ) in enumerate(argtypes)
                if GPUCompiler.isghosttype(source_typ) || Core.Compiler.isconstType(source_typ)
                    continue
                end

                argexpr = :(args[$source_i])

                isboxed = GPUCompiler.deserves_argbox(source_typ)
                et = isboxed ? T_prjlvalue : convert(LLVMType, source_typ; ctx)

                if isboxed
                    push!(ccall_types, Any)
                elseif isa(et, LLVM.SequentialType) # et->isAggregateType
                    push!(ccall_types, Ptr{source_typ})
                    argexpr = Expr(:call, GlobalRef(Base, :Ref), argexpr)
                else
                    push!(ccall_types, source_typ)
                end
                push!(argexprs, argexpr)
            end

            if GPUCompiler.isghosttype(rettype) || Core.Compiler.isconstType(rettype)
                # Do nothing...
                # In theory we could set `rettype` to `T_void`, but ccall will do that for us
            # elseif jl_is_uniontype?
            elseif !GPUCompiler.deserves_retbox(rettype)
                rt = convert(LLVMType, rettype; ctx)
                if !isa(rt, LLVM.VoidType) && GPUCompiler.deserves_sret(rettype, rt)
                    before = :(sret = Ref{$rettype}())
                    pushfirst!(argexprs, :(sret))
                    pushfirst!(ccall_types, Ptr{rettype})
                    rettype = Nothing
                    after = :(sret[])
                end
            else
                # rt = T_prjlvalue
            end
        end

        quote
            $before
            ret = ccall(f, $rettype, ($(ccall_types...),), $(argexprs...))
            $after
        end
    end

    @inline function call_delayed(f::F, args...) where F
        tt = Tuple{map(Core.Typeof, args)...}
        rt = Core.Compiler.return_type(f, tt)
        ptr = deferred_codegen(Val(f), Val(tt))
        abi_call(ptr, rt, tt, args...)
    end
end


text:
using Documenter
using Literate
using LatticeTools

example_directory = joinpath(@__DIR__, "..", "examples")
output_directory = joinpath(@__DIR__, "src/generated")
for filename in [
    "example_group_isomorphism.jl",
    "example_honeycomb_symmetry.jl",
    "example_hypercube.jl",
    "example_kagome_symmetry.jl",
    "example_little_group_2d.jl",
    "example_little_group_3d.jl",
    "example_little_symmetry_kspace_honeycomb.jl",
    "example_point_group_4mm.jl",
    "example_point_group_elements.jl",
]
    example = joinpath(example_directory, filename)
    Literate.markdown(example, output_directory)
end

makedocs(
    modules=[LatticeTools],
    doctest=true,
    sitename="LatticeTools.jl",
    format=Documenter.HTML(prettyurls=!("local" in ARGS)),
    authors="Kyungmin Lee",
    checkdocs=:all,
    pages = [
        "Home" => "index.md",
        "lattice.md",
        "Symmetry Analysis" => [
            "Symmetry/space-symmetry.md",
            "Symmetry/space-operation.md",
            "Symmetry/irrep.md",
            "Symmetry/space-symmetry-database.md",
        ],
        "Examples" => [
            "Group Isomorphism" => "generated/example_group_isomorphism.md",
            "Hypercube" => "generated/example_hypercube.md",
            "Point Group Elements" => "generated/example_point_group_elements.md",
            "Point Group 4mm (C₄ᵥ)" => "generated/example_point_group_4mm.md",
            "Little Group in 2D" => "generated/example_little_group_2d.md",
            "Little Group in 3D" => "generated/example_little_group_3d.md",
            "Honeycomb Lattice" => "generated/example_honeycomb_symmetry.md",
            "Honeycomb Lattice in k-space" => "generated/example_little_symmetry_kspace_honeycomb.md",
            "Kagome Lattice" => "generated/example_kagome_symmetry.md",
        ],
        "API Reference" => [
            "Basic" => "API/basic.md",
            "Symmetry Operation" => "API/symmetry-operation.md",
            "Symmetry" => "API/symmetry.md",
            "Symmetry Embedding" => "API/symmetry-embedding.md",
            "Symmetry Irrep Component" => "API/symmetry-irrep-component.md",
        ]
    ]
)

deploydocs(;
    repo="github.com/kyungminlee/LatticeTools.jl.git",
    devbranch="dev",
)


text:
const imaginary_unit_str =  split(string(Complex(0,0)), "0")[end]

const ARB_STR_FLAGS                = UInt64(0)

const ARB_STR_MORE                 = UInt64(1)        # extra, maybe incorrect digs
const ARB_STR_NO_RADIUS            = UInt64(2)        # omit radius
const ARB_STR_RADIUS               = UInt64(4)
const ARB_STR_CONDENSE             = UInt64(16)
const ARB_STR_CONDENSE_36_Digits   = UInt64(50*16)
const ARB_STR_CONDENSE_72_Digits   = UInt64(72*16)

const NO_FLAGS  = UInt64(0)
const NO_RADIUS = ARB_STR_NO_RADIUS

@inline function digit_precision(bitprecision::Int)
    bitprecision = evincedbits(bitprecision)
    return maximin_digits(bitprecision)
end

function trimzeros(str::String)
    m =  n = length(str)
    n === 0 && return str
    while n>1 && str[n] === '0'
        n -= 1
    end
    m === n ? str : str[1:n+1]
end

function string(x::Mag, maxdigits::Int = maximin_digits(30), flags::UInt = NO_FLAGS)
    y = ArbFloat{precision(x)}()
    ccall(@libarb(arf_set_mag), Cvoid, (Ref{ArbFloat}, Ref{Mag}), y, x)
    z = ArbReal{precision(x)}()
    ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), z, y)
    unsafestr = ccall(@libarb(arb_get_str), Cstring,
                      (Ref{ArbReal}, Clong, Culong), z, maxdigits, flags)
    str = deepcopy( unsafe_string(pointer(unsafestr)) )
    ccall(@libflint(flint_free), Cvoid, (Cstring,), unsafestr)
    return str
end

function string(x::ArbFloat{P}; midpoint::Bool=false) where {P}
    prec = midpoint ? digits4bits(P) : digit_precision(P)
    flags = ARB_STR_NO_RADIUS
    return arbstring(x, prec, flags=flags)
end

stringall(x::ArbFloat{P}) where {P} = string(x, midpoint=true)

function arbstring(x::ArbFloat{P}, maxdigits::Int=digit_precision(P); flags::UInt = ARB_STR_NO_RADIUS) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), z, x)
    unsafestr = ccall(@libarb(arb_get_str), Cstring,
                      (Ref{ArbReal}, Clong, Culong), z, maxdigits, flags)
    str = deepcopy( unsafe_string(pointer(unsafestr)) )
    str = trimzeros(str)
    ccall(@libflint(flint_free), Cvoid, (Cstring,), unsafestr)
    return str
end

function string(x::ArbReal{P}; midpoint::Bool=false, radius::Bool=false) where {P}
    prec = midpoint ? digits4bits(P) : digit_precision(P)
    flags = radius ? ARB_STR_RADIUS : ARB_STR_NO_RADIUS
    return arbstring(x, prec, flags=flags)
end

stringall(x::ArbReal{P}; radius::Bool=false) where {P} = string(x, midpoint=true, radius=radius)

function arbstring(x::ArbReal{P}, maxdigits::Int=digit_precision(P); flags::UInt = NO_FLAGS) where {P}
    unsafestr = ccall(@libarb(arb_get_str), Cstring,
                      (Ref{ArbReal}, Clong, Culong), x, maxdigits, flags)
    str = deepcopy( unsafe_string(pointer(unsafestr)) )
    str = trimzeros(str)
    ccall(@libflint(flint_free), Cvoid, (Cstring,), unsafestr)
    return str
end

function string(x::ArbComplex{P}; midpoint::Bool=false, radius::Bool=false) where {P}
    prec = midpoint ? digits4bits(P) : digit_precision(P)
    flags = radius ? ARB_STR_RADIUS : ARB_STR_NO_RADIUS
    return arbstring(x, prec, flags=flags)
end

stringall(x::ArbComplex{P}; radius::Bool=false) where {P} = string(x, midpoint=true, radius=radius)

function arbstring(x::ArbComplex{P}, maxdigits::Int=digit_precision(P); flags::UInt = NO_FLAGS) where {P}
    # rea, ima = real(x), imag(x)
    rea = ArbReal{P}()
    ima = ArbReal{P}()
    ccall(@libarb(acb_get_real), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}), rea, x)
    ccall(@libarb(acb_get_imag), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}), ima, x)

    ima_isneg = signbit(ima)
    ima_abs = ima_isneg ? -ima : ima
    connection = ima_isneg ? " - " : " + "
    rea_str = arbstring(rea, maxdigits, flags=flags)
    ima_str = arbstring(abs(ima), maxdigits, flags=flags)
    str = string(rea_str, connection, ima_str, imaginary_unit_str)
    return str
end

#= from Nemo.jl
function show(io::IO, x::arb)
  d = ceil(parent(x).prec * 0.30102999566398119521)
  cstr = ccall((:arb_get_str, :libarb), Ptr{UInt8}, (Ref{arb}, Int, UInt),
                                                  x, Int(d), UInt(0))
  print(io, unsafe_string(cstr))
  ccall((:flint_free, :libflint), Void, (Ptr{UInt8},), cstr)
end
=#


text:
#arrowhead.jl

using JuMP
using Ipopt

N = parse(Int, ARGS[1])
K = parse(Int, ARGS[2])

if parse(Int, ARGS[3]) == 1
    m = Model(solver=IpoptSolver())
elseif parse(Int, ARGS[3]) == 2
    include("interface_tape.jl")
    using TapeInterface
    m = Model(solver=TapeSolver(IpoptSolver()))
elseif parse(Int, ARGS[3]) == 3
    include("interface_tape_dict.jl")
    using TapeInterface
    m = Model(solver=TapeSolver(IpoptSolver()))
else
    @assert false
end

@variable(m, -100<=x[1:N+K]<=100)

@NLobjective(m, Min, sum{ cos(sum{ x[i+j], j=1:K}) + sum{ (x[i] + x[j])^2, j=1:K} , i=1:N} )

solve(m)


text:
@testset "testing interoperability, matrix-based net" begin

# tree, some edge lengths missing
tree1 = readTopology("(A,(B:1.0,(C:1.0,D:1.0):1.0):1.0);");
@test_nowarn PhyloNetworks.resetNodeNumbers!(tree1);
@test PhyloNetworks.majoredgematrix(tree1) == [5 1; 5 6; 6 2; 6 7; 7 3; 7 4]
@test all(PhyloNetworks.majoredgelength(tree1) .===  [missing, 1.0,1.0,1.0,1.0,1.0])
@test PhyloNetworks.minorreticulationmatrix(tree1) == Array{Int64,2}(0,2)
@test size(PhyloNetworks.minorreticulationlength(tree1)) == (0,)
@test size(PhyloNetworks.minorreticulationgamma(tree1)) == (0,)

# network, h=1, some missing gamma values
net1 = readTopology("(((A:4.0,(B:1.0)#H1:1.1::0.9):0.5,(C:0.6,#H1:1.0):1.0):3.0,D:5.0);");
@test_nowarn PhyloNetworks.resetNodeNumbers!(net1);
@test PhyloNetworks.majoredgematrix(net1) == [5 6; 5 4; 6 8; 6 7; 7 3; 8 1; 8 9; 9 2]
@test PhyloNetworks.majoredgelength(net1) ==  [3.,5.,.5,1.,.6,4.,1.1,1.]
@test PhyloNetworks.minorreticulationmatrix(net1) == [7 9]
@test PhyloNetworks.minorreticulationlength(net1) == [1.]
@test all(PhyloNetworks.minorreticulationgamma(net1) .=== [missing])

# network, h=2 hybridizations
s = "(((Ag,(#H1:7.159::0.056,((Ak,(E:0.08,#H2:0.0::0.004):0.023):0.078,(M:0.0)#H2:::0.996):2.49):2.214):0.026,
      (((((Az:0.002,Ag2:0.023):2.11,As:2.027):1.697)#H1:0.0::0.944,Ap):0.187,Ar):0.723):5.943,(P,20):1.863,165);";
net2 = readTopology(s);
@test_nowarn PhyloNetworks.resetNodeNumbers!(net2);
@test PhyloNetworks.majoredgematrix(net2) == [13 15; 13 14; 13 12; 14 10; 14 11; 15 18;
  15 16; 16 17; 16 9; 17 24; 17 8; 18 1; 18 19; 19 20; 20 21; 20 23; 21 2; 21 22; 22 3;
  23 4; 24 25; 25 26; 25 7; 26 5; 26 6]
@test all(PhyloNetworks.majoredgelength(net2) .===
  [5.943,1.863,missing,missing,missing,.026,.723,.187,missing,0,missing,missing,
  2.214,2.490,.078,missing,missing,.023,.08,0,1.697,2.11,2.027,.002,.023])

#a.isnull[[3,4,5,9,11,12,16,17]] == [true for i in 1:8]
#@test convert(Array, a[[1,2,6,7,8,10,13,14,15]]) == [5.943,1.863, .026,.723,.187, 0, 2.214,2.490,.078]
#@test convert(Array, a[18:end]) == [.023,.08,0,1.697,2.11,2.027,.002,.023]
@test PhyloNetworks.minorreticulationmatrix(net2) == [19 24; 22 23]
@test PhyloNetworks.minorreticulationlength(net2) == [7.159,0.]
@test PhyloNetworks.minorreticulationgamma(net2) == [0.056,0.004]

end # of testset for interop


text:
module ElasticSearch

include("search.jl")

export ES, Query, search, test
export RangeBody, MultiMatchBody, SortBody, TermsBody
export MustQueryBlock, HighlightBlock, SourceBlock, FromBlock, SizeBlock,
        SortBlock, Block

end


text:
"""
`PDBFile <: Format`

Protein Data Bank (PDB) format.
It provides a standard representation for macromolecular
structure data derived from X-ray diffraction and NMR studies.
"""
immutable PDBFile <: Format end

"""
`parse(io, ::Type{PDBFile}; chain="all", model="all", group="all", atomname="all", onlyheavy=false, occupancyfilter=false)`

Reads a text file of a PDB entry.
Returns a list of `PDBResidue` (view `MIToS.PDB.PDBResidues`).
Setting `chain`, `model`, `group`, `atomname` and `onlyheavy` values
can be used to select of a subset of all residues. Group can be `ATOM`
or `HETATM`. If not set, all residues are returned.
If the keyword argument `occupancyfilter` (default: `false`) is `true`,
only the atoms with the best occupancy are returned.
"""
function parse(io::Union{IO, ASCIIString}, ::Type{PDBFile}; chain::ASCIIString = "all",
               model::ASCIIString = "all", group::ASCIIString = "all", atomname::ASCIIString="all",
               onlyheavy::Bool=false, occupancyfilter::Bool=false)
    residue_dict = OrderedDict{PDBResidueIdentifier, Vector{PDBAtom}}()
    atom_model = 0
    for line in eachline(io)
        line_id = length(line) < 6 ? replace(line, ' ', "") : replace(line[1:6], ' ', "") # i.e. "END\n"
        if line_id == "MODEL"
            atom_model += 1
        end
        if (group=="all" && (line_id=="ATOM" || line_id=="HETATM")) || (line_id == group)
            atom_chain = string(line[22])
            atom_name = replace(line[13:16],' ',"")
            element = replace(line[77:78],' ',"")
            if  (chain=="all" || chain==atom_chain) && ((model==0 || model=="all") || model==string(atom_model+1)) &&
                    (atomname=="all" || atomname==atom_name) && (!onlyheavy || element!="H")

                # 23 - 26        Integer         Residue sequence number.
                # 27             AChar           Code for insertion of residues.
                PDB_number = replace(line[23:27],' ',"")

                name = replace(line[18:20],' ',"")
                x = float(replace(line[31:38],' ',""))
                y = float(replace(line[39:46],' ',""))
                z = float(replace(line[47:54],' ',""))
                occupancy = float(replace(line[55:60],' ',""))
                B = replace(line[61:66],' ',"")

                mdl = atom_model == 0 ? "1" : string(atom_model)

                residue_id = PDBResidueIdentifier("", PDB_number, name, line_id, mdl, atom_chain)
                atom_data  = PDBAtom(Coordinates(x,y,z), atom_name, element, occupancy, B)

                value = get!(residue_dict, residue_id, PDBAtom[])
                push!(value, atom_data)

            end
        end
    end
    _generate_residues(residue_dict, occupancyfilter)
end

# Print PDB
# =========

# ATOM & HETATM
# COLUMNS        DATA TYPE       CONTENTS
# --------------------------------------------------------------------------------
#  1 -  6        Record name     "ATOM  "
#  7 - 11        Integer         Atom serial number.
# 13 - 16        Atom            Atom name.
# 17             Character       Alternate location indicator.
# 18 - 20        Residue name    Residue name.
# 22             Character       Chain identifier.
# 23 - 26        Integer         Residue sequence number.
# 27             AChar           Code for insertion of residues.
# 31 - 38        Real(8.3)       Orthogonal coordinates for X in Angstroms.
# 39 - 46        Real(8.3)       Orthogonal coordinates for Y in Angstroms.
# 47 - 54        Real(8.3)       Orthogonal coordinates for Z in Angstroms.
# 55 - 60        Real(6.2)       Occupancy.
# 61 - 66        Real(6.2)       Temperature factor (Default = 0.0).
# 73 - 76        LString(4)      Segment identifier, left-justified.
# 77 - 78        LString(2)      Element symbol, right-justified.
# 79 - 80        LString(2)      Charge on the atom.
# Example:
#          1         2         3         4         5         6         7         8
# 12345678901234567890123456789012345678901234567890123456789012345678901234567890
# ATOM    145  N   VAL A  25      32.433  16.336  57.540  1.00 11.92      A1   N
# ATOM    146  CA  VAL A  25      31.132  16.439  58.160  1.00 11.85      A1   C
# ATOM    147  C   VAL A  25      30.447  15.105  58.363  1.00 12.34      A1   C
# ATOM    148  O   VAL A  25      29.520  15.059  59.174  1.00 15.65      A1   O
# ATOM    149  CB AVAL A  25      30.385  17.437  57.230  0.28 13.88      A1   C
# ATOM    150  CB BVAL A  25      30.166  17.399  57.373  0.72 15.41      A1   C
# ATOM    151  CG1AVAL A  25      28.870  17.401  57.336  0.28 12.64      A1   C
# ATOM    152  CG1BVAL A  25      30.805  18.788  57.449  0.72 15.11      A1   C
# ATOM    153  CG2AVAL A  25      30.835  18.826  57.661  0.28 13.58      A1   C
# ATOM    154  CG2BVAL A  25      29.909  16.996  55.922  0.72 13.25      A1   C
#
# HETATM 1357 MG    MG   168       4.669  34.118  19.123  1.00  3.16          MG2+
# HETATM 3835 FE   HEM     1      17.140   3.115  15.066  1.00 14.14          FE3+

const _Format_PDB_ATOM = FormatExpr(
    #          1         2         3         4         5         6         7         8
    # 12345678901234567890123456789012345678901234567890123456789012345678901234567890
    # >    <>   < >  <|> < |>  <|   >      <>      <>      <>    <>    <      >  <><><
    "{:<6}{:>5d} {:<4}{:>1}{:>3} {:>1}{:>4}{:>1}   {:>8.3f}{:>8.3f}{:>8.3f}{:>6.2f}{:>6.2f}      {:<4}{:>2}{:>2}\n"
    )

# Models are numbered sequentially beginning with 1.
# Each MODEL must have a corresponding ENDMDL record.
# COLUMNS        DATA TYPE       CONTENTS
# --------------------------------------------------------------------------------
# 1 -  6       Record name    "MODEL "
# 11 - 14       Integer        Model serial number
# Example:
#          1         2         3         4         5         6         7         8
# 12345678901234567890123456789012345678901234567890123456789012345678901234567890
# MODEL        1
# ATOM      1  N   ALA     1      11.104   6.134  -6.504  1.00  0.00           N
# ATOM    294 2HG  GLU    18     -13.630  -3.769   0.160  1.00  0.00           H
# TER     295      GLU    18
# ENDMDL

const _Format_PDB_MODEL = FormatExpr(
    #          1
    # 12345678901234
    # MODEL        1
     "MODEL     {:>4}\n"
    )

# TER
# Indicates the end of a list of ATOM/HETATM records for a chain
# The TER records occur in the coordinate section of the entry, and indicate
# the last residue presented for each polypeptide and/or nucleic acid chain for
# which there are coordinates.
# The TER record has the same residue name, chain identifier, sequence number
# and insertion code as the terminal residue. The serial number of the TER
# record is one number greater than the serial number of the ATOM/HETATM
# preceding the TER.
# The residue name appearing on the TER record must be the same as the residue name
# of the immediately preceding ATOM or non-water HETATM record.
# COLUMNS         DATA TYPE         CONTENTS
# --------------------------------------------------------------------------------
#  1 -  6         Record name       "TER   "
#  7 - 11         Integer           Serial number
# 18 - 20         Residue name      Residue name
# 22              Character         Chain identifier
# 23 - 26         Integer           Residue sequence number
# 27              AChar             Insertion code
# Example:
#          1         2         3         4         5         6         7         8
# 12345678901234567890123456789012345678901234567890123456789012345678901234567890
# ATOM   4150  H   ALA A 431       8.674  16.036  12.858  1.00  0.00           H
# TER    4151      ALA A 431

# HETATM 1415  O2  BLE P   1      13.775  30.147  14.862  1.09 20.95           O
# TER    1416      BLE P   1

const _Format_PDB_TER = FormatExpr(
    # TER    4151      ALA A 431
    #          1         2         3         4         5         6         7         8
    # 12345678901234567890123456789012345678901234567890123456789012345678901234567890
    "TER   {:>5d}      {:>3} {:>1}{:>4}{:>1}\n"
    )

function print(io::IO, res::PDBResidue, format::Type{PDBFile}, atom_index::Int, serial_number::Int)
    number = match(r"(\d+)(\D?)", res.id.number)
    atomname = res.atoms[atom_index].atom
    printfmt(io, _Format_PDB_ATOM,
             res.id.group,
             serial_number,
             length(atomname) <= 3 ? string(" ", atomname) : atomname, # It works with NACCESS
             " ",
             res.id.name,
             res.id.chain,
             number[1],
             number[2],
             res.atoms[atom_index].coordinates.x,
             res.atoms[atom_index].coordinates.y,
             res.atoms[atom_index].coordinates.z,
             res.atoms[atom_index].occupancy,
             res.atoms[atom_index].B,
             " ",
             res.atoms[atom_index].element,
             " ")
    serial_number + 1
end

function print(io::IO, res::PDBResidue, format::Type{PDBFile}, start::Int=1)
    next = start
    for i in 1:length(res.atoms)
        next = print(io, res, format, i, next)
    end
    nothing
end

function print(io::IO, reslist::AbstractVector{PDBResidue}, format::Type{PDBFile}, start::Int=1)
    next = start

    use_model = length(unique(map(res->res.id.model, reslist))) > 1
    if use_model
        model = "START"
    end

    for resindex in 1:length(reslist)
        res = reslist[resindex]

        # MODEL
        if use_model
            if model != res.id.model
                if model != "START"
                   println(io, "ENDMDL")
                end
                printfmt(io, _Format_PDB_MODEL, res.id.model)
            end
            model = res.id.model
        end

        # TER

        # MIToS only prints TER for the ATOM group if the chain changes.
        # Some modified residues are annotated as HETATM in the middle of the ATOM chain:
        # TER can not be printed from ATOM to HETATM if the chain doesn’t change.
        if resindex > 1
            previous_res = reslist[ resindex - 1 ]
            if (previous_res.id.group == "ATOM") && (previous_res.id.chain != res.id.chain)
                number = match(r"(\d+)(\D?)", previous_res.id.number)
                printfmt(io, _Format_PDB_TER,
                    next,
                    previous_res.id.name,
                    previous_res.id.chain,
                    number[1],
                    number[2]
                    )
                next += 1
            end
        end

        # ATOM/HETATM
        for i in 1:length(res.atoms)
            next = print(io, res, format, i, next)
        end

    end

    if use_model
        println(io, "ENDMDL")
    end

    println(io, "END")
    nothing
end

print(reslist::AbstractVector{PDBResidue}, format::Type{PDBFile}) = print(STDOUT, reslist, format)
print(res::PDBResidue, format::Type{PDBFile}) = print(STDOUT, res, format)

@doc """
`print(io, res, format::Type{PDBFile})`
`print(res, format::Type{PDBFile})`

Print a `PDBResidue` or a vector of `PDBResidue`s in PDB format.
""" print


text:
using DSP, Random
export AutoRegressive

"""
    AutoRegressive(n, method = LPCLevinson())

Autoregressive surrogates of order-`n`. The autoregressive coefficients `φ` are estimated
using `DSP.lpc(x, n, method)`, and thus see the documentation of DSP.jl for possible
`method`s.

While these surrogates are obviously suited to test the null hypothesis whether the data
are coming from a autoregressive process, the Fourier Transform-based surrogates are
probably a better option. The current method is more like an explicit way to
produce surrogates for the same hypothesis by fitting a model.
It can be used as convient way to estimate
autoregressive coefficients and automatically generate surrogates based on them.

The coefficients φ of the autoregressive fit can be found by doing
```julia
sg = surrogenerator(x, AutoRegressive(n))
φ = sg.init.φ
```
"""
struct AutoRegressive{M} <: Surrogate
    n::Int
    m::M
end
AutoRegressive(n::Int) = AutoRegressive(n, LPCLevinson())

function surrogenerator(x, method::AutoRegressive, rng = Random.default_rng())
    φ, e = lpc(x, method.n, method.m)
    N, f = length(x), length(φ)
    d = Normal(0, std(x))

    init = (
        d = d,
        φ = φ,
        N = N,
    )

    # The surrogate
    T = typeof(d).parameters[1]
    s = zeros(T, N + f)

    return SurrogateGenerator(method, x, s, init, rng)
end

function (sg::SurrogateGenerator{<:AutoRegressive})()
    N, φ, d = sg.init.N, sg.init.φ, sg.init.d
    return autoregressive!(sg.s, d, N, φ, sg.rng)
end

function autoregressive!(s, d::Normal{T}, N, φ, rng) where {T}
    f = length(φ)
    @inbounds for i in 1:f; s[i] = rand(rng, d); end
    @inbounds for i in f+1:length(s)
        s̃ = zero(T)
        for j in 1:f
            s̃ += φ[j] * s[i - j]
        end
        s̃ += rand(rng, d)
        s[i] = s̃
    end
    return view(s, f+1:N+f)
end


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

@test integrate(F^(c*(a+b*x))*sin(d+e*x)^3, x) == :(-6*F^(c*(a+b*x))*e^3*(9*e^4+b^4*c^4*log(F)^4+10*b^2*c^2*e^2*log(F)^2)^-1*cos(d+e*x)+-3*e*F^(c*(a+b*x))*(9*e^2+b^2*c^2*log(F)^2)^-1*sin(d+e*x)^2*cos(d+e*x)+b*c*F^(c*(a+b*x))*(9*e^2+b^2*c^2*log(F)^2)^-1*sin(d+e*x)^3*log(F)+6*b*c*F^(c*(a+b*x))*e^2*(9*e^4+b^4*c^4*log(F)^4+10*b^2*c^2*e^2*log(F)^2)^-1*log(F)*sin(d+e*x))
@test integrate(F^(c*(a+b*x))*sin(d+e*x)^2, x) == :(-2*e*F^(c*(a+b*x))*(4*e^2+b^2*c^2*log(F)^2)^-1*cos(d+e*x)*sin(d+e*x)+b*c*F^(c*(a+b*x))*(4*e^2+b^2*c^2*log(F)^2)^-1*sin(d+e*x)^2*log(F)+2*F^(c*(a+b*x))*b^-1*c^-1*e^2*(4*e^2+b^2*c^2*log(F)^2)^-1*log(F)^-1)
@test integrate(F^(c*(a+b*x))*sin(d+e*x), x) == :(-1*e*F^(c*(a+b*x))*(e^2+b^2*c^2*log(F)^2)^-1*cos(d+e*x)+b*c*F^(c*(a+b*x))*(e^2+b^2*c^2*log(F)^2)^-1*log(F)*sin(d+e*x))
@test integrate(sin(x)^4*exp(x), x) == :(24//85*exp(x)+1//17*sin(x)^4*exp(x)+12//85*sin(x)^2*exp(x)+-24//85*cos(x)*exp(x)*sin(x)+-4//17*sin(x)^3*cos(x)*exp(x))
@test integrate(F^(c*(a+b*x))*cos(d+e*x)^3, x) == :(6*F^(c*(a+b*x))*e^3*(9*e^4+b^4*c^4*log(F)^4+10*b^2*c^2*e^2*log(F)^2)^-1*sin(d+e*x)+3*e*F^(c*(a+b*x))*(9*e^2+b^2*c^2*log(F)^2)^-1*cos(d+e*x)^2*sin(d+e*x)+b*c*F^(c*(a+b*x))*(9*e^2+b^2*c^2*log(F)^2)^-1*cos(d+e*x)^3*log(F)+6*b*c*F^(c*(a+b*x))*e^2*(9*e^4+b^4*c^4*log(F)^4+10*b^2*c^2*e^2*log(F)^2)^-1*cos(d+e*x)*log(F))
@test integrate(F^(c*(a+b*x))*cos(d+e*x)^2, x) == :(2*e*F^(c*(a+b*x))*(4*e^2+b^2*c^2*log(F)^2)^-1*cos(d+e*x)*sin(d+e*x)+b*c*F^(c*(a+b*x))*(4*e^2+b^2*c^2*log(F)^2)^-1*cos(d+e*x)^2*log(F)+2*F^(c*(a+b*x))*b^-1*c^-1*e^2*(4*e^2+b^2*c^2*log(F)^2)^-1*log(F)^-1)
@test integrate(F^(c*(a+b*x))*cos(d+e*x), x) == :(e*F^(c*(a+b*x))*(e^2+b^2*c^2*log(F)^2)^-1*sin(d+e*x)+b*c*F^(c*(a+b*x))*(e^2+b^2*c^2*log(F)^2)^-1*cos(d+e*x)*log(F))
@test integrate(cos(x)^4*exp(x), x) == :(24//85*exp(x)+1//17*cos(x)^4*exp(x)+12//85*cos(x)^2*exp(x)+4//17*cos(x)^3*exp(x)*sin(x)+24//85*cos(x)*exp(x)*sin(x))
@test integrate(F^(c*(a+b*x))*(f*x)^m*sin(d+e*x)^2, x) == :(F^(a*c)*(f*x)^m*(x*(2*im*e+-1*b*c*log(F)))^(-1m)*(8*im*e+-4*b*c*log(F))^-1*SpecialFunctions.gamma(1+m,x*(2*im*e+-1*b*c*log(F)))*exp(-2*im*d)+-1*F^(a*c)*(f*x)^m*(-1*x*(2*im*e+b*c*log(F)))^(-1m)*(8*im*e+4*b*c*log(F))^-1*SpecialFunctions.gamma(1+m,-1*x*(2*im*e+b*c*log(F)))*exp(2*im*d)+(1/2)*F^(a*c)*b^-1*c^-1*(f*x)^m*(-1*b*c*x*log(F))^(-1m)*log(F)^-1*SpecialFunctions.gamma(1+m,-1*b*c*x*log(F)))
@test integrate(F^(c*(a+b*x))*(f*x)^m*sin(d+e*x), x) == :(-1*F^(a*c)*(f*x)^m*(x*(im*e+-1*b*c*log(F)))^(-1m)*(2*e+2*im*b*c*log(F))^-1*SpecialFunctions.gamma(1+m,x*(im*e+-1*b*c*log(F)))*exp(-1*im*d)+-1*F^(a*c)*(f*x)^m*(-1*x*(im*e+b*c*log(F)))^(-1m)*(2*e+-2*im*b*c*log(F))^-1*SpecialFunctions.gamma(1+m,-1*x*(im*e+b*c*log(F)))*exp(im*d))
@test integrate(f*F^(c*(a+b*x))*(f*x)^(-2+m)*((-1+m+b*c*x*log(F))*sin(d+e*x)+e*x*cos(d+e*x)), x) == :(F^(a*c+b*c*x)*(f*x)^(-1+m)*sin(d+e*x))
@test integrate(f*F^(c*(a+b*x))*(f*x)^m*((1+m+b*c*x*log(F))*sin(d+e*x)+e*x*cos(d+e*x)), x) == :(f*x*F^(c*(a+b*x))*(f*x)^m*sin(d+e*x))
@test integrate(F^(c*(a+b*x))*x^-1*(f*x)^m*((m+b*c*x*log(F))*sin(d+e*x)+e*x*cos(d+e*x)), x) == :(F^(a*c+b*c*x)*(f*x)^m*sin(d+e*x))
@test integrate(F^(c*(a+b*x))*(e*cos(d+e*x)+b*c*log(F)*sin(d+e*x)), x) == :(F^(c*(a+b*x))*sin(d+e*x))
@test integrate(F^(c*(a+b*x))*x^-2*((-1+b*c*x*log(F))*sin(d+e*x)+e*x*cos(d+e*x)), x) == :(F^(a*c+b*c*x)*x^-1*sin(d+e*x))
@test integrate(F^(c*(a+b*x))*x^-3*((-2+b*c*x*log(F))*sin(d+e*x)+e*x*cos(d+e*x)), x) == :(F^(a*c+b*c*x)*x^-2*sin(d+e*x))
@test integrate(cos(c+d*x)*exp(a+b*x)*sin(c+d*x), x) == :(b*(2*b^2+8*d^2)^-1*exp(a+b*x)*sin(2c+2*d*x)+-1*d*(b^2+4*d^2)^-1*cos(2c+2*d*x)*exp(a+b*x))
@test integrate(sin(c+d*x)^2*cos(c+d*x)*exp(a+b*x), x) == :(b*(4*b^2+4*d^2)^-1*cos(c+d*x)*exp(a+b*x)+d*(4*b^2+4*d^2)^-1*exp(a+b*x)*sin(c+d*x)+-1*b*(4*b^2+36*d^2)^-1*cos(3c+3*d*x)*exp(a+b*x)+-3*d*(4*b^2+36*d^2)^-1*exp(a+b*x)*sin(3c+3*d*x))
@test integrate(sin(c+d*x)^3*cos(c+d*x)*exp(a+b*x), x) == :(b*(4*b^2+16*d^2)^-1*exp(a+b*x)*sin(2c+2*d*x)+d*(2*b^2+32*d^2)^-1*cos(4c+4*d*x)*exp(a+b*x)+-1*b*(8*b^2+128*d^2)^-1*exp(a+b*x)*sin(4c+4*d*x)+-1*d*(2*b^2+8*d^2)^-1*cos(2c+2*d*x)*exp(a+b*x))
@test integrate(cos(c+d*x)^2*exp(a+b*x)*sin(c+d*x), x) == :(b*(4*b^2+4*d^2)^-1*exp(a+b*x)*sin(c+d*x)+b*(4*b^2+36*d^2)^-1*exp(a+b*x)*sin(3c+3*d*x)+-1*d*(4*b^2+4*d^2)^-1*cos(c+d*x)*exp(a+b*x)+-3*d*(4*b^2+36*d^2)^-1*cos(3c+3*d*x)*exp(a+b*x))
@test integrate(cos(c+d*x)^2*sin(c+d*x)^2*exp(a+b*x), x) == :(1//8*b^-1*exp(a+b*x)+-1*b*(8*b^2+128*d^2)^-1*cos(4c+4*d*x)*exp(a+b*x)+-1*d*(2*b^2+32*d^2)^-1*exp(a+b*x)*sin(4c+4*d*x))
@test integrate(cos(c+d*x)^2*sin(c+d*x)^3*exp(a+b*x), x) == :(b*(8*b^2+8*d^2)^-1*exp(a+b*x)*sin(c+d*x)+b*(16*b^2+144*d^2)^-1*exp(a+b*x)*sin(3c+3*d*x)+-1*b*(16*b^2+400*d^2)^-1*exp(a+b*x)*sin(5c+5*d*x)+-1*d*(8*b^2+8*d^2)^-1*cos(c+d*x)*exp(a+b*x)+-3*d*(16*b^2+144*d^2)^-1*cos(3c+3*d*x)*exp(a+b*x)+5*d*(16*b^2+400*d^2)^-1*cos(5c+5*d*x)*exp(a+b*x))
@test integrate(cos(c+d*x)^3*exp(a+b*x)*sin(c+d*x), x) == :(b*(4*b^2+16*d^2)^-1*exp(a+b*x)*sin(2c+2*d*x)+b*(8*b^2+128*d^2)^-1*exp(a+b*x)*sin(4c+4*d*x)+-1*d*(2*b^2+8*d^2)^-1*cos(2c+2*d*x)*exp(a+b*x)+-1*d*(2*b^2+32*d^2)^-1*cos(4c+4*d*x)*exp(a+b*x))
@test integrate(cos(c+d*x)^3*sin(c+d*x)^2*exp(a+b*x), x) == :(b*(8*b^2+8*d^2)^-1*cos(c+d*x)*exp(a+b*x)+d*(8*b^2+8*d^2)^-1*exp(a+b*x)*sin(c+d*x)+-1*b*(16*b^2+144*d^2)^-1*cos(3c+3*d*x)*exp(a+b*x)+-1*b*(16*b^2+400*d^2)^-1*cos(5c+5*d*x)*exp(a+b*x)+-5*d*(16*b^2+400*d^2)^-1*exp(a+b*x)*sin(5c+5*d*x)+-3*d*(16*b^2+144*d^2)^-1*exp(a+b*x)*sin(3c+3*d*x))
@test integrate(cos(c+d*x)^3*sin(c+d*x)^3*exp(a+b*x), x) == :(-1*b*(32*b^2+1152*d^2)^-1*exp(a+b*x)*sin(6c+6*d*x)+-3*d*(16*b^2+64*d^2)^-1*cos(2c+2*d*x)*exp(a+b*x)+3*b*(32*b^2+128*d^2)^-1*exp(a+b*x)*sin(2c+2*d*x)+3*d*(16*b^2+576*d^2)^-1*cos(6c+6*d*x)*exp(a+b*x))
@test integrate(x*exp(x)*sin(x), x) == :((1/2)*cos(x)*exp(x)+(1/2)*x*exp(x)*sin(x)+-1//2*x*cos(x)*exp(x))
@test integrate(x^2*exp(x)*sin(x), x) == :(-1//2*cos(x)*exp(x)+-1//2*exp(x)*sin(x)+x*cos(x)*exp(x)+(1/2)*x^2*exp(x)*sin(x)+-1//2*x^2*cos(x)*exp(x))
@test integrate(x*cos(x)*exp(x), x) == :(-1//2*exp(x)*sin(x)+(1/2)*x*cos(x)*exp(x)+(1/2)*x*exp(x)*sin(x))
@test integrate(x^2*cos(x)*exp(x), x) == :((1/2)*exp(x)*sin(x)+-1//2*cos(x)*exp(x)+(1/2)*x^2*cos(x)*exp(x)+(1/2)*x^2*exp(x)*sin(x)+-1*x*exp(x)*sin(x))
@test integrate((-5*cos(4x)+2*sin(4x))*exp(3x), x) == :(-23//25*cos(4x)*exp(3x)+-14//25*exp(3x)*sin(4x))
@test integrate(exp(x)*sin(x)+exp(-1x)*sin(x), x) == :((1/2)*exp(x)*sin(x)+-1//2*cos(x)*exp(x)+-1//2*cos(x)*exp(-1x)+-1//2*exp(-1x)*sin(x))
@test integrate(exp(x^2)*sin(b*x), x) == :(-1//4*im*pi^(1/2)*erfi(x+(1/2)*im*b)*exp(1//4*b^2)+-1//4*im*pi^(1/2)*erfi(-1x+(1/2)*im*b)*exp(1//4*b^2))
@test integrate(cos(b*x)*exp(x^2), x) == :(-1//4*pi^(1/2)*erfi(-1x+(1/2)*im*b)*exp(1//4*b^2)+1//4*pi^(1/2)*erfi(x+(1/2)*im*b)*exp(1//4*b^2))
@test integrate(exp(x^2)*sin(a+b*x), x) == :(-1//4*im*pi^(1/2)*erfi(x+(1/2)*im*b)*exp(1//4*b^2+im*a)+-1//4*im*pi^(1/2)*erfi(-1x+(1/2)*im*b)*exp(1//4*b^2+-1*im*a))
@test integrate(cos(a+b*x)*exp(x^2), x) == :(-1//4*pi^(1/2)*erfi(-1x+(1/2)*im*b)*exp(1//4*b^2+-1*im*a)+1//4*pi^(1/2)*erfi(x+(1/2)*im*b)*exp(1//4*b^2+im*a))
@test integrate(x*cos(2*x^2)*exp(2*x^2), x) == :(1//8*cos(2*x^2)*exp(2*x^2)+1//8*exp(2*x^2)*sin(2*x^2))
@test integrate(exp(x)*sin(exp(x)), x) == :(-1*cos(exp(x)))
@test integrate(csc(exp(x))*exp(x)*sec(exp(x)), x) == :(log(tan(exp(x))))
@test integrate(cos(exp(x))*exp(x), x) == :(sin(exp(x)))
@test integrate(cos(exp(2x))*exp(2x), x) == :((1/2)*sin(exp(2x)))
@test integrate(cos(exp(-2x))*exp(-2x), x) == :(-1//2*sin(exp(-2x)))
@test integrate(cos(exp(x))*exp(2x), x) == :(exp(x)*sin(exp(x))+cos(exp(x)))
@test integrate(cos(exp(-1+3x))*exp(-1+3x)*sin(1+exp(-1+3x)), x) == :(-1//12*cos(1+2*exp(-1+3x))+1//6*exp(-1+3x)*sin(1))
@test integrate(exp(x)*tan(exp(x)), x) == :(-1*log(cos(exp(x))))
@test integrate(exp(x)*sec(exp(x)), x) == :(arctanh(sin(exp(x))))
@test integrate(exp(x)*sec(exp(x))*tan(exp(x)), x) == :(sec(exp(x)))
@test integrate(csc(exp(x))^2*exp(x), x) == :(-1*cot(exp(x)))
@test integrate(exp(x)*sin(a+b*x), x) == :((1+b^2)^-1*exp(x)*sin(a+b*x)+-1*b*(1+b^2)^-1*cos(a+b*x)*exp(x))
@test integrate(exp(x)*sin(a+c*x^2), x) == :(1//4*(-1)^3//4*pi^(1/2)*c^-1//2*erf((1/2)*(-1)^1//4*c^-1//2*(1+2*im*c*x))*exp(1//4*im*(c^-1+4a))+1//4*(-1)^3//4*pi^(1/2)*c^-1//2*erfi((1/2)*(-1)^1//4*c^-1//2*(1+-2*im*c*x))*exp(-1//4*im*(c^-1+4a)))
@test integrate(exp(x)*sin(a+b*x+c*x^2), x) == :(1//4*(-1)^3//4*pi^(1/2)*c^-1//2*erf((1/2)*(-1)^1//4*c^-1//2*(1+im*b+2*im*c*x))*exp(1//4*im*(4a+c^-1*(1+im*b)^2))+1//4*(-1)^3//4*pi^(1/2)*c^-1//2*erfi((1/2)*(-1)^1//4*c^-1//2*(1+-1*im*b+-2*im*c*x))*exp(-1*im*a+1//4*im*c^-1*(im+b)^2))
@test integrate(exp(x^2)*sin(a+b*x), x) == :(-1//4*im*pi^(1/2)*erfi(x+(1/2)*im*b)*exp(1//4*b^2+im*a)+-1//4*im*pi^(1/2)*erfi(-1x+(1/2)*im*b)*exp(1//4*b^2+-1*im*a))
@test integrate(exp(x^2)*sin(a+c*x^2), x) == :(-1//4*im*pi^(1/2)*(1+im*c)^-1//2*erfi(x*(1+im*c)^(1/2))*exp(im*a)+1//4*im*pi^(1/2)*(1+-1*im*c)^-1//2*erfi(x*(1+-1*im*c)^(1/2))*exp(-1*im*a))
@test integrate(exp(x^2)*sin(a+b*x+c*x^2), x) == :(-1//4*im*pi^(1/2)*(1+im*c)^-1//2*erfi((1/2)*(1+im*c)^-1//2*(im*b+x*(2+2*im*c)))*exp(im*a+b^2*(4+4*im*c)^-1)+-1//4*im*pi^(1/2)*(1+-1*im*c)^-1//2*erfi((1/2)*(1+-1*im*c)^-1//2*(im*b+-1*x*(2+-2*im*c)))*exp(-1*im*(a+-1*b^2*(4im+4c)^-1)))
@test integrate(f^(a+b*x)*sin(d+f*x^2), x) == :(-1//4*(-1)^3//4*pi^(1/2)*f^(-1//2+a)*erfi((1/2)*(-1)^1//4*f^-1//2*(-1*b*log(f)+2*im*f*x))*exp(-1//4*im*(4d+b^2*f^-1*log(f)^2))+1//4*(-1)^3//4*pi^(1/2)*f^(-1//2+a)*erf((1/2)*(-1)^1//4*f^-1//2*(b*log(f)+2*im*f*x))*exp(1//4*im*(4d+b^2*f^-1*log(f)^2)))
@test integrate(f^(a+b*x)*sin(d+f*x^2)^2, x) == :((1/2)*b^-1*f^(a+b*x)*log(f)^-1+pi^(1/2)*f^(-1//2+a)*(1//16+1//16*im)*erf(f^-1//2*(1//4+1//4*im)*(b*log(f)+4*im*f*x))*exp(2*im*d+1//8*im*b^2*f^-1*log(f)^2)+pi^(1/2)*f^(-1//2+a)*(1//16+1//16*im)*erfi(f^-1//2*(1//4+1//4*im)*(-1*b*log(f)+4*im*f*x))*exp(-1//8*im*(16d+b^2*f^-1*log(f)^2)))
@test integrate(f^(a+b*x)*sin(d+f*x^2)^3, x) == :(-3//16*(-1)^3//4*pi^(1/2)*f^(-1//2+a)*erfi((1/2)*(-1)^1//4*f^-1//2*(-1*b*log(f)+2*im*f*x))*exp(-1//4*im*(4d+b^2*f^-1*log(f)^2))+3//16*(-1)^3//4*pi^(1/2)*f^(-1//2+a)*erf((1/2)*(-1)^1//4*f^-1//2*(b*log(f)+2*im*f*x))*exp(1//4*im*(4d+b^2*f^-1*log(f)^2))+-1//6*6^(1/2)*pi^(1/2)*f^(-1//2+a)*(1//16+-1//16*im)*erfi(1//6*6^(1/2)*f^-1//2*(1/2+(1/2)*im)*(-1*b*log(f)+6*im*f*x))*exp(-1//12*im*(36d+b^2*f^-1*log(f)^2))+1//6*6^(1/2)*pi^(1/2)*f^(-1//2+a)*(1//16+-1//16*im)*erf(1//6*6^(1/2)*f^-1//2*(1/2+(1/2)*im)*(b*log(f)+6*im*f*x))*exp(3*im*d+1//12*im*b^2*f^-1*log(f)^2))
@test integrate(f^(a+b*x)*sin(d+e*x+f*x^2), x) == :(-1//4*(-1)^3//4*pi^(1/2)*f^(-1//2+a)*erfi((1/2)*(-1)^1//4*f^-1//2*(im*e+-1*b*log(f)+2*im*f*x))*exp(-1*im*d+1//4*im*f^-1*(e+im*b*log(f))^2)+1//4*(-1)^3//4*pi^(1/2)*f^(-1//2+a)*erf((1/2)*(-1)^1//4*f^-1//2*(im*e+b*log(f)+2*im*f*x))*exp(1//4*im*(4d+f^-1*(im*e+b*log(f))^2)))
@test integrate(f^(a+b*x)*sin(d+e*x+f*x^2)^2, x) == :((1/2)*b^-1*f^(a+b*x)*log(f)^-1+pi^(1/2)*f^(-1//2+a)*(1//16+1//16*im)*erf(f^-1//2*(1//4+1//4*im)*(b*log(f)+2*im*e+4*im*f*x))*exp(2*im*d+1//8*im*f^-1*(b*log(f)+2*im*e)^2)+pi^(1/2)*f^(-1//2+a)*(1//16+1//16*im)*erfi(f^-1//2*(1//4+1//4*im)*(-1*b*log(f)+2*im*e+4*im*f*x))*exp(-2*im*d+1//8*im*f^-1*(2*e+im*b*log(f))^2))
@test integrate(f^(a+b*x)*sin(d+e*x+f*x^2)^3, x) == :(-3//16*(-1)^3//4*pi^(1/2)*f^(-1//2+a)*erfi((1/2)*(-1)^1//4*f^-1//2*(im*e+-1*b*log(f)+2*im*f*x))*exp(-1*im*d+1//4*im*f^-1*(e+im*b*log(f))^2)+3//16*(-1)^3//4*pi^(1/2)*f^(-1//2+a)*erf((1/2)*(-1)^1//4*f^-1//2*(im*e+b*log(f)+2*im*f*x))*exp(1//4*im*(4d+f^-1*(im*e+b*log(f))^2))+-1//6*6^(1/2)*pi^(1/2)*f^(-1//2+a)*(1//16+-1//16*im)*erfi(1//6*6^(1/2)*f^-1//2*(1/2+(1/2)*im)*(-1*b*log(f)+3*im*e+6*im*f*x))*exp(-3*im*d+1//12*im*f^-1*(3*e+im*b*log(f))^2)+1//6*6^(1/2)*pi^(1/2)*f^(-1//2+a)*(1//16+-1//16*im)*erf(1//6*6^(1/2)*f^-1//2*(1/2+(1/2)*im)*(b*log(f)+3*im*e+6*im*f*x))*exp(3*im*d+1//12*im*f^-1*(b*log(f)+3*im*e)^2))
@test integrate(f^(a+c*x^2)*sin(d+f*x^2), x) == :(-1//4*im*pi^(1/2)*f^a*(im*f+c*log(f))^-1//2*erfi(x*(im*f+c*log(f))^(1/2))*exp(im*d)+1//4*im*pi^(1/2)*f^a*(im*f+-1*c*log(f))^-1//2*erf(x*(im*f+-1*c*log(f))^(1/2))*exp(-1*im*d))
@test integrate(f^(a+c*x^2)*sin(d+f*x^2)^2, x) == :(-1//8*pi^(1/2)*f^a*(c*log(f)+2*im*f)^-1//2*erfi(x*(c*log(f)+2*im*f)^(1/2))*exp(2*im*d)+-1//8*pi^(1/2)*f^a*(-1*c*log(f)+2*im*f)^-1//2*erf(x*(-1*c*log(f)+2*im*f)^(1/2))*exp(-2*im*d)+1//4*pi^(1/2)*c^-1//2*f^a*log(f)^-1//2*erfi(x*c^(1/2)*log(f)^(1/2)))
@test integrate(f^(a+c*x^2)*sin(d+f*x^2)^3, x) == :(-3//16*im*pi^(1/2)*f^a*(im*f+c*log(f))^-1//2*erfi(x*(im*f+c*log(f))^(1/2))*exp(im*d)+-1//16*im*pi^(1/2)*f^a*(-1*c*log(f)+3*im*f)^-1//2*erf(x*(-1*c*log(f)+3*im*f)^(1/2))*exp(-3*im*d)+1//16*im*pi^(1/2)*f^a*(c*log(f)+3*im*f)^-1//2*erfi(x*(c*log(f)+3*im*f)^(1/2))*exp(3*im*d)+3//16*im*pi^(1/2)*f^a*(im*f+-1*c*log(f))^-1//2*erf(x*(im*f+-1*c*log(f))^(1/2))*exp(-1*im*d))
@test integrate(f^(a+c*x^2)*sin(d+e*x+f*x^2), x) == :(-1//4*im*pi^(1/2)*f^a*(im*f+c*log(f))^-1//2*erfi((1/2)*(im*f+c*log(f))^-1//2*(im*e+2*x*(im*f+c*log(f))))*exp(im*d+e^2*(4*im*f+4*c*log(f))^-1)+1//4*im*pi^(1/2)*f^a*(im*f+-1*c*log(f))^-1//2*erf((1/2)*(im*f+-1*c*log(f))^-1//2*(im*e+2*x*(im*f+-1*c*log(f))))*exp(-1*im*d+-1*e^2*(-4*c*log(f)+4*im*f)^-1))
@test integrate(f^(a+c*x^2)*sin(d+e*x+f*x^2)^2, x) == :(-1//8*pi^(1/2)*f^a*(c*log(f)+2*im*f)^-1//2*erfi((c*log(f)+2*im*f)^-1//2*(im*e+x*(c*log(f)+2*im*f)))*exp(e^2*(c*log(f)+2*im*f)^-1+2*im*d)+-1//8*pi^(1/2)*f^a*(-1*c*log(f)+2*im*f)^-1//2*erf((-1*c*log(f)+2*im*f)^-1//2*(im*e+x*(-1*c*log(f)+2*im*f)))*exp(-1*e^2*(-1*c*log(f)+2*im*f)^-1+-2*im*d)+1//4*pi^(1/2)*c^-1//2*f^a*log(f)^-1//2*erfi(x*c^(1/2)*log(f)^(1/2)))
@test integrate(f^(a+c*x^2)*sin(d+e*x+f*x^2)^3, x) == :(-3//16*im*pi^(1/2)*f^a*(im*f+c*log(f))^-1//2*erfi((1/2)*(im*f+c*log(f))^-1//2*(im*e+2*x*(im*f+c*log(f))))*exp(im*d+e^2*(4*im*f+4*c*log(f))^-1)+-1//16*im*pi^(1/2)*f^a*(-1*c*log(f)+3*im*f)^-1//2*erf((1/2)*(-1*c*log(f)+3*im*f)^-1//2*(2*x*(-1*c*log(f)+3*im*f)+3*im*e))*exp(-9*e^2*(-4*c*log(f)+12*im*f)^-1+-3*im*d)+1//16*im*pi^(1/2)*f^a*(c*log(f)+3*im*f)^-1//2*erfi((1/2)*(c*log(f)+3*im*f)^-1//2*(2*x*(c*log(f)+3*im*f)+3*im*e))*exp(3*im*d+9*e^2*(4*c*log(f)+12*im*f)^-1)+3//16*im*pi^(1/2)*f^a*(im*f+-1*c*log(f))^-1//2*erf((1/2)*(im*f+-1*c*log(f))^-1//2*(im*e+2*x*(im*f+-1*c*log(f))))*exp(-1*im*d+-1*e^2*(-4*c*log(f)+4*im*f)^-1))
@test integrate(f^(a+b*x+c*x^2)*sin(d+f*x^2), x) == :(-1//4*im*pi^(1/2)*f^a*(im*f+c*log(f))^-1//2*erfi((1/2)*(im*f+c*log(f))^-1//2*(b*log(f)+2*x*(im*f+c*log(f))))*exp(im*d+-1*b^2*(4*im*f+4*c*log(f))^-1*log(f)^2)+-1//4*im*pi^(1/2)*f^a*(im*f+-1*c*log(f))^-1//2*erf((1/2)*(im*f+-1*c*log(f))^-1//2*(b*log(f)+-2*x*(im*f+-1*c*log(f))))*exp(-1*im*d+b^2*(-4*c*log(f)+4*im*f)^-1*log(f)^2))
@test integrate(f^(a+b*x+c*x^2)*sin(d+f*x^2)^2, x) == :(-1//8*pi^(1/2)*f^a*(c*log(f)+2*im*f)^-1//2*erfi((1/2)*(c*log(f)+2*im*f)^-1//2*(b*log(f)+2*x*(c*log(f)+2*im*f)))*exp(2*im*d+-1*b^2*(4*c*log(f)+8*im*f)^-1*log(f)^2)+1//4*pi^(1/2)*c^-1//2*f^(a+-1//4*b^2*c^-1)*log(f)^-1//2*erfi((1/2)*c^-1//2*log(f)^(1/2)*(b+2*c*x))+1//8*pi^(1/2)*f^a*(-1*c*log(f)+2*im*f)^-1//2*erf((1/2)*(-1*c*log(f)+2*im*f)^-1//2*(b*log(f)+-2*x*(-1*c*log(f)+2*im*f)))*exp(-2*im*d+b^2*(-4*c*log(f)+8*im*f)^-1*log(f)^2))
@test integrate(f^(a+b*x+c*x^2)*sin(d+f*x^2)^3, x) == :(-3//16*im*pi^(1/2)*f^a*(im*f+c*log(f))^-1//2*erfi((1/2)*(im*f+c*log(f))^-1//2*(b*log(f)+2*x*(im*f+c*log(f))))*exp(im*d+-1*b^2*(4*im*f+4*c*log(f))^-1*log(f)^2)+-3//16*im*pi^(1/2)*f^a*(im*f+-1*c*log(f))^-1//2*erf((1/2)*(im*f+-1*c*log(f))^-1//2*(b*log(f)+-2*x*(im*f+-1*c*log(f))))*exp(-1*im*d+b^2*(-4*c*log(f)+4*im*f)^-1*log(f)^2)+1//16*im*pi^(1/2)*f^a*(c*log(f)+3*im*f)^-1//2*erfi((1/2)*(c*log(f)+3*im*f)^-1//2*(b*log(f)+2*x*(c*log(f)+3*im*f)))*exp(3*im*d+-1*b^2*(4*c*log(f)+12*im*f)^-1*log(f)^2)+1//16*im*pi^(1/2)*f^a*(-1*c*log(f)+3*im*f)^-1//2*erf((1/2)*(-1*c*log(f)+3*im*f)^-1//2*(b*log(f)+-2*x*(-1*c*log(f)+3*im*f)))*exp(-3*im*d+b^2*(-4*c*log(f)+12*im*f)^-1*log(f)^2))
@test integrate(f^(a+b*x+c*x^2)*sin(d+e*x+f*x^2), x) == :(-1//4*im*pi^(1/2)*f^a*(im*f+c*log(f))^-1//2*erfi((1/2)*(im*f+c*log(f))^-1//2*(im*e+b*log(f)+2*x*(im*f+c*log(f))))*exp(im*d+(e+-1*im*b*log(f))^2*(4*im*f+4*c*log(f))^-1)+1//4*im*pi^(1/2)*f^a*(im*f+-1*c*log(f))^-1//2*erf((1/2)*(im*f+-1*c*log(f))^-1//2*(im*e+-1*b*log(f)+2*x*(im*f+-1*c*log(f))))*exp(-1*im*d+-1*(e+im*b*log(f))^2*(-4*c*log(f)+4*im*f)^-1))
@test integrate(f^(a+b*x+c*x^2)*sin(d+e*x+f*x^2)^2, x) == :(-1//8*pi^(1/2)*f^a*(c*log(f)+2*im*f)^-1//2*erfi((1/2)*(c*log(f)+2*im*f)^-1//2*(b*log(f)+2*im*e+2*x*(c*log(f)+2*im*f)))*exp((2*e+-1*im*b*log(f))^2*(4*c*log(f)+8*im*f)^-1+2*im*d)+-1//8*pi^(1/2)*f^a*(-1*c*log(f)+2*im*f)^-1//2*erf((1/2)*(-1*c*log(f)+2*im*f)^-1//2*(-1*b*log(f)+2*im*e+2*x*(-1*c*log(f)+2*im*f)))*exp(-1*(2*e+im*b*log(f))^2*(-4*c*log(f)+8*im*f)^-1+-2*im*d)+1//4*pi^(1/2)*c^-1//2*f^(a+-1//4*b^2*c^-1)*log(f)^-1//2*erfi((1/2)*c^-1//2*log(f)^(1/2)*(b+2*c*x)))
@test integrate(f^(a+b*x+c*x^2)*sin(d+e*x+f*x^2)^3, x) == :(-3//16*im*pi^(1/2)*f^a*(im*f+c*log(f))^-1//2*erfi((1/2)*(im*f+c*log(f))^-1//2*(im*e+b*log(f)+2*x*(im*f+c*log(f))))*exp(im*d+(e+-1*im*b*log(f))^2*(4*im*f+4*c*log(f))^-1)+-1//16*im*pi^(1/2)*f^a*(-1*c*log(f)+3*im*f)^-1//2*erf((1/2)*(-1*c*log(f)+3*im*f)^-1//2*(-1*b*log(f)+2*x*(-1*c*log(f)+3*im*f)+3*im*e))*exp(-1*(3*e+im*b*log(f))^2*(-4*c*log(f)+12*im*f)^-1+-3*im*d)+1//16*im*pi^(1/2)*f^a*(c*log(f)+3*im*f)^-1//2*erfi((1/2)*(c*log(f)+3*im*f)^-1//2*(b*log(f)+2*x*(c*log(f)+3*im*f)+3*im*e))*exp(-1*(b*log(f)+3*im*e)^2*(4*c*log(f)+12*im*f)^-1+3*im*d)+3//16*im*pi^(1/2)*f^a*(im*f+-1*c*log(f))^-1//2*erf((1/2)*(im*f+-1*c*log(f))^-1//2*(im*e+-1*b*log(f)+2*x*(im*f+-1*c*log(f))))*exp(-1*im*d+-1*(e+im*b*log(f))^2*(-4*c*log(f)+4*im*f)^-1))
@test integrate(cos(a+b*x)*exp(x), x) == :((1+b^2)^-1*cos(a+b*x)*exp(x)+b*(1+b^2)^-1*exp(x)*sin(a+b*x))
@test integrate(cos(a+c*x^2)*exp(x), x) == :(-1//4*(-1)^1//4*pi^(1/2)*c^-1//2*erf((1/2)*(-1)^1//4*c^-1//2*(1+2*im*c*x))*exp(1//4*im*(c^-1+4a))+1//4*(-1)^1//4*pi^(1/2)*c^-1//2*erfi((1/2)*(-1)^1//4*c^-1//2*(1+-2*im*c*x))*exp(-1//4*im*(c^-1+4a)))
@test integrate(cos(a+b*x+c*x^2)*exp(x), x) == :(-1//4*(-1)^1//4*pi^(1/2)*c^-1//2*erf((1/2)*(-1)^1//4*c^-1//2*(1+im*b+2*im*c*x))*exp(1//4*im*(4a+c^-1*(1+im*b)^2))+1//4*(-1)^1//4*pi^(1/2)*c^-1//2*erfi((1/2)*(-1)^1//4*c^-1//2*(1+-1*im*b+-2*im*c*x))*exp(-1*im*a+1//4*im*c^-1*(im+b)^2))
@test integrate(cos(a+b*x)*exp(x^2), x) == :(-1//4*pi^(1/2)*erfi(-1x+(1/2)*im*b)*exp(1//4*b^2+-1*im*a)+1//4*pi^(1/2)*erfi(x+(1/2)*im*b)*exp(1//4*b^2+im*a))
@test integrate(cos(a+c*x^2)*exp(x^2), x) == :(1//4*pi^(1/2)*(1+im*c)^-1//2*erfi(x*(1+im*c)^(1/2))*exp(im*a)+1//4*pi^(1/2)*(1+-1*im*c)^-1//2*erfi(x*(1+-1*im*c)^(1/2))*exp(-1*im*a))
@test integrate(cos(a+b*x+c*x^2)*exp(x^2), x) == :(-1//4*pi^(1/2)*(1+-1*im*c)^-1//2*erfi((1/2)*(1+-1*im*c)^-1//2*(im*b+-1*x*(2+-2*im*c)))*exp(-1*im*(a+-1*b^2*(4im+4c)^-1))+1//4*pi^(1/2)*(1+im*c)^-1//2*erfi((1/2)*(1+im*c)^-1//2*(im*b+x*(2+2*im*c)))*exp(im*a+b^2*(4+4*im*c)^-1))
@test integrate(f^(a+b*x)*cos(d+f*x^2), x) == :(-1//4*(-1)^1//4*pi^(1/2)*f^(-1//2+a)*erf((1/2)*(-1)^1//4*f^-1//2*(b*log(f)+2*im*f*x))*exp(1//4*im*(4d+b^2*f^-1*log(f)^2))+-1//4*(-1)^1//4*pi^(1/2)*f^(-1//2+a)*erfi((1/2)*(-1)^1//4*f^-1//2*(-1*b*log(f)+2*im*f*x))*exp(-1//4*im*(4d+b^2*f^-1*log(f)^2)))
@test integrate(f^(a+b*x)*cos(d+f*x^2)^2, x) == :((1/2)*b^-1*f^(a+b*x)*log(f)^-1+pi^(1/2)*f^(-1//2+a)*(-1//16+-1//16*im)*erf(f^-1//2*(1//4+1//4*im)*(b*log(f)+4*im*f*x))*exp(2*im*d+1//8*im*b^2*f^-1*log(f)^2)+-1*pi^(1/2)*f^(-1//2+a)*(1//16+1//16*im)*erfi(f^-1//2*(1//4+1//4*im)*(-1*b*log(f)+4*im*f*x))*exp(-1//8*im*(16d+b^2*f^-1*log(f)^2)))
@test integrate(f^(a+b*x)*cos(d+f*x^2)^3, x) == :(-3//16*(-1)^1//4*pi^(1/2)*f^(-1//2+a)*erf((1/2)*(-1)^1//4*f^-1//2*(b*log(f)+2*im*f*x))*exp(1//4*im*(4d+b^2*f^-1*log(f)^2))+-3//16*(-1)^1//4*pi^(1/2)*f^(-1//2+a)*erfi((1/2)*(-1)^1//4*f^-1//2*(-1*b*log(f)+2*im*f*x))*exp(-1//4*im*(4d+b^2*f^-1*log(f)^2))+-1//6*6^(1/2)*pi^(1/2)*f^(-1//2+a)*(1//16+1//16*im)*erf(1//6*6^(1/2)*f^-1//2*(1/2+(1/2)*im)*(b*log(f)+6*im*f*x))*exp(3*im*d+1//12*im*b^2*f^-1*log(f)^2)+-1//6*6^(1/2)*pi^(1/2)*f^(-1//2+a)*(1//16+1//16*im)*erfi(1//6*6^(1/2)*f^-1//2*(1/2+(1/2)*im)*(-1*b*log(f)+6*im*f*x))*exp(-1//12*im*(36d+b^2*f^-1*log(f)^2)))
@test integrate(f^(a+b*x)*cos(d+e*x+f*x^2), x) == :(-1//4*(-1)^1//4*pi^(1/2)*f^(-1//2+a)*erf((1/2)*(-1)^1//4*f^-1//2*(im*e+b*log(f)+2*im*f*x))*exp(1//4*im*(4d+f^-1*(im*e+b*log(f))^2))+-1//4*(-1)^1//4*pi^(1/2)*f^(-1//2+a)*erfi((1/2)*(-1)^1//4*f^-1//2*(im*e+-1*b*log(f)+2*im*f*x))*exp(-1*im*d+1//4*im*f^-1*(e+im*b*log(f))^2))
@test integrate(f^(a+b*x)*cos(d+e*x+f*x^2)^2, x) == :((1/2)*b^-1*f^(a+b*x)*log(f)^-1+pi^(1/2)*f^(-1//2+a)*(-1//16+-1//16*im)*erf(f^-1//2*(1//4+1//4*im)*(b*log(f)+2*im*e+4*im*f*x))*exp(2*im*d+1//8*im*f^-1*(b*log(f)+2*im*e)^2)+-1*pi^(1/2)*f^(-1//2+a)*(1//16+1//16*im)*erfi(f^-1//2*(1//4+1//4*im)*(-1*b*log(f)+2*im*e+4*im*f*x))*exp(-2*im*d+1//8*im*f^-1*(2*e+im*b*log(f))^2))
@test integrate(f^(a+b*x)*cos(d+e*x+f*x^2)^3, x) == :(-3//16*(-1)^1//4*pi^(1/2)*f^(-1//2+a)*erf((1/2)*(-1)^1//4*f^-1//2*(im*e+b*log(f)+2*im*f*x))*exp(1//4*im*(4d+f^-1*(im*e+b*log(f))^2))+-3//16*(-1)^1//4*pi^(1/2)*f^(-1//2+a)*erfi((1/2)*(-1)^1//4*f^-1//2*(im*e+-1*b*log(f)+2*im*f*x))*exp(-1*im*d+1//4*im*f^-1*(e+im*b*log(f))^2)+-1//6*6^(1/2)*pi^(1/2)*f^(-1//2+a)*(1//16+1//16*im)*erf(1//6*6^(1/2)*f^-1//2*(1/2+(1/2)*im)*(b*log(f)+3*im*e+6*im*f*x))*exp(3*im*d+1//12*im*f^-1*(b*log(f)+3*im*e)^2)+-1//6*6^(1/2)*pi^(1/2)*f^(-1//2+a)*(1//16+1//16*im)*erfi(1//6*6^(1/2)*f^-1//2*(1/2+(1/2)*im)*(-1*b*log(f)+3*im*e+6*im*f*x))*exp(-3*im*d+1//12*im*f^-1*(3*e+im*b*log(f))^2))
@test integrate(f^(a+c*x^2)*cos(d+f*x^2), x) == :(1//4*pi^(1/2)*f^a*(im*f+c*log(f))^-1//2*erfi(x*(im*f+c*log(f))^(1/2))*exp(im*d)+1//4*pi^(1/2)*f^a*(im*f+-1*c*log(f))^-1//2*erf(x*(im*f+-1*c*log(f))^(1/2))*exp(-1*im*d))
@test integrate(f^(a+c*x^2)*cos(d+f*x^2)^2, x) == :(1//4*pi^(1/2)*c^-1//2*f^a*log(f)^-1//2*erfi(x*c^(1/2)*log(f)^(1/2))+1//8*pi^(1/2)*f^a*(c*log(f)+2*im*f)^-1//2*erfi(x*(c*log(f)+2*im*f)^(1/2))*exp(2*im*d)+1//8*pi^(1/2)*f^a*(-1*c*log(f)+2*im*f)^-1//2*erf(x*(-1*c*log(f)+2*im*f)^(1/2))*exp(-2*im*d))
@test integrate(f^(a+c*x^2)*cos(d+f*x^2)^3, x) == :(1//16*pi^(1/2)*f^a*(c*log(f)+3*im*f)^-1//2*erfi(x*(c*log(f)+3*im*f)^(1/2))*exp(3*im*d)+1//16*pi^(1/2)*f^a*(-1*c*log(f)+3*im*f)^-1//2*erf(x*(-1*c*log(f)+3*im*f)^(1/2))*exp(-3*im*d)+3//16*pi^(1/2)*f^a*(im*f+c*log(f))^-1//2*erfi(x*(im*f+c*log(f))^(1/2))*exp(im*d)+3//16*pi^(1/2)*f^a*(im*f+-1*c*log(f))^-1//2*erf(x*(im*f+-1*c*log(f))^(1/2))*exp(-1*im*d))
@test integrate(f^(a+c*x^2)*cos(d+e*x+f*x^2), x) == :(1//4*pi^(1/2)*f^a*(im*f+c*log(f))^-1//2*erfi((1/2)*(im*f+c*log(f))^-1//2*(im*e+2*x*(im*f+c*log(f))))*exp(im*d+e^2*(4*im*f+4*c*log(f))^-1)+1//4*pi^(1/2)*f^a*(im*f+-1*c*log(f))^-1//2*erf((1/2)*(im*f+-1*c*log(f))^-1//2*(im*e+2*x*(im*f+-1*c*log(f))))*exp(-1*im*d+-1*e^2*(-4*c*log(f)+4*im*f)^-1))
@test integrate(f^(a+c*x^2)*cos(d+e*x+f*x^2)^2, x) == :(1//4*pi^(1/2)*c^-1//2*f^a*log(f)^-1//2*erfi(x*c^(1/2)*log(f)^(1/2))+1//8*pi^(1/2)*f^a*(c*log(f)+2*im*f)^-1//2*erfi((c*log(f)+2*im*f)^-1//2*(im*e+x*(c*log(f)+2*im*f)))*exp(e^2*(c*log(f)+2*im*f)^-1+2*im*d)+1//8*pi^(1/2)*f^a*(-1*c*log(f)+2*im*f)^-1//2*erf((-1*c*log(f)+2*im*f)^-1//2*(im*e+x*(-1*c*log(f)+2*im*f)))*exp(-1*e^2*(-1*c*log(f)+2*im*f)^-1+-2*im*d))
@test integrate(f^(a+c*x^2)*cos(d+e*x+f*x^2)^3, x) == :(1//16*pi^(1/2)*f^a*(c*log(f)+3*im*f)^-1//2*erfi((1/2)*(c*log(f)+3*im*f)^-1//2*(2*x*(c*log(f)+3*im*f)+3*im*e))*exp(3*im*d+9*e^2*(4*c*log(f)+12*im*f)^-1)+1//16*pi^(1/2)*f^a*(-1*c*log(f)+3*im*f)^-1//2*erf((1/2)*(-1*c*log(f)+3*im*f)^-1//2*(2*x*(-1*c*log(f)+3*im*f)+3*im*e))*exp(-9*e^2*(-4*c*log(f)+12*im*f)^-1+-3*im*d)+3//16*pi^(1/2)*f^a*(im*f+c*log(f))^-1//2*erfi((1/2)*(im*f+c*log(f))^-1//2*(im*e+2*x*(im*f+c*log(f))))*exp(im*d+e^2*(4*im*f+4*c*log(f))^-1)+3//16*pi^(1/2)*f^a*(im*f+-1*c*log(f))^-1//2*erf((1/2)*(im*f+-1*c*log(f))^-1//2*(im*e+2*x*(im*f+-1*c*log(f))))*exp(-1*im*d+-1*e^2*(-4*c*log(f)+4*im*f)^-1))
@test integrate(f^(a+b*x+c*x^2)*cos(d+f*x^2), x) == :(-1//4*pi^(1/2)*f^a*(im*f+-1*c*log(f))^-1//2*erf((1/2)*(im*f+-1*c*log(f))^-1//2*(b*log(f)+-2*x*(im*f+-1*c*log(f))))*exp(-1*im*d+b^2*(-4*c*log(f)+4*im*f)^-1*log(f)^2)+1//4*pi^(1/2)*f^a*(im*f+c*log(f))^-1//2*erfi((1/2)*(im*f+c*log(f))^-1//2*(b*log(f)+2*x*(im*f+c*log(f))))*exp(im*d+-1*b^2*(4*im*f+4*c*log(f))^-1*log(f)^2))
@test integrate(f^(a+b*x+c*x^2)*cos(d+f*x^2)^2, x) == :(-1//8*pi^(1/2)*f^a*(-1*c*log(f)+2*im*f)^-1//2*erf((1/2)*(-1*c*log(f)+2*im*f)^-1//2*(b*log(f)+-2*x*(-1*c*log(f)+2*im*f)))*exp(-2*im*d+b^2*(-4*c*log(f)+8*im*f)^-1*log(f)^2)+1//4*pi^(1/2)*c^-1//2*f^(a+-1//4*b^2*c^-1)*log(f)^-1//2*erfi((1/2)*c^-1//2*log(f)^(1/2)*(b+2*c*x))+1//8*pi^(1/2)*f^a*(c*log(f)+2*im*f)^-1//2*erfi((1/2)*(c*log(f)+2*im*f)^-1//2*(b*log(f)+2*x*(c*log(f)+2*im*f)))*exp(2*im*d+-1*b^2*(4*c*log(f)+8*im*f)^-1*log(f)^2))
@test integrate(f^(a+b*x+c*x^2)*cos(d+f*x^2)^3, x) == :(-3//16*pi^(1/2)*f^a*(im*f+-1*c*log(f))^-1//2*erf((1/2)*(im*f+-1*c*log(f))^-1//2*(b*log(f)+-2*x*(im*f+-1*c*log(f))))*exp(-1*im*d+b^2*(-4*c*log(f)+4*im*f)^-1*log(f)^2)+-1//16*pi^(1/2)*f^a*(-1*c*log(f)+3*im*f)^-1//2*erf((1/2)*(-1*c*log(f)+3*im*f)^-1//2*(b*log(f)+-2*x*(-1*c*log(f)+3*im*f)))*exp(-3*im*d+b^2*(-4*c*log(f)+12*im*f)^-1*log(f)^2)+1//16*pi^(1/2)*f^a*(c*log(f)+3*im*f)^-1//2*erfi((1/2)*(c*log(f)+3*im*f)^-1//2*(b*log(f)+2*x*(c*log(f)+3*im*f)))*exp(3*im*d+-1*b^2*(4*c*log(f)+12*im*f)^-1*log(f)^2)+3//16*pi^(1/2)*f^a*(im*f+c*log(f))^-1//2*erfi((1/2)*(im*f+c*log(f))^-1//2*(b*log(f)+2*x*(im*f+c*log(f))))*exp(im*d+-1*b^2*(4*im*f+4*c*log(f))^-1*log(f)^2))
@test integrate(f^(a+b*x+c*x^2)*cos(d+e*x+f*x^2), x) == :(1//4*pi^(1/2)*f^a*(im*f+c*log(f))^-1//2*erfi((1/2)*(im*f+c*log(f))^-1//2*(im*e+b*log(f)+2*x*(im*f+c*log(f))))*exp(im*d+(e+-1*im*b*log(f))^2*(4*im*f+4*c*log(f))^-1)+1//4*pi^(1/2)*f^a*(im*f+-1*c*log(f))^-1//2*erf((1/2)*(im*f+-1*c*log(f))^-1//2*(im*e+-1*b*log(f)+2*x*(im*f+-1*c*log(f))))*exp(-1*im*d+-1*(e+im*b*log(f))^2*(-4*c*log(f)+4*im*f)^-1))
@test integrate(f^(a+b*x+c*x^2)*cos(d+e*x+f*x^2)^2, x) == :(1//4*pi^(1/2)*c^-1//2*f^(a+-1//4*b^2*c^-1)*log(f)^-1//2*erfi((1/2)*c^-1//2*log(f)^(1/2)*(b+2*c*x))+1//8*pi^(1/2)*f^a*(c*log(f)+2*im*f)^-1//2*erfi((1/2)*(c*log(f)+2*im*f)^-1//2*(b*log(f)+2*im*e+2*x*(c*log(f)+2*im*f)))*exp((2*e+-1*im*b*log(f))^2*(4*c*log(f)+8*im*f)^-1+2*im*d)+1//8*pi^(1/2)*f^a*(-1*c*log(f)+2*im*f)^-1//2*erf((1/2)*(-1*c*log(f)+2*im*f)^-1//2*(-1*b*log(f)+2*im*e+2*x*(-1*c*log(f)+2*im*f)))*exp(-1*(2*e+im*b*log(f))^2*(-4*c*log(f)+8*im*f)^-1+-2*im*d))
@test integrate(f^(a+b*x+c*x^2)*cos(d+e*x+f*x^2)^3, x) == :(1//16*pi^(1/2)*f^a*(c*log(f)+3*im*f)^-1//2*erfi((1/2)*(c*log(f)+3*im*f)^-1//2*(b*log(f)+2*x*(c*log(f)+3*im*f)+3*im*e))*exp(-1*(b*log(f)+3*im*e)^2*(4*c*log(f)+12*im*f)^-1+3*im*d)+1//16*pi^(1/2)*f^a*(-1*c*log(f)+3*im*f)^-1//2*erf((1/2)*(-1*c*log(f)+3*im*f)^-1//2*(-1*b*log(f)+2*x*(-1*c*log(f)+3*im*f)+3*im*e))*exp(-1*(3*e+im*b*log(f))^2*(-4*c*log(f)+12*im*f)^-1+-3*im*d)+3//16*pi^(1/2)*f^a*(im*f+c*log(f))^-1//2*erfi((1/2)*(im*f+c*log(f))^-1//2*(im*e+b*log(f)+2*x*(im*f+c*log(f))))*exp(im*d+(e+-1*im*b*log(f))^2*(4*im*f+4*c*log(f))^-1)+3//16*pi^(1/2)*f^a*(im*f+-1*c*log(f))^-1//2*erf((1/2)*(im*f+-1*c*log(f))^-1//2*(im*e+-1*b*log(f)+2*x*(im*f+-1*c*log(f))))*exp(-1*im*d+-1*(e+im*b*log(f))^2*(-4*c*log(f)+4*im*f)^-1))
@test integrate(F^(c*(a+b*x))*(f+f*sin(d+e*x))^2, x) == :(F^(a*c+b*c*x)*b^-1*c^-1*f^2*log(F)^-1+-2*e*F^(a*c+b*c*x)*f^2*(e^2+b^2*c^2*log(F)^2)^-1*cos(d+e*x)+-2*e*F^(a*c+b*c*x)*f^2*(4*e^2+b^2*c^2*log(F)^2)^-1*cos(d+e*x)*sin(d+e*x)+b*c*F^(a*c+b*c*x)*f^2*(4*e^2+b^2*c^2*log(F)^2)^-1*sin(d+e*x)^2*log(F)+2*b*c*F^(a*c+b*c*x)*f^2*(e^2+b^2*c^2*log(F)^2)^-1*log(F)*sin(d+e*x)+2*F^(a*c+b*c*x)*b^-1*c^-1*e^2*f^2*(4*e^2+b^2*c^2*log(F)^2)^-1*log(F)^-1)
@test integrate(F^(c*(a+b*x))*(f+f*sin(d+e*x)), x) == :(f*F^(a*c+b*c*x)*b^-1*c^-1*log(F)^-1+-1*e*f*F^(a*c+b*c*x)*(e^2+b^2*c^2*log(F)^2)^-1*cos(d+e*x)+b*c*f*F^(a*c+b*c*x)*(e^2+b^2*c^2*log(F)^2)^-1*log(F)*sin(d+e*x))
@test integrate(F^(c*(a+b*x))*(f+f*cos(d+e*x))^2, x) == :(F^(a*c+b*c*x)*b^-1*c^-1*f^2*log(F)^-1+2*e*F^(a*c+b*c*x)*f^2*(e^2+b^2*c^2*log(F)^2)^-1*sin(d+e*x)+2*e*F^(a*c+b*c*x)*f^2*(4*e^2+b^2*c^2*log(F)^2)^-1*cos(d+e*x)*sin(d+e*x)+b*c*F^(a*c+b*c*x)*f^2*(4*e^2+b^2*c^2*log(F)^2)^-1*cos(d+e*x)^2*log(F)+2*b*c*F^(a*c+b*c*x)*f^2*(e^2+b^2*c^2*log(F)^2)^-1*cos(d+e*x)*log(F)+2*F^(a*c+b*c*x)*b^-1*c^-1*e^2*f^2*(4*e^2+b^2*c^2*log(F)^2)^-1*log(F)^-1)
@test integrate(F^(c*(a+b*x))*(f+f*cos(d+e*x)), x) == :(e*f*F^(a*c+b*c*x)*(e^2+b^2*c^2*log(F)^2)^-1*sin(d+e*x)+f*F^(a*c+b*c*x)*b^-1*c^-1*log(F)^-1+b*c*f*F^(a*c+b*c*x)*(e^2+b^2*c^2*log(F)^2)^-1*cos(d+e*x)*log(F))

max_stars_repo_path: test/4 Trig functions/4.7 Miscellaneous/4.7.6 f^(a+b x+c x^2) trig(d+e x+f x^2)^n.jl