text:
### Struct variable
struct PowerFlowAC
    main::Array{Float64,2}
    flow::Array{Float64,2}
    generation::Array{Float64,2}
    iterations::Int64
end

### AC power flow
@inbounds function runacpf(system, num, settings, info)

    ########## Pre-processing ##########
    busi, typei, Pload, Qload, Gshunt, Bshunt, Vini, Tini,
    Pgeni, Qgeni, Qmaxi, Qmini, Vgen, genOn,
    resistance, reactance, charging, transTap, transShift, branchOn = read_acsystem(system)

    main = fill(0.0, num.Nbus, 11)
    flow =  fill(0.0, num.Nbranch, 14)
    gene = fill(0.0, num.Ngen, 3)
    Va, Vm, Pinj, Qinj, Pbus, Qbus, Pshunt, Qshunt,
    Pij, Qij, Pji, Qji, Qcharging, Ploss, Qloss, Imij, Iaij, Imji, Iaji,
    Pgen, Qgen = write_acsystem(system, num, main, flow, gene)

  algtime = @elapsed begin
    ########## Convert in integers ##########
    geni = convert(Array{Int64,1}, system.generator[:, 1])
    fromi = convert(Array{Int64,1}, system.branch[:, 2])
    toi = convert(Array{Int64,1}, system.branch[:, 3])

    ########## Numbering and slack bus ##########
    numbering = false
    bus = collect(1:num.Nbus)
    slack = 0
    for i = 1:num.Nbus
        if bus[i] != busi[i]
            numbering = true
        end
        if typei[i] == 3
            slack = bus[i]
        end
    end
    if numbering
        println("The new bus numbering is running.")
    end
    if slack == 0
        slack = 1
        println("The slack bus is not found. Slack bus is the first bus.")
    end

    gen_bus = renumber(geni, busi, bus, numbering)
    from = renumber(fromi, busi, bus, numbering)
    to = renumber(toi, busi, bus, numbering)

    ########## Generation ##########
    slackLimit = copy(slack)
    multiple = fill(0, num.Nbus)
    isMultiple = false
    limit = fill(1, num.Ngen)
    Qmin = fill(0.0, num.Ngen)
    Qmax = fill(0.0, num.Ngen)
    type = fill(1, num.Nbus)

    for (k, i) in enumerate(gen_bus)
        if genOn[k] == 1
            Pbus[i] += Pgeni[k] / system.basePower
            Qbus[i] += Qgeni[k] / system.basePower
            Pgen[k] = Pgeni[k] / system.basePower
            Qgen[k] = Qgeni[k] / system.basePower
            Qmin[k] = Qmini[k] / system.basePower
            Qmax[k] = Qmaxi[k] / system.basePower
            Vini[i] = Vgen[k]
            type[i] = 2

            multiple[i] += 1
            if multiple[i] != 1
                isMultiple = true
            end
        end
    end
    type[slack] = 3

    ########## Ybus matrix ##########
    Ybus, YbusT, Ytt, Yff, Yft, Ytf, admittance, tap, shunt = ybusac(system, num, bus, from, to, branchOn,
        Gshunt, Bshunt, resistance, reactance, charging, transTap, transShift)

    ########## Solve the system ##########
    Vc = Vini .* exp.(im * (pi / 180) * Tini)
    iter = 0
    while settings.reactive[2]
        if  settings.reactive[1] && settings.reactive[2]
            Vc = Vini .* exp.(im * (pi / 180) * Tini)
        end

        if settings.algorithm == "gs"
            Vc, iter = gauss_seidel(system, num, settings, Ybus, YbusT, slack, Vc, Pbus, Qbus, Pload, Qload, Vini, type, iter)
        end
        if settings.algorithm == "nr"
            Vc, iter = newton_raphson(system, num, settings, Ybus, YbusT, slack, Vc, Pbus, Qbus, Pload, Qload, type, iter)
        end
        if settings.algorithm == "fnrbx" || settings.algorithm == "fnrxb"
            Vc, iter = fast_newton_raphson(system, num, settings, branchOn, Ybus, YbusT, slack, Vc, Pbus, Qbus,
                Pload, Qload, type, resistance, reactance, transShift, Gshunt, Bshunt, charging, transTap, from, to, iter)
        end

        for i = 1:num.Nbus
            Sshunt = Vc[i] * conj(Vc[i] * shunt[i])
            Pshunt[i] = real(Sshunt) * system.basePower
            Qshunt[i] = imag(Sshunt) * system.basePower

            I = 0.0 + im * 0.0
            for j in Ybus.colptr[i]:(Ybus.colptr[i + 1] - 1)
                row = Ybus.rowval[j]
                I += conj(YbusT[row, i]) * conj(Vc[row])
            end
            Si::ComplexF64 = I * Vc[i]
            Pinj[i] = real(Si) * system.basePower
            Qinj[i] = imag(Si) * system.basePower
            if type[i] != 1
                Qbus[i] = (Qinj[i] + Qload[i]) / system.basePower
            end
        end
        Pbus[slack] = (Pinj[slack] + Pload[slack]) / system.basePower

        if !isMultiple
            for i = 1:num.Ngen
                if genOn[i] == 1
                    Pgen[i] = Pbus[gen_bus[i]]
                    Qgen[i] = Qbus[gen_bus[i]]
                end
            end
        else
            Qmintotal = fill(0.0, num.Nbus); Qmaxtotal = fill(0.0, num.Nbus)
            QminInf = fill(0.0, num.Nbus); QmaxInf = fill(0.0, num.Nbus)
            QminNew = copy(Qmin);  QmaxNew = copy(Qmax)
            Qgentotal = fill(0.0, num.Nbus)
            for i = 1:num.Ngen
                if genOn[i] == 1
                    j = gen_bus[i]
                    if !isinf(Qmin[i])
                        Qmintotal[j] += Qmin[i]
                    end
                    if !isinf(Qmax[i])
                        Qmaxtotal[j] += Qmax[i]
                    end
                    Qgentotal[j] += ((Qinj[j] + Qload[j]) / system.basePower) / multiple[j]
                end
            end
            for i = 1:num.Ngen
                if genOn[i] == 1
                    j = gen_bus[i]
                    if Qmin[i] == Inf
                        QminInf[i] = abs(Qgentotal[j]) + abs(Qmintotal[j]) + abs(Qmaxtotal[j])
                    end
                    if Qmin[i] == -Inf
                        QminInf[i] = -abs(Qgentotal[j]) - abs(Qmintotal[j]) - abs(Qmaxtotal[j])
                    end
                    if Qmax[i] == Inf
                        QmaxInf[i] = abs(Qgentotal[j]) + abs(Qmintotal[j]) + abs(Qmaxtotal[j])
                    end
                    if Qmax[i] == -Inf
                        QmaxInf[i] = -abs(Qgentotal[j]) - abs(Qmintotal[j]) - abs(Qmaxtotal[j])
                    end
                end
            end

            for i = 1:num.Ngen
                if genOn[i] == 1
                    j = gen_bus[i]
                    if isinf(Qmin[i])
                        Qmintotal[j] += QminInf[i]
                        QminNew[i] = QminInf[i]
                    end
                    if isinf(Qmax[i])
                        Qmaxtotal[j] += QmaxInf[i]
                        QmaxNew[i] =  QmaxInf[i]
                    end
                end
            end
            flag = true
            tempslack = 0
            for i = 1:num.Ngen
                if genOn[i] == 1
                    j = gen_bus[i]
                    if system.basePower * abs(Qmintotal[j] - Qmaxtotal[j]) > 10 * eps(Float64)
                        Qgen[i] = QminNew[i] + ((Qgentotal[j] - Qmintotal[j]) / (Qmaxtotal[j] - Qmintotal[j])) * (QmaxNew[i] - QminNew[i])
                    else
                        Qgen[i] = QminNew[i] + (Qgentotal[j] - Qmintotal[j]) / multiple[j]
                    end
                    if gen_bus[i] == slack && flag == false
                        Pgen[tempslack] -= Pgen[i]
                    end
                    if gen_bus[i] == slack && flag == true
                        Pgen[i] = Pbus[slack]
                        tempslack = i
                        flag = false
                    end
                end
            end
        end

        if settings.reactive[1] && settings.reactive[2]
            settings.reactive[2] = false
            for i = 1:num.Ngen
                if genOn[i] == 1
                    j = gen_bus[i]
                    if limit[i] == 1 && (Qgen[i] < Qmin[i] || Qgen[i] > Qmax[i])
                        settings.reactive[2] = true
                        j = gen_bus[i]
                        if Qgen[i] < Qmin[i]
                            limit[i] = 2
                            Qbus[j] -= Qgen[i]
                            Qgen[i] = Qmin[i]
                            Qbus[j] += Qgen[i]
                        end
                        if Qgen[i] > Qmax[i]
                            limit[i] = 3
                            Qbus[j] -= Qgen[i]
                            Qgen[i] = Qmax[i]
                            Qbus[j] += Qgen[i]
                        end

                        type[j] = 1
                        if j == slack
                            for k = 1:num.Nbus
                                if type[k] == 2
                                    println("Bus $(trunc(Int, system.bus[k, 1])) is the new slack bus.")
                                    slack = bus[k]
                                    type[k] = 3
                                    settings.reactive[3] = true
                                    break
                                end
                            end
                        end
                    end
                end
            end
        else
            settings.reactive[2] = false
        end
    end

    if settings.reactive[3]
        T = angle(Vc[slackLimit])
        for i = 1:num.Nbus
            Vc[i] = abs(Vc[i]) * exp(im * (angle(Vc[i]) - T +  (pi / 180) * Tini[slackLimit]))
        end
    end

    ########## Post-processing ##########
    for i = 1:num.Nbus
        Vm[i] = abs(Vc[i])
        Va[i] = (180 / pi) * angle(Vc[i])
        Pbus[i] = Pbus[i] * system.basePower
        Qbus[i] = Qbus[i] * system.basePower
    end

    for i = 1:num.Ngen
        Pgen[i] = Pgen[i] * system.basePower
        Qgen[i] = Qgen[i] * system.basePower
    end

    acflow(system.basePower, num, from, to, branchOn, Vc, Yff, Yft, Ytt, Ytf,
        admittance, tap, resistance, reactance, charging,
        Pij, Qij, Pji, Qji, Qcharging, Ploss, Qloss, Imij, Iaij, Imji, Iaji)
  end #algtime

    ########## Results ##########
    results = PowerFlowAC(main, flow, gene, iter)
    header, group = results_flowac(system, num, settings, results, slack, limit, algtime)
    if !isempty(settings.save)
        savedata(results, system; info = info, group = group, header = header, path = settings.save)
    end

    return results
end

### Read data
function read_acsystem(system)
    busi = @view(system.bus[:, 1])
    typei = @view(system.bus[:, 2])
    Pload = @view(system.bus[:, 3])
    Qload = @view(system.bus[:, 4])
    Gshunt = @view(system.bus[:, 5])
    Bshunt =  @view(system.bus[:, 6])
    Vini =  @view(system.bus[:, 8])
    Tini = @view(system.bus[:, 9])

    Pgeni = @view(system.generator[:, 2])
    Qgeni = @view(system.generator[:, 3])
    Qmaxi = @view(system.generator[:, 4])
    Qmini = @view(system.generator[:, 5])
    Vgen = @view(system.generator[:, 6])
    genOn = @view(system.generator[:, 8])

    resistance = @view(system.branch[:, 4])
    reactance = @view(system.branch[:, 5])
    charging = @view(system.branch[:, 6])
    transTap = @view(system.branch[:, 10])
    transShift = @view(system.branch[:, 11])
    branchOn = @view(system.branch[:, 12])

    return busi, typei, Pload, Qload, Gshunt, Bshunt, Vini, Tini,
        Pgeni, Qgeni, Qmaxi, Qmini, Vgen, genOn,
        resistance, reactance, charging, transTap, transShift, branchOn
end

### Write data
function write_acsystem(system, num, main, flow, gene)
    for i = 1:num.Nbus
        main[i, 1] = system.bus[i, 1]
        main[i, 8] = system.bus[i, 3]
        main[i, 9] = system.bus[i, 4]
    end
    Vm = @view(main[:, 2])
    Va = @view(main[:, 3])
    Pinj = @view(main[:, 4])
    Qinj = @view(main[:, 5])
    Pbus = @view(main[:, 6])
    Qbus = @view(main[:, 7])
    Pshunt = @view(main[:, 10])
    Qshunt = @view(main[:, 11])

    for i = 1:num.Nbranch
        flow[i, 1] = system.branch[i, 1]
        flow[i, 2] = system.branch[i, 2]
        flow[i, 3] = system.branch[i, 3]
    end
    Pij = @view(flow[:, 4])
    Qij = @view(flow[:, 5])
    Pji = @view(flow[:, 6])
    Qji = @view(flow[:, 7])
    Qcharging = @view(flow[:, 8])
    Ploss = @view(flow[:, 9])
    Qloss = @view(flow[:, 10])
    Imij = @view(flow[:, 11])
    Iaij = @view(flow[:, 12])
    Imji = @view(flow[:, 13])
    Iaji = @view(flow[:, 14])

    for i = 1:num.Ngen
        gene[i, 1] = system.generator[i, 1]
    end
    Pgen = @view(gene[:, 2])
    Qgen = @view(gene[:, 3])

    return Va, Vm, Pinj, Qinj, Pbus, Qbus, Pshunt, Qshunt,
        Pij, Qij, Pji, Qji, Qcharging, Ploss, Qloss, Imij, Iaij, Imji, Iaji,
        Pgen, Qgen
end

max_stars_repo_name: mcosovic/JuliaGrid

max_issues_repo_name: mcosovic/JuliaGrid

max_forks_repo_name: mcosovic/JuliaGrid


text:
using NBodySimulator, StaticArrays

function generate_bodies_in_line(n::Int, m::Real, v_dev::Real, L::Real)
    dL = L / (ceil(n^(1 / 3)))
    n_line = floor(Int, L / dL)
    rng = MersenneTwister(n);
    velocities = v_dev * randn(rng, Float64, (3, n_line))
    bodies = MassBody[]
    x = y = L / 2
    for i = 1:n_line
        r = SVector(x, y, i * dL)
        v = SVector{3}(velocities[:,i])
        body = MassBody(r, v, m)
        push!(bodies, body)
    end
    return bodies
end

function generate_random_directions(n::Int)
    theta = acos.(1 - 2 * rand(n));
    phi = 2 * pi * rand(n);
    directions = [@SVector [sin(theta[i]) .* cos(phi[i]), sin(theta[i]) .* sin(phi[i]), cos(theta[i])] for i = 1:n]
end

units = :real
units = :reduced

const T = 120.0 # °K
const T0 = 90.0
const kb = 1.38e-23 # J/K
const ϵ = T * kb
const σ = 3.4e-10 # m
const ρ = 1374 # kg/m^3
const m = 39.95 * 1.6747 * 1e-27 # kg
const N = 216#floor(Int, ρ * L^3 / m)
const L = (m*N/ρ)^(1/3)#10.229σ
const R = 0.5*L
const v_dev = sqrt(kb * T / m)
const τ = 0.5e-15 # σ/v
const t1 = 0τ
const t2 = 30000τ
#bodies = generate_bodies_randomly(N, m, v_dev, L)
bodies = generate_bodies_in_cell_nodes(N, m, v_dev, L)
#bodies = generate_bodies_in_line(N, m, v_dev, L)
jl_parameters = LennardJonesParameters(ϵ, σ, R)
pbc = CubicPeriodicBoundaryConditions(L)
#thermostat = AndersenThermostat(0.02, T, kb)
lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => jl_parameters));
simulation = NBodySimulation(lj_system, (t1, t2), pbc, kb);
#result = run_simulation(simulation, Tsit5())
result = @time run_simulation(simulation, VelocityVerlet(), dt=τ)

#=
using Plots
import GR
(rs, grf) = rdf(result)
(ts, dr2) = msd(result)
plot(rs/σ, grf, xlim=[0, 0.4999L/σ], label=["Radial distribution function"],ylabel="g(r)", xlabel="r/sigma")

using JLD
time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
Nactual = length(bodies)
timesteps = round(length(result.solution.t))
#save("D:/water $Nactual molecules $timesteps steps.jld", "rs", rs, "grf", grf, "ts", ts, "dr2", dr2)
save("D:/liquid argon $Nactual molecules $timesteps steps $time_now.jld", "rs", rs, "grf", grf, "ts", ts, "dr2", dr2, "e_tot", e_tot, "e_kin", e_kin, "e_pot", e_pot)
=#
#=
using Plots
import GR
time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
Nactual = length(bodies)
timesteps = round(length(result.solution.t))
@time animate(result, "D:/$Nactual liquid argon particles with $timesteps timesteps $time_now.gif")

#plot(simulation)
=#

#=
(rs, grf) = rdf(result)
(ts, dr2) = msd(result)

t = t1:τ:result.solution.t[end-1]
temper = temperature.(result, t)

time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
Nactual = length(bodies)
timesteps = round(length(result.solution.t))

using MAT
file = matopen("d:/liquid argon SI units and $timesteps timesteps_$time_now.mat", "w")
write(file, "t", collect(t))
write(file, "temper", temper)
write(file, "rs", rs)
write(file, "ts", ts)
write(file, "grf", grf)
write(file, "dr2", dr2)
close(file)
=#

ext: jl


text:
module OceanParameterizations

if VERSION < v"1.5"
    error("OceanParameterizations.jl requires Julia v1.5 or newer.")
end

export
    # Utils
    coarse_grain, Dᶠ, Dᶜ,
    ZeroMeanUnitVarianceScaling, MinMaxScaling, scale, unscale,
    predict, animate_gif,

    # GaussianProcesses
    gp_model, best_kernel, get_kernel, euclidean_distance, derivative_distance, antiderivative_distance, error_per_gamma,

    # NeuralNetworks
    nn_model,

    # PhysicalParameterizations
    closure_kpp_full_evolution,
    closure_tke_full_evolution

using LinearAlgebra
using Printf
using Statistics
using Logging

using Flux
using JLD2
using NCDatasets
using Plots
using Oceananigans.Utils

using Oceananigans.Grids: Center, Face

import Base.inv

include("differentiation_operators.jl")
include("predict.jl")

# submodules
include("DataWrangling/DataWrangling.jl")
include("GaussianProcesses/GaussianProcesses.jl")
include("NeuralNetworks/NeuralNetworks.jl")
include("PhysicalParameterizations/PhysicalParameterizations.jl")

using .DataWrangling
using .GaussianProcesses
using .NeuralNetworks
using .PhysicalParameterizations

end # module


text:
@testset "Symmetric TSP: pr76, by coordinates" begin
    tsp = TSPLIB.readTSPLIB(:pr76)
    ap = AlgorithmParameters(timeLimit=1)
    result = solve_tsp(tsp, ap)
    @show result.cost
    @show result.routes
    @test result.cost <= 108159 * 1.02
end

@testset "Asymmetric TSP: br17, by dist_mtx" begin
    tsp = TSPLIB.readTSP("br17.atsp")
    ap = AlgorithmParameters(timeLimit=1)
    result = solve_tsp(tsp, ap; use_dist_mtx=true)
    @show result.cost
    @show result.routes
    @test result.cost <= 39 * 1.02
end

@testset "Asymmetric TSP: ftv64, by dist_mtx" begin
    tsp = TSPLIB.readTSP("ftv64.atsp")
    ap = AlgorithmParameters(timeLimit=1)
    result = solve_tsp(tsp, ap; use_dist_mtx=true)
    @show result.cost
    @show result.routes
    @test result.cost <= 1839 * 1.02
end

@testset "x, y, dist_mtx TSP" begin
    tsp = TSPLIB.readTSPLIB(:pr76)
    dist_mtx = tsp.weights
    x = tsp.nodes[:, 1]
    y = tsp.nodes[:, 2]

    ap = AlgorithmParameters(timeLimit=1)
    result1 = solve_tsp(x, y, ap)
    result2 = solve_tsp(dist_mtx, ap)
    result3 = solve_tsp(dist_mtx, ap; x_coordinates=x, y_coordinates=y)

    @test result1.cost == result2.cost == result3.cost
end


text:
"""
Baysian Structure learning seeks to maximize P(G|D)
In the Bayesian fashion, we can provide a prior over the parameters in our learning network.
This is described using a Dirichlet Prior.
"""
abstract type DirichletPrior end

"""
A uniform Dirichlet prior such that all α are the same

Defaults to the popular K2 prior, α = 1, which is similar to Laplace Smoothing

    https://en.wikipedia.org/wiki/Additive_smoothing
"""
struct UniformPrior <: DirichletPrior
    α::Float64
    UniformPrior(α::Float64=1.0) = new(α)
end

Base.print(io::IO, p::UniformPrior) = Base.print(io, "UniformPrior(%.2f)", p.α)

Base.get(p::UniformPrior, ncategories::Integer) = fill(p.α, ncategories)
function Base.get{I<:Integer, J<:Integer}(p::UniformPrior,
    var_index::Integer,
    ncategories::AbstractVector{I}, # [nvars]
    parents::AbstractVector{J},    # [nvars]
    )

    r = ncategories[var_index]
    q = isempty(parents) ? 1 : prod(ncategories[parents])
    α = p.α

    fill(α, r, q)
end

"""
Assigns equal scores to Markov equivalent structures

    α_ijk = x/{q_i * r_i} for each j, k and some given x

see DMU section 2.4.3
"""
struct BDeuPrior <: DirichletPrior
    x::Float64
    BDeuPrior(x::Float64=1.0) = new(x)
end

Base.print(io::IO, p::BDeuPrior) = @printf(io, "BDeuPrior(%.2f)", p.x)

Base.get(p::BDeuPrior, ncategories::Integer) = fill(p.x/ncategories, ncategories)
function Base.get{I<:Integer, J<:Integer}(p::BDeuPrior,
    var_index::Integer,
    ncategories::AbstractVector{I}, # [nvars]
    parents::AbstractVector{J},    # [nvars]
    )

    x = p.x
    r = ncategories[var_index]
    q = isempty(parents) ? 1 : prod(ncategories[parents])

    α = x / (r*q)

    fill(α, r, q)
end


text:
using ConcatOperator
using Base.Test

@test "abc" ++ "def" == "abcdef"
@test "ab" ++ "cd" ++ "ef" == "abcdef"
@test "a" ++ 'b' ++ "cd" ++ 'e' ++ "f" == "abcdef"
@test 'π' ++ " ≠" ++ utf16(" α ±") ++ utf32(" β") == "π ≠ α ± β"

@test [1,2] ++ [3,4] == [1,2,3,4]
@test [1,2] ++ [3.0,4.0] == [1.0,2.0,3.0,4.0]


text:
# Tutorial 2.6. "Superconductors": orbitals, conservation laws and symmetries
# ===========================================================================
#
# Physics background
# ------------------
# - conductance of a NS-junction (Andreev reflection, superconducting gap)
#
# Kwant features highlighted
# --------------------------
# - Implementing electron and hole ("orbital") degrees of freedom
#   using conservation laws.
# - Use of discrete symmetries to relate scattering states.

tau_x = [0 1; 1 0]
tau_y = [0  -1im; 1im  0]
tau_z = [1   0  ;  0  -1]

function make_system(;a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
                mu=0.4, Delta=0.1, Deltapos=4, t=1.0)
    # Start with an empty tight-binding system. On each site, there
    # are now electron and hole orbitals, so we must specify the
    # number of orbitals per site. The orbital structure is the same
    # as in the Hamiltonian.
    lat = kwant.lattice.square(norbs=2)
    syst = kwant.Builder()

    #### Define the scattering region. ####
    # The superconducting order parameter couples electron and hole orbitals
    # on each site, and hence enters as an onsite potential.
    # The pairing is only included beyond the point 'Deltapos' in the scattering region.
    # syst[(lat(x, y) for x in range(Deltapos) for y in range(W))] = (4 * t - mu) * tau_z
    # syst[(lat(x, y) for x in range(Deltapos, L) for y in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x
    #
    # The tunnel barrier
    syst[[lat(x, y) for x in range(barrierpos[1], stop=barrierpos[2])
         for y in range(0,length=W)]] = (4 * t + barrier - mu) * tau_z

    # Hoppings
    syst[lat.neighbors()] = -t * tau_z
    #### Define the leads. ####
    # Left lead - normal, so the order parameter is zero.
    sym_left = kwant.TranslationalSymmetry((-a, 0))
    # Specify the conservation law used to treat electrons and holes separately.
    # We only do this in the left lead, where the pairing is zero.
    lead0 = kwant.Builder(sym_left, conservation_law=-tau_z, particle_hole=tau_y)
    lead0[[lat(0, j) for j in range(0,length=W)]] = (4 * t - mu) * tau_z
    lead0[lat.neighbors()] = -t * tau_z
    # Right lead - superconducting, so the order parameter is included.
    sym_right = kwant.TranslationalSymmetry((a, 0))
    lead1 = kwant.Builder(sym_right)
    lead1[[lat(0, j) for j in range(0,length=W)]] = (4 * t - mu) * tau_z + Delta * tau_x
    lead1[lat.neighbors()] = -t * tau_z

    #### Attach the leads and return the system. ####
    syst.attach_lead(lead0)
    syst.attach_lead(lead1)

    return syst
end

function plot_conductance(syst; energies)
    # Compute conductance
    data = []
    for energy in energies
        smatrix = kwant.smatrix(syst, energy)
        # Conductance is N - R_ee + R_he
        append!(data,smatrix.submatrix((0, 0), (0, 0)).shape[1] -
                    smatrix.transmission((0, 0), (0, 0)) +
                    smatrix.transmission((0, 1), (0, 0)))
    end
    pyplot.figure()
    pyplot.plot(energies, data)
    pyplot.xlabel("energy [t]")
    pyplot.ylabel("conductance [e^2/h]")
    pyplot.gcf()
end

function check_PHS(syst)
    # Scattering matrix
    s = kwant.smatrix(syst, 0)
    # Electron to electron block
    s_ee = convert(Array,s.submatrix((0,0), (0,0)))
    # Hole to hole block
    s_hh = convert(Array,s.submatrix((0,1), (0,1)))
    println("s_ee: \n", round.(s_ee; digits=3))
    println("s_hh: \n", round.(s_hh[end:-1:1, end:-1:1], digits=3))
    println("s_ee - s_hh^*: \n",
          round.((s_ee - conj(s_hh[end:-1:1, end:-1:1])), digits=3), "\n")
    # Electron to hole block
    s_he = convert(Array,s.submatrix((0,1), (0,0)))
    # Hole to electron block
    s_eh = convert(Array,s.submatrix((0,0), (0,1)))
    println("s_he: \n", round.(s_he, digits=3))
    println("s_eh: \n", round.(s_eh[end:-1:1,end:-1:1], digits=3))
    println("s_he + s_eh^*: \n",
          round.(s_he + conj(s_eh[end:-1:1,end:-1:1]), digits=3))
end


text:
# Function to obtain the conditional means for given time series,
# parameters, and model specification

# Input:
# y:        Time series
# model:    Model specification
# θ:        Parameters
# X:        Matrix of regressors

# Output:
# λ:        Conditional mean

# Core operations are put in functions for speedup

function CoreIID(β0::Float64, logl::Bool, T::Int64)::Vector{Float64}
    if logl
        return fill(exp(β0), T)
    else
        return fill(Float64(β0), T)
    end
end

function CoreIID(β0::Float64, X::Matrix{Float64}, η::Vector{Float64}, logl::Bool, T::Int64)::Vector{Float64}
    out = fill(β0, T)
    out += (η'X)[1, :]
    if logl
        return exp.(out)
    else
        return out
    end
end

function LinPred(y::Array{T,1} where T<:Integer,
    model::IIDModel,
    θ::parameter,
    initiate::String = "first")::Vector{Float64}
    if !(initiate in ["first", "intercept", "marginal"])
        println("Initiation not specified correctly - Changed to 'fist'")
        initiate = "first"
    end

    T = length(y)
    r = length(model.external)

    logl = (model.link == "Log")

    β0 = Float64(θ.β0)

    if r > 0
        X = Matrix{Float64}(model.X)
        η = Vector{Float64}(θ.η)
        return CoreIID(β0, X, η, logl, T)
    else
        return CoreIID(β0, logl, T)
    end
end

function CoreARCH(y::Vector{Int64},
    β0::Float64,
    α::Vector{Float64},
    pastObs::Vector{Int64},
    P::Int64,
    logl::Bool,
    T::Int64)::Vector{Float64}

    out = fill(β0, T)
    if logl
        y = @. log(y + 1)
    end
    @inbounds for i = pastObs
        out[P+1:end] .+= α[i] .* y[P+1-i:end-i]
    end
    if logl
        return exp.(out)
    else
        return out
    end
end

function CoreARCH(y::Vector{Int64},
    β0::Float64,
    α::Vector{Float64},
    pastObs::Vector{Int64},
    P::Int64,
    η::Vector{Float64},
    X::Matrix{Float64},
    logl::Bool,
    T::Int64)::Vector{Float64}

    out = fill(β0, T)
    if logl
        y = @. log(y + 1)
    end
    @inbounds for i = pastObs
        out[P+1:end] .+= α[i] .* y[P+1-i:end-i]
    end
    out += (η'*X)[1, :]

    if logl
        return exp.(out)
    else
        return out
    end
end

function LinPred(y::Array{T,1} where T<: Integer,
    model::INARCHModel,
    θ::parameter,
    initiate::String = "first")::Vector{Float64}

    if !(initiate in ["first", "intercept", "marginal"])
        println("Initiation not specified correctly - Changed to 'fist'")
        initiate = "first"
    end

    T = length(y)

    p = length(model.pastObs)
    if p == 0
        P = 0
    else
        P = maximum(model.pastObs)
    end
    r = length(model.external)

    logl = (model.link == "Log")
    lin = !logl

    β0 = Float64(θ.β0)
    α = Vector{Float64}(θ.α)

    if r > 0
        X = Matrix{Float64}(model.X)
        η = Vector{Float64}(θ.η)
        return CoreARCH(y, β0, α, model.pastObs, P, η, X, logl, T)
    else
        return CoreARCH(y, β0, α, model.pastObs, P, logl, T)
    end
end

# No regressors
function CoreGARCH(y::Vector{Int64},
    λinit::Float64,
    β0::Float64,
    α::Vector{Float64},
    pastObs::Vector{Int64},
    p::Int64,
    β::Vector{Float64},
    pastMean::Vector{Int64},
    q::Int64,
    M::Int64,
    logl::Bool,
    T::Int64)::Vector{Float64}

    out = fill(β0, T)
    out[1:M] .= λinit

    if logl
        y = @. log(y + 1)
    end
    for t = M+1:T
        for i = 1:p
            out[t] += α[i] * y[t-pastObs[i]]
        end

        for i = 1:q
            out[t] += β[i] * out[t - pastMean[i]]
        end
    end

    if logl
        return exp.(out)
    else
        return out
    end
end

# Only Internal regressors
function CoreGARCH(y::Vector{Int64},
    λinit::Float64,
    β0::Float64,
    α::Vector{Float64},
    pastObs::Vector{Int64},
    p::Int64,
    β::Vector{Float64},
    pastMean::Vector{Int64},
    q::Int64,
    M::Int64,
    ηI::Vector{Float64},
    rI::Int64,
    XI::Matrix{Float64},
    logl::Bool,
    T::Int64)::Vector{Float64}

    out = fill(β0, T)
    out[1:M] .= λinit

    if logl
        y = @. log(y + 1)
    end
    @inbounds for t = M+1:T
        for i = 1:p
            out[t] += α[i] * y[t-pastObs[i]]
        end

        for i = 1:q
            out[t] += β[i] * out[t - pastMean[i]]
        end

        for i = 1:rI
            out[t] += ηI[i] * XI[i, t]
        end
    end

    if logl
        return exp.(out)
    else
        return out
    end
end

# Both, internal and external regressors
function CoreGARCH(y::Vector{Int64},
    λinit::Float64,
    β0::Float64,
    α::Vector{Float64},
    pastObs::Vector{Int64},
    p::Int64,
    β::Vector{Float64},
    pastMean::Vector{Int64},
    q::Int64,
    M::Int64,
    ηI::Vector{Float64},
    rI::Int64,
    XI::Matrix{Float64},
    logl::Bool,
    T::Int64,
    ηE::Vector{Float64},
    XE::Matrix{Float64})::Vector{Float64}

    out = fill(β0, T)
    out[1:M] .= λinit

    if logl
        y = @. log(y + 1)
    end
    @inbounds for t = M+1:T
        for i = 1:p
            out[t] += α[i] * y[t-pastObs[i]]
        end

        for i = 1:q
            out[t] += β[i] * out[t - pastMean[i]]
        end

        for i = 1:rI
            out[t] += ηI[i] * XI[i, t]
        end
    end

    if logl
        out = exp.(out)
    end

    return out + (ηE'XE)[1, :]
end

# Only external regressors
function CoreGARCH(y::Vector{Int64},
    λinit::Float64,
    β0::Float64,
    α::Vector{Float64},
    pastObs::Vector{Int64},
    p::Int64,
    β::Vector{Float64},
    pastMean::Vector{Int64},
    q::Int64,
    M::Int64,
    logl::Bool,
    T::Int64,
    ηE::Vector{Float64},
    XE::Matrix{Float64})::Vector{Float64}

    out = fill(β0, T)
    out[1:M] .= λinit

    if logl
        y = @. log(y + 1)
    end
    @inbounds for t = M+1:T
        for i = 1:p
            out[t] += α[i] * y[t-pastObs[i]]
        end

        for i = 1:q
            out[t] += β[i] * out[t - pastMean[i]]
        end
    end

    if logl
        out = exp.(out)
    end

    return out + (ηE'XE)[1, :]
end

function LinPred(y::Array{T,1} where T<: Integer,
    model::INGARCHModel,
    θ::parameter,
    initiate::String = "first")::Vector{Float64}

    if !(initiate in ["first", "intercept", "marginal"])
        println("Initiation not specified correctly - Changed to 'fist'")
        initiate = "first"
    end

    T = length(y)

    p = length(model.pastObs)
    if p == 0
        P = 0
    else
        P = Int64(maximum(model.pastObs))
    end
    q = length(model.pastMean)
    if q == 0
        Q = 0
    else
        Q = Int64(maximum(model.pastMean))
    end

    M = Int64(maximum([P, Q]))
    r = length(model.external)

    logl = (model.link == "Log")
    lin = !logl

    β0 = Float64(θ.β0)
    α = Vector{Float64}(θ.α)
    β = Vector{Float64}(θ.β)

    if initiate == "first"
        λinit = Float64(y[1])
        if logl
            λinit = log(λinit + 1)
        end
    end
    if initiate == "marginal"
        λinit = β0/(1 - sum(α) - sum(β))
        if logl
            λinit = log(λinit)
        end
    end
    if initiate == "intercept"
        λinit = β0
        if logl
            λinit = log(λinit)
        end
    end

    if r > 0
        iE = findall(model.external)
        iI = setdiff(1:r, iE)
        rE = length(iE)
        rI = length(iI)
        X = Matrix{Float64}(model.X)
        η = Vector{Float64}(θ.η)
        if rE == r
            return CoreGARCH(y, λinit, β0, α, model.pastObs, p, β, model.pastMean, q, M, logl, T, η, X)
        end
        if rI == r
            return CoreGARCH(y, λinit, β0, α, model.pastObs, p, β, model.pastMean, q, M, η, r, X, logl, T)
        end

        return CoreGARCH(y, λinit, β0, α, model.pastObs, p, β, model.pastMean, q, M, η[iI], rI, X[iI, :], logl, T, η[iE], X[iE, :])
    else
        return CoreGARCH(y, λinit, β0, α, model.pastObs, p, β, model.pastMean, q, M, logl, T)
    end
end


text:
using .Plots

"""
    animate_draws(uc::UncertainCurve, bnds::Vector{UncertainBound}; n=20, fps=5, filepath=nothing, kw...)

Create a gif animation showing the first `n` random draws.
Requires the Plots.jl package.

Create an animation of the Monte-Carlo iterations when integrating the
`UncertainCurve` `uc` using the `UncertainBound`s `bnds`.

# Keyword arguments

`n`: number of draws that are included in animation

`fps`: frames-per-second

`filepath`: if provided, the gif-animation will be stored using this filepath,
otherwise a temporary file will be created.

Further keyword arguments are passed on the `Plots.plot()` function.
"""
function animate_draws(
    uc::T, bnds::Vector{S}; n=20, fps=5, filepath=nothing, label="", kw...
) where {T<:UncertainCurve, S<:UncertainBound}
    # get minimum and maximum y value of all curve draws
    # in the selected `n` samples, used to determine the
    # y-span when plotting
    miny, maxy = [
        mm(curve.y for curve ∈ [NoisySignalIntegration.get_draw(i, uc) for i ∈ 1:n]) |> mm
        for mm in (minimum, maximum)
    ]
    anim = @animate for i in 1:n
        plot(NoisySignalIntegration.get_draw(i, uc), bnds, i; label=label, ylim=(miny*0.9, maxy*1.1), kw...)
    end
    fp = filepath === nothing ? tempname()*".gif" : filepath
    return gif(anim, fp, fps=fps)
end

export animate_draws


text:
using ..Fed: initialize_stats
using GD
using HTTP

"""
    build_router(node::Node)

Build the routes to the node endpoints.
"""
function build_router(node::Node)
    router = HTTP.Router()

    HTTP.@register(
        router,
        "POST",
        node.config.base.fit_node,
        (request::HTTP.Request) -> fit_service(node, request)
    )

    return router
end

"""
    start(node::Node)

Start the given `node` by setting up its router, and registering it to the
server.
"""
function start(node::Node)
    router = build_router(node)

    # hackish way to prevent PayloadSerde to fail while recording stats
    initialize_stats(node.config, 0)

    response = register_to_server(node)
    @assert response.status == 200
    @info "[$(node.host):$(node.port)] $(String(response.body))"

    HTTP.serve(router, node.host, node.port)
end


text:
int_rules_1_2_2_3 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.2.2.3*(d+e*x^2)^q*(a+b*x^2+c*x^4)^p =#
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / (~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ (3 / 4), ~x) => (-2 * (~c * ~d - ~b * ~e) * (~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (1 / 4)) / (~b * ~c * ~x) + (~e / ~c) * Antiderivative((~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (1 / 4) / (~x) ^ 2, ~x) <-- FreeQ([~b, ~c, ~d, ~e], ~x)
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) * (~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => (~e * (~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1)) / (~c * (4 * ~p + 3) * ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~p], ~x) && (Not(IntegerQ(~p)) && (NeQ(4 * ~p + 3, 0) && EqQ(~b * ~e * (2 * ~p + 1) - ~c * ~d * (4 * ~p + 3), 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) * (~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => (~e * (~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1)) / (~c * (4 * ~p + 3) * ~x) - ((~b * ~e * (2 * ~p + 1) - ~c * ~d * (4 * ~p + 3)) / (~c * (4 * ~p + 3))) * Antiderivative((~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~p], ~x) && (Not(IntegerQ(~p)) && (NeQ(4 * ~p + 3, 0) && NeQ(~b * ~e * (2 * ~p + 1) - ~c * ~d * (4 * ~p + 3), 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(q') * (~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => ((~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ FracPart(~p) / ((~x) ^ (2 * FracPart(~p)) * (~b + ~c * (~x) ^ 2) ^ FracPart(~p))) * Antiderivative((~x) ^ (2 * ~p) * (~d + ~e * (~x) ^ 2) ^ ~q * (~b + ~c * (~x) ^ 2) ^ ~p, ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~p, ~q], ~x) && Not(IntegerQ(~p))
    #= Antiderivative(((~d)+(~e')*(~x)^2)^(~q')*((~a)+(~b')*(~x)^2+(~c')*(~x)^4)^(~p'),~x) :=  1/c^p*Antiderivative((d+e*x^2)^q*(b/2+c*x^2)^(2*p),x) <--  FreeQ([a,b,c,d,e,p,q],x) && EqQ(b^2-4*a*c,0) && IntegerQ(p) =#
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(q') * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => ((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p / (~d + ~e * (~x) ^ 2) ^ (2 * ~p)) * Antiderivative((~d + ~e * (~x) ^ 2) ^ (~q + 2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p, ~q], ~x) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (Not(IntegerQ(~p)) && EqQ(2 * ~c * ~d - ~b * ~e, 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(q') * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => ((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ FracPart(~p) / ((~c) ^ IntPart(~p) * (~b / 2 + ~c * (~x) ^ 2) ^ (2 * FracPart(~p)))) * Antiderivative((~d + ~e * (~x) ^ 2) ^ ~q * (~b / 2 + ~c * (~x) ^ 2) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p, ~q], ~x) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~p)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(q') * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => Antiderivative((~d + ~e * (~x) ^ 2) ^ (~p + ~q) * (~a / ~d + (~c / ~e) * (~x) ^ 2) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~q], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (EqQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && IntegerQ(~p)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(q') * (~a + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => Antiderivative((~d + ~e * (~x) ^ 2) ^ (~p + ~q) * (~a / ~d + (~c / ~e) * (~x) ^ 2) ^ ~p, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~q], ~x) && (EqQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && IntegerQ(~p))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => ((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ FracPart(~p) / ((~d + ~e * (~x) ^ 2) ^ FracPart(~p) * (~a / ~d + (~c * (~x) ^ 2) / ~e) ^ FracPart(~p))) * Antiderivative((~d + ~e * (~x) ^ 2) ^ (~p + ~q) * (~a / ~d + (~c / ~e) * (~x) ^ 2) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p, ~q], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (EqQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q * (~a + ~(c') * (~x) ^ 4) ^ ~p, ~x) => ((~a + ~c * (~x) ^ 4) ^ FracPart(~p) / ((~d + ~e * (~x) ^ 2) ^ FracPart(~p) * (~a / ~d + (~c * (~x) ^ 2) / ~e) ^ FracPart(~p))) * Antiderivative((~d + ~e * (~x) ^ 2) ^ (~p + ~q) * (~a / ~d + (~c / ~e) * (~x) ^ 2) ^ ~p, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~p, ~q], ~x) && (EqQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && Not(IntegerQ(~p)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(q') * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ 2) ^ ~q * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (IGtQ(~p, 0) && IGtQ(~q, -2))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(q') * (~a + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ 2) ^ ~q * (~a + ~c * (~x) ^ 4) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (IGtQ(~p, 0) && IGtQ(~q, -2)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => ((~a) ^ ~p * ~x * (~d + ~e * (~x) ^ 2) ^ (~q + 1)) / ~d + (1 / ~d) * Antiderivative((~x) ^ 2 * (~d + ~e * (~x) ^ 2) ^ ~q * (~d * PolynomialQuotient((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p - (~a) ^ ~p, (~x) ^ 2, ~x) - ~e * (~a) ^ ~p * (2 * ~q + 3)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (IGtQ(~p, 0) && (ILtQ(~q + 1 / 2, 0) && LtQ(4 * ~p + 2 * ~q + 1, 0)))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q * (~a + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => ((~a) ^ ~p * ~x * (~d + ~e * (~x) ^ 2) ^ (~q + 1)) / ~d + (1 / ~d) * Antiderivative((~x) ^ 2 * (~d + ~e * (~x) ^ 2) ^ ~q * (~d * PolynomialQuotient((~a + ~c * (~x) ^ 4) ^ ~p - (~a) ^ ~p, (~x) ^ 2, ~x) - ~e * (~a) ^ ~p * (2 * ~q + 3)), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (IGtQ(~p, 0) && (ILtQ(~q + 1 / 2, 0) && LtQ(4 * ~p + 2 * ~q + 1, 0))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => With([Qx = PolynomialQuotient((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~d + ~e * (~x) ^ 2, ~x), R = Coeff(PolynomialRemainder((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~d + ~e * (~x) ^ 2, ~x), ~x, 0)], (-R * ~x * (~d + ~e * (~x) ^ 2) ^ (~q + 1)) / (2 * ~d * (~q + 1)) + (1 / (2 * ~d * (~q + 1))) * Antiderivative((~d + ~e * (~x) ^ 2) ^ (~q + 1) * ExpandToSum(2 * ~d * (~q + 1) * Qx + R * (2 * ~q + 3), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (IGtQ(~p, 0) && LtQ(~q, -1))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q * (~a + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => With([Qx = PolynomialQuotient((~a + ~c * (~x) ^ 4) ^ ~p, ~d + ~e * (~x) ^ 2, ~x), R = Coeff(PolynomialRemainder((~a + ~c * (~x) ^ 4) ^ ~p, ~d + ~e * (~x) ^ 2, ~x), ~x, 0)], (-R * ~x * (~d + ~e * (~x) ^ 2) ^ (~q + 1)) / (2 * ~d * (~q + 1)) + (1 / (2 * ~d * (~q + 1))) * Antiderivative((~d + ~e * (~x) ^ 2) ^ (~q + 1) * ExpandToSum(2 * ~d * (~q + 1) * Qx + R * (2 * ~q + 3), ~x), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (IGtQ(~p, 0) && LtQ(~q, -1)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => ((~c) ^ ~p * (~x) ^ (4 * ~p - 1) * (~d + ~e * (~x) ^ 2) ^ (~q + 1)) / (~e * (4 * ~p + 2 * ~q + 1)) + (1 / (~e * (4 * ~p + 2 * ~q + 1))) * Antiderivative((~d + ~e * (~x) ^ 2) ^ ~q * ExpandToSum((~e * (4 * ~p + 2 * ~q + 1) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p - ~d * (~c) ^ ~p * (4 * ~p - 1) * (~x) ^ (4 * ~p - 2)) - ~e * (~c) ^ ~p * (4 * ~p + 2 * ~q + 1) * (~x) ^ (4 * ~p), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~q], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (IGtQ(~p, 0) && Not(LtQ(~q, -1)))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q * (~a + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => ((~c) ^ ~p * (~x) ^ (4 * ~p - 1) * (~d + ~e * (~x) ^ 2) ^ (~q + 1)) / (~e * (4 * ~p + 2 * ~q + 1)) + (1 / (~e * (4 * ~p + 2 * ~q + 1))) * Antiderivative((~d + ~e * (~x) ^ 2) ^ ~q * ExpandToSum((~e * (4 * ~p + 2 * ~q + 1) * (~a + ~c * (~x) ^ 4) ^ ~p - ~d * (~c) ^ ~p * (4 * ~p - 1) * (~x) ^ (4 * ~p - 2)) - ~e * (~c) ^ ~p * (4 * ~p + 2 * ~q + 1) * (~x) ^ (4 * ~p), ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~q], ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (IGtQ(~p, 0) && Not(LtQ(~q, -1))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((2 * ~d) / ~e - ~b / ~c, 2)], (~e / (2 * ~c)) * Antiderivative(1 / Simp(~d / ~e + q * ~x + (~x) ^ 2, ~x), ~x) + (~e / (2 * ~c)) * Antiderivative(1 / Simp((~d / ~e - q * ~x) + (~x) ^ 2, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (EqQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0) && (GtQ((2 * ~d) / ~e - ~b / ~c, 0) || Not(LtQ((2 * ~d) / ~e - ~b / ~c, 0)) && EqQ(~d - ~e * Rt(~a / ~c, 2), 0))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / (~a + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((2 * ~d) / ~e, 2)], (~e / (2 * ~c)) * Antiderivative(1 / Simp(~d / ~e + q * ~x + (~x) ^ 2, ~x), ~x) + (~e / (2 * ~c)) * Antiderivative(1 / Simp((~d / ~e - q * ~x) + (~x) ^ 2, ~x), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (EqQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0) && PosQ(~d * ~e))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], (~e / 2 + (2 * ~c * ~d - ~b * ~e) / (2q)) * Antiderivative(1 / ((~b / 2 - q / 2) + ~c * (~x) ^ 2), ~x) + (~e / 2 - (2 * ~c * ~d - ~b * ~e) / (2q)) * Antiderivative(1 / (~b / 2 + q / 2 + ~c * (~x) ^ 2), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (EqQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0) && GtQ((~b) ^ 2 - 4 * ~a * ~c, 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((-2 * ~d) / ~e - ~b / ~c, 2)], (~e / (2 * ~c * q)) * Antiderivative((q - 2 * ~x) / Simp((~d / ~e + q * ~x) - (~x) ^ 2, ~x), ~x) + (~e / (2 * ~c * q)) * Antiderivative((q + 2 * ~x) / Simp((~d / ~e - q * ~x) - (~x) ^ 2, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (EqQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0) && Not(GtQ((~b) ^ 2 - 4 * ~a * ~c, 0))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / (~a + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((-2 * ~d) / ~e, 2)], (~e / (2 * ~c * q)) * Antiderivative((q - 2 * ~x) / Simp((~d / ~e + q * ~x) - (~x) ^ 2, ~x), ~x) + (~e / (2 * ~c * q)) * Antiderivative((q + 2 * ~x) / Simp((~d / ~e - q * ~x) - (~x) ^ 2, ~x), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (EqQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0) && NegQ(~d * ~e))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], (~e / 2 + (2 * ~c * ~d - ~b * ~e) / (2q)) * Antiderivative(1 / ((~b / 2 - q / 2) + ~c * (~x) ^ 2), ~x) + (~e / 2 - (2 * ~c * ~d - ~b * ~e) / (2q)) * Antiderivative(1 / (~b / 2 + q / 2 + ~c * (~x) ^ 2), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0) && PosQ((~b) ^ 2 - 4 * ~a * ~c)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / (~a + ~(c') * (~x) ^ 4), ~x) => With([q = Rt(-(~a) * ~c, 2)], (~e / 2 + (~c * ~d) / (2q)) * Antiderivative(1 / (-q + ~c * (~x) ^ 2), ~x) + (~e / 2 - (~c * ~d) / (2q)) * Antiderivative(1 / (q + ~c * (~x) ^ 2), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NeQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0) && PosQ(-(~a) * ~c))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / (~a + ~(c') * (~x) ^ 4), ~x) => With([q = Rt(~a * ~c, 2)], ((~d * q + ~a * ~e) / (2 * ~a * ~c)) * Antiderivative((q + ~c * (~x) ^ 2) / (~a + ~c * (~x) ^ 4), ~x) + ((~d * q - ~a * ~e) / (2 * ~a * ~c)) * Antiderivative((q - ~c * (~x) ^ 2) / (~a + ~c * (~x) ^ 4), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (NeQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0) && NegQ(-(~a) * ~c)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt(~a / ~c, 2)], With([r = Rt(2q - ~b / ~c, 2)], (1 / (2 * ~c * q * r)) * Antiderivative((~d * r - (~d - ~e * q) * ~x) / ((q - r * ~x) + (~x) ^ 2), ~x) + (1 / (2 * ~c * q * r)) * Antiderivative((~d * r + (~d - ~e * q) * ~x) / (q + r * ~x + (~x) ^ 2), ~x))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && NegQ((~b) ^ 2 - 4 * ~a * ~c)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q / (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ 2) ^ ~q / (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && IntegerQ(~q)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q / (~a + ~(c') * (~x) ^ 4), ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ 2) ^ ~q / (~a + ~c * (~x) ^ 4), ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && IntegerQ(~q))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q / (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => ((~e) ^ 2 / ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2)) * Antiderivative((~d + ~e * (~x) ^ 2) ^ ~q, ~x) + (1 / ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2)) * Antiderivative(((~d + ~e * (~x) ^ 2) ^ (~q + 1) * ((~c * ~d - ~b * ~e) - ~c * ~e * (~x) ^ 2)) / (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (Not(IntegerQ(~q)) && LtQ(~q, -1))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q / (~a + ~(c') * (~x) ^ 4), ~x) => ((~e) ^ 2 / (~c * (~d) ^ 2 + ~a * (~e) ^ 2)) * Antiderivative((~d + ~e * (~x) ^ 2) ^ ~q, ~x) + (~c / (~c * (~d) ^ 2 + ~a * (~e) ^ 2)) * Antiderivative(((~d + ~e * (~x) ^ 2) ^ (~q + 1) * (~d - ~e * (~x) ^ 2)) / (~a + ~c * (~x) ^ 4), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (Not(IntegerQ(~q)) && LtQ(~q, -1)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q / (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([r = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((2 * ~c) / r) * Antiderivative((~d + ~e * (~x) ^ 2) ^ ~q / ((~b - r) + 2 * ~c * (~x) ^ 2), ~x) - ((2 * ~c) / r) * Antiderivative((~d + ~e * (~x) ^ 2) ^ ~q / (~b + r + 2 * ~c * (~x) ^ 2), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~q], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && Not(IntegerQ(~q))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q / (~a + ~(c') * (~x) ^ 4), ~x) => With([r = Rt(-(~a) * ~c, 2)], (-(~c) / (2r)) * Antiderivative((~d + ~e * (~x) ^ 2) ^ ~q / (r - ~c * (~x) ^ 2), ~x) - (~c / (2r)) * Antiderivative((~d + ~e * (~x) ^ 2) ^ ~q / (r + ~c * (~x) ^ 2), ~x)) <-- FreeQ([~a, ~c, ~d, ~e, ~q], ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && Not(IntegerQ(~q)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => (~x * (2 * ~b * ~e * ~p + ~c * ~d * (4 * ~p + 3) + ~c * ~e * (4 * ~p + 1) * (~x) ^ 2) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p) / (~c * (4 * ~p + 1) * (4 * ~p + 3)) + ((2 * ~p) / (~c * (4 * ~p + 1) * (4 * ~p + 3))) * Antiderivative(Simp((2 * ~a * ~c * ~d * (4 * ~p + 3) - ~a * ~b * ~e) + ((2 * ~a * ~c * ~e * (4 * ~p + 1) + ~b * ~c * ~d * (4 * ~p + 3)) - (~b) ^ 2 * ~e * (2 * ~p + 1)) * (~x) ^ 2, ~x) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (GtQ(~p, 0) && (FractionQ(~p) && IntegerQ(2 * ~p)))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) * (~a + ~(c') * (~x) ^ 4) ^ ~p, ~x) => (~x * (~d * (4 * ~p + 3) + ~e * (4 * ~p + 1) * (~x) ^ 2) * (~a + ~c * (~x) ^ 4) ^ ~p) / ((4 * ~p + 1) * (4 * ~p + 3)) + ((2 * ~p) / ((4 * ~p + 1) * (4 * ~p + 3))) * Antiderivative(Simp(2 * ~a * ~d * (4 * ~p + 3) + (2 * ~a * ~e * (4 * ~p + 1)) * (~x) ^ 2, ~x) * (~a + ~c * (~x) ^ 4) ^ (~p - 1), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (GtQ(~p, 0) && (FractionQ(~p) && IntegerQ(2 * ~p))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => (~x * ((~a * ~b * ~e - ~d * ((~b) ^ 2 - 2 * ~a * ~c)) - ~c * (~b * ~d - 2 * ~a * ~e) * (~x) ^ 2) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1)) / (2 * ~a * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c)) + (1 / (2 * ~a * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c))) * Antiderivative(Simp((((2 * ~p + 3) * ~d * (~b) ^ 2 - ~a * ~b * ~e) - 2 * ~a * ~c * ~d * (4 * ~p + 5)) + (4 * ~p + 7) * (~d * ~b - 2 * ~a * ~e) * ~c * (~x) ^ 2, ~x) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (LtQ(~p, -1) && IntegerQ(2 * ~p))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) * (~a + ~(c') * (~x) ^ 4) ^ ~p, ~x) => (-(~x) * (~d + ~e * (~x) ^ 2) * (~a + ~c * (~x) ^ 4) ^ (~p + 1)) / (4 * ~a * (~p + 1)) + (1 / (4 * ~a * (~p + 1))) * Antiderivative(Simp(~d * (4 * ~p + 5) + ~e * (4 * ~p + 7) * (~x) ^ 2, ~x) * (~a + ~c * (~x) ^ 4) ^ (~p + 1), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (LtQ(~p, -1) && IntegerQ(2 * ~p)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], 2 * sqrt(-(~c)) * Antiderivative((~d + ~e * (~x) ^ 2) / (sqrt(~b + q + 2 * ~c * (~x) ^ 2) * sqrt((-(~b) + q) - 2 * ~c * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (GtQ((~b) ^ 2 - 4 * ~a * ~c, 0) && LtQ(~c, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(c') * (~x) ^ 4), ~x) => With([q = Rt(-(~a) * ~c, 2)], sqrt(-(~c)) * Antiderivative((~d + ~e * (~x) ^ 2) / (sqrt(q + ~c * (~x) ^ 2) * sqrt(q - ~c * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (GtQ(~a, 0) && LtQ(~c, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt(~c / ~a, 4)], (-(~d) * ~x * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)) / (~a * (1 + q ^ 2 * (~x) ^ 2)) + ((~d * (1 + q ^ 2 * (~x) ^ 2) * sqrt((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) / (~a * (1 + q ^ 2 * (~x) ^ 2) ^ 2))) / (q * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4))) * Elliptic.E(2 * atan(q * ~x), 1 / 2 - (~b * q ^ 2) / (4 * ~c)) <-- EqQ(~e + ~d * q ^ 2, 0)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (GtQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (GtQ(~c / ~a, 0) && LtQ(~b / ~a, 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt(~c / ~a, 2)], ((~e + ~d * q) / q) * Antiderivative(1 / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x) - (~e / q) * Antiderivative((1 - q * (~x) ^ 2) / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x) <-- NeQ(~e + ~d * q, 0)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (GtQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (GtQ(~c / ~a, 0) && LtQ(~b / ~a, 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], (~e * ~x * (~b + q + 2 * ~c * (~x) ^ 2)) / (2 * ~c * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)) - ((~e * q * sqrt((2 * ~a + (~b - q) * (~x) ^ 2) / (2 * ~a + (~b + q) * (~x) ^ 2)) * sqrt((2 * ~a + (~b + q) * (~x) ^ 2) / q)) / (2 * ~c * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) * sqrt(~a / (2 * ~a + (~b + q) * (~x) ^ 2)))) * Elliptic.E(asin(~x / sqrt((2 * ~a + (~b + q) * (~x) ^ 2) / (2q))), (~b + q) / (2q)) <-- EqQ(2 * ~c * ~d - ~e * (~b - q), 0)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (GtQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (LtQ(~a, 0) && GtQ(~c, 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(c') * (~x) ^ 4), ~x) => With([q = Rt(-(~a) * ~c, 2)], (~e * ~x * (q + ~c * (~x) ^ 2)) / (~c * sqrt(~a + ~c * (~x) ^ 4)) - ((sqrt(2) * ~e * q * sqrt(-(~a) + q * (~x) ^ 2) * sqrt((~a + q * (~x) ^ 2) / q)) / (sqrt(-(~a)) * ~c * sqrt(~a + ~c * (~x) ^ 4))) * Elliptic.E(asin(~x / sqrt((~a + q * (~x) ^ 2) / (2q))), 1 / 2) <-- EqQ(~c * ~d + ~e * q, 0) && IntegerQ(q)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (LtQ(~a, 0) && GtQ(~c, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(c') * (~x) ^ 4), ~x) => With([q = Rt(-(~a) * ~c, 2)], (~e * ~x * (q + ~c * (~x) ^ 2)) / (~c * sqrt(~a + ~c * (~x) ^ 4)) - ((sqrt(2) * ~e * q * sqrt((~a - q * (~x) ^ 2) / (~a + q * (~x) ^ 2)) * sqrt((~a + q * (~x) ^ 2) / q)) / (~c * sqrt(~a + ~c * (~x) ^ 4) * sqrt(~a / (~a + q * (~x) ^ 2)))) * Elliptic.E(asin(~x / sqrt((~a + q * (~x) ^ 2) / (2q))), 1 / 2) <-- EqQ(~c * ~d + ~e * q, 0)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (LtQ(~a, 0) && GtQ(~c, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((2 * ~c * ~d - ~e * (~b - q)) / (2 * ~c)) * Antiderivative(1 / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x) + (~e / (2 * ~c)) * Antiderivative(((~b - q) + 2 * ~c * (~x) ^ 2) / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x) <-- NeQ(2 * ~c * ~d - ~e * (~b - q), 0)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (GtQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (LtQ(~a, 0) && GtQ(~c, 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(c') * (~x) ^ 4), ~x) => With([q = Rt(-(~a) * ~c, 2)], ((~c * ~d + ~e * q) / ~c) * Antiderivative(1 / sqrt(~a + ~c * (~x) ^ 4), ~x) - (~e / ~c) * Antiderivative((q - ~c * (~x) ^ 2) / sqrt(~a + ~c * (~x) ^ 4), ~x) <-- NeQ(~c * ~d + ~e * q, 0)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (LtQ(~a, 0) && GtQ(~c, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ~d * Antiderivative(1 / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x) + ~e * Antiderivative((~x) ^ 2 / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x) <-- PosQ((~b + q) / ~a) || PosQ((~b - q) / ~a)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && GtQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(c') * (~x) ^ 4), ~x) => ~d * Antiderivative(1 / sqrt(~a + ~c * (~x) ^ 4), ~x) + ~e * Antiderivative((~x) ^ 2 / sqrt(~a + ~c * (~x) ^ 4), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && GtQ(-(~a) * ~c, 0)
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((-(~a) * ~e * Rt(-((~b + q)) / (2 * ~a), 2) * sqrt(1 + ((~b + q) * (~x) ^ 2) / (2 * ~a)) * sqrt(1 + ((~b - q) * (~x) ^ 2) / (2 * ~a))) / (~c * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4))) * Elliptic.E(asin(Rt(-((~b + q)) / (2 * ~a), 2) * ~x), (~b - q) / (~b + q)) <-- NegQ((~b + q) / ~a) && (EqQ(2 * ~c * ~d - ~e * (~b + q), 0) && Not(SimplerSqrtQ(-((~b - q)) / (2 * ~a), -((~b + q)) / (2 * ~a))))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && GtQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((2 * ~c * ~d - ~e * (~b + q)) / (2 * ~c)) * Antiderivative(1 / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x) + (~e / (2 * ~c)) * Antiderivative((~b + q + 2 * ~c * (~x) ^ 2) / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x) <-- NegQ((~b + q) / ~a) && (NeQ(2 * ~c * ~d - ~e * (~b + q), 0) && Not(SimplerSqrtQ(-((~b - q)) / (2 * ~a), -((~b + q)) / (2 * ~a))))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && GtQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((-(~a) * ~e * Rt(-((~b - q)) / (2 * ~a), 2) * sqrt(1 + ((~b - q) * (~x) ^ 2) / (2 * ~a)) * sqrt(1 + ((~b + q) * (~x) ^ 2) / (2 * ~a))) / (~c * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4))) * Elliptic.E(asin(Rt(-((~b - q)) / (2 * ~a), 2) * ~x), (~b + q) / (~b - q)) <-- NegQ((~b - q) / ~a) && EqQ(2 * ~c * ~d - ~e * (~b - q), 0)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && GtQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((2 * ~c * ~d - ~e * (~b - q)) / (2 * ~c)) * Antiderivative(1 / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x) + (~e / (2 * ~c)) * Antiderivative(((~b - q) + 2 * ~c * (~x) ^ 2) / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x) <-- NegQ((~b - q) / ~a) && NeQ(2 * ~c * ~d - ~e * (~b - q), 0)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && GtQ((~b) ^ 2 - 4 * ~a * ~c, 0)
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt(~c / ~a, 4)], (-(~d) * ~x * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)) / (~a * (1 + q ^ 2 * (~x) ^ 2)) + ((~d * (1 + q ^ 2 * (~x) ^ 2) * sqrt((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) / (~a * (1 + q ^ 2 * (~x) ^ 2) ^ 2))) / (q * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4))) * Elliptic.E(2 * atan(q * ~x), 1 / 2 - (~b * q ^ 2) / (4 * ~c)) <-- EqQ(~e + ~d * q ^ 2, 0)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && PosQ(~c / ~a))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(c') * (~x) ^ 4), ~x) => With([q = Rt(~c / ~a, 4)], (-(~d) * ~x * sqrt(~a + ~c * (~x) ^ 4)) / (~a * (1 + q ^ 2 * (~x) ^ 2)) + ((~d * (1 + q ^ 2 * (~x) ^ 2) * sqrt((~a + ~c * (~x) ^ 4) / (~a * (1 + q ^ 2 * (~x) ^ 2) ^ 2))) / (q * sqrt(~a + ~c * (~x) ^ 4))) * Elliptic.E(2 * atan(q * ~x), 1 / 2) <-- EqQ(~e + ~d * q ^ 2, 0)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && PosQ(~c / ~a)
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt(~c / ~a, 2)], ((~e + ~d * q) / q) * Antiderivative(1 / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x) - (~e / q) * Antiderivative((1 - q * (~x) ^ 2) / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x) <-- NeQ(~e + ~d * q, 0)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && PosQ(~c / ~a))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(c') * (~x) ^ 4), ~x) => With([q = Rt(~c / ~a, 2)], ((~e + ~d * q) / q) * Antiderivative(1 / sqrt(~a + ~c * (~x) ^ 4), ~x) - (~e / q) * Antiderivative((1 - q * (~x) ^ 2) / sqrt(~a + ~c * (~x) ^ 4), ~x) <-- NeQ(~e + ~d * q, 0)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && PosQ(~c / ~a)
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(c') * (~x) ^ 4), ~x) => (~d / sqrt(~a)) * Antiderivative(sqrt(1 + (~e * (~x) ^ 2) / ~d) / sqrt(1 - (~e * (~x) ^ 2) / ~d), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NegQ(~c / ~a) && (EqQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && GtQ(~a, 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(c') * (~x) ^ 4), ~x) => (sqrt(1 + (~c * (~x) ^ 4) / ~a) / sqrt(~a + ~c * (~x) ^ 4)) * Antiderivative((~d + ~e * (~x) ^ 2) / sqrt(1 + (~c * (~x) ^ 4) / ~a), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NegQ(~c / ~a) && (EqQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && Not(GtQ(~a, 0))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(c') * (~x) ^ 4), ~x) => With([q = Rt(-(~c) / ~a, 2)], ((~d * q - ~e) / q) * Antiderivative(1 / sqrt(~a + ~c * (~x) ^ 4), ~x) + (~e / q) * Antiderivative((1 + q * (~x) ^ 2) / sqrt(~a + ~c * (~x) ^ 4), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NegQ(~c / ~a) && NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((sqrt(1 + (2 * ~c * (~x) ^ 2) / (~b - q)) * sqrt(1 + (2 * ~c * (~x) ^ 2) / (~b + q))) / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)) * Antiderivative((~d + ~e * (~x) ^ 2) / (sqrt(1 + (2 * ~c * (~x) ^ 2) / (~b - q)) * sqrt(1 + (2 * ~c * (~x) ^ 2) / (~b + q))), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && NegQ(~c / ~a))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ 2) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) * (~a + ~(c') * (~x) ^ 4) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ 2) * (~a + ~c * (~x) ^ 4) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0)
    #= Antiderivative(((~d)+(~e')*(~x)^2)^2/sqrt((~a)+(~b')*(~x)^2+(~c')*(~x)^4),~x) :=  e^2*x*sqrt(a+b*x^2+c*x^4)/(3*c) + 2*(3*c*d-b*e)/(3*c)*Antiderivative((d+e*x^2)/sqrt(a+b*x^2+c*x^4),x) - (3*c*d^2-2*b*d*e+a*e^2)/(3*c)*Antiderivative(1/sqrt(a+b*x^2+c*x^4),x) <--  FreeQ([a,b,c,d,e],x) && NeQ(b^2-4*a*c,0) && NeQ(c*d^2-b*d*e+a*e^2,0) =#
    #= Antiderivative(((~d)+(~e')*(~x)^2)^2/sqrt((~a)+(~c')*(~x)^4),~x) :=  e^2*x*sqrt(a+c*x^4)/(3*c) + 2*d*Antiderivative((d+e*x^2)/sqrt(a+c*x^4),x) - (3*c*d^2+a*e^2)/(3*c)*Antiderivative(1/sqrt(a+c*x^4),x) <-- FreeQ([a,c,d,e],x) &&  NeQ(c*d^2+a*e^2,0) =#
    #= Antiderivative(((~d)+(~e')*(~x)^2)^(~q)/sqrt((~a)+(~b')*(~x)^2+(~c')*(~x)^4),~x) :=  e^2*x*(d+e*x^2)^(q-2)*sqrt(a+b*x^2+c*x^4)/(c*(2*q-1)) + 2*(q-1)*(3*c*d-b*e)/(c*(2*q-1))*Antiderivative((d+e*x^2)^(q-1)/sqrt(a+b*x^2+c* x^4),x) - (2*q-3)*(3*c*d^2-2*b*d*e+a*e^2)/(c*(2*q-1))*Antiderivative((d+e*x^2)^(q-2)/ sqrt(a+b*x^2+c*x^4),x) + 2*d*(q-2)*(c*d^2-b*d*e+a*e^2)/(c*(2*q-1))*Antiderivative((d+e*x^2)^(q-3)/sqrt( a+b*x^2+c*x^4),x) <-- FreeQ([a,b,c,d,e],x) && NeQ(b^2-4*a*c,0) &&  IGtQ(q,2) =#
    #= Antiderivative(((~d)+(~e')*(~x)^2)^(~q)/sqrt((~a)+(~c')*(~x)^4),~x) :=  e^2*x*(d+e*x^2)^(q-2)*sqrt(a+c*x^4)/(c*(2*q-1)) + 6*d*(q-1)/(2*q-1)*Antiderivative((d+e*x^2)^(q-1)/sqrt(a+c*x^4),x) - (2*q-3)*(3*c*d^2+a*e^2)/(c*(2*q-1))*Antiderivative((d+e*x^2)^(q-2)/sqrt(a+c*x^ 4),x) + 2*d*(q-2)*(c*d^2+a*e^2)/(c*(2*q-1))*Antiderivative((d+e*x^2)^(q-3)/sqrt(a+c*x^ 4),x) <-- FreeQ([a,c,d,e],x) && IGtQ(q,2) =#
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => With([f = Coeff(PolynomialRemainder((~d + ~e * (~x) ^ 2) ^ ~q, ~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4, ~x), ~x, 0), g = Coeff(PolynomialRemainder((~d + ~e * (~x) ^ 2) ^ ~q, ~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4, ~x), ~x, 2)], (~x * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1) * ((~a * ~b * g - f * ((~b) ^ 2 - 2 * ~a * ~c)) - ~c * (~b * f - 2 * ~a * g) * (~x) ^ 2)) / (2 * ~a * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c)) + (1 / (2 * ~a * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c))) * Antiderivative((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1) * ExpandToSum((((2 * ~a * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c) * PolynomialQuotient((~d + ~e * (~x) ^ 2) ^ ~q, ~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4, ~x) + (~b) ^ 2 * f * (2 * ~p + 3)) - 2 * ~a * ~c * f * (4 * ~p + 5)) - ~a * ~b * g) + ~c * (4 * ~p + 7) * (~b * f - 2 * ~a * g) * (~x) ^ 2, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (IGtQ(~q, 1) && LtQ(~p, -1))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => ((~e) ^ ~q * (~x) ^ (2 * ~q - 3) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1)) / (~c * (4 * ~p + 2 * ~q + 1)) + (1 / (~c * (4 * ~p + 2 * ~q + 1))) * Antiderivative((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p * ExpandToSum(((~c * (4 * ~p + 2 * ~q + 1) * (~d + ~e * (~x) ^ 2) ^ ~q - ~a * (2 * ~q - 3) * (~e) ^ ~q * (~x) ^ (2 * ~q - 4)) - ~b * ((2 * ~p + 2 * ~q) - 1) * (~e) ^ ~q * (~x) ^ (2 * ~q - 2)) - ~c * (4 * ~p + 2 * ~q + 1) * (~e) ^ ~q * (~x) ^ (2 * ~q), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && IGtQ(~q, 1)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q * (~a + ~(c') * (~x) ^ 4) ^ ~p, ~x) => ((~e) ^ ~q * (~x) ^ (2 * ~q - 3) * (~a + ~c * (~x) ^ 4) ^ (~p + 1)) / (~c * (4 * ~p + 2 * ~q + 1)) + (1 / (~c * (4 * ~p + 2 * ~q + 1))) * Antiderivative((~a + ~c * (~x) ^ 4) ^ ~p * ExpandToSum((~c * (4 * ~p + 2 * ~q + 1) * (~d + ~e * (~x) ^ 2) ^ ~q - ~a * (2 * ~q - 3) * (~e) ^ ~q * (~x) ^ (2 * ~q - 4)) - ~c * (4 * ~p + 2 * ~q + 1) * (~e) ^ ~q * (~x) ^ (2 * ~q), ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~p], ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && IGtQ(~q, 1))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p / (~d + ~(e') * (~x) ^ 2), ~x) => (-1 / (~e) ^ 2) * Antiderivative(((~c * ~d - ~b * ~e) - ~c * ~e * (~x) ^ 2) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p - 1), ~x) + (((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2) / (~e) ^ 2) * Antiderivative((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p - 1) / (~d + ~e * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && IGtQ(~p + 1 / 2, 0)))
    @apply_utils Antiderivative((~a + ~(c') * (~x) ^ 4) ^ ~p / (~d + ~(e') * (~x) ^ 2), ~x) => (-1 / (~e) ^ 2) * Antiderivative((~c * ~d - ~c * ~e * (~x) ^ 2) * (~a + ~c * (~x) ^ 4) ^ (~p - 1), ~x) + ((~c * (~d) ^ 2 + ~a * (~e) ^ 2) / (~e) ^ 2) * Antiderivative((~a + ~c * (~x) ^ 4) ^ (~p - 1) / (~d + ~e * (~x) ^ 2), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && IGtQ(~p + 1 / 2, 0))
    @apply_utils Antiderivative(1 / ((~d + ~(e') * (~x) ^ 2) * sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4)), ~x) => (1 / (2 * ~d)) * Antiderivative(1 / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x) + (1 / (2 * ~d)) * Antiderivative((~d - ~e * (~x) ^ 2) / ((~d + ~e * (~x) ^ 2) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && EqQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0)))
    @apply_utils Antiderivative(1 / ((~d + ~(e') * (~x) ^ 2) * sqrt(~a + ~(c') * (~x) ^ 4)), ~x) => (1 / (2 * ~d)) * Antiderivative(1 / sqrt(~a + ~c * (~x) ^ 4), ~x) + (1 / (2 * ~d)) * Antiderivative((~d - ~e * (~x) ^ 2) / ((~d + ~e * (~x) ^ 2) * sqrt(~a + ~c * (~x) ^ 4)), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && EqQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0))
    @apply_utils Antiderivative(1 / ((~d + ~(e') * (~x) ^ 2) * sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4)), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], 2 * sqrt(-(~c)) * Antiderivative(1 / ((~d + ~e * (~x) ^ 2) * sqrt(~b + q + 2 * ~c * (~x) ^ 2) * sqrt((-(~b) + q) - 2 * ~c * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (GtQ((~b) ^ 2 - 4 * ~a * ~c, 0) && LtQ(~c, 0))
    @apply_utils Antiderivative(1 / ((~d + ~(e') * (~x) ^ 2) * sqrt(~a + ~(c') * (~x) ^ 4)), ~x) => With([q = Rt(-(~a) * ~c, 2)], sqrt(-(~c)) * Antiderivative(1 / ((~d + ~e * (~x) ^ 2) * sqrt(q + ~c * (~x) ^ 2) * sqrt(q - ~c * (~x) ^ 2)), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (GtQ(~a, 0) && LtQ(~c, 0))
    @apply_utils Antiderivative(1 / ((~d + ~(e') * (~x) ^ 2) * sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4)), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((2 * ~c) / (2 * ~c * ~d - ~e * (~b - q))) * Antiderivative(1 / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x) - (~e / (2 * ~c * ~d - ~e * (~b - q))) * Antiderivative(((~b - q) + 2 * ~c * (~x) ^ 2) / ((~d + ~e * (~x) ^ 2) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (GtQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(LtQ(~c, 0)))
    @apply_utils Antiderivative(1 / ((~d + ~(e') * (~x) ^ 2) * sqrt(~a + ~(c') * (~x) ^ 4)), ~x) => With([q = Rt(-(~a) * ~c, 2)], (~c / (~c * ~d + ~e * q)) * Antiderivative(1 / sqrt(~a + ~c * (~x) ^ 4), ~x) + (~e / (~c * ~d + ~e * q)) * Antiderivative((q - ~c * (~x) ^ 2) / ((~d + ~e * (~x) ^ 2) * sqrt(~a + ~c * (~x) ^ 4)), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (GtQ(-(~a) * ~c, 0) && Not(LtQ(~c, 0)))
    @apply_utils Antiderivative(1 / ((~d + ~(e') * (~x) ^ 2) * sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4)), ~x) => With([q = Rt(~c / ~a, 2)], ((~c * ~d + ~a * ~e * q) / (~c * (~d) ^ 2 - ~a * (~e) ^ 2)) * Antiderivative(1 / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x) - ((~a * ~e * (~e + ~d * q)) / (~c * (~d) ^ 2 - ~a * (~e) ^ 2)) * Antiderivative((1 + q * (~x) ^ 2) / ((~d + ~e * (~x) ^ 2) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (NeQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0) && PosQ(~c / ~a))))
    @apply_utils Antiderivative(1 / ((~d + ~(e') * (~x) ^ 2) * sqrt(~a + ~(c') * (~x) ^ 4)), ~x) => With([q = Rt(~c / ~a, 2)], ((~c * ~d + ~a * ~e * q) / (~c * (~d) ^ 2 - ~a * (~e) ^ 2)) * Antiderivative(1 / sqrt(~a + ~c * (~x) ^ 4), ~x) - ((~a * ~e * (~e + ~d * q)) / (~c * (~d) ^ 2 - ~a * (~e) ^ 2)) * Antiderivative((1 + q * (~x) ^ 2) / ((~d + ~e * (~x) ^ 2) * sqrt(~a + ~c * (~x) ^ 4)), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (NeQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0) && PosQ(~c / ~a)))
    @apply_utils Antiderivative(1 / ((~d + ~(e') * (~x) ^ 2) * sqrt(~a + ~(c') * (~x) ^ 4)), ~x) => With([q = Rt(-(~c) / ~a, 4)], (1 / (~d * sqrt(~a) * q)) * Elliptic.Pi(-(~e) / (~d * q ^ 2), asin(q * ~x), -1)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NegQ(~c / ~a) && GtQ(~a, 0))
    @apply_utils Antiderivative(1 / ((~d + ~(e') * (~x) ^ 2) * sqrt(~a + ~(c') * (~x) ^ 4)), ~x) => (sqrt(1 + (~c * (~x) ^ 4) / ~a) / sqrt(~a + ~c * (~x) ^ 4)) * Antiderivative(1 / ((~d + ~e * (~x) ^ 2) * sqrt(1 + (~c * (~x) ^ 4) / ~a)), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NegQ(~c / ~a) && Not(GtQ(~a, 0)))
    @apply_utils Antiderivative(1 / ((~d + ~(e') * (~x) ^ 2) * sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4)), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((sqrt(1 + (2 * ~c * (~x) ^ 2) / (~b - q)) * sqrt(1 + (2 * ~c * (~x) ^ 2) / (~b + q))) / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)) * Antiderivative(1 / ((~d + ~e * (~x) ^ 2) * sqrt(1 + (2 * ~c * (~x) ^ 2) / (~b - q)) * sqrt(1 + (2 * ~c * (~x) ^ 2) / (~b + q))), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && NegQ(~c / ~a))
    @apply_utils Antiderivative((~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p / (~d + ~(e') * (~x) ^ 2), ~x) => (1 / ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2)) * Antiderivative(((~c * ~d - ~b * ~e) - ~c * ~e * (~x) ^ 2) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~x) + ((~e) ^ 2 / ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2)) * Antiderivative((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1) / (~d + ~e * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && ILtQ(~p + 1 / 2, 0)))
    @apply_utils Antiderivative((~a + ~(c') * (~x) ^ 4) ^ ~p / (~d + ~(e') * (~x) ^ 2), ~x) => (1 / (~c * (~d) ^ 2 + ~a * (~e) ^ 2)) * Antiderivative((~c * ~d - ~c * ~e * (~x) ^ 2) * (~a + ~c * (~x) ^ 4) ^ ~p, ~x) + ((~e) ^ 2 / (~c * (~d) ^ 2 + ~a * (~e) ^ 2)) * Antiderivative((~a + ~c * (~x) ^ 4) ^ (~p + 1) / (~d + ~e * (~x) ^ 2), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && ILtQ(~p + 1 / 2, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q / sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4), ~x) => (-((~e) ^ 2) * ~x * (~d + ~e * (~x) ^ 2) ^ (~q + 1) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)) / (2 * ~d * (~q + 1) * ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2)) + (1 / (2 * ~d * (~q + 1) * ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2))) * Antiderivative(((~d + ~e * (~x) ^ 2) ^ (~q + 1) / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)) * Simp(((~a * (~e) ^ 2 * (2 * ~q + 3) + 2 * ~d * (~c * ~d - ~b * ~e) * (~q + 1)) - 2 * ~e * (~c * ~d * (~q + 1) - ~b * ~e * (~q + 2)) * (~x) ^ 2) + ~c * (~e) ^ 2 * (2 * ~q + 5) * (~x) ^ 4, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && ILtQ(~q, -1))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q / sqrt(~a + ~(c') * (~x) ^ 4), ~x) => (-((~e) ^ 2) * ~x * (~d + ~e * (~x) ^ 2) ^ (~q + 1) * sqrt(~a + ~c * (~x) ^ 4)) / (2 * ~d * (~q + 1) * (~c * (~d) ^ 2 + ~a * (~e) ^ 2)) + (1 / (2 * ~d * (~q + 1) * (~c * (~d) ^ 2 + ~a * (~e) ^ 2))) * Antiderivative(((~d + ~e * (~x) ^ 2) ^ (~q + 1) / sqrt(~a + ~c * (~x) ^ 4)) * Simp(((~a * (~e) ^ 2 * (2 * ~q + 3) + 2 * ~c * (~d) ^ 2 * (~q + 1)) - 2 * ~e * ~c * ~d * (~q + 1) * (~x) ^ 2) + ~c * (~e) ^ 2 * (2 * ~q + 5) * (~x) ^ 4, ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && ILtQ(~q, -1)
    @apply_utils Antiderivative(sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) / (~d + ~(e') * (~x) ^ 2) ^ 2, ~x) => With([q = Rt(~e / ~d, 2)], ((~c * (~d + ~e * (~x) ^ 2) * sqrt(((~e) ^ 2 * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)) / (~c * (~d + ~e * (~x) ^ 2) ^ 2))) / (2 * ~d * (~e) ^ 2 * q * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4))) * Elliptic.E(2 * atan(q * ~x), (2 * ~c * ~d - ~b * ~e) / (4 * ~c * ~d))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (EqQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0) && PosQ(~e / ~d))))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) / (~d + ~(e') * (~x) ^ 2) ^ 2, ~x) => ((~x * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)) / (2 * ~d * (~d + ~e * (~x) ^ 2)) + (~c / (2 * ~d * (~e) ^ 2)) * Antiderivative((~d - ~e * (~x) ^ 2) / sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4), ~x)) - ((~c * (~d) ^ 2 - ~a * (~e) ^ 2) / (2 * ~d * (~e) ^ 2)) * Antiderivative(1 / ((~d + ~e * (~x) ^ 2) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0))
    @apply_utils Antiderivative(sqrt(~a + ~(c') * (~x) ^ 4) / (~d + ~(e') * (~x) ^ 2) ^ 2, ~x) => ((~x * sqrt(~a + ~c * (~x) ^ 4)) / (2 * ~d * (~d + ~e * (~x) ^ 2)) + (~c / (2 * ~d * (~e) ^ 2)) * Antiderivative((~d - ~e * (~x) ^ 2) / sqrt(~a + ~c * (~x) ^ 4), ~x)) - ((~c * (~d) ^ 2 - ~a * (~e) ^ 2) / (2 * ~d * (~e) ^ 2)) * Antiderivative(1 / ((~d + ~e * (~x) ^ 2) * sqrt(~a + ~c * (~x) ^ 4)), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0)
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => Module([aa, bb, cc], Antiderivative(ReplaceAll(ExpandIntegrand(1 / sqrt(aa + bb * (~x) ^ 2 + cc * (~x) ^ 4), (~d + ~e * (~x) ^ 2) ^ ~q * (aa + bb * (~x) ^ 2 + cc * (~x) ^ 4) ^ (~p + 1 / 2), ~x), [(aa->begin
                                    #= none:1 =#
                                    #= none:3 =#
                                    ~a
                                end), (bb->begin
                                    #= none:4 =#
                                    #= none:6 =#
                                    ~b
                                end), (cc->begin
                                    #= none:7 =#
                                    #= none:9 =#
                                    ~c
                                end)]), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (ILtQ(~q, 0) && IntegerQ(~p + 1 / 2))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q * (~a + ~(c') * (~x) ^ 4) ^ ~p, ~x) => Module([aa, cc], Antiderivative(ReplaceAll(ExpandIntegrand(1 / sqrt(aa + cc * (~x) ^ 4), (~d + ~e * (~x) ^ 2) ^ ~q * (aa + cc * (~x) ^ 4) ^ (~p + 1 / 2), ~x), [(aa->begin
                                    #= none:1 =#
                                    #= none:3 =#
                                    ~a
                                end), (cc->begin
                                    #= none:4 =#
                                    #= none:6 =#
                                    ~c
                                end)]), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (ILtQ(~q, 0) && IntegerQ(~p + 1 / 2)))
    @apply_utils Antiderivative(1 / (sqrt(~d + ~(e') * (~x) ^ 2) * sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4)), ~x) => (1 / (2 * sqrt(~a) * sqrt(~d) * Rt(-(~e) / ~d, 2))) * Elliptic.F(2 * asin(Rt(-(~e) / ~d, 2) * ~x), (~b * ~d) / (4 * ~a * ~e)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(~c * ~d - ~b * ~e, 0) && (GtQ(~a, 0) && GtQ(~d, 0)))
    @apply_utils Antiderivative(1 / (sqrt(~d + ~(e') * (~x) ^ 2) * sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4)), ~x) => ((sqrt((~d + ~e * (~x) ^ 2) / ~d) * sqrt((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) / ~a)) / (sqrt(~d + ~e * (~x) ^ 2) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4))) * Antiderivative(1 / (sqrt(1 + (~e / ~d) * (~x) ^ 2) * sqrt(1 + (~b / ~a) * (~x) ^ 2 + (~c / ~a) * (~x) ^ 4)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(~c * ~d - ~b * ~e, 0) && Not(GtQ(~a, 0) && GtQ(~d, 0)))
    @apply_utils Antiderivative(1 / (sqrt(~d + ~(e') * (~x) ^ 2) * sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4)), ~x) => (((~x) ^ 3 * sqrt(~e + ~d / (~x) ^ 2) * sqrt(~c + ~b / (~x) ^ 2 + ~a / (~x) ^ 4)) / (sqrt(~d + ~e * (~x) ^ 2) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4))) * Antiderivative(1 / ((~x) ^ 3 * sqrt(~e + ~d / (~x) ^ 2) * sqrt(~c + ~b / (~x) ^ 2 + ~a / (~x) ^ 4)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0))
    @apply_utils Antiderivative(1 / (sqrt(~d + ~(e') * (~x) ^ 2) * sqrt(~a + ~(c') * (~x) ^ 4)), ~x) => (((~x) ^ 3 * sqrt(~e + ~d / (~x) ^ 2) * sqrt(~c + ~a / (~x) ^ 4)) / (sqrt(~d + ~e * (~x) ^ 2) * sqrt(~a + ~c * (~x) ^ 4))) * Antiderivative(1 / ((~x) ^ 3 * sqrt(~e + ~d / (~x) ^ 2) * sqrt(~c + ~a / (~x) ^ 4)), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0)
    @apply_utils Antiderivative(sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) / sqrt(~d + ~(e') * (~x) ^ 2), ~x) => (sqrt(~a) / (2 * sqrt(~d) * Rt(-(~e) / ~d, 2))) * Elliptic.E(2 * asin(Rt(-(~e) / ~d, 2) * ~x), (~b * ~d) / (4 * ~a * ~e)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(~c * ~d - ~b * ~e, 0) && (GtQ(~a, 0) && GtQ(~d, 0)))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) / sqrt(~d + ~(e') * (~x) ^ 2), ~x) => ((sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) * sqrt((~d + ~e * (~x) ^ 2) / ~d)) / (sqrt(~d + ~e * (~x) ^ 2) * sqrt((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) / ~a))) * Antiderivative(sqrt(1 + (~b / ~a) * (~x) ^ 2 + (~c / ~a) * (~x) ^ 4) / sqrt(1 + (~e / ~d) * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(~c * ~d - ~b * ~e, 0) && Not(GtQ(~a, 0) && GtQ(~d, 0)))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) / sqrt(~d + ~(e') * (~x) ^ 2), ~x) => ((sqrt(~e + ~d / (~x) ^ 2) * sqrt(~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4)) / (~x * sqrt(~d + ~e * (~x) ^ 2) * sqrt(~c + ~b / (~x) ^ 2 + ~a / (~x) ^ 4))) * Antiderivative((~x * sqrt(~c + ~b / (~x) ^ 2 + ~a / (~x) ^ 4)) / sqrt(~e + ~d / (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0))
    @apply_utils Antiderivative(sqrt(~a + ~(c') * (~x) ^ 4) / sqrt(~d + ~(e') * (~x) ^ 2), ~x) => ((sqrt(~e + ~d / (~x) ^ 2) * sqrt(~a + ~c * (~x) ^ 4)) / (~x * sqrt(~d + ~e * (~x) ^ 2) * sqrt(~c + ~a / (~x) ^ 4))) * Antiderivative((~x * sqrt(~c + ~a / (~x) ^ 4)) / sqrt(~e + ~d / (~x) ^ 2), ~x) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0)
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ 2) ^ ~q * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p, ~q], ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IntegerQ(~p) && IntegerQ(~q) || (IGtQ(~p, 0) || IGtQ(~q, 0))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q * (~a + ~(c') * (~x) ^ 4) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ 2) ^ ~q * (~a + ~c * (~x) ^ 4) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~p, ~q], ~x) && (IntegerQ(~p) && IntegerQ(~q) || IGtQ(~p, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~q * (~a + ~(c') * (~x) ^ 4) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~a + ~c * (~x) ^ 4) ^ ~p, (~d / ((~d) ^ 2 - (~e) ^ 2 * (~x) ^ 4) - (~e * (~x) ^ 2) / ((~d) ^ 2 - (~e) ^ 2 * (~x) ^ 4)) ^ -(~q), ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~p], ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (Not(IntegerQ(~p)) && ILtQ(~q, 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(q') * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => Unintegrable((~d + ~e * (~x) ^ 2) ^ ~q * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p, ~q], ~x)
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(q') * (~a + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => Unintegrable((~d + ~e * (~x) ^ 2) ^ ~q * (~a + ~c * (~x) ^ 4) ^ ~p, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~p, ~q], ~x)
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.2 Trinomial products/1.2.2 Quartic/1.2.2.3 (d+e x^2)^q (a+b x^2+c x^4)^p.jl


text:
encode(metadata, x) = x

encode(metadata, x::Dict) =
    [(encode(metadata, k), encode(metadata, v)) for (k, v) in x] |> sort

encode(metadata, x::Union{AbstractArray, Tuple}) =
    map(y -> encode(metadata, y), x)

function encode(metadata, x::Real)
    precision = get(metadata, "precision", 3)
    round(x, precision)
end

function encode(metadata, x::SymTridiagonal)
    "symtridiag(dv:" * JSON.string(x.dv) * ",ev:" *JSON.string(x.ev)
end

@require SymPy begin
    Homework.encode(metadata, x::SymPy.Sym) = "sym-" * string(x)
end


text:
function Vnk_NGG(n, k, β, σ, prec)
    RF = RealField(prec)
    CF = ComplexField(prec)
    β_arb = RF(β)
    σ_arb = RF(σ)
    return exp(β_arb) * σ_arb^(k-1) // RF(gamma(n)) * sum([binom(n-1, i, RF) * (-1)^i * β_arb ^(i//σ_arb) * real(gamma(k-i//CF(σ_arb), CF(β_arb))) for i in 0:(n-1)])
end

function Vnk_2PD(n, k, θ, σ, prec)
    RF = RealField(prec)
    θ_arb = RF(θ)
    σ_arb = RF(σ)
    if k==1
        num::arb = RF(1)
    else
        num = prod([θ_arb + i * σ_arb for i in 1:(k-1)])
    end
    return num // risingfac(1+θ_arb,n-1)
end

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


text:
## for parameters

function _posterior(
    prior::Distributions.DiscreteNonParametric,
    parameters::Vector{<:AbstractParameters},
    responses::Vector{Float64},
    )
    return map( (x, w) ->
            mapreduce( (par, r) ->
            __likelihood(r, x, par)
            ,
            *,
            parameters,
            responses,
            )*w,
        prior.support,
        prior.p
        )
end

function _update_posterior(
    latent::Latent1D,
    parameters::Vector{<:AbstractParameters},
    responses::Vector{Float64}
    )
    likelihood = _posterior(latent.prior, parameters, responses)
    normalizer = sum(likelihood)
    if normalizer > typemin(Float64)
        latent.posterior = Distributions.DiscreteNonParametric(latent.prior.support, likelihood ./ normalizer; check_args = false);
        latent.likelihood = normalizer
    else
        latent.posterior = Distributions.DiscreteNonParametric(latent.prior.support, likelihood; check_args = false);
    end
    return latent::Latent1D
end

## for items

function _posterior(
    prior::Distributions.DiscreteNonParametric,
    items::Vector{<:AbstractItem},
    responses::Vector{Float64}
    )
    return map( (x, w) ->
            mapreduce( (i, r) ->
            _likelihood(r, x, i)
            ,
            *,
            items,
            responses,
            )*w,
        prior.support,
        prior.p
        )
end

function _update_posterior!(
    latent::Latent1D,
    items::Vector{<:AbstractItem},
    responses::Vector{Float64}
    )
    likelihood = _posterior(latent.prior, items, responses)
    normalizer = sum(likelihood)
    if normalizer > typemin(Float64)
        latent.posterior = Distributions.DiscreteNonParametric(latent.prior.support, likelihood ./ normalizer; check_args = false);
        latent.likelihood = normalizer
    else
        latent.posterior = Distributions.DiscreteNonParametric(latent.prior.support, likelihood; check_args = false);
    end
    return nothing
end

function update_posterior!(
    examinee::AbstractExaminee,
    items::Vector{<:AbstractItem},
    responses::Matrix{Union{Missing, Float64}};
    already_sorted::Bool = false,
    kwargs...
    )
    if !already_sorted
        responses_filtered = responses[:, examinee.idx]
        items_filtered = items[.!ismissing.(responses_filtered)]
    else
        responses_filtered = responses
        items_filtered = items
    end
    _update_posterior!(examinee.latent, items_filtered, responses_filtered)
    return nothing
end

function update_posterior!(
    examinee::AbstractExaminee,
    items::Vector{<:AbstractItem},
    responses::Vector{Union{Missing, Float64}};
    kwargs...
    )
    _update_posterior!(examinee.latent, items, responses)
    return nothing
end

function update_posterior!(
    examinees::Vector{<:AbstractExaminee},
    items::Vector{<:AbstractItem},
    responses::Matrix{Union{Missing, Float64}};
    already_sorted::Bool = false,
    kwargs...
)
    map( e -> update_posterior!(e, items, responses; already_sorted = already_sorted, kwargs...), examinees)
    return nothing
end


text:
abstract AbstractAction
botid(a::AbstractAction) = a.bot_id
name(a::AbstractAction) = a.name
position(a::AbstractAction) = a.pos

#########################################################
#  actions that are send to the server
#########################################################

const VALID_ACTIONS = ["move", "radar", "cannon"]

immutable PosAction <: AbstractAction
	bot_id::Int
	pos::Position
	name::ASCIIString
	function PosAction(id, pos, name)
		@assert name ∈ VALID_ACTIONS
		new(id, pos, name)
	end
end

function to_dict(a::PosAction)
	return Dict("botId" => botid(a), "type" => name(a), "pos" => position(a))
end

function make_action(bot_id::Int, x::Int, y::Int, name)
	return PosAction(bot_id, Position(x, y), name)
end

function make_action(bot_id::Int, p::Position, name)
	return PosAction(bot_id, p, name)
end

function MoveAction(args...)
	return make_action(args..., "move")
end

function RadarAction(args...)
	return make_action(args..., "radar")
end

function CannonAction(args...)
	return make_action(args..., "cannon")
end

#########################################################
#  action type for planning that has a associated
#  rating field for its rating
#########################################################

immutable ActionPlan <: AbstractAction
  name::ASCIIString
  pos::Position
  rating::Float64 # rating of this action
  function ActionPlan(name, pos, rating)
    @assert name ∈ VALID_ACTIONS
    new(name, pos, rating)
  end
end
botid(a::ActionPlan) = error("Action plans do not yet have an associated bot!")

# convert ActionPlan to an actual action
make_action(a::ActionPlan, bot::Int) = make_action(bot, position(a), name(a))
make_action(a::ActionPlan, bot::AbstractBot) = make_action(a, botid(bot))

# generate a bunch of actions plans, all actions of type name, at positions pos with ratings rating
plan_actions(name::ASCIIString, pos::Vector{Position}, rating::Vector{Float64}) = [ActionPlan(name, p, w) for (p,w) in zip(pos, rating)]

# returns a randomized copy of the actions plans, were ratings are shifted by [0, rmax].
randomize(old::Vector{ActionPlan}, rmax::Real) =  map(t->ActionPlan(t.name, t.pos, t.rating + rand() * rmax), old)

# new action plans based on old ones, just shifting weights
plan_actions(old::Vector{ActionPlan}, shift::Real) = map(t->ActionPlan(t.name, t.pos, t.rating + shift), old)
function plan_actions(old::Vector{ActionPlan}, shift::Vector{Float64})
  result =  map(x->ActionPlan(x[1].name, x[1].pos, x[1].rating + x[2]), zip(old, shift))
  return convert(Vector{ActionPlan}, result)
end

# gets the best action
best_action(actions::Vector{ActionPlan}, N::Integer = 1) = sort(actions, by = a->a.rating, rev=true)[1]

max_stars_repo_name: Condimentum/spacetyckiting

max_issues_repo_name: Condimentum/spacetyckiting

max_forks_repo_name: Condimentum/spacetyckiting


text:
using JuLIP, QMMM2, SHIPs, PrettyTables, LinearAlgebra, Plots,
      DataFrames, JuLIPMaterials
using JuLIP.MLIPs
include(@__DIR__() * "/swqmmm.jl")

## reference model and MM potential(s) (cf. mmpot_test.jl)
r0 = rnn(:Si)
sw = StillingerWeber()
ships = [ SWqmmm.train_ship(3, deg, wL=1.0, rinfact=0.6, rtol=5e-4)
          for deg in [3,4,5] ]
# all ships have the same cutoff as sw:
rcut = cutoff(sw)

## Solve the reference problem
# domain radius of reference domain ("exact" solution)
Rmax = 16.0 * r0
atmax0 = SWqmmm.int_config(Rmax, sw)
atmax  = SWqmmm.int_config(Rmax, sw)
X0 = positions(atmax0)
E0_max =  energy(atmax0)
minimise!(atmax, precond = FF(atmax, sw))
# add two Newton iterations to properly converge this!
H = lu( hessian(atmax) )
for i = 1:2
   g = gradient(atmax)
   x = dofs(atmax)
   x -= H \ g
   set_dofs!(atmax, x)
end
@show norm(gradient(atmax), Inf)
E1_max = energy(atmax)
dE_max = E0_max - E1_max

## Setup sequence of QM-regions and MM potentials
RQM = [2.0, 3.0, 5.0, 7.0] * r0

errE = zeros(length(RQM), length(ships))
err2 = zeros(length(RQM), length(ships))
errinf = zeros(length(RQM), length(ships))

for n = length(RQM):-1:1, imm = length(ships):-1:1
   @show RQM[n], imm
   at = deepcopy(atmax)
   xc = at["xcore"]
   r = [ norm(x - xc) for x in positions(at) ]
   Iqm = findall(r .< RQM[n])
   at = prepare_qmmm!(at, EnergyMixing;
                      Vqm = sw, Vmm = ships[imm][1], Iqm = Iqm)
   at0 = deepcopy(at)
   set_positions!(at0, X0)
   E0 = energy(at0)
   success = false
   try
      optresult = minimise!(at; precond = FF(at, sw), verbose=2,
                                g_calls_limit=20)
      if optresult.g_residual < 1e-3
         success = true
      end
   catch
      @info("`minimise!` crashed!")
   end
   if !success
      errE[n, imm] = NaN
      err2[n, imm] = NaN
      errinf[n, imm] = NaN
      @info("Optimisation did not converge")
   else
      @info("Optimisation converged; improve with Newton step")
      x = dofs(at) - H \ gradient(at)
      set_dofs!(at, x)
      @show norm(gradient(at), Inf)
      @info("Compute errors")
      E1 = energy(at)
      dE = E0 - E1
      errE[n, imm] = abs(dE - dE_max)

      U = positions(at) - positions(atmax)
      ee = JuLIPMaterials.strains(U, atmax0)
      err2[n, imm] =  norm(ee, 2)
      errinf[n, imm] = norm(ee, Inf)
      @show abs(dE - dE_max)
      @show norm(ee, 2)
   end
end

##

using Printf
println("-"^85)
println("       |" * prod("|         Model$n          " for n = 1:3))
println("  Rqm  |" * "|  E [eV] 2-norm  ∞-norm  "^3)
println("-"^85)
for n =  1:length(RQM)
   @printf(" %5.1f |", RQM[n])
   for imm = 1:length(ships)
      if isnan(errE[n, imm])
         print("|" * "    x   "^3 * " ")
      else
         @printf("| %3.1e %3.1e %3.1e ",
                 errE[n, imm], err2[n, imm], errinf[n, imm])
      end
   end
   println()
end

##
#
# n = 1; imm = 1
# @show RQM[n], imm
# at = deepcopy(atmax)
# xc = at["xcore"]
# r = [ norm(x - xc) for x in positions(at) ]
# Iqm = findall(r .< RQM[n])
# at = prepare_qmmm!(at, EnergyMixing;
#                    Vqm = sw, Vmm = ships[imm][1], Iqm = Iqm)
# at0 = deepcopy(at)
# set_positions!(at0, X0)
# E0 = energy(at0)
# optresult = minimise!(at, precond = FF(at, sw), verbose=2, g_calls_limit=30)
# if optresult.g_residual > 1e-3
#
# E1 = energy(at)
# dE = E0 - E1
# errE[n, imm] = abs(dE - dE_max)
#
# U = positions(at) - positions(atmax)
# ee = JuLIPMaterials.strains(U, atmax0)
# err2[n, imm] =  norm(ee, 2)
# errinf[n, imm] = norm(ee, Inf)
# @show norm(ee, 2)
# @show abs(dE - dE_max)


text:
d1 = domain([4,3,2,1])
d2 = domain(1:4)
domains = Dictionary(1:2, [d1, d2])
@testset "Internals: Domains" begin
    for d in domains
        # constructors and ∈
        for x in [1,2,3,4]
            @test x ∈ d
        end
        # length
        @test length(d) == 4
        # draw and ∈
        @test rand(d) ∈ d
    end
    # add!
    ConstraintDomains.add!(d1, 5)
    @test 5 ∈ d1
    # delete!
    delete!(d1, 5)
    @test 5 ∉ d1
end

x1 = variable([4,3,2,1])
x2 = variable(d2)
x3 = variable() # TODO: tailored test for free variable
vars = Dictionary(1:2, [x1, x2])
@testset "Internals: variables" begin
    for x in vars
        # add and delete from constraint
        LS._add_to_constraint!(x, 1)
        LS._add_to_constraint!(x, 2)
        LS._delete_from_constraint!(x, 2)
        @test x ∈ 1
        @test x ∉ 2
        @test LS._constriction(x) == 1
        @test length(x) == 4
        for y in [1,2,3,4]
            @test y ∈ x
        end
        @test rand(x) ∈ x
    end
    add!(x1, 5)
    @test 5 ∈ x1
    delete!(x1, 5)
    @test 5 ∉ x1
end

values = [1, 2, 3]
inds = [1, 2]
err = error_f(usual_constraints[:all_different])
c1 = constraint(err, inds)
c2 = constraint(err, inds)
cons = Dictionary(1:2, [c1, c2])
@testset "Internals: constraints" begin
    for c in cons
        LS._add!(c, 3)
        @test 3 ∈ c
        LS._delete!(c, 3)
        @test 3 ∉ c
        @test LS._length(c) == 2
        c.f(values, Matrix{Float64}(undef, 3, CompositionalNetworks.max_icn_length()); dom_size=3)
    end
end

o1 = objective(sum, "Objective 1: sum")
o2 = objective(prod, "Objective 2: product")
objs = Dictionary(1:2, [o1, o2])
@testset "Internals: objectives" begin
    for o in objs
        @test o.f(values) == 6
    end
end

m = model()
describe(m)

x1 = variable([4,3,2,1])
x2 = variable(d2)
vars = Dictionary(1:2, [x1, x2])

values = [1, 2, 3]
inds = [1, 2]
c1 = constraint(err, inds)
c2 = constraint(err, inds)
cons = Dictionary(1:2, [c1, c2])

o1 = objective(sum, "Objective 1: sum")
o2 = objective(prod, "Objective 2: product")
objs = Dictionary(1:2, [o1, o2])

@testset "Internals: model" begin
    for x in vars
        add!(m, x)
    end
    variable!(m, d1)

    for c in cons
        add!(m, c)
    end
    constraint!(m, err, [1,2])

    for o in objs
        add!(m, o)
    end

    objective!(m, max)
    length_var(m, 1)
    length_cons(m, 1)
    constriction(m, 1)
    draw(m, 1)
    get_objective(m, 1)
    delete_value!(m, 1, 1)
    add_value!(m, 1, 1)
    delete_var_from_cons!(m, 1, 1)
    add_var_to_cons!(m, 1, 1)
    describe(m)
end

## Test Solver
s1 = solver()

LS.get_error(LS.EmptyState())