text:
using Flux, Gym
using Flux.Optimise: Optimiser
using Flux.Tracker: data
using Statistics: mean
using DataStructures: CircularBuffer
using Distributions: sample
using Printf
#using CuArrays

# Load game environment
env = make("CartPole-v0", :human_pane)
reset!(env)

# ----------------------------- Parameters -------------------------------------

STATE_SIZE = length(env._env.state)    # 4
ACTION_SIZE = length(env._env.action_space) # 2
MEM_SIZE = 100_000
BATCH_SIZE = 64
γ = 1f0   			  # discount rate

# Exploration params
ϵ = 1f0       # Initial exploration rate
ϵ_MIN = 1f-2    # Final exploratin rate
ϵ_DECAY = 995f-3

# Optimiser params
η = 1f-2   # Learning rate
η_decay = 1f-2

memory = CircularBuffer{Any}(MEM_SIZE) # Used to remember past results

# ------------------------------ Model Architecture ----------------------------

model = Chain(Dense(STATE_SIZE, 24, tanh),
              Dense(24, 48, tanh),
              Dense(48, ACTION_SIZE)) |> gpu

loss(x, y) = Flux.mse(model(x |> gpu), y)

opt = Optimiser(ADAM(η), InvDecay(η_decay))

# ----------------------------- Helper Functions -------------------------------

get_ϵ(e) = max(ϵ_MIN, min(ϵ, 1f0 - log10(e * ϵ_DECAY)))

remember(state, action, reward, next_state, done) =
  push!(memory, (data(state), action, reward, data(next_state), done))

function action(state, train=true)
  train && rand() <= get_ϵ(e) && (return Gym.sample(env.action_space))
  act_values = model(state |> gpu)
  return Flux.onecold(act_values)
end

function replay()
  global ϵ
  batch_size = min(BATCH_SIZE, length(memory))
  minibatch = sample(memory, batch_size, replace = false)

  x = []
  y = []
  for (iter, (state, action, reward, next_state, done)) in enumerate(minibatch)
    target = reward
    if !done
      target += γ * maximum(data(model(next_state |> gpu)))
    end

    target_f = data(model(state |> gpu))
    target_f[action] = target

    push!(x, state)
    push!(y, target_f)
  end
  x = hcat(x...) |> gpu
  y = hcat(y...) |> gpu
  Flux.train!(loss, params(model), [(x, y)], opt)

  ϵ *= ϵ > ϵ_MIN ? ϵ_DECAY : 1.0f0
end

function episode!(env, train=true, draw=false)
  done = false
  total_reward = 0f0
  while !done
    draw && render!(env)
    s = env._env.state
    a = action(s, train)
    s′, r, done, _ = step!(env, a)
    total_reward += r
    train && remember(s, a, r, s′, done)
    draw && sleep(0.01)
  end

  total_reward
end

# ------------------------------ Training --------------------------------------

e = 1
scores = CircularBuffer{Float32}(100)

while true
  global e
  reset!(env)
  total_reward = episode!(env)
  push!(scores, total_reward)
  print("Episode: $e | Score: $total_reward ")
  last_100_mean = mean(scores)
  print("Last 100 episodes mean score: $(@sprintf "%6.2f" last_100_mean)")
  if last_100_mean > 195
    println("\nCartPole-v0 solved!")
    break
  end
  println()
  replay()
  e += 1
end

# -------------------------------- Testing -------------------------------------

ee = 1

while true
  global ee
  reset!(env)
  total_reward = episode!(env, false, true)
  println("Episode: $ee | Score: $total_reward")
  ee += 1
end


text:
"""

    lad(setting; starting_betas = nothing)

Perform Least Absolute Deviations regression for a given regression setting.

# Arguments
- `setting::RegressionSetting`: RegressionSetting object with a formula and dataset.
- `starting_betas::Array{Float64,1}`: Starting values of parameter estimations that fed to local search optimizer.

# Examples
```julia-repl
julia> reg0001 = createRegressionSetting(@formula(calls ~ year), phones);
julia> lad(reg0001)
Dict{Any,Any} with 2 entries:
  "betas"     => [-57.3269, 1.19155]
  "residuals" => [2.14958, 1.25803, 0.0664872, 0.0749413, -0.416605, -0.90815, -1.2997, -1.79124,…

```

"""
function lad(setting::RegressionSetting; starting_betas=nothing)
    X, y = @extractRegressionSetting setting
    return lad(X, y, starting_betas=starting_betas)
end

"""

    lad(X, y; starting_betas = nothing)

Perform Least Absolute Deviations regression for a given regression setting.

# Arguments
- `X::Array{Float64, 2}`: Design matrix of the linear model.
- `y::Array{Float64, 1}`: Response vector of the linear model.
- `starting_betas::Array{Float64,1}`: Starting values of parameter estimations that fed to local search optimizer.

"""
function lad(X::Array{Float64,2}, y::Array{Float64,1}; starting_betas=nothing)
    n, p = size(X)

    if starting_betas isa Nothing
        starting_betas = zeros(Float64, p)
    end

    function goal(betas::Array{Float64,1})::Float64
        sum(abs.(y .- X * betas))
    end

    optim_result = optimize(goal, starting_betas, NelderMead())
    betas = optim_result.minimizer
    residuals = y .- X * betas

    result = Dict()
    result["betas"] = optim_result.minimizer
    result["residuals"] = residuals

    return result
end

ext: jl

max_stars_repo_name: angeris/LinRegOutliers

max_issues_repo_name: angeris/LinRegOutliers

max_forks_repo_name: angeris/LinRegOutliers


text:
# kuramoto_8-19-19.jl
# I've worked my way up to tasks/coroutines in the manual.
#
# https://docs.julialang.org/en/v1/manual/control-flow/#man-tasks-1
#
# Pretty wild stuff, but I need to play around with them
# a bit to really get a handle on how they work.

function producer(c::Channel)
    put!(c, "Channel initialized.")
    for n = 1:100
        sleep(5) # Hang for 5 seconds, arbitrarily.
        put!(c, n)
    end
    put!(c, "Channel completed.")
end

# How long does this code take to run? are the producer
# tasks always running in the background, generating and
# storing values, or is it more like lazy evaluation? I'll
# need to do some experiments to find out.

#----------------------------------------------------------#

# First, the most basic test.

channel_1 = Channel(producer)

@time begin
    @time println(take!(channel_1))
    @time println(take!(channel_1))
    @time println(take!(channel_1))
end

#----------------------------------------------------------#

channel_2 = Channel(producer)

@time begin
    @time println(take!(channel_2))
    @time sleep(2)
    @time println(take!(channel_2))
    @time sleep(2)
    @time println(take!(channel_2))
end

# Oh, holy shit, no. The producers DO run in the background,
# even when other functions are working. Both of these run
# at around 10 seconds apiece.

#----------------------------------------------------------#

channel_3 = Channel(producer)

@time begin
    @time println(take!(channel_3))
    @time sleep(10)
    @time println(take!(channel_3))
    @time sleep(10)
    @time println(take!(channel_3))
end

# Final test. If this runs at 20 seconds, it's evidence that
# the producers are running even when I'm not take!ing from
# them. If this runs at 30 seconds, it's evidence that they
# only fire up to produce the next value once I call the
# take! function. If it's something else... I don't know.

# 20 seconds. So the producers can run in the background.

#----------------------------------------------------------#

# Okay, next question. We know it runs in the background
# until the next call. But does it run until it has a value
# to return, and *then* hang? Or does it keep running so it
# has a stockpile of values just in case? My guess is
# actually the latter, but I wonder if there's some way to
# force the behavior of the former. Let's test it.

# If it stores multiple values, then a sleep value larger
# than all of the 5-second delays should lead to the values
# all spilling out at once. If it's only storing one value,
# this should take about 55 seconds, instead of 45.

channel_4 = Channel(producer)

@time begin
    @time println(take!(channel_4))
    @time sleep(45)
    @time println(take!(channel_4))
    @time println(take!(channel_4))
    @time println(take!(channel_4))
end

# 55 seconds! Aha! So it only works up until it put!s the
# next value, then it hands off to the next lad.

#----------------------------------------------------------#

# Okay, one last test. I think this one is pretty
# straightforward.

function producer_slow(c::Channel)
    put!(c, "Slow producer online!")
    while true
        sleep(7)
        put!(c, "x")
    end
end

function producer_fast(c::Channel)
    put!(c, "Fast producer online!")
    while true
        sleep(5)
        put!(c, "y")
    end
end

channel_x = Channel(producer_slow)
channel_y = Channel(producer_fast)

@time begin
    for i = 1:5
        @time println(take!(channel_y))
        @time println(take!(channel_x))
    end
end

#----------------------------------------------------------#

# Let's finally make some code for practicing actually
# making tasks.

producer(c::Channel, t) = while true put!(c, t) end

channel = Channel(producer)

take!(channel)

# MethodError, of course.

#----------------------------------------------------------#

producer(c::Channel, t) = while true put!(c, t) end

channel = Channel((c::Channel) -> producer(c, 5))

take!(channel)
take!(channel)
take!(channel)

#----------------------------------------------------------#

producer(c::Channel, t) = while true put!(c, t) end

channel = Channel(() -> producer(5))

take!(channel)
take!(channel)
take!(channel)

#----------------------------------------------------------#

producer(c::Channel, t) = while true put!(c, t) end

channel = Channel(producer(c,5))

take!(channel)
take!(channel)
take!(channel)

#----------------------------------------------------------#

# Okay, now let's get some practice with creating Tasks.

f(x) = sleep(x);
g(y) = sleep(y);

t1 = Task(() -> f(2));
t2 = @task g(3);

#----------------------------------------------------------#


text:
using Gadfly

plot(x=rand(20), y=rand(20), color=vcat(fill("A", 10), fill("B", 10)),
     Guide.colorkey("Species"))


text:
#=
object.jl
Types for making (often) piece-wise constant phantoms
out of elementary objects having shapes like
ellipses, rectangles, and gaussian bumps.
=#

export AbstractObject
export Object, Object2d, Object3d
#export rotate, scale, translate # no: to avoid conflicts with Plots
export phantom, radon, spectrum

atuple(x::Any, n::Int) = ntuple(i -> x, n)

abstract type AbstractObject end

"""
    Object{S, D, V, ...}(center, width, angle, value, param) <: AbstractObject
General container for 2D and 3D objects for defining image phantoms.

* `center::NTuple{D,C}` coordinates of "center" of this object
* `width::NTuple{W,C}` "width" along axis; usually `W=D` or `W=1`
  (e.g., FWHM for Gauss, radii for Ellipse, radius of Circle)
* `angle::NTuple{D-1,A}` angle of x' axis relative to x axis,
  in radians (or with units)
* `value::V` "intensity" value for this object
* `param` optional additional shape parameters (often `nothing`)

# Example

```jldoctest
julia> Object(Ellipse(), (0,0), (1,2), 0.0, 1//2, nothing)
Object2d{Ellipse, Rational{Int64}, 2, Int64, Float64, Nothing} (S, D, V, ...)
 shape::Ellipse Ellipse()
 center::NTuple{2,Int64} (0, 0)
 width::NTuple{2,Int64} (1, 2)
 angle::Tuple{Float64} (0.0,)
 value::Rational{Int64} 1//2
 param::Nothing nothing
```
"""
struct Object{S, D, V, W, C, A, P} <: AbstractObject
    shape::S
    "x,y center coordinates"
    center::NTuple{D,C}
    "'width' along x',y' axes (FWHM for Gauss, radii for Ellipse)"
    width::NTuple{W,C}
    "angle of x' axis relative to x axis, in radians (or with units)"
    angle::NTuple{Da,A} where Da
    "'intensity' value for this shape"
    value::V
    "optional additional shape parameters"
    param::P # often `nothing`

    """
        shape = Object(shape, center, width, angle, value, param)
    Top-level shape constructor where `angle` must be tuple.
    """
    function Object{S}(
        shape::S,
        center::NTuple{D,C},
        width::NTuple{W,C},
        angle::NTuple{Da,A},
        value::V,
        param::P,
    ) where {S <: AbstractShape, D, V <: Number, W, C <: RealU, Da, A <: RealU, P}
        1 ≤ Da == D-1 || throw(ArgumentError("Dϕ=$Dϕ != D-1, where D=$D"))
        all(width .> zero(eltype(width))) || throw(ArgumentError("widths must be positive"))
        new{S,D,V,W,C,A,P}(shape, center, width, angle, value, param)
    end
end

"""
    Object2d = Object{S,2} where S <: AbstractObject
For 2D objects
"""
const Object2d = Object{S,2} where S

"""
    Object3d = Object{S,3} where S <: AbstractObject
For 3D objects
"""
const Object3d = Object{S,3} where S

# Constructors

"""
    ob = Object(shape, center, width, angle, value, param)
General outer object constructor where `angle` is a tuple, including 3D case.
"""
function Object(
    shape::S,
    center::NTuple{D,C},
    width::NTuple{W,C} where W = atuple(zero(eltype(center)), length(center)),
    angle::NTuple{Da,RealU} where Da = atuple(0, length(center)-1),
    value::Number = 1,
    param::Any = nothing,
) where {S <: AbstractShape, D, C <: RealU}
    Object{S}(shape, center, width, angle, value, param)
end

"""
    ob = Object(shape, center, width, angle=0, value=1, param=nothing)
2D object constructor.
"""
function Object(
    shape::AbstractShape2,
    center::NTuple{2,C},
    width::NTuple{W,C} where W,
    angle::RealU = 0,
    value::Number = 1,
    param::Any = nothing,
) where {C <: RealU}
    Object(shape, center, width, (angle,), value, param)
end

"""
    ob = Object(shape ; center=(0,0), width=(1,1), angle=0, value=1, param=nothing)
2D object constructor using keywords.
"""
function Object(
    shape::AbstractShape2 ;
    center::NTuple{2,C} = (0,0),
    width::NTuple{W,C} where W = atuple(one(eltype(center)), 2),
    angle::RealU = 0,
    value::Number = 1,
    param = nothing,
) where {C <: RealU}
    Object(shape, center, width, angle, value, param)
end

#=
# not needed due to general case above
"""
    ob = Object(shape, center, width, angle=(0,0), value=1, param=nothing)
3D object constructor.
"""
function Object(
    shape::Shape3d,
    center::NTuple{3,C},
    width::NTuple{W,C} where W,
    angle::NTuple{2,RealU} = (0,0),
    value::Number = 1,
    param::Any = nothing,
) where {C <: RealU}
    Object(shape, center, width, angle, value, param)
end
=#

"""
    ob = Object(shape ; center=(0,0,0), width=(1,1,1), angle=(0,0), value=1, param=nothing)
3D object constructor using keywords.
"""
function Object(
    shape::AbstractShape3 ;
    center::NTuple{3,C} = (0,0,0),
    width::NTuple{W,C} = atuple(one(eltype(center)), 3),
    angle::NTuple{2,RealU} = (0,0),
    value::Number = 1,
    param = nothing,
) where {W, C <: RealU}
    Object(shape, center, width, angle, value, param)
end

#=
"""
    Object(shape ; x, y, rx, ry, ϕ=0, value=1 ; param=nothing)
2D object constructor without tuples.
"""
function Object(
    shape::AbstractShape,
    cx::C, cy::C,
    rx::C, ry::C,
    ϕ::RealU,
    value::Number = 1 ;
    param = nothing,
) where {C <: RealU}
# todo: might want to use promote
    Object(shape, (cx, cy), (rx, ry), ϕ, value, param)
end
=#

"""
    Object(ob::Object ; center, width, angle, value, param)
Make a copy of Object `ob`, optionally modifying some values.
"""
function Object(ob::Object{S,D} ;
    center::NTuple{D} = ob.center,
    width::NTuple{D} = ob.width,
    angle::NTuple{Da,<:RealU} = ob.angle,
    value::Number = ob.value,
    param = ob.param,
) where {S, D, Da}
    Da == D-1 || throw(ArgumentError("Dϕ=$Dϕ != D-1, where D=$D"))
    Object(ob.shape, center, width, angle, value, param)
end

# Methods for objects

Base.eltype(::Object{S,D,V}) where {S,D,V} = V
Base.ndims(::Object{S,D}) where {S,D} = D

"""
    show(io::IO, ::MIME"text/plain", ob::Object)
"""
function Base.show(io::IO, ::MIME"text/plain", ob::Object{S,D}) where {S,D}
    println(io, typeof(ob), " (S, D, V, ...)")
    for f in (:shape, :center, :width, :angle, :value, :param)
        p = getproperty(ob, f)
        t = typeof(p)
        t = t == NTuple{D,eltype(t)} ? "NTuple{$D,$(eltype(t))}" : "$t"
        println(io, " ", f, "::", t, " ", p)
    end
end

# rotate
"""
    rotate(ob::Object2d, θ::RealU)
In-plane rotation of a 2D object.
"""
rotate(ob::Object2d, θ::RealU) =
    Object(ob.shape, ob.center, ob.width, ob.angle .+ (θ,), ob.value, ob.param)
#rotate(ob::Object3d, θ::RealU) =
#Object(ob.shape, ob.center, ob.width, ob.angle .+ (θ,0), ob.value, ob.param)

"""
    rotate(ob::Object3d, (α,β))
    rotate(ob::Object3d, α, β=0)
Rotation of a 3D object.
"""
rotate(ob::Object3d, θ::NTuple{2,RealU}) =
    Object(ob.shape, ob.center, ob.width, ob.angle .+ θ, ob.value, ob.param)
rotate(ob::Object3d, α::RealU, β::RealU=0) = rotate(ob, (α,β))

# scale width (tricky because length(width) < D is possible
"""
    scale(ob::Object, factor::RealU)
    scale(ob::Object, factor::NTuple{W,RealU})
Scale the width(s) by `factor`.
"""
scale(ob::Object{S,D,V,W}, factor::NTuple{W,RealU}) where {S,D,V,W} =
    Object(ob.shape, ob.center, ob.width .* factor, ob.angle, ob.value, ob.param)
scale(ob::Object{S,D,V,W}, factor::RealU) where {S,D,V,W} =
    scale(ob, atuple(factor,W))

# scale value
"""
    (*)(ob::Object, x::Number)
    (*)(x::Number, ob::Object)
Scale object `value` by `x`.
"""
Base.:(*)(ob::AbstractObject, x::Number) =
    Object(ob.shape, ob.center, ob.width, ob.angle, ob.value * x, ob.param)
Base.:(/)(ob::Object, x::Number) = ob * (1 / x)
Base.:(*)(x::Number, ob::Object) = ob * x

# translate
"""
    translate(ob::Object, shift::NTuple{D,RealU})
    translate(ob::Object2d, xshift, yshift)
    translate(ob::Object3d, xshift, yshift, zshift)
Translate the center coordinates of an object by `shift`
"""
translate(ob::Object{S,D}, shift::NTuple{D,RealU}) where {S, D} =
    Object(ob.shape, ob.center .+ shift, ob.width, ob.angle, ob.value, ob.param)
translate(ob::Object3d, x::RealU, y::RealU, z::RealU) = translate(ob, (x,y,z))
translate(ob::Object2d, x::RealU, y::RealU) = translate(ob, (x,y))

"""
    coords(object::Object2d, x::RealU, y::RealU)
Put coordinates `(x,y)` in canonical axes associated with `object`.
"""
function coords(ob::Object2d, x::RealU, y::RealU)
    (x, y) = rotate2d(x - ob.center[1], y - ob.center[2], ob.angle[1])
    return (x / ob.width[1], y / ob.width[2]) # unitless
end

# methods for phantoms: an array of objects

"""
    image = phantom(oa::Array{<:Object2d})::Function
Return function `image(x,y)` that user can sample at any `(x,y)` locations
to make a phantom image.
"""
function phantom(oa::Array{<:Object2d})
    return (x,y) -> sum(ob -> phantom(ob)(x,y), oa)
end

"""
    image = phantom(x, y, oa::Array{<:Object2d}, oversample::Int; T)
Return a digital image of the phantom sampled at `(x,y)` locations,
with over-sampling factor `oversample` and element type `T`.
"""
function phantom(
    x::AbstractVector,
    y::AbstractVector,
    oa::Array{<:Object2d},
    oversample::Int;
    T::DataType = promote_type(eltype.(oa)..., Float32),
)
    oversample < 1 && throw(ArgumentError("oversample $oversample"))
    dx = x[2] - x[1]
    dy = y[2] - y[1]
    all(≈(dx), diff(x)) || throw("oversample requires uniform x")
    all(≈(dy), diff(y)) || throw("oversample requires uniform y")
    tmp = ((1:oversample) .- (oversample+1)/2) / oversample
    ophantom = ob ->
       (x,y) -> T(sum(phantom(ob).(x .+ dx*tmp, y .+ dy*tmp')) / abs2(oversample))
    return sum(ob -> ophantom(ob).(x,y'), oa)
end

"""
    image = phantom(x, y, oa::Array{<:Object2d})
Return a digital image of the phantom sampled at `(x,y)` locations.
"""
function phantom(
    x::AbstractVector,
    y::AbstractVector,
    oa::Array{<:Object2d},
)

    return sum(ob -> phantom(ob).(x,y'), oa)
end

function phantom(xx::AbstractArray, yy::AbstractArray, oa::Array{<:Object2d})
    return sum(ob -> phantom(ob).(xx,yy), oa)
end

"""
    sino = radon(oa::Array{<:Object2d})::Function
Return function `sino(r,ϕ)` that user can sample at any `(r,ϕ)` locations
to make a phantom 2D sinogram.

The coordinate system used here is such that `ϕ=0` corresponds to
line integrals along the ``y`` axis for an object ``f(x,y)``.
Then as `ϕ` increases, the line integrals rotate counter-clockwise.
"""
function radon(oa::Array{<:Object2d})
    return (r,ϕ) -> sum(ob -> radon(ob)(r,ϕ), oa)
end

"""
    sino = radon(r, ϕ, oa::Array{<:Object2d})
Return parallel-beam 2D sinogram `sino` sampled at given `(r,ϕ)` locations.
"""
function radon(r::AbstractVector, ϕ::AbstractVector, oa::Array{<:Object2d})
    return sum(ob -> radon(ob).(r,ϕ'), oa)
end

function radon(rr::AbstractArray, ϕϕ::AbstractArray, oa::Array{<:Object2d})
    return sum(ob -> radon(ob).(rr,ϕϕ), oa)
end

"""
    kspace = spectrum(oa::Array{<:Object2d})::Function
Return function `kspace(fx,fy)` that user can sample at any `(fx,fy)` locations
to make phantom 2D k-space data.
"""
function spectrum(oa::Array{<:Object2d})
    return (fx,fy) -> sum(ob -> spectrum(ob)(fx,fy), oa)
end

"""
    kspace = spectrum(fx, fy, oa::Array{<:Object2d})::Function
Return k-space matrix `kspace` sampled at given `(fx,fy)` locations.
"""
function spectrum(fx::AbstractVector, fy::AbstractVector, oa::Array{<:Object2d})
    return sum(ob -> spectrum(ob).(fx,fy'), oa)
end

function spectrum(fx::AbstractArray, fy::AbstractArray, oa::Array{<:Object2d})
    return sum(ob -> spectrum(ob).(fx,fy), oa)
end

# helpers

function rotate2d(x::RealU, y::RealU, θ::RealU)
    (s, c) = sincos(θ)
    return (c * x + s * y, -s * x + c * y)
end

rotate2d(xy::NTuple{2,RealU}, θ::RealU) = rotate2d(xy..., θ)


text:
# Cholesky decomposition

L = zeros(A)

for k=1:m
    lk = L[1:k-1,1:k-1]\A[1:k-1,k]
    L[k,1:k-1] = lk ’
    L[k:k,k] = sqrt(A[k:k,k] - lk ’*lk)
end

max_stars_repo_path: Code/Numerical Analysis & Machine Learning/Numerical Linear Algebra/CholeskyDecompostion.jl

max_stars_repo_name: BambooFlower/Math-Scripts

max_issues_repo_name: BambooFlower/Math-Scripts

max_forks_repo_name: BambooFlower/Math-Scripts


text:
module Diffusion
using StaticArrays
using CoordinateTransformations

export absorb_boundary
export reflect_boundary
export simmolreflect
export simmolabsorb
export nmolreflect
export nmolabsorb

###    boundary conditions     ###

function absorb_boundary(pos,r) #models an absorbing boundary-pos inside the boundary stay, but pos outside disappear
    if(sum(abs2, pos) < r^2)
        return pos
    else
        return nothing
    end
end

function topolar(pos) #turns a cartesian coordinate into a polar coordinate
    if pos[1] >= 0 && pos[2] >= 0
        theta = atan(pos[2]/pos[1])
    elseif pos[1] >=0 && pos[2] <= 0
        theta = atan(pos[2]/pos[1]) + 2π
    else
        theta = atan(pos[2]/pos[1]) + π
    end
    r = sqrt(sum(abs2, pos))
    return SA[r,theta]
end

function tocart(pos) #turns a polar coordinate into a cartesian coordinate
    x = pos[1]*cos(pos[2])
    y = pos[1]*sin(pos[2])
    return SA[x,y]
end

function reflect(pos,r) #reflects a point radially over a circle
    pos1 = topolar(pos)
    dist = pos1[1] - r
    pos2 = SA[r - dist,pos1[2]]
    posf = tocart(pos2)
    return posf
end

function reflect_boundary(pos,r) #models a reflective boundary. pos inside the boundary stay, pos outside the boundary are reflected back in
    if(sum(abs2, pos) <= r^2)
        return pos
    else
        pos = reflect(pos,r)
        return reflect_boundary(pos, r)
    end
end

###   simulating the movement of molecules inside a particular boundary   ###

function movemol(pos, σ) #moves the molecule once
    return pos + SA[σ*randn(),σ*randn()]
end

#=
    simmol functions
    inputs: σ aka sqrt(2*d*deltat) aka the sd of the displacement, number of steps, initial position, radius
    outputs: vector of static arrays showing where the molecule has been
=#

function simmolreflect(σ, steps, pos, r) #simmol under refelctive boundary conditions
    molhistory = [pos]
    for x in 1:steps
        pos = movemol(pos, σ)
        pos = reflect_boundary(pos, r)
        push!(molhistory, pos)
    end
    return molhistory
end

function simmolabsorb(σ, steps, pos, r) #simmol under absorbing boundary conditions
    molhistory = [pos]
    for x in 1:steps
        pos = movemol(pos, σ)
        pos = absorb_boundary(pos, r)
        if pos === nothing
            return molhistory
        else
        push!(molhistory, pos)
        end
    end
    return molhistory
end

function nmolreflect(σ, steps, pos, r, molecules) #simulates multiple molecules at once-reflective
    molhistories = Vector{typeof(pos)}[]
    for x in 1:molecules
        push!(molhistories, simmolreflect(σ, steps,pos,r))
    end
    return molhistories
end

function nmolabsorb(σ, steps, pos, r, molecules) #simulates multiple molecules at once-absorbing
    molhistories = Vector{typeof(pos)}[]
    for x in 1:molecules
        push!(molhistories, simmolabsorb(σ, steps,pos,r))
    end
    return molhistories
end

end


text:
# Apply default faceting look to a grid of `AxisEntries`

function facet_wrap!(fig, aes::AbstractMatrix{AxisEntries})
    scale = get(aes[1].scales, :layout, nothing)
    isnothing(scale) && return aes
    linkaxes!(aes...)
    for ae in aes
        ax = Axis(ae)
        entries = ae.entries
        vs = Iterators.filter(
            !isnothing,
            (get(entry.mappings, :layout, nothing) for entry in entries)
        )
        it = iterate(vs)
        if isnothing(it)
            delete!(ax)
        else
            v, _ = it
            ax.title[] = string(first(v))
        end
    end
    return nothing
end

# F. Greimel implementation from AlgebraOfGraphics

# TODO: add configuration options here
function facet_grid!(fig, aes::AbstractMatrix{AxisEntries})
    M, N = size(aes)
    row_scale, col_scale = map(sym -> get(aes[1].scales, sym, nothing), (:row, :col))
    all(isnothing, (row_scale, col_scale)) || hideinnerdecorations!(aes)
    if !isnothing(row_scale)
        linkxaxes!(aes...)
        for ae in aes
            Axis(ae).ylabelvisible[] = false
        end
        row_dict = Dict(zip(row_scale.plot, row_scale.data))
        for m in 1:M
            Box(fig[m, N, Right()], color=:gray85, strokevisible=true)
            Label(fig[m, N, Right()], string(row_dict[m]); rotation=-π/2)
        end
        protrusion = lift(
            (xs...) -> maximum(x -> x.left, xs),
            (MakieLayout.protrusionsobservable(Axis(ae)) for ae in aes[:, 1])...
        )
        # TODO: here and below, set in such a way that one can change padding after the fact?
        padding = lift(protrusion, Axis(aes[1]).ylabelpadding) do val, p
            return (0f0, val + p, 0f0, 0f0)
        end
        Label(fig[:, 1, Left()], Axis(aes[1]).ylabel; rotation=π/2, padding)
    end
    if !isnothing(col_scale)
        linkyaxes!(aes...)
        for ae in aes
            Axis(ae).xlabelvisible[] = false
        end
        col_dict = Dict(zip(col_scale.plot, col_scale.data))
        for n in 1:N
            Box(fig[1, n, Top()], color=:gray85, strokevisible=true)
            Label(fig[1, n, Top()], string(col_dict[n]))
        end
        protrusion = lift(
            (xs...) -> maximum(x -> x.bottom, xs),
            (MakieLayout.protrusionsobservable(Axis(ae)) for ae in aes[M, :])...
        )
        padding = lift(protrusion, Axis(aes[1]).xlabelpadding) do val, p
            return (0f0, 0f0, 0f0, val + p)
        end
        Label(fig[M, :, Bottom()], Axis(aes[1]).xlabel; padding)
    end
    return nothing
end

function facet!(fig, aes::AbstractMatrix{AxisEntries})
    facet_wrap!(fig, aes)
    facet_grid!(fig, aes)
    return nothing
end

function facet!(fg::FigureGrid)
    facet!(fg.figure, fg.grid)
    return fg
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   The Jacchia-Roberts 1971 Atmospheric Model.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Roberts, C. R (1971). An analytic model for upper atmosphere densities
#       based upon Jacchia's 1970 models.
#
#   [2] Jacchia, L. G (1970). New static models of the thermosphere and
#       exosphere with empirical temperature profiles. SAO Special Report #313.
#
#   [3] Long, A. C., Cappellari Jr., J. O., Velez, C. E., Fuchs, A. J (editors)
#       (1989). Goddard Trajectory Determination System (GTDS) Mathematical
#       Theory (Revision 1). FDD/552-89/0001 and CSC/TR-89/6001.
#
#   [4] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
#   [5] GMT-2668 - http://bugs.gmatcentral.org/browse/GMT-2668
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export jr1971

include("./jr1971_priv.jl")

"""
    jr1971(JD::Number, glat::Number, glon::Number, h::Number, F10::Number, F10ₐ::Number, Kp::Number)

Compute the atmospheric density using the Jacchia-Roberts 1971 model.

# Args

* `JD`: Julian day.
* `glat`: Geodetic latitude [rad].
* `glon`: Geodetic longitude [rad].
* `h`: Altitude [m].
* `F10`: 10.7-cm solar flux \\[10⁻²² W/(M² Hz)].
* `F10ₐ`: 10.7-cm averaged solar flux, 81-day centered on input time.
* `Kp`: Kp geomagnetic index (with a delay of 3 hours).

# Returns

An instance of the structure `JR1971_Output` with the computed values.

"""
function jr1971(JD::Number, glat::Number, glon::Number, h::Number, F10::Number,
                F10ₐ::Number, Kp::Number)

    # Constants
    # ==========================================================================

    @unpack_JR1971_CONSTANTS _jr1971_constants

    # Auxiliary variables
    # ==========================================================================

    Ra² = Ra*Ra       # Mean Earth radius squared [km²].

    # Preliminaries
    # ==========================================================================

    # Convert the altitude from [m] to [km].
    h /= 1000

    # Compute the Sun position represented in the inertial reference frame
    # (MOD).
    s_i = sun_position_i(JD)

    # Compute the Sun declination [rad].
    δs = atan( s_i[3], sqrt(s_i[1]*s_i[1] + s_i[2]*s_i[2]) )

    # Compute the Sun right ascension [rad].
    Ωs = atan( s_i[2], s_i[1] )

    # Compute the right ascension of the selected location w.r.t. the inertial
    # reference frame.
    Ωp = glon + JDtoGMST(JD)

    # Compute the hour angle at the selected location, which is the angle
    # measured at the XY plane between the right ascension of the selected
    # position and the right ascension of the Sun.
    H = Ωp - Ωs

    # Algorithm
    # ==========================================================================

    # Exospheric temperature
    # ==========================================================================

    # Diurnal variation
    # =================

    # Eq. 14 [2], eq. 4-89 [3]
    # --------------------------------------------------------------------------
    #
    # Nighttime minimum of the global exospheric temperature distribution when
    # the planetary geomagnetic index Kp is zero.
    ΔF10 =  F10 - F10ₐ
    Tc   = 379 + 3.24F10 + 1.3ΔF10

    # Eq. 15 [2], eq. 4-91 [3]
    # --------------------------------------------------------------------------
    η = abs(glat - δs)/2
    θ = abs(glat + δs)/2

    # Eq. 16 [2], eq. 4-91 [3]
    # --------------------------------------------------------------------------
    τ = H + ( -37 + 6sin(H + 43*π/180) )*π/180

    # Eq. 17 [2], eq. 4.90 [3]
    # --------------------------------------------------------------------------
    m = 2.2
    n = 3.0
    R = 0.3
    C = cos(η)^m
    S = sin(θ)^m

    Tl = Tc*( 1 + R*( S + (C-S)*cos(τ/2)^n ) )

    # Variations with geomagnetic activity
    # ====================================

    # Eq. 18 or Eq. 20 [2], eq. 4-94 [3]
    # --------------------------------------------------------------------------

    ΔT∞ = (h < 200) ? 14Kp + 0.02exp(Kp) : 28Kp + 0.03exp(Kp)

    # Eq. 4-95 [3]
    # --------------------------------------------------------------------------
    #
    # Compute the local exospheric temperature with the geomagnetic storm
    # effect.
    T∞ = Tl + ΔT∞

    # Temperature at desired altitude
    # ==========================================================================

    # Eq. 4-96 [3]
    # --------------------------------------------------------------------------
    #
    # Compute the temperature at inflection point `zx`.
    #
    # The values at [1, p. 369] are from an old version of Jacchia 1971 model.
    # We will use the new values available at [2].
    a  =  371.6678
    b  =    0.0518806
    c  = -294.3505
    d  =   -0.00216222
    Tx = a + b*T∞ + c*exp(d*T∞)

    # Compute the temperature at desired point.
    Tz = _jr1971_T(h, Tx, T∞)

    # Compute corrections to density, eqs. 4-96 to 4-101 [3]
    # ==========================================================================

    # If the user wants, then compute the corrections of geomegnetic effect,
    # semi-annual variation, and seasonal latitudinal variation. Otherwise, this
    # function will return the standard density [4].

    Δlog₁₀ρ_g  = 0.0
    Δlog₁₀ρ_sa = 0.0
    Δlog₁₀ρ_lt = 0.0
    Δlog₁₀ρ_c  = 0.0

    # Geomagnetic effect, eq. 4-97 [3]
    # --------------------------------------------------------------------------
    (h < 200) && ( Δlog₁₀ρ_g = 0.012Kp + 1.2e-5exp(Kp) )

    # Semi-annual variation, eqs. 4-98 to 4-99 [3]
    # --------------------------------------------------------------------------
    Φ    = (JD - 2436204.5)/365.2422  # Number of days since January 1, 1958.

    τ_sa = Φ + 0.09544( ( 1/2 + 1/2*sin( 2π*Φ + 6.035 ) )^(1.65) - 1/2 )
    f_z  = ( 5.876e-7h^2.331 + 0.06328 )*exp( -0.002868h )
    g_t  = 0.02835 + ( 0.3817 + 0.17829sin( 2π*τ_sa + 4.137 ) )*
    sin( 4π*τ_sa + 4.259 )

    Δlog₁₀ρ_sa = f_z*g_t

    # Seasonal latitudinal variation, eq. 4-100 [3]
    # --------------------------------------------------------------------------
    Δlog₁₀ρ_lt = 0.014( h - 90 )*exp( -0.0013( h - 90 )^2 )*
    sin( 2π*Φ + 1.72 )*sin( glat )*abs( sin(glat) )

    # Total correction, eq. B-10 [4]
    # --------------------------------------------------------------------------
    Δlog₁₀ρ_c = Δlog₁₀ρ_g + Δlog₁₀ρ_lt + Δlog₁₀ρ_sa

    # Density
    # ==========================================================================

    if h == z₁
        ρ = ρ₁ * 10^Δlog₁₀ρ_c

        # Convert to SI and return.
        return JR1971_Output(nN2   = (ρ * μi.N₂) * Av/Mi.N₂ * 1e6,
                             nO2   = (ρ * μi.O₂) * Av/Mi.O₂ * 1e6,
                             nO    = (ρ * μi.O)  * Av/Mi.O  * 1e6,
                             nAr   = (ρ * μi.Ar) * Av/Mi.Ar * 1e6,
                             nHe   = (ρ * μi.He) * Av/Mi.He * 1e6,
                             nH    = (ρ * μi.H)  * Av/Mi.H  * 1e6,
                             rho   = ρ*1e3,
                             T_exo = T∞,
                             Tz    = Tz )

    elseif z₁ < h <= zx

        # First, we need to find the roots of the polynomial:
        #
        #   P(Z) = c₀ + c₁ ⋅ z + c₂ ⋅ z² + c₃ ⋅ z³ + c₄ ⋅ z⁴
        c₀ = ( 35^4*Tx/(Tx - T₁) + Ca[1] )/Ca[5]
        c₁ = Ca[2]/Ca[5]
        c₂ = Ca[3]/Ca[5]
        c₃ = Ca[4]/Ca[5]
        c₄ = Ca[5]/Ca[5]
        r₁, r₂, x, y = _jr1971_roots([c₀, c₁, c₂, c₃, c₄])

        # f and k, [1. p. 371]
        # ----------------------------------------------------------------------
        f = 35^4*Ra²/Ca[5]
        k = -g₀/( Rstar*(Tx - T₁) )

        # U(ν), V(ν), and W(ν) functions, [1, p. 372]
        # ----------------------------------------------------------------------

        # TODO: Create copies of the roots, which are the variables used in the
        # functions U(ν), V(ν), and W(ν). Otherwise, the code is not type
        # stable. This seems related to:
        #
        #   https://github.com/JuliaLang/julia/issues/15276
        #
        cr₁, cr₂, cx, cy = r₁, r₂, x, y

        U(ν) = (ν + Ra)^2*(ν^2 - 2cx*ν + cx^2 + cy^2)*(cr₁ - cr₂)
        V(ν) = (ν^2 - 2cx*ν + cx^2 + cy^2)*(ν - cr₁)*(ν - cr₂)

        # The equation for `W(ν)` in [1] was:
        #
        #   W(ν) = r₁*r₂*Ra²*(ν + Ra) + (x^2 + y^2)*Ra*(Ra*ν + r₁*r₂)
        #
        # However, [3,4] mention that this is not correct, and must be replaced
        # by:
        W(ν) = cr₁*cr₂*Ra*(Ra + ν)*( Ra + ( cx^2 + cy^2 )/ν )

        # X, [1, p. 372]
        # ------------------------------------------------------------------
        X = -2r₁*r₂*Ra*( Ra² + 2x*Ra + x^2 + y^2)

        # The original paper [1] divided into two sections: 90 km to 105 km, and
        # 105 km to 125 km. However, [3] divided into 90 to 100 km, and 100 km
        # to 125 km.

        if h <= z₂

            # Altitudes between 90 km and 100 km
            # ==================================================================

            # S(z) polynomial, [1, p. 371]
            # ------------------------------------------------------------------
            B₀, B₁, B₂, B₃, B₄, B₅ = [ α[k] + β[k]*Tx/(Tx - T₁) for k = 1:6 ]

            S(z) = @evalpoly(z, B₀, B₁, B₂, B₃, B₄, B₅)

            # Auxiliary variables, [1. p. 372]
            # ------------------------------------------------------------------
            p₂ =   S(r₁)/U(r₁)
            p₃ =  -S(r₂)/U(r₂)
            p₅ =  S(-Ra)/V(-Ra)

            # There is a typo in the fourth term in [1] that was corrected in
            # [3].
            p₄ = ( B₀ - r₁*r₂*Ra²*( B₄ + ( 2x + r₁ + r₂ - Ra )*B₅ ) -
                  r₁*r₂*Ra*( x^2 + y^2 )*B₅ + r₁*r₂*( Ra² - ( x^2 + y^2 ) )*p₅ +
                  W(r₁)*p₂ + W(r₂)*p₃ )/X

            p₆ = B₄ + ( 2x + r₁ + r₂ - Ra )*B₅ - p₅ - 2( x + Ra )*p₄ -
                 ( r₂ + Ra )*p₃ - (r₁ + Ra)*p₂

            p₁ = B₅ - 2p₄ - p₃ - p₂

            # F₁ and F₂, [1, p. 372-373]
            # ------------------------------------------------------------------
            log_F₁ = p₁ * log( ( h + Ra )/( z₁ + Ra ) ) +
                     p₂ * log( ( h - r₁ )/( z₁ - r₁ ) ) +
                     p₃ * log( ( h - r₂ )/( z₁ - r₂ ) ) +
                     p₄ * log( (  h^2 - 2x*h  + x^2 + y^2 )/
                               ( z₁^2 - 2x*z₁ + x^2 + y^2 ) )

            # This equation in [4] is wrong, since `f` is multiplying `A₆`. We
            # will use the one in [3].
            F₂ = ( h - z₁ )*( Aa[7] + p₅/( (h + Ra)*(z₁ + Ra) ) ) +
                 p₆/y*atan( y*(h - z₁) / ( y^2 + (h - x)*(z₁ - x) ) )

            # Compute the density, eq. 13 [1]
            # ------------------------------------------------------------------
            Mz = _jr1971_M(h)
            ρ  = ρ₁*Mz*T₁/( M₁*Tz )*exp( k*( log_F₁ + F₂ ) ) * 10^Δlog₁₀ρ_c

            # Convert to SI and return.
            return JR1971_Output(nN2   = (ρ * μi.N₂) * Av/Mi.N₂ * 1e6,
                                 nO2   = (ρ * μi.O₂) * Av/Mi.O₂ * 1e6,
                                 nO    = (ρ * μi.O)  * Av/Mi.O  * 1e6,
                                 nAr   = (ρ * μi.Ar) * Av/Mi.Ar * 1e6,
                                 nHe   = (ρ * μi.He) * Av/Mi.He * 1e6,
                                 nH    = (ρ * μi.H)  * Av/Mi.H  * 1e6,
                                 rho   = ρ*1e3,
                                 T_exo = T∞,
                                 Tz    = Tz )
        else

            # Altitudes between 100 km and 125 km
            # ==================================================================

            # First, we need to compute the temperature and density at 100 km.
            T_₁₀₀ = _jr1971_T(z₂, Tx, T∞)

            # References [3,4] suggest to compute the density using a polynomial
            # fit, so that the computational burden can be reduced:
            #
            ρ_₁₀₀ = @evalpoly(T∞, ζ[1], ζ[2], ζ[3], ζ[4], ζ[5], ζ[6], ζ[7])*M₀

            # However, it turns out that this approach leads to discontinuity at
            # 100 km. This was also seen by GMAT [5].
            #
            # TODO: Check if we need to fix this.

            # Auxiliary variables, [1, p. 374]
            # ------------------------------------------------------------------

            q₂ =  1/U(r₁)
            q₃ = -1/U(r₂)
            q₅ =  1/V(-Ra)
            q₄ = (1 + r₁*r₂*(Ra² - (x^2 + y^2))*q₅ + W(r₁)*q₂ + W(r₂)*q₃)/X
            q₆ = -q₅ - 2(x + Ra)*q₄ - (r₂ + Ra)*q₃ - (r₁ + Ra)*q₂
            q₁ = -2q₄ - q₃ - q₂

            # F₃ and F₄, [1, p. 374]
            # ------------------------------------------------------------------
            log_F₃ = q₁ * log( ( h + Ra )/( z₂ + Ra ) ) +
                     q₂ * log( ( h - r₁ )/( z₂ - r₁ ) ) +
                     q₃ * log( ( h - r₂ )/( z₂ - r₂ ) ) +
                     q₄ * log( ( h^2 - 2x*h + x^2 + y^2 )/
                               ( z₂^2 - 2x*z₂ + x^2 + y^2 ) )

            F₄ = ( q₅*( h - z₂ ) )/( ( h + Ra )*( Ra + z₂ ) ) +
                 q₆/y*atan( y*(h - z₂)/( y^2 + (h - x)*(z₂ - x) ) )

            # Compute the density of each specie [3]
            # ------------------------------------------------------------------

            ρi = zeros(5)

            for i = 1:5
                ρi[i] = ρ_₁₀₀*Mi[i]/M₀*μi[i]*( T_₁₀₀/Tz )^(1 + αi[i])*
                        exp( Mi[i]*k*f*(log_F₃ + F₄ ) )
            end

            # Apply the corrections.
            ρi .*= 10^Δlog₁₀ρ_c

            # Convert to SI and return.
            return JR1971_Output(nN2   = ρi[1] * Av/Mi.N₂ * 1e6,
                                 nO2   = ρi[2] * Av/Mi.O₂ * 1e6,
                                 nO    = ρi[3] * Av/Mi.O  * 1e6,
                                 nAr   = ρi[4] * Av/Mi.Ar * 1e6,
                                 nHe   = ρi[5] * Av/Mi.He * 1e6,
                                 nH    = 0.0,
                                 rho   = sum(ρi)*1e3,
                                 T_exo = T∞,
                                 Tz    = Tz )
        end

    else

        # Altitudes higher than 125 km
        # ======================================================================

        # First, we need to compute the density at 125 km.

        # References [3,4] suggest to compute the density using a polynomial
        # fit, so that the computational burden can be reduced:

        ρi_₁₂₅ = zeros(5)

        for i = 1:5
            aux = @evalpoly(T∞, δij[i][1], δij[i][2], δij[i][3], δij[i][4],
                                δij[i][5], δij[i][6], δij[i][7])

            ρi_₁₂₅[i] = Mi[i]*10^aux/Av
        end

        # However, it turns out that this approach leads to discontinuity at 100
        # km. This was also seen by GMAT [5].
        #
        # TODO: Check if we need to fix this.

        # Compute `l` according to eq. 4-136 [3]
        # ----------------------------------------------------------------------
        l = @evalpoly(T∞, la[1], la[2], la[3], la[4], la[5])

        # Eq. 25' [1]
        # ----------------------------------------------------------------------
        γi = values(Mi).*(g₀*Ra²/( Rstar*l*T∞ )*(T∞ - Tx)/(Tx - T₁)*(zx - z₁)/(Ra + zx))

        # Eq. 25 [1]
        # ----------------------------------------------------------------------
        ρi = zeros(5)

        for i = 1:5
            ρi[i] = ρi_₁₂₅[i]*(Tx/Tz)^(1 + αi[i] + γi[i])* ( (T∞ - Tz)/
                                                             (T∞ - Tx) )^γi[i]
        end

        # Correction of seasonal variations of helium by latitude, eq. 4-101 [3]
        # ----------------------------------------------------------------------

        Δlog₁₀ρ_He = 0.65/(23.439291*π/180)*abs(δs)*
        ( sin( π/4 - glat*δs/( 2abs(δs) ) )^3 - 0.35355)

        ρi[_jr1971_id.He] *= 10^(Δlog₁₀ρ_He)

        # For altitude higher than 500 km, we must account for H.
        # ----------------------------------------------------------------------

        ρH = 0.0

        if h > 500
            # Compute the temperature and the H density at 500 km.
            T_₅₀₀       = _jr1971_T(500.0, Tx, T∞)
            log₁₀_T_₅₀₀ = log10(T_₅₀₀)
            ρH_₅₀₀      = Mi.H/Av*10^( 73.13 -
                                      ( 39.4 - 5.5log₁₀_T_₅₀₀ )*log₁₀_T_₅₀₀ )

            # Compute the H density at desired altitude.
            γH = Mi.H*g₀*Ra²/( Rstar*l*T∞ )*
                   (T∞ - Tx)/(Tx - T₁)*
                   (zx - z₁)/(Ra + zx)

            ρH = ρH_₅₀₀*(T_₅₀₀/Tz)^( 1 + αi.H + γH )*( (T∞ - Tz   )/
                                                       (T∞ - T_₅₀₀) )^γH
        end

        # Apply the corrections.
        ρi .*= 10^Δlog₁₀ρ_c
        ρH  *= 10^Δlog₁₀ρ_c

        # Convert to SI and return.
        return JR1971_Output(nN2   = ρi[1] * Av/Mi.N₂ * 1e6,
                             nO2   = ρi[2] * Av/Mi.O₂ * 1e6,
                             nO    = ρi[3] * Av/Mi.O  * 1e6,
                             nAr   = ρi[4] * Av/Mi.Ar * 1e6,
                             nHe   = ρi[5] * Av/Mi.He * 1e6,
                             nH    = ρH    * Av/Mi.H  * 1e6,
                             rho   = (sum(ρi) + ρH)*1e3,
                             T_exo = T∞,
                             Tz    = Tz )
    end
end


text:
"""
    SysConst{T}

Struct that contains all run time constants, e.g. lattice size, surface tension `γ` and so on.
"""
@with_kw struct SysConst{T}
    # Lattice
    Lx :: Int = 256
    Ly :: Int = 256
    Tmax :: Int = 1000
    tdump :: Int = Tmax÷10
    # Collision related
    τ :: T = 1.0
    cₛ :: T = 1/sqrt(3.0)
    μ :: T =  cₛ^2 / 2 *( τ - 0.5)
    # Force related
    δ :: T = 1.0
    kbt :: T = 0.0
    γ :: T = 0.01
    n :: Int = 9
    m :: Int = 3
    hmin :: T = 0.1
    hcrit :: T = 0.05
    g :: T = 0.0
end

"""
    SysConst_1D{T}

Struct that contains all run time constants, e.g. lattice size, surface tension `γ` and so on.
"""
@with_kw struct SysConst_1D{T}
    # Lattice
    L :: Int = 256
    Tmax :: Int = 1000
    tdump :: Int = Tmax÷10
    # Collision related
    τ :: T = 1.0
    cₛ :: T = 1/sqrt(3.0)
    μ :: T =  cₛ^2 / 2 *( τ - 0.5)
    # Force related
    δ :: T = 1.0
    kbt :: T = 0.0
    γ :: T = 0.01
    n :: Int = 9
    m :: Int = 3
    hmin :: T = 0.1
    hcrit :: T = 0.05
    g :: T = 0.0
end

"""
    Sys(sysc, device, exotic)

Mostly allocations of arrays used to run a simulation, but all within one function :)
"""
function Sys(sysc::SysConst, device::String, exotic::Bool, T)
    if device == "CPU"
        # Meso
        fout = zeros(sysc.Lx, sysc.Ly, 9)
        ftemp = zeros(sysc.Lx, sysc.Ly, 9)
        feq = zeros(sysc.Lx, sysc.Ly, 9)

        # Macro
        height = ones(sysc.Lx, sysc.Ly)
        velx = zeros(sysc.Lx, sysc.Ly)
        vely = zeros(sysc.Lx, sysc.Ly)
        vsq = zeros(sysc.Lx, sysc.Ly)
        pressure = zeros(sysc.Lx, sysc.Ly)
        dgrad = zeros(sysc.Lx, sysc.Ly, 8)

        # Forces
        Fx = zeros(sysc.Lx, sysc.Ly)
        Fy = zeros(sysc.Lx, sysc.Ly)
        slipx = zeros(sysc.Lx, sysc.Ly)
        slipy = zeros(sysc.Lx, sysc.Ly)
        h∇px = zeros(sysc.Lx, sysc.Ly)
        h∇py = zeros(sysc.Lx, sysc.Ly)
        if exotic
            # Probably many more forces here in the future
            fthermalx = zeros(sysc.Lx, sysc.Ly)
            fthermaly = zeros(sysc.Lx, sysc.Ly)
            return fout, ftemp, feq, height, velx, vely, vsq, pressure, dgrad, Fx, Fy, slipx, slipy, h∇px, h∇py, fthermalx, fthermaly
        end
        return fout, ftemp, feq, height, velx, vely, vsq, pressure, dgrad, Fx, Fy, slipx, slipy, h∇px, h∇py
    elseif device == "GPU"
        # Meso
        fout = CUDA.zeros(T, sysc.Lx, sysc.Ly, 9)
        ftemp = CUDA.zeros(T, sysc.Lx, sysc.Ly, 9)
        feq = CUDA.zeros(T, sysc.Lx, sysc.Ly, 9)

        # Macro
        height = CUDA.ones(T, sysc.Lx, sysc.Ly)
        velx = CUDA.zeros(T, sysc.Lx, sysc.Ly)
        vely = CUDA.zeros(T, sysc.Lx, sysc.Ly)
        vsq = CUDA.zeros(T, sysc.Lx, sysc.Ly)
        pressure = CUDA.zeros(T, sysc.Lx, sysc.Ly)
        dgrad = CUDA.zeros(T, sysc.Lx, sysc.Ly, 8)

        # Forces
        Fx = CUDA.zeros(T, sysc.Lx, sysc.Ly)
        Fy = CUDA.zeros(T, sysc.Lx, sysc.Ly)
        slipx = CUDA.zeros(T, sysc.Lx, sysc.Ly)
        slipy = CUDA.zeros(T, sysc.Lx, sysc.Ly)
        h∇px = CUDA.zeros(T, sysc.Lx, sysc.Ly)
        h∇py = CUDA.zeros(T, sysc.Lx, sysc.Ly)
        if exotic
            fthermalx = CUDA.zeros(T, sysc.Lx, sysc.Ly)
            fthermaly = CUDA.zeros(T, sysc.Lx, sysc.Ly)
            return fout, ftemp, feq, height, velx, vely, vsq, pressure, dgrad, Fx, Fy, slipx, slipy, h∇px, h∇py, fthermalx, fthermaly
        end
        return fout, ftemp, feq, height, velx, vely, vsq, pressure, dgrad, Fx, Fy, slipx, slipy, h∇px, h∇py
    end
end

function Sys(sysc::SysConst_1D, exotic::Bool, T)
    # Meso
    fout = zeros(sysc.L, 3)
    ftemp = zeros(sysc.L, 3)
    feq = zeros(sysc.L, 3)

    # Macro
    height = ones(sysc.L)
    vel = zeros(sysc.L)
    pressure = zeros(sysc.L)
    dgrad = zeros(sysc.L, 2)
    # Forces
    F = zeros(sysc.L)
    slip = zeros(sysc.L)
    h∇p = zeros(sysc.L)
    if exotic
        # Probably many more forces here in the future
        fthermal = zeros(sysc.L)

        return fout, ftemp, feq, height, vel, pressure, dgrad, F, slip, h∇p, fthermal
    end

    return fout, ftemp, feq, height, vel, pressure, dgrad, F, slip, h∇p

end


text:
@testset "Fitness" begin

    struct MockFitness <: AbstractFitness f end
    NaiveGAflux._fitness(s::MockFitness, c) = s.f

    struct MockCandidate{F} <: AbstractCandidate
        f::F
    end
    IdCand() = MockCandidate(identity)
    NaiveGAflux.model(c::MockCandidate) = c.f

    @testset "LogFitness" begin
        import NaiveGAflux: FittedCandidate
        iv = inputvertex("in", 2)
        c = CompGraph(iv, fluxvertex(Dense(2,3), iv)) |> CandidateModel
        lf = LogFitness(MockFitness(3))
        @test @test_logs (:info, " Candidate:   1\tvertices:   2\tparams:  0.01k\tfitness: 3") fitness(lf, c) == 3
        @test @test_logs (:info, " Candidate:   2\tvertices:   2\tparams:  0.01k\tfitness: 3") fitness(lf, c) == 3
        fc = FittedCandidate(1, 3, c)
        @test @test_logs (:info, " Candidate:   1\tvertices:   2\tparams:  0.01k\tfitness: 3") fitness(lf, fc) == 3
        @test @test_logs (:info, " Candidate:   2\tvertices:   2\tparams:  0.01k\tfitness: 3") fitness(lf, fc) == 3
    end

    @testset "GpuFitness" begin
        label = Val(:GpuTest)
        wasgpumapped = false
        function Flux.gpu(x::MockCandidate{Val{:GpuTest}})
            wasgpumapped = true
            return x
        end

        wascpumapped = false
        function NaiveGAflux.transferstate!(::MockCandidate{Val{:GpuTest}}, ::MockCandidate{Val{:GpuTest}})
            wascpumapped = true
        end

        @test fitness(GpuFitness(MockFitness(4)), MockCandidate(label)) == 4
        @test wasgpumapped
        @test wascpumapped
    end

    @testset "AccuracyFitness" begin
        struct DummyIter end
        Base.iterate(::DummyIter) = (([0 1 0 0; 1 0 0 1; 0 0 1 0], [1 0 0 0; 1 0 0 1; 0 0 1 0]), 1)
        Base.iterate(::DummyIter, state) = state == 1 ? (([1 0; 0 0; 0 1], [0 0; 1 1; 0 0]), 2) : nothing

        @test fitness(AccuracyFitness(DummyIter()), IdCand()) == 0.5
    end

    @testset "TrainThenFitness" begin
        iv = inputvertex("in", 2)
        ov = fluxvertex(Dense(2, 2), iv)
        cand = CandidateModel(CompGraph(iv, ov))

        ttf = TrainThenFitness(
            dataiter = [(ones(Float32, 2, 1), ones(Float32, 2, 1))],
            defaultloss=Flux.mse,
            defaultopt = Descent(0.0001),
            fitstrat = MockFitness(17),
            invalidfitness = 123.456)

        @test fitness(ttf, cand) === 17

        struct NaNCandidateModel <: AbstractCandidate end
        NaiveGAflux.ninputs(::NaNCandidateModel) = 1
        NaiveGAflux.model(f::NaNCandidateModel) = f
        (::NaNCandidateModel)(x) = x
        NaiveGAflux.lossfun(::NaNCandidateModel; default) = (args...) -> NaN32

        @test @test_logs (:warn, "NaN loss detected when training!") fitness(ttf, NaNCandidateModel()) == 123.456
    end

    @testset "TrainAccuracyFitness" begin
        dataiter  = [([0 1 0; 1 0 0; 0 0 1], [0 0 0; 1 1 0; 0 0 1])]
        iv = inputvertex("in", 3)
        idgraph = CompGraph(iv, iv)

        taf = TrainAccuracyFitness(;drop=0.5,dataiter, defaultloss = Flux.Losses.mse, defaultopt = Descent(0.001))

        # Note: First example will be dropped
        @test fitness(taf, CandidateModel(idgraph))  == 0.5
    end

    @testset "MapFitness" begin
        @test fitness(MapFitness(f -> 2f, MockFitness(3)), IdCand()) == 6
    end

    @testset "EwmaFitness" begin
        import NaiveGAflux: FittedCandidate

        @test fitness(EwmaFitness(MockFitness(10)), IdCand()) == 10
        @test fitness(EwmaFitness(MockFitness(20)), FittedCandidate(0, 10, IdCand())) == 15
        @test fitness(EwmaFitness(MockFitness(20)), FittedCandidate(0, 15, IdCand())) == 17.5
    end

    @testset "TimeFitness" begin

        function busysleep(t)
            t0 = time()
            # Busy wait to avoid yielding since this causes sporadic failures in CI
            while time() - t0 < t
                1+1
            end
        end
        tf = TimeFitness(MockFitness(13))
        NaiveGAflux._fitness(s::MockFitness, c::MockCandidate{typeof(busysleep)}) = (busysleep(0.1); s.f)
        ftime, f = fitness(tf, MockCandidate(busysleep))
        @test ftime ≈ 0.1 atol=0.01
        @test f == 13
    end

    @testset "SizeFitness" begin
        @test fitness(SizeFitness(), MockCandidate(Dense(2,3))) == 9
    end

    @testset "AggFitness" begin
        af = AggFitness(+, MockFitness(3), MockFitness(2))

        @test fitness(af, IdCand()) == 5
        @test_throws ErrorException AggFitness(+)
    end
end


text:
using StackOverflow
using Test

qHolder = getrecentquestionsfortag()
@test_nowarn qHolder[1].answer_count

@test_nowarn getrecentbadgerecipients()
@test_nowarn geterrors()

siteinfoholder = getsiteinfo()
@test_nowarn dump(siteinfoholder)

@test_nowarn searchtag("julia")
@test_nowarn getusers()
@test_nowarn getcomments()
@test_nowarn getsuggested_edits()


text:
#### z parsing.

function countz(Gs::Vector{NMRSpectraSimulator.zFIDModelType{T,SST}}) where {T,SST}

    N_z = 0
    N_z_singlets = 0
    for n = 1:length(Gs)
        # for i = 1:length(Gs[n].zs)
        #     #for l = 1:length(Gs[n].zs[i])
        #     N_z += length(Gs[n].zs[i])
        #
        #     #end
        # end
        N_z += sum(length.(Gs[n].zs))

        N_z_singlets += length(Gs[n].zs_singlets)
    end

    return N_z, N_z_singlets
end

function parsez!(Gs::Vector{NMRSpectraSimulator.zFIDModelType{T,SST}},
    z_vec::Vector{Complex{T}}) where {T,SST}

    @assert length(z_vec) == sum(countz(Gs))

    j = 0
    for n = 1:length(Gs)
        for i = 1:length(Gs[n].zs)
            for l = 1:length(Gs[n].zs[i])
                j += 1
                Gs[n].zs[i][l] = z_vec[j]
            end
        end

        for i = 1:length(Gs[n].zs_singlets)
            j += 1
            Gs[n].zs_singlets[i] = z_vec[j]
        end
    end

    return j
end

function resetz!(Gs::Vector{NMRSpectraSimulator.zFIDModelType{T,SST}}) where {T,SST}

    j = 0
    for n = 1:length(Gs)
        for i = 1:length(Gs[n].zs)
            for l = 1:length(Gs[n].zs[i])
                j += 1
                Gs[n].zs[i][l] = one(T)
            end
        end

        for i = 1:length(Gs[n].zs_singlets)
            j += 1
            Gs[n].zs_singlets[i] = one(T)
        end
    end

    return j
end

######### z, compensation for α at the (i,k) level. (spin system index, partition element index).
# include singlets too.

function updatez!(  A::Matrix{Complex{T}},
    b::Vector{T},
    z::Vector{Complex{T}},
    U_rad_LS,
    Gs::Vector{NMRSpectraSimulator.zFIDModelType{T,SST}},
    w::Vector{T}) where {T <: Real,SST}

    #@assert size(A) == (length(b), length(αS))
    @assert length(b) == 2*length(U_rad_LS)

    evaldesignmatrixz!(A, U_rad_LS, Gs, w)

    W = reinterpret(T, A)
    z[:] = W\b

    return nothing
end

function evaldesignmatrixz!(B::Matrix{Complex{T}},
    U_rad,
    Es::Vector{NMRSpectraSimulator.zFIDModelType{T,NMRSpectraSimulator.SpinSysParamsType1{T}}},
    w::Vector{T}) where T <: Real

    #
    M = length(U_rad)
    N = length(Es)

    #N_κ, N_κ_singlets = countκ(Es)
    #@assert size(B) == (M, N_κ + N_κ_singlets)
    #fill!(B, Inf) # debug.

    #resetκ!(Es)
    j = 0

    # loop over each κ partition element in Es.
    for n = 1:N
        A = Es[n]
        @assert length(A.zs) == length(A.core.qs)

        #spin system.
        #for i = 1:length(A.zs)
        for i in eachindex(A.zs)

            # partition
            for k in eachindex(A.zs[i])
                j += 1

                B[:,j] = collect( w[n]*A.core.gs[i][k](U_rad[m] - A.core.ss_params.d[i], A.core.ss_params.κs_λ[i]) for m = 1:M )

            end
        end

        # singlets.
        for k = 1:length(A.zs_singlets)
            j += 1

            B[:,j] = collect( w[n]*evalκsingletsz(U_rad[m], A.core.d_singlets,
                A.core.αs_singlets, A.core.Ωs_singlets,
                A.core.β_singlets, A.core.λ0, A.core.κs_λ_singlets) for m = 1:M )
        end

    end

    return nothing
end

# modified from NMRSpectraSimulator.evalκsinglets().
function evalκsingletsz(u_rad::T, d::Vector{T},
    αs_singlets::Vector{T}, Ωs_singlets,
    λ0::T, λ_multipliers::Vector{T})::Complex{T} where T <: Real

    out = zero(Complex{T})
    for i = 1:length(αs_singlets)
        τ = u_rad - d[i]

        λ = λ0*λ_multipliers[i]
        Ω = Ωs_singlets[i]
        #out += αs_singlets[i]*exp(im*βs_singlets[i])/(λ+im*(τ-Ω))
        out += αs_singlets[i]/(λ+im*(τ-Ω))
    end
    return out
end


text:
mutable struct MissingVector <: AbstractVector{Missing}
    len::Int
end

Base.IndexStyle(::Type{MissingVector}) = Base.IndexLinear()
Base.size(x::MissingVector) = (x.len,)

MissingVector(::UndefInitializer, dims::Dims) = MissingVector(dims[1])

Base.@propagate_inbounds function Base.getindex(x::MissingVector, i::Int)
    @boundscheck checkbounds(x, i)
    missing
end

Base.@propagate_inbounds function Base.setindex!(x::MissingVector, ::Missing, i::Int)
    @boundscheck checkbounds(x, i)
    missing
end

Base.similar(x::MissingVector, ::Type{Missing}, dims::Dims{N}) where {N} = MissingVector(dims[1])

function Base.empty!(x::MissingVector)
    x.len = 0
    return x
end

Base.empty(m::MissingVector, ::Missing) = MissingVector(0)

function Base.resize!(x::MissingVector, len)
    len >= 0 || throw(ArgumentError("`len` must be >= 0 when resizing MissingVector"))
    x.len = len
    return x
end

function Base.push!(x::MissingVector, ::Missing)
    x.len += 1
    return x
end

function Base.pushfirst!(x::MissingVector, ::Missing)
    x.len += 1
    return x
end

Base.@propagate_inbounds function Base.deleteat!(x::MissingVector, i::Integer)
    @boundscheck checkbounds(x, i)
    x.len -= 1
    return x
end

Base.@propagate_inbounds function Base.deleteat!(x::MissingVector, inds)
    for i in reverse(inds)
        deleteat!(x, i)
    end
    return x
end

Base.@propagate_inbounds function Base.deleteat!(x::MissingVector, inds::AbstractVector{Bool})
    length(inds) == length(x) || throw(BoundsError(x, inds))
    for i = length(x):-1:1
        if inds[i]
            deleteat!(x, i)
        end
    end
    return x
end

function Base.pop!(x::MissingVector)
    if isempty(x)
        throw(ArgumentError("array must be non-empty"))
    end
    x.len -= 1
    return missing
end

function Base.popfirst!(x::MissingVector)
    if isempty(x)
        throw(ArgumentError("array must be non-empty"))
    end
    x.len -= 1
    return missing
end

Base.@propagate_inbounds function Base.insert!(x::MissingVector, i::Integer, ::Missing)
    if !(i == 1 && isempty(x))
        @boundscheck checkbounds(x, i)
    end
    x.len += 1
    return x
end

function Base.vcat(x::MissingVector, arrays::MissingVector...)
    A = MissingVector(x.len)
    for array in arrays
        A.len += array.len
    end
    return A
end

for f in (:append!, :prepend!)
    @eval begin
        function Base.$f(x::MissingVector, y::AbstractVector{Missing})
            x.len += length(y)
            return x
        end

        function Base.$f(x::MissingVector, itr)
            for i in itr
                push!(x, i)
            end
            return x
        end
    end
end


text:
#=
The power set of a set is the set of all its subsets. Write a function that, given a set, generates its power set.

For example, given the set {1, 2, 3}, it should return {{}, {1}, {2}, {3}, {1, 2}, {1, 3}, {2, 3}, {1, 2, 3}}.

You may also use a list or array to represent a set.
=#

using Test
include("Solutions/problem37_generate_powerset.jl")

@test get_powerset([1,2,3]) == [[], [1], [2], [3], [1, 2], [1, 3], [2, 3], [1, 2, 3]]

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
# Octahedron Stein Distribution
#
# Represents a uniform distribution over an octahedron on the integer lattice
# in three dimensions, i.e., a scaled l1 ball in R^3 intersected with Z^3.
# The vertices of the octahedron are at (0,0,ulim), (0,0,-ulim), (0,ulim,0), etc.

using StatsBase: wsample

type SteinOctahedron <: SteinDistribution
    ulim::Int64
    c1::Array{Float64}
    c2::Array{Float64}
    c3::Array{Float64}
end

SteinOctahedron(ulim::Int64) = SteinOctahedron(
    ulim,
    [1.0 for i=1:3],
    [1.0 for i=1:3],
    [1.0 for i=1:3]
);

# Draw n independent samples from distribution via rejection sampling
function rand(d::SteinOctahedron, n::Int64)
    x = zeros(n,3);
    i = 0;
    while i < n
        y = rand(-d.ulim:d.ulim,3);
        if sum(abs(y)) <= d.ulim
            i += 1;
            x[i,:] = y;
        end
    end
    x
end

# Joint distribution evaluated at a single point (unnormalized)
function jointdistribution(d::SteinOctahedron, y::Array{Float64})
    sum(abs(y)) <= d.ulim ? 1.0 : 0.0
end

# Conditional distribution of jth coordinate of y, given all other coordinates,
# represented as a probability vector
function conddistribution(d::SteinOctahedron, y::Array{Float64,1}, j::Int64)
    conddist = [jointdistribution(d, [y[1:(j-1)], yj, y[(j+1):end]])
                for yj in -d.ulim:d.ulim]
    conddist / sum(conddist)
end

function supportlowerbound(d::SteinOctahedron, j::Int64)
    -d.ulim
end

function supportupperbound(d::SteinOctahedron, j::Int64)
    d.ulim
end

function numdimensions(d::SteinOctahedron)
    3
end


text:
int_rules_1_2_3_3 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.2.3.3*(d+e*x^n)^q*(a+b*x^n+c*x^(2*n))^p =#
    #= Antiderivative(((~d)+(~e')*(~x)^(~n'))^(~q')*((~a)+(~b')*(~x)^(~n')+(~c')*(~x)^(~n2'))^(~p),~x) :=   1/c^p*Antiderivative((d+e*x^n)^q*(b/2+c*x^n)^(2*p),x) <--  FreeQ([a,b,c,d,e,n,p,q],x) && EqQ(n2,2*n) && EqQ(b^2-4*a*c,0) &&  IntegerQ(p) =#
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~(n')) ^ ~(q') * (~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => ((~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p / (~d + ~e * (~x) ^ ~n) ^ (2 * ~p)) * Antiderivative((~d + ~e * (~x) ^ ~n) ^ (~q + 2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p, ~q], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (Not(IntegerQ(~p)) && EqQ(2 * ~c * ~d - ~b * ~e, 0))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~(n')) ^ ~(q') * (~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => ((~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ FracPart(~p) / ((~c) ^ IntPart(~p) * (~b / 2 + ~c * (~x) ^ ~n) ^ (2 * FracPart(~p)))) * Antiderivative((~d + ~e * (~x) ^ ~n) ^ ~q * (~b / 2 + ~c * (~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p, ~q], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~(q') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => Antiderivative((~x) ^ (~n * (2 * ~p + ~q)) * (~e + ~d * (~x) ^ -(~n)) ^ ~q * (~c + ~b * (~x) ^ -(~n) + ~a * (~x) ^ (-2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (IntegersQ(~p, ~q) && NegQ(~n)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~(q') * (~a + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => Antiderivative((~x) ^ (~n * (2 * ~p + ~q)) * (~e + ~d * (~x) ^ -(~n)) ^ ~q * (~c + ~a * (~x) ^ (-2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (IntegersQ(~p, ~q) && NegQ(~n)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~(q') * (~(a') + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => -(Subst(Antiderivative(((~d + ~e * (~x) ^ -(~n)) ^ ~q * (~a + ~b * (~x) ^ -(~n) + ~c * (~x) ^ (-2 * ~n)) ^ ~p) / (~x) ^ 2, ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p, ~q], ~x) && (EqQ(~n2, 2 * ~n) && ILtQ(~n, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~(q') * (~a + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => -(Subst(Antiderivative(((~d + ~e * (~x) ^ -(~n)) ^ ~q * (~a + ~c * (~x) ^ (-2 * ~n)) ^ ~p) / (~x) ^ 2, ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~c, ~d, ~e, ~p, ~q], ~x) && (EqQ(~n2, 2 * ~n) && ILtQ(~n, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~(q') * (~(a') + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => With([g = Denominator(~n)], g * Subst(Antiderivative((~x) ^ (g - 1) * (~d + ~e * (~x) ^ (g * ~n)) ^ ~q * (~a + ~b * (~x) ^ (g * ~n) + ~c * (~x) ^ (2 * g * ~n)) ^ ~p, ~x), ~x, (~x) ^ (1 / g))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p, ~q], ~x) && (EqQ(~n2, 2 * ~n) && FractionQ(~n))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~(q') * (~a + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => With([g = Denominator(~n)], g * Subst(Antiderivative((~x) ^ (g - 1) * (~d + ~e * (~x) ^ (g * ~n)) ^ ~q * (~a + ~c * (~x) ^ (2 * g * ~n)) ^ ~p, ~x), ~x, (~x) ^ (1 / g))) <-- FreeQ([~a, ~c, ~d, ~e, ~p, ~q], ~x) && (EqQ(~n2, 2 * ~n) && FractionQ(~n))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) * (~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => ((~b * ~e - ~d * ~c) * (~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1)) / (~b * ~c * ~n * (~p + 1) * (~x) ^ (2 * ~n * (~p + 1))) + (~e / ~c) * Antiderivative((~x) ^ -(~n) * (~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1), ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (Not(IntegerQ(~p)) && EqQ(~n * (2 * ~p + 1) + 1, 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) * (~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => (~e * (~x) ^ (-(~n) + 1) * (~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1)) / (~c * (~n * (2 * ~p + 1) + 1)) <-- FreeQ([~b, ~c, ~d, ~e, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (Not(IntegerQ(~p)) && (NeQ(~n * (2 * ~p + 1) + 1, 0) && EqQ(~b * ~e * (~n * ~p + 1) - ~c * ~d * (~n * (2 * ~p + 1) + 1), 0))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) * (~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => (~e * (~x) ^ (-(~n) + 1) * (~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1)) / (~c * (~n * (2 * ~p + 1) + 1)) - ((~b * ~e * (~n * ~p + 1) - ~c * ~d * (~n * (2 * ~p + 1) + 1)) / (~c * (~n * (2 * ~p + 1) + 1))) * Antiderivative((~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (Not(IntegerQ(~p)) && (NeQ(~n * (2 * ~p + 1) + 1, 0) && NeQ(~b * ~e * (~n * ~p + 1) - ~c * ~d * (~n * (2 * ~p + 1) + 1), 0))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~(q') * (~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => ((~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ FracPart(~p) / ((~x) ^ (~n * FracPart(~p)) * (~b + ~c * (~x) ^ ~n) ^ FracPart(~p))) * Antiderivative((~x) ^ (~n * ~p) * (~d + ~e * (~x) ^ ~n) ^ ~q * (~b + ~c * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~b, ~c, ~d, ~e, ~n, ~p, ~q], ~x) && (EqQ(~n2, 2 * ~n) && Not(IntegerQ(~p)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~(q') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~(p'), ~x) => Antiderivative((~d + ~e * (~x) ^ ~n) ^ (~p + ~q) * (~a / ~d + (~c / ~e) * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~q], ~x) && (EqQ(~n2, 2 * ~n) && (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) ^ ~n) ^ ~(q') * (~a + ~(c') * (~x) ^ ~n2) ^ ~(p'), ~x) => Antiderivative((~d + ~e * (~x) ^ ~n) ^ (~p + ~q) * (~a / ~d + (~c / ~e) * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~n, ~q], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && IntegerQ(~p)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~q * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => ((~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ FracPart(~p) / ((~d + ~e * (~x) ^ ~n) ^ FracPart(~p) * (~a / ~d + (~c * (~x) ^ ~n) / ~e) ^ FracPart(~p))) * Antiderivative((~d + ~e * (~x) ^ ~n) ^ (~p + ~q) * (~a / ~d + (~c / ~e) * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p, ~q], ~x) && (EqQ(~n2, 2 * ~n) && (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) ^ ~n) ^ ~q * (~a + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => ((~a + ~c * (~x) ^ (2 * ~n)) ^ FracPart(~p) / ((~d + ~e * (~x) ^ ~n) ^ FracPart(~p) * (~a / ~d + (~c * (~x) ^ ~n) / ~e) ^ FracPart(~p))) * Antiderivative((~d + ~e * (~x) ^ ~n) ^ (~p + ~q) * (~a / ~d + (~c / ~e) * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~n, ~p, ~q], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~(q') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2), ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ ~n) ^ ~q * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && IGtQ(~q, 0))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~(q') * (~a + ~(c') * (~x) ^ ~n2), ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ ~n) ^ ~q * (~a + ~c * (~x) ^ (2 * ~n)), ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && IGtQ(~q, 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~q * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2), ~x) => (-(((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2)) * ~x * (~d + ~e * (~x) ^ ~n) ^ (~q + 1)) / (~d * (~e) ^ 2 * ~n * (~q + 1)) + (1 / (~n * (~q + 1) * ~d * (~e) ^ 2)) * Antiderivative((~d + ~e * (~x) ^ ~n) ^ (~q + 1) * Simp((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2 * (~n * (~q + 1) + 1) + ~c * ~d * ~e * ~n * (~q + 1) * (~x) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && LtQ(~q, -1))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~q * (~a + ~(c') * (~x) ^ ~n2), ~x) => (-((~c * (~d) ^ 2 + ~a * (~e) ^ 2)) * ~x * (~d + ~e * (~x) ^ ~n) ^ (~q + 1)) / (~d * (~e) ^ 2 * ~n * (~q + 1)) + (1 / (~n * (~q + 1) * ~d * (~e) ^ 2)) * Antiderivative((~d + ~e * (~x) ^ ~n) ^ (~q + 1) * Simp(~c * (~d) ^ 2 + ~a * (~e) ^ 2 * (~n * (~q + 1) + 1) + ~c * ~d * ~e * ~n * (~q + 1) * (~x) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && LtQ(~q, -1)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~q * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2), ~x) => (~c * (~x) ^ (~n + 1) * (~d + ~e * (~x) ^ ~n) ^ (~q + 1)) / (~e * (~n * (~q + 2) + 1)) + (1 / (~e * (~n * (~q + 2) + 1))) * Antiderivative((~d + ~e * (~x) ^ ~n) ^ ~q * (~a * ~e * (~n * (~q + 2) + 1) - (~c * ~d * (~n + 1) - ~b * ~e * (~n * (~q + 2) + 1)) * (~x) ^ ~n), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~q], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~q * (~a + ~(c') * (~x) ^ ~n2), ~x) => (~c * (~x) ^ (~n + 1) * (~d + ~e * (~x) ^ ~n) ^ (~q + 1)) / (~e * (~n * (~q + 2) + 1)) + (1 / (~e * (~n * (~q + 2) + 1))) * Antiderivative((~d + ~e * (~x) ^ ~n) ^ ~q * (~a * ~e * (~n * (~q + 2) + 1) - ~c * ~d * (~n + 1) * (~x) ^ ~n), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~n, ~q], ~x) && (EqQ(~n2, 2 * ~n) && NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) / (~a + ~(c') * (~x) ^ ~n2), ~x) => With([q = Rt(2 * ~d * ~e, 2)], ((~e) ^ 2 / (2 * ~c)) * Antiderivative(1 / (~d + q * (~x) ^ (~n / 2) + ~e * (~x) ^ ~n), ~x) + ((~e) ^ 2 / (2 * ~c)) * Antiderivative(1 / ((~d - q * (~x) ^ (~n / 2)) + ~e * (~x) ^ ~n), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0) && (IGtQ(~n / 2, 0) && PosQ(~d * ~e))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) / (~a + ~(c') * (~x) ^ ~n2), ~x) => With([q = Rt(-2 * ~d * ~e, 2)], (~d / (2 * ~a)) * Antiderivative((~d - q * (~x) ^ (~n / 2)) / ((~d - q * (~x) ^ (~n / 2)) - ~e * (~x) ^ ~n), ~x) + (~d / (2 * ~a)) * Antiderivative((~d + q * (~x) ^ (~n / 2)) / ((~d + q * (~x) ^ (~n / 2)) - ~e * (~x) ^ ~n), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0) && (IGtQ(~n / 2, 0) && NegQ(~d * ~e))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) / (~a + ~(c') * (~x) ^ ~n2), ~x) => With([q = Rt(~a / ~c, 4)], (1 / (2 * sqrt(2) * ~c * q ^ 3)) * Antiderivative((sqrt(2) * ~d * q - (~d - ~e * q ^ 2) * (~x) ^ (~n / 2)) / ((q ^ 2 - sqrt(2) * q * (~x) ^ (~n / 2)) + (~x) ^ ~n), ~x) + (1 / (2 * sqrt(2) * ~c * q ^ 3)) * Antiderivative((sqrt(2) * ~d * q + (~d - ~e * q ^ 2) * (~x) ^ (~n / 2)) / (q ^ 2 + sqrt(2) * q * (~x) ^ (~n / 2) + (~x) ^ ~n), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (NeQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0) && (IGtQ(~n / 2, 0) && PosQ(~a * ~c)))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 3) / (~a + ~(c') * (~x) ^ 6), ~x) => With([q = Rt(~c / ~a, 6)], (1 / (3 * ~a * q ^ 2)) * Antiderivative((q ^ 2 * ~d - ~e * ~x) / (1 + q ^ 2 * (~x) ^ 2), ~x) + (1 / (6 * ~a * q ^ 2)) * Antiderivative((2 * q ^ 2 * ~d - (sqrt(3) * q ^ 3 * ~d - ~e) * ~x) / ((1 - sqrt(3) * q * ~x) + q ^ 2 * (~x) ^ 2), ~x) + (1 / (6 * ~a * q ^ 2)) * Antiderivative((2 * q ^ 2 * ~d + (sqrt(3) * q ^ 3 * ~d + ~e) * ~x) / (1 + sqrt(3) * q * ~x + q ^ 2 * (~x) ^ 2), ~x)) <-- FreeQ([~a, ~c, ~d, ~e], ~x) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && PosQ(~c / ~a))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) / (~a + ~(c') * (~x) ^ ~n2), ~x) => With([q = Rt(-(~a) / ~c, 2)], ((~d + ~e * q) / 2) * Antiderivative(1 / (~a + ~c * q * (~x) ^ ~n), ~x) + ((~d - ~e * q) / 2) * Antiderivative(1 / (~a - ~c * q * (~x) ^ ~n), ~x)) <-- FreeQ([~a, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (NegQ(~a * ~c) && IntegerQ(~n))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) / (~a + ~(c') * (~x) ^ ~n2), ~x) => ~d * Antiderivative(1 / (~a + ~c * (~x) ^ (2 * ~n)), ~x) + ~e * Antiderivative((~x) ^ ~n / (~a + ~c * (~x) ^ (2 * ~n)), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (PosQ(~a * ~c) || Not(IntegerQ(~n)))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) / (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2), ~x) => With([q = Rt((2 * ~d) / ~e - ~b / ~c, 2)], (~e / (2 * ~c)) * Antiderivative(1 / Simp(~d / ~e + q * (~x) ^ (~n / 2) + (~x) ^ ~n, ~x), ~x) + (~e / (2 * ~c)) * Antiderivative(1 / Simp((~d / ~e - q * (~x) ^ (~n / 2)) + (~x) ^ ~n, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (EqQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0) && (IGtQ(~n / 2, 0) && (GtQ((2 * ~d) / ~e - ~b / ~c, 0) || Not(LtQ((2 * ~d) / ~e - ~b / ~c, 0)) && EqQ(~d, ~e * Rt(~a / ~c, 2)))))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) / (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2), ~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) ^ ~n), ~x) + (~e / 2 - (2 * ~c * ~d - ~b * ~e) / (2q)) * Antiderivative(1 / (~b / 2 + q / 2 + ~c * (~x) ^ ~n), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (EqQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0) && (IGtQ(~n / 2, 0) && GtQ((~b) ^ 2 - 4 * ~a * ~c, 0)))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) / (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2), ~x) => With([q = Rt((-2 * ~d) / ~e - ~b / ~c, 2)], (~e / (2 * ~c * q)) * Antiderivative((q - 2 * (~x) ^ (~n / 2)) / Simp((~d / ~e + q * (~x) ^ (~n / 2)) - (~x) ^ ~n, ~x), ~x) + (~e / (2 * ~c * q)) * Antiderivative((q + 2 * (~x) ^ (~n / 2)) / Simp((~d / ~e - q * (~x) ^ (~n / 2)) - (~x) ^ ~n, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (EqQ(~c * (~d) ^ 2 - ~a * (~e) ^ 2, 0) && (IGtQ(~n / 2, 0) && Not(GtQ((~b) ^ 2 - 4 * ~a * ~c, 0))))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) / (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2), ~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) ^ ~n), ~x) + (~e / 2 - (2 * ~c * ~d - ~b * ~e) / (2q)) * Antiderivative(1 / (~b / 2 + q / 2 + ~c * (~x) ^ ~n), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (PosQ((~b) ^ 2 - 4 * ~a * ~c) || Not(IGtQ(~n / 2, 0))))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) / (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2), ~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) ^ (~n / 2)) / ((q - r * (~x) ^ (~n / 2)) + (~x) ^ ~n), ~x) + (1 / (2 * ~c * q * r)) * Antiderivative((~d * r + (~d - ~e * q) * (~x) ^ (~n / 2)) / (q + r * (~x) ^ (~n / 2) + (~x) ^ ~n), ~x))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (IGtQ(~n / 2, 0) && NegQ((~b) ^ 2 - 4 * ~a * ~c)))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~q / (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2), ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ ~n) ^ ~q / (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (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) ^ ~n) ^ ~q / (~a + ~(c') * (~x) ^ ~n2), ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ ~n) ^ ~q / (~a + ~c * (~x) ^ (2 * ~n)), ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && IntegerQ(~q)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~q / (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2), ~x) => ((~e) ^ 2 / ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2)) * Antiderivative((~d + ~e * (~x) ^ ~n) ^ ~q, ~x) + (1 / ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2)) * Antiderivative(((~d + ~e * (~x) ^ ~n) ^ (~q + 1) * ((~c * ~d - ~b * ~e) - ~c * ~e * (~x) ^ ~n)) / (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (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) ^ ~n) ^ ~q / (~a + ~(c') * (~x) ^ ~n2), ~x) => ((~e) ^ 2 / (~c * (~d) ^ 2 + ~a * (~e) ^ 2)) * Antiderivative((~d + ~e * (~x) ^ ~n) ^ ~q, ~x) + (~c / (~c * (~d) ^ 2 + ~a * (~e) ^ 2)) * Antiderivative(((~d + ~e * (~x) ^ ~n) ^ (~q + 1) * (~d - ~e * (~x) ^ ~n)) / (~a + ~c * (~x) ^ (2 * ~n)), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (Not(IntegerQ(~q)) && LtQ(~q, -1))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~q / (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2), ~x) => With([r = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((2 * ~c) / r) * Antiderivative((~d + ~e * (~x) ^ ~n) ^ ~q / ((~b - r) + 2 * ~c * (~x) ^ ~n), ~x) - ((2 * ~c) / r) * Antiderivative((~d + ~e * (~x) ^ ~n) ^ ~q / (~b + r + 2 * ~c * (~x) ^ ~n), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~q], ~x) && (EqQ(~n2, 2 * ~n) && (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) ^ ~n) ^ ~q / (~a + ~(c') * (~x) ^ ~n2), ~x) => With([r = Rt(-(~a) * ~c, 2)], (-(~c) / (2r)) * Antiderivative((~d + ~e * (~x) ^ ~n) ^ ~q / (r - ~c * (~x) ^ ~n), ~x) - (~c / (2r)) * Antiderivative((~d + ~e * (~x) ^ ~n) ^ ~q / (r + ~c * (~x) ^ ~n), ~x)) <-- FreeQ([~a, ~c, ~d, ~e, ~n, ~q], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && Not(IntegerQ(~q))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => (-(~x) * (((~d * (~b) ^ 2 - ~a * ~b * ~e) - 2 * ~a * ~c * ~d) + (~b * ~d - 2 * ~a * ~e) * ~c * (~x) ^ ~n) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1)) / (~a * ~n * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c)) + (1 / (~a * ~n * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c))) * Antiderivative(Simp((((~n * ~p + ~n + 1) * ~d * (~b) ^ 2 - ~a * ~b * ~e) - 2 * ~a * ~c * ~d * (2 * ~n * ~p + 2 * ~n + 1)) + (2 * ~n * ~p + 3 * ~n + 1) * (~d * ~b - 2 * ~a * ~e) * ~c * (~x) ^ ~n, ~x) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && ILtQ(~p, -1)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) * (~a + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => (-(~x) * (~d + ~e * (~x) ^ ~n) * (~a + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1)) / (2 * ~a * ~n * (~p + 1)) + (1 / (2 * ~a * ~n * (~p + 1))) * Antiderivative((~d * (2 * ~n * ~p + 2 * ~n + 1) + ~e * (2 * ~n * ~p + 3 * ~n + 1) * (~x) ^ ~n) * (~a + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && ILtQ(~p, -1))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ ~n) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ(~n2, 2 * ~n) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) * (~a + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ ~n) * (~a + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~n], ~x) && EqQ(~n2, 2 * ~n)
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~q * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => ((~c) ^ ~p * (~x) ^ ((2 * ~n * ~p - ~n) + 1) * (~d + ~e * (~x) ^ ~n) ^ (~q + 1)) / (~e * (2 * ~n * ~p + ~n * ~q + 1)) + Antiderivative((~d + ~e * (~x) ^ ~n) ^ ~q * ExpandToSum(((~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p - (~c) ^ ~p * (~x) ^ (2 * ~n * ~p)) - (~d * (~c) ^ ~p * ((2 * ~n * ~p - ~n) + 1) * (~x) ^ (2 * ~n * ~p - ~n)) / (~e * (2 * ~n * ~p + ~n * ~q + 1)), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~q], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~p, 0) && (NeQ(2 * ~n * ~p + ~n * ~q + 1, 0) && (IGtQ(~n, 0) && Not(IGtQ(~q, 0)))))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~q * (~a + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => ((~c) ^ ~p * (~x) ^ ((2 * ~n * ~p - ~n) + 1) * (~d + ~e * (~x) ^ ~n) ^ (~q + 1)) / (~e * (2 * ~n * ~p + ~n * ~q + 1)) + Antiderivative((~d + ~e * (~x) ^ ~n) ^ ~q * ExpandToSum(((~a + ~c * (~x) ^ (2 * ~n)) ^ ~p - (~c) ^ ~p * (~x) ^ (2 * ~n * ~p)) - (~d * (~c) ^ ~p * ((2 * ~n * ~p - ~n) + 1) * (~x) ^ (2 * ~n * ~p - ~n)) / (~e * (2 * ~n * ~p + ~n * ~q + 1)), ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~n, ~q], ~x) && (EqQ(~n2, 2 * ~n) && (IGtQ(~p, 0) && (NeQ(2 * ~n * ~p + ~n * ~q + 1, 0) && (IGtQ(~n, 0) && Not(IGtQ(~q, 0))))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~q * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ ~n) ^ ~q * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p, ~q], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (NeQ((~c * (~d) ^ 2 - ~b * ~d * ~e) + ~a * (~e) ^ 2, 0) && (IntegersQ(~p, ~q) && Not(IntegerQ(~n)) || (IGtQ(~p, 0) || IGtQ(~q, 0) && Not(IntegerQ(~n)))))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~q * (~a + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~d + ~e * (~x) ^ ~n) ^ ~q * (~a + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~n, ~p, ~q], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (IntegersQ(~p, ~q) && Not(IntegerQ(~n)) || (IGtQ(~p, 0) || IGtQ(~q, 0) && Not(IntegerQ(~n))))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~q * (~a + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~a + ~c * (~x) ^ (2 * ~n)) ^ ~p, (~d / ((~d) ^ 2 - (~e) ^ 2 * (~x) ^ (2 * ~n)) - (~e * (~x) ^ ~n) / ((~d) ^ 2 - (~e) ^ 2 * (~x) ^ (2 * ~n))) ^ -(~q), ~x), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (NeQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && (Not(IntegerQ(~p)) && ILtQ(~q, 0))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~q * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => Unintegrable((~d + ~e * (~x) ^ ~n) ^ ~q * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p, ~q], ~x) && EqQ(~n2, 2 * ~n)
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~q * (~a + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => Unintegrable((~d + ~e * (~x) ^ ~n) ^ ~q * (~a + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~n, ~p, ~q], ~x) && EqQ(~n2, 2 * ~n)
    @apply_utils Antiderivative((~d + ~(e') * (~u) ^ ~n) ^ ~(q') * (~a + ~(b') * (~u) ^ ~n + ~(c') * (~u) ^ ~n2) ^ ~(p'), ~x) => (1 / Coefficient(~u, ~x, 1)) * Subst(Antiderivative((~d + ~e * (~x) ^ ~n) ^ ~q * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x, ~u) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p, ~q], ~x) && (EqQ(~n2, 2 * ~n) && (LinearQ(~u, ~x) && NeQ(~u, ~x)))
    @apply_utils Antiderivative((~d + ~(e') * (~u) ^ ~n) ^ ~(q') * (~a + ~(c') * (~u) ^ ~n2) ^ ~(p'), ~x) => (1 / Coefficient(~u, ~x, 1)) * Subst(Antiderivative((~d + ~e * (~x) ^ ~n) ^ ~q * (~a + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x, ~u) <-- FreeQ([~a, ~c, ~d, ~e, ~n, ~p, ~q], ~x) && (EqQ(~n2, 2 * ~n) && (LinearQ(~u, ~x) && NeQ(~u, ~x)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~(mn')) ^ ~(q') * (~(a') + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => Antiderivative((~x) ^ (-(~n) * ~q) * (~e + ~d * (~x) ^ ~n) ^ ~q * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~mn, -(~n)) && (IntegerQ(~q) && (PosQ(~n) || Not(IntegerQ(~p))))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~(mn')) ^ ~(q') * (~a + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => Antiderivative((~x) ^ (~mn * ~q) * (~e + ~d * (~x) ^ -(~mn)) ^ ~q * (~a + ~c * (~x) ^ ~n2) ^ ~p, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~mn, ~p], ~x) && (EqQ(~n2, -2 * ~mn) && (IntegerQ(~q) && (PosQ(~n2) || Not(IntegerQ(~p)))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~(n')) ^ ~(q') * (~(a') + ~(b') * (~x) ^ ~(mn') + ~(c') * (~x) ^ ~(mn2')) ^ ~(p'), ~x) => Antiderivative((~x) ^ (-2 * ~n * ~p) * (~d + ~e * (~x) ^ ~n) ^ ~q * (~c + ~b * (~x) ^ ~n + ~a * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~q], ~x) && (EqQ(~mn, -(~n)) && (EqQ(~mn2, 2 * ~mn) && IntegerQ(~p)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~(n')) ^ ~(q') * (~(a') + ~(c') * (~x) ^ ~(mn2')) ^ ~(p'), ~x) => Antiderivative((~x) ^ (-2 * ~n * ~p) * (~d + ~e * (~x) ^ ~n) ^ ~q * (~c + ~a * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~n, ~q], ~x) && (EqQ(~mn2, -2 * ~n) && IntegerQ(~p))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~(mn')) ^ ~q * (~(a') + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => (((~e) ^ IntPart(~q) * (~x) ^ (~n * FracPart(~q)) * (~d + ~e * (~x) ^ -(~n)) ^ FracPart(~q)) / (1 + (~d * (~x) ^ ~n) / ~e) ^ FracPart(~q)) * Antiderivative((~x) ^ (-(~n) * ~q) * (1 + (~d * (~x) ^ ~n) / ~e) ^ ~q * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p, ~q], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~mn, -(~n)) && (Not(IntegerQ(~p)) && (Not(IntegerQ(~q)) && PosQ(~n)))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~(mn')) ^ ~q * (~a + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => (((~e) ^ IntPart(~q) * (~x) ^ (-(~mn) * FracPart(~q)) * (~d + ~e * (~x) ^ ~mn) ^ FracPart(~q)) / (1 + (~d * (~x) ^ -(~mn)) / ~e) ^ FracPart(~q)) * Antiderivative((~x) ^ (~mn * ~q) * (1 + (~d * (~x) ^ -(~mn)) / ~e) ^ ~q * (~a + ~c * (~x) ^ ~n2) ^ ~p, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~mn, ~p, ~q], ~x) && (EqQ(~n2, -2 * ~mn) && (Not(IntegerQ(~p)) && (Not(IntegerQ(~q)) && PosQ(~n2))))
    #= Antiderivative(((~d)+(~e')*(~x)^(~mn'))^(~q)*((~a')+(~b')*(~x)^(~n')+(~c')*(~x)^(~n2'))^(~p'),~x)  := x^(n*FracPart(q))*(d+e*x^(-n))^FracPart(q)/(e+d*x^n)^FracPart(q)* Antiderivative(x^(-n*q)*(e+d*x^n)^q*(a+b*x^n+c*x^(2*n))^p,x) <--  FreeQ([a,b,c,d,e,n,p,q],x) && EqQ(n2,2*n) && EqQ(mn,-n) &&  Not(IntegerQ(p)) && Not(IntegerQ(q)) && PosQ(n) =#
    #= Antiderivative(((~d)+(~e')*(~x)^(~mn'))^(~q)*((~a)+(~c')*(~x)^(~n2'))^(~p'),~x) :=  x^(-mn*FracPart(q))*(d+e*x^mn)^FracPart(q)/(e+d*x^(-mn))^FracPart(q)* Antiderivative(x^(mn*q)*(e+d*x^(-mn))^q*(a+c*x^n2)^p,x) <--  FreeQ([a,c,d,e,mn,p,q],x) && EqQ(n2,-2*mn) && Not(IntegerQ(p)) &&  Not(IntegerQ(q)) && PosQ(n2) =#
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~(n')) ^ ~(q') * (~(a') + ~(b') * (~x) ^ ~(mn') + ~(c') * (~x) ^ ~(mn2')) ^ ~p, ~x) => (((~x) ^ (2 * ~n * FracPart(~p)) * (~a + ~b * (~x) ^ -(~n) + ~c * (~x) ^ (-2 * ~n)) ^ FracPart(~p)) / (~c + ~b * (~x) ^ ~n + ~a * (~x) ^ (2 * ~n)) ^ FracPart(~p)) * Antiderivative((~x) ^ (-2 * ~n * ~p) * (~d + ~e * (~x) ^ ~n) ^ ~q * (~c + ~b * (~x) ^ ~n + ~a * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p, ~q], ~x) && (EqQ(~mn, -(~n)) && (EqQ(~mn2, 2 * ~mn) && (Not(IntegerQ(~p)) && (Not(IntegerQ(~q)) && PosQ(~n)))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~(n')) ^ ~(q') * (~(a') + ~(c') * (~x) ^ ~(mn2')) ^ ~p, ~x) => (((~x) ^ (2 * ~n * FracPart(~p)) * (~a + ~c * (~x) ^ (-2 * ~n)) ^ FracPart(~p)) / (~c + ~a * (~x) ^ (2 * ~n)) ^ FracPart(~p)) * Antiderivative((~x) ^ (-2 * ~n * ~p) * (~d + ~e * (~x) ^ ~n) ^ ~q * (~c + ~a * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~n, ~p, ~q], ~x) && (EqQ(~mn2, -2 * ~n) && (Not(IntegerQ(~p)) && (Not(IntegerQ(~q)) && PosQ(~n))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~(q') * (~a + ~(b') * (~x) ^ ~mn + ~(c') * (~x) ^ ~(n')) ^ ~(p'), ~x) => Antiderivative((~x) ^ (-(~n) * ~p) * (~d + ~e * (~x) ^ ~n) ^ ~q * (~b + ~a * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~q], ~x) && (EqQ(~mn, -(~n)) && IntegerQ(~p))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~(q') * (~a + ~(b') * (~x) ^ ~mn + ~(c') * (~x) ^ ~(n')) ^ ~(p'), ~x) => (((~x) ^ (~n * FracPart(~p)) * (~a + ~b / (~x) ^ ~n + ~c * (~x) ^ ~n) ^ FracPart(~p)) / (~b + ~a * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ FracPart(~p)) * Antiderivative((~x) ^ (-(~n) * ~p) * (~d + ~e * (~x) ^ ~n) ^ ~q * (~b + ~a * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p, ~q], ~x) && (EqQ(~mn, -(~n)) && Not(IntegerQ(~p)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~(q') * (~f + ~(g') * (~x) ^ ~n) ^ ~(r') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => ((~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ FracPart(~p) / ((4 * ~c) ^ IntPart(~p) * (~b + 2 * ~c * (~x) ^ ~n) ^ (2 * FracPart(~p)))) * Antiderivative((~d + ~e * (~x) ^ ~n) ^ ~q * (~f + ~g * (~x) ^ ~n) ^ ~r * (~b + 2 * ~c * (~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n, ~p, ~q, ~r], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~(q') * (~f + ~(g') * (~x) ^ ~n) ^ ~(r') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~(p'), ~x) => Antiderivative((~d + ~e * (~x) ^ ~n) ^ (~p + ~q) * (~f + ~g * (~x) ^ ~n) ^ ~r * (~a / ~d + (~c / ~e) * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n, ~q, ~r], ~x) && (EqQ(~n2, 2 * ~n) && (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) ^ ~n) ^ ~(q') * (~f + ~(g') * (~x) ^ ~n) ^ ~(r') * (~a + ~(c') * (~x) ^ ~n2) ^ ~(p'), ~x) => Antiderivative((~d + ~e * (~x) ^ ~n) ^ (~p + ~q) * (~f + ~g * (~x) ^ ~n) ^ ~r * (~a / ~d + (~c / ~e) * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g, ~n, ~q, ~r], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && IntegerQ(~p)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ ~n) ^ ~(q') * (~f + ~(g') * (~x) ^ ~n) ^ ~(r') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => ((~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ FracPart(~p) / ((~d + ~e * (~x) ^ ~n) ^ FracPart(~p) * (~a / ~d + (~c * (~x) ^ ~n) / ~e) ^ FracPart(~p))) * Antiderivative((~d + ~e * (~x) ^ ~n) ^ (~p + ~q) * (~f + ~g * (~x) ^ ~n) ^ ~r * (~a / ~d + (~c / ~e) * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~n, ~p, ~q, ~r], ~x) && (EqQ(~n2, 2 * ~n) && (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) ^ ~n) ^ ~(q') * (~f + ~(g') * (~x) ^ ~n) ^ ~(r') * (~a + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => ((~a + ~c * (~x) ^ (2 * ~n)) ^ FracPart(~p) / ((~d + ~e * (~x) ^ ~n) ^ FracPart(~p) * (~a / ~d + (~c * (~x) ^ ~n) / ~e) ^ FracPart(~p))) * Antiderivative((~d + ~e * (~x) ^ ~n) ^ (~p + ~q) * (~f + ~g * (~x) ^ ~n) ^ ~r * (~a / ~d + (~c / ~e) * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g, ~n, ~p, ~q, ~r], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~c * (~d) ^ 2 + ~a * (~e) ^ 2, 0) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative((~d1 + ~(e1') * (~x) ^ ~(non2')) ^ ~(q') * (~d2 + ~(e2') * (~x) ^ ~(non2')) ^ ~(q') * (~(a') + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~(p'), ~x) => Antiderivative((~d1 * ~d2 + ~e1 * ~e2 * (~x) ^ ~n) ^ ~q * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~n, ~p, ~q], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~non2, ~n / 2) && (EqQ(~d2 * ~e1 + ~d1 * ~e2, 0) && (IntegerQ(~q) || GtQ(~d1, 0) && GtQ(~d2, 0)))))
    @apply_utils Antiderivative((~d1 + ~(e1') * (~x) ^ ~(non2')) ^ ~(q') * (~d2 + ~(e2') * (~x) ^ ~(non2')) ^ ~(q') * (~(a') + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~(p'), ~x) => (((~d1 + ~e1 * (~x) ^ (~n / 2)) ^ FracPart(~q) * (~d2 + ~e2 * (~x) ^ (~n / 2)) ^ FracPart(~q)) / (~d1 * ~d2 + ~e1 * ~e2 * (~x) ^ ~n) ^ FracPart(~q)) * Antiderivative((~d1 * ~d2 + ~e1 * ~e2 * (~x) ^ ~n) ^ ~q * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~n, ~p, ~q], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~non2, ~n / 2) && EqQ(~d2 * ~e1 + ~d1 * ~e2, 0)))
    @apply_utils Antiderivative((~A + ~(B') * (~x) ^ ~(m')) * (~d + ~(e') * (~x) ^ ~n) ^ ~(q') * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~(p'), ~x) => ~A * Antiderivative((~d + ~e * (~x) ^ ~n) ^ ~q * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) + ~B * Antiderivative((~x) ^ ~m * (~d + ~e * (~x) ^ ~n) ^ ~q * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~A, ~B, ~m, ~n, ~p, ~q], ~x) && (EqQ(~n2, 2 * ~n) && EqQ((~m - ~n) + 1, 0))
    @apply_utils Antiderivative((~A + ~(B') * (~x) ^ ~(m')) * (~d + ~(e') * (~x) ^ ~n) ^ ~(q') * (~a + ~(c') * (~x) ^ ~n2) ^ ~(p'), ~x) => ~A * Antiderivative((~d + ~e * (~x) ^ ~n) ^ ~q * (~a + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) + ~B * Antiderivative((~x) ^ ~m * (~d + ~e * (~x) ^ ~n) ^ ~q * (~a + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~A, ~B, ~m, ~n, ~p, ~q], ~x) && (EqQ(~n2, 2 * ~n) && EqQ((~m - ~n) + 1, 0))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.2 Trinomial products/1.2.3 General/1.2.3.3 (d+e x^n)^q (a+b x^n+c x^(2 n))^p.jl


text:
using Knet: relu, gpu

@testset "layer" begin
    atype = gpu() >= 0 ? KnetArray{Float64} : Array{Float64}

    input = rand(1:10)
    output = rand(1:10)
    layer = AutoML.LinearLayer(input, output)

    batchsize = map(x->2^x, rand(3:10))
    x = atype(randn(input, batchsize))
    y = atype(randn(output, batchsize))

    output = layer(x)
    @test typeof(output) <: atype

    error = layer(x, y)
    @test typeof(error) == typeof(float(rand(1:2)))
end

@testset "linearchain" begin
    atype = gpu() >= 0 ? KnetArray{Float64} : Array{Float64}

    input = rand(1:10)
    hidden = rand(1:10)
    output = rand(1:10)

    layer = AutoML.LinearLayer(input, hidden)
    layer2 = AutoML.LinearLayer(hidden, output)
    chain = AutoML.LinearChain(layer, layer2)

    batchsize = map(x->2^x, rand(3:10))
    x = atype(randn(input, batchsize))
    y = atype(randn(output, batchsize))

    output = chain(x)
    @test typeof(output) <: atype

    error = chain(x, y)
    @test typeof(error) == typeof(float(rand(1:2)))
end

@testset "classificationchain" begin
    atype = gpu() >= 0 ? KnetArray{Float64} : Array{Float64}

    input = rand(1:10)
    hidden = rand(1:10)
    output = rand(1:10)

    layer = AutoML.LinearLayer(input, hidden)
    layer2 = AutoML.LinearLayer(hidden, output)
    chain = AutoML.CategoricalChain(layer, layer2)

    batchsize = map(x->2^x, rand(3:10))
    numclass = rand(2:10)

    x = atype(randn(input, batchsize))
    y = rand([1:numclass;], 1, batchsize)

    output = chain(x)
    @test typeof(output) <: atype

    error = chain(x, y)
    @test typeof(error) == typeof(float(rand(1:2)))
end

max_stars_repo_name: alperkesen/AutoML

max_issues_repo_name: alperkesen/AutoML


text:
module SDSS

include("io.jl")
include("db.jl")

end # module


text:
function lgd05()
#
#  test polynomial suggested by Goedecker
#   Legendre polynomial of degree 5
#
    p = lgd(5);
    z = [0;
         0.90617984593866;
         0.53846931010568;
         -0.90617984593866;
         -0.53846931010568
        ];
    z = [z ones(5)]
    p, PolyZeros(z)
end