text:
export Mesh

using Combinatorics
using Compat.Iterators

export vertexarray, cellarray

"""
    `U::Int`: indicating dimension of the embedding space
    `D1::Int`: one plus the manifold dimension
    `T<:AbstractFloat`: the type of a vertex coordinate
"""
abstract type AbstractMesh{U,D1,T} end

mutable struct Mesh{U,D1,T} <: AbstractMesh{U,D1,T}
    vertices::Vector{SVector{U,T}}
    faces::Vector{SVector{D1,Int}}
    """
    maps a face on its index in enumeration
    """
    dict::Dict{SVector{D1,Int},Int}
end

function Mesh(vertices, faces)
    dict = Dict((f,i) for (i,f) in enumerate(faces))
    T = eltype(eltype(vertices))
    U = length(eltype(vertices))
    D1 = length(eltype(faces))
    Mesh{U,D1,T}(vertices, faces, dict)
end

vertexarray(m::Mesh) = [ v[i] for v in m.vertices, i in 1:universedimension(m) ]
cellarray(m::Mesh) = [ k[i] for k in m.faces, i in 1:dimension(m)+1 ]

"""
    mesh(type, mdim, udim=mdim+1)

Returns an empty mesh with `coordtype` equal to `type`, of dimension `mdim`
and embedded in a universe of dimension `udim`
"""
mesh(T, mdim, udim=mdim+1) = Mesh(Pt{udim,T}[], SVector{mdim+1,Int}[])

"""
    vt = vertextype(mesh)

Returns type of the vertices used to define the cells of the mesh.
"""
vertextype(m::Mesh) = eltype(m.vertices)

"""
    celltype(mesh)

Returns the type of the index tuples stored in the mesh.
"""
celltype(m::Mesh) = eltype(m.faces)

"""
    coordtype(mesh)

Returns `eltype(vertextype(mesh))`
"""
coordtype(m::AbstractMesh{U,D1,T}) where {U,D1,T} = T
# coordtype(m::Mesh) = eltype(vertextype(m))

"""
    dim = dimension(mesh)

Returns the dimension of the mesh. Note that this is
the dimension of the cells, not of the surrounding space.
"""
dimension(m::AbstractMesh{U,D1}) where {U,D1} = D1-1
# dimension(m::Mesh) = length(celltype(m)) - 1

"""
    universedimension(mesh)

Returns the dimension of the surrounding space. Equals
the number of coordinates required to describe a vertex.
"""
universedimension(m::Mesh) = length(vertextype(m))

"""
    vertices(mesh)

Returns an indexable iterable to the vertices of the mesh
"""
vertices(m::Mesh) = m.vertices

#     vertices(mesh, i)
#     vertices(mesh, I)
#
# Select one or multiple vertices from the mesh. In the second
# form, only statically sized arrays are allowed to discourage
# memory allocation. The returned vector in that case will also
# be statically typed and of the same size as `I`.
vertices(m::Mesh, i::Number) = m.vertices[i]
@generated function vertices(m::Mesh, I::SVector)
    N = length(I)
    xp = :(())
    for i in 1:N
        push!(xp.args, :(m.vertices[I[$i]]))
    end
    :(SVector($xp))
end

"""
    numvertices(mesh)

Returns the number of vertices in the mesh.

*Note*: this is the number of vertices in the vertex buffer and might include floatin vertices
or vertices not appearing in any cell. In other words the following is not necessarily true:

```julia
    numvertices(mesh) == numcells(skeleton(mesh,0))
```
"""
numvertices(m::Mesh) = length(m.vertices)

"""
    numcells(mesh)

Returns the number of cells in the mesh.
"""
numcells(m::Mesh) = length(m.faces)

"""
    cells(mesh)

Return an iterable collection containing the cells making up the mesh.
"""
cells(mesh) = mesh.faces

# """
#     cellvertices(mesh, i)
#
# Return an indexable collection containing the vertices of cell `i`.
# Shorthand for `vertices(mesh, cells(mesh, i))`.
# """
# cellvertices(mesh,i) = vertices(mesh, cells(mesh, i))

"""
    translate(mesh, v)

Creates a new mesh by translating `mesh` over vector `v`
"""
translate(Γ::Mesh, v) = Mesh([w + v for w in Γ.vertices], deepcopy(Γ.faces))

"""
    translate!(mesh, v)

Translates `mesh` over vector `v` inplace.
"""
function translate!(Γ::Mesh, v)
    for i in 1:length(Γ.vertices)
        Γ.vertices[i] += v
    end
    Γ
end

"""
    flip(cell)

Change the orientation of a cell by interchanging the first to indices.
"""
@generated function flip(cell)
    # generate `T(cell[2],cell[1],cell[3],...)`, with `T = typeof(cell)`
    N = length(cell)
    xp = :($cell(cell[2], cell[1]))
    for i in 3:N
        push!(xp.args, :(cell[$i]))
    end
    xp
end

"""
    flipmesh!(mesh)

Change the orientation of a mesh
"""
function flipmesh!(mesh)
    mesh.faces .= flip.(mesh.faces)
    mesh
end

flipmesh(mesh) = flipmesh!(deepcopy(mesh))

export mirrormesh, mirrormesh!

"""
    mirror(vertex, normal, anchor)

Mirror vertex across a plane defined by its normal and a containing point.
"""
function mirror(vertex, normal, anchor)
    h = dot(vertex - anchor, normal)
    vertex - 2*h * normal
end

function mirrormesh!(mesh, normal, anchor)
    #mesh.vertices .= mirror.(mesh.vertices, normal, anchor)
    for i in eachindex(mesh.vertices)
        mesh.vertices[i] = mirror(mesh.vertices[i], normal, anchor)
    end
    mesh
end

mirrormesh(mesh, args...) = mirrormesh!(deepcopy(mesh), args...)

function rotate!(Γ::Mesh{3}, v)
    α = norm(v)
    u = v/α

    cα = cos(α/2)
    sα = sin(α/2)

    p = point(cα, sα*u[1], sα*u[2], sα*u[3])
    q = point(cα, -sα*u[1], -sα*u[2], -sα*u[3])

    P = @SMatrix [
        +p[1] -p[2] -p[3] -p[4];
        +p[2] +p[1] -p[4] +p[3];
        +p[3] +p[4] +p[1] -p[2];
        +p[4] -p[3] +p[2] +p[1]]

    Q = @SMatrix [
        +q[1] -q[2] -q[3] -q[4];
        +q[2] +q[1] +q[4] -q[3];
        +q[3] -q[4] +q[1] +q[2];
        +q[4] +q[3] -q[2] +q[1]]

    for i in 1:numvertices(Γ)
        r = vertices(Γ,i)
        t = point(0, r[1], r[2], r[3])
        s = P * Q * t
        Γ.vertices[i] = point(s[2], s[3], s[4])
    end
    Γ
end

function rotate(Γ::Mesh{3}, v)
    R = deepcopy(Γ)
    rotate!(R,v)
    R
end

"""
    fliporientation(mesh)

Changes the mesh orientation inplace. If non-orientatble, undefined.
"""
function fliporientation!(m::Mesh)
    for i in 1:numcells(m)
        m.faces[i] = fliporientation(m.faces[i])
    end
    return m
end

"""
    fliporientation(mesh)

Returns a mesh of opposite orientation.
"""
function fliporientation(m::Mesh)
    n = deepcopy(m)
    fliporientation!(n)
end

@generated function fliporientation(I::SVector{N,T}) where {N,T}
    @assert N >= 2
    xp = :(SVector{N,T}(I[2],I[1]))
    for i in 3:N
        push!(xp.args, :(I[$i]))
    end
    return xp
end

"""
    -mesh -> flipped_mesh

Create a mesh with opposite orientation.
"""
Base.:-(m::Mesh) = fliporientation(m)

"""
    boundary(mesh)

Returns the boundary of `mesh` as a mesh of lower dimension.
"""
function boundary(mesh)

    D = dimension(mesh)

    # vertices have no boundary
    @assert 0 < D

    # build a list of D-1 cells
    edges = skeleton(mesh, D-1)
    faces = skeleton(mesh, D)

    # get the edge-face connection matrix
    conn = connectivity(edges, faces)

    # find the edges that only have one adjacent face
    #i = find(x -> x < 2, sum(abs.(conn), dims=1))
    sums = sum(abs.(conn), dims=1)
    i = (LinearIndices(sums))[findall(x->x<2, sums)]

    # create a mesh out of these
    bnd = Mesh(mesh.vertices, edges.faces[i])
end

"""
Complement to boundary. This function selects those edges that have at
least two faces adjacent. The case with more than two neighboring faces
occurs on non-manifold structures (e.g. containing junctions)
"""
function interior(mesh::Mesh, edges=skeleton(mesh,1))
    @assert dimension(mesh) == 2
    @assert vertices(mesh) === vertices(edges)

    C = connectivity(edges, mesh)
    @assert size(C) == (numcells(mesh), numcells(edges))

    nn = vec(sum(abs.(C), dims=1))
    T = CompScienceMeshes.celltype(edges)
    interior_edges = Vector{T}()
    for (i,edge) in pairs(cells(edges))
        nn[i] > 1 && push!(interior_edges, edge)
    end
    Mesh(vertices(mesh), interior_edges)
end

"""
    vertextocellmap(mesh) -> vertextocells, numneighbors

Computed an V×M array `vertextocells` where V is the number of vertices
and M is the maximum number of cells adjacent to any given vertex such
that `vertextocells[v,i]` is the index in the cells of `mesh` of the `i`th
cell adjacent to teh `v`-th vertex. `numneighbors[v]` contains the number
of cells adjacent to the `v`-th vertex.

This method allows e.g. for the efficient computation of the connectivity
matrix of the mesh.
"""
function vertextocellmap(mesh)

    # cells = mesh.faces
    C = cells(mesh)

    numverts = length(vertices(mesh)) #numvertices(mesh)
    numcells = length(C)
    numneighbors = zeros(Int, numverts)
    for i = 1 : numcells
        for m in C[i]
            numneighbors[m] += 1
        end
        #numneighbors[ cells[i] ] += 1
    end

    npos = -1
    vertstocells = fill(npos, numverts, maximum(numneighbors))
    numneighbors = zeros(Int, numverts)
    for i = 1 : numcells
        cell = C[i]
        for j = 1 : length(cell)
            v = cell[j]
            k = (numneighbors[v] += 1)
            vertstocells[v,k] = i
        end
    end

    vertstocells, numneighbors
end

"""
    skeleton(mesh, dim)

Returns a mesh comprising the `dim`-dimensional sub cells of `mesh`. For example to retrieve
the edges of a given surface `mesh`,

```julia
edges = skelton(mesh, 1)
```
"""
function skeleton(mesh, dim::Integer)

    meshdim = dimension(mesh)
    @assert 0 <= dim <= meshdim

    if dim == meshdim
        return mesh
    end

    nc = numcells(mesh)
    C = SVector{dim+1,Int}
    simplices = zeros(C, nc*binomial(meshdim+1,dim+1))

    n = 1
    for c = 1 : nc

        cell = mesh.faces[c]
        for simplex in combinations(cell,dim+1)
            simplices[n] = sort(simplex)
            n += 1
        end
    end

    simplices = unique(simplices)
    Mesh(mesh.vertices, simplices)
end

"""
    skeleton(pred, mesh, dim)

Like `skeleton(mesh, dim)`, but only cells for which `pred(cell)`
returns true are withheld.
"""
function skeleton(pred, mesh, dim)

    meshdim = dimension(mesh)
    @assert 0 <= dim <= meshdim

    nc = numcells(mesh)
    C = SVector{dim+1,Int}
    simplices = zeros(C, nc*binomial(meshdim+1,dim+1))

    n = 1
    for c = 1 : nc

        cell = mesh.faces[c]
        for simplex in combinations(cell,dim+1)
            if pred(SVector{dim+1,Int}(simplex))
                simplices[n] = sort(simplex)
                n += 1
            end
        end
    end

    simplices = simplices[1:n-1]
    simplices = unique(simplices)

    Mesh(mesh.vertices, simplices)
end

"""
    connectivity(faces, cells, op=sign)

Create a sparse matrix `D` of size `numcells(cells)` by `numcells(faces)` that
contiains the connectivity info of the mesh. In particular `D[m,k]` is `op(r)`
where `r` is the local index of face `k` in cell `m`. The sign of `r` is
positive or negative depending on the relative orientation of face `k` in cell
`m`.

For `op=sign`, the matrix returned is the classic connectivity matrix, i.e.
the graph version of the exterior derivative.
"""
function connectivity(kcells::Mesh, mcells::Mesh, op = sign)

    vtok, _ = vertextocellmap(kcells)
    vtom, _ = vertextocellmap(mcells)

    npos = -1

    dimk = numcells(kcells)
    dimm = numcells(mcells)

    D = spzeros(Int, dimm, dimk)

    for v in 1:numvertices(kcells)
        for i in vtok[v,:]
            i == npos && break
            #kcell = cells(kcells, i)
            kcell = kcells.faces[i]
            for j in vtom[v,:]
                j == npos && break
                #mcell = cells(mcells, j)
                mcell = mcells.faces[j]
                D[j,i] = op(relorientation(kcell, mcell))
            end
        end
    end

    return D
end

"""
    pairs = cellpairs(mesh, edges, dropjunctionpair=false)

Given a mesh and set of oriented edges from that mesh (as generated by `skeleton`),
    `cellpairs` will generate a 2 x K matrix, where K is the number of pairs
    and each column contains a pair of indices in the cell array of `mesh` that have
    one of the supplied edges in common.

Returns an array of pairs of indices, each pair corresponding to a pair of adjacent faces.

(If the mesh is oriented, the first row of `facepairs` will contain indices to the cell
    for which the corresponding edge has a positive relative orientation.

If a edge lies on the boundary of the mesh, and only has one neighboring cell, the
    second row of `facepairs` will contain `-k` with `k` the local index of the corresponding
    edge in its neighboring triangle.

If an edge has more than two neighboring cells (i.e. the edge is on a junction),
    all possible pairs of cells that have the junction edge in common are supplied. if
    `dropjunctionpair == false` then one of the possible pairs of cells is not recorded.
    This is done to avoid the creation of linearly dependent basis functions in the
    construction of boundary element methods for Maxwell's equations.)
"""
function cellpairs(mesh, edges; dropjunctionpair=false)

    ndrops = dropjunctionpair ? 1 : 0
    @assert dimension(edges)+1 == dimension(mesh)

    numedges = numcells(edges)

    # perform a dry run to determine the number of cellpairs
    v2e, nn = vertextocellmap(mesh)
    k = 0
    for edge = edges.faces

        v = edge[1]; n = nn[v]; nbd1 = v2e[v,1:n]
        v = edge[2]; n = nn[v]; nbd2 = v2e[v,1:n]

        nbd = intersect(nbd1, nbd2)
        n = length(nbd)
        if n < 3
            k += 1
        else
            k += (n-ndrops)
        end

    end
    facepairs = zeros(Int, 2, k)

    k = 1
    for edge in edges.faces

        # neighborhood of startvertex
        v = edge[1]
        n = nn[v]
        nbd1 = v2e[v,1:n]

        # neighborhood of endvertex
        v = edge[2]
        n = nn[v]
        nbd2 = v2e[v,1:n]

        # get the intersection
        nbd = intersect(nbd1, nbd2)
        n = length(nbd)
        @assert 0 < n

        if n == 1 # boundary edge

            c = nbd[1]
            cell = mesh.faces[c]
            s = relorientation(edge, cell)
            facepairs[1,k] = c
            facepairs[2,k] = -abs(s)

            k += 1
        elseif n == 2

            c = nbd

            cell1 = mesh.faces[c[1]]
            cell2 = mesh.faces[c[2]]

            r1 = relorientation(edge, cell1)
            r2 = relorientation(edge, cell2)

            if r1 > 0 && r2 < 0
                c1 = c[2]
                c2 = c[1]
            else
                c1 = c[1]
                c2 = c[2]
            end

            facepairs[1,k] = c1
            facepairs[2,k] = c2

            k += 1
        else

            for c in drop(combinations(nbd,2), ndrops)

                cell1 = mesh.faces[c[1]]
                cell2 = mesh.faces[c[2]]

                r1 = relorientation(edge, cell1)
                r2 = relorientation(edge, cell2)

                if r1 > 0 && r2 < 0
                    c1 = c[2]
                    c2 = c[1]
                else
                    c1 = c[1]
                    c2 = c[2]
                end

                facepairs[1,k] = c1
                facepairs[2,k] = c2

                k += 1
            end
        end
    end

    facepairs
end

"""
    chart(mesh, cell) -> cell_chart

Return a chart describing the supplied cell of `mesh`.
"""
chart(mesh::Mesh, cell) = simplex(vertices(mesh,cell))

"""
    isoriented(mesh) -> Bool

Returns true is all cells are consistently oriented, false otherwise.
"""
function isoriented(m::AbstractMesh)

    @assert dimension(m) >= 0
    edges = skeleton(m, dimension(m)-1)

    D = connectivity(edges, m)
    S = (abs.(sum(D,dims=1)) .<= 1)
    return all(S)
end

max_stars_repo_name: shakisparki/CompScienceMeshes

max_issues_repo_name: shakisparki/CompScienceMeshes

max_forks_repo_name: shakisparki/CompScienceMeshes


text:
export RandomIdGenerator

using Random

abstract type AbstractIdGenerator end

Base.@kwdef struct RandomIdGenerator{R} <: AbstractIdGenerator
    rng::R = Random.GLOBAL_RNG
end

Base.rand(g::RandomIdGenerator, T) = rand(g.rng, T)

generate_trace_id(rng::RandomIdGenerator) = rand(rng.rng, API.TraceIdType)
generate_span_id(rng::RandomIdGenerator) = rand(rng.rng, API.SpanIdType)

ext: jl


text:
@doc raw"""
    simple_estimator(model::AshkinTeller, T::Real, Jsigma, Jtau, K)

Returns the following observables as `Dict{String, Any}`

# Observables
- `"Energy"`
    - energy density
- `"Energy^2"`
    - square of energy density
- `"|Magnetization|"`
    - absolute value of magnetization density, ``|m| = \sqrt{ (\sum_i \sigma_i )^2 + (\sum_i \tau_i)^2 } / N``
- `"|Magnetization|^2"`
    - square of magnetization density
- `"|Magnetization|^4"`
    - quadruple of magnetization density
- `"Magnetization sigma"`
    - magnetization density (sigma spin)
- `"|Magnetization sigma|"`
- `"Magnetization sigma^2"`
- `"Magnetization sigma^4"`
- `"Magnetization tau"`
    - magnetization density (tau spin)
- `"|Magnetization tau|"`
- `"Magnetization tau^2"`
- `"Magnetization tau^4"`

"""
function simple_estimator(model::AshkinTeller, T::Real, Jsigma, Jtau, K, _=nothing)
    nsites = numsites(model)
    nbonds = numbonds(model)

    Ms = mean(model.spins, dims=2)
    E = 0.0
    @inbounds for b in bonds(model)
        s1, s2 = source(b), target(b)
        bt = bondtype(b)
        E += model.spins[1,s1] * model.spins[1,s2] * Jsigma[bt]
        E += model.spins[2,s1] * model.spins[2,s2] * Jtau[bt]
        E += model.spins[1,s1] * model.spins[2,s1] * model.spins[1,s2] * model.spins[2,s2] * K[bt]
    end
    E /= nsites

    res = Measurement()
    res["|Magnetization|^2"] = Ms[1]^2 + Ms[2]^2
    res["|Magnetization|^4"] = res["|Magnetization|^2"]^2
    res["|Magnetization|"] = sqrt(res["|Magnetization|^2"])

    res["Magnetization sigma"] = Ms[1]
    res["|Magnetization sigma|"] = abs(Ms[1])
    res["Magnetization sigma^2"] = Ms[1]^2
    res["Magnetization sigma^4"] = Ms[1]^4

    res["Magnetization tau"] = Ms[2]
    res["|Magnetization tau|"] = abs(Ms[2])
    res["Magnetization tau^2"] = Ms[2]^2
    res["Magnetization tau^4"] = Ms[2]^4

    res["Energy"] = E
    res["Energy^2"] = E^2
    return res
end

default_estimator(model::AshkinTeller, update) = simple_estimator


text:
using NewtonsMethod
using Test

# Define functions to be tested
m(x) = (x - 4)^3
m′(x) = 3*(x - 4)^2
h(x) = log(x)
h′(x) = 1/x
b(x) = 2*x + x^4
b′(x) = 2 + 4*x^3
g(x) = 2 + x^2
g′(x) = 2*x

@testset "NewtonsMethod.jl" begin
    # Run tests on newton root with analytical derivatives
    @test newtonroot(m, m′, x_0=0.5).value ≈ 4.0 rtol=1E-7
    @test newtonroot(h, h′, x_0=0.5).value ≈ 1.0 rtol=1E-7
    @test newtonroot(b, b′, x_0=0.5).value ≈ 0.0 rtol=1E-7
    @test newtonroot(b, b′, x_0=-2.5).value ≈ -2^(1/3) rtol=1E-7

    # Run tests on newton root with nummerical derivatives
    @test newtonroot(m, x_0=0.5).value ≈ 4.0 rtol=1E-7
    @test newtonroot(h, x_0=0.5).value ≈ 1.0 rtol=1E-7
    @test newtonroot(b, x_0=0.5).value ≈ 0.0 rtol=1E-7
    @test newtonroot(b, x_0=-2.5).value ≈ -2^(1/3) rtol=1E-7

    # Run tests on newton root with analytical derivatives, BigFloat
    @test newtonroot(m, m′, x_0=BigFloat(0.5)).value ≈ 4.0 rtol=1E-7
    @test newtonroot(h, h′, x_0=BigFloat(0.5)).value ≈ 1.0 rtol=1E-7
    @test newtonroot(b, b′, x_0=BigFloat(-2.5)).value ≈ -2^(1/3) rtol=1E-7

    # Run tests on newton root with nummerical derivatives, BigFloat
    @test newtonroot(m, x_0=BigFloat(0.5)).value ≈ 4.0 rtol=1E-7
    @test newtonroot(h, x_0=BigFloat(0.5)).value ≈ 1.0 rtol=1E-7
    @test newtonroot(b, x_0=BigFloat(-2.5)).value ≈ -2^(1/3) rtol=1E-7

    # Run tests with non-convergence function
    @test newtonroot(g,g′, x_0=0.5) == nothing
    @test newtonroot(g, x_0=0.5) == nothing

    # Run tests with changed maxiter and tol for some of the functions
    @test newtonroot(h, x_0=0.5, maxiter=5) == nothing
    @test newtonroot(b, b′, x_0=BigFloat(-2.5), maxiter=10_000).value ≈ -2^(1/3)
    @test newtonroot(m, x_0=0.5, tol=10E-9).value ≈ 4.0 rtol=1E-7
    @test !isapprox(newtonroot(m, m′, x_0=BigFloat(0.5), tol=10E-2).value,4.0)

end


text:
struct PORollout
    solver::Union{POMDPs.Solver,POMDPs.Policy,Function}
    updater::POMDPs.Updater
end

struct SolvedPORollout{P<:POMDPs.Policy,U<:POMDPs.Updater,RNG<:AbstractRNG}
    policy::P
    updater::U
    rng::RNG
end

struct FORollout # fully observable rollout
    solver::Union{POMDPs.Solver,POMDPs.Policy}
end

struct SolvedFORollout{P<:POMDPs.Policy,RNG<:AbstractRNG}
    policy::P
    rng::RNG
end

struct FOValue
    solver::Union{POMDPs.Solver, POMDPs.Policy}
end

struct SolvedFOValue{P<:POMDPs.Policy}
    policy::P
end

"""
    estimate_value(estimator, problem::POMDPs.POMDP, start_state, h::BeliefNode, steps::Int)

Return an initial unbiased estimate of the value at belief node h.

By default this runs a rollout simulation
"""
function estimate_value end
estimate_value(f::Function, pomdp::POMDPs.POMDP, start_state, h::BeliefNode, steps::Int) = f(pomdp, start_state, h, steps)
estimate_value(n::Number, pomdp::POMDPs.POMDP, start_state, h::BeliefNode, steps::Int) = convert(Float64, n)

function estimate_value(estimator::Union{SolvedPORollout,SolvedFORollout}, pomdp::POMDPs.POMDP, start_state, h::BeliefNode, steps::Int)
    rollout(estimator, pomdp, start_state, h, steps)
end

@POMDP_require estimate_value(estimator::Union{SolvedPORollout,SolvedFORollout}, pomdp::POMDPs.POMDP, start_state, h::BeliefNode, steps::Int) begin
    @subreq rollout(estimator, pomdp, start_state, h, steps)
end

function estimate_value(estimator::SolvedFOValue, pomdp::POMDPs.POMDP, start_state, h::BeliefNode, steps::Int)
    POMDPs.value(estimator.policy, start_state)
end

function convert_estimator(ev::RolloutEstimator, solver, pomdp)
    policy = MCTS.convert_to_policy(ev.solver, pomdp)
    SolvedPORollout(policy, updater(policy), solver.rng)
end

function convert_estimator(ev::PORollout, solver, pomdp)
    policy = MCTS.convert_to_policy(ev.solver, pomdp)
    SolvedPORollout(policy, ev.updater, solver.rng)
end

function convert_estimator(est::FORollout, solver, pomdp)
    policy = MCTS.convert_to_policy(est.solver, pomdp)
    SolvedFORollout(policy, solver.rng)
end

function convert_estimator(est::FOValue, solver::AbstractPOMCPSolver, pomdp::POMDPs.POMDP)
    policy = MCTS.convert_to_policy(est.solver, pomdp)
    SolvedFOValue(policy)
end

"""
Perform a rollout simulation to estimate the value.
"""
function rollout(est::SolvedPORollout, pomdp::POMDPs.POMDP, start_state, h::BeliefNode, steps::Int)
    b = extract_belief(est.updater, h)
    sim = RolloutSimulator(est.rng,
                           steps)
    return POMDPs.simulate(sim, pomdp, est.policy, est.updater, b, start_state)
end

@POMDP_require rollout(est::SolvedPORollout, pomdp::POMDPs.POMDP, start_state, h::BeliefNode, steps::Int) begin
    @req extract_belief(::typeof(est.updater), ::typeof(h))
    b = extract_belief(est.updater, h)
    sim = RolloutSimulator(est.rng,
                           steps)
    @subreq POMDPs.simulate(sim, pomdp, est.policy, est.updater, b, start_state)
end

function rollout(est::SolvedFORollout, pomdp::POMDPs.POMDP, start_state, h::BeliefNode, steps::Int)
    sim = RolloutSimulator(est.rng,
                                        steps)
    return POMDPs.simulate(sim, pomdp, est.policy, start_state)
end

@POMDP_require rollout(est::SolvedFORollout, pomdp::POMDPs.POMDP, start_state, h::BeliefNode, steps::Int) begin
    sim = RolloutSimulator(est.rng,
                                        steps)
    @subreq POMDPs.simulate(sim, pomdp, est.policy, start_state)
end

"""
    extract_belief(rollout_updater::POMDPs.Updater, node::BeliefNode)

Return a belief compatible with the `rollout_updater` from the belief in `node`.

When a rollout simulation is started, this function is used to create the initial belief (compatible with `rollout_updater`) based on the appropriate `BeliefNode` at the edge of the tree. By overriding this, a belief can be constructed based on the entire tree or entire observation-action history.
"""
function extract_belief end

# some defaults are provided
extract_belief(::NothingUpdater, node::BeliefNode) = nothing

function extract_belief(::PreviousObservationUpdater, node::BeliefNode)
    if node.node==1 && !isdefined(node.tree.o_labels, node.node)
        missing
    else
        node.tree.o_labels[node.node]
    end
end


text:
my_f(x,y) = 2x+3y
derivative_of_my_f(x,y) = ForwardDiff.derivative(x->my_f(x,y),x)


text:
module test_julia_axes

using Test

A = [1 2; 3 4; 5 6; 7 8;]
#=
4×2 Array{Int64,2}:
 1  2
 3  4
 5  6
 7  8
=#

@test size(A) == (4, 2)
@test axes(A, 1) == Base.OneTo(4)
@test axes(A, 2) == Base.OneTo(2)

end # module test_julia_axes


text:
using Test, Astro

# references are to a book by Meeus

print("22.a Nutation (delta psi)")
deltaPsi = nut_in_lon(2446895.5)
d, m, s = d_to_dms(rad2deg(deltaPsi))
@test(d == 0)
@test(m == 0)
@test isapprox(s, -3.788, atol=0.001)
println(" passed")

print("22.a Nutation (delta epsilon)")
deltaEps = nut_in_obl(2446895.5)
d, m, s = d_to_dms(rad2deg(deltaEps))
@test(d == 0)
@test(m == 0)
@test isapprox(s, 9.443, atol=0.001)
println(" passed")

print("22.a Nutation (epsilon)")
epsi = obliquity(2446895.5)
d, m, s = d_to_dms(rad2deg(epsi))
@test(d == 23)
@test(m == 26)
@test isapprox(s, 27.407, atol=0.001)
println(" passed")

print("22.a Nutation (epsilon - high precision)")
epsi = obliquity_high(2446895.5)
d, m, s = d_to_dms(rad2deg(epsi))
@test(d == 23)
@test(m == 26)
@test isapprox(s, 27.407, atol=0.001)
println(" passed")


text:
# ps3000_get_timebase.jl

export ps3000_get_timebase

function ps3000_get_timebase(handle::Int16, timebase, no_of_samples, oversample)
  time_interval = Ref{Int32}()
  time_units = Ref{Int16}()
  max_samples = Ref{Int32}()
  ps_status = ccall((:ps3000_get_timebase, ps3000driver),Int16,
    (Int16, Int16, Int32, Ref{Int32}, Ref{Int16}, Int16, Ref{Int32}),
    handle, timebase, no_of_samples, time_interval,
    time_units, oversample, max_samples)
  if ps_status == 0
    error("ps3000 Error: Parameter Out Of Range")
  end
  return (time_interval[], time_units[], max_samples[])
end


text:
using Base.Threads
using Base.Test
using Base.Threads.Atomic
using LightGraphs
using GraphBenchmarks
using BenchmarkTools
using UnsafeAtomics

import Base: push!, shift!, isempty, getindex
import GraphBenchmarks: preparebench

include("lg.jl")

immutable ThreadQueue{T}
    data::Vector{T}
    head::Atomic{Int}
    tail::Atomic{Int}
end

abstract LGBFSAlgs <: AbstractBFS

immutable NaiveSerialBFS <: LGBFSAlgs end
immutable ParallBFS <: LGBFSAlgs end
immutable LevelSynchronous <: LGBFSAlgs end

function ThreadQueue(T::Type, maxlength::Int)
    q = ThreadQueue(Vector{T}(maxlength), Atomic{Int}(1), Atomic{Int}(1))
    return q
end

function push!{T}(q::ThreadQueue{T}, val::T)
    # TODO: check that head > tail
    offset = atomic_add!(q.tail, 1)
    q.data[offset] = val
    return offset
end

function shift!{T}(q::ThreadQueue{T})
    # TODO: check that head < tail
    offset = atomic_add!(q.head, 1)
    return q.data[offset]
end

function isempty(q::ThreadQueue)
    return ( q.head[] == q.tail[] ) && q.head != 1
    # return q.head == length(q.data)
end

function getindex{T}(q::ThreadQueue{T}, iter)
    return q.data[iter]
end

function bfs(alg::NaiveSerialBFS, next::Vector{Int}, g::SimpleGraph, source::Int64)
    parents = fill(-2, nv(g)) #Initialize parents array with -2's.
    parents[source]=-1 #Set source to -1
    while !isempty(next)
        src = shift!(next) #Get first element
        vertexneighbors = neighbors(g, src)
        for vertex in vertexneighbors
            #If not already set, and is not found in the queue.
            if parents[vertex]==-2
                push!(next, vertex) #Push onto queue
                parents[vertex] = src
            end
        end
    end
    return parents
end

function bfskernel(alg::ParallBFS, next::ThreadQueue, src::Int64, parents::Array{Int64}, vertexneighbors::Array{Int64})
    @threads for vertex in vertexneighbors
        #If not already set, and is not found in the queue.
        if parents[vertex]==-2
            # println("pushing: $vertex")
            push!(next, vertex) #Push onto queue
            # println("pushed: $vertex")
            parents[vertex] = src
        end
    end
end

function bfs(alg::ParallBFS, next::ThreadQueue, g::SimpleGraph, source::Int64, parents::Array{Int64})
    parents[source]=-1 #Set source to -1
    push!(next, source)
    while !isempty(next)
        src = shift!(next) #Get first element
        vertexneighbors = neighbors(g, src)
        bfskernel(alg, next, src, parents, vertexneighbors)
    end
    return parents
end

function bfskernel(alg::LevelSynchronous, next::ThreadQueue, g::SimpleGraph, parents::Array{Int64}, level::Array{Int64})
    @threads for src in level
        vertexneighbors = neighbors(g, src)
        for vertex in vertexneighbors
            #Set parent value if not set yet.
            parent = UnsafeAtomics.unsafe_atomic_cas!(parents, vertex, -2, src)
            if parent==-2
                push!(next, vertex) #Push onto queue
            end
        end
    end
end

function bfs(
        alg::LevelSynchronous, next::ThreadQueue, g::SimpleGraph, source::Int64,
        parents::Array{Int64}
    )
    parents[source]=-1 #Set source to -1
    push!(next, source)
    while !isempty(next)
        level = next[next.head[]:next.tail[]-1]
        next.head[] = next.tail[] #reset the queue
        bfskernel(alg, next, g, parents, level)
    end
    return parents
end

function bfs(alg::NaiveSerialBFS, g::SimpleGraph, source::Int, nv::Int)
    next = Vector{Int64}([source])
    sizehint!(next, nv)
    return bfs(NaiveSerialBFS(), next, g, source)
end

function bfs(alg::LevelSynchronous, g::SimpleGraph, source::Int64, nv::Int64)
    next = ThreadQueue(Int, nv)
    parents = fill(-2, nv)
    bfs(alg, next, g, source, parents)
end

function bfs(alg::ParallBFS, g::SimpleGraph, source::Int64, nv::Int64)
    next = ThreadQueue(Int, nv)
    parents = fill(-2, nv)
    bfs(alg, next, g, source, parents)
end

function preparebench(benchmark::GraphBenchmarkSpec, alg::LGBFSAlgs, lg::LGGraph)
    n = nv(lg.g)
    sources = picksources(benchmark, n)
    bfsbench = @benchmarkable begin
        for src in $sources
            bfs($alg, $(lg.g), src, $n)
        end
    end seconds=6000 samples=10
    bfsbench
end

function bfsbenchsuite{T <: SimpleGraph}(t::Type{T}, bench::GraphBenchmarkSpec, scale::Int, edgefactor::Int)
    srand(0)
    gen = Kronecker(scale, edgefactor)
    suite = BenchmarkGroup(["bfs", "lg"])
    suite["serial"] = BenchmarkGroup(["serialbfs"])
    suite["serial"]["lgbfs"] = preparebenchmark(bench, LGGraph{T}, BFS(), gen)
    suite["serial"]["naivebfs"] = preparebenchmark(bench, LGGraph{T}, NaiveSerialBFS(), gen)
    suite["parallel"] = BenchmarkGroup(["parallelbfs"])
    suite["parallel"]["parallbfs"] = preparebenchmark(bench, LGGraph{T}, ParallBFS(), gen)
    suite["parallel"]["levelsyncbfs"] = preparebenchmark(bench, LGGraph{T}, LevelSynchronous(), gen)
    suite
end

function comparebfs{T <: SimpleGraph}(t::Type{T}, bench::GraphBenchmarkSpec, scale::Int, edgefactor::Int)
    suite = bfsbenchsuite(T, bench, scale, edgefactor)
    run(suite, verbose=true)
end
#Taken from the LightGraphs BFS test
function istree(parents::Vector{Int}, maxdepth)
    flag = true
    for i in 1:maxdepth
        s = i
        depth = 0
        while parents[s] > 0 && parents[s] != s
            s = parents[s]
            depth += 1
            if depth > maxdepth
                return false
            end
        end
    end
    return flag
end

function testbfsalgs()
    gen = Kronecker(10, 16)
    lg = construct(LGGraph{DiGraph}, generate(gen))
    n = nv(lg.g)
    visitor = LightGraphs.TreeBFSVisitorVector(zeros(Int,n))
    LightGraphs.bfs_tree!(visitor, lg.g, 1)
    @test istree(visitor.tree, n)
    @test istree(bfs(NaiveSerialBFS(), lg.g, 1, n), n)
    @test istree(bfs(LevelSynchronous(), lg.g, 1, n), n)
    @test istree(bfs(ParallBFS(), lg.g, 1, n), n)
end


text:
__precompile__()

module Interpolation

    export Interpolator, HermiteInterpolation, evaluate

    include("interpolation/interpolation.jl")
    include("interpolation/hermite_interpolation.jl")

end


text:
@testset "traversal" begin
    tree = ParseNewick("((B:5,A:5)C:9,(D:5,E:5)F:5)G:5;")
    po_list = post_order(tree)
    retarray = []
    for x in po_list
        push!(retarray,x.name)
    end
    @test retarray == ["B","A","C","D","E","F","G"]

    retarray2 = GeneralNode[]
    post_order(tree,retarray2)
    bla = []
    for x in retarray2
        push!(bla,x.name)
    end
    @test bla == ["B","A","C","D","E","F","G"]

    prelist = pre_order(tree)
    namelist = []
    for x in prelist
        push!(namelist,x.name)
    end
    @test namelist == ["G", "C", "B", "A", "F", "D", "E"]

    retarray2 = GeneralNode[]
    pre_order(tree,retarray2)
    bla = []
    for x in retarray2
        push!(bla,x.name)
    end
    @test bla == ["G", "C", "B", "A", "F", "D", "E"]

    lo_list = level_order(tree)
    bla = []
    for x in lo_list
        push!(bla,x.name)
    end
    @test bla == ["G", "C", "F", "B", "A", "D", "E"]

    lvl1 = GeneralNode[]
    lvl2 = GeneralNode[]
    lvl3 = GeneralNode[]

    MCPhyloTree.level_traverse(tree,1,lvl1)
    MCPhyloTree.level_traverse(tree,2,lvl2)
    MCPhyloTree.level_traverse(tree,3,lvl3)

    namelist = []
    for x in lvl1
        push!(namelist,x.name)
    end
    @test namelist == ["G"]

    namelist = []
    for x in lvl2
        push!(namelist,x.name)
    end
    @test namelist == ["C", "F"]

    namelist = []
    for x in lvl3
        push!(namelist,x.name)
    end
    @test namelist == ["B", "A", "D", "E"]

    leafret = GeneralNode[]
    get_leaves(tree,leafret)
    leafnames = []
    for x in leafret
        push!(leafnames,x.name)
    end
    @test leafnames==["B","A","D","E"]

    leafret = GeneralNode[]
    get_leaves(tree,leafret)
    leafnames = []
    for x in leafret
        push!(leafnames,x.name)
    end
    @test leafnames == ["B","A","D","E"]
end


text:
mutable struct ClassicPatch{T<:IO,C<:Codec} <: Patch
    io::T
    new_size::Int64
    ctrl::TranscodingStream{C,IOBuffer}
    diff::TranscodingStream{C,IOBuffer}
    data::TranscodingStream{C,IOBuffer}
end
function ClassicPatch(patch_io::IO, new_size::Int64 = typemax(Int64))
    ctrl = TranscodingStream(compressor(), IOBuffer())
    diff = TranscodingStream(compressor(), IOBuffer())
    data = TranscodingStream(compressor(), IOBuffer())
    patch = ClassicPatch(patch_io, new_size, ctrl, diff, data)
    finalizer(finalize_patch, patch)
    return patch
end

format_magic(::Type{ClassicPatch}) = "BSDIFF40"

function write_start(
    ::Type{ClassicPatch},
    patch_io::IO,
    old_data::AbstractVector{UInt8},
    new_data::AbstractVector{UInt8},
)
    ClassicPatch(patch_io, length(new_data))
end

function read_start(::Type{ClassicPatch}, patch_io::IO)
    ctrl_size = read_int(patch_io)
    diff_size = read_int(patch_io)
    new_size  = read_int(patch_io)
    ctrl_io = IOBuffer(read(patch_io, ctrl_size))
    diff_io = IOBuffer(read(patch_io, diff_size))
    data_io = IOBuffer(read(patch_io))
    ctrl = TranscodingStream(decompressor(), ctrl_io)
    diff = TranscodingStream(decompressor(), diff_io)
    data = TranscodingStream(decompressor(), data_io)
    patch = ClassicPatch(patch_io, new_size, ctrl, diff, data)
    finalizer(finalize_patch, patch)
    return patch
end

function finalize_patch(patch::ClassicPatch)
    for stream in (patch.ctrl, patch.diff, patch.data)
        # must be called to avoid leaking memory
        TranscodingStreams.changemode!(stream, :close)
    end
end

function write_finish(patch::ClassicPatch)
    for stream in (patch.ctrl, patch.diff, patch.data)
        write(stream, TranscodingStreams.TOKEN_END)
    end
    write_int(patch.io, patch.ctrl.stream.size)
    write_int(patch.io, patch.diff.stream.size)
    write_int(patch.io, patch.new_size)
    for stream in (patch.ctrl, patch.diff, patch.data)
        write(patch.io, resize!(stream.stream.data, stream.stream.size))
    end
    flush(patch.io)
    finalize(patch)
end

function encode_control(
    patch::ClassicPatch,
    diff_size::Int,
    copy_size::Int,
    skip_size::Int,
)
    write_int(patch.ctrl, diff_size)
    write_int(patch.ctrl, copy_size)
    write_int(patch.ctrl, skip_size)
end

function decode_control(patch::ClassicPatch)
    eof(patch.ctrl) && return nothing
    diff_size = read_int(patch.ctrl)
    copy_size = read_int(patch.ctrl)
    skip_size = read_int(patch.ctrl)
    return diff_size, copy_size, skip_size
end

function encode_diff(
    patch::ClassicPatch,
    diff_size::Int,
    new::AbstractVector{UInt8}, new_pos::Int,
    old::AbstractVector{UInt8}, old_pos::Int,
)
    for i = 1:diff_size
        write(patch.diff, new[new_pos + i] - old[old_pos + i])
    end
end

function decode_diff(
    patch::ClassicPatch,
    diff_size::Int,
    new::IO,
    old::AbstractVector{UInt8},
    old_pos::Int,
)
    for i = 1:diff_size
        write(new, old[old_pos + i] + read(patch.diff, UInt8))
    end
end

function encode_data(
    patch::ClassicPatch,
    copy_size::Int,
    new::AbstractVector{UInt8}, pos::Int,
)
    for i = 1:copy_size
        write(patch.data, new[pos + i])
    end
end

function decode_data(
    patch::ClassicPatch,
    copy_size::Int,
    new::IO,
)
    for i = 1:copy_size
        write(new, read(patch.data, UInt8))
    end
end


text:
## Data-driven SMC algorithm

@gen function migrate_proposal(prev_trace, t, obs)
    model = get_retval(prev_trace)
    # Propose to migration locations for each observed location
    for i in 1:nagents(model)
        loc_addr = :step => t => :obs => :location => i
        if !has_value(obs, loc_addr) continue end
        probs = [c == obs[loc_addr] ? .99 : .01/(model.C-1) for c in 1:model.C]
        {:step => t => :agents => :step! => :agent_step! => 1 =>
         i => :migrate! => :m} ~ Gen.categorical(probs)
    end
end

function data_driven_smc(T::Int, observations::ChoiceMap, n_particles::Int;
                         drift::Bool = false, obs_noise::Float64=5.0)
    # Initialize filter
    state = pf_initialize(bayesian_sir, (0, obs_noise), choicemap(), n_particles)
    model = get_retval(state.traces[1])
    println("Running SMC with data-driven migration proposal...")
    @showprogress for t=1:T
        obs = get_selected(observations, Gen.select(:step => t))
        # Apply drift rejuvenation.
        if drift && effective_sample_size(state) < 0.25 * n_particles
            pf_resample!(state, :residual)
            pf_rejuvenate!(state, mh, (gaussian_drift, ()))
        end
        # Upsample transmission counts based on number newly infected
        pf_update!(state, (t, obs_noise), (UnknownChange(), NoChange()),
                   obs, migrate_proposal, (t, obs))
    end
    # Return particles and log ML estimate
    return state
end


text:
module dispatch

function fn1(dt::Array{Float64,1})
    info("In the fn1")
end

function fn2(dt::Array{Float64,2})
    info("In the fn2")
end

end


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

include("model_Location.jl")
include("model_LocationListResult.jl")
include("model_Subscription.jl")
include("model_SubscriptionListResult.jl")
include("model_SubscriptionPolicies.jl")
include("model_TenantIdDescription.jl")
include("model_TenantListResult.jl")


text:
import Distributions: Exponential, Normal
import JSON: json

count = 1000

"""
generate data that is designed to look a certain way
"""
samples = map(_ -> begin
  x0 = rand(Exponential(1.0))
  x1 = x0 + 32.0
  k = rand() > 0.5
  x2 = k ? x0 + rand(Normal(1.0, 3.0)) : -2 + sin(x0) * x1
  x3 = x0 + sin(x1^2 - x1 + x2) + rand(Normal(1.0, 4.0))
  [x0, x1, x2, x3]
end, 1:count)

write("data.json", json(samples))

max_stars_repo_name: mvarble/react-marginals

max_issues_repo_name: mvarble/react-marginals

max_forks_repo_name: mvarble/react-marginals


text:
using JuLIP
using Test
using LinearAlgebra: I

h3("check that `bulk` evaluates ok...")
at = bulk(:Si)
println(@test typeof(at) == Atoms{Float64})

h3("... and that we can repeat it.")
println(@test length(at * (1,2,3)) == 6 * length(at))

h3("check deepcopy and == ...")
at = bulk(:Si)
at1 = deepcopy(at)
println(@test at == at1)

h3("Check setindex! and getindex ...")
at = bulk(:Si, cubic=true)
x = at[2]
println(@test x == JVec(1.3575, 1.3575, 1.3575))
x += 0.1
at[2] = x
X = positions(at)
println(@test !(X === at.X))
println(@test X[2] == x)

h3("set_positions ...")
at = bulk(:Si, cubic=true) * 2
X = positions(at)
println(@test !(X === at.X))
println(@test X == at.X)
X += 0.1 * rand(JVecF, length(at))
println(@test X != at.X)
set_positions!(at, X)
println(@test X == at.X)
println(@test !(X === at.X))

# TODO: set_momenta, set_masses, set_numbers

println(@test chemical_symbols(at) == fill(:Si, 64))
println(@test chemical_symbols(bulk(:W)) == [:W])

h3("test set_positions!")
Y = copy(X)
Y[3] = JVec(rand(3))
println(@test Y != positions(at))
set_positions!(at, Y)
println(@test Y == positions(at))

h3("test set_momenta!")
Nat = length(at)
Ifree = rand(collect(1:Nat), Nat * 4 ÷ 5) |> unique |> sort # prepare for test below
Iclamp = setdiff(1:Nat, Ifree)
P = rand(size(Y |> mat)...)
P[:, Iclamp] .= 0.0
P = vecs(P)
set_momenta!(at, P)
println(@test P == momenta(at))

h3("test set_dofs!, etc")
# this is making an assumptions on the ordering of dofs; since a new
# implementation of the DOF manager could change this, this test needs to be
# re-implemented if that happens.
set_free!(at, Ifree)
println(@test dofs(at) == position_dofs(at) == mat(Y)[:, Ifree][:])
println(@test momentum_dofs(at) == mat(P)[:, Ifree][:])
q = position_dofs(at)
q = rand(length(q))
set_position_dofs!(at, q)
println(@test q == position_dofs(at))
p = rand(length(q))
set_momentum_dofs!(at, p)
println(@test p == momentum_dofs(at))


text:
xRange, yRange = [parse.(Int, (split(x[3:end], "..")...,)) for x in split.(readline("input.txt")[14:end], ", ")]

#! Part 1

vy = -yRange[1]-1
yMax = vy * (vy + 1)/2

#! Part 2

# Sorta brute force 😅

function doesItHit(velocityX, velocityY)
    positionX = 0
    positionY = 0
    while true
        positionX += velocityX
        positionY += velocityY

        if xRange[1] <= positionX <= xRange[2] && yRange[1] <= positionY <= yRange[2]
            return true
        end

        if positionY < yRange[1] || positionX > xRange[2]
            return false
        end

        velocityY -= 1
        velocityX -= velocityX == 0 ? 0 : abs(velocityX)/velocityX
    end
end

function countVelocities()
    a = max(abs.(xRange)...)
    b = max(abs.(yRange)...)
    return count([doesItHit(x, y) for x in -a:a for y in -b:b])
end

countVelocities()

max_stars_repo_path: Day 17/main.jl

max_issues_repo_path: Day 17/main.jl

max_forks_repo_path: Day 17/main.jl


text:
# Interaction with other languages

mycompare(a, b)::Cint = sign(a - b)

#+

mycompare_c = @cfunction(mycompare, Cint, (Ref{Cdouble}, Ref{Cdouble}))

#+

A = [1.3, -2.7, 4.4, 3.1];

ccall(:qsort, Cvoid, (Ptr{Cdouble}, Csize_t, Csize_t, Ptr{Cvoid}),
      A, length(A), sizeof(eltype(A)), mycompare_c)

#+

A

# ## RCall.jl

using RCall, RDatasets
mtcars = dataset("datasets", "mtcars");

#+

R"""
library(ggplot2)
ggplot($mtcars, aes(x = WT, y = MPG)) + geom_point()
"""


text:
using Zeros

struct Complex2{R,I} <: Number
  re::R
  im::I
end

Complex2(::Zero, ::Zero) = Zero()

Complex2(x::Real) = Complex2(x, Zero())
Complex2(x::Number) = Complex2(real(x), imag(x))
Complex2{R,I}(x::Number) where {R,I} = Complex2(R(real(x)), I(imag(x)))

Base.Complex(x::Complex2) = Complex(x.re, x.im)

Base.convert(::Type{T}, x::Complex2) where {T<:Real} = x.im == 0 ? T(x.re) : throw(InexactError(Symbol(T), T, x))
Base.convert(::Type{Complex{T}}, x::Complex2) where {T<:Real} = Complex(T(x.re), T(x.im))

Base.promote_rule(::Type{Complex2{R,I}}, ::Type{Complex2{S,J}}) where {R<:Real,I<:Real,S<:Real,J<:Real} =
    Complex2{promote_type(R,S), promote_type(I,J)}
Base.promote_rule(::Type{Complex2{R,I}}, ::Type{S}) where {R<:Real,I<:Real,S<:Real} =
    Complex2{promote_type(R,S), promote_type(I,Zero)}
Base.promote_rule(::Type{Complex2{R,I}}, ::Type{Complex{S}}) where {R<:Real,I<:Real,S<:Real} =
    Complex2{promote_type(R,S), promote_type(I,S)}

# Avoid unnecessary/incorrect promotion
for op in (:+, :-, :*, :/)
    @eval begin
        Base.$op(x::Real, y::Complex2) = $op(Complex2(x, Zero()), y)
        Base.$op(x::Complex2, y::Real) = $op(x, Complex2(y, Zero()))
    end
end

Base.real(x::Complex2) = x.re
Base.imag(x::Complex2) = x.im
Base.conj(x::Complex2) = Complex2(x.re, -x.im)
Base.adjoint(x::Complex2) = conj(x)
Base.abs(x::Complex2) = hypot(x.re, x.im)
Base.abs2(x::Complex2) = x.re*x.re + x.im*x.im
Base.angle(x::Complex2) = atan(z.im, z.re)
# Many more functions could be copied from base/complex.jl with no or minimal modifications

Base.:+(x::Complex2, y::Complex2) = Complex2(x.re + y.re, x.im + y.im)
Base.:-(x::Complex2, y::Complex2) = Complex2(x.re - y.re, x.im - y.im)
Base.:-(x::Complex2) = Complex2(-x.re, -x.im)

Base.:*(x::Complex2, y::Complex2) = Complex2(x.re*y.re - x.im*y.im, x.re*y.im + x.im*y.re)

const MyReal = Complex2{<:Real,Zero}
const Imag2 = Complex2{Zero,<:Real}

Base.:/(x::Complex2, y::MyReal) = Complex2(x.re/y.re, x.im/y.re)
Base.:/(x::Complex2, y::Imag2) = Complex2(x.im/y.im, -x.re/y.im)
Base.:/(x::Complex2, y::Complex2) = (x*y')/abs2(y)

const i = Complex2(Zero(), One())

using Test
@test sizeof(i) == 0
@test sizeof(im) == 2
@test sizeof(1*i) == sizeof(1*im)/2
@test i*i isa MyReal
@test 1/i isa Imag2
@test abs(-7i) === 7
@test abs(-7im) === 7.0

using BenchmarkTools
f(x,y) = x + y*im
g(x,y) = x + y*i
a = one(BigFloat)
b = one(BigFloat)
@btime f($a, $b)
@btime g($a, $b)


text:
# make BandwidthBenchmark loadable
push!(LOAD_PATH, joinpath(@__DIR__, "../.."))

using BandwidthBenchmark
using LIKWID

bwbench(; verbose=true)

max_stars_repo_name: carstenbauer/BandwidthBenchmark

max_issues_repo_name: carstenbauer/BandwidthBenchmark

max_forks_repo_name: carstenbauer/BandwidthBenchmark


text:
#jl #! format: off
# # 5-bus Market simulation with [PowerSimulations.jl](https://github.com/NREL-SIIP/PowerSimulations.jl)

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# PowerSimulations.jl supports simulations that consist of sequential optimization problems
# where results from previous problems inform subsequent problems in a variety of ways. This
# example demonstrates some of these capabilities to represent electricity market clearing.

# ## Dependencies and Data
# First, let's create `System`s to represent the Day-Ahead and Real-Time market clearing
# process with hourly, and 5-minute time series data, respectively.

# ### Modeling Packages
using PowerSystems
using PowerSimulations
using PowerSystemCaseBuilder

# ### Data management packages
using Dates
using DataFrames

# ### Optimization packages
using Cbc # mip solver
solver = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 1, "ratioGap" => 0.5)
using Ipopt # solver that supports duals
ipopt_solver = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)

# ### 5-bus Data
# The five bus system data here includes hourly day-ahead data, 5-minute real-time market
# data, and 6-second actual data. We'll only use the hourly and 5-minute data for the
# example simulations below, but the 6-second data is included for future development.
sys_DA = build_system(SIIPExampleSystems, "5_bus_matpower_DA")
sys_RT = build_system(SIIPExampleSystems, "5_bus_matpower_RT")

# ## `OperationsProblemTemplate`s

template_uc = template_unit_commitment()
template_ed = template_economic_dispatch()

# ### Define the Simulation Sequence

problems = SimulationProblems(
    UC = OperationsProblem(
        template_uc,
        sys_DA,
        optimizer = solver,
        balance_slack_variables = true,
    ),
    ED = OperationsProblem(
        template_ed,
        sys_RT,
        optimizer = ipopt_solver,
        constraint_duals = [:CopperPlateBalance],
    ),
)

feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24))

feedforward = Dict(
    ("ED", :devices, :ThermalStandard) => SemiContinuousFF(
        binary_source_problem = ON,
        affected_variables = [ACTIVE_POWER],
    ),
)

#cache = Dict("UC" => [TimeStatusChange(ThermalStandard, PSI.ON)])
intervals = Dict("UC" => (Hour(24), Consecutive()), "ED" => (Hour(1), Consecutive()))

DA_RT_sequence = SimulationSequence(
    problems = problems,
    intervals = intervals,
    ini_cond_chronology = InterProblemChronology(),
    feedforward_chronologies = feedforward_chronologies,
    feedforward = feedforward,
)

# ## `Simulation`
file_path = mkpath(joinpath(".", "5-bus-simulation"))
sim = Simulation(
    name = "5bus-test",
    steps = 1,
    problems = problems,
    sequence = DA_RT_sequence,
    simulation_folder = file_path,
)

# ### Build simulation

build!(sim)

# ### Execute simulation
# ```julia
execute!(sim, enable_progress_bar = false)
# ```

## Results
# First we can load the result metadata
# ```julia
# results = SimulationResults(sim);
# uc_results = get_problem_results(results, "UC")
# ed_results = get_problem_results(results, "ED");
# ```

# Then we can read and examine the results of interest
# ```julia
# prices = read_dual(ed_results, :CopperPlateBalance)
# ```

# or if we want to look at the realized values
# ```julia
# read_realized_duals(ed_results)[:CopperPlateBalance]
# ```

# *note that in this simulation the prices are all equal to the balance slack
# penalty value of $100000/MWh because there is unserved energy in the result*


text:
# Fixes GR warnings in Travis
ENV["GKSwstype"] = "100"

using Documenter, Omniscape

const formats = Any[
    Documenter.HTML(
    	assets = [
    		"assets/custom.css"
    	],
	edit_link = :commit,
    ),
]

makedocs(
	format = formats,
    modules = [Omniscape],
    authors = "Vincent A. Landau",
    sitename = "Omniscape.jl",
    pages = ["About" => "index.md",
			 "How It Works" => "algorithm.md",
    		 "User Guide" => "usage.md",
			 "Examples" => "examples.md",
			 "API Documentation" => "apidocs.md"],
)

deploydocs(
    repo = "github.com/Circuitscape/Omniscape.jl.git",
    devbranch = "main",
    devurl = "latest"
)


text:
function mean_of_zvmean(m)
  nsets = length(m)
  lzv = isnan((m[1].means)[1, 2]) ? false : true
  qzv = isnan((m[1].means)[1, 3]) ? false : true
  cols = (lzv ? (qzv ? (1:3) : (1:2)) : (qzv ? (1:2:3) : (1:1)))
  nmethods = length(cols)
  mm = zeros(Float64, size(m[1].means, 1), nmethods)

  for i in 1:nsets
    mm += (m[i].means)[:, cols]
  end

  mm/nsets
end

function var_of_zvmean(m; means::Union(Matrix{Float64}, Nothing)=nothing)
  nsets = length(m)
  lzv = isnan((m[1].means)[1, 2]) ? false : true
  qzv = isnan((m[1].means)[1, 3]) ? false : true
  cols = (lzv ? (qzv ? (1:3) : (1:2)) : (qzv ? (1:2:3) : (1:1)))
  nmethods = length(cols)
  v = zeros(Float64, size(m[1].means, 1), nmethods)

  mm = means == nothing ? mean_of_zvmean(m) : means
  for i in 1:nsets
    v += ((m[i].means)[:, cols]-mm).^2
  end

  v/(nsets-1)
end

function vrf(v)
  npars, nmethods = size(v)
  vrfs = Array(Float64, npars, nmethods-1)

  for i in 1:nmethods-1
    vrfs[:, i] = v[:, 1]./v[:, i+1]
  end

  vrfs
end

collect_diagnostic(m; diagnostic::ASCIIString="rate") = [(m[i].diagnostics)[diagnostic] for i in 1:length(m)]


text:
macro prog_print(verbose, args...)
    quote
        if $(esc(verbose))
            printstyled(stderr, $(esc.(args)...); color=:green)
        end
    end
end
macro prog_println(verbose, args...)
    quote
        if $(esc(verbose))
            printstyled(stderr, $(esc.(args)...), '\n'; color=:green)
        end
    end
end

simplex_name(::Type{<:Simplex{1}}) = "edges"
simplex_name(::Type{<:Simplex{2}}) = "triangles"
simplex_name(::Type{<:Simplex{3}}) = "tetrahedra"
simplex_name(::Type{<:AbstractSimplex{D}}) where {D} = "$D-simplices"
simplex_name(::Type{<:AbstractCell{D}}) where {D} = "$D-cells"

function fmt_number(i)
    # Stolen from Humanize.jl
    value = string(i)
    group_ends = reverse(collect(length(value):-3:1))
    groups = [value[max(end_index - 2, 1):end_index] for end_index in group_ends]
    return join(groups, ",")
end


text:
is_xlsxfile(f)::Bool = endswith(f, r".xlsx|.xlsm")
is_jsonfile(f)::Bool = endswith(f, ".json")

isnull(x) = ismissing(x) | isnothing(x)
function skipnothing(itr)
    [x for x in itr if !isnothing(x)]
end
function skipnull(itr)
    skipnothing(skipmissing(itr))
end
function issamedata(a, b)::Bool
    isequal(hash(a), hash(b))
end

function namedtuple(d::AbstractDict{String,T}) where T
    key = (Symbol.(collect(keys(d)))...,)
    val = (collect(values(d))...,)
    NamedTuple{key}(val)
end
function namedtuple(d::AbstractDict{Symbol,T}) where T
    key = (collect(keys(d))...,)
    val = (collect(values(d))...,)
    NamedTuple{key}(val)
end

function print_write_result(path, msg="다음과 같습니다")
    print_section("$(msg)\n   SAVED => $(normpath(path))", "연산결과"; color=:green)

    nothing
end

function print_section(message, title="NOTE"; color=:normal)
    msglines = split(chomp(string(message)), '\n')

    for (i, el) in enumerate(msglines)
        prefix = length(msglines) == 1 ? "[ $title: " :
                                i == 1 ? "┌ $title: " :
                                el == last(msglines) ? "└ " : "│ "

        printstyled(stderr, prefix; color=color)
        print(stderr, el)
        print(stderr,  '\n')
    end
    nothing
end

function throw_fuzzylookupname(keyset, idx; kwargs...)
    throw_fuzzylookupname(collect(keyset), idx; kwargs...)
end

function throw_fuzzylookupname(names::AbstractArray, idx::AbstractString;
                                msg="'$idx' not found")
    l = Dict{AbstractString,Int}(zip(names, eachindex(names)))
    candidates = XLSXasJSON.fuzzymatch(l, idx)
    if isempty(candidates)
        throw(ArgumentError(msg))
    end
    candidatesstr = join(string.("\"", candidates, "\""), ", ", " and ")
    throw(ArgumentError(msg * "\nexisting most similar names are: $candidatesstr"))
end


text:
using Documenter
using LightOSM

makedocs(
    sitename="LightOSM.jl Documentation",
    # format=Documenter.HTML(prettyurls=false),
    pages=[
        "Home" => "index.md",
        "Interface" => [
            "types.md",
            "download_network.md",
            "create_graph.md",
            "shortest_path.md",
            "nearest_node.md",
            "download_buildings.md",
            "create_buildings.md",
            "geolocation.md"
        ],
        "Unit Test Use" => "testing_use.md",
    ]
)

deploydocs(
    repo="github.com/DeloitteDigitalAPAC/LightOSM.jl.git",
    devurl="docs",
    push_preview=true,
)


text:
@testset "rule2" begin
    @syms u x z y

    # Example from Section 3.4.3 (Pearl, 2009).
    edges = [
        :U => :X,
        :U => :Y,
        :X => :Z,
        :Z => :Y
    ]
    G = graph(edges)
    @test [name2node(G, s) for s in [:U, :X, :Z, :Y]] == [1, 2, 3, 4]
    r = rule2(G)
    @eqtest r(P(y, x + Do(z))) == P(y, x + z)
    @eqtest r(P(z, Do(x))) == P(z, x)

end


text:
############################################################################################
"""
$(TYPEDEF)

Utility struct to help assign boundaries to parameter - keeps scalar parameter constrained.

# Fields
$(TYPEDFIELDS)
"""
struct Constrained{T<:Real} <: AbstractConstraint
    min::T
    max::T
    function Constrained(min::S, max::T, ϵ = 1e-6) where {S<:AbstractFloat, T<:AbstractFloat}
        ArgCheck.@argcheck min < max
        vals = promote(min, max)
        return new{eltype(vals)}(vals[1] + ϵ, vals[2] - ϵ)
    end
end
#=
#!NOTE: Bijectors seems to constrain Int to Float -> so we will remove this option
function Constrained(min::I, max::I) where {I<:Integer}
    ArgCheck.@argcheck min < max
    return Constrained{I}(min, max)
end
=#

############################################################################################
function construct_transform(info::Constrained, val)
    transform = Bijectors.TruncatedBijector(info.min, info.max)
    return transform, Bijectors.inverse(transform)
end

############################################################################################
function _check(
    _rng::Random.AbstractRNG,
    constraint::Constrained,
    val::R,
) where {R<:Real}
    ArgCheck.@argcheck typeof(val) == typeof(constraint.min) == typeof(constraint.max) "Type of Constrained boundaries must match type of val"
    ArgCheck.@argcheck constraint.min < val < constraint.max "val must be between boundaries"
    return true
end

############################################################################################
# Constrained constraint -> just flatten
function construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::U,
    constraint::Constrained,
    x
) where {T<:AbstractFloat,F<:FlattenTypes,U<:UnflattenTypes}
    return construct_flatten(T, flattentype, unflattentype, x)
end

############################################################################################
#=
!NOTE: For this constraint, we use a Bijector as Transformer, so we do not need to add a new functors
1. MyTransformer <: AbstractTransformer, MyInverseTransformer <: AbstractTransformer
2. define a function construct_transform(MyConstraint, val) -> MyTransformer, MyInverseTransformer
3. overload unconstrain and log_abs_det_jac on MyTransformer.
4. overload constrain on MyInverseTransformer.
=#

############################################################################################
#Export
export
    Constrained,
    construct_flatten,
    construct_transform,
    constrain,
    unconstrain,
    log_abs_det_jac


text:
int_rules_7_3_7 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 7.3.7*Miscellaneous*inverse*hyperbolic*tangent =#
    @apply_utils Antiderivative(atanh(~a + ~(b') * (~x) ^ ~n), ~x) => ~x * atanh(~a + ~b * (~x) ^ ~n) - ~b * ~n * Antiderivative((~x) ^ ~n / (((1 - (~a) ^ 2) - 2 * ~a * ~b * (~x) ^ ~n) - (~b) ^ 2 * (~x) ^ (2 * ~n)), ~x) <-- FreeQ([~a, ~b, ~n], ~x)
    @apply_utils Antiderivative(acoth(~a + ~(b') * (~x) ^ ~n), ~x) => ~x * acoth(~a + ~b * (~x) ^ ~n) - ~b * ~n * Antiderivative((~x) ^ ~n / (((1 - (~a) ^ 2) - 2 * ~a * ~b * (~x) ^ ~n) - (~b) ^ 2 * (~x) ^ (2 * ~n)), ~x) <-- FreeQ([~a, ~b, ~n], ~x)
    @apply_utils Antiderivative(atanh(~(a') + ~(b') * (~x) ^ ~(n')) / ~x, ~x) => (1 / 2) * Antiderivative(log(1 + ~a + ~b * (~x) ^ ~n) / ~x, ~x) - (1 / 2) * Antiderivative(log((1 - ~a) - ~b * (~x) ^ ~n) / ~x, ~x) <-- FreeQ([~a, ~b, ~n], ~x)
    @apply_utils Antiderivative(acoth(~(a') + ~(b') * (~x) ^ ~(n')) / ~x, ~x) => (1 / 2) * Antiderivative(log(1 + 1 / (~a + ~b * (~x) ^ ~n)) / ~x, ~x) - (1 / 2) * Antiderivative(log(1 - 1 / (~a + ~b * (~x) ^ ~n)) / ~x, ~x) <-- FreeQ([~a, ~b, ~n], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * atanh(~a + ~(b') * (~x) ^ ~n), ~x) => ((~x) ^ (~m + 1) * atanh(~a + ~b * (~x) ^ ~n)) / (~m + 1) - ((~b * ~n) / (~m + 1)) * Antiderivative((~x) ^ (~m + ~n) / (((1 - (~a) ^ 2) - 2 * ~a * ~b * (~x) ^ ~n) - (~b) ^ 2 * (~x) ^ (2 * ~n)), ~x) <-- FreeQ([~a, ~b], ~x) && (RationalQ(~m, ~n) && (NeQ(~m, -1) && NeQ(~m + 1, ~n)))
    @apply_utils Antiderivative((~x) ^ ~(m') * acoth(~a + ~(b') * (~x) ^ ~n), ~x) => ((~x) ^ (~m + 1) * acoth(~a + ~b * (~x) ^ ~n)) / (~m + 1) - ((~b * ~n) / (~m + 1)) * Antiderivative((~x) ^ (~m + ~n) / (((1 - (~a) ^ 2) - 2 * ~a * ~b * (~x) ^ ~n) - (~b) ^ 2 * (~x) ^ (2 * ~n)), ~x) <-- FreeQ([~a, ~b], ~x) && (RationalQ(~m, ~n) && (NeQ(~m, -1) && NeQ(~m + 1, ~n)))
    @apply_utils Antiderivative(atanh(~(a') + ~(b') * (~f) ^ (~(c') + ~(d') * ~x)), ~x) => (1 / 2) * Antiderivative(log(1 + ~a + ~b * (~f) ^ (~c + ~d * ~x)), ~x) - (1 / 2) * Antiderivative(log((1 - ~a) - ~b * (~f) ^ (~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~f], ~x)
    @apply_utils Antiderivative(acoth(~(a') + ~(b') * (~f) ^ (~(c') + ~(d') * ~x)), ~x) => (1 / 2) * Antiderivative(log(1 + 1 / (~a + ~b * (~f) ^ (~c + ~d * ~x))), ~x) - (1 / 2) * Antiderivative(log(1 - 1 / (~a + ~b * (~f) ^ (~c + ~d * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~f], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * atanh(~(a') + ~(b') * (~f) ^ (~(c') + ~(d') * ~x)), ~x) => (1 / 2) * Antiderivative((~x) ^ ~m * log(1 + ~a + ~b * (~f) ^ (~c + ~d * ~x)), ~x) - (1 / 2) * Antiderivative((~x) ^ ~m * log((1 - ~a) - ~b * (~f) ^ (~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~f], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~x) ^ ~(m') * acoth(~(a') + ~(b') * (~f) ^ (~(c') + ~(d') * ~x)), ~x) => (1 / 2) * Antiderivative((~x) ^ ~m * log(1 + 1 / (~a + ~b * (~f) ^ (~c + ~d * ~x))), ~x) - (1 / 2) * Antiderivative((~x) ^ ~m * log(1 - 1 / (~a + ~b * (~f) ^ (~c + ~d * ~x))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~f], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative(~(u') * atanh(~(c') / (~(a') + ~(b') * (~x) ^ ~(n'))) ^ ~(m'), ~x) => Antiderivative(~u * acoth(~a / ~c + (~b * (~x) ^ ~n) / ~c) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~m], ~x)
    @apply_utils Antiderivative(~(u') * acoth(~(c') / (~(a') + ~(b') * (~x) ^ ~(n'))) ^ ~(m'), ~x) => Antiderivative(~u * atanh(~a / ~c + (~b * (~x) ^ ~n) / ~c) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~c, ~n, ~m], ~x)
    @apply_utils Antiderivative(atanh((~(c') * ~x) / sqrt(~(a') + ~(b') * (~x) ^ 2)), ~x) => ~x * atanh((~c * ~x) / sqrt(~a + ~b * (~x) ^ 2)) - ~c * Antiderivative(~x / sqrt(~a + ~b * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && EqQ(~b, (~c) ^ 2)
    @apply_utils Antiderivative(acoth((~(c') * ~x) / sqrt(~(a') + ~(b') * (~x) ^ 2)), ~x) => ~x * acoth((~c * ~x) / sqrt(~a + ~b * (~x) ^ 2)) - ~c * Antiderivative(~x / sqrt(~a + ~b * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && EqQ(~b, (~c) ^ 2)
    @apply_utils Antiderivative(atanh((~(c') * ~x) / sqrt(~(a') + ~(b') * (~x) ^ 2)) / ~x, ~x) => atanh((~c * ~x) / sqrt(~a + ~b * (~x) ^ 2)) * log(~x) - ~c * Antiderivative(log(~x) / sqrt(~a + ~b * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && EqQ(~b, (~c) ^ 2)
    @apply_utils Antiderivative(acoth((~(c') * ~x) / sqrt(~(a') + ~(b') * (~x) ^ 2)) / ~x, ~x) => acoth((~c * ~x) / sqrt(~a + ~b * (~x) ^ 2)) * log(~x) - ~c * Antiderivative(log(~x) / sqrt(~a + ~b * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && EqQ(~b, (~c) ^ 2)
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * atanh((~(c') * ~x) / sqrt(~(a') + ~(b') * (~x) ^ 2)), ~x) => ((~d * ~x) ^ (~m + 1) * atanh((~c * ~x) / sqrt(~a + ~b * (~x) ^ 2))) / (~d * (~m + 1)) - (~c / (~d * (~m + 1))) * Antiderivative((~d * ~x) ^ (~m + 1) / sqrt(~a + ~b * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (EqQ(~b, (~c) ^ 2) && NeQ(~m, -1))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * acoth((~(c') * ~x) / sqrt(~(a') + ~(b') * (~x) ^ 2)), ~x) => ((~d * ~x) ^ (~m + 1) * acoth((~c * ~x) / sqrt(~a + ~b * (~x) ^ 2))) / (~d * (~m + 1)) - (~c / (~d * (~m + 1))) * Antiderivative((~d * ~x) ^ (~m + 1) / sqrt(~a + ~b * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (EqQ(~b, (~c) ^ 2) && NeQ(~m, -1))
    @apply_utils Antiderivative(1 / (sqrt(~(a') + ~(b') * (~x) ^ 2) * atanh((~(c') * ~x) / sqrt(~(a') + ~(b') * (~x) ^ 2))), ~x) => (1 / ~c) * log(atanh((~c * ~x) / sqrt(~a + ~b * (~x) ^ 2))) <-- FreeQ([~a, ~b, ~c], ~x) && EqQ(~b, (~c) ^ 2)
    @apply_utils Antiderivative(1 / (sqrt(~(a') + ~(b') * (~x) ^ 2) * acoth((~(c') * ~x) / sqrt(~(a') + ~(b') * (~x) ^ 2))), ~x) => (-1 / ~c) * log(acoth((~c * ~x) / sqrt(~a + ~b * (~x) ^ 2))) <-- FreeQ([~a, ~b, ~c], ~x) && EqQ(~b, (~c) ^ 2)
    @apply_utils Antiderivative(atanh((~(c') * ~x) / sqrt(~(a') + ~(b') * (~x) ^ 2)) ^ ~(m') / sqrt(~(a') + ~(b') * (~x) ^ 2), ~x) => atanh((~c * ~x) / sqrt(~a + ~b * (~x) ^ 2)) ^ (~m + 1) / (~c * (~m + 1)) <-- FreeQ([~a, ~b, ~c, ~m], ~x) && (EqQ(~b, (~c) ^ 2) && NeQ(~m, -1))
    @apply_utils Antiderivative(acoth((~(c') * ~x) / sqrt(~(a') + ~(b') * (~x) ^ 2)) ^ ~(m') / sqrt(~(a') + ~(b') * (~x) ^ 2), ~x) => -(acoth((~c * ~x) / sqrt(~a + ~b * (~x) ^ 2)) ^ (~m + 1)) / (~c * (~m + 1)) <-- FreeQ([~a, ~b, ~c, ~m], ~x) && (EqQ(~b, (~c) ^ 2) && NeQ(~m, -1))
    @apply_utils Antiderivative(atanh((~(c') * ~x) / sqrt(~(a') + ~(b') * (~x) ^ 2)) ^ ~(m') / sqrt(~(d') + ~(e') * (~x) ^ 2), ~x) => (sqrt(~a + ~b * (~x) ^ 2) / sqrt(~d + ~e * (~x) ^ 2)) * Antiderivative(atanh((~c * ~x) / sqrt(~a + ~b * (~x) ^ 2)) ^ ~m / sqrt(~a + ~b * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (EqQ(~b, (~c) ^ 2) && EqQ(~b * ~d - ~a * ~e, 0))
    @apply_utils Antiderivative(acoth((~(c') * ~x) / sqrt(~(a') + ~(b') * (~x) ^ 2)) ^ ~(m') / sqrt(~(d') + ~(e') * (~x) ^ 2), ~x) => (sqrt(~a + ~b * (~x) ^ 2) / sqrt(~d + ~e * (~x) ^ 2)) * Antiderivative(acoth((~c * ~x) / sqrt(~a + ~b * (~x) ^ 2)) ^ ~m / sqrt(~a + ~b * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (EqQ(~b, (~c) ^ 2) && EqQ(~b * ~d - ~a * ~e, 0))
    If(TrueQ($LoadShowSteps), Antiderivative((~u)*(~v)^(~n'), ~x) := With([tmp = InverseFunctionOfLinear(u, x)], ShowStep("", "Antiderivative(f(x,atanh(a+b*x))/(1-(a+b*x)^2),x)", "Subst(Antiderivative(f(-a/b+tanh(x)/b,x),x),x,atanh(a+b*x))/b", Hold( (-Discriminant(v, x)/(4*Coefficient(v, x, 2)))^n/ Coefficient(tmp((1)), x, 1)* Subst( Antiderivative(SimplifyIntegrand( SubstForInverseFunction(u, tmp, x)*sech(x)^(2*(n + 1)), x), x), x, tmp))) <-- Not(FalseQ(tmp)) && EqQ(Head(tmp), atanh) && EqQ(Discriminant(v, x)*tmp((1))^2 - D(v, x)^2, 0)) <-- SimplifyFlag && QuadraticQ(v, x) && ILtQ(n, 0) && PosQ(Discriminant(v, x)) && MatchQ(u, (~r')*(~f)^(~w) <-- FreeQ(f, x)), Antiderivative((~u)*(~v)^(~n'), ~x) := With([tmp = InverseFunctionOfLinear(u, x)], (-Discriminant(v, x)/(4*Coefficient(v, x, 2)))^n/ Coefficient(tmp((1)), x, 1)* Subst( Antiderivative(SimplifyIntegrand( SubstForInverseFunction(u, tmp, x)*sech(x)^(2*(n + 1)), x), x), x, tmp) <-- Not(FalseQ(tmp)) && EqQ(Head(tmp), atanh) && EqQ(Discriminant(v, x)*tmp((1))^2 - D(v, x)^2, 0)) <-- QuadraticQ(v, x) && ILtQ(n, 0) && PosQ(Discriminant(v, x)) && MatchQ(u, (~r')*(~f)^(~w) <-- FreeQ(f, x)))
    If(TrueQ($LoadShowSteps), Antiderivative((~u)*(~v)^(~n'), ~x) := With([tmp = InverseFunctionOfLinear(u, x)], ShowStep("", "Antiderivative(f(x,acoth(a+b*x))/(1-(a+b*x)^2),x)", "Subst(Antiderivative(f(-a/b+coth(x)/b,x),x),x,acoth(a+b*x))/b", Hold( (-Discriminant(v, x)/(4*Coefficient(v, x, 2)))^n/ Coefficient(tmp((1)), x, 1)* Subst( Antiderivative(SimplifyIntegrand( SubstForInverseFunction(u, tmp, x)*(-csch(x)^2)^(n + 1), x), x), x, tmp))) <-- Not(FalseQ(tmp)) && EqQ(Head(tmp), acoth) && EqQ(Discriminant(v, x)*tmp((1))^2 - D(v, x)^2, 0)) <-- SimplifyFlag && QuadraticQ(v, x) && ILtQ(n, 0) && PosQ(Discriminant(v, x)) && MatchQ(u, (~r')*(~f)^(~w) <-- FreeQ(f, x)), Antiderivative((~u)*(~v)^(~n'), ~x) := With([tmp = InverseFunctionOfLinear(u, x)], (-Discriminant(v, x)/(4*Coefficient(v, x, 2)))^n/ Coefficient(tmp((1)), x, 1)* Subst( Antiderivative(SimplifyIntegrand( SubstForInverseFunction(u, tmp, x)*(-csch(x)^2)^(n + 1), x), x), x, tmp) <-- Not(FalseQ(tmp)) && EqQ(Head(tmp), acoth) && EqQ(Discriminant(v, x)*tmp((1))^2 - D(v, x)^2, 0)) <-- QuadraticQ(v, x) && ILtQ(n, 0) && PosQ(Discriminant(v, x)) && MatchQ(u, (~r')*(~f)^(~w) <-- FreeQ(f, x)))
    @apply_utils Antiderivative(atanh(~(c') + ~(d') * tanh(~(a') + ~(b') * ~x)), ~x) => ~x * atanh(~c + ~d * tanh(~a + ~b * ~x)) + ~b * Antiderivative(~x / ((~c - ~d) + ~c * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c - ~d) ^ 2, 1)
    @apply_utils Antiderivative(acoth(~(c') + ~(d') * tanh(~(a') + ~(b') * ~x)), ~x) => ~x * acoth(~c + ~d * tanh(~a + ~b * ~x)) + ~b * Antiderivative(~x / ((~c - ~d) + ~c * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c - ~d) ^ 2, 1)
    @apply_utils Antiderivative(atanh(~(c') + ~(d') * coth(~(a') + ~(b') * ~x)), ~x) => ~x * atanh(~c + ~d * coth(~a + ~b * ~x)) + ~b * Antiderivative(~x / ((~c - ~d) - ~c * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c - ~d) ^ 2, 1)
    @apply_utils Antiderivative(acoth(~(c') + ~(d') * coth(~(a') + ~(b') * ~x)), ~x) => ~x * acoth(~c + ~d * coth(~a + ~b * ~x)) + ~b * Antiderivative(~x / ((~c - ~d) - ~c * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c - ~d) ^ 2, 1)
    @apply_utils Antiderivative(atanh(~(c') + ~(d') * tanh(~(a') + ~(b') * ~x)), ~x) => (~x * atanh(~c + ~d * tanh(~a + ~b * ~x)) + ~b * ((1 - ~c) - ~d) * Antiderivative((~x * ℯ ^ (2 * ~a + 2 * ~b * ~x)) / ((1 - ~c) + ~d + ((1 - ~c) - ~d) * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x)) - ~b * (1 + ~c + ~d) * Antiderivative((~x * ℯ ^ (2 * ~a + 2 * ~b * ~x)) / (((1 + ~c) - ~d) + (1 + ~c + ~d) * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~c - ~d) ^ 2, 1)
    @apply_utils Antiderivative(acoth(~(c') + ~(d') * tanh(~(a') + ~(b') * ~x)), ~x) => (~x * acoth(~c + ~d * tanh(~a + ~b * ~x)) + ~b * ((1 - ~c) - ~d) * Antiderivative((~x * ℯ ^ (2 * ~a + 2 * ~b * ~x)) / ((1 - ~c) + ~d + ((1 - ~c) - ~d) * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x)) - ~b * (1 + ~c + ~d) * Antiderivative((~x * ℯ ^ (2 * ~a + 2 * ~b * ~x)) / (((1 + ~c) - ~d) + (1 + ~c + ~d) * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~c - ~d) ^ 2, 1)
    @apply_utils Antiderivative(atanh(~(c') + ~(d') * coth(~(a') + ~(b') * ~x)), ~x) => (~x * atanh(~c + ~d * coth(~a + ~b * ~x)) + ~b * (1 + ~c + ~d) * Antiderivative((~x * ℯ ^ (2 * ~a + 2 * ~b * ~x)) / (((1 + ~c) - ~d) - (1 + ~c + ~d) * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x)) - ~b * ((1 - ~c) - ~d) * Antiderivative((~x * ℯ ^ (2 * ~a + 2 * ~b * ~x)) / (((1 - ~c) + ~d) - ((1 - ~c) - ~d) * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~c - ~d) ^ 2, 1)
    @apply_utils Antiderivative(acoth(~(c') + ~(d') * coth(~(a') + ~(b') * ~x)), ~x) => (~x * acoth(~c + ~d * coth(~a + ~b * ~x)) + ~b * (1 + ~c + ~d) * Antiderivative((~x * ℯ ^ (2 * ~a + 2 * ~b * ~x)) / (((1 + ~c) - ~d) - (1 + ~c + ~d) * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x)) - ~b * ((1 - ~c) - ~d) * Antiderivative((~x * ℯ ^ (2 * ~a + 2 * ~b * ~x)) / (((1 - ~c) + ~d) - ((1 - ~c) - ~d) * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~c - ~d) ^ 2, 1)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * atanh(~(c') + ~(d') * tanh(~(a') + ~(b') * ~x)), ~x) => ((~e + ~f * ~x) ^ (~m + 1) * atanh(~c + ~d * tanh(~a + ~b * ~x))) / (~f * (~m + 1)) + (~b / (~f * (~m + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m + 1) / ((~c - ~d) + ~c * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~m, 0) && EqQ((~c - ~d) ^ 2, 1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * acoth(~(c') + ~(d') * tanh(~(a') + ~(b') * ~x)), ~x) => ((~e + ~f * ~x) ^ (~m + 1) * acoth(~c + ~d * tanh(~a + ~b * ~x))) / (~f * (~m + 1)) + (~b / (~f * (~m + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m + 1) / ((~c - ~d) + ~c * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~m, 0) && EqQ((~c - ~d) ^ 2, 1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * atanh(~(c') + ~(d') * coth(~(a') + ~(b') * ~x)), ~x) => ((~e + ~f * ~x) ^ (~m + 1) * atanh(~c + ~d * coth(~a + ~b * ~x))) / (~f * (~m + 1)) + (~b / (~f * (~m + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m + 1) / ((~c - ~d) - ~c * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~m, 0) && EqQ((~c - ~d) ^ 2, 1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * acoth(~(c') + ~(d') * coth(~(a') + ~(b') * ~x)), ~x) => ((~e + ~f * ~x) ^ (~m + 1) * acoth(~c + ~d * coth(~a + ~b * ~x))) / (~f * (~m + 1)) + (~b / (~f * (~m + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m + 1) / ((~c - ~d) - ~c * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~m, 0) && EqQ((~c - ~d) ^ 2, 1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * atanh(~(c') + ~(d') * tanh(~(a') + ~(b') * ~x)), ~x) => (((~e + ~f * ~x) ^ (~m + 1) * atanh(~c + ~d * tanh(~a + ~b * ~x))) / (~f * (~m + 1)) + ((~b * ((1 - ~c) - ~d)) / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * ℯ ^ (2 * ~a + 2 * ~b * ~x)) / ((1 - ~c) + ~d + ((1 - ~c) - ~d) * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x)) - ((~b * (1 + ~c + ~d)) / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * ℯ ^ (2 * ~a + 2 * ~b * ~x)) / (((1 + ~c) - ~d) + (1 + ~c + ~d) * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~m, 0) && NeQ((~c - ~d) ^ 2, 1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * acoth(~(c') + ~(d') * tanh(~(a') + ~(b') * ~x)), ~x) => (((~e + ~f * ~x) ^ (~m + 1) * acoth(~c + ~d * tanh(~a + ~b * ~x))) / (~f * (~m + 1)) + ((~b * ((1 - ~c) - ~d)) / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * ℯ ^ (2 * ~a + 2 * ~b * ~x)) / ((1 - ~c) + ~d + ((1 - ~c) - ~d) * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x)) - ((~b * (1 + ~c + ~d)) / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * ℯ ^ (2 * ~a + 2 * ~b * ~x)) / (((1 + ~c) - ~d) + (1 + ~c + ~d) * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~m, 0) && NeQ((~c - ~d) ^ 2, 1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * atanh(~(c') + ~(d') * coth(~(a') + ~(b') * ~x)), ~x) => (((~e + ~f * ~x) ^ (~m + 1) * atanh(~c + ~d * coth(~a + ~b * ~x))) / (~f * (~m + 1)) + ((~b * (1 + ~c + ~d)) / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * ℯ ^ (2 * ~a + 2 * ~b * ~x)) / (((1 + ~c) - ~d) - (1 + ~c + ~d) * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x)) - ((~b * ((1 - ~c) - ~d)) / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * ℯ ^ (2 * ~a + 2 * ~b * ~x)) / (((1 - ~c) + ~d) - ((1 - ~c) - ~d) * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~m, 0) && NeQ((~c - ~d) ^ 2, 1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * acoth(~(c') + ~(d') * coth(~(a') + ~(b') * ~x)), ~x) => (((~e + ~f * ~x) ^ (~m + 1) * acoth(~c + ~d * coth(~a + ~b * ~x))) / (~f * (~m + 1)) + ((~b * (1 + ~c + ~d)) / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * ℯ ^ (2 * ~a + 2 * ~b * ~x)) / (((1 + ~c) - ~d) - (1 + ~c + ~d) * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x)) - ((~b * ((1 - ~c) - ~d)) / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * ℯ ^ (2 * ~a + 2 * ~b * ~x)) / (((1 - ~c) + ~d) - ((1 - ~c) - ~d) * ℯ ^ (2 * ~a + 2 * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~m, 0) && NeQ((~c - ~d) ^ 2, 1))
    @apply_utils Antiderivative(atanh(tan(~(a') + ~(b') * ~x)), ~x) => ~x * atanh(tan(~a + ~b * ~x)) - ~b * Antiderivative(~x * sec(2 * ~a + 2 * ~b * ~x), ~x) <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(acoth(tan(~(a') + ~(b') * ~x)), ~x) => ~x * acoth(tan(~a + ~b * ~x)) - ~b * Antiderivative(~x * sec(2 * ~a + 2 * ~b * ~x), ~x) <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(atanh(cot(~(a') + ~(b') * ~x)), ~x) => ~x * atanh(cot(~a + ~b * ~x)) - ~b * Antiderivative(~x * sec(2 * ~a + 2 * ~b * ~x), ~x) <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(acoth(cot(~(a') + ~(b') * ~x)), ~x) => ~x * acoth(cot(~a + ~b * ~x)) - ~b * Antiderivative(~x * sec(2 * ~a + 2 * ~b * ~x), ~x) <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * atanh(tan(~(a') + ~(b') * ~x)), ~x) => ((~e + ~f * ~x) ^ (~m + 1) * atanh(tan(~a + ~b * ~x))) / (~f * (~m + 1)) - (~b / (~f * (~m + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m + 1) * sec(2 * ~a + 2 * ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * acoth(tan(~(a') + ~(b') * ~x)), ~x) => ((~e + ~f * ~x) ^ (~m + 1) * acoth(tan(~a + ~b * ~x))) / (~f * (~m + 1)) - (~b / (~f * (~m + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m + 1) * sec(2 * ~a + 2 * ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * atanh(cot(~(a') + ~(b') * ~x)), ~x) => ((~e + ~f * ~x) ^ (~m + 1) * atanh(cot(~a + ~b * ~x))) / (~f * (~m + 1)) - (~b / (~f * (~m + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m + 1) * sec(2 * ~a + 2 * ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * acoth(cot(~(a') + ~(b') * ~x)), ~x) => ((~e + ~f * ~x) ^ (~m + 1) * acoth(cot(~a + ~b * ~x))) / (~f * (~m + 1)) - (~b / (~f * (~m + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m + 1) * sec(2 * ~a + 2 * ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative(atanh(~(c') + ~(d') * tan(~(a') + ~(b') * ~x)), ~x) => ~x * atanh(~c + ~d * tan(~a + ~b * ~x)) + im * ~b * Antiderivative(~x / (~c + im * ~d + ~c * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c + im * ~d) ^ 2, 1)
    @apply_utils Antiderivative(acoth(~(c') + ~(d') * tan(~(a') + ~(b') * ~x)), ~x) => ~x * acoth(~c + ~d * tan(~a + ~b * ~x)) + im * ~b * Antiderivative(~x / (~c + im * ~d + ~c * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c + im * ~d) ^ 2, 1)
    @apply_utils Antiderivative(atanh(~(c') + ~(d') * cot(~(a') + ~(b') * ~x)), ~x) => ~x * atanh(~c + ~d * cot(~a + ~b * ~x)) + im * ~b * Antiderivative(~x / ((~c - im * ~d) - ~c * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c - im * ~d) ^ 2, 1)
    @apply_utils Antiderivative(acoth(~(c') + ~(d') * cot(~(a') + ~(b') * ~x)), ~x) => ~x * acoth(~c + ~d * cot(~a + ~b * ~x)) + im * ~b * Antiderivative(~x / ((~c - im * ~d) - ~c * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~c - im * ~d) ^ 2, 1)
    @apply_utils Antiderivative(atanh(~(c') + ~(d') * tan(~(a') + ~(b') * ~x)), ~x) => (~x * atanh(~c + ~d * tan(~a + ~b * ~x)) + im * ~b * ((1 - ~c) + im * ~d) * Antiderivative((~x * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)) / (((1 - ~c) - im * ~d) + ((1 - ~c) + im * ~d) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x)) - im * ~b * ((1 + ~c) - im * ~d) * Antiderivative((~x * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)) / (1 + ~c + im * ~d + ((1 + ~c) - im * ~d) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~c + im * ~d) ^ 2, 1)
    @apply_utils Antiderivative(acoth(~(c') + ~(d') * tan(~(a') + ~(b') * ~x)), ~x) => (~x * acoth(~c + ~d * tan(~a + ~b * ~x)) + im * ~b * ((1 - ~c) + im * ~d) * Antiderivative((~x * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)) / (((1 - ~c) - im * ~d) + ((1 - ~c) + im * ~d) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x)) - im * ~b * ((1 + ~c) - im * ~d) * Antiderivative((~x * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)) / (1 + ~c + im * ~d + ((1 + ~c) - im * ~d) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~c + im * ~d) ^ 2, 1)
    @apply_utils Antiderivative(atanh(~(c') + ~(d') * cot(~(a') + ~(b') * ~x)), ~x) => (~x * atanh(~c + ~d * cot(~a + ~b * ~x)) - im * ~b * ((1 - ~c) - im * ~d) * Antiderivative((~x * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)) / (((1 - ~c) + im * ~d) - ((1 - ~c) - im * ~d) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x)) + im * ~b * (1 + ~c + im * ~d) * Antiderivative((~x * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)) / (((1 + ~c) - im * ~d) - (1 + ~c + im * ~d) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~c - im * ~d) ^ 2, 1)
    @apply_utils Antiderivative(acoth(~(c') + ~(d') * cot(~(a') + ~(b') * ~x)), ~x) => (~x * acoth(~c + ~d * cot(~a + ~b * ~x)) - im * ~b * ((1 - ~c) - im * ~d) * Antiderivative((~x * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)) / (((1 - ~c) + im * ~d) - ((1 - ~c) - im * ~d) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x)) + im * ~b * (1 + ~c + im * ~d) * Antiderivative((~x * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)) / (((1 + ~c) - im * ~d) - (1 + ~c + im * ~d) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~c - im * ~d) ^ 2, 1)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * atanh(~(c') + ~(d') * tan(~(a') + ~(b') * ~x)), ~x) => ((~e + ~f * ~x) ^ (~m + 1) * atanh(~c + ~d * tan(~a + ~b * ~x))) / (~f * (~m + 1)) + ((im * ~b) / (~f * (~m + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m + 1) / (~c + im * ~d + ~c * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~m, 0) && EqQ((~c + im * ~d) ^ 2, 1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * acoth(~(c') + ~(d') * tan(~(a') + ~(b') * ~x)), ~x) => ((~e + ~f * ~x) ^ (~m + 1) * acoth(~c + ~d * tan(~a + ~b * ~x))) / (~f * (~m + 1)) + ((im * ~b) / (~f * (~m + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m + 1) / (~c + im * ~d + ~c * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~m, 0) && EqQ((~c + im * ~d) ^ 2, 1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * atanh(~(c') + ~(d') * cot(~(a') + ~(b') * ~x)), ~x) => ((~e + ~f * ~x) ^ (~m + 1) * atanh(~c + ~d * cot(~a + ~b * ~x))) / (~f * (~m + 1)) + ((im * ~b) / (~f * (~m + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m + 1) / ((~c - im * ~d) - ~c * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~m, 0) && EqQ((~c - im * ~d) ^ 2, 1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * acoth(~(c') + ~(d') * cot(~(a') + ~(b') * ~x)), ~x) => ((~e + ~f * ~x) ^ (~m + 1) * acoth(~c + ~d * cot(~a + ~b * ~x))) / (~f * (~m + 1)) + ((im * ~b) / (~f * (~m + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m + 1) / ((~c - im * ~d) - ~c * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~m, 0) && EqQ((~c - im * ~d) ^ 2, 1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * atanh(~(c') + ~(d') * tan(~(a') + ~(b') * ~x)), ~x) => (((~e + ~f * ~x) ^ (~m + 1) * atanh(~c + ~d * tan(~a + ~b * ~x))) / (~f * (~m + 1)) + ((im * ~b * ((1 - ~c) + im * ~d)) / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)) / (((1 - ~c) - im * ~d) + ((1 - ~c) + im * ~d) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x)) - ((im * ~b * ((1 + ~c) - im * ~d)) / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)) / (1 + ~c + im * ~d + ((1 + ~c) - im * ~d) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~m, 0) && NeQ((~c + im * ~d) ^ 2, 1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * acoth(~(c') + ~(d') * tan(~(a') + ~(b') * ~x)), ~x) => (((~e + ~f * ~x) ^ (~m + 1) * acoth(~c + ~d * tan(~a + ~b * ~x))) / (~f * (~m + 1)) + ((im * ~b * ((1 - ~c) + im * ~d)) / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)) / (((1 - ~c) - im * ~d) + ((1 - ~c) + im * ~d) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x)) - ((im * ~b * ((1 + ~c) - im * ~d)) / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)) / (1 + ~c + im * ~d + ((1 + ~c) - im * ~d) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~m, 0) && NeQ((~c + im * ~d) ^ 2, 1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * atanh(~(c') + ~(d') * cot(~(a') + ~(b') * ~x)), ~x) => (((~e + ~f * ~x) ^ (~m + 1) * atanh(~c + ~d * cot(~a + ~b * ~x))) / (~f * (~m + 1)) - ((im * ~b * ((1 - ~c) - im * ~d)) / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)) / (((1 - ~c) + im * ~d) - ((1 - ~c) - im * ~d) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x)) + ((im * ~b * (1 + ~c + im * ~d)) / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)) / (((1 + ~c) - im * ~d) - (1 + ~c + im * ~d) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~m, 0) && NeQ((~c - im * ~d) ^ 2, 1))
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * acoth(~(c') + ~(d') * cot(~(a') + ~(b') * ~x)), ~x) => (((~e + ~f * ~x) ^ (~m + 1) * acoth(~c + ~d * cot(~a + ~b * ~x))) / (~f * (~m + 1)) - ((im * ~b * ((1 - ~c) - im * ~d)) / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)) / (((1 - ~c) + im * ~d) - ((1 - ~c) - im * ~d) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x)) + ((im * ~b * (1 + ~c + im * ~d)) / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)) / (((1 + ~c) - im * ~d) - (1 + ~c + im * ~d) * ℯ ^ (2 * im * ~a + 2 * im * ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~m, 0) && NeQ((~c - im * ~d) ^ 2, 1))
    @apply_utils Antiderivative(atanh(~u), ~x) => ~x * atanh(~u) - Antiderivative(SimplifyIntegrand((~x * D(~u, ~x)) / (1 - (~u) ^ 2), ~x), ~x) <-- InverseFunctionFreeQ(~u, ~x)
    @apply_utils Antiderivative(acoth(~u), ~x) => ~x * acoth(~u) - Antiderivative(SimplifyIntegrand((~x * D(~u, ~x)) / (1 - (~u) ^ 2), ~x), ~x) <-- InverseFunctionFreeQ(~u, ~x)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * atanh(~u)), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * (~a + ~b * atanh(~u))) / (~d * (~m + 1)) - (~b / (~d * (~m + 1))) * Antiderivative(SimplifyIntegrand(((~c + ~d * ~x) ^ (~m + 1) * D(~u, ~x)) / (1 - (~u) ^ 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (NeQ(~m, -1) && (InverseFunctionFreeQ(~u, ~x) && (Not(FunctionOfQ((~c + ~d * ~x) ^ (~m + 1), ~u, ~x)) && FalseQ(PowerVariableExpn(~u, ~m + 1, ~x)))))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * acoth(~u)), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * (~a + ~b * acoth(~u))) / (~d * (~m + 1)) - (~b / (~d * (~m + 1))) * Antiderivative(SimplifyIntegrand(((~c + ~d * ~x) ^ (~m + 1) * D(~u, ~x)) / (1 - (~u) ^ 2), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (NeQ(~m, -1) && (InverseFunctionFreeQ(~u, ~x) && (Not(FunctionOfQ((~c + ~d * ~x) ^ (~m + 1), ~u, ~x)) && FalseQ(PowerVariableExpn(~u, ~m + 1, ~x)))))
    @apply_utils Antiderivative(~v * (~(a') + ~(b') * atanh(~u)), ~x) => With([w = IntHide(~v, ~x)], Dist(~a + ~b * atanh(~u), w, ~x) - ~b * Antiderivative(SimplifyIntegrand((w * D(~u, ~x)) / (1 - (~u) ^ 2), ~x), ~x) <-- InverseFunctionFreeQ(w, ~x)) <-- FreeQ([~a, ~b], ~x) && (InverseFunctionFreeQ(~u, ~x) && (Not(MatchQ(~v, (~(c') + ~(d') * ~x) ^ ~(m') <-- FreeQ([c, d, m], ~x))) && FalseQ(FunctionOfLinear(~v * (~a + ~b * atanh(~u)), ~x))))
    @apply_utils Antiderivative(~v * (~(a') + ~(b') * acoth(~u)), ~x) => With([w = IntHide(~v, ~x)], Dist(~a + ~b * acoth(~u), w, ~x) - ~b * Antiderivative(SimplifyIntegrand((w * D(~u, ~x)) / (1 - (~u) ^ 2), ~x), ~x) <-- InverseFunctionFreeQ(w, ~x)) <-- FreeQ([~a, ~b], ~x) && (InverseFunctionFreeQ(~u, ~x) && (Not(MatchQ(~v, (~(c') + ~(d') * ~x) ^ ~(m') <-- FreeQ([c, d, m], ~x))) && FalseQ(FunctionOfLinear(~v * (~a + ~b * acoth(~u)), ~x))))
end

max_stars_repo_path: src/rules/7 Inverse hyperbolic functions/7.3 Inverse hyperbolic tangent/7.3.7 Miscellaneous inverse hyperbolic tangent.jl


text:
"""
    FEMMHeatDiffSurfModule

Module for operations on boundaries of domains to construct system matrices and
system vectors for linear heat diffusion/conduction.
"""
module FEMMHeatDiffSurfModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.FESetModule: AbstractFESet, nodesperelem, manifdim
import FinEtools.IntegDomainModule: IntegDomain, integrationdata, Jacobiansurface
import FinEtools.FieldModule: ndofs, gatherdofnums!, gatherfixedvalues_asvec!, gathervalues_asvec!, gathervalues_asmat!
import FinEtools.NodalFieldModule: NodalField
import FinEtools.AssemblyModule: AbstractSysvecAssembler, AbstractSysmatAssembler, SysmatAssemblerSparseSymm, startassembly!, assemble!, makematrix!, makevector!, SysvecAssembler
import FinEtools.FEMMBaseModule: AbstractFEMM
import FinEtools.MatrixUtilityModule: add_gkgt_ut_only!, add_nnt_ut_only!, complete_lt!, locjac!
import LinearAlgebra: norm, dot, cross

"""
    FEMMHeatDiffSurf{S<:AbstractFESet, F<:Function} <: AbstractFEMM

    Type for heat diffusion finite element modeling machine for boundary integrals.
"""
mutable struct FEMMHeatDiffSurf{S<:AbstractFESet, F<:Function} <: AbstractFEMM
    integdomain::IntegDomain{S, F} # geometry data
    surfacetransfercoeff::FFlt # material object
end

"""
    surfacetransfer(self::FEMMHeatDiffSurf,  assembler::A,
      geom::NodalField{FFlt}, temp::NodalField{FFlt}) where {A<:AbstractSysmatAssembler}

Compute the surface heat transfer matrix.

# Arguments
- `self` = model machine,
- `assembler` = matrix assembler
- `geom` = geometry field,
- `temp` = temperature field
"""
function surfacetransfer(self::FEMMHeatDiffSurf,  assembler::A, geom::NodalField{FFlt}, temp::NodalField{FFlt}) where {A<:AbstractSysmatAssembler}
    fes = self.integdomain.fes
    # Constants
    nfes = count(fes); # number of finite elements in the set
    ndn = ndofs(temp); # number of degrees of freedom per node
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    Hedim = ndn*nne;             # dimension of the element matrix
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    # Prepare assembler and temporaries
    ecoords = fill(zero(FFlt), nne, ndofs(geom)); # array of Element coordinates
    He = fill(zero(FFlt), Hedim, Hedim);                # element matrix -- used as a buffer
    dofnums = zeros(FInt, Hedim); # degree of freedom array -- used as a buffer
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- used as a buffer
    J = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- used as a buffer
    startassembly!(assembler, Hedim, Hedim, nfes, temp.nfreedofs, temp.nfreedofs);
    for i = 1:nfes # Loop over elements
        gathervalues_asmat!(geom, ecoords, fes.conn[i]);
        fill!(He,  0.0); # Initialize element matrix
        for j=1:npts # Loop over quadrature points
            locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
            Jac = Jacobiansurface(self.integdomain, J, loc, fes.conn[i],  Ns[j]);
            add_nnt_ut_only!(He, Ns[j], self.surfacetransfercoeff*Jac*w[j])
        end # Loop over quadrature points
        complete_lt!(He)
        gatherdofnums!(temp, dofnums, fes.conn[i]);# retrieve degrees of freedom
        assemble!(assembler, He, dofnums, dofnums);# assemble symmetric matrix
    end # Loop over elements
    return makematrix!(assembler);
end

function surfacetransfer(self::FEMMHeatDiffSurf{S},
                                   geom::NodalField{FFlt},
                                   temp::NodalField{FFlt}) where {S<:AbstractFESet}
    assembler = SysmatAssemblerSparseSymm()
    return  surfacetransfer(self, assembler, geom, temp);
end

"""
    surfacetransferloads(self::FEMMHeatDiffSurf,  assembler::A,
      geom::NodalField{FFlt}, temp::NodalField{FFlt},
      ambtemp::NodalField{FFlt}) where {A<:AbstractSysvecAssembler}

Compute the load vector corresponding to surface heat transfer.

# Arguments
- `self` = model machine,
- `assembler` = matrix assembler
- `geom` = geometry field,
- `temp` = temperature field
- `ambtemp` = ambient temperature field on the surface
"""
function surfacetransferloads(self::FEMMHeatDiffSurf,  assembler::A,  geom::NodalField{FFlt}, temp::NodalField{FFlt},  ambtemp::NodalField{FFlt}) where {A<:AbstractSysvecAssembler}
    fes = self.integdomain.fes
    # Constants
    nfes = count(fes); # number of finite elements in the set
    ndn = ndofs(temp); # number of degrees of freedom per node
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    Hedim = ndn*nne;             # dimension of the element matrix
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    # Prepare assembler and temporaries
    ecoords = fill(zero(FFlt), nne, ndofs(geom)); # array of Element coordinates
    Fe = fill(zero(FFlt), Hedim, 1); # element matrix -- used as a buffer
    dofnums = zeros(FInt, Hedim); # degree of freedom array -- used as a buffer
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- used as a buffer
    J = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- used as a buffer
    pT = fill(zero(FFlt), Hedim);
    startassembly!(assembler,  temp.nfreedofs);
    for i = 1:nfes # Loop over elements
        gathervalues_asvec!(ambtemp, pT, fes.conn[i]);# retrieve ambient temp
        if norm(pT, Inf) != 0.0    # Is the load nonzero?
            gathervalues_asmat!(geom, ecoords, fes.conn[i]);
            fill!(Fe,  0.0); # Initialize element matrix
            for j=1:npts # Loop over quadrature points
                locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
                Jac = Jacobiansurface(self.integdomain, J, loc, fes.conn[i],  Ns[j]);
                Ta = dot(vec(pT), vec(Ns[j]))
                factor = Ta*self.surfacetransfercoeff*Jac*w[j]
                Fe .+= factor*Ns[j]
            end # Loop over quadrature points
            gatherdofnums!(temp, dofnums, fes.conn[i]); # retrieve degrees of freedom
            assemble!(assembler,  Fe,  dofnums); # assemble element load vector
        end
    end
    F = makevector!(assembler);
    return F
end

function surfacetransferloads(self::FEMMHeatDiffSurf,
                                        geom::NodalField{FFlt},
                                        temp::NodalField{FFlt},
                                        ambtemp::NodalField{FFlt})
    assembler = SysvecAssembler()
    return  surfacetransferloads(self, assembler, geom, temp, ambtemp);
end

"""
    nzebcsurfacetransferloads(self::FEMMHeatDiffSurf, assembler::A,
      geom::NodalField{FFlt}, temp::NodalField{FFlt}) where {A<:AbstractSysvecAssembler}

Compute load vector for nonzero EBC for fixed temperature.

# Arguments
- `self` = model machine,
- `assembler` = matrix assembler
- `geom` = geometry field,
- `temp` = temperature field; the `fixed_values` attribute needs to list
      the prescribed values of the temperature. This is performed with
      `applyebc!`.
"""
function nzebcsurfacetransferloads(self::FEMMHeatDiffSurf, assembler::A,  geom::NodalField{FFlt}, temp::NodalField{FFlt}) where {A<:AbstractSysvecAssembler}
    fes = self.integdomain.fes
    # Constants
    nfes = count(fes); # number of finite elements in the set
    ndn = ndofs(temp); # number of degrees of freedom per node
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    Hedim = ndn*nne;             # dimension of the element matrix
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    # Prepare assembler and temporaries
    ecoords = fill(zero(FFlt), nne, ndofs(geom)); # array of Element coordinates
    He = fill(zero(FFlt), Hedim, Hedim);                # element matrix -- used as a buffer
    dofnums = zeros(FInt, Hedim); # degree of freedom array -- used as a buffer
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- used as a buffer
    J = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- used as a buffer
    pT = fill(zero(FFlt), Hedim);
    startassembly!(assembler,  temp.nfreedofs);
    # Now loop over all finite elements in the set
    for i=1:nfes # Loop over elements
        gatherfixedvalues_asvec!(temp, pT, fes.conn[i]);# retrieve element temp
        if norm(pT, Inf) != 0.0    # Is the load nonzero?
            gathervalues_asmat!(geom, ecoords, fes.conn[i]);
            fill!(He,  0.0);
            for j=1:npts # Loop over quadrature points
                locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
                Jac = Jacobiansurface(self.integdomain, J, loc, fes.conn[i],  Ns[j]);
                add_nnt_ut_only!(He, Ns[j], self.surfacetransfercoeff*Jac*w[j])
            end # Loop over quadrature points
            complete_lt!(He)
            gatherdofnums!(temp, dofnums, fes.conn[i]); # retrieve degrees of freedom
            assemble!(assembler,  vec(He*pT),  dofnums); # assemble element load vector
        end
    end
    F= makevector!(assembler);
    return F
end

function nzebcsurfacetransferloads(self::FEMMHeatDiffSurf,
                                          geom::NodalField{FFlt},
                                          temp::NodalField{FFlt})
    assembler = SysvecAssembler()
    return  nzebcsurfacetransferloads(self, assembler, geom, temp);
end

end

max_stars_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff

max_issues_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff

max_forks_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff


text:
@testset "Custom function evaluation" begin

funcs = Dict()
funcs[:pi] = pi
funcs[:zero] = () -> 0
funcs[:sin] = sin
funcs[:cos] = cos
funcs[:square] = x -> x * x
funcs[:dup] = x -> (x, x)
funcs[:pair] = (x, y) -> (x, y)
funcs[:fst] = tup -> tup[1]
funcs[:snd] = tup -> tup[2]
funcs[:fakesum] = Dict((1, 1) => 2, (2, 2) => 4)

@test resolve(@julog(zero == 0), Clause[], funcs=funcs)[1] == true
@test resolve(@julog(sin(pi / 2) == 1), Clause[], funcs=funcs)[1] == true
@test resolve(@julog(cos(pi) == -1), Clause[], funcs=funcs)[1] == true
@test resolve(@julog(square(5) == 25), Clause[], funcs=funcs)[1] == true
@test resolve(@julog(dup(6) == pair(6, 6)), Clause[], funcs=funcs)[1] == true
@test resolve(@julog(fakesum(1, 1) == 2), Clause[], funcs=funcs)[1] == true
@test resolve(@julog(fakesum(2, 2) == 4), Clause[], funcs=funcs)[1] == true

clauses = @julog [
    even(X) <<= mod(X, 2) == 0,
    is_dup(X) <<= fst(X) == snd(X),
    is_square(5, 5, 25) <<= true
]

@test resolve(@julog(even(fakesum(1, 1))), clauses, funcs=funcs)[1] == true
@test resolve(@julog(is_dup(dup(8))), clauses, funcs=funcs)[1] == true
@test resolve(@julog(is_square(5, X, square(5))), clauses, funcs=funcs)[1] == true

clauses = @julog [
    on_circ(Rad, Pt) <<= square(fst(Pt)) + square(snd(Pt)) == square(Rad),
    on_diag(X, Y) <<= dup(X) == pair(X, Y)
]

# Is the point (3, -4) on the circle of radius 5?
@test resolve(@julog(on_circ(5, (3, -4))), clauses, funcs=funcs)[1] == true
# How about the point (5 * sin(1), 5 * cos(1))?
@test resolve(@julog(on_circ(5, pair(5*sin(1), 5*cos(1)))), clauses, funcs=funcs)[1] == true
# Is the point (10, 10) on the line X=Y?
@test resolve(@julog(on_diag(10, 10)), clauses, funcs=funcs)[1] == true

end


text:
include("euler/euler.jl")
using .BigIntegers: BigInteger
using .Numbers: digits_sum

function compute(n::Int)::Int
    current, previous = BigInteger(1), BigInteger(1)
    for i ∈ 2:n + 1
        current, previous = (i % 3 == 1 ? 2i ÷ 3 : 1)current + previous, current
    end
    return digits_sum(current)
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler