text:
"""
Log of the InverseGamma PDF. Equivalent to using Distributions; logpdf(InverseGamma(α, β), x)
https://en.wikipedia.org/wiki/Inverse-gamma_distribution
"""
function log_inverse_gamma(x::Real, α::Real=1., β::Real=1.; d::Integer=0)
    @assert 0 <= d <= 2
    if d == 0
        x > 0 ? val = -(β / x) - (1 + α) * log(x) + α * log(β) - loggamma(α) : val = -Inf
    elseif d == 1
        x > 0 ? val = (β / x - (1 + α)) / x : val = 0
    else
        x > 0 ? val = (-2 * β / x + (1 + α)) / (x * x) : val = 0
    end
    return val
end

function gamma_mode_std_2_alpha_theta(m::Real, s::Real)
    θ = (sqrt(m ^ 2 + 4 * s ^ 2) - m) / 2
    α = m / θ + 1
    return [α, θ]
end

"""
Log of the Gamma PDF. Equivalent to using Distributions; logpdf(Gamma(α, β), x)
https://en.wikipedia.org/wiki/Gamma_distribution
"""
function log_gamma(x::Real, parameters::Vector{<:Real}; d::Integer=0, passed_mode_std::Bool=false)
    @assert 0 <= d <= 2
    @assert length(parameters) == 2
    assert_positive(parameters)
    if passed_mode_std
        parameters = gamma_mode_std_2_alpha_theta(parameters[1], parameters[2])
    end
    α = parameters[1]
    θ = parameters[2]
    if d == 0
        x > 0 ? val = -(x / θ) + (α - 1) * log(x) - α * log(θ) - loggamma(α) : val = -Inf
    elseif d == 1
        x > 0 ? val = (α - 1) / x - 1 / θ : val = 0
    else
        x > 0 ? val = -(α - 1) / (x * x) : val = 0
    end
    return val
end

gauss_cdf(x::Real) = (1 + erf(x))/2

"log of the Gaussian PDF. Equivalent to using Distributions; logpdf(Gaussian(μ, σ), x)"
function log_gaussian(x::Real, parameters::Vector{<:Real}; d::Integer=0, min::Real=-Inf, max::Real=Inf)
    @assert 0 <= d <= 2
    @assert length(parameters) == 2
    μ = parameters[1]
    σ = parameters[2]
    assert_positive(σ)
    @assert min < max
    normalization = 1 - gauss_cdf(min - μ) - gauss_cdf(μ - max)
    if d == 0
        (min < x < max) ? val = -((x - μ)^2/(2 * σ * σ)) - log(sqrt(2 * π) * σ) - log(normalization) : val = -Inf
    elseif d == 1
        (min < x < max) ? val = -(x - μ)/(σ * σ) : val = 0
    else
        (min < x < max) ? val = -1 / (σ * σ) : val = 0
    end
    return val
end

"Log of the Uniform PDF."
function log_uniform(x::Real, min_max::Vector{<:Real}=[0,1]; d::Integer=0)
    @assert 0 <= d <= 2
    @assert length(min_max) == 2
    min, max = min_max
    @assert min < max
    if d == 0
        min <= x <= max ? -log(max - min) : -Inf
    else
        return 0
    end
end

"""
Log of the log-Uniform PDF.
Flattens out in log space starting at shift
Also known as a (modified in shifted case) Jeffrey's prior
"""
function log_loguniform(x::Real, min_max::Vector{<:Real}; d::Integer=0, shift::Real=0)
    @assert 0 <= d <= 2
    @assert length(min_max) == 2
    min, max = min_max
    @assert 0 < min + shift < max + shift
    xpshift = x + shift
    if d == 0
        min <= x <= max ? val = -log(xpshift) - log(log((max + shift)/(min + shift))) : val = -Inf
    elseif d == 1
        min <= x <= max ? val = -1 / xpshift : val = 0
    elseif d == 2
        min <= x <= max ? val = 1 / (xpshift * xpshift) : val = 0
    end
    return val
end

"""
Log of the Rayleigh PDF.
"""
function log_Rayleigh(x::Real, σ::Real; d::Integer=0)
    @assert 0 <= d <= 2
    if d == 0
        0 <= x <= 1 ? val = -(x * x / (2 * σ * σ)) + log(x / σ / σ) -  log(1 - exp(-1 / (2 * σ * σ))) : val = -Inf
    elseif d == 1
        0 <= x <= 1 ? val = 1 / x - x / σ / σ : val = 0
    elseif d == 2
        0 <= x <= 1 ? val = -1 / x / x - 1 / σ / σ : val = 0
    end
    return val
end

"""
Log of the 2D circle PDF
"""
function log_circle(x::Vector{<:Real}, min_max_r::Vector{<:Real}; d::Vector{<:Integer}=[0,0])
    @assert minimum(d) == 0
    @assert maximum(d) <= 2
    @assert sum(d) <= 2

    @assert length(x) == length(min_max_r) == length(d) == 2
    min_r, max_r = min_max_r

    # @assert min_r < sqrt(dot(x, x)) < max_r
    # @assert min_θ < atan(x[1], x[2]) < max_θ
    if all(d .== 0)
        min_r < sqrt(dot(x, x)) < max_r ? val = -log(2 * π * (max_r ^ 2 - min_r ^ 2)) : val = -Inf
    else
        return 0
    end
    return val
end

"""
Log of the 2D unit cone PDF
"""
function log_cone(x::Vector{<:Real}; d::Vector{<:Integer}=[0,0])
    @assert minimum(d) >= 0
    @assert maximum(d) <= 2
    @assert sum(d) <= 2

    @assert length(x) == length(d) == 2

    r_sq = dot(x, x)  # x^2 + y^2
    r = sqrt(r_sq)
    if d == [0,0]
        0 <= r < 1 ? val = log(3 / π * (1 - r)) : val = -Inf
    elseif d == [0,1]
        0 <= r < 1 ? val = x[2] / (r_sq - r) : val = 0
    elseif d == [0,2]
        0 <= r < 1 ? val = (-x[2] ^ 2 * r + x[1] ^ 2 * (r - 1)) /
            (r ^ 3 * (1 - 2 * r + r_sq)) : val = 0
    elseif d == [1,0]
        0 <= r < 1 ? val = x[1] / (r_sq - r) : val = 0
    elseif d == [1,1]
        0 <= r < 1 ? val = (x[1] * x[2] * (1 - 2 * r)) /
            (r ^ 3 * (1 - 2 * r + r_sq)) : val = 0
    elseif d == [2,0]
        0 <= r < 1 ? val = (-x[1] ^ 2 * r + x[2] ^ 2 * (r - 1)) /
            (r ^ 3 * (1 - 2 * r + r_sq)) : val = 0
    end
    return val
end

"""
Log of the 2D unit quadratic cone PDF
"""
function log_quad_cone(x::Vector{<:Real}; d::Vector{<:Integer}=[0,0])
    @assert minimum(d) >= 0
    @assert maximum(d) <= 2
    @assert sum(d) <= 2

    @assert length(x) == length(d) == 2

    r_sq = dot(x, x)  # x^2 + y^2
    r = sqrt(r_sq)
    if d == [0,0]
        0 <= r < 1 ? val = log(6 / π * (1 - 2 * r + r_sq)) : val = -Inf
    elseif d == [0,1]
        0 <= r < 1 ? val = 2 * x[2] / (r_sq - r) : val = 0
    elseif d == [0,2]
        0 <= r < 1 ? val = 2 * (-x[2] ^ 2 * r + x[1] ^ 2 * (r - 1)) /
            (r ^ 3 * (1 - 2 * r + r_sq)) : val = 0
    elseif d == [1,0]
        0 <= r < 1 ? val = 2 * x[1] / (r_sq - r) : val = 0
    elseif d == [1,1]
        0 <= r < 1 ? val = (2 * x[1] * x[2] * (1 - 2 * r)) /
            (r ^ 3 * (1 - 2 * r + r_sq)) : val = 0
    elseif d == [2,0]
        0 <= r < 1 ? val = 2 * (-x[1] ^ 2 * r + x[2] ^ 2 * (r - 1)) /
            (r ^ 3 * (1 - 2 * r + r_sq)) : val = 0
    end
    return val
end

"""
Log of the 2D unit cubic cone PDF
"""
function log_cubic_cone(x::Vector{<:Real}; d::Vector{<:Integer}=[0,0])
    @assert minimum(d) >= 0
    @assert maximum(d) <= 2
    @assert sum(d) <= 2

    @assert length(x) == length(d) == 2

    r_sq = dot(x, x)  # x^2 + y^2
    r = sqrt(r_sq)
    if d == [0,0]
        0 <= r < 1 ? val = log(10 / π * (1 - r)^3) : val = -Inf
    elseif d == [0,1]
        0 <= r < 1 ? val = 3 * x[2] / (r_sq - r) : val = 0
    elseif d == [0,2]
        0 <= r < 1 ? val = 3 * (-x[2] ^ 2 * r + x[1] ^ 2 * (r - 1)) /
            (r ^ 3 * (1 - 2 * r + r_sq)) : val = 0
    elseif d == [1,0]
        0 <= r < 1 ? val = 3 * x[1] / (r_sq - r) : val = 0
    elseif d == [1,1]
        0 <= r < 1 ? val = (3 * x[1] * x[2] * (1 - 2 * r)) /
            (r ^ 3 * (1 - 2 * r + r_sq)) : val = 0
    elseif d == [2,0]
        0 <= r < 1 ? val = 3 * (-x[1] ^ 2 * r + x[2] ^ 2 * (r - 1)) /
            (r ^ 3 * (1 - 2 * r + r_sq)) : val = 0
    end
    return val
end

"""
Log of the 2D rotated Rayleigh PDF that is cutoff at r=1
ONLY ROUGHLY NORMALIZED according to σ = 1/5
"""
function log_rot_Rayleigh(x::Vector{<:Real}; d::Vector{<:Integer}=[0,0], σ=1/5)
    @assert minimum(d) >= 0
    @assert maximum(d) <= 2
    @assert sum(d) <= 2

    @assert length(x) == length(d) == 2
    r_sq = dot(x, x)  # x^2 + y^2
    r = sqrt(r_sq)
    σ_sq = σ ^ 2
    log_norm = -2 * log(σ) - 0.454215
    if d == [0,0]
        0 <= r < 1 ? val = -r_sq / (2 * σ_sq) + log(r) + log_norm : val = -Inf
    elseif d == [0,1]
        0 <= r < 1 ? val = -x[2] * (r_sq - σ_sq) / (r_sq * σ_sq) : val = 0
    elseif d == [0,2]
        0 <= r < 1 ? val = -(x[1] ^ 4 + x[1] ^ 2 * (2 * x[2] ^ 2 - σ_sq) + x[2] ^ 2 * (x[2] ^ 2 + σ_sq)) / (r_sq ^ 2 * σ_sq) : val = 0
    elseif d == [1,0]
        0 <= r < 1 ? val = -x[1] * (r_sq - σ_sq) / (r_sq * σ_sq) : val = 0
    elseif d == [1,1]
        0 <= r < 1 ? val = -2 * x[1] * x[2] / r_sq ^ 2 : val = 0
    elseif d == [2,0]
        0 <= r < 1 ? val = -(x[2] ^ 4 + x[2] ^ 2 * (2 * x[1] ^ 2 - σ_sq) + x[1] ^ 2 * (x[1] ^ 2 + σ_sq)) / (r_sq ^ 2 * σ_sq) : val = 0
    end
    return val
end

"""
Log of the bivariate normal PDF
NOTE THAT THAT WHEN USING lows!=[-∞,...], THIS IS NOT PROPERLY NORMALIZED
"""
function log_bvnormal(x::Vector{T}, Σ::Cholesky{T,Matrix{T}}; μ::Vector{T}=zeros(T, length(x)), d::Vector{<:Integer}=[0,0], lows::Vector{T}=zeros(T, length(x)) .- Inf) where {T<:Real}
    @assert minimum(d) >= 0
    @assert maximum(d) <= 2
    @assert sum(d) <= 2

    @assert length(x) == length(d) == length(μ) == 2
    y = x - μ

    if sum(d) == 0
        all(x .>= lows) ? val = -nlogL(Σ, y) : val = -Inf
    elseif sum(d) == 1
        y1 = Float64.(d)
        all(x .>= lows) ? val = -dnlogLdθ(y1, Σ \ y) : val = 0
    elseif d == [0,2]
        y1 = y2 = [0, 1.]
        all(x .>= lows) ? val = -d2nlogLdθ(y2, [0,0.], Σ \ y, Σ \ y1) : val = 0
    elseif d == [1,1]
        y1 = [1, 0.]
        y2 = [0, 1.]
        all(x .>= lows) ? val = -d2nlogLdθ(y2, [0,0.], Σ \ y, Σ \ y1) : val = 0
    elseif d == [2,0]
        y1 = y2 = [1, 0.]
        all(x .>= lows) ? val = -d2nlogLdθ(y2, [0,0.], Σ \ y, Σ \ y1) : val = 0
    end
    return val
end
function bvnormal_covariance(σ11::Real, σ22::Real, ρ::Real)
    assert_positive(σ11, σ22)
    @assert 0 <= ρ <= 1
    v12 = σ11*σ22*ρ
    return ridge_chol([σ11^2 v12;v12 σ22^2])
end


text:
const Counts = Clonglong

ContiguousSubArray{T,N,P,I<:Tuple{Union{Base.Slice, AbstractUnitRange}, Vararg{Any}}} = SubArray{T,N,P,I,true}
ContiguousArray{T} = Union{Vector{T}, ContiguousSubArray{T}}

"""
	start_counters(events)

Start counting hardware events.  This function cannot be called if the counters have already been started.
"""
function start_counters(evtset::EventSet)
    nevts = length(evtset)
    if nevts == 0
        throw(ArgumentError("one or more PAPI.Events required"))
	end

    ncounters = num_counters()
    if nevts > ncounters
        throw(ArgumentError("number of PAPI.Events must be ≤ PAPI.num_counters(), got $nevts"))
	end

    @papichk ccall((:PAPI_start, libpapi), Cint, (Cint, ), evtset)
    evtset
end

start_counters(evts::Vector{Event}) = start_counters(EventSet(evts))
start_counters(evts::Event...) = start_counters(collect(Event, evts))

"""
	reset_counters!(evtset::EventSet)

Reset counters and leaves them running after the call.
"""
function reset_counters!(evtset::EventSet)
	@papichk ccall((:PAPI_reset, libpapi), Cint, (Cint,), evtset)
	evtset
end

"""
	read_counters!(evtset::EventSet, values::Vector{Counts})

Read and reset counters.
`read_counters!` copies the event counters into values. The counters are reset and left running after the call.

The user must provide a vector of the correct size (equal to the number of events)
"""
function read_counters!(evtset::EventSet, values::ContiguousArray{Counts})
	@assert length(evtset) == length(values)
	@papichk ccall((:PAPI_read, libpapi), Cint, (Cint, Ptr{Counts}), evtset, values)
	@papichk ccall((:PAPI_reset, libpapi), Cint, (Cint,), evtset)
	values
end

"""
	accum_counters!(evtset::EventSet, values::Vector{Counts})

Accumulate and reset counters.
`accum_counters!` accumulates the event counters into values. The counters are reset and left running after the call.

The user must provide a vector of the correct size (equal to the number of events)
"""
function accum_counters!(evtset::EventSet, values::ContiguousArray{Counts})
	@assert length(evtset) == length(values)
  @papichk ccall((:PAPI_accum, libpapi), Cint, (Cint, Ptr{Counts}), evtset, values)
  values
end

"""
	stop_counters!(evtset::EventSet, values::Vector{Counts})

Stop counters and return current counts.
The counters must have been started by a previous call to `start_counters`

The user must provide a vector of the correct size (equal to the number of events)
"""
function stop_counters!(evtset::EventSet, values::ContiguousArray{Counts})
	@assert length(evtset) == length(values)
	@papichk ccall((:PAPI_stop, libpapi), Cint, (Cint, Ptr{Counts}), evtset, values)
	values
end

"""
	stop_counters(evtset::EventSet)

Stop counters and returns counts
The counters must have been started by a previous call to `start_counters`
"""
function stop_counters(evtset::EventSet)
	values = Vector{Counts}(undef, length(evtset))
	stop_counters!(evtset, values)
end

ext: jl


text:
"""
`Residual`: perform residual sampling

### Arguments
`n` - number of trials
`ω` - probabilities
"""
function Residual(n,ω)

    m = length(ω);

    R = sum(floor.(Int, n * ω));

    if (R < n)
        ω̄ = (n * ω .- floor.(n * ω)) / (n - R);
        N̄vals = rand(Multinomial(n-R, ω̄));
        Nvals = floor.(Int, n * ω) .+ N̄vals;
    else
        Nvals = floor.(Int, n * ω);
    end

    return Nvals

end

"""
`stratified`: perform stratified sampling

### Arguments
`n` - number of trials
`ω` - probabilities
"""
function stratified(n,ω)
    U = range(0,stop=n-1)/n .+ rand(n)/n;

    Nvals = counts(quantile.(Categorical(ω), U), 1:length(ω));

    return Nvals
end

"""
`systematic`: perform systematic sampling

### Arguments
`n` - number of trials
`ω` - probabilities
"""
function systematic(n,ω)

    U = range(0,stop=n-1)/n .+ rand()/n;

    Nvals = counts(quantile.(Categorical(ω), U), 1:length(ω));

    return Nvals
end

"""
`multinomial`: perform multinomial sampling

### Arguments
`n` - number of trials
`ω` - probabilities
"""
function multinomial(n, ω)

    multinomial = Multinomial(n, ω);
    Nvals = rand(multinomial);

    return Nvals
end


text:
include("helpers.jl")

function laff_axpy(alpha::Number, x::T, y::T) where T <: Matrix{Number}

    my, ny = verifydimensions(x, y)

    # Redimensionar o x para o formato de y
    x = reshape(x, my, ny)

    x * alpha .+ y
end

## Testar ##
x = Matrix([1,2,3], 3, 1)
y = Matrix([4,5,6], 1, 3)

laff_axpy(0.5, x, y)

max_stars_repo_path: LAFF/Semana 1/laff_axpy.jl

max_issues_repo_path: LAFF/Semana 1/laff_axpy.jl

max_forks_repo_path: LAFF/Semana 1/laff_axpy.jl


text:
module LegHemodynamics

importall NumericalIntegration
importall MAT
importall Dierckx
importall CSV
importall Missings
importall Distributions

# conversion factors and solver parameters
include("constants.jl")
include("solverparams.jl")

# file loader
include("loadtexttree.jl")

# timer
include("buildtimer.jl")

# type definitions
include("buildbranches.jl")
include("buildall.jl")

# initialization/discretization
include("calcbranchprops.jl")
include("discretizebranches.jl")
include("assignterminals.jl")
include("discretizeperiphery.jl")

# initial conditions for each submodel
include("applybranchics.jl")
include("applyperipheryics.jl")
include("updatediscretization.jl")

# cardiac cycle counter
include("setnumbeats.jl")

# TVD RK3 + WENO3 for artery interiors
include("tvdrk3.jl")
include("invariants.jl")
include("invariantodes.jl")
include("weno3.jl")
include("smoothinds.jl")
include("weights.jl")
include("reconstruct.jl")
include("F1d.jl")
include("J1d.jl")
include("abseigs.jl")

# 0D-1D coupling and 0D updates
include("interiorbcs.jl")
include("coupleproximal.jl")
include("coupledistal.jl")
include("coupling.jl")
include("fdist.jl")
include("Jdist.jl")
include("newtondist.jl")
include("linedist.jl")

# 1D interior junction updates
include("solvesplits.jl")
include("newton.jl")
include("fsingle.jl")
include("fdouble.jl")
include("ftriple.jl")
include("Jsingle.jl")
include("Jdouble.jl")
include("Jtriple.jl")

# 1D pressure updates
include("arterialpressure.jl")

# data assimilation
include("advancetime.jl")
include("builderrors.jl")
include("paramwalk.jl")
include("applycustomics.jl")
include("rediscretizet.jl")

# main
include("main.jl")

end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
using Tranformers.Basic.Embed
using Transformers

abstract type abstractSEM end
struct SEM <: abstractSEM
    tok_embed1::Embed
    tok_embed2::Dense
    pos_embed::PositionEmbedding
    seg_embed::Embed
end

@functor SEM

"just a wrapper for two dense layer."
SEM(vocab_size::Int, embed::Int,hs::Int, seg::Int,maxlen::int) = SEM(
    Embed(vocab_size,embed),
    Dense(embed, hs),
    PositionEmbedding(hs,maxlen)
    Embed(seg,hidden)

)

function (s::SEM)(x::AbstractArray{T, N},seg::AbstractArray{T, N}) where {T, N}
    # size(x) == (dims, seq_len)
    e =s.tok_embed1(x)
    e=s.tok_embed2(e)
    e=e .+ s.postionEmbedding(x) .+ s.seg_embed(seg)

end


text:
using HTTP
using JSON

const ROUTER = HTTP.Router()

function surface(req::HTTP.Request)
    headers = [
        "Access-Control-Allow-Origin" => "*",
        "Access-Control-Allow-Methods" => "GET, OPTIONS"
    ]
    if HTTP.method(req) == "OPTIONS"
        return HTTP.Response(200, headers)
    end

    uri = HTTP.URIs.splitpath(req.target)
    if length(uri) != 5
        return HTTP.Response(400, headers)
    end

    (offset, step, stop) = parse.(Float64, uri[3:5])
    r = offset:step:stop

    return HTTP.Response(200, headers; body=JSON.json([r, surface(r, σ=1, μ=0)]))
end

HTTP.@register(ROUTER, "GET", "/api/surface", surface)

function Base.run()
    @info "service started at http://127.0.0.1:8080/api/surface/-10/0.1/10"
    HTTP.serve(ROUTER, "127.0.0.1", 8080)
end

max_stars_repo_name: foldfelis/Replotter

max_issues_repo_name: foldfelis/Replotter

max_forks_repo_name: foldfelis/Replotter


text:
int_rules_1_2_3_6 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.2.3.6*P(x)*(d*x)^m*(a+b*x^n+c*x^(2*n))^p =#
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => (1 / ~n) * Subst(Antiderivative(SubstFor((~x) ^ ~n, ~Pq, ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, (~x) ^ ~n) && EqQ(Simplify((~m - ~n) + 1), 0)))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => Antiderivative(ExpandIntegrand((~d * ~x) ^ ~m * ~Pq * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, ~x) && IGtQ(~p, 0)))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~d + ~(e') * (~x) ^ ~(n') + ~(f') * (~x) ^ ~(n2')) * (~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => (~d * (~g * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1)) / (~a * ~g * (~m + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~a * ~e * (~m + 1) - ~b * ~d * (~m + ~n * (~p + 1) + 1), 0) && (EqQ(~a * ~f * (~m + 1) - ~c * ~d * (~m + 2 * ~n * (~p + 1) + 1), 0) && NeQ(~m, -1))))
    @apply_utils Antiderivative((~(g') * ~x) ^ ~(m') * (~d + ~(f') * (~x) ^ ~(n2')) * (~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => (~d * (~g * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1)) / (~a * ~g * (~m + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~f, ~g, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~m + ~n * (~p + 1) + 1, 0) && (EqQ(~c * ~d + ~a * ~f, 0) && NeQ(~m, -1))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * ~Pq * (~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 * ~x) ^ ~m * ~Pq * (~b + 2 * ~c * (~x) ^ ~n) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, ~x) && (EqQ((~b) ^ 2 - 4 * ~a * ~c, 0) && Not(IntegerQ(2 * ~p)))))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (Simplify((~m + 1) / ~n) - 1) * SubstFor((~x) ^ ~n, ~Pq, ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, (~x) ^ ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(Simplify((~m + 1) / ~n)))))
    @apply_utils Antiderivative((~d * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => ((~d * ~x) ^ ~m / (~x) ^ ~m) * Antiderivative((~x) ^ ~m * ~Pq * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, (~x) ^ ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IntegerQ(Simplify((~m + 1) / ~n)))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => (1 / ~d) * Antiderivative((~d * ~x) ^ (~m + 1) * PolynomialQuotient(~Pq, ~x, ~x) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, ~x) && EqQ(Coeff(~Pq, ~x, 0), 0)))
    @apply_utils Antiderivative(((~x) ^ ~(m') * (~e + ~(f') * (~x) ^ ~(q') + ~(g') * (~x) ^ ~(r') + ~(h') * (~x) ^ ~(s'))) / (~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')) ^ (3 / 2), ~x) => -((2 * ~c * (~b * ~f - 2 * ~a * ~g) + 2 * ~h * ((~b) ^ 2 - 4 * ~a * ~c) * (~x) ^ (~n / 2) + 2 * ~c * (2 * ~c * ~f - ~b * ~g) * (~x) ^ ~n)) / (~c * ~n * ((~b) ^ 2 - 4 * ~a * ~c) * sqrt(~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n))) <-- FreeQ([~a, ~b, ~c, ~e, ~f, ~g, ~h, ~m, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~q, ~n / 2) && (EqQ(~r, (3 * ~n) / 2) && (EqQ(~s, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (EqQ((2 * ~m - ~n) + 2, 0) && EqQ(~c * ~e + ~a * ~h, 0)))))))
    @apply_utils Antiderivative(((~d * ~x) ^ ~(m') * (~e + ~(f') * (~x) ^ ~(q') + ~(g') * (~x) ^ ~(r') + ~(h') * (~x) ^ ~(s'))) / (~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')) ^ (3 / 2), ~x) => ((~d * ~x) ^ ~m / (~x) ^ ~m) * Antiderivative(((~x) ^ ~m * (~e + ~f * (~x) ^ (~n / 2) + ~g * (~x) ^ ((3 * ~n) / 2) + ~h * (~x) ^ (2 * ~n))) / (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (3 / 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~q, ~n / 2) && (EqQ(~r, (3 * ~n) / 2) && (EqQ(~s, 2 * ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (EqQ((2 * ~m - ~n) + 2, 0) && EqQ(~c * ~e + ~a * ~h, 0)))))))
    @apply_utils Antiderivative((~x) ^ ~m * ~Pq * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => Module([q = Expon(~Pq, ~x)], Module([Q = PolynomialQuotient(~a * (~b * ~c) ^ (Floor((q - 1) / ~n) + 1) * (~x) ^ ~m * ~Pq, ~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n), ~x), R = PolynomialRemainder(~a * (~b * ~c) ^ (Floor((q - 1) / ~n) + 1) * (~x) ^ ~m * ~Pq, ~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n), ~x), i], ((-(~x) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1)) / ((~a) ^ 2 * ~n * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c) * (~b * ~c) ^ (Floor((q - 1) / ~n) + 1))) * Sum((((~b) ^ 2 - 2 * ~a * ~c) * Coeff(R, ~x, i) - ~a * ~b * Coeff(R, ~x, ~n + i)) * (~x) ^ i + ~c * (~b * Coeff(R, ~x, i) - 2 * ~a * Coeff(R, ~x, ~n + i)) * (~x) ^ (~n + i), [i, 0, ~n - 1]) + (1 / (~a * ~n * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c) * (~b * ~c) ^ (Floor((q - 1) / ~n) + 1))) * Antiderivative((~x) ^ ~m * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1) * ExpandToSum(~n * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c) * (~x) ^ -(~m) * Q + Sum(((((~b) ^ 2 * (~n * (~p + 1) + i + 1)) / ~a - 2 * ~c * (2 * ~n * (~p + 1) + i + 1)) * Coeff(R, ~x, i) - ~b * (i + 1) * Coeff(R, ~x, ~n + i)) * (~x) ^ (i - ~m) + ~c * (~n * (2 * ~p + 3) + i + 1) * ((~b / ~a) * Coeff(R, ~x, i) - 2 * Coeff(R, ~x, ~n + i)) * (~x) ^ ((~n + i) - ~m), [i, 0, ~n - 1]), ~x), ~x)) <-- GeQ(q, 2 * ~n)) <-- FreeQ([~a, ~b, ~c], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n, 0) && (LtQ(~p, -1) && ILtQ(~m, 0))))))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => With([g = GCD(~m + 1, ~n)], (1 / g) * Subst(Antiderivative((~x) ^ ((~m + 1) / g - 1) * ReplaceAll(~Pq, (~x->begin
                                                #= none:1 =#
                                                #= none:3 =#
                                                (~x) ^ (1 / g)
                                            end)) * (~a + ~b * (~x) ^ (~n / g) + ~c * (~x) ^ ((2 * ~n) / g)) ^ ~p, ~x), ~x, (~x) ^ g) <-- NeQ(g, 1)) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, (~x) ^ ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n, 0) && IntegerQ(~m)))))
    @apply_utils Antiderivative(((~(d') * ~x) ^ ~(m') * ~Pq) / (~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~n2), ~x) => Antiderivative(ExpandIntegrand(((~d * ~x) ^ ~m * ~Pq) / (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, (~x) ^ ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n, 0) && NiceSqrtQ((~b) ^ 2 - 4 * ~a * ~c)))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => With([q = Expon(~Pq, ~x)], With([Pqq = Coeff(~Pq, ~x, q)], (Pqq * (~d * ~x) ^ (((~m + q) - 2 * ~n) + 1) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ (~p + 1)) / (~c * (~d) ^ ((q - 2 * ~n) + 1) * (~m + q + 2 * ~n * ~p + 1)) + Antiderivative((~d * ~x) ^ ~m * ExpandToSum((~Pq - Pqq * (~x) ^ q) - (Pqq * (~a * (((~m + q) - 2 * ~n) + 1) * (~x) ^ (q - 2 * ~n) + ~b * (~m + q + ~n * (~p - 1) + 1) * (~x) ^ (q - ~n))) / (~c * (~m + q + 2 * ~n * ~p + 1)), ~x) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x)) <-- GeQ(q, 2 * ~n) && (NeQ(~m + q + 2 * ~n * ~p + 1, 0) && (IntegerQ(2 * ~p) || (EqQ(~n, 1) && IntegerQ(4 * ~p) || IntegerQ(~p + (q + 1) / (2 * ~n)))))) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, (~x) ^ ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IGtQ(~n, 0))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~n2) ^ ~p, ~x) => Module([q = Expon(~Pq, ~x), j, k], Antiderivative(Sum((1 / (~d) ^ j) * (~d * ~x) ^ (~m + j) * Sum(Coeff(~Pq, ~x, j + k * ~n) * (~x) ^ (k * ~n), [k, 0, (q - j) / ~n + 1]) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, [j, 0, ~n - 1]), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IGtQ(~n, 0) && Not(PolyQ(~Pq, (~x) ^ ~n))))))
    @apply_utils Antiderivative(((~(d') * ~x) ^ ~(m') * ~Pq) / (~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')), ~x) => Antiderivative(RationalFunctionExpand(((~d * ~x) ^ ~m * ~Pq) / (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && IGtQ(~n, 0))))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => With([q = Expon(~Pq, ~x)], -(Subst(Antiderivative((ExpandToSum((~x) ^ q * ReplaceAll(~Pq, (~x->begin
                                                        #= none:1 =#
                                                        #= none:3 =#
                                                        (~x) ^ -1
                                                    end)), ~x) * (~a + ~b * (~x) ^ -(~n) + ~c * (~x) ^ (-2 * ~n)) ^ ~p) / (~x) ^ (~m + q + 2), ~x), ~x, 1 / ~x))) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (ILtQ(~n, 0) && IntegerQ(~m)))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => With([g = Denominator(~m), q = Expon(~Pq, ~x)], (-g / ~d) * Subst(Antiderivative((ExpandToSum((~x) ^ (g * q) * ReplaceAll(~Pq, (~x->begin
                                                        #= none:1 =#
                                                        #= none:3 =#
                                                        (~d) ^ -1 * (~x) ^ -g
                                                    end)), ~x) * (~a + ~b * (~d) ^ -(~n) * (~x) ^ (-g * ~n) + ~c * (~d) ^ (-2 * ~n) * (~x) ^ (-2 * g * ~n)) ^ ~p) / (~x) ^ (g * (~m + q + 1) + 1), ~x), ~x, 1 / (~d * ~x) ^ (1 / g))) <-- FreeQ([~a, ~b, ~c, ~d, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (ILtQ(~n, 0) && FractionQ(~m)))))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~m * ~Pq * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => With([q = Expon(~Pq, ~x)], -((~d * ~x) ^ ~m) * ((~x) ^ -1) ^ ~m * Subst(Antiderivative((ExpandToSum((~x) ^ q * ReplaceAll(~Pq, (~x->begin
                                                        #= none:1 =#
                                                        #= none:3 =#
                                                        (~x) ^ -1
                                                    end)), ~x) * (~a + ~b * (~x) ^ -(~n) + ~c * (~x) ^ (-2 * ~n)) ^ ~p) / (~x) ^ (~m + q + 2), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (ILtQ(~n, 0) && Not(RationalQ(~m))))))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => With([g = Denominator(~n)], g * Subst(Antiderivative((~x) ^ (g * (~m + 1) - 1) * ReplaceAll(~Pq, (~x->begin
                                            #= none:1 =#
                                            #= none:3 =#
                                            (~x) ^ g
                                        end)) * (~a + ~b * (~x) ^ (g * ~n) + ~c * (~x) ^ (2 * g * ~n)) ^ ~p, ~x), ~x, (~x) ^ (1 / g))) <-- FreeQ([~a, ~b, ~c, ~m, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && FractionQ(~n))))
    @apply_utils Antiderivative((~d * ~x) ^ ~m * ~Pq * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => (((~d) ^ (~m - 1 / 2) * sqrt(~d * ~x)) / sqrt(~x)) * Antiderivative((~x) ^ ~m * ~Pq * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (FractionQ(~n) && IGtQ(~m + 1 / 2, 0)))))
    @apply_utils Antiderivative((~d * ~x) ^ ~m * ~Pq * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => (((~d) ^ (~m + 1 / 2) * sqrt(~x)) / sqrt(~d * ~x)) * Antiderivative((~x) ^ ~m * ~Pq * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (FractionQ(~n) && ILtQ(~m - 1 / 2, 0)))))
    @apply_utils Antiderivative((~d * ~x) ^ ~m * ~Pq * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => ((~d * ~x) ^ ~m / (~x) ^ ~m) * Antiderivative((~x) ^ ~m * ~Pq * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, ~x) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && FractionQ(~n))))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => (1 / (~m + 1)) * Subst(Antiderivative(ReplaceAll(SubstFor((~x) ^ ~n, ~Pq, ~x), (~x->begin
                                        #= none:1 =#
                                        #= none:3 =#
                                        (~x) ^ Simplify(~n / (~m + 1))
                                    end)) * (~a + ~b * (~x) ^ Simplify(~n / (~m + 1)) + ~c * (~x) ^ Simplify((2 * ~n) / (~m + 1))) ^ ~p, ~x), ~x, (~x) ^ (~m + 1)) <-- FreeQ([~a, ~b, ~c, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, (~x) ^ ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IntegerQ(Simplify(~n / (~m + 1))) && Not(IntegerQ(~n))))))
    @apply_utils Antiderivative((~d * ~x) ^ ~m * ~Pq * (~a + ~(b') * (~x) ^ ~n + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => ((~d * ~x) ^ ~m / (~x) ^ ~m) * Antiderivative((~x) ^ ~m * ~Pq * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, (~x) ^ ~n) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (IntegerQ(Simplify(~n / (~m + 1))) && Not(IntegerQ(~n))))))
    @apply_utils Antiderivative(((~(d') * ~x) ^ ~(m') * ~Pq) / (~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')), ~x) => With([q = Rt((~b) ^ 2 - 4 * ~a * ~c, 2)], ((2 * ~c) / q) * Antiderivative(((~d * ~x) ^ ~m * ~Pq) / ((~b - q) + 2 * ~c * (~x) ^ ~n), ~x) - ((2 * ~c) / q) * Antiderivative(((~d * ~x) ^ ~m * ~Pq) / (~b + q + 2 * ~c * (~x) ^ ~n), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, ~x) && NeQ((~b) ^ 2 - 4 * ~a * ~c, 0)))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')) ^ ~p, ~x) => Antiderivative(ExpandIntegrand((~d * ~x) ^ ~m * ~Pq * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, ~x) && ILtQ(~p + 1, 0)))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~(n') + ~(c') * (~x) ^ ~(n2')) ^ ~(p'), ~x) => Unintegrable((~d * ~x) ^ ~m * ~Pq * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (PolyQ(~Pq, ~x) || PolyQ(~Pq, (~x) ^ ~n)))
    @apply_utils Antiderivative((~u) ^ ~(m') * ~Pq * (~a + ~(b') * (~v) ^ ~n + ~(c') * (~v) ^ ~(n2')) ^ ~(p'), ~x) => ((~u) ^ ~m / (Coefficient(~v, ~x, 1) * (~v) ^ ~m)) * Subst(Antiderivative((~x) ^ ~m * SubstFor(~v, ~Pq, ~x) * (~a + ~b * (~x) ^ ~n + ~c * (~x) ^ (2 * ~n)) ^ ~p, ~x), ~x, ~v) <-- FreeQ([~a, ~b, ~c, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (LinearPairQ(~u, ~v, ~x) && PolyQ(~Pq, (~v) ^ ~n)))
end

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


text:
include("../benchmark/bench_morris_lecar.jl")

@time flc_solver = pre_solve!(make_morris_lecar_flc_solver())
@time sweep!(flc_solver)

@time flc_solver = pre_solve!(make_morris_lecar_flc_solver())
Profile.clear()
@profile sweep!(flc_solver)


text:
#=
# CG solver
=#
module CGSolver

export solve, nonlinear_solve

using LinearAlgebra, SparseArrays, CUDA

# See finch_import_symbols.jl for a list of all imported symbols.
import ..Finch: @import_finch_symbols
@import_finch_symbols()

include("fe_boundary.jl");
include("nonlinear.jl")
include("cg_matrixfree.jl");
include("level_benchmark.jl");

# Things to do before and after each step(or stage) ########
function default_pre_step() end
function default_post_step() end

pre_step_function = default_pre_step;
post_step_function = default_post_step;

function set_pre_step(fun) global pre_step_function = fun; end
function set_post_step(fun) global post_step_function = fun; end
############################################################

function init_solver()
    global pre_step_function = default_pre_step;
    global post_step_function = default_post_step;
end

function linear_solve(var, bilinear, linear, stepper=nothing; assemble_func=nothing)
    if config.linalg_matrixfree
        return solve_matrix_free_sym(var, bilinear, linear, stepper, assemble_func=assemble_func);
        #return solve_matrix_free_asym(var, bilinear, linear, stepper);
    end

    # If more than one variable
    if typeof(var) <: Array
        # multiple variables being solved for simultaneously
        dofs_per_node = 0;
        dofs_per_loop = 0;
        for vi=1:length(var)
            dofs_per_loop += length(var[vi].symvar);
            dofs_per_node += var[vi].total_components;
        end
    else
        # one variable
        dofs_per_loop = length(var.symvar);
        dofs_per_node = var.total_components;
    end
    N1 = size(grid_data.allnodes,2);
    Nn = dofs_per_node * N1;
    Np = refel.Np;
    nel = size(grid_data.loc2glb,2);

    # For partitioned meshes, keep these numbers handy
    (b_order, b_sizes) = get_partitioned_ordering(N1, dofs_per_node);

    # Allocate arrays that will be used by assemble
    rhsvec = zeros(Nn);
    lhsmatI = zeros(Int, nel*dofs_per_node*Np*dofs_per_node*Np);
    lhsmatJ = zeros(Int, nel*dofs_per_node*Np*dofs_per_node*Np);
    lhsmatV = zeros(nel*dofs_per_node*Np*dofs_per_node*Np);
    allocated_vecs = [rhsvec, lhsmatI, lhsmatJ, lhsmatV];

    if prob.time_dependent && !(stepper === nothing)
        # First assemble both lhs and rhs
        assemble_t = @elapsed begin
                        (A, b) = assemble(var, bilinear, linear, allocated_vecs, dofs_per_node, dofs_per_loop, 0, stepper.dt, assemble_loops=assemble_func);
                        (A, b) = gather_system(A, b, N1, dofs_per_node, b_order, b_sizes);
                    end

        log_entry("Initial assembly took "*string(assemble_t)*" seconds");

        log_entry("Beginning "*string(stepper.Nsteps)*" time steps.");
        t = 0;
        sol = get_var_vals(var);

        start_t = Base.Libc.time();
        assemble_t = 0;
        linsolve_t = 0;
        last2update = 0;
        last10update = 0;

        # allocate any temporary storage needed by steppers
        if stepper.type == LSRK4
            tmppi = zeros(size(b));
            tmpki = zeros(size(b));
        elseif stepper.stages > 1
            last_result = zeros(size(b));
            tmpresult = zeros(size(b));
            tmpki = zeros(length(b), stepper.stages);
        end
        print("Time stepping progress(%): 0");
        for i=1:stepper.Nsteps
            if stepper.stages > 1
                # LSRK4 is a special case, low storage
                if stepper.type == LSRK4
                    # Low storage RK4:
                    # p0 = u
                    #   ki = ai*k(i-1) + dt*f(p(i-1), t+ci*dt)
                    #   pi = p(i-1) + bi*ki
                    # u = p5

                    tmppi = get_var_vals(var, tmppi);
                    for rki=1:stepper.stages
                        rktime = t + stepper.c[rki]*stepper.dt;
                        # p(i-1) is currently in u
                        pre_step_function();

                        assemble_t += @elapsed begin
                                            b = assemble(var, nothing, linear, allocated_vecs, dofs_per_node, dofs_per_loop, rktime, stepper.dt; rhs_only = true, assemble_loops=assemble_func);
                                            b = gather_system(nothing, b, N1, dofs_per_node, b_order, b_sizes);
                                        end

                        linsolve_t += @elapsed(sol = distribute_solution( linear_system_solve(A,b) , N1, dofs_per_node, b_order, b_sizes));

                        # At this point sol holds the boundary values
                        # directly write them to the variable values and zero sol.
                        copy_bdry_vals_to_variables(var, sol, grid_data, dofs_per_node, zero_vals=true);

                        if rki == 1 # because a1 == 0
                            tmpki = stepper.dt .* sol;
                        else
                            tmpki = stepper.a[rki].*tmpki + stepper.dt.*sol;
                        end
                        tmppi = tmppi + stepper.b[rki].*tmpki

                        copy_bdry_vals_to_vector(var, tmppi, grid_data, dofs_per_node);
                        place_sol_in_vars(var, tmppi, stepper);

                        post_step_function();
                    end

                else
                    # Explicit multi-stage methods:
                    # x = x + dt*sum(bi*ki)
                    # ki = rhs(t+ci*dt, x+dt*sum(aij*kj)))   j < i

                    # will hold the final result
                    last_result = get_var_vals(var, last_result);
                    # will be placed in var.values for each stage
                    # Storage for each stage
                    # tmpki = zeros(length(b), stepper.stages);
                    for stage=1:stepper.stages
                        stime = t + stepper.c[stage]*stepper.dt;

                        # Update the values in vars to be used in this stage
                        if stage > 1
                            initialized_tmpresult = false;
                            for j=1:stage
                                if stepper.a[stage, j] > 0
                                    if !initialized_tmpresult
                                        initialized_tmpresult = true;
                                        for k=1:length(last_result)
                                            tmpresult[k] = last_result[k] + stepper.dt * stepper.a[stage, j] * tmpki[k,j];
                                        end
                                    else
                                        for k=1:length(last_result)
                                            tmpresult[k] += stepper.dt * stepper.a[stage, j] * tmpki[k,j];
                                        end
                                    end
                                end
                            end

                            if initialized_tmpresult
                                copy_bdry_vals_to_vector(var, tmpresult, grid_data, dofs_per_node);
                                place_sol_in_vars(var, tmpresult, stepper);
                            end
                            post_step_function(); # seems weird, but imagine this is happening after stage-1
                        end

                        pre_step_function();

                        assemble_t += @elapsed begin
                                            b = assemble(var, nothing, linear, allocated_vecs, dofs_per_node, dofs_per_loop, stime, stepper.dt; rhs_only = true, assemble_loops=assemble_func)
                                            b = gather_system(nothing, b, N1, dofs_per_node, b_order, b_sizes);
                                        end

                        linsolve_t += @elapsed(tmpki[:,stage] = distribute_solution( linear_system_solve(A,b) , N1, dofs_per_node, b_order, b_sizes));

                        # At this point tmpki[:,stage] holds the boundary values
                        # directly write them to the variable values and zero sol.
                        copy_bdry_vals_to_variables(var, tmpki[:,stage], grid_data, dofs_per_node, zero_vals=true);

                    end
                    for stage=1:stepper.stages
                        last_result += stepper.dt * stepper.b[stage] .* tmpki[:, stage];
                    end
                    copy_bdry_vals_to_vector(var, last_result, grid_data, dofs_per_node);
                    place_sol_in_vars(var, last_result, stepper);
                    post_step_function();
                end

            elseif stepper.type == EULER_EXPLICIT
                pre_step_function();
                assemble_t += @elapsed begin
                                            b = assemble(var, nothing, linear, allocated_vecs, dofs_per_node, dofs_per_loop, t, stepper.dt; rhs_only = true, assemble_loops=assemble_func)
                                            b = gather_system(nothing, b, N1, dofs_per_node, b_order, b_sizes);
                                        end

                linsolve_t += @elapsed(tmpvec = distribute_solution( linear_system_solve(A,b) , N1, dofs_per_node, b_order, b_sizes));

                # At this point tmpvec holds the boundary values
                # directly write them to the variable values and zero sol.
                copy_bdry_vals_to_variables(var, tmpvec, grid_data, dofs_per_node, zero_vals=true);

                sol = sol .+ stepper.dt .* tmpvec;

                copy_bdry_vals_to_vector(var, sol, grid_data, dofs_per_node);
                place_sol_in_vars(var, sol, stepper);
                post_step_function();

            elseif stepper.type == PECE
                # Predictor (explicit Euler)
                pre_step_function();
                assemble_t += @elapsed begin
                                            b = assemble(var, nothing, linear, allocated_vecs, dofs_per_node, dofs_per_loop, t, stepper.dt; rhs_only = true, assemble_loops=assemble_func)
                                            b = gather_system(nothing, b, N1, dofs_per_node, b_order, b_sizes);
                                        end

                linsolve_t += @elapsed(tmpvec = distribute_solution( linear_system_solve(A,b) , N1, dofs_per_node, b_order, b_sizes));

                # At this point tmpvec holds the boundary values
                # directly write them to the variable values and zero sol.
                copy_bdry_vals_to_variables(var, tmpvec, grid_data, dofs_per_node, zero_vals=true);

                tmpvec = sol .+ stepper.dt .* tmpvec;

                copy_bdry_vals_to_vector(var, tmpvec, grid_data, dofs_per_node);
                place_sol_in_vars(var, tmpvec, stepper);
                post_step_function();

                # Corrector (implicit Euler)
                pre_step_function();
                assemble_t += @elapsed begin
                                            b = assemble(var, nothing, linear, allocated_vecs, dofs_per_node, dofs_per_loop, t+stepper.dt, stepper.dt; rhs_only = true, assemble_loops=assemble_func)
                                            b = gather_system(nothing, b, N1, dofs_per_node, b_order, b_sizes);
                                        end

                linsolve_t += @elapsed(tmpvec = distribute_solution( linear_system_solve(A,b) , N1, dofs_per_node, b_order, b_sizes));

                # At this point tmpvec holds the boundary values
                # directly write them to the variable values and zero sol.
                copy_bdry_vals_to_variables(var, tmpvec, grid_data, dofs_per_node, zero_vals=true);

                sol = sol .+ stepper.dt .* tmpvec;

                copy_bdry_vals_to_vector(var, sol, grid_data, dofs_per_node);
                place_sol_in_vars(var, sol, stepper);
                post_step_function();

            else # implicit methods include the dt part
                pre_step_function();
                assemble_t += @elapsed begin
                                            b = assemble(var, nothing, linear, allocated_vecs, dofs_per_node, dofs_per_loop, t, stepper.dt; rhs_only = true, assemble_loops=assemble_func)
                                            b = gather_system(nothing, b, N1, dofs_per_node, b_order, b_sizes);
                                        end

                linsolve_t += @elapsed(sol = distribute_solution( linear_system_solve(A,b) , N1, dofs_per_node, b_order, b_sizes));

                place_sol_in_vars(var, sol, stepper);
                post_step_function();
            end

            t += stepper.dt;

            progressPercent = Int(floor(i*100.0/stepper.Nsteps));
            if progressPercent - last2update >= 2
                last2update = progressPercent;
                if progressPercent - last10update >= 10
                    print(progressPercent);
                    last10update = progressPercent;
                else
                    print(".");
                end
            end
        end
        println("");
        end_t = Base.Libc.time();

        log_entry("Stepping took "*string(end_t-start_t)*" seconds. ("*string(assemble_t)*" for assembly, "*string(linsolve_t)*" for linear solve)");
        #display(sol);
		# outfile = "linear_sol.txt"
		# open(outfile, "w") do f
  		# 	for ii in sol
    	# 		println(f, ii)
  		# 	end
		# end # the file f is automatically closed after this block finishes
        return sol;

    else
        assemble_t = @elapsed begin
                        (A, b) = assemble(var, bilinear, linear, allocated_vecs, dofs_per_node, dofs_per_loop, 0, 0, assemble_loops=assemble_func);
                        (A, b) = gather_system(A, b, N1, dofs_per_node, b_order, b_sizes);
                    end
        # uncomment to look at A
        # global Amat = A;

        sol_t = @elapsed(sol = distribute_solution( linear_system_solve(A,b) , N1, dofs_per_node, b_order, b_sizes));

        log_entry("Assembly took "*string(assemble_t)*" seconds");
        log_entry("Linear solve took "*string(sol_t)*" seconds");
        # display(A);
        # display(b);
        # display(sol);
        return sol;
    end
end

function nonlinear_solve(var, nlvar, bilinear, linear, stepper=nothing; assemble_loops=nothing)
    if config.num_pertitions > 1
        printerr("nonlinear solver is not ready for partitioned meshes. sorry.", fatal=true);
    end
    if prob.time_dependent && !(stepper === nothing)
        #TODO time dependent coefficients

        log_entry("Beginning "*string(stepper.Nsteps)*" time steps.");
        t = 0;
        start_t = Base.Libc.time();
        nl = nonlinear(100, 1e-9, 1e-9);
        init_nonlinear(nl, var, nlvar, bilinear, linear);
        for i=1:stepper.Nsteps
			println("solve for time step ", i);
			newton(nl, assemble, assemble_rhs_only, nlvar, t, stepper.dt);
            t += stepper.dt;
			nlvar[4].values = copy(nlvar[1].values);
			nlvar[5].values = copy(nlvar[2].values);
			if (i % 10 ==0)
				outfile = string("LD_u_",i);
				open(outfile, "w") do f
					for ii=1:length(nlvar[1].values)
						println(f, nlvar[1].values[ii])
					end
					close(f)
				end
				outfile = string("LD_v_",i);
				open(outfile, "w") do f
					for ii=1:length(nlvar[2].values)
						println(f, nlvar[2].values[ii])
					end
					close(f)
				end
			end
        end
        end_t = Base.Libc.time();

        log_entry("Solve took "*string(end_t-start_t)*" seconds");
        #display(sol);
        #return nlvar.values;
        return [];

    else
        start_t = Base.Libc.time();
        nl = nonlinear(100, 1e-12, 1e-12);
        init_nonlinear(nl, var, nlvar, bilinear, linear);
        newton(nl, assemble, assemble_rhs_only, nlvar);
        end_t = Base.Libc.time();

        log_entry("Solve took "*string(end_t-start_t)*" seconds");

        #return nlvar.values;
        return [];
    end
end

# assembles the A and b in Au=b
function assemble(var, bilinear, linear, allocated_vecs, dofs_per_node=1, dofs_per_loop=1, t=0, dt=0; assemble_loops=nothing, rhs_only = false)
    # If an assembly loop function was provided, use it
    if !(assemble_loops === nothing)
        return assemble_loops.func(var, bilinear, linear, allocated_vecs, dofs_per_node, dofs_per_loop, t, dt; rhs_only=rhs_only);
    end

    # Label things that were allocated externally
    b = allocated_vecs[1];
    if !rhs_only
        AI = allocated_vecs[2];
        AJ = allocated_vecs[3];
        AV = allocated_vecs[4];
    end
    # zero b
    b .= 0;

    Np = refel.Np;
    nel = size(grid_data.loc2glb,2);

    # Stiffness and mass are precomputed for uniform grid meshes
    precomputed_mass_stiffness = config.mesh_type == UNIFORM_GRID && config.geometry == SQUARE
    if precomputed_mass_stiffness
        wdetj = refel.wg .* geo_factors.detJ[1];
        J = geo_factors.J[1];

        if config.dimension == 1
            (RQ1, RD1) = @level_bench Level1 build_deriv_matrix(refel, J);
            TRQ1 = RQ1';
            stiffness = [(TRQ1 * diagm(wdetj) * RQ1)];
        elseif config.dimension == 2
            (RQ1, RQ2, RD1, RD2) = build_deriv_matrix(refel, J);
            (TRQ1, TRQ2) = (RQ1', RQ2');
            stiffness = [(TRQ1 * diagm(wdetj) * RQ1) , (TRQ2 * diagm(wdetj) * RQ2)];
        else
            (RQ1, RQ2, RQ3, RD1, RD2, RD3) = build_deriv_matrix(refel, J);
            (TRQ1, TRQ2, TRQ3) = (RQ1', RQ2', RQ3');
            stiffness = [(TRQ1 * diagm(wdetj) * RQ1) , (TRQ2 * diagm(wdetj) * RQ2) , (TRQ3 * diagm(wdetj) * RQ3)];
        end
        mass = (refel.Q)' * diagm(wdetj) * refel.Q;
    else
        stiffness = 0;
        mass = 0;
    end

    loop_time = Base.Libc.time();
    # Elemental loop follows elemental ordering
    for ei=1:nel
        eid = elemental_order[ei];
        loc2glb = grid_data.loc2glb[:,eid]; # global indices of this element's nodes

        if !rhs_only
            Astart = (eid-1)*Np*dofs_per_node*Np*dofs_per_node + 1; # The segment of AI, AJ, AV for this element
        end

        volargs = (var, eid, 0, grid_data, geo_factors, refel, t, dt, stiffness, mass);

        if dofs_per_node == 1
            linchunk = @level_bench Level1 linear.func(volargs);  # get the elemental linear part
            b[loc2glb] += linchunk;

            if !rhs_only
                bilinchunk = bilinear.func(volargs); # the elemental bilinear part
                #A[glb, glb] .+= bilinchunk;         # This will be very inefficient for sparse A
                for jj=1:Np
                    offset = Astart - 1 + (jj-1)*Np;
                    for ii=1:Np
                        AI[offset + ii] = loc2glb[ii];
                        AJ[offset + ii] = loc2glb[jj];
                        AV[offset + ii] = bilinchunk[ii, jj];
                    end
                end
            end

        else # more than one DOF per node
            linchunk = linear.func(volargs);
            insert_linear!(b, linchunk, loc2glb, 1:dofs_per_node, dofs_per_node);

            if !rhs_only
                bilinchunk = bilinear.func(volargs);
                @level_bench Level1  insert_bilinear!(AI, AJ, AV, Astart, bilinchunk, loc2glb, 1:dofs_per_node, dofs_per_node);
            end
        end
    end
    loop_time = Base.Libc.time() - loop_time;

    if !rhs_only
        # Build the sparse A. Uses default + to combine overlaps
        A = @level_bench Level1 sparse(AI, AJ, AV);
    end

    # Boundary conditions
    bc_time = Base.Libc.time();
    if rhs_only
        b =@level_bench Level1  apply_boundary_conditions_rhs_only(var, b, t);
    else
        (A, b) = @level_bench Level1 apply_boundary_conditions_lhs_rhs(var, A, b, t);
    end
    bc_time = Base.Libc.time() - bc_time;

    if rhs_only
        return b;

    else
        log_entry("Elemental loop time:     "*string(loop_time));
        log_entry("Boundary condition time: "*string(bc_time));
        return (A, b);
    end
end

function assemble_rhs_only(var, linear, allocated_vecs, dofs_per_node=1, dofs_per_loop=1, t=0, dt=0; assemble_loops=nothing)
    return assemble(var, nothing, linear, allocated_vecs, dofs_per_node, dofs_per_loop, t, dt; assemble_loops=assemble_loops, rhs_only=true)
end

# ######################################################
# # To be romoved. Fix nonlinear solve first.
# # assembles the A and b in Au=b
# function assemble_rhs_only(var, linear, t=0.0, dt=0.0; keep_geometric_factors = false, saved_geometric_factors = nothing)
#     Np = refel.Np;
#     nel = size(grid_data.loc2glb,2);
#     N1 = size(grid_data.allnodes,2);
#     multivar = typeof(var) <: Array;
#     maxvarindex = 0;
#     if multivar
#         # multiple variables being solved for simultaneously
#         dofs_per_node = 0;
#         var_to_dofs = [];
#         for vi=1:length(var)
#             tmp = dofs_per_node;
#             dofs_per_node += length(var[vi].symvar);
#             push!(var_to_dofs, (tmp+1):dofs_per_node);
#             maxvarindex = max(maxvarindex,var[vi].index);
#         end
#     else
#         # one variable
#         dofs_per_node = length(var.symvar);
#         maxvarindex = var.index;
#     end
#     Nn = dofs_per_node * N1;

#     b = zeros(Nn);

#     # Save geometric factors if desired
#     if keep_geometric_factors
#         saved_wdetj = [];
#         saved_J = [];
#     end

#     # Stiffness and mass are precomputed for uniform grid meshes
#     precomputed_mass_stiffness = config.mesh_type == UNIFORM_GRID && config.geometry == SQUARE
#     if precomputed_mass_stiffness
#         glb = grid_data.loc2glb[:,1];
#         xe = grid_data.allnodes[:,glb[:]];
#         (detJ, J) = geometric_factors(refel, xe);
#         wdetj = refel.wg .* detJ;
#         if keep_geometric_factors
#             saved_wdetj = wdetj;
#             saved_J = J;
#         end
#         if config.dimension == 1
#             (RQ1, RD1) = build_deriv_matrix(refel, J);
#             TRQ1 = RQ1';
#             stiffness = [(TRQ1 * diagm(wdetj) * RQ1)];
#         elseif config.dimension == 2
#             (RQ1, RQ2, RD1, RD2) = build_deriv_matrix(refel, J);
#             (TRQ1, TRQ2) = (RQ1', RQ2');
#             stiffness = [(TRQ1 * diagm(wdetj) * RQ1) , (TRQ2 * diagm(wdetj) * RQ2)];
#         else
#             (RQ1, RQ2, RQ3, RD1, RD2, RD3) = build_deriv_matrix(refel, J);
#             (TRQ1, TRQ2, TRQ3) = (RQ1', RQ2', RQ3');
#             stiffness = [(TRQ1 * diagm(wdetj) * RQ1) , (TRQ2 * diagm(wdetj) * RQ2) , (TRQ3 * diagm(wdetj) * RQ3)];
#         end
#         mass = (refel.Q)' * diagm(wdetj) * refel.Q;
#     else
#         stiffness = 0;
#         mass = 0;
#     end

#     # Elemental loop follows elemental ordering
#     for e=elemental_order;
#         glb = grid_data.loc2glb[:,e];       # global indices of this element's nodes for extracting values from var arrays
#         xe = grid_data.allnodes[:,glb[:]];  # coordinates of this element's nodes for evaluating coefficient functions

#         if !precomputed_mass_stiffness
#             if saved_geometric_factors === nothing
#                 (detJ, J) = geometric_factors(refel, xe);
#                 wdetj = refel.wg .* detJ;

#                 if keep_geometric_factors
#                     push!(saved_wdetj, wdetj);
#                     push!(saved_J, J);
#                 end
#             else
#                 wdetj = saved_geometric_factors[1][e];
#                 J = saved_geometric_factors[2][e];
#             end
#         end

#         rhsargs = (var, xe, glb, refel, wdetj, J, RHS, t, dt, stiffness, mass);

#         #linchunk = linear.func(args);  # get the elemental linear part
#         if dofs_per_node == 1
#             linchunk = linear.func(rhsargs);  # get the elemental linear part
#             b[glb] .+= linchunk;

#         elseif typeof(var) == Variable
#             # only one variable, but more than one dof
#             linchunk = linear.func(rhsargs);
#             insert_linear!(b, linchunk, glb, 1:dofs_per_node, dofs_per_node);

#         else
#             linchunk = linear.func(rhsargs);
#             insert_linear!(b, linchunk, glb, 1:dofs_per_node, dofs_per_node);

#         end
#     end

#     # Boundary conditions
#     bidcount = length(grid_data.bids); # the number of BIDs
#     if dofs_per_node > 1
#         if multivar
#             dofind = 0;
#             for vi=1:length(var)
#                 for compo=1:length(var[vi].symvar)
#                     dofind = dofind + 1;
#                     for bid=1:bidcount
#                         if prob.bc_type[var[vi].index, bid] == NO_BC
#                             # do nothing
#                         else
#                             b = dirichlet_bc_rhs_only(b, prob.bc_func[var[vi].index, bid][compo], grid_data.bdry[bid], t, dofind, dofs_per_node);
#                         end
#                     end
#                 end
#             end
#         else
#             for d=1:dofs_per_node
#                 #rows = ((d-1)*length(glb)+1):(d*length(glb));
#                 dofind = d;
#                 for bid=1:bidcount
#                     if prob.bc_type[var.index, bid] == NO_BC
#                         # do nothing
#                     else
#                         b = dirichlet_bc_rhs_only(b, prob.bc_func[var.index, bid][d], grid_data.bdry[bid], t, d, dofs_per_node);
#                     end
#                 end
#             end
#         end
#     else
#         for bid=1:bidcount
#             if prob.bc_type[var.index, bid] == NO_BC
#                 # do nothing
#             else
#                 b = dirichlet_bc_rhs_only(b, prob.bc_func[var.index, bid][1], grid_data.bdry[bid], t);
#             end
#         end
#     end

#     # Reference points
#     if size(prob.ref_point,1) >= maxvarindex
#         if multivar
#             posind = zeros(Int,0);
#             vals = zeros(0);
#             for vi=1:length(var)
#                 if prob.ref_point[var[vi].index,1]
#                     eii = prob.ref_point[var[vi].index, 2];
#                     tmp = (grid_data.glbvertex[eii[1], eii[2]] - 1)*dofs_per_node + var_to_dofs[vi][1];
#                     if length(prob.ref_point[var[vi].index, 3]) > 1
#                         tmp = tmp:(tmp+length(prob.ref_point[var[vi].index, 3])-1);
#                     end
#                     posind = [posind; tmp];
#                     vals = [vals; prob.ref_point[var[vi].index, 3]];
#                 end
#             end
#             if length(vals) > 0
#                 b[posind] = vals;
#             end

#         else
#             if prob.ref_point[var.index,1]
#                 eii = prob.ref_point[var.index, 2];
#                 posind = (grid_data.glbvertex[eii[1], eii[2]] - 1)*dofs_per_node + 1;
#                 if length(prob.ref_point[var.index, 3]) > 1
#                     posind = posind:(posind+length(prob.ref_point[var[vi].index, 3])-1);
#                 end
#                 b[posind] = prob.ref_point[var.index, 3];
#             end
#         end
#     end

#     return b;
# end
# ##########################################################

# Inset the single dof into the greater construct
function insert_linear!(b, bel, glb, dof, Ndofs)
    # group nodal dofs
    for d=1:length(dof)
        ind = glb.*Ndofs .- (Ndofs-dof[d]);
        ind2 = ((d-1)*length(glb)+1):(d*length(glb));

        b[ind] = b[ind] + bel[ind2];
    end
end

function insert_bilinear!(AI, AJ, AV, Astart, ael, glb, dof, Ndofs)
    Np = length(glb);
    # group nodal dofs
    for dj=1:length(dof)
        indj = glb.*Ndofs .- (Ndofs-dof[dj]);
        indj2 = ((dj-1)*Np+1):(dj*Np);
        for di=1:length(dof)
            indi = glb.*Ndofs .- (Ndofs-dof[di]);
            indi2 = ((di-1)*Np+1):(di*Np);

            #a[indi, indj] = a[indi, indj] + ael[indi2, indj2];

            for jj=1:Np
                offset = Astart + (jj-1 + Np*(dj-1))*Np*Ndofs + Np*(di-1) - 1;
                for ii=1:Np
                    AI[offset + ii] = indi[ii];
                    AJ[offset + ii] = indj[jj];
                    AV[offset + ii] = ael[indi2[ii], indj2[jj]];
                end
            end
        end
    end

end

function place_sol_in_vars(var, sol, stepper)
    # place the values in the variable value arrays
    if typeof(var) <: Array
        tmp = 0;
        totalcomponents = 0;
        for vi=1:length(var)
            totalcomponents = totalcomponents + var[vi].total_components;
        end
        for vi=1:length(var)
            components = var[vi].total_components;
            for compi=1:components
                if stepper.type == EULER_EXPLICIT
                    var[vi].values[compi,:] = sol[(compi+tmp):totalcomponents:end];
                else
                    var[vi].values[compi,:] = sol[(compi+tmp):totalcomponents:end];
                end
                tmp = tmp + 1;
            end
        end
    else
        components = var.total_components;
        for compi=1:components
            if stepper.type == EULER_EXPLICIT
                var.values[compi,:] = sol[compi:components:end];
            else
                var.values[compi,:] = sol[compi:components:end];
            end
        end
    end
end

function get_var_vals(var, vect=nothing)
    # place the variable values in a vector
    if typeof(var) <: Array
        tmp = 0;
        totalcomponents = 0;
        for vi=1:length(var)
            totalcomponents = totalcomponents + var[vi].total_components;
        end
        if vect === nothing
            vect = zeros(totalcomponents * length(var[1].values[1,:]));
        end
        for vi=1:length(var)
            components = var[vi].total_components;
            for compi=1:components
                vect[(compi+tmp):totalcomponents:end] = var[vi].values[compi,:];
                tmp = tmp + 1;
            end
        end
    else
        components = var.total_components;
        if vect === nothing
            vect = zeros(components * length(var.values[1,:]));
        end
        for compi=1:components
            vect[compi:components:end] = var.values[compi,:];
        end
    end

    return vect;
end

######################################################################################
# For partitioned meshes
function get_partitioned_ordering(nnodes, dofs_per_node)
    if config.num_procs > 1
        # The global ordering of this partition's b vector entries
        b_order_loc = zeros(Int, nnodes*dofs_per_node);
        for ni=1:nnodes
            for di=1:dofs_per_node
                b_order_loc[(ni-1)*dofs_per_node + di] = (grid_data.partition2global[ni]-1)*dofs_per_node + di;
            end
        end

        # only proc 0 needs this info for now
        if config.proc_rank == 0
            # The b_sizes
            send_buf = [nnodes*dofs_per_node];
            recv_buf = zeros(Int, config.num_procs);
            MPI.Gather!(send_buf, recv_buf, 0, MPI.COMM_WORLD);

            # b_order
            chunk_sizes = recv_buf;
            displacements = zeros(Int, config.num_procs); # for the irregular gatherv
            total_length = chunk_sizes[1];
            for proc_i=2:config.num_procs
                displacements[proc_i] = displacements[proc_i-1] + chunk_sizes[proc_i-1];
                total_length += chunk_sizes[proc_i];
            end
            full_b_order = zeros(Int, total_length);
            b_order_buf = MPI.VBuffer(full_b_order, chunk_sizes, displacements, MPI.Datatype(Int));
            MPI.Gatherv!(b_order_loc, b_order_buf, 0, MPI.COMM_WORLD);

            return (full_b_order, chunk_sizes);

        else
            send_buf = [nnodes*dofs_per_node];
            MPI.Gather!(send_buf, nothing, 0, MPI.COMM_WORLD);
            MPI.Gatherv!(b_order_loc, nothing, 0, MPI.COMM_WORLD);

            return ([1], [1]);
        end

    else
        return (nothing, nothing);
    end
end

# For multiple processes, gather the system, distribute the solution
# Note: rescatter_b only applies when rhs_only
function gather_system(A, b, nnodes, dofs_per_node, b_order, b_sizes; rescatter_b=false)
    rhs_only = (A===nothing);
    if config.num_procs > 1
        if !rhs_only
            # For now just gather all of A in proc 0 to assemble.
            # The row and column indices have to be changed according to partition2global.
            # Also, b must be reordered on proc 0, so send the needed indices as well.
            (AI, AJ, AV) = findnz(A);
            for i=1:length(AI)
                dof = mod(AI[i]-1,dofs_per_node)+1;
                node = Int(floor((AI[i]-dof) / dofs_per_node) + 1);
                AI[i] = (grid_data.partition2global[node] - 1) * dofs_per_node + dof;

                dof = mod(AJ[i]-1,dofs_per_node)+1;
                node = Int(floor((AJ[i]-dof) / dofs_per_node) + 1);
                AJ[i] = (grid_data.partition2global[node] - 1) * dofs_per_node + dof;
            end

            if config.proc_rank == 0
                # First figure out how long each proc's arrays are
                send_buf = [length(AI)];
                recv_buf = zeros(Int, config.num_procs);
                MPI.Gather!(send_buf, recv_buf, 0, MPI.COMM_WORLD);

                # Use gatherv to accumulate A
                chunk_sizes = recv_buf;
                displacements = zeros(Int, config.num_procs); # for the irregular gatherv
                total_length = chunk_sizes[1];
                for proc_i=2:config.num_procs
                    displacements[proc_i] = displacements[proc_i-1] + chunk_sizes[proc_i-1];
                    total_length += chunk_sizes[proc_i];
                end
                full_AI = zeros(Int, total_length);
                full_AJ = zeros(Int, total_length);
                full_AV = zeros(Float64, total_length);
                AI_buf = MPI.VBuffer(full_AI, chunk_sizes, displacements, MPI.Datatype(Int));
                AJ_buf = MPI.VBuffer(full_AJ, chunk_sizes, displacements, MPI.Datatype(Int));
                AV_buf = MPI.VBuffer(full_AV, chunk_sizes, displacements, MPI.Datatype(Float64));

                MPI.Gatherv!(AI, AI_buf, 0, MPI.COMM_WORLD);
                MPI.Gatherv!(AJ, AJ_buf, 0, MPI.COMM_WORLD);
                MPI.Gatherv!(AV, AV_buf, 0, MPI.COMM_WORLD);

                # # Modify AI and AJ to global indices using b_order
                # b_start = 0;
                # for proc_i=1:config.num_procs
                #     for ai=(displacements[proc_i]+1):(displacements[proc_i] + chunk_sizes[proc_i])
                #         full_AI[ai] = b_order[b_start + full_AI[ai]];
                #         full_AJ[ai] = b_order[b_start + full_AJ[ai]];
                #     end
                #     b_start += b_sizes[proc_i];
                # end

                # Assemble A
                full_A = sparse(full_AI, full_AJ, full_AV);

                # Next gather b
                chunk_sizes = b_sizes;
                displacements = zeros(Int, config.num_procs); # for the irregular gatherv
                total_length = chunk_sizes[1];
                for proc_i=2:config.num_procs
                    displacements[proc_i] = displacements[proc_i-1] + chunk_sizes[proc_i-1];
                    total_length += chunk_sizes[proc_i];
                end
                full_b = zeros(total_length);
                b_buf = MPI.VBuffer(full_b, chunk_sizes, displacements, MPI.Datatype(Float64));
                MPI.Gatherv!(b, b_buf, 0, MPI.COMM_WORLD);

                # Overlapping values will be added.
                new_b = zeros(grid_data.nnodes_global * dofs_per_node);
                for i=1:total_length
                    new_b[b_order[i]] += full_b[i];
                end

                return (full_A, new_b);

            else # other procs just send their data
                send_buf = [length(AI)];
                MPI.Gather!(send_buf, nothing, 0, MPI.COMM_WORLD);
                MPI.Gatherv!(AI, nothing, 0, MPI.COMM_WORLD);
                MPI.Gatherv!(AJ, nothing, 0, MPI.COMM_WORLD);
                MPI.Gatherv!(AV, nothing, 0, MPI.COMM_WORLD);
                MPI.Gatherv!(b, nothing, 0, MPI.COMM_WORLD);

                return (ones(1,1),[1]);
            end

        else # RHS only\
            if config.proc_rank == 0
                # gather b
                chunk_sizes = b_sizes;
                displacements = zeros(Int, config.num_procs); # for the irregular gatherv
                total_length = chunk_sizes[1];
                for proc_i=2:config.num_procs
                    displacements[proc_i] = displacements[proc_i-1] + chunk_sizes[proc_i-1];
                    total_length += chunk_sizes[proc_i];
                end
                full_b = zeros(total_length);
                b_buf = MPI.VBuffer(full_b, chunk_sizes, displacements, MPI.Datatype(Float64));
                MPI.Gatherv!(b, b_buf, 0, MPI.COMM_WORLD);

                # Overlapping values will be added.
                new_b = zeros(grid_data.nnodes_global * dofs_per_node);
                for i=1:total_length
                    new_b[b_order[i]] += full_b[i];
                end

                if rescatter_b
                    return distribute_solution(new_b, nnodes, dofs_per_node, b_order, b_sizes)
                else
                    return new_b;
                end

            else # other procs just send their data
                MPI.Gatherv!(b, nothing, 0, MPI.COMM_WORLD);

                if rescatter_b
                    return distribute_solution(nothing, nnodes, dofs_per_node, b_order, b_sizes)
                else
                    return [1];
                end
            end
        end

    else # one process
        if rhs_only
            return b;
        else
            return (A, b);
        end
    end
end

function distribute_solution(sol, nnodes, dofs_per_node, b_order, b_sizes)
    if config.num_procs > 1
        my_sol = zeros(nnodes*dofs_per_node);
        if config.proc_rank == 0 # 0 has the full sol
            # Need to reorder b and put in a larger array according to b_order
            total_length = sum(b_sizes);
            full_sol = zeros(total_length);
            for i=1:total_length
                full_sol[i] = sol[b_order[i]];
            end

            # scatter b
            chunk_sizes = b_sizes;
            displacements = zeros(Int, config.num_procs); # for the irregular scatterv
            for proc_i=2:config.num_procs
                displacements[proc_i] = displacements[proc_i-1] + chunk_sizes[proc_i-1];
            end
            sol_buf = MPI.VBuffer(full_sol, chunk_sizes, displacements, MPI.Datatype(Float64));

            # Scatter it amongst the little ones
            MPI.Scatterv!(sol_buf, my_sol, 0, MPI.COMM_WORLD);

        else # Others don't have sol
            MPI.Scatterv!(nothing, my_sol, 0, MPI.COMM_WORLD);
        end
        return my_sol;

    else
        return sol;
    end
end

# Simply does a reduction.
# Works for scalar values or vectors, but vectors are not themselves reduced.
# 1, 2, 3 -> 6
# [1,10,100], [2,20,200], [3,30,300] -> [6,60,600]
function combine_values(val; combine_op = +)
    if config.num_procs > 1
        rval = MPI.Allreduce(val, combine_op, MPI.COMM_WORLD);
        return rval;

    else
        return val;
    end
end

# This reduces a global vector to one value
# [1,10,100], [2,20,200], [3,30,300] -> 666
function reduce_vector(vec::Array)
    if config.num_procs > 1
        rval = MPI.Allreduce(sum(vec), +, MPI.COMM_WORLD);
        return rval;

    else
        return sum(vec);
    end
end

################################################################
## Options for solving the assembled linear system
################################################################
function linear_system_solve(A,b)
    if config.num_procs > 1 && config.proc_rank > 0
        # Other procs don't have the full A, just return b
        return b;
    end

    if config.linalg_backend == DEFAULT_SOLVER
        return A\b;
    elseif config.linalg_backend == PETSC_SOLVER
        # For now try this simple setup.
        # There will be many options that need to be available to the user. TODO
        # The matrix should really be constructed from the begninning as a PETSc one. TODO

        petsclib = PETSc.petsclibs[1]           #
        inttype = PETSc.inttype(petsclib);      # These should maybe be kept in config
        scalartype = PETSc.scalartype(petsclib);#

        (I, J, V) = findnz(A);
        (n,n) = size(A);
        nnz = zeros(inttype, n); # number of non-zeros per row
        for i=1:length(I)
            nnz[I[i]] += 1;
        end

        petscA = PETSc.MatSeqAIJ{scalartype}(n,n,nnz);
        for i=1:length(I)
            petscA[I[i],J[i]] = V[i];
        end
        PETSc.assemble(petscA);

        ksp = PETSc.KSP(petscA; ksp_rtol=1e-8, pc_type="jacobi", ksp_monitor=false); # Options should be available to user

        return ksp\b;
    elseif config.linalg_backend == CUDA_SOLVER
        elty = typeof(b[1])
        x = zeros(elty ,length(b))
        tol = convert(real(elty),1e-6)
        x = CUDA.CUSOLVER.csrlsvlu!(A, b, x, tol ,one(Cint),'O')
        return x
    end
end

end #module

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
module PolaronMakie

export viewfunction, frohlich_α, feynman_variation, osaka_free_energy, ℜD, ℑD, D, ℜS, ℑS, S, ℜχ, ℑχ, χ

using QuadGK
using AbstractPlotting.MakieLayout
using AbstractPlotting
# using GLMakie
using Optim

include("../src/function.jl")
include("../src/Functions/D(u)_function.jl")
include("../src/Functions/S(u)_function.jl")
include("../src/Functions/X(u)_function.jl")
include("../src/Polaron/frohlich.jl")
include("../src/Polaron/feynmantheory.jl")
include("../src/Polaron/osaka.jl")
include("../src/Polaron/hellwarththeory.jl")

# Physical constants
const ħ = 1.05457162825e-34; # kg m^2 s^{-1
const eV = 1.05457162825e-34; # kg m^2 s^{-2}
const m_e = 9.10938188e-31; # kg
const k_B = 1.3806504e-23; # kg m^2 K^{-1} s^2
const ϵ_0 = 8.854e-12 # C^2 N^{-1} m^{-2}

struct Polaron
    α::Float64      # Frohlich alpha
    v::Float64      # Variational parameter
    w::Float64      # Variational parameter
    F::Float64      # Free energy
    E::Float64      # Ground-state energy
    m_b::Float64    # Effective mass
    μ::Float64      # Mobility
end

function Polaron(α, v, w, F, E, m_b, μ)
    new(α, v, w, F, E, m_b, μ)
end

function Base.show(io::IO, quantities::Polaron)
    print(io, "$α, $v, $w, $F, $E, $m_b, $μ")
end

end # PolaronMakie


text:
# module Ecto

import Base: in
import .MirrorTypes: schema_to_mirrorstruct, mirrorstruct_to_mirror
import .MirrorTypes: MirrorModel
import .Schema: Assoc

function in(modul::Module)
    struct = schema_to_mirrorstruct(modul)
    mirrorstruct_to_mirror(struct)
end

function from{Model<:MirrorModel}(model::Model; kw...)::Query.t
    vector = Assoc(kw)
    prefix = haskey(vector, :prefix) ? vector[:prefix] : nothing
    sources = [v for (k,v) in vector if :source==k]
    from = model
    joins = [v for (k,v) in vector if :join==k]
    wheres = [v for (k,v) in vector if :where==k]
    select = haskey(vector, :select) ? vector[:select] : nothing
    order_bys = [v for (k,v) in vector if :order_by==k]
    limit = haskey(vector, :limit) ? vector[:limit] : nothing
    offset = haskey(vector, :offset) ? vector[:offset] : nothing
    group_bys = [v for (k,v) in vector if :group_by==k]
    updates = [v for (k,v) in vector if :update==k]
    havings = [v for (k,v) in vector if :having==k]
    preloads = [v for (k,v) in vector if :preload==k]
    assocs = [v for (k,v) in vector if :assoc==k]
    distinct = haskey(vector, :distinct) ? vector[:distinct] : nothing
    Query.t(prefix, sources, from, joins, wheres, select, order_bys, limit, offset, group_bys, updates, havings, preloads, assocs, distinct, lock)
end

function and
end

function set
end

## Errors
type QueryError
    message
end

module Query

import Base: |>, isempty

abstract Queryable

type Applicable
end

type t <: Queryable
    prefix
    sources::Vector
    from
    joins::Vector
    wheres::Vector
    select
    order_bys::Vector
    limit
    offset
    group_bys::Vector
    updates::Vector
    havings::Vector
    preloads::Vector
    assocs::Vector
    distinct
    lock
    t() = new(nothing,[],nothing,[],[],nothing,[],nothing,nothing,[],[],[],[],[],nothing,nothing)
    t(prefix, sources, from, joins, wheres, select, order_bys, limit, offset, group_bys, updates, havings, preloads, assocs, distinct, lock) = new(prefix, sources, from, joins, wheres, select, order_bys, limit, offset, group_bys, updates, havings, preloads, assocs, distinct, lock)
end

## Errors
type CastError
    message
end

function |>(queryable::Query.t, applicable::Query.Applicable)::Query.t
    queryable
end

function isempty(queryable::Query.t, key::Symbol)::Bool
    value = getfield(queryable, key)
    isa(value, Vector) ? isempty(value) : isa(value, Void)
end

end # module Ecto.Query


text:
@testset "bsuite" begin
    @testset "bsuite_cartpole_shuffled" begin
        ds = rl_unplugged_bsuite_dataset(
            "cartpole",
            [1, 2],
            "full";
            is_shuffle = true,
            stochasticity = 0.0,
            shuffle_buffer_size = 10_000,
            tf_reader_bufsize = 10_000,
            tf_reader_sz = 10_000,
            batch_size = 256,
            n_preallocations = Threads.nthreads() * 12
        )

        @test typeof(ds)<:RingBuffer

        s_size = 6

        data_1 = take!(ds)

        @test size(data_1.state) == (s_size, batch_size)
        @test size(data_1.next_state) == (s_size, batch_size)
        @test size(data_1.action) == (batch_size,)
        @test size(data_1.reward) == (batch_size,)
        @test size(data_1.terminal) == (batch_size,)

        @test typeof(data_1.state) == Array{Float32, 2}
        @test typeof(data_1.next_state) == Array{Float32, 2}
        @test typeof(data_1.action) == Array{Int, 1}
        @test typeof(data_1.reward) == Array{Float32, 1}
        @test typeof(data_1.terminal) == Array{Bool, 1}

    end

    @testset "bsuite_cartpole" begin
        ds = rl_unplugged_bsuite_dataset(
            "cartpole",
            [1, 2],
            "full";
            is_shuffle = false,
            stochasticity = 0.0,
            shuffle_buffer_size = 10_000,
            tf_reader_bufsize = 10_000,
            tf_reader_sz = 10_000,
            batch_size = 256,
            n_preallocations = Threads.nthreads() * 12
        )

        @test typeof(ds)<:RingBuffer

        s_size = 6

        data_1 = take!(ds)

        @test size(data_1.state) == (s_size, batch_size)
        @test size(data_1.next_state) == (s_size, batch_size)
        @test size(data_1.action) == (batch_size,)
        @test size(data_1.reward) == (batch_size,)
        @test size(data_1.terminal) == (batch_size,)

        @test typeof(data_1.state) == Array{Float32, 2}
        @test typeof(data_1.next_state) == Array{Float32, 2}
        @test typeof(data_1.action) == Array{Int, 1}
        @test typeof(data_1.reward) == Array{Float32, 1}
        @test typeof(data_1.terminal) == Array{Bool, 1}

    end

    @testset "bsuite_catch_shuffle" begin
        ds = rl_unplugged_bsuite_dataset(
            "catch",
            [1, 2],
            "full";
            is_shuffle = true,
            stochasticity = 0.0,
            shuffle_buffer_size = 10_000,
            tf_reader_bufsize = 10_000,
            tf_reader_sz = 10_000,
            batch_size = 256,
            n_preallocations = Threads.nthreads() * 12
        )

        @test typeof(ds)<:RingBuffer

        s_size = (10, 5)

        data_1 = take!(ds)

        @test size(data_1.state) == (s_size[1], s_size[2], batch_size)
        @test size(data_1.next_state) == (s_size[1], s_size[2], batch_size)
        @test size(data_1.action) == (batch_size,)
        @test size(data_1.reward) == (batch_size,)
        @test size(data_1.terminal) == (batch_size,)

        @test typeof(data_1.state) == Array{Float32, 3}
        @test typeof(data_1.next_state) == Array{Float32, 3}
        @test typeof(data_1.action) == Array{Int, 1}
        @test typeof(data_1.reward) == Array{Float32, 1}
        @test typeof(data_1.terminal) == Array{Bool, 1}

    end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Hexahedron(p1, p2, ..., p8)

A hexahedron with points `p1`, `p2`, ..., `p8`.
"""
struct Hexahedron{Dim,T,V<:AbstractVector{Point{Dim,T}}} <: Polyhedron{Dim,T}
  vertices::V
end