text:
"""
Compute the power of a test without checking arguments. Useful inside of loops
that need to evaluate the power of many proposed designs.
"""
function analytic_power(
    ::Type{T},
    ns::Union{Real, Tuple},
    δ::Real,
    σs::Union{Real, Tuple},
    α::Real,
    one_sided::Bool,
) where {T}
    null, alt = hypotheses(T, ns, δ, σs)
    if one_sided
        c_r = cquantile(null, α)
        return ccdf(alt, c_r)
    else
        c_l = quantile(null, α / 2)
        c_r = cquantile(null, α / 2)
        return cdf(alt, c_l) + ccdf(alt, c_r)
    end
end

"""
Compute the power of a test. Users must specify:

* `ns`: Per-group sample size(s)
* `δ`: Deviation from the null
* `σs`: Per-group standard deviations
* `α`: Significance level (defaults to 0.05)
* `one_sided`: Is the test one-sided? (defaults to false)
"""
function power(
    ::Type{T},
    ns::Union{Real, Tuple},
    δ::Real,
    σs::Union{Real, Tuple},
    α::Real = 0.05,
    one_sided::Bool = false,
) where {T}
    check_args(ns, δ, σs, 0.8, α)
    return analytic_power(T, ns, δ, σs, α, one_sided)
end

"""
Compute the smallest effect size measurable using a test. Users must specify:

* `ns`: Per-group sample size(s)
* `p`: Desired power
* `α`: Significance level (defaults to 0.05)
* `one_sided`: Is the test one-sided? (defaults to false)
"""
function effect_size(
    ::Type{T},
    ns::Union{Real, Tuple},
    p::Real,
    α::Real = 0.05,
    one_sided::Bool = false,
) where {T}
    check_args(ns, 1.0, 1.0, p, α)
    gap(δ::Real) = analytic_power(T, ns, δ, 1.0, α, one_sided) - p
    return fzero(gap, 1e-8, 100.0)
end

"""
Compute the smallest sample size that allows using a test to measure a
specified effect. Users must specify:

* `δ`: Deviation from the null
* `σs`: Per-group standard deviations
* `p`: Desired power
* `α`: Significance level (defaults to 0.05)
* `one_sided`: Is the test one-sided? (defaults to false)
"""
function sample_size(
    ::Type{T},
    δ::Real,
    σs::Union{Real, Tuple},
    p::Real,
    α::Real = 0.05,
    one_sided::Bool = false,
) where {T}
    check_args(2, δ, σs, p, α)
    gap(n::Real) = analytic_power(T, n, δ, σs, α, one_sided) - p
    return fzero(gap, 2.0, 1.0e12)
end

"""
Compute the type 1 error rate of a test. Users must specify:

* `ns`: Per-group sample size(s)
* `δ`: Deviation from the null
* `σs`: Per-group standard deviations
* `α`: Significance level (defaults to 0.05)
* `one_sided`: Is the test one-sided? (defaults to false)
"""
function type_1(
    ::Type{T},
    ns::Union{Real, Tuple},
    δ::Real,
    σs::Union{Real, Tuple},
    α::Real = 0.05,
    one_sided::Bool = false,
) where {T}
    check_args(ns, δ, σs, 0.8, α)
    return α
end

"""
Compute the type 2 error rate of a test. Users must specify:

* `ns`: Per-group sample size(s)
* `δ`: Deviation from the null
* `σs`: Per-group standard deviations
* `α`: Significance level (defaults to 0.05)
* `one_sided`: Is the test one-sided? (defaults to false)
"""
function type_2(
    ::Type{T},
    ns::Union{Real, Tuple},
    δ::Real,
    σs::Union{Real, Tuple},
    α::Real = 0.05,
    one_sided::Bool = false,
) where {T}
    check_args(ns, δ, σs, 0.8, α)
    return 1 - power(T, ns, δ, σs, α, one_sided)
end

"""
Compute the type S error rate of a test. Users must specify:

* `ns`: Per-group sample size(s)
* `δ`: Deviation from the null
* `σs`: Per-group standard deviations
* `α`: Significance level (defaults to 0.05)
* `one_sided`: Is the test one-sided? (defaults to false)
"""
function type_s(
    ::Type{T},
    ns::Union{Real, Tuple},
    δ::Real,
    σs::Union{Real, Tuple},
    α::Real = 0.05,
) where {T}
    check_args(ns, δ, σs, 0.8, α)
    null, alt = hypotheses(T, ns, δ, σs)
    c_l, c_r = quantile(null, α / 2), cquantile(null, α / 2)
    p_l, p_r = cdf(alt, c_l), ccdf(alt, c_r)
    return p_l / (p_l + p_r)
end

"""
Compute the type M error rate of a test. Users must specify:

* `ns`: Per-group sample size(s)
* `δ`: Deviation from the null
* `σs`: Per-group standard deviations
* `α`: Significance level (defaults to 0.05)
* `one_sided`: Is the test one-sided? (defaults to false)
"""
function type_m(
    ::Type{T},
    ns::Union{Real, Tuple},
    δ::Real,
    σs::Union{Real, Tuple},
    α::Real = 0.05,
    one_sided::Bool = false,
) where {T}
    check_args(ns, δ, σs, 0.8, α)
    null, alt = hypotheses(T, ns, δ, σs)
    if one_sided
        c_r = cquantile(null, α)
        t_l, t_r = thresholds(alt)
        p_significant = ccdf(alt, c_r)
        p_overestimate = ccdf(alt, max(c_r, t_r))
        return p_overestimate / p_significant
    else
        c_l, c_r = quantile(null, α / 2), cquantile(null, α / 2)
        t_l, t_r = thresholds(alt)
        p_significant = cdf(alt, c_l) + ccdf(alt, c_r)
        p_overestimate = cdf(alt, min(c_l, t_l)) + ccdf(alt, max(c_r, t_r))
        return p_overestimate / p_significant
    end
end

"""
Compute the average exaggeration factor of a test. Users must specify:

* `ns`: Per-group sample size(s)
* `δ`: Deviation from the null
* `σs`: Per-group standard deviations
* `α`: Significance level (defaults to 0.05)
* `one_sided`: Is the test one-sided? (defaults to false)
"""
function exaggeration_factor(
    ::Type{T},
    ns::Union{Real, Tuple},
    δ::Real,
    σs::Union{Real, Tuple},
    α::Real = 0.05,
    one_sided::Bool = false,
    bounds::Real = 1e-8,
    tol::Real = 1e-8,
) where {T}
    check_args(ns, δ, σs, 0.8, α)
    null, alt = hypotheses(T, ns, δ, σs)
    if one_sided
        c_r = cquantile(null, α)
        t_l, t_r = thresholds(alt)
        p_significant = ccdf(alt, c_r)
        f0(t) = (t / t_r) * (pdf(alt, t) / p_significant)
        i_r = cquantile(alt, bounds)
        e_r = gauss(f0, c_r, i_r, rtol = tol)[1]
        return e_r
    else
        c_l, c_r = quantile(null, α / 2), cquantile(null, α / 2)
        t_l, t_r = thresholds(alt)
        p_significant = cdf(alt, c_l) + ccdf(alt, c_r)
        f1(t) = (t / t_l) * (pdf(alt, t) / p_significant)
        f2(t) = (t / t_r) * (pdf(alt, t) / p_significant)
        i_l, i_r = quantile(alt, bounds), cquantile(alt, bounds)
        e_l = gauss(f1, i_l, c_l, rtol = tol)[1]
        e_r = gauss(f2, c_r, i_r, rtol = tol)[1]
        return e_l + e_r
    end
end


text:
"""
    bab_sqrt(S::Real; tolerance = 1e-6, guess = nothing)

The Babylonian Method of calculating a square root is a simple iterative method to determine square roots. A full description of the algorithm can be found on [Wikipedia](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method)

# Arguments:
- `S`: The number to calculate the square root for.

# Positional Arguments
- `tolerance`: How close the square of the square root needs to be from the input value. `abs(S - xn^2) < tolerance`
- `guess`: The initial value to use for `xn`

# Examples/Tests
```julia
julia> bab_sqrt(100)
10.000000000107445

julia> bab_sqrt(100, guess = 15)
10.000000000131072

julia> bab_sqrt(π, guess = 1)
1.7724538555800293

julia> bab_sqrt(π, guess = 1, tolerance = 2)
2.0707963267948966
```

# Algorithm:

```julia
while tolerance <= abs(xn^2 - S)
    xn = (1 / 2) * (xn + S / xn)
end
```

# References:
[Methods of computing square roots](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method)

```

# Contributed by:- [Anson Biggs](https://ansonbiggs.com)
"""
function bab_sqrt(S::Real; tolerance = 1e-6, guess = nothing)

    # Can only calculate for positive numbers
    if S < 0
        throw(DomainError("Input must be greater than zero"))
    end

    xn = guess
    if xn === nothing
        xn = S / 2
    elseif xn < 0
        xn = abs(xn)
    end

    while tolerance <= abs(xn^2 - S)
        xn = (1 / 2) * (xn + S / xn)
    end
    return xn

end

ext: jl

max_stars_repo_name: standardgalactic/Julia-algorithms


text:
using Singularity
using Test

@testset "Singularity.jl" begin
    # Write your own tests here.
end


text:
#=
read_rdb_hdr.jl
based on
https://gitlab.com/fMRI/toppe/blob/master/+toppe/+utils/read_rdb_hdr.m
=#

export read_rdb_hdr

include("rdb-26_002.jl")

"""
`s = read_rdb_hdr(fid::IOStream)`

Read GE raw (RDB) header for MRI scans

Returns NamedTuple `s` with header values accessible by `s.key`

Matlab version is:
Copyright (c) 2012 by General Electric Company. All rights reserved.

2019-05-22 Julia version by Jeff Fessler
"""
function read_rdb_hdr(fid::IOStream)
	seek(fid, 0)
	rdbm_rev = read(fid, Float32) # raw header (RDBM) revision number

	if rdbm_rev == Float32(26.002)
		s = read_rdb_hdr_26_002(fid)
	else
		throw("unknown RDBM rev $rdbm_rev")
	end

	if (s.rdbm_rev != rdbm_rev)
		throw("rev mismatch: $(s.rdbm_rev) != $rdbm_rev")
	end

	return s
end

"""
`s = read_rdb_hdr(file::String)`
read from `file`
"""
function read_rdb_hdr(file::String)
	return open(read_rdb_hdr, file)
end

"""
`read_rdb_hdr(:test)`
self test
"""
function read_rdb_hdr(test::Symbol)
	file = "/n/ir71/d3/fessler/fmri-data-michelle-L+S/P97792.7"
	if isfile(file)
		return read_rdb_hdr(file).dab[2] == 31
	end
	true
end


text:
module Data

using MLUtils
export DataLoader

end#module


text:
# Header parsing logic, as in https://github.com/cjdoris/ARFFFiles.jl

const _SPACE = UInt8(' ')
const _AT = UInt8('@')
const _COMMENT = UInt8('%')
const _DQUO = UInt8('"')
const _SQUO = UInt8('\'')
const _SEP = UInt8(',')
const _MISSING = UInt8('?')
const _ESC = UInt8('\\')
const _NOMSTART = UInt8('{')
const _NOMEND = UInt8('}')
const _ILLEGAL_RAWCHAR = (_SPACE, _AT, _COMMENT, _DQUO, _SQUO, _SEP, _MISSING, _ESC, _NOMSTART, _NOMEND)
const _A = (UInt8('a'), UInt8('A'))
const _B = (UInt8('b'), UInt8('B'))
const _C = (UInt8('c'), UInt8('C'))
const _D = (UInt8('d'), UInt8('D'))
const _E = (UInt8('e'), UInt8('E'))
const _G = (UInt8('g'), UInt8('G'))
const _I = (UInt8('i'), UInt8('I'))
const _L = (UInt8('l'), UInt8('L'))
const _M = (UInt8('m'), UInt8('M'))
const _N = (UInt8('n'), UInt8('N'))
const _O = (UInt8('o'), UInt8('O'))
const _R = (UInt8('r'), UInt8('R'))
const _S = (UInt8('s'), UInt8('S'))
const _T = (UInt8('t'), UInt8('T'))
const _U = (UInt8('u'), UInt8('U'))

mutable struct State
    src :: String
    pos :: Int
    lineno :: Int
end

Base.length(s::State) = sizeof(s.src)
Base.@propagate_inbounds Base.getindex(s::State, i::Int) = codeunit(s.src, i)
@inline hasmore(s::State) = s.pos ≤ length(s)
@inline inc!(s::State) = (s.pos += 1; nothing)
@inline curbyte(s::State) = s[s.pos]
@inline maybeskip!(s::State, c::UInt8) = @inbounds (hasmore(s) && curbyte(s) == c ? (inc!(s); true) : false)
@inline maybeskip!(s::State, cs::Tuple{Vararg{UInt8}}) = @inbounds (hasmore(s) && curbyte(s) in cs ? (inc!(s); true) : false)
@inline skip!(s::State, c::UInt8) = @inbounds (hasmore(s) && curbyte(s) == c ? (inc!(s); nothing) : _error_expecting(s, c))
@inline skip!(s::State, cs::Tuple{Vararg{UInt8}}) = @inbounds (hasmore(s) && curbyte(s) in cs ? (inc!(s); nothing) : _error_expecting(s, cs))
@inline skip!(s::State, c, cs...) = (skip!(s, c); skip!(s, cs...))

@noinline _error_expecting(s::State, c::UInt8) = _error(s, "expecting $(repr(Char(c)))")
@noinline _error_expecting(s::State, cs::Tuple{Vararg{UInt8}}) = _error(s, "expecting one of $(join(map(repr∘Char, cs), " "))")
@noinline _error_expecting(s::State, what::String) = _error(s, "expecting $what")
@noinline _error(s::State, msg) = error("line $(s.lineno):$(s.pos): parsing error: $msg")

skipspace!(s::State) =
    @inbounds while maybeskip!(s, _SPACE); end

skipcomment!(s::State) =
    @inbounds if maybeskip!(s, _COMMENT); s.pos = length(s)+1; end

function parse_type(s::State)
    if maybeskip!(s, _N)
        skip!(s, _U, _M, _E, _R, _I, _C)
        return ARFFNumericType()
    elseif maybeskip!(s, _R)
        skip!(s, _E, _A, _L)
        return ARFFNumericType()
    elseif maybeskip!(s, _I)
        skip!(s, _N, _T, _E, _G, _E, _R)
        return ARFFNumericType()
    elseif maybeskip!(s, _S)
        skip!(s, _T, _R, _I, _N, _G)
        return ARFFStringType()
    elseif maybeskip!(s, _D)
        skip!(s, _A, _T, _E)
        skipspace!(s)
        if hasmore(s) && (curbyte(s) in (_SQUO, _DQUO) || curbyte(s) ∉ _ILLEGAL_RAWCHAR)
            fmt = parse_string(s)
            return ARFFDateType(fmt)
        else
            return ARFFDateType()
        end
    elseif maybeskip!(s, _NOMSTART)
        skipspace!(s)
        xs = String[]
        while true
            x = parse_string(s)
            push!(xs, x)
            skipspace!(s)
            if maybeskip!(s, _SEP)
                skipspace!(s)
            elseif maybeskip!(s, _NOMEND)
                break
            else
                _error_expecting(s, (_SEP, _NOMEND))
            end
        end
        return ARFFNominalType(xs)
    else
        _error_expecting(s, "a type")
    end
end

function parse_string(s::State)
    io = IOBuffer()
    if maybeskip!(s, _DQUO)
        @inbounds while hasmore(s)
            c = curbyte(s)
            if c == _DQUO
                inc!(s)
                break
            elseif c == _ESC
                error("escape sequences not supported yet")
            else
                inc!(s)
                write(io, c)
            end
        end
    elseif maybeskip!(s, _SQUO)
        @inbounds while hasmore(s)
            c = curbyte(s)
            if c == _SQUO
                inc!(s)
                break
            elseif c == _ESC
                error("escapes sequences not supported yet")
            else
                inc!(s)
                write(io, c)
            end
        end
    else
        @inbounds while hasmore(s)
            c = curbyte(s)
            if c in _ILLEGAL_RAWCHAR
                break
            else
                write(io, c)
                inc!(s)
            end
        end
        io.size > 0 || _error_expecting(s, "a string")
    end
    String(take!(io))
end

function parse_header_line(s::State)
    skipspace!(s)
    r = nothing
    if maybeskip!(s, _AT)
        if maybeskip!(s, _D)
            skip!(s, _A, _T, _A)
            r = ARFFData()
        elseif maybeskip!(s, _A)
            skip!(s, _T, _T, _R, _I, _B, _U, _T, _E)
            skipspace!(s)
            name = parse_string(s)
            skipspace!(s)
            tp = parse_type(s)
            r = ARFFAttribute(name, tp)
        elseif maybeskip!(s, _R)
            skip!(s, _E, _L, _A, _T, _I, _O, _N)
            skipspace!(s)
            name = parse_string(s)
            r = ARFFRelation(name)
        else
            _error(s, "invalid attribute")
        end
        skipspace!(s)
    end
    skipcomment!(s)
    hasmore(s) && _error_expecting(s, r===nothing ? "command or end of string" : "end of string")
    return r
end

# Parse the ARFF header from `io`, stopping after `@data` is seen.
function parse_header(io::IO)
    relation = missing
    attributes = ARFFAttribute[]
    lineno = 0
    while !eof(io)
        lineno += 1
        h = parse_header_line(State(readline(io), 1, lineno))
        if h === nothing
            # blank line
        elseif h isa ARFFRelation
            relation === missing || error("line $lineno: @relation seen twice")
            relation = h.name
        elseif h isa ARFFAttribute
            relation === missing && error("line $lineno: @relation required before @attribute")
            push!(attributes, h)
        elseif h isa ARFFData
            relation === missing && error("line $lineno: @relation required before @data")
            return ARFFHeader(relation, attributes), lineno
        else
            error()
        end
    end
    return error("reached end of file before seeing @data")
end


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

@doc raw"""ResourceFieldSelector represents container resources (cpu, memory) and their output format

    IoK8sApiCoreV1ResourceFieldSelector(;
        containerName=nothing,
        divisor=nothing,
        resource=nothing,
    )

    - containerName::String : Container name: required for volumes, optional for env vars
    - divisor::IoK8sApimachineryPkgApiResourceQuantity : Specifies the output format of the exposed resources, defaults to \"1\"
    - resource::String : Required: resource to select
"""
mutable struct IoK8sApiCoreV1ResourceFieldSelector <: SwaggerModel
    containerName::Any # spec type: Union{ Nothing, String } # spec name: containerName
    divisor::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApiResourceQuantity } # spec name: divisor
    resource::Any # spec type: Union{ Nothing, String } # spec name: resource

    function IoK8sApiCoreV1ResourceFieldSelector(;containerName=nothing, divisor=nothing, resource=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ResourceFieldSelector, Symbol("containerName"), containerName)
        setfield!(o, Symbol("containerName"), containerName)
        validate_property(IoK8sApiCoreV1ResourceFieldSelector, Symbol("divisor"), divisor)
        setfield!(o, Symbol("divisor"), divisor)
        validate_property(IoK8sApiCoreV1ResourceFieldSelector, Symbol("resource"), resource)
        setfield!(o, Symbol("resource"), resource)
        o
    end
end # type IoK8sApiCoreV1ResourceFieldSelector

const _property_map_IoK8sApiCoreV1ResourceFieldSelector = Dict{Symbol,Symbol}(Symbol("containerName")=>Symbol("containerName"), Symbol("divisor")=>Symbol("divisor"), Symbol("resource")=>Symbol("resource"))
const _property_types_IoK8sApiCoreV1ResourceFieldSelector = Dict{Symbol,String}(Symbol("containerName")=>"String", Symbol("divisor")=>"IoK8sApimachineryPkgApiResourceQuantity", Symbol("resource")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1ResourceFieldSelector }) = collect(keys(_property_map_IoK8sApiCoreV1ResourceFieldSelector))
Swagger.property_type(::Type{ IoK8sApiCoreV1ResourceFieldSelector }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ResourceFieldSelector[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ResourceFieldSelector }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ResourceFieldSelector[property_name]

function check_required(o::IoK8sApiCoreV1ResourceFieldSelector)
    (getproperty(o, Symbol("resource")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ResourceFieldSelector }, name::Symbol, val)
end


text:
using CSV, DataFrames, CategoricalArrays

data = CSV.read("../data/purchaseData.csv", copycols = true)
println("Levels of Grade: ", levels(data.Grade))
println("Data points: ", nrow(data))

n = sum(.!(ismissing.(data.Grade)))
println("Non-missing data points: ", n)
data2 = dropmissing(data[:,[:Grade]],:Grade)

gradeInQuestion = "E"
indicatorVector = data2.Grade .== gradeInQuestion
numSuccess = sum(indicatorVector)
phat = numSuccess/n
serr = sqrt(phat*(1-phat)/n)

alpha = 0.05
confidencePercent = 100*(1-alpha)
zVal = quantile(Normal(),1-alpha/2)
confInt = (phat - zVal*serr, phat + zVal*serr)

println("\nOut of $n non-missing observations, "*
        "$numSuccess are at level $gradeInQuestion.")
println("Hence a point estimate for the proportion "*
        "of grades at level $gradeInQuestion is $phat.")
println("A $confidencePercent% confidence interval for "*
        "the proprotion of level $gradeInQuestion is:\n$confInt.")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#    API implementation for two body orbit propagator.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

get_epoch(orbp::OrbitPropagatorTwoBody) = orbp.tbd.epoch

function get_mean_elements(orbp::OrbitPropagatorTwoBody)
    tbd = orbp.tbd
    orb = KeplerianElements(
        tbd.epoch + tbd.Δt / 86400,
        tbd.a_0,
        tbd.e_0,
        tbd.i_0,
        tbd.Ω_0,
        tbd.ω_0,
        tbd.f_k
    )

    return orb
end

"""
    init_orbit_propagator(::Val{:twobody}, epoch::Number, a_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, f_0::Number; μ::T = m0) where T

Initialize the two body orbit propagator.

# Args

- `epoch::Number`: Initial orbit epoch [Julian Day].
- `a_0::Number`: Initial mean semi-major axis [m].
- `e_0::Number`: Initial mean eccentricity.
- `i_0::Number`: Initial mean inclination [rad].
- `Ω_0::Number`: Initial mean right ascension of the ascending node [rad].
- `ω_0::Number`: Initial mean argument of perigee [rad].
- `f_0::Number`: Initial mean true anomaly [rad].
- `orb_0::Orbit`: Object of type [`Orbit`](@ref) with the initial mean orbital
    elements [SI].

## Keywords

* `μ::Number`: (OPTIONAL) Standard gravitational parameter of the central body
    \\[m^3/s^2] (**Default** = `m0`).
"""
function init_orbit_propagator(
    ::Val{:twobody},
    epoch::Number,
    a_0::Number,
    e_0::Number,
    i_0::Number,
    Ω_0::Number,
    ω_0::Number,
    f_0::Number;
    μ::T = m0
) where T
    # Create the new Two Body propagator structure.
    tbd = twobody_init(epoch, a_0, e_0, i_0, Ω_0, ω_0, f_0; μ = μ)

    # Create and return the orbit propagator structure.
    return OrbitPropagatorTwoBody(tbd)
end

function init_orbit_propagator(::Val{:twobody}, orb_0::Orbit; μ::Number = m0)
    # Convert the orbit representation to Keplerian elements.
    k_0 = convert(KeplerianElements, orb_0)

    return init_orbit_propagator(
        Val(:twobody),
        k_0.t,
        k_0.a,
        k_0.e,
        k_0.i,
        k_0.Ω,
        k_0.ω,
        k_0.f;
        μ = μ
    )
end

function propagate!(orbp::OrbitPropagatorTwoBody, t::Number)
    # Auxiliary variables.
    tbd = orbp.tbd

    # Propagate the orbit.
    r_i, v_i = twobody!(tbd, t)

    return r_i, v_i
end

function step!(orbp::OrbitPropagatorTwoBody, Δt::Number)
    # Auxiliary variables.
    tbd = orbp.tbd

    # Propagate the orbit.
    r_i, v_i = twobody!(tbd, tbd.Δt + Δt)

    # Return the information about the step.
    return r_i, v_i
end


text:
using Random

xm = rand(10, 2)
ym = rand(10, 2)
v = Dataset(rand(10, 2))
@testset "Deprecations" begin
    # rqa types
    R = RecurrenceMatrix(xm[:,1], 0.5)
    rqa_params = rqa(R, onlydiagonal=true)
    @test_logs (:warn, "x.RR is deprecated for results of `rqa`; use x[:RR] instead") rqa_params.RR
    rqa_tuple = @test_logs (:warn, "RQA with `NamedTuple` output is deprecated, and will be removed in later versions") rqa(NamedTuple, R, onlydiagonal=true)
    @test rqa_params[:RR] == rqa_tuple.RR
    @test rqa_params[:DET] == rqa_tuple.DET
    @test rqa_params[:ENTR] == rqa_tuple.ENTR
    @test rqa_params[:L] == rqa_tuple.L
    # Matrix inputs
    @test_logs (:warn, "`RecurrenceMatrix{WithinRange}(x::AbstractMatrix, ε; kwargs...)` is deprecated, use `RecurrenceMatrix{WithinRange}(Dataset(x), ε; kwargs...)`") RecurrenceMatrix(xm, 0.5)
    @test_logs (:warn, "`CrossRecurrenceMatrix{WithinRange}(x::AbstractMatrix, y, ε; kwargs...)` is deprecated, use `CrossRecurrenceMatrix{WithinRange}(Dataset(x), y, ε; kwargs...)`") CrossRecurrenceMatrix(xm, v, 0.5)
    @test_logs (:warn, "`CrossRecurrenceMatrix{WithinRange}(x, y::AbstractMatrix, ε; kwargs...)` is deprecated, use `CrossRecurrenceMatrix{WithinRange}(x, Dataset(y), ε; kwargs...)`") CrossRecurrenceMatrix(v, ym, 0.5)
    @test_logs (:warn, "`CrossRecurrenceMatrix{WithinRange}(x::AbstractMatrix, y::AbstractMatrix, ε; kwargs...)` is deprecated, use `CrossRecurrenceMatrix{WithinRange}(Dataset(x), Dataset(y), ε; kwargs...)`") CrossRecurrenceMatrix(xm, ym, 0.5)
    @test_logs (:warn, "`JointRecurrenceMatrix{WithinRange}(x::AbstractMatrix, y, ε; kwargs...)` is deprecated, use `JointRecurrenceMatrix{WithinRange}(Dataset(x), y, ε; kwargs...)`") JointRecurrenceMatrix(xm, v, 0.5)
    @test_logs (:warn, "`JointRecurrenceMatrix{WithinRange}(x, y::AbstractMatrix, ε; kwargs...)` is deprecated, use `JointRecurrenceMatrix{WithinRange}(x, Dataset(y), ε; kwargs...)`") JointRecurrenceMatrix(v, ym, 0.5)
    @test_logs (:warn, "`JointRecurrenceMatrix{WithinRange}(x::AbstractMatrix, y::AbstractMatrix, ε; kwargs...)` is deprecated, use `JointRecurrenceMatrix{WithinRange}(Dataset(x), Dataset(y), ε; kwargs...)`") JointRecurrenceMatrix(xm, ym, 0.5)
    # transitivity
    @test_logs (:warn, "`transitivity(x::AbstractRecurrenceMatrix)` is deprecated`, use `rna` to analyse network parameters") transitivity(R)
end


text:
#=
breadth_first_search:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-10-02

search the nearest node from a given node

Tree definition
     1
    / \
   2   3
 / \  / \
4  5 6  7
=#

module BFS
    # node class
    mutable struct Node
        val
        left
        right
    end

    function search(root, target)
        # define empty queue
        node_que = []

        # add root node to queue
        push!(node_que, root)

        # iterate search until empty
        while length(node_que) > 0
            # get node from head of queue
            node = node_que[1]
            deleteat!(node_que, 1)

            if node != nothing
                println("Search Node $(node.val)")

                # found target
                if node.val == target
                    return true
                end

                # add child nodes to queue
                push!(node_que, node.left)
                push!(node_que, node.right)
            end
        end

        # not found target
        return false
    end

    function main()
        # difine tree
        root = Node(1, nothing, nothing)
        root.left = Node(2, nothing, nothing)
        root.right = Node(3, nothing, nothing)
        root.left.left = Node(4, nothing, nothing)
        root.left.right = Node(5, nothing, nothing)
        root.right.left = Node(6, nothing, nothing)
        root.right.right = Node(7, nothing, nothing)

        # target node
        target = 4

        # search target node
        result = search(root, target)
        if result == true
            println("Found target node $(target)")
        else
            println("Not found")
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .BFS
    BFS.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
# Gun: variant R2 (fully rational case; leja + repeated nodes + freeze)

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test

include(joinpath("..", "rk_helper", "gun_test_utils.jl"))

@bench @testset "NLEIGS: Gun variant R2" begin
    verbose = displaylevel

    nep, Σ, Ξ, v, nodes, funres = gun_init()

    # solve nlep
    lambda, X, res, solution_info = nleigs(nep, Σ, Ξ=Ξ, logger=verbose > 0 ? 1 : 0, minit=60, maxit=100, v=v, nodes=nodes, errmeasure=funres, return_details=verbose > 1)

    verify_lambdas(21, nep, lambda, X)

    if verbose > 1
        include("nleigs_residual_plot.jl")
        nleigs_residual_plot("Gun: variant R2", solution_info, Σ; ylims=[1e-17, 1e-1])
    end
end


text:
immutable GiteaError <: Exception
	code::Int64
	msg::Union{String,Vector{Dict{String,Any}},Dict{String,Any}}
end

function Base.show(io::IO,err::GiteaError)
	if isa(err.msg,String)
		msg = err.msg
	elseif isa(err.msg,Vector)
		msg = err.msg[1]["message"]
	else
		msg = err.msg["message"]
	end
	print(io, "Gitea Error:\tError Code: $(err.code)\tError Msg: \"$(JSON.json(err.msg))\"")

end


text:
@testset "6.1.3 (e x)^m (a+b sinh(c+d x^n))^p" begin
    (a, b, c, d, e, m, n, p, x, ) = @variables a b c d e m n p x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*sinh(c+d*x^n))^p=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*sinh(c+d*x^2))^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*(a+b*sinh(c+d*x^2))^p=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [x^3*sinh(a + b*x^2), x, 3, (x^2*cosh(a + b*x^2))/(2*b) - sinh(a + b*x^2)/(2*b^2)]
    @test_int [x^2*sinh(a + b*x^2), x, 4, (x*cosh(a + b*x^2))/(2*b) - (sqrt(pi)*erf(sqrt(b)*x))/(ℯ^a*(8*b^(3/2))) - (ℯ^a*sqrt(pi)*erfi(sqrt(b)*x))/(8*b^(3/2))]
    @test_int [x*sinh(a + b*x^2), x, 2, cosh(a + b*x^2)/(2*b)]
    @test_int [sinh(a + b*x^2), x, 3, -((sqrt(pi)*erf(sqrt(b)*x))/(ℯ^a*(4*sqrt(b)))) + (ℯ^a*sqrt(pi)*erfi(sqrt(b)*x))/(4*sqrt(b))]
    @test_int [sinh(a + b*x^2)/x, x, 3, (1/2)*CoshIntegral(b*x^2)*sinh(a) + (1/2)*cosh(a)*SinhIntegral(b*x^2)]
    @test_int [sinh(a + b*x^2)/x^2, x, 4, ((1/2)*sqrt(b)*sqrt(pi)*erf(sqrt(b)*x))/ℯ^a + (1/2)*sqrt(b)*ℯ^a*sqrt(pi)*erfi(sqrt(b)*x) - sinh(a + b*x^2)/x]
    @test_int [sinh(a + b*x^2)/x^3, x, 5, (1/2)*b*cosh(a)*CoshIntegral(b*x^2) - sinh(a + b*x^2)/(2*x^2) + (1/2)*b*sinh(a)*SinhIntegral(b*x^2)]

    @test_int [x^3*sinh(a + b*x^2)^2, x, 3, -(x^4/8) + (x^2*cosh(a + b*x^2)*sinh(a + b*x^2))/(4*b) - sinh(a + b*x^2)^2/(8*b^2)]
    @test_int [x^2*sinh(a + b*x^2)^2, x, 6, -(x^3/6) + (sqrt(pi/2)*erf(sqrt(2)*sqrt(b)*x))/(ℯ^(2*a)*(32*b^(3/2))) - (ℯ^(2*a)*sqrt(pi/2)*erfi(sqrt(2)*sqrt(b)*x))/(32*b^(3/2)) + (x*sinh(2*a + 2*b*x^2))/(8*b)]
    @test_int [x*sinh(a + b*x^2)^2, x, 3, -(x^2/4) + (cosh(a + b*x^2)*sinh(a + b*x^2))/(4*b)]
    @test_int [sinh(a + b*x^2)^2, x, 5, -(x/2) + (sqrt(pi/2)*erf(sqrt(2)*sqrt(b)*x))/(ℯ^(2*a)*(8*sqrt(b))) + (ℯ^(2*a)*sqrt(pi/2)*erfi(sqrt(2)*sqrt(b)*x))/(8*sqrt(b))]
    @test_int [sinh(a + b*x^2)^2/x, x, 5, (1/4)*cosh(2*a)*CoshIntegral(2*b*x^2) - log(x)/2 + (1/4)*sinh(2*a)*SinhIntegral(2*b*x^2)]
    @test_int [sinh(a + b*x^2)^2/x^2, x, 6, ((-(1/2))*sqrt(b)*sqrt(pi/2)*erf(sqrt(2)*sqrt(b)*x))/ℯ^(2*a) + (1/2)*sqrt(b)*ℯ^(2*a)*sqrt(pi/2)*erfi(sqrt(2)*sqrt(b)*x) - sinh(a + b*x^2)^2/x]
    @test_int [sinh(a + b*x^2)^2/x^3, x, 7, 1/(4*x^2) - cosh(2*(a + b*x^2))/(4*x^2) + (1/2)*b*CoshIntegral(2*b*x^2)*sinh(2*a) + (1/2)*b*cosh(2*a)*SinhIntegral(2*b*x^2)]

    @test_int [x^3*sinh(a + b*x^2)^3, x, 4, -((x^2*cosh(a + b*x^2))/(3*b)) + sinh(a + b*x^2)/(3*b^2) + (x^2*cosh(a + b*x^2)*sinh(a + b*x^2)^2)/(6*b) - sinh(a + b*x^2)^3/(18*b^2)]
    @test_int [x^2*sinh(a + b*x^2)^3, x, 10, -((3*x*cosh(a + b*x^2))/(8*b)) + (x*cosh(3*a + 3*b*x^2))/(24*b) + (3*sqrt(pi)*erf(sqrt(b)*x))/(ℯ^a*(32*b^(3/2))) - (sqrt(pi/3)*erf(sqrt(3)*sqrt(b)*x))/(ℯ^(3*a)*(96*b^(3/2))) + (3*ℯ^a*sqrt(pi)*erfi(sqrt(b)*x))/(32*b^(3/2)) - (ℯ^(3*a)*sqrt(pi/3)*erfi(sqrt(3)*sqrt(b)*x))/(96*b^(3/2))]
    @test_int [x*sinh(a + b*x^2)^3, x, 3, -(cosh(a + b*x^2)/(2*b)) + cosh(a + b*x^2)^3/(6*b)]
    @test_int [sinh(a + b*x^2)^3, x, 8, (3*sqrt(pi)*erf(sqrt(b)*x))/(ℯ^a*(16*sqrt(b))) - (sqrt(pi/3)*erf(sqrt(3)*sqrt(b)*x))/(ℯ^(3*a)*(16*sqrt(b))) - (3*ℯ^a*sqrt(pi)*erfi(sqrt(b)*x))/(16*sqrt(b)) + (ℯ^(3*a)*sqrt(pi/3)*erfi(sqrt(3)*sqrt(b)*x))/(16*sqrt(b))]
    @test_int [sinh(a + b*x^2)^3/x, x, 8, (-(3/8))*CoshIntegral(b*x^2)*sinh(a) + (1/8)*CoshIntegral(3*b*x^2)*sinh(3*a) - (3/8)*cosh(a)*SinhIntegral(b*x^2) + (1/8)*cosh(3*a)*SinhIntegral(3*b*x^2)]
    @test_int [sinh(a + b*x^2)^3/x^2, x, 9, ((-(3/8))*sqrt(b)*sqrt(pi)*erf(sqrt(b)*x))/ℯ^a + ((1/8)*sqrt(b)*sqrt(3*pi)*erf(sqrt(3)*sqrt(b)*x))/ℯ^(3*a) - (3/8)*sqrt(b)*ℯ^a*sqrt(pi)*erfi(sqrt(b)*x) + (1/8)*sqrt(b)*ℯ^(3*a)*sqrt(3*pi)*erfi(sqrt(3)*sqrt(b)*x) - sinh(a + b*x^2)^3/x]
    @test_int [sinh(a + b*x^2)^3/x^3, x, 12, (-(3/8))*b*cosh(a)*CoshIntegral(b*x^2) + (3/8)*b*cosh(3*a)*CoshIntegral(3*b*x^2) + (3*sinh(a + b*x^2))/(8*x^2) - sinh(3*(a + b*x^2))/(8*x^2) - (3/8)*b*sinh(a)*SinhIntegral(b*x^2) + (3/8)*b*sinh(3*a)*SinhIntegral(3*b*x^2)]

    @test_int [x*sinh(a + b*x^2)^7, x, 3, -(cosh(a + b*x^2)/(2*b)) + cosh(a + b*x^2)^3/(2*b) - (3*cosh(a + b*x^2)^5)/(10*b) + cosh(a + b*x^2)^7/(14*b)]

    #= ::Subsubsection:: =#
    #=p<0=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*sinh(c+d*x^2))^p*with*m*symbolic=#

    @test_int [(e*x)^m*sinh(a + b*x^2)^p, x, 0, Unintegrable((e*x)^m*sinh(a + b*x^2)^p, x)]

    @test_int [(e*x)^m*sinh(a + b*x^2)^3, x, 8, -((3^(-(1/2) - m/2)*ℯ^(3*a)*(e*x)^(1 + m)*((-b)*x^2)^((1/2)*(-1 - m))*Gamma((1 + m)/2, -3*b*x^2))/(16*e)) + (3*ℯ^a*(e*x)^(1 + m)*((-b)*x^2)^((1/2)*(-1 - m))*Gamma((1 + m)/2, (-b)*x^2))/(16*e) - (3*(e*x)^(1 + m)*(b*x^2)^((1/2)*(-1 - m))*Gamma((1 + m)/2, b*x^2))/(ℯ^a*(16*e)) + (3^(-(1/2) - m/2)*(e*x)^(1 + m)*(b*x^2)^((1/2)*(-1 - m))*Gamma((1 + m)/2, 3*b*x^2))/(ℯ^(3*a)*(16*e))]
    @test_int [(e*x)^m*sinh(a + b*x^2)^2, x, 5, -((e*x)^(1 + m)/(2*e*(1 + m))) - (2^(-(7/2) - m/2)*ℯ^(2*a)*(e*x)^(1 + m)*((-b)*x^2)^((1/2)*(-1 - m))*Gamma((1 + m)/2, -2*b*x^2))/e - (2^(-(7/2) - m/2)*(e*x)^(1 + m)*(b*x^2)^((1/2)*(-1 - m))*Gamma((1 + m)/2, 2*b*x^2))/(ℯ^(2*a)*e)]
    @test_int [(e*x)^m*sinh(a + b*x^2)^1, x, 3, -((ℯ^a*(e*x)^(1 + m)*((-b)*x^2)^((1/2)*(-1 - m))*Gamma((1 + m)/2, (-b)*x^2))/(4*e)) + ((e*x)^(1 + m)*(b*x^2)^((1/2)*(-1 - m))*Gamma((1 + m)/2, b*x^2))/(ℯ^a*(4*e))]
    @test_int [(e*x)^m/sinh(a + b*x^2)^1, x, 1, ((e*x)^m*Unintegrable(x^m*csch(a + b*x^2), x))/x^m]

    #= ::Section:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*sinh(c+d*x^3))^p=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*sinh(c+d*x^4))^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*(a+b*sinh(c+d*x^2))^p=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [x^3*sinh(a + b*x^4), x, 2, cosh(a + b*x^4)/(4*b)]

    #= ::Subsubsection:: =#
    #=p<0=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*sinh(c+d / x^1))^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*(a+b*sinh(c+d / x))^p=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [sinh(a + b/x)*x^2, x, 7, (1/6)*b*x^2*cosh(a + b/x) - (1/6)*b^3*cosh(a)*CoshIntegral(b/x) + (1/6)*b^2*x*sinh(a + b/x) + (1/3)*x^3*sinh(a + b/x) - (1/6)*b^3*sinh(a)*SinhIntegral(b/x)]
    @test_int [sinh(a + b/x)*x^1, x, 6, (1/2)*b*x*cosh(a + b/x) - (1/2)*b^2*CoshIntegral(b/x)*sinh(a) + (1/2)*x^2*sinh(a + b/x) - (1/2)*b^2*cosh(a)*SinhIntegral(b/x)]
    @test_int [sinh(a + b/x)*x^0, x, 5, (-b)*cosh(a)*CoshIntegral(b/x) + x*sinh(a + b/x) - b*sinh(a)*SinhIntegral(b/x)]
    @test_int [sinh(a + b/x)/x^1, x, 3, (-CoshIntegral(b/x))*sinh(a) - cosh(a)*SinhIntegral(b/x)]
    @test_int [sinh(a + b/x)/x^2, x, 2, -(cosh(a + b/x)/b)]
    @test_int [sinh(a + b/x)/x^3, x, 3, -(cosh(a + b/x)/(b*x)) + sinh(a + b/x)/b^2]
    @test_int [sinh(a + b/x)/x^4, x, 4, -((2*cosh(a + b/x))/b^3) - cosh(a + b/x)/(b*x^2) + (2*sinh(a + b/x))/(b^2*x)]
    @test_int [sinh(a + b/x)/x^5, x, 5, -(cosh(a + b/x)/(b*x^3)) - (6*cosh(a + b/x))/(b^3*x) + (6*sinh(a + b/x))/b^4 + (3*sinh(a + b/x))/(b^2*x^2)]

    #= ::Subsubsection:: =#
    #=p<0=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*sinh(c+d / x))^p*with*m*symbolic=#

    @test_int [(e*x)^m*sinh(a + b/x)^3, x, 9, (-(1/8))*3^(1 + m)*b*ℯ^(3*a)*(-(b/x))^m*(e*x)^m*Gamma(-1 - m, -((3*b)/x)) + (3/8)*b*ℯ^a*(-(b/x))^m*(e*x)^m*Gamma(-1 - m, -(b/x)) + ((3/8)*b*(b/x)^m*(e*x)^m*Gamma(-1 - m, b/x))/ℯ^a - ((1/8)*3^(1 + m)*b*(b/x)^m*(e*x)^m*Gamma(-1 - m, (3*b)/x))/ℯ^(3*a)]
    @test_int [(e*x)^m*sinh(a + b/x)^2, x, 6, -((x*(e*x)^m)/(2*(1 + m))) - 2^(-1 + m)*b*ℯ^(2*a)*(-(b/x))^m*(e*x)^m*Gamma(-1 - m, -((2*b)/x)) + (2^(-1 + m)*b*(b/x)^m*(e*x)^m*Gamma(-1 - m, (2*b)/x))/ℯ^(2*a)]
    @test_int [(e*x)^m*sinh(a + b/x)^1, x, 4, (-(1/2))*b*ℯ^a*(-(b/x))^m*(e*x)^m*Gamma(-1 - m, -(b/x)) - ((1/2)*b*(b/x)^m*(e*x)^m*Gamma(-1 - m, b/x))/ℯ^a]
    @test_int [(e*x)^m/sinh(a + b/x)^1, x, 1, ((e*x)^m*Unintegrable(x^m*csch(a + b/x), x))/x^m]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*sinh(c+d / x^2))^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*(a+b*sinh(c+d / x^2))^p=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [sinh(a + b/x^2)*x^4, x, 7, (2/15)*b*x^3*cosh(a + b/x^2) - ((2/15)*b^(5/2)*sqrt(pi)*erf(sqrt(b)/x))/ℯ^a - (2/15)*b^(5/2)*ℯ^a*sqrt(pi)*erfi(sqrt(b)/x) + (4/15)*b^2*x*sinh(a + b/x^2) + (1/5)*x^5*sinh(a + b/x^2)]
    @test_int [sinh(a + b/x^2)*x^3, x, 6, (1/4)*b*x^2*cosh(a + b/x^2) - (1/4)*b^2*CoshIntegral(b/x^2)*sinh(a) + (1/4)*x^4*sinh(a + b/x^2) - (1/4)*b^2*cosh(a)*SinhIntegral(b/x^2)]
    @test_int [sinh(a + b/x^2)*x^2, x, 6, (2/3)*b*x*cosh(a + b/x^2) + ((1/3)*b^(3/2)*sqrt(pi)*erf(sqrt(b)/x))/ℯ^a - (1/3)*b^(3/2)*ℯ^a*sqrt(pi)*erfi(sqrt(b)/x) + (1/3)*x^3*sinh(a + b/x^2)]
    @test_int [sinh(a + b/x^2)*x^1, x, 5, (-(1/2))*b*cosh(a)*CoshIntegral(b/x^2) + (1/2)*x^2*sinh(a + b/x^2) - (1/2)*b*sinh(a)*SinhIntegral(b/x^2)]
    @test_int [sinh(a + b/x^2)*x^0, x, 5, ((-(1/2))*sqrt(b)*sqrt(pi)*erf(sqrt(b)/x))/ℯ^a - (1/2)*sqrt(b)*ℯ^a*sqrt(pi)*erfi(sqrt(b)/x) + x*sinh(a + b/x^2)]
    @test_int [sinh(a + b/x^2)/x^1, x, 3, (-(1/2))*CoshIntegral(b/x^2)*sinh(a) - (1/2)*cosh(a)*SinhIntegral(b/x^2)]
    @test_int [sinh(a + b/x^2)/x^2, x, 4, (sqrt(pi)*erf(sqrt(b)/x))/(ℯ^a*(4*sqrt(b))) - (ℯ^a*sqrt(pi)*erfi(sqrt(b)/x))/(4*sqrt(b))]
    @test_int [sinh(a + b/x^2)/x^3, x, 2, -(cosh(a + b/x^2)/(2*b))]
    @test_int [sinh(a + b/x^2)/x^4, x, 5, -(cosh(a + b/x^2)/(2*b*x)) + (sqrt(pi)*erf(sqrt(b)/x))/(ℯ^a*(8*b^(3/2))) + (ℯ^a*sqrt(pi)*erfi(sqrt(b)/x))/(8*b^(3/2))]
    @test_int [sinh(a + b/x^2)/x^5, x, 3, -(cosh(a + b/x^2)/(2*b*x^2)) + sinh(a + b/x^2)/(2*b^2)]
    @test_int [sinh(a + b/x^2)/x^6, x, 6, -(cosh(a + b/x^2)/(2*b*x^3)) + (3*sqrt(pi)*erf(sqrt(b)/x))/(ℯ^a*(16*b^(5/2))) - (3*ℯ^a*sqrt(pi)*erfi(sqrt(b)/x))/(16*b^(5/2)) + (3*sinh(a + b/x^2))/(4*b^2*x)]
    @test_int [sinh(a + b/x^2)/x^7, x, 4, -(cosh(a + b/x^2)/b^3) - cosh(a + b/x^2)/(2*b*x^4) + sinh(a + b/x^2)/(b^2*x^2)]

    #= ::Subsubsection:: =#
    #=p<0=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*sinh(c+d / x))^p*with*m*symbolic=#

    @test_int [(e*x)^m*sinh(a + b/x^2)^3, x, 9, (1/16)*3^((1 + m)/2)*ℯ^(3*a)*(-(b/x^2))^((1 + m)/2)*x*(e*x)^m*Gamma((1/2)*(-1 - m), -((3*b)/x^2)) - (3/16)*ℯ^a*(-(b/x^2))^((1 + m)/2)*x*(e*x)^m*Gamma((1/2)*(-1 - m), -(b/x^2)) + ((3/16)*(b/x^2)^((1 + m)/2)*x*(e*x)^m*Gamma((1/2)*(-1 - m), b/x^2))/ℯ^a - ((1/16)*3^((1 + m)/2)*(b/x^2)^((1 + m)/2)*x*(e*x)^m*Gamma((1/2)*(-1 - m), (3*b)/x^2))/ℯ^(3*a)]
    @test_int [(e*x)^m*sinh(a + b/x^2)^2, x, 6, -((x*(e*x)^m)/(2*(1 + m))) + 2^((1/2)*(-5 + m))*ℯ^(2*a)*(-(b/x^2))^((1 + m)/2)*x*(e*x)^m*Gamma((1/2)*(-1 - m), -((2*b)/x^2)) + (2^((1/2)*(-5 + m))*(b/x^2)^((1 + m)/2)*x*(e*x)^m*Gamma((1/2)*(-1 - m), (2*b)/x^2))/ℯ^(2*a)]
    @test_int [(e*x)^m*sinh(a + b/x^2)^1, x, 4, (1/4)*ℯ^a*(-(b/x^2))^((1 + m)/2)*x*(e*x)^m*Gamma((1/2)*(-1 - m), -(b/x^2)) - ((1/4)*(b/x^2)^((1 + m)/2)*x*(e*x)^m*Gamma((1/2)*(-1 - m), b/x^2))/ℯ^a]
    @test_int [(e*x)^m/sinh(a + b/x^2)^1, x, 1, ((e*x)^m*Unintegrable(x^m*csch(a + b/x^2), x))/x^m]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*sinh(c+d*x^(1/2)))^p=#

    @test_int [sinh(sqrt(x))/sqrt(x), x, 2, 2*cosh(sqrt(x))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*sinh(c+d*x^n))^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*sinh(c+d*x^n))^p=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [x^2*sinh(a + b*x^n), x, 3, -((ℯ^a*x^3*Gamma(3/n, (-b)*x^n))/(((-b)*x^n)^(3/n)*(2*n))) + (x^3*Gamma(3/n, b*x^n))/(ℯ^a*(b*x^n)^(3/n)*(2*n))]
    @test_int [x^1*sinh(a + b*x^n), x, 3, -((ℯ^a*x^2*Gamma(2/n, (-b)*x^n))/(((-b)*x^n)^(2/n)*(2*n))) + (x^2*Gamma(2/n, b*x^n))/(ℯ^a*(b*x^n)^(2/n)*(2*n))]
    @test_int [x^0*sinh(a + b*x^n), x, 3, -((ℯ^a*x*Gamma(1/n, (-b)*x^n))/(((-b)*x^n)^n^(-1)*(2*n))) + (x*Gamma(1/n, b*x^n))/(ℯ^a*(b*x^n)^n^(-1)*(2*n))]
    @test_int [sinh(a + b*x^n)/x^1, x, 3, (CoshIntegral(b*x^n)*sinh(a))/n + (cosh(a)*SinhIntegral(b*x^n))/n]
    @test_int [sinh(a + b*x^n)/x^2, x, 3, -((ℯ^a*((-b)*x^n)^(1/n)*Gamma(-(1/n), (-b)*x^n))/(2*n*x)) + ((b*x^n)^(1/n)*Gamma(-(1/n), b*x^n))/(ℯ^a*(2*n*x))]
    @test_int [sinh(a + b*x^n)/x^3, x, 3, -((ℯ^a*((-b)*x^n)^(2/n)*Gamma(-(2/n), (-b)*x^n))/(2*n*x^2)) + ((b*x^n)^(2/n)*Gamma(-(2/n), b*x^n))/(ℯ^a*(2*n*x^2))]

    @test_int [x^2*sinh(a + b*x^n)^2, x, 5, -(x^3/6) - (2^(-2 - 3/n)*ℯ^(2*a)*x^3*Gamma(3/n, -2*b*x^n))/(((-b)*x^n)^(3/n)*n) - (2^(-2 - 3/n)*x^3*Gamma(3/n, 2*b*x^n))/(ℯ^(2*a)*(b*x^n)^(3/n)*n)]
    @test_int [x^1*sinh(a + b*x^n)^2, x, 5, -(x^2/4) - (4^(-1 - 1/n)*ℯ^(2*a)*x^2*Gamma(2/n, -2*b*x^n))/(((-b)*x^n)^(2/n)*n) - (4^(-1 - 1/n)*x^2*Gamma(2/n, 2*b*x^n))/(ℯ^(2*a)*(b*x^n)^(2/n)*n)]
    @test_int [x^0*sinh(a + b*x^n)^2, x, 5, -(x/2) - (2^(-2 - 1/n)*ℯ^(2*a)*x*Gamma(1/n, -2*b*x^n))/(((-b)*x^n)^n^(-1)*n) - (2^(-2 - 1/n)*x*Gamma(1/n, 2*b*x^n))/(ℯ^(2*a)*(b*x^n)^n^(-1)*n)]
    @test_int [sinh(a + b*x^n)^2/x^1, x, 5, (cosh(2*a)*CoshIntegral(2*b*x^n))/(2*n) - log(x)/2 + (sinh(2*a)*SinhIntegral(2*b*x^n))/(2*n)]
    @test_int [sinh(a + b*x^n)^2/x^2, x, 5, 1/(2*x) - (2^(-2 + 1/n)*ℯ^(2*a)*((-b)*x^n)^(1/n)*Gamma(-(1/n), -2*b*x^n))/(n*x) - (2^(-2 + 1/n)*(b*x^n)^(1/n)*Gamma(-(1/n), 2*b*x^n))/(ℯ^(2*a)*(n*x))]

    @test_int [x^2*sinh(a + b*x^n)^3, x, 8, -((ℯ^(3*a)*x^3*Gamma(3/n, -3*b*x^n))/(3^(3/n)*((-b)*x^n)^(3/n)*(8*n))) + (3*ℯ^a*x^3*Gamma(3/n, (-b)*x^n))/(((-b)*x^n)^(3/n)*(8*n)) - (3*x^3*Gamma(3/n, b*x^n))/(ℯ^a*(b*x^n)^(3/n)*(8*n)) + (x^3*Gamma(3/n, 3*b*x^n))/(3^(3/n)*ℯ^(3*a)*(b*x^n)^(3/n)*(8*n))]
    @test_int [x^1*sinh(a + b*x^n)^3, x, 8, -((ℯ^(3*a)*x^2*Gamma(2/n, -3*b*x^n))/(9^n^(-1)*((-b)*x^n)^(2/n)*(8*n))) + (3*ℯ^a*x^2*Gamma(2/n, (-b)*x^n))/(((-b)*x^n)^(2/n)*(8*n)) - (3*x^2*Gamma(2/n, b*x^n))/(ℯ^a*(b*x^n)^(2/n)*(8*n)) + (x^2*Gamma(2/n, 3*b*x^n))/(9^n^(-1)*ℯ^(3*a)*(b*x^n)^(2/n)*(8*n))]
    @test_int [x^0*sinh(a + b*x^n)^3, x, 8, -((ℯ^(3*a)*x*Gamma(1/n, -3*b*x^n))/(3^n^(-1)*((-b)*x^n)^n^(-1)*(8*n))) + (3*ℯ^a*x*Gamma(1/n, (-b)*x^n))/(((-b)*x^n)^n^(-1)*(8*n)) - (3*x*Gamma(1/n, b*x^n))/(ℯ^a*(b*x^n)^n^(-1)*(8*n)) + (x*Gamma(1/n, 3*b*x^n))/(3^n^(-1)*ℯ^(3*a)*(b*x^n)^n^(-1)*(8*n))]
    @test_int [sinh(a + b*x^n)^3/x^1, x, 8, -((3*CoshIntegral(b*x^n)*sinh(a))/(4*n)) + (CoshIntegral(3*b*x^n)*sinh(3*a))/(4*n) - (3*cosh(a)*SinhIntegral(b*x^n))/(4*n) + (cosh(3*a)*SinhIntegral(3*b*x^n))/(4*n)]
    @test_int [sinh(a + b*x^n)^3/x^2, x, 8, -((3^(1/n)*ℯ^(3*a)*((-b)*x^n)^(1/n)*Gamma(-(1/n), -3*b*x^n))/(8*n*x)) + (3*ℯ^a*((-b)*x^n)^(1/n)*Gamma(-(1/n), (-b)*x^n))/(8*n*x) - (3*(b*x^n)^(1/n)*Gamma(-(1/n), b*x^n))/(ℯ^a*(8*n*x)) + (3^(1/n)*(b*x^n)^(1/n)*Gamma(-(1/n), 3*b*x^n))/(ℯ^(3*a)*(8*n*x))]

    #= ::Subsubsection:: =#
    #=p<0=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*(a+b*sinh(c+d*x^n))^p*with*m*symbolic=#

    @test_int [(e*x)^m*(b*sinh(c + d*x^n))^p, x, 0, Unintegrable((e*x)^m*(b*sinh(c + d*x^n))^p, x)]
    @test_int [(e*x)^m*(a + b*sinh(c + d*x^n))^p, x, 0, Unintegrable((e*x)^m*(a + b*sinh(c + d*x^n))^p, x)]

    @test_int [(e*x)^(n - 1)*(b*sinh(c + d*x^n))^p, x, 3, ((e*x)^n*cosh(c + d*x^n)*HypergeometricFunctions._₂F₁(1/2, (1 + p)/2, (3 + p)/2, -sinh(c + d*x^n)^2)*(b*sinh(c + d*x^n))^(1 + p))/(x^n*(b*d*e*n*(1 + p)*sqrt(cosh(c + d*x^n)^2)))]
    @test_int [(e*x)^(2*n - 1)*(b*sinh(c + d*x^n))^p, x, 1, ((e*x)^(2*n)*Unintegrable(x^(-1 + 2*n)*(b*sinh(c + d*x^n))^p, x))/(x^(2*n)*e)]

    @test_int [(e*x)^(n - 1)*(a + b*sinh(c + d*x^n))^p, x, 5, (im*sqrt(2)*(e*x)^n*AppellF1(1/2, 1/2, -p, 3/2, (1/2)*(1 - im*sinh(c + d*x^n)), (b*(1 - im*sinh(c + d*x^n)))/(im*a + b))*cosh(c + d*x^n)*(a + b*sinh(c + d*x^n))^p)/(x^n*((a + b*sinh(c + d*x^n))/(a - im*b))^p*(d*e*n*sqrt(1 + im*sinh(c + d*x^n))))]
    @test_int [(e*x)^(2*n - 1)*(a + b*sinh(c + d*x^n))^p, x, 1, ((e*x)^(2*n)*Unintegrable(x^(-1 + 2*n)*(a + b*sinh(c + d*x^n))^p, x))/(x^(2*n)*e)]

    @test_int [(e*x)^m*sinh(a + b*x^n)^3, x, 8, -((ℯ^(3*a)*(e*x)^(1 + m)*Gamma((1 + m)/n, -3*b*x^n))/(3^((1 + m)/n)*((-b)*x^n)^((1 + m)/n)*(8*e*n))) + (3*ℯ^a*(e*x)^(1 + m)*Gamma((1 + m)/n, (-b)*x^n))/(((-b)*x^n)^((1 + m)/n)*(8*e*n)) - (3*(e*x)^(1 + m)*Gamma((1 + m)/n, b*x^n))/(ℯ^a*(b*x^n)^((1 + m)/n)*(8*e*n)) + ((e*x)^(1 + m)*Gamma((1 + m)/n, 3*b*x^n))/(3^((1 + m)/n)*ℯ^(3*a)*(b*x^n)^((1 + m)/n)*(8*e*n))]
    @test_int [(e*x)^m*sinh(a + b*x^n)^2, x, 5, -((e*x)^(1 + m)/(2*e*(1 + m))) - (ℯ^(2*a)*(e*x)^(1 + m)*Gamma((1 + m)/n, -2*b*x^n))/(2^((1 + m + 2*n)/n)*((-b)*x^n)^((1 + m)/n)*(e*n)) - ((e*x)^(1 + m)*Gamma((1 + m)/n, 2*b*x^n))/(2^((1 + m + 2*n)/n)*ℯ^(2*a)*(b*x^n)^((1 + m)/n)*(e*n))]
    @test_int [(e*x)^m*sinh(a + b*x^n)^1, x, 3, -((ℯ^a*(e*x)^(1 + m)*Gamma((1 + m)/n, (-b)*x^n))/(((-b)*x^n)^((1 + m)/n)*(2*e*n))) + ((e*x)^(1 + m)*Gamma((1 + m)/n, b*x^n))/(ℯ^a*(b*x^n)^((1 + m)/n)*(2*e*n))]
    @test_int [(e*x)^m/sinh(a + b*x^n)^2, x, 1, ((e*x)^m*Unintegrable(x^m*csch(a + b*x^n)^2, x))/x^m]

    @test_int [sinh(a + b*x^n)^1/x^(n + 1), x, 5, (b*cosh(a)*CoshIntegral(b*x^n))/n - sinh(a + b*x^n)/(x^n*n) + (b*sinh(a)*SinhIntegral(b*x^n))/n]
    @test_int [sinh(a + b*x^n)^2/x^(n + 1), x, 7, 1/(x^n*(2*n)) - cosh(2*(a + b*x^n))/(x^n*(2*n)) + (b*CoshIntegral(2*b*x^n)*sinh(2*a))/n + (b*cosh(2*a)*SinhIntegral(2*b*x^n))/n]
    @test_int [sinh(a + b*x^n)^3/x^(n + 1), x, 12, -((3*b*cosh(a)*CoshIntegral(b*x^n))/(4*n)) + (3*b*cosh(3*a)*CoshIntegral(3*b*x^n))/(4*n) + (3*sinh(a + b*x^n))/(x^n*(4*n)) - sinh(3*(a + b*x^n))/(x^n*(4*n)) - (3*b*sinh(a)*SinhIntegral(b*x^n))/(4*n) + (3*b*sinh(3*a)*SinhIntegral(3*b*x^n))/(4*n)]

    @test_int [x^(n/2 - 1)*sinh(a + b*x^n), x, 4, -((sqrt(pi)*erf(sqrt(b)*x^(n/2)))/(ℯ^a*(2*sqrt(b)*n))) + (ℯ^a*sqrt(pi)*erfi(sqrt(b)*x^(n/2)))/(2*sqrt(b)*n)]

    #= ::Title:: =#
    #=Integrands*of*the*form*(e*x)^m*sinh(a+b*(c+d*x)^n)=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*sinh(a+b*(c+d*x)^n)=#

    @test_int [x^2*sinh((a + b*x)^2), x, 12, -((a*cosh((a + b*x)^2))/b^3) + ((a + b*x)*cosh((a + b*x)^2))/(2*b^3) - (sqrt(pi)*erf(a + b*x))/(8*b^3) - (a^2*sqrt(pi)*erf(a + b*x))/(4*b^3) - (sqrt(pi)*erfi(a + b*x))/(8*b^3) + (a^2*sqrt(pi)*erfi(a + b*x))/(4*b^3)]
    @test_int [x^1*sinh((a + b*x)^2), x, 8, cosh((a + b*x)^2)/(2*b^2) + (a*sqrt(pi)*erf(a + b*x))/(4*b^2) - (a*sqrt(pi)*erfi(a + b*x))/(4*b^2)]
    @test_int [x^0*sinh((a + b*x)^2), x, 4, -((sqrt(pi)*erf(a + b*x))/(4*b)) + (sqrt(pi)*erfi(a + b*x))/(4*b)]
    @test_int [sinh((a + b*x)^2)/x^1, x, 1, b*CannotIntegrate(sinh((a + b*x)^2)/(b*x), x)]
    @test_int [sinh((a + b*x)^2)/x^2, x, 1, Unintegrable(sinh((a + b*x)^2)/x^2, x), b^2*CannotIntegrate(sinh((a + b*x)^2)/(b^2*x^2), x)]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*sinh(a+b*(c+d*x)^(n/2))=#

    @test_int [sinh(a + b*sqrt(c + d*x))*x^2, x, 16, (240*sqrt(c + d*x)*cosh(a + b*sqrt(c + d*x)))/(b^5*d^3) - (24*c*sqrt(c + d*x)*cosh(a + b*sqrt(c + d*x)))/(b^3*d^3) + (2*c^2*sqrt(c + d*x)*cosh(a + b*sqrt(c + d*x)))/(b*d^3) + (40*(c + d*x)^(3/2)*cosh(a + b*sqrt(c + d*x)))/(b^3*d^3) - (4*c*(c + d*x)^(3/2)*cosh(a + b*sqrt(c + d*x)))/(b*d^3) + (2*(c + d*x)^(5/2)*cosh(a + b*sqrt(c + d*x)))/(b*d^3) - (240*sinh(a + b*sqrt(c + d*x)))/(b^6*d^3) + (24*c*sinh(a + b*sqrt(c + d*x)))/(b^4*d^3) - (2*c^2*sinh(a + b*sqrt(c + d*x)))/(b^2*d^3) - (120*(c + d*x)*sinh(a + b*sqrt(c + d*x)))/(b^4*d^3) + (12*c*(c + d*x)*sinh(a + b*sqrt(c + d*x)))/(b^2*d^3) - (10*(c + d*x)^2*sinh(a + b*sqrt(c + d*x)))/(b^2*d^3)]
    @test_int [sinh(a + b*sqrt(c + d*x))*x^1, x, 10, (12*sqrt(c + d*x)*cosh(a + b*sqrt(c + d*x)))/(b^3*d^2) - (2*c*sqrt(c + d*x)*cosh(a + b*sqrt(c + d*x)))/(b*d^2) + (2*(c + d*x)^(3/2)*cosh(a + b*sqrt(c + d*x)))/(b*d^2) - (12*sinh(a + b*sqrt(c + d*x)))/(b^4*d^2) + (2*c*sinh(a + b*sqrt(c + d*x)))/(b^2*d^2) - (6*(c + d*x)*sinh(a + b*sqrt(c + d*x)))/(b^2*d^2)]
    @test_int [sinh(a + b*sqrt(c + d*x))*x^0, x, 4, (2*sqrt(c + d*x)*cosh(a + b*sqrt(c + d*x)))/(b*d) - (2*sinh(a + b*sqrt(c + d*x)))/(b^2*d)]
    @test_int [sinh(a + b*sqrt(c + d*x))/x^1, x, 10, CoshIntegral(b*(sqrt(c) + sqrt(c + d*x)))*sinh(a - b*sqrt(c)) + CoshIntegral(b*(sqrt(c) - sqrt(c + d*x)))*sinh(a + b*sqrt(c)) - cosh(a + b*sqrt(c))*SinhIntegral(b*(sqrt(c) - sqrt(c + d*x))) + cosh(a - b*sqrt(c))*SinhIntegral(b*(sqrt(c) + sqrt(c + d*x)))]
    @test_int [sinh(a + b*sqrt(c + d*x))/x^2, x, 11, (b*d*cosh(a + b*sqrt(c))*CoshIntegral(b*(sqrt(c) - sqrt(c + d*x))))/(2*sqrt(c)) - (b*d*cosh(a - b*sqrt(c))*CoshIntegral(b*(sqrt(c) + sqrt(c + d*x))))/(2*sqrt(c)) - sinh(a + b*sqrt(c + d*x))/x - (b*d*sinh(a + b*sqrt(c))*SinhIntegral(b*(sqrt(c) - sqrt(c + d*x))))/(2*sqrt(c)) - (b*d*sinh(a - b*sqrt(c))*SinhIntegral(b*(sqrt(c) + sqrt(c + d*x))))/(2*sqrt(c))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(e*x)^m*sinh(a+b*(c+d*x)^(n/3))=#

    @test_int [sinh(a + b*(c + d*x)^(1/3))*x^2, x, 23, (120960*cosh(a + b*(c + d*x)^(1/3)))/(b^9*d^3) + (6*c^2*cosh(a + b*(c + d*x)^(1/3)))/(b^3*d^3) - (720*c*(c + d*x)^(1/3)*cosh(a + b*(c + d*x)^(1/3)))/(b^5*d^3) + (60480*(c + d*x)^(2/3)*cosh(a + b*(c + d*x)^(1/3)))/(b^7*d^3) + (3*c^2*(c + d*x)^(2/3)*cosh(a + b*(c + d*x)^(1/3)))/(b*d^3) - (120*c*(c + d*x)*cosh(a + b*(c + d*x)^(1/3)))/(b^3*d^3) + (5040*(c + d*x)^(4/3)*cosh(a + b*(c + d*x)^(1/3)))/(b^5*d^3) - (6*c*(c + d*x)^(5/3)*cosh(a + b*(c + d*x)^(1/3)))/(b*d^3) + (168*(c + d*x)^2*cosh(a + b*(c + d*x)^(1/3)))/(b^3*d^3) + (3*(c + d*x)^(8/3)*cosh(a + b*(c + d*x)^(1/3)))/(b*d^3) + (720*c*sinh(a + b*(c + d*x)^(1/3)))/(b^6*d^3) - (120960*(c + d*x)^(1/3)*sinh(a + b*(c + d*x)^(1/3)))/(b^8*d^3) - (6*c^2*(c + d*x)^(1/3)*sinh(a + b*(c + d*x)^(1/3)))/(b^2*d^3) + (360*c*(c + d*x)^(2/3)*sinh(a + b*(c + d*x)^(1/3)))/(b^4*d^3) - (20160*(c + d*x)*sinh(a + b*(c + d*x)^(1/3)))/(b^6*d^3) + (30*c*(c + d*x)^(4/3)*sinh(a + b*(c + d*x)^(1/3)))/(b^2*d^3) - (1008*(c + d*x)^(5/3)*sinh(a + b*(c + d*x)^(1/3)))/(b^4*d^3) - (24*(c + d*x)^(7/3)*sinh(a + b*(c + d*x)^(1/3)))/(b^2*d^3)]
    @test_int [sinh(a + b*(c + d*x)^(1/3))*x^1, x, 13, -((6*c*cosh(a + b*(c + d*x)^(1/3)))/(b^3*d^2)) + (360*(c + d*x)^(1/3)*cosh(a + b*(c + d*x)^(1/3)))/(b^5*d^2) - (3*c*(c + d*x)^(2/3)*cosh(a + b*(c + d*x)^(1/3)))/(b*d^2) + (60*(c + d*x)*cosh(a + b*(c + d*x)^(1/3)))/(b^3*d^2) + (3*(c + d*x)^(5/3)*cosh(a + b*(c + d*x)^(1/3)))/(b*d^2) - (360*sinh(a + b*(c + d*x)^(1/3)))/(b^6*d^2) + (6*c*(c + d*x)^(1/3)*sinh(a + b*(c + d*x)^(1/3)))/(b^2*d^2) - (180*(c + d*x)^(2/3)*sinh(a + b*(c + d*x)^(1/3)))/(b^4*d^2) - (15*(c + d*x)^(4/3)*sinh(a + b*(c + d*x)^(1/3)))/(b^2*d^2)]
    @test_int [sinh(a + b*(c + d*x)^(1/3))*x^0, x, 5, (6*cosh(a + b*(c + d*x)^(1/3)))/(b^3*d) + (3*(c + d*x)^(2/3)*cosh(a + b*(c + d*x)^(1/3)))/(b*d) - (6*(c + d*x)^(1/3)*sinh(a + b*(c + d*x)^(1/3)))/(b^2*d)]
    @test_int [sinh(a + b*(c + d*x)^(1/3))/x^1, x, 13, CoshIntegral(b*(c^(1/3) - (c + d*x)^(1/3)))*sinh(a + b*c^(1/3)) + CoshIntegral(b*((-1)^(1/3)*c^(1/3) + (c + d*x)^(1/3)))*sinh(a - (-1)^(1/3)*b*c^(1/3)) + CoshIntegral((-b)*((-1)^(2/3)*c^(1/3) - (c + d*x)^(1/3)))*sinh(a + (-1)^(2/3)*b*c^(1/3)) - cosh(a + b*c^(1/3))*SinhIntegral(b*(c^(1/3) - (c + d*x)^(1/3))) - cosh(a + (-1)^(2/3)*b*c^(1/3))*SinhIntegral(b*((-1)^(2/3)*c^(1/3) - (c + d*x)^(1/3))) + cosh(a - (-1)^(1/3)*b*c^(1/3))*SinhIntegral(b*((-1)^(1/3)*c^(1/3) + (c + d*x)^(1/3)))]
    @test_int [sinh(a + b*(c + d*x)^(1/3))/x^2, x, 14, (b*d*cosh(a + b*c^(1/3))*CoshIntegral(b*(c^(1/3) - (c + d*x)^(1/3))))/(3*c^(2/3)) + ((-1)^(2/3)*b*d*cosh(a + (-1)^(2/3)*b*c^(1/3))*CoshIntegral((-b)*((-1)^(2/3)*c^(1/3) - (c + d*x)^(1/3))))/(3*c^(2/3)) - ((-1)^(1/3)*b*d*cosh(a - (-1)^(1/3)*b*c^(1/3))*CoshIntegral(b*((-1)^(1/3)*c^(1/3) + (c + d*x)^(1/3))))/(3*c^(2/3)) - sinh(a + b*(c + d*x)^(1/3))/x - (b*d*sinh(a + b*c^(1/3))*SinhIntegral(b*(c^(1/3) - (c + d*x)^(1/3))))/(3*c^(2/3)) - ((-1)^(2/3)*b*d*sinh(a + (-1)^(2/3)*b*c^(1/3))*SinhIntegral(b*((-1)^(2/3)*c^(1/3) - (c + d*x)^(1/3))))/(3*c^(2/3)) - ((-1)^(1/3)*b*d*sinh(a - (-1)^(1/3)*b*c^(1/3))*SinhIntegral(b*((-1)^(1/3)*c^(1/3) + (c + d*x)^(1/3))))/(3*c^(2/3))]
    end

max_stars_repo_path: test/6 Hyperbolic functions/6.1 Hyperbolic sine/6.1.3 (e x)^m (a+b sinh(c+d x^n))^p.jl


text:
module PaddleChainRules

using DLPack
using Requires

import FillArrays
import Adapt
import Functors: fmap
import ChainRulesCore

struct PyAdaptor{T} end
Adapt.adapt_storage(to::PyAdaptor{T}, x::AbstractArray) where {T} = convert(Array, x)
Adapt.adapt_storage(to::PyAdaptor{T}, x::StridedArray) where {T} = x
Adapt.adapt_storage(to::PyAdaptor{<:AbstractArray}, x::FillArrays.AbstractFill) = collect(x)
# Handle views
Adapt.adapt_structure(to::PyAdaptor{<:AbstractArray}, x::A) where {A <: SubArray}  = collect(x)

### XXX: what's a little piracy between us
fmap(f, x::ChainRulesCore.Tangent) = fmap(f, x.backing)

include("paddle.jl")

function __init__()
    @require CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" begin
        include("cuda.jl")
    end
end

end # module

max_stars_repo_name: PangHua/InfiniteScience

max_issues_repo_name: PangHua/InfiniteScience

max_forks_repo_name: PangHua/InfiniteScience


text:
module UStatsTest

using Test
using Random, LinearAlgebra, Statistics
using DiffPrecTest: compute, variance,
                    createVarianceStatistic, createCovarianceStatistic, createSecondMomentStatistic

@testset "one sample" begin
  x = randn(10)
  @test compute(createVarianceStatistic(), x) ≈ var(x)

  s = createSecondMomentStatistic()
  @test compute(s, x) ≈ sum(abs2, x) / 10
  @test variance(s, x) ≈ var( x .* x)

  X = randn(10, 2)
  @test compute(createCovarianceStatistic(1, 2), X) ≈ cov(X[:, 1], X[:, 2])

  s = createSecondMomentStatistic(1, 2)
  @test compute(s, X) ≈ mean(X[:, 1] .* X[:, 2])
  @test variance(s, X) ≈ var(X[:, 1] .* X[:, 2])
end

end


text:
module CryptoDashApp

export run_app

using AlphaVantage, DataFrames, Dates, PlotlyJS, Dash, DashHtmlComponents,
      DashCoreComponents, CSV, Statistics, StatsBase, GLM, HTTP, JSON, Query,
      DelimitedFiles

include("CryptoFunctions.jl")
include("GetDataFunctions.jl")
include("PlotFunctions.jl")
include("ConfigureApp.jl")

end # module


text:
# High level model construction

function cplex_model(env::Env;    # solver environment
    name::ASCIIString="",          # model name
    sense::Symbol=:Min,       # :minimize or :maximize
    H::CoeffMat=emptyfmat,         # quadratic coefficient matrix
    f::FVec=emptyfvec,             # linear coefficient vector
    A::CoeffMat=emptyfmat,         # LHS of inequality constraints
    b::FVec=emptyfvec,             # RHS of inequality constraints
    Aeq::CoeffMat=emptyfmat,       # LHS of equality constraints
    beq::FVec=emptyfvec,           # RHS of equality constraints
    lb::Bounds=-Inf,               # upper bounds
    ub::Bounds=Inf)                # lower bounds

  # check f
  if isempty(f)
    error("f must be specified.")
  end

    # create model
    model = Model(env, name)

    # set sense
    set_sense!(model, sense)

    # add variables
    add_vars!(model, f, lb, ub)

    # add qpterms
    if !isempty(H)
        add_qpterms!(model, H)
    end

    # add constraints
    if !isempty(A) && !isempty(b)
        add_constrs!(model, A, '<', b)
    end

    if !isempty(Aeq) && !isempty(beq)
        add_constrs!(model, Aeq, '=', beq)
    end
    return model
end


text:
using KFEstimate
using LinearAlgebra, Plots, Revise, ProgressBars
using Flux, Flux.Optimise
using Flux: Params, gradient
pathof(KFEstimate)
global dt = 0.001
g = 9.81

# non linear transition functions must always take x and u as variables regardless if they are used
function f(x, u)
    return [x[1] + x[2]*dt; x[2] - g*dt*sin(x[1])]
end

# non linear observation function only takes x as an input variable
function h(x)
    return [sin(x[1]), sin(x[2])]
end

Q = 0.1*[dt^3/3 dt^2/2; dt^2/2 dt]
R = 1.0*Matrix{Float64}(I, 2, 2)

ekf = ExtendedKalmanFilter(f, Q, h, R);
# run simulation
time_step = 0.0:5000
x0 = [pi/2, 0.0]
P0 = Matrix{Float64}(I, 2, 2)
action_sequence = [[0.0; 0.0] for t in time_step]
sim_states, sim_measurements = run_simulation(ekf, x0, action_sequence)
# run kalman filter
s0 = State(x0, P0) # initial state belief
filtered_states = run_filter(ekf, s0, action_sequence, sim_measurements)

# unpack sim and filtered states
μ, Σ = unpack(filtered_states)
θ = hcat([x[1] for x in sim_states], [x[2] for x in sim_states])
# plot
plot(time_step, θ[2:end, 1], label = "simulated θ", legend=:bottomleft)
plot!(time_step, θ[2:end, 2], label = "simulated dθ", legend=:bottomleft)
plot!(time_step, μ[2:end, 1], label = "filterd θ", legend=:bottomleft)
plot!(time_step, μ[2:end, 2], label = "filterd dθ", legend=:bottomleft)
xlabel!("time step (t)")

##

function fhat(θ, x, u)
    return [x[1] + x[2]*dt; x[2] - θ[1]*dt*sin(x[1])]
end

function hhat(θ, x)
    return [sin(x[1]), sin(x[2])]
end

Rhat(θ) = R
Qhat(θ) = Q

param_ekf = ExtendedParamKalmanFilter(fhat, Qhat, hhat, Rhat)

opt = ADAM(0.001)
epochs = 400

θ0 = [g+randn()]
println(θ0)

newθ, loss = run_ekf_gradient(θ0, param_ekf, s0, action_sequence, sim_measurements, opt, epochs)

grad_states = run_param_filter(newθ, param_ekf, s0, action_sequence, sim_measurements)
μgrad, Σgrad = unpack(grad_states)

l = @layout [a{0.7h};grid(1, 3)]
p1 = plot(time_step, [θ[2:end, 1] θ[2:end, 2]], label = ["simulated θ" "simulated dθ"], legend=:bottomright, xlabel="time step (t)")
p1 = plot!(time_step, [μ[2:end, 1] μ[2:end, 2] ], label = ["filtered θ" "filtered dθ"], legend=:bottomright, xlabel="time step (t)")
p1 = plot!(time_step, [μgrad[2:end, 1] μgrad[2:end, 2]], label = ["learned θ" "learned dθ"], legend=:bottomright, xlabel="time step (t)")
p2 = plot(1:epochs, loss, title="loss", xlabel="number of epochs")
p3 = plot(time_step, (θ[2:end, :]-μ[2:end, :]).^2, title="KF vs. sim error", xlabel="time step (t)")
p4 = plot(time_step, (θ[2:end, :]-μgrad[2:end, :]).^2, title="grad vs. sim error", xlabel="time step (t)")
plot(p1, p2, p3, p4, layout=l, titlefont = font(12), size=(1000, 700))


text:
function upair(r, data :: Data)

	@unpack sig6, eps4 = data

	r6 = (r^2)^3

	return eps4*(sig6/r6 - 1.0)/r6
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct CategoricalMapping <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CategoricalMapping(; kwargs...)
        obj = new(meta(CategoricalMapping), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct CategoricalMapping
const __meta_CategoricalMapping = Ref{ProtoMeta}()
function meta(::Type{CategoricalMapping})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CategoricalMapping)
            __meta_CategoricalMapping[] = target = ProtoMeta(CategoricalMapping)
            fnum = Int[1,2,101,102]
            allflds = Pair{Symbol,Union{Type,String}}[:stringToInt64Map => StringToInt64Map, :int64ToStringMap => Int64ToStringMap, :strValue => AbstractString, :int64Value => Int64]
            oneofs = Int[1,1,2,2]
            oneof_names = Symbol[Symbol("MappingType"),Symbol("ValueOnUnknown")]
            meta(target, CategoricalMapping, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_CategoricalMapping[]
    end
end
function Base.getproperty(obj::CategoricalMapping, name::Symbol)
    if name === :stringToInt64Map
        return (obj.__protobuf_jl_internal_values[name])::StringToInt64Map
    elseif name === :int64ToStringMap
        return (obj.__protobuf_jl_internal_values[name])::Int64ToStringMap
    elseif name === :strValue
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :int64Value
        return (obj.__protobuf_jl_internal_values[name])::Int64
    else
        getfield(obj, name)
    end
end

export CategoricalMapping


text:
using Decimals
using Compat.Test

@testset "Equality" begin

@testset "isequal" begin
    @test isequal(Decimal(0, 2, -3), Decimal(0, 2, -3))
    @test ~isequal(Decimal(0, 2, -3), Decimal(0, 2, 3))
    @test isequal(Decimal(0, 2, -3), 0.002)
    @test isequal(Decimal(1, 2, 0), -2)
    @test ~isequal(Decimal(1, 2, 0), 2)
end

@testset "==" begin
    @test Decimal(0, 2, -3) == Decimal(0, 2, -3)
    @test Decimal(0, 2, -3) != Decimal(0, 2, 3)
    @test Decimal(0, 2, -3) == 0.002

    @test -2 == Decimal(1, 2, 0)
    @test 2 != Decimal(1, 2, 0)

    @test Decimal(1, 2, 0) == -2
    @test Decimal(1, 2, 0) != 2

    bf_pi = BigFloat(pi)
    @test Decimal(bf_pi) == bf_pi
    @test bf_pi == Decimal(bf_pi)

    bi = big"4608230166434464229556241992703"
    @test Decimal(bi) == bi
    @test bi == Decimal(bi)
end

@testset "<" begin
    @test Decimal(1, 1, 1) < Decimal(0, 1, 1)
    @test !(Decimal(0, 1, 1) < Decimal(1, 1, 1))
    @test Decimal(1, 1, 1) < Decimal(1, 0, 1)
    @test !(Decimal(1, 0, 1) < Decimal(1, 1, 1))
    @test Decimal(0, 2, -3) < Decimal(0, 2, 3)
    @test !(Decimal(0, 2, 3) < Decimal(0, 2, -3))
end

end


text:
using KernelAbstractions: NoneEvent
using OffsetArrays: OffsetArray
using CUDA

using Oceananigans.Operators: Δzᵃᵃᶜ
using Oceananigans.BoundaryConditions: left_gradient, right_gradient, linearly_extrapolate, FBC, VBC, GBC
using Oceananigans.BoundaryConditions: fill_bottom_halo!, fill_top_halo!, apply_z_bottom_bc!, apply_z_top_bc!
using Oceananigans.Grids: Flat, Bounded
using Oceananigans.Coriolis: AbstractRotation
using Oceananigans.Architectures: device_event
using Oceananigans.TurbulenceClosures: AbstractTurbulenceClosure
using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: _top_tke_flux, CATKEVDArray

import Oceananigans.Utils: launch!
import Oceananigans.Grids: validate_size, validate_halo
import Oceananigans.BoundaryConditions: fill_halo_regions!
import Oceananigans.TurbulenceClosures: time_discretization, calculate_diffusivities!
import Oceananigans.TurbulenceClosures: ∂ⱼ_τ₁ⱼ, ∂ⱼ_τ₂ⱼ, ∂ⱼ_τ₃ⱼ, ∇_dot_qᶜ
import Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: top_tke_flux
import Oceananigans.Coriolis: x_f_cross_U, y_f_cross_U, z_f_cross_U

#####
##### Implements a "single column model mode" for HydrostaticFreeSurfaceModel
#####

const SingleColumnGrid = AbstractGrid{<:AbstractFloat, <:Flat, <:Flat, <:Bounded}

#####
##### Model constructor utils
#####

PressureField(::SingleColumnGrid) = (; pHY′ = nothing)
FreeSurface(free_surface::ExplicitFreeSurface{Nothing}, velocities, arch, ::SingleColumnGrid) = nothing
FreeSurface(free_surface::ImplicitFreeSurface{Nothing}, velocities, arch, ::SingleColumnGrid) = nothing

validate_momentum_advection(momentum_advection, ::SingleColumnGrid) = nothing
validate_tracer_advection(tracer_advection::AbstractAdvectionScheme, ::SingleColumnGrid) = nothing, NamedTuple()
validate_tracer_advection(tracer_advection::Nothing, ::SingleColumnGrid) = nothing, NamedTuple()

#####
##### Time-step optimizations
#####

calculate_free_surface_tendency!(arch, ::SingleColumnGrid, args...) = NoneEvent()

# Fast state update and halo filling

function update_state!(model::HydrostaticFreeSurfaceModel, grid::SingleColumnGrid)

    fill_halo_regions!(prognostic_fields(model), model.architecture, model.clock, fields(model))

    compute_auxiliary_fields!(model.auxiliary_fields)

    # Calculate diffusivities
    calculate_diffusivities!(model.diffusivity_fields, model.closure, model)

    fill_halo_regions!(model.diffusivity_fields,
                       model.architecture,
                       model.clock,
                       fields(model))

    return nothing
end

const ClosureArray = AbstractArray{<:AbstractTurbulenceClosure}

@inline function ∂ⱼ_τ₁ⱼ(i, j, k, grid::SingleColumnGrid, closure_array::ClosureArray, args...)
    @inbounds closure = closure_array[i, j]
    return ∂ⱼ_τ₁ⱼ(i, j, k, grid, closure, args...)
end

@inline function ∂ⱼ_τ₂ⱼ(i, j, k, grid::SingleColumnGrid, closure_array::ClosureArray, args...)
    @inbounds closure = closure_array[i, j]
    return ∂ⱼ_τ₂ⱼ(i, j, k, grid, closure, args...)
end

@inline function ∇_dot_qᶜ(i, j, k, grid::SingleColumnGrid, closure_array::ClosureArray, c, tracer_index, args...)
    @inbounds closure = closure_array[i, j]
    return ∇_dot_qᶜ(i, j, k, grid, closure, c, tracer_index, args...)
end

struct ColumnEnsembleSize{C<:Tuple{Int, Int}}
    ensemble :: C
    Nz :: Int
    Hz :: Int
end

ColumnEnsembleSize(; Nz, ensemble=(0, 0), Hz=1) = ColumnEnsembleSize(ensemble, Nz, Hz)

validate_size(TX, TY, TZ, e::ColumnEnsembleSize) = tuple(e.ensemble[1], e.ensemble[2], e.Nz)
validate_halo(TX, TY, TZ, e::ColumnEnsembleSize) = tuple(0, 0, e.Hz)

@inline function time_discretization(closure_array::AbstractArray)
    first_closure = CUDA.@allowscalar first(closure_array) # assumes all closures have same time-discretization
    return time_discretization(first_closure)
end

#####
##### CATKEVerticalDiffusivity helpers
#####

@inline tracer_tendency_kernel_function(model::HydrostaticFreeSurfaceModel, closure::CATKEVDArray, ::Val{:e}) =
    hydrostatic_turbulent_kinetic_energy_tendency

""" Compute the flux of TKE through the surface / top boundary. """
@inline function top_tke_flux(i, j, grid::SingleColumnGrid, clock, fields, parameters, closure_array::CATKEVDArray, buoyancy)
    top_tracer_bcs = parameters.top_tracer_boundary_conditions
    top_velocity_bcs = parameters.top_velocity_boundary_conditions
    @inbounds closure = closure_array[i, j]

    return _top_tke_flux(i, j, grid, closure.surface_TKE_flux, closure,
                         buoyancy, fields, top_tracer_bcs, top_velocity_bcs, clock)
end

@inline function hydrostatic_turbulent_kinetic_energy_tendency(i, j, k, grid::SingleColumnGrid,
                                                               val_tracer_index::Val{tracer_index},
                                                               advection,
                                                               closure_array::CATKEVDArray, args...) where tracer_index

    @inbounds closure = closure_array[i, j]
    return hydrostatic_turbulent_kinetic_energy_tendency(i, j, k, grid, val_tracer_index, advection, closure, args...)
end

#####
##### Arrays of Coriolises
#####

const CoriolisMatrix = AbstractMatrix{<:AbstractRotation}

@inline function x_f_cross_U(i, j, k, grid::SingleColumnGrid, coriolis_array::CoriolisMatrix, U)
    @inbounds coriolis = coriolis_array[i, j]
    return x_f_cross_U(i, j, k, grid, coriolis, U)
end

@inline function y_f_cross_U(i, j, k, grid::SingleColumnGrid, coriolis_array::CoriolisMatrix, U)
    @inbounds coriolis = coriolis_array[i, j]
    return y_f_cross_U(i, j, k, grid, coriolis, U)
end

@inline function z_f_cross_U(i, j, k, grid::SingleColumnGrid, coriolis_array::CoriolisMatrix, U)
    @inbounds coriolis = coriolis_array[i, j]
    return z_f_cross_U(i, j, k, grid, coriolis, U)
end


text:
@doc raw"""
    struct RadialBasis

Wrapper type for the `helfem::atomic::basis::RadialBasis` class, representing a finite
element basis for radial functions on the domain ``r \in [0, r_{\rm{max}}]``.
"""
struct RadialBasis
    b :: helfem.RadialBasis
    nnodes :: Int
    primbas :: Int
    rmax :: Float64
    igrid :: Int
    zexp :: Float64
    nquad :: Int

    function RadialBasis(nnodes, nelem; primbas = 4, rmax = 40.0, igrid = 4, zexp = 2.0, nquad = nothing)
        nquad = isnothing(nquad) ? 0 : nquad
        new(
            helfem.basis(nnodes, nelem, primbas, rmax, igrid, zexp, nquad),
            nnodes, primbas, rmax, igrid, zexp, nquad,
        )
    end

    function RadialBasis(pb::PolynomialBasis, grid::AbstractVector, nquad::Integer)
        b = helfem.RadialBasis(CxxPtr(pb.pb), nquad, helfem.ArmaVector(grid))
        nnodes = helfem.pb_order(pb.pb)
        nelem = length(grid) - 1
        primbas = pb.primbas
        rmax = last(grid)
        igrid = 0 # let's declare the igrid == 0 means a custom grid
        zexp = NaN
        new(b, nnodes, primbas, rmax, igrid, zexp, nquad)
    end
end

function Base.getproperty(b::RadialBasis, s::Symbol)
    if s === :nelem
        helfem.nel(b.b)
    else
        getfield(b, s)
    end
end

Base.propertynames(::RadialBasis) = (fieldnames(RadialBasis)..., :nelem)

# Special RadialBasis constructor that allows you to construct a new RadialBasis based on an
# existing one, where you have added (or, also, removed) elements without changing the
# values of the first elements.
#
# TODO: add tests for this
function HelFEM.RadialBasis(b::RadialBasis; nelem)
    b_new = HelFEM.RadialBasis(b.nnodes, nelem;
        primbas = b.primbas, igrid = b.igrid, zexp = b.zexp, nquad = b.nquad,
        # calculated from the new nelem value:
        rmax = ^(1 + b.rmax, (nelem/b.nelem)^b.zexp) - 1,
    )
    imax = min(nelem, b.nelem) + 1
    @assert HelFEM.boundaries(b)[1:imax] ≈ HelFEM.boundaries(b_new)[1:imax]
    return b_new
end

Base.length(b::RadialBasis) = Int(helfem.nbf(b.b))

function Base.show(io::IO, b::RadialBasis)
    nbf = helfem.nbf(b.b)
    nel = helfem.nel(b.b)
    write(io, "RadialBasis: $nbf basis functions ($nel elements)")
end

nquad(b::RadialBasis) = helfem.get_nquad(b.b)

function overlap(b::RadialBasis; invh=false)
    S = helfem.overlap(b.b, b.b)
    if invh
        collect(S), collect(HelFEM.invh(S))
    else
        collect(S)
    end
end

function radial_integral(b1::RadialBasis, n::Integer, b2::RadialBasis = b1; lderivative=false, rderivative=false)
    I = helfem.radial_integral(b1.b, b2.b, n, lderivative, rderivative)
    collect(I)
end

function potential(b1::RadialBasis, model::Symbol, rms::Real = 0.0, b2::RadialBasis=b1)
    if model === :point && rms != 0.0
        @warn "Point nuclear model should pass rms = 0.0. RMS value ignored." rms
    elseif model !== :point && rms <= 0.0
        throw(ArgumentError("Must pass rms >= 0 for nuclear model $model (passed $rms)"))
    end
    modelobj = (model === :point) ? helfem.PointNucleus(1) :
        (model === :gaussian)     ? helfem.GaussianNucleus(1, rms) :
        (model === :spherical)    ? helfem.SphericalNucleus(1, rms) :
        (model === :hollow)       ? helfem.HollowNucleus(1, rms) :
        throw(ArgumentError("Invalid nucler model: $(model)"))
    I = helfem.model_potential(b1.b, b2.b, CxxPtr(modelobj), false, false)
    collect(I)
end

boundaries(b::RadialBasis) = collect(helfem.get_bval(b.b))
add_boundary!(b::RadialBasis, r::Real) = helfem.add_boundary(b.b, r)

function quadraturepoints(b::RadialBasis)
    nq, nel = nquad(b), Int(helfem.nel(b.b))
    rs = Vector{Float64}(undef, nel*nq)
    for i = 0:nel-1
        el_rs = collect(helfem.get_r(b.b, i))
        @assert length(el_rs) == nq
        rs[nq*i+1:nq*(i+1)] .= el_rs
    end
    return rs
end

function quadratureweights(b::RadialBasis)
    nq, nel = nquad(b), Int(helfem.nel(b.b))
    rs = Vector{Float64}(undef, nel*nq)
    for i = 0:nel-1
        el_rs = collect(helfem.get_wrad(b.b, i))
        @assert length(el_rs) == nq
        rs[nq*i+1:nq*(i+1)] .= el_rs
    end
    return rs
end

function basisvalues(b::RadialBasis)
    nq, nel, nbf = nquad(b), Int(helfem.nel(b.b)), length(b)
    ys = zeros(nel*nq, nbf)
    nbf_count = 0
    for i = 0:nel-1
        el_ys = collect(helfem.get_bf(b.b, i))
        @assert size(el_ys, 1) == nq
        #@info "el=$i" size(el_ys)
        nbf_in_element = size(el_ys, 2)
        rs_range = nq*i+1 : nq*(i+1)
        bf_range_start = nbf_count == 0 ? 1 : nbf_count
        bf_range = bf_range_start:(bf_range_start+nbf_in_element-1)
        #@show bf_range
        ys[rs_range, bf_range] .= el_ys
        nbf_count = last(bf_range)
    end
    @assert nbf_count == nbf
    return ys
end

"""
Returns a vector of grid boundaries.
"""
function radialgrid(gridtype, nelem, rmax; zexp=nothing)
    igrid = (gridtype == :linear) ? 1 :
        (gridtype == :quadratic) ? 2 :
        (gridtype == :polynomial) ? 3 :
        (gridtype == :exponential) ? 4 :
        throw(ArgumentError("Invalid gridtype $gridtype"))
    if (gridtype == :linear || gridtype == :quadratic) && !isnothing(zexp)
        throw(ArgumentError("zexp has no effect on linear or quadratic grids"))
    end
    if (gridtype == :polynomial || gridtype == :exponential) && isnothing(zexp)
        throw(ArgumentError("Must set zexp for polynomial and exponential grids"))
    end
    collect(helfem.get_grid(rmax, nelem, igrid, isnothing(zexp) ? 0.0 : zexp))
end

function (b::RadialBasis)(rs)
    element_boundaries = boundaries(b)
    @assert minimum(element_boundaries) == 0
    @assert minimum(element_boundaries) <= minimum(rs)
    @assert maximum(element_boundaries) >= maximum(rs)
    # Group the rs values by by element. Note that we'll use the C++ indexing convention,
    # i.e. first element is labelled with 0.
    bf = zeros(Float64, (length(rs), length(b)))
    pb_ptr = helfem.get_poly(b.b)
    poly_nbf = helfem.get_nbf(pb_ptr)
    for iel = 0:(b.nelem - 1)
        rmin, rmax = element_boundaries[iel+1], element_boundaries[iel+2]
        r0, rλ = (rmax + rmin) / 2, (rmax - rmin) / 2
        # We assign the r values on the boundaries to the element on the left. In other words,
        # we'll define the element to be r ∈ (rmin, rmax]. This excludes r=0, but that is
        # intentional --  we need to special case that because of the 1/r scaling.
        idxs = findall(r -> rmin < r <= rmax, rs)
        isempty(idxs) && continue
        # Scale the rs values in this element down to the [-1, 1] range
        xs = (rs[idxs] .- r0) ./ rλ
        # Calculate the basis values in the element
        ys = collect(helfem.pb_eval(pb_ptr, helfem.ArmaVector(xs)))
        # Assign the calculated basis values to the correct place in the basis matrix
        bfrange_start = (iel == 0) ? 1 : (poly_nbf - 1) * iel
        bfrange_end = (iel == b.nelem - 1) ? length(b) : (poly_nbf - 1) * (iel + 1)
        ysrange_start = (iel == 0) ? 2 : 1
        ysrange_end = (iel == b.nelem - 1) ? (poly_nbf - 1) : poly_nbf
        # We also divide all the basis values by r, because that is how RadialBasis is defined
        # in HelFEM.
        bf[idxs, bfrange_start:bfrange_end] .= ys[:, ysrange_start:ysrange_end] ./ rs[idxs]
    end
    # Handle the r = 0 case, for which we need to look at the derivative of the polynomial
    # on the edge. That is, if f(0) = 0, then f(x)/x|_{x=0} = f'(0).
    idxs = findall(isequal(0), rs)
    if length(idxs) > 0
        zero_r = helfem.ArmaVector([-1]) # r = 0, therefore x = -1 in the first element
        ys, dys = helfem.ArmaMatrix(1, poly_nbf), helfem.ArmaMatrix(1, poly_nbf)
        helfem.pb_eval(pb_ptr, zero_r, ys, dys)
        dys = collect(dys)
        # If any of the r values are zero, we know that they will belong to the first
        # element, so we need to remove the leftmost basis function.
        for i in idxs
            # The 2/element_boundaries[2] factor comes from the scaling of the derivative
            # due to the scaling of the x-axis when going from [-1, 1] to [0, r_1].
            bf[i, 1:(poly_nbf - 1)] .= dys[1, 2:end] .* (2/element_boundaries[2])
        end
    end
    return bf
end

"""
    elementrange(b::RadialBasis, k::Integer) -> (r_k, r_{k+1})

Returns a tuple with the start and end boundary of the `k`-th element.
"""
function elementrange(b::RadialBasis, k::Integer)
    element_boundaries = boundaries(b)
    @assert 1 <= k < length(element_boundaries) # TODO: throw proper error
    element_boundaries[k], element_boundaries[k+1]
end

"""
    scale_to_element(b::RadialBasis, k::Integer, xs)

Scales the ``x`` values within the element `k` to the corresponding ``r`` coordinates.
"""
function scale_to_element(b::RadialBasis, k::Integer, xs)
    rmin, rmax = elementrange(b, k)
    r0, rλ = (rmax + rmin) / 2, (rmax - rmin) / 2
    xs .* rλ .+ r0
end

function controlpoints(b::RadialBasis)
    # We only know how to calculate control points for the LIP basis currently
    b.primbas == 4 || throw(ArgumentError("Can only calculate control points for LIP basis"))
    # The LIPBasis uses Gauss-Lobatto nodes, so the following only applies to LIP
    xs, _ = lobatto(b.nnodes) # we only need the nodes, can discard the weights
    rs = Vector{Float64}(undef, (b.nnodes - 1)*b.nelem + 1)
    for k in 1:b.nelem
        rs[(b.nnodes-1)*(k-1)+1:(b.nnodes-1)*k+1] .= scale_to_element(b, k, xs)
    end
    return rs
end


text:
using Random
using LinearAlgebra
using Statistics

mutable struct CMAES
    N::Int
    μ::Int
    λ::Int
    τ::Float64
    τ_c::Float64
    τ_σ::Float64
    population::Array{Array{Float64}}
    offspring::Array{Array{Float64}}
    F_μ::Array{Float64}
    F_λ::Array{Float64}
    C::Array{Float64}
    s::Array{Float64}
    s_σ::Array{Float64}
    σ::Float64
    E::Array{Float64}
    W::Array{Float64}
    x::Array{Float64}
end

function CMAES(;N=2, μ=1, λ=10, τ=sqrt(N), τ_c=N^2, τ_σ=sqrt(N))
    x = randn(N)
    population = fill(x, µ)
    offspring = Array{Array{Float64}}(undef, λ)
    F_µ = Inf .* ones(µ)
    F_λ = Inf .* ones(λ)
    C = Array(Diagonal{Float64}(I, N))
    s = zeros(N)
    s_σ = zeros(N)
    σ = 1.0
    E = zeros(N, λ)
    W = zeros(N, λ);
    CMAES(N, μ, λ, τ, τ_c, τ_σ, population, offspring, F_µ, F_λ, C, s, s_σ, σ, E, W, x)
end

function step!(c::CMAES, obj::Function)
    # L1
    sqrt_c = cholesky((c.C + c.C') / 2.0).U
    for i in 1:c.λ
        c.E[:,i] = randn(c.N)
        c.W[:,i] = c.σ * (sqrt_c * c.E[:,i])
        c.offspring[i] = c.x + c.W[:,i]
        c.F_λ[i] = obj(c.offspring[i])
    end
    # Select new parent population
    idx = sortperm(c.F_λ)[1:c.μ]
    for i in 1:c.μ
        c.population[i] = c.offspring[idx[i]]
        c.F_μ[i] = c.F_λ[idx[i]]
    end
    # L2
    w = vec(mean(c.W[:,idx], dims=2))
    c.x += w
    # L3
    c.s = (1.0 - 1.0/c.τ)*c.s + (sqrt(c.μ/c.τ * (2.0 - 1.0/c.τ))/c.σ)*w
    # L4
    c.C = (1.0 - 1.0/c.τ_c).*c.C + (c.s./c.τ_c)*c.s'
    # L5
    ɛ = vec(mean(c.E[:,idx], dims=2))
    c.s_σ = (1.0 - 1.0/c.τ_σ)*c.s_σ + sqrt(c.μ/c.τ_σ*(2.0 - 1.0/c.τ_σ))*ɛ
    # L6
    c.σ = c.σ*exp(((c.s_σ'*c.s_σ)[1] - c.N)/(2*c.N*sqrt(c.N)))
    c
end

max_stars_repo_name: alban-deBourges/evolution

max_issues_repo_name: alban-deBourges/evolution

max_forks_repo_name: alban-deBourges/evolution


text:
module HPhiJulia
    export HPhi
    currentd = pwd()
    const models = ["Fermion Hubbard", "Spin", "Kondo Lattice", "Fermion HubbardGC","SpinGC", "Kondo LatticeGC", "SpinGCCMA"]
    const methods = ["Lanczos", "TPQ", "Full Diag", "CG", "Time-Evolution"]
    const lattices = ["Chain Lattice", "Square Lattice", "Triangular Lattice", "Honeycomb Lattice", "Ladder", "Kagome"]

    cd(homedir()*"/.julia/packages")
    if isdir("HPhiJulia")
    else
        mkdir("HPhiJulia")
    end
    cd("HPhiJulia")

    const HPhipath = homedir()*"/.julia/packages/HPhiJulia/HPhi.build/src/HPhi"
    const HPhiversion = "3.2.0"
#    const HPhiversion = "3.1.2"
    const HPhitar = "HPhi-"*HPhiversion*".tar.gz"
    const HPhiurl = "https://github.com/issp-center-dev/HPhi/releases/download/v"*HPhiversion*"/"

    if isfile(HPhipath) == false
        using CMake
        if isfile(HPhitar) == false
            println("Downloading HPhi...")
            download( HPhiurl*HPhitar,HPhitar)
	        println("done.")

	    end
        run(`tar xzvf $HPhitar`)
        if isdir("HPhi.build") == false
       	    mkdir( "HPhi.build" )
        end
        cd("HPhi.build")
        run(`$cmake ../HPhi-$HPhiversion/`)
        run(`make`)
    end
    cd(currentd)

    function create_stanin(model,lattice,W,L,twoSz,method,
        mu,t,tp,U,V,Vp,J,Jp,nelec)
       fp = open("stan.in","w")
       println(fp,"model = \"$model\"")
       println(fp,"method = \"$method\"")
       println(fp,"lattice = \"$lattice\"")
       if lattice != "Chain Lattice"
              println(fp,"W= $W")
       end
       println(fp,"L = $L")

       if mu != 0 println(fp,"mu = $mu") end
       if t != 0 println(fp,"t = $t") end
       if tp != 0 println(fp,"t` = $tp") end
       if U != 0 println(fp,"U = $U") end
       if V != 0 println(fp,"V = $V") end
       if Vp != 0 println(fp,"V` = $Vp") end
       if J != 0 println(fp,"J = $J") end
       if Jp != 0 println(fp,"J` = $Jp") end
       if nelec != 0 println(fp,"nelec = $nelec")
              if model == "Fermion HubbardGC" || model == "Spin" || model == "SpinGC"
                     error("nelec should not be used in model: $model .")
              end
       end
       if model == "Fermion HubbardGC" || model == "SpinGC"
       else
              println(fp,"2Sz = $twoSz")
       end
       close(fp)
    end

    function read_calcmod()
       calcmod = readlines("calcmod.def")
       Calcmod = Dict(
       "CalcType" => parse(Int8,split(calcmod[5])[2]),
       "CalcModel" => parse(Int8,split(calcmod[6])[2]),
       "ReStart" => parse(Int8,split(calcmod[7])[2]),
       "CalcSpec" => parse(Int8,split(calcmod[8])[2]),
       "CalcEigenVec" => parse(Int8,split(calcmod[9])[2]),
       "InitialVecType" => parse(Int8,split(calcmod[10])[2]),
       "InputEigenVec" => parse(Int8,split(calcmod[11])[2]),
       "OutputEigenVec" => parse(Int8,split(calcmod[12])[2])
       )
       return Calcmod
    end

    function write_calcmod(Calcmod)
       fp = open("calcmod.def","w")
       println(fp,"#CalcType = 0:Lanczos, 1:TPQCalc, 2:FullDiag, 3:CG, 4:Time-evolution")
       println(fp,"#CalcModel = 0:Hubbard, 1:Spin, 2:Kondo, 3:HubbardGC, 4:SpinGC, 5:KondoGC")
       println(fp,"#Restart = 0:None, 1:Save, 2:Restart&Save, 3:Restart")
       println(fp,"#CalcSpec = 0:None, 1:Normal, 2:No H*Phi, 3:Save, 4:Restart, 5:Restart&Save")
       println(fp,"CalcType   $(Calcmod["CalcType"])")
       println(fp,"CalcModel   $(Calcmod["CalcModel"])")
       println(fp,"CalcSpec   $(Calcmod["CalcSpec"])")
       println(fp,"CalcEigenVec   $(Calcmod["CalcEigenVec"])")
       println(fp,"InitialVecType   $(Calcmod["InitialVecType"])")
       println(fp,"InputEigenVec   $(Calcmod["InputEigenVec"])")
       println(fp,"OutputEigenVec   $(Calcmod["OutputEigenVec"])")
       close(fp)

       return
    end

    function HPhi(model,lattice,W,L;expart=false,twoSz=0,method="CG",mpinum = 1,
        mu = 0,t=0,tp=0,U=0,V=0,Vp=0,J=0,Jp=0,nelec=0,OutputEigenVec=false,usenamelist=false)

        if expart == false
              if  OutputEigenVec
                     error("export should be true when OutputEigenVec = true")
              end
              create_stanin(model,lattice,W,L,twoSz,method,mu,t,tp,U,V,Vp,J,Jp,nelec)
              if mpinum == 1
                     run(`$HPhipath -s stan.in`)
              else
                     run(`mpirun -np $mpinum $HPhipath -s stan.in`)
              end
        else
              if usenamelist == false
                     create_stanin(model,lattice,W,L,twoSz,method,mu,t,tp,U,V,Vp,J,Jp,nelec)
                     run(`$HPhipath -sdry stan.in`)
              end
              Calcmod = read_calcmod()
              if  OutputEigenVec
                     Calcmod["OutputEigenVec"] = 1
              end
              write_calcmod(Calcmod)
              if mpinum == 1
                     run(`$HPhipath -e namelist.def`)
              else
                     run(`mpirun -np $mpinum $HPhipath -e namelist.def`)
              end
        end

        zvo_energy = readlines("./output/zvo_energy.dat")

        Energy = parse(Float64,split(zvo_energy[2])[2])
        Doublon = parse(Float64,split(zvo_energy[3])[2])
        Sz = parse(Float64,split(zvo_energy[4])[2])
        println("Energy is $Energy")
        println("Doublon is $Doublon")
        println("Sz is $Sz")
        results = Dict("Energy" => Energy,"Doublon" => Doublon,"Sz" => Sz)

        if OutputEigenVec
            fname = "output/zvo_eigenvec_0_rank_0.dat"
            println("loading $fname")
            evec,Hilbert = load_evec(fname)
            return results,evec
       else
              return results
        end

    end

    function load_evec(fname)
       io = open(fname)
       Lanczos = read(io, Int32)
       Hilbert= read(io, Int64)
       discard= read(io, ComplexF64)
       evec=zeros(ComplexF64, Hilbert)
       for i in 1:Hilbert
              v = read(io, ComplexF64)
              evec[i]= v
       end
       println("Eigenvector was loaded")
       return evec,Hilbert
    end

end

max_stars_repo_name: cometscome/HPhiJulia


text:
@testset "Ambiguity Resolution" begin
    @testset "fill" begin
        k = fill(Leaf(Tape(), 1.2), 3, 4)
        @test k.pullback !== nothing

        x = fill(Leaf(Tape(), 1.2), (3, 4))
        @test x.pullback !== nothing
    end
end


text:
__precompile__()

module PDE

using DataStructures: OrderedDict

# types
export Model, ModelState
export DiagnosticEquation, PrognosticEquation

# model functions
export iterate, iterate_heaps

# time stepping routines
export euler_forward

# operator

include("equations.jl")
include("model.jl")
include("integrator.jl")
include("operator.jl")

end


text:
# The policy type
mutable struct ValueIterationPolicy{F} <: Policy
    qmat::AbstractMatrix{F} # Q matrix storing Q(s,a) values
    util::AbstractVector{F} # The value function V(s)
    policy::AbstractVector{Int64} # Policy array, maps state index to action index
    action_map::Vector # Maps the action index to the concrete action type
    include_Q::Bool # Flag for including the Q-matrix
    mdp::Union{MDP,POMDP} # uses the model for indexing in the action function
end

# constructor with an optinal initial value function argument
function ValueIterationPolicy(mdp::Union{MDP,POMDP};
                              utility::AbstractVector{Float64}=zeros(n_states(mdp)),
                              policy::AbstractVector{Int64}=zeros(Int64, n_states(mdp)),
                              include_Q::Bool=true)
    ns = n_states(mdp)
    na = n_actions(mdp)
    @assert length(utility) == ns "Input utility dimension mismatch"
    @assert length(policy) == ns "Input policy dimension mismatch"
    action_map = ordered_actions(mdp)
    include_Q ? qmat = zeros(ns,na) : qmat = zeros(0,0)
    return ValueIterationPolicy(qmat, utility, policy, action_map, include_Q, mdp)
end

# constructor for solved q, util and policy
function ValueIterationPolicy(mdp::Union{MDP,POMDP}, q::AbstractMatrix{F}, util::AbstractVector{F}, policy::Vector{Int64}) where {F}
    action_map = ordered_actions(mdp)
    include_Q = true
    return ValueIterationPolicy(q, util, policy, action_map, include_Q, mdp)
end

# constructor for default Q-matrix
function ValueIterationPolicy(mdp::Union{MDP,POMDP}, q::AbstractMatrix{F}) where {F}
    (ns, na) = size(q)
    p = zeros(Int64, ns)
    u = zeros(ns)
    for i = 1:ns
        p[i] = argmax(q[i,:])
        u[i] = maximum(q[i,:])
    end
    action_map = ordered_actions(mdp)
    include_Q = true
    return ValueIterationPolicy(q, u, p, action_map, include_Q, mdp)
end

# returns the fields of the policy type
function locals(p::ValueIterationPolicy)
    return (p.qmat,p.util,p.policy,p.action_map)
end

function action(policy::ValueIterationPolicy, s::S) where S
    sidx = stateindex(policy.mdp, s)
    aidx = policy.policy[sidx]
    return policy.action_map[aidx]
end

function value(policy::ValueIterationPolicy, s::S) where S
    sidx = stateindex(policy.mdp, s)
    policy.util[sidx]
end

function POMDPPolicies.actionvalues(policy::ValueIterationPolicy, s::S) where S
    if !policy.include_Q
        error("ValueIterationPolicyError: the policy does not contain the Q function!")
    else
        sidx = stateindex(policy.mdp, s)
        return policy.qmat[sidx,:]
    end
end


text:
const InDims = Union{AllMetadata,Type,UnionAll,Dimension,IndexMode,ModeComponent,Symbol,Nothing}
const DimArrayOrStack = Union{AbstractDimArray,AbstractDimStack}

"""
    set(x, val)

Set the field matching the supertypes of values in xs and return a new object.

As DimensionalData is so strongly typed you do not need to specify what field
to `set` - there is no ambiguity.

To set fields of dimensions you need to specify the dimension. This can be done using
`Dimension => x` pairs, `X = x` keyword arguments, `Dimension` wrapped arguments,
or a `NamedTuple`.

When dimensions or IndexModes are passed to `set` to replace the existing ones,
fields that are not set will keep their original values.

## Notes:

Changing the dimension index range will set the `Sampled` mode
component `Regular` with a new step size, and set the dimension order.

Setting [`Order`](@ref) will *not* reverse the array or dimension to match. Use
[`reverse`](@ref) and [`reorder`](@ref) to do this.

## Examples

```julia
da = DimArray(rand(3, 4), (Dim{:custom}(10.0:010.0:30.0), Z(-20:010.0:10.0)))

# Set the array values
set(da, zeros(3, 4))

# Set the array name
set(da, "newname") # Swap dimension type # Using Pairs set(da, :Z => Ti, :custom => Z) set(da, :custom => X, Z => Y)
# Using keyword arguments
set(da, custom = X, Z = :a)

# Using Dimension wrappers
set(da, Dim{:custom}(X), Z(Dim{:a}))

# Set the dimension index

# To an `AbstractArray` set(da, Z => [:a, :b, :c, :d], :custom => Val((4, 5, 6)))

# To a `Val` tuple index (compile time indexing)
set(da, Z(Val((:a, :b, :c, :d))), custom = 4:6)

# Set dim modes
set(da, Z=NoIndex(), custom=Sampled())
set(da, :custom => Irregular(10, 12), Z => Regular(9.9))
set(da, (Z=NoIndex(), custom=Sampled()))
set(da, custom=Ordered(array=Reverse()), Z=Unordered())
```
"""
function set end
set(A::DimArrayOrStack, name::T) where {T<:Union{Mode,ModeComponent}} = _onlydimerror(T)
set(x::DimArrayOrStack, ::Type{T}) where T = set(x, T())
set(A, x) = _cantseterror(A, x)

"""
    set(x, args::Pairs...) => x with updated field/s
    set(x, args...; kw...) => x with updated field/s
    set(x, args::Tuple{Vararg{<:Dimension}}; kw...) => x with updated field/s

Set the dimensions or any properties of the dimensions for `AbstractDimArray`
or `AbstractDimStack`.

Set can be passed Keyword arguments or arguments of Pairs using dimension names,
tuples of values wrapped in the intended dimensions. Or fully or partially
constructed dimensions with val, mode or metadata fields set to intended
values. Dimension fields not assigned a value will be ignored, and the orginals kept.
"""
set(A::DimArrayOrStack, args::Union{Dimension,DimTuple,Pair}...; kw...) =
    rebuild(A, data(A), _set(dims(A), args...; kw...))
"""
    set(A::AbstractDimArray, data::AbstractArray) => AbstractDimArray

`AbstractArray` is always data, and update the `data` field of the array.
This is what is returned by `parent(A)`. It must be the same size as the
original value to match the `Dimension`s in the dims field.
"""
set(A::AbstractDimArray, newdata::AbstractArray) = begin
    axes(A) == axes(newdata) || _axiserr(A, newdata)
    rebuild(A; data=newdata)
end
"""
    set(A::AbstractDimArray, metadata::Union{AbstractMetadata,AbstractDict}) => AbstractDimArray

Update the `metadata` field of the array.
"""
set(A::AbstractDimArray, metadata::AllMetadata) = rebuild(A; metadata=metadata)
"""
    set(A::AbstractDimArray, name::AbstractName) => AbstractDimArray

Symbols are always names, and update the `name` field of the array.
"""
set(A::AbstractDimArray, name::Union{Symbol,AbstractName}) = rebuild(A; name=name)
"""
    set(s::AbstractDimStack, data::NamedTuple) => AbstractDimStack

`NamedTuple`s are always data, and update the `data` field of the dataset.
The values must be `AbstractArray of the same size as the original data, to
match the `Dimension`s in the dims field.
"""
set(s::AbstractDimStack, newdata::NamedTuple) = begin
    dat = data(s)
    keys(dat) === keys(newdata) || _keyerr(keys(dat), keys(newdata))
    map(dat, newdata) do d, nd
        axes(d) == axes(nd) || _axiserr(d, nd)
    end
    rebuild(s; data=newdata)
end
"""
    set(s::AbstractDimStack, metadata::AbstractMetadata) => AbstractDimStack

Update the `metadata` field of the stack.
"""
set(s::AbstractDimStack, metadata::AbstractMetadata) = rebuild(s; metadata=metadata)
"""
    set(dim::Dimension, index::Unioon{AbstractArray,Val}) => Dimension
    set(dim::Dimension, mode::Mode) => Dimension
    set(dim::Dimension, modecomponent::ModeComponent) => Dimension
    set(dim::Dimension, metadata::AbstractMetadata) => Dimension

Set fields of the dimension
"""
set(dim::Dimension, x::InDims) = _set(dim, x)

# Array or Stack
_set(A, x) = _cantseterror(A, x)

# Dimension
# Convert args/kw to dims and set
_set(dims_::DimTuple, args::Dimension...; kw...) = _set(dims_, (args..., _kwdims(kw)...))
# Convert pairs to wrapped dims and set
_set(dims_::DimTuple, p::Pair, ps::Vararg{<:Pair}) = _set(dims_, (p, ps...))
_set(dims_::DimTuple, ps::Tuple{Vararg{<:Pair}}) = _set(dims_, _pairdims(ps...))
# Set dims with (possibly unsorted) wrapper vals
_set(dims::DimTuple, wrappers::DimTuple) = begin
    # Check the dimension types match
    map(wrappers) do w
        hasdim(dims, w) || _wrongdimserr(dims, w)
    end
    # Missing dims return `nothing` from sortdims
    newdims = map(_set, dims, sortdims(wrappers, dims))
    # Swaps existing dims with non-nothing new dims
    swapdims(dims, newdims)
end
# Set things wrapped in dims
_set(dim::Dimension, wrapper::Dimension{<:InDims}) = _set(dim::Dimension, val(wrapper))
# Set the index
_set(dim::Dimension, index::Val) = rebuild(dim; val=index)
_set(dim::Dimension, index::Colon) = dim
_set(dim::Dimension, index::AbstractArray) = rebuild(dim; val=index)
_set(dim::Dimension, index::AbstractRange) = begin
    dim = _set(dim, _orderof(index))
    # We might need to update the index mode
    _set(mode(dim), dim, index)
end
# Update the Sampling mode of Sampled dims - it must match the range.
_set(mode::AbstractSampled, dim::Dimension, index::AbstractRange) =
    rebuild(dim; val=index, mode=_set(mode, Regular(step(index))))
_set(mode::IndexMode, dim::Dimension, index::AbstractRange) = rebuild(dim; val=index)
# Set the dim, checking the mode
_set(dim::Dimension, newdim::Dimension) = begin
    # Get new metadata and val
    dim = _set(dim, val(newdim))
    dim = _set(dim, metadata(newdim))
    dim = _set(dim, mode(newdim))
    # then wrap the updated dim in the new type
    basetypeof(newdim)(val(dim), mode(dim), metadata(dim))
end
# Construct types
_set(dim::Dimension, ::Type{T}) where T = _set(dim, T())
_set(dim::Dimension, key::Symbol) = _set(dim, key2dim(key))
_set(dim::Dimension, dt::DimType) = basetypeof(dt)(val(dim), mode(dim), metadata(dim))
# Set the mode
_set(dim::Dimension, newmode::IndexMode) = rebuild(dim; mode=_set(mode(dim), newmode))
# Otherwise pass this on to set fields on the mode
_set(dim::Dimension, x::ModeComponent) = _set(mode(dim), dim, x)
# Fallback dispatch without the dimension
_set(mode::IndexMode, dim::Dimension, x::ModeComponent) = rebuild(dim; mode=_set(mode, x))

# IndexMode
_set(mode::IndexMode, newmode::AutoMode) = mode
_set(mode::IndexMode, newmode::Categorical) =
    rebuild(newmode; order=_set(order(mode), order(newmode)))
_set(mode::IndexMode, newmode::AbstractSampled) = begin
    # Update each field separately. The old mode may not have these fields, or may have
    # a subset with the rest being traits. The new mode may have some auto fields.
    o = _set(order(mode), order(newmode))
    sp = _set(span(mode), span(newmode))
    sa = _set(sampling(mode), sampling(newmode))
    # Rebuild the new mode with the merged fields
    rebuild(newmode; order=o, span=sp, sampling=sa)
end
_set(mode::IndexMode, newmode::NoIndex) = newmode

# Order
_set(mode::IndexMode, neworder::Order) = rebuild(mode; order=_set(order(mode), neworder))
_set(mode::NoIndex, neworder::Order) = mode
_set(order::Union{Ordered,Unordered}, neworder::Ordered) = begin
    index = _set(indexorder(order), indexorder(neworder))
    array = _set(arrayorder(order), arrayorder(neworder))
    rel = _set(relation(order), relation(neworder))
    Ordered(index, array, rel)
end
_set(order::Union{Ordered,Unordered}, neworder::Unordered) =
    Unordered(_set(relation(order), relation(neworder)))
# AutoOrder
_set(order::Order, neworder::AutoOrder) = order
_set(order::AutoOrder, neworder::Order) = neworder
# Just for resolving ambiguity
_set(order::AutoOrder, neworder::AutoOrder) = AutoOrder()
# SubOrder
_set(order::Unordered, suborder::Relation) =
    rebuild(order; relation=_set(relation(order), suborder))
_set(order::Ordered, suborder::IndexOrder) =
    rebuild(order; index=_set(indexorder(order), suborder))
_set(order::Ordered, suborder::ArrayOrder) =
    rebuild(order; array=_set(arrayorder(order), suborder))
_set(order::Ordered, suborder::Relation) =
    rebuild(order; relation=_set(relation(order), suborder))
_set(suborder::ArrayOrder, newsuborder::ArrayOrder) = newsuborder
_set(suborder::IndexOrder, newsuborder::IndexOrder) = newsuborder
_set(suborder::Relation, newsuborder::Relation) = newsuborder

# Span
_set(mode::AbstractSampled, dim::Dimension, ::Irregular{AutoBounds}) =
    _set(dim, Irregular(bounds(dim)))
_set(mode::AbstractSampled, dim::Dimension, ::Explicit{AutoBounds}) = begin
    bnds = dim2boundsmatrix(locus(mode), span(mode), dim)
    rebuild(dim; mode=rebuild(mode; span=Explicit(bnds), sampling=Intervals(locus(mode))))
end
_set(span::Regular, dim::Dimension, ::Regular{AutoStep}) = dim
_set(span::Irregular, dim::Dimension, ::Irregular{AutoBounds}) = dim
_set(span::Explicit, dim::Dimension, ::Explicit{AutoBounds}) = dim

_set(mode::AbstractSampled, span::Span) = rebuild(mode; span=span)
_set(mode::AbstractSampled, span::AutoSpan) = mode
_set(span::Span, newspan::Span) = newspan
_set(span::Span, newspan::AutoSpan) = span

# Sampling
_set(mode::AbstractSampled, newsampling::Sampling) =
    rebuild(mode; sampling=_set(sampling(mode), newsampling))
_set(mode::AbstractSampled, sampling::AutoSampling) = mode
_set(sampling::Sampling, newsampling::Sampling) = newsampling
_set(sampling::Sampling, newsampling::AutoSampling) = sampling

# Locus
_set(mode::AbstractSampled, locus::Locus) =
    rebuild(mode; sampling=_set(sampling(mode), locus))
_set(sampling::Points, locus::Union{AutoLocus,Center}) = Points()
_set(sampling::Points, locus::Locus) = _locuserror()
_set(sampling::Intervals, locus::Locus) = Intervals(locus)
_set(sampling::Intervals, locus::AutoLocus) = sampling

# Metadata
_set(dim::Dimension, newmetadata::AllMetadata) =
    rebuild(dim, val(dim), mode(dim), newmetadata)

_set(x, ::Nothing) = x
_set(::Nothing, x) = x
_set(::Nothing, ::Nothing) = nothing

@noinline _locuserror() = throw(ArgumentError("Can't set a locus for `Points` sampling other than `Center` - the index values are the exact points"))
@noinline _cantseterror(a, b) = throw(ArgumentError("Can not set any fields of $(typeof(a)) to $(typeof(b))"))
@noinline _onlydimerror(T) = throw(ArgumentError("Can only set $(typeof(T)) for a dimension. Specify which dimension you want to set it for"))
@noinline _axiserr(a, b) = throw(ArgumentError("passed in axes $(axes(b)) do not match the currect axes $(axes(a))"))
@noinline _wrongdimserr(dims, w) = throw(ArgumentError("dim $(basetypeof(w))) not in $(map(basetypeof, dims))"))
@noinline _keyerr(ka, kb) = throw(ArgumentError("keys $ka and $kb do not match"))