text:
# bc.jl

@doc """

  This function calculates the boundary flux for the portion of the boundary
  with a particular boundary condition.  The eqn.q are converted to
  conservative variables if needed.  For the DG version, eqn.q_bndry must
  already be populated with the q variables interpolated to the boundary

  Inputs:

  mesh : AbstractMesh
  sbp : AbstractOperator
  eqn : AdvectionEquation
  functor : a callable object that calculates the boundary flux at a node
  idx_range: the Range describing which Boundaries have the current BC
  bndry_facenums:  An array with elements of type Boundary that tell which
                   element faces have the boundary condition
  Outputs:

  bndryflux : the array to store the boundary flux, corresponds to
              bndry_facenums

  note that bndry_facenums and bndryflux must be only the portion of the
  their parent arrays that correspond to the Boundaries that have the
  current boundary condition applied.

  The functor must have the signature:
  functor( q, aux_vars, x, dxidx, nrm, bndryflux_i, eqn.params)
  where q are the *conservative* variables.
  where all arguments (except params and nrm) are vectors of values at a node.

  params is the ParamType associated with the the EulerEquation object
  nrm = mesh.sbpface.normal[:, current_node]

  This is a mid level function.
"""->
# mid level function
function calcBoundaryFlux( mesh::AbstractCGMesh{Tmsh},
       sbp::AbstractOperator, eqn::AdvectionData{Tsol},
       functor::BCType, idx_range::UnitRange,
       bndry_facenums::AbstractArray{Boundary,1},
       bndryflux::AbstractArray{Tres, 3}) where {Tmsh,  Tsol, Tres}
# calculate the boundary flux for the boundary condition evaluated by the functor

#  println("enterted calcBoundaryFlux CG")
#TODO: update this to use mesh.nrm_bndry

  t = eqn.t
  nfaces = length(bndry_facenums)
  nrm_scaled = zeros(Tmsh, mesh.dim)
  for i=1:nfaces  # loop over faces with this BC
    bndry_i = bndry_facenums[i]
    for j = 1:mesh.numNodesPerFace
      k = mesh.facenodes[j, bndry_i.face]

      # get components
      q = eqn.q[ 1, k, bndry_i.element]
      alpha_x = eqn.params.alpha_x
      alpha_y = eqn.params.alpha_y
      coords = ro_sview(mesh.coords, :, k, bndry_i.element)
      dxidx = ro_sview(mesh.dxidx, :, :, k, bndry_i.element)
      nrm = ro_sview(mesh.sbpface.normal, :, bndry_i.face)
      calcBCNormal(eqn.params, dxidx, nrm, nrm_scaled)
      bndryflux[1, j, i] = -functor(eqn.params, q, coords, nrm_scaled, t)
    end
  end

  return nothing
end

# DG version
function calcBoundaryFlux( mesh::AbstractDGMesh{Tmsh},
       sbp::AbstractOperator, eqn::AdvectionData{Tsol},
       functor::BCType, idx_range::UnitRange,
       bndry_facenums::AbstractArray{Boundary,1},
       bndryflux::AbstractArray{Tres, 3}) where {Tmsh,  Tsol, Tres}

# calculate the boundary flux for the boundary condition evaluated by the functor

#  println("entered calcBoundaryFlux DG")
  t = eqn.t
  nfaces = length(bndry_facenums)
  for i=1:nfaces  # loop over faces with this BC
    bndry_i = bndry_facenums[i]
    global_facenum = idx_range[i]
    for j = 1:mesh.numNodesPerFace

      # get components
      q = eqn.q_bndry[ 1, j, global_facenum]
      alpha_x = eqn.params.alpha_x
      alpha_y = eqn.params.alpha_y
      coords = ro_sview(mesh.coords_bndry, :, j, global_facenum)
      nrm_scaled = ro_sview(mesh.nrm_bndry, :, j, global_facenum)
      bndryflux[1, j, i] = -functor(eqn.params, q, coords, nrm_scaled, t)
    end
  end

  return nothing
end

"""
  This function computes the boundary integrals (and should probably be renamed)
  without using eqn.q_bndry of eqn.bndryflux.  eqn.res is updated with
  the results.

  See calcBoundaryFlux for the meaning of the arguments
"""
function calcBoundaryFlux_nopre( mesh::AbstractDGMesh{Tmsh},
                          sbp::AbstractOperator, eqn::AdvectionData{Tsol},
                          functor::BCType, idx_range::UnitRange,
                          bndry_facenums::AbstractArray{Boundary,1},
                          bndryflux::AbstractArray{Tres, 3}) where {Tmsh,  Tsol, Tres}

  t = eqn.t
  nfaces = length(bndry_facenums)
  q_face = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  flux_face = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

  for i=1:nfaces  # loop over faces with this BC
    bndry_i = bndry_facenums[i]
    global_facenum = idx_range[i]

    # interpolate to face
    q_el = ro_sview(eqn.q, :, :, bndry_i.element)
    boundaryFaceInterpolate!(mesh.sbpface, bndry_i.face, q_el, q_face)

    for j = 1:mesh.numNodesPerFace
      # get components
      q = q_face[ 1, j]
      coords = ro_sview(mesh.coords_bndry, :, j, global_facenum)
      nrm_scaled = ro_sview(mesh.nrm_bndry, :, j, i)
      flux_face[1, j] = -functor(eqn.params, q, coords, nrm_scaled, t)
    end

    res_i = sview(eqn.res, :, :, bndry_i.element)
    boundaryFaceIntegrate!(mesh.sbpface, bndry_i.face, flux_face, res_i)
  end

  return nothing
end

@doc """
### AdvectionEquationMod.x5plusy5BC

Calculates q at the boundary which is equal to x^5 + y^5. It is a nodal
level function.

**Inputs**

*  `u` : Advection variable (eqn.q)
*  `params`: the equation ParamType
*  `coords` : Nodal coordinates
*  `nrm_scaled`    : scaled face normal vector in x-y space
*  `t`:  current time value

**Outputs**
*  `bndryflux` : Flux at the boundary

*  None

"""->

mutable struct x5plusy5BC <: BCType
end

function (obj::x5plusy5BC)(params::ParamType2, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_x5plusy5(params, coords, t) # Calculate the actual analytic value of u at the bondary
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

mutable struct constantBC <: BCType
end

function (obj::constantBC)(params::ParamTypes, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = 2
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)
  return bndryflux
end

@doc """
### AdvectionEquationMod.exp_xplusyBC

Calculates q at the boundary which is equal to exp(x+y). It is a nodal
level function.

**Inputs**

*  `u` : Advection variable (eqn.q)
*  `alpha_x` & `alpha_y` : velocities in the X & Y directions
*  `coords` : Nodal coordinates
*  `dxidx`  : Mapping Jacobian
*  `nrm`    : SBP face-normal vectors
*  `bndryflux` : Flux at the boundary

**Outputs**

*  None

"""->
mutable struct exp_xplusyBC <: BCType
end

function (obj::exp_xplusyBC)(params::ParamType2, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_exp_xplusy(params, coords, t)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

@doc """
### AdvectionEquationMod.sinwave_BC

  Uses the Roe solver to calculate the boundary flux using calc_sinewave to
  get the boundary state
"""->
mutable struct sinwave_BC <: BCType
end

function (obj::sinwave_BC)(params::ParamType2, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_sinwave(params, coords, t)
#  println("  u_bc = ", u_bc)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

@doc """
### AdvectionEquationMod.sinwavey_BC

  Uses the Roe solver to calculate the boundary flux using calc_sinewavey to
  get the boundary state
"""->
mutable struct sinwavey_BC <: BCType
end

function (obj::sinwavey_BC)(params::ParamType2, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_sinwavey(params, coords, t)
#  println("  u_bc = ", u_bc)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

@doc """
### AdvectionEquationMod.sinwavey_pertBC

  Uses the Roe solver to calculate the boundary flux using calc_sinewave_pert to
  get the boundary state
"""->

mutable struct sinwavey_pertBC <: BCType
end

function (obj::sinwavey_pertBC)(params::ParamType2, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_sinwavey_pert(params, coords, t)
#  println("  u_bc = ", u_bc)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

@doc """
### AdvectionEquationMod.sinwave_ampl_BC

  Uses the Roe solver to calculate the boundary flux using calc_sinewave_ampl to
  get the boundary state
"""->
mutable struct sinwave_ampl_BC <: BCType
end

function (obj::sinwave_ampl_BC)(params::ParamType2, u::Tsol,
              coords::AbstractArray{Tmsh,1}, nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_sinwave_ampl(params, coords, t)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

@doc """
### AdvectionEquationMod.mms1BC

  Uses the Roe solver to calculate the boundary flux using calc_mms1 to get
  the boundary state.
"""->
mutable struct mms1BC <: BCType
end

function (obj::mms1BC)(params::ParamType2, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_mms1(params, coords, t)
#  println("  u_bc = ", u_bc)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

@doc """
### AdvectionEquationMod.x4BC

  Uses the Roe solver to calculate the boundary flux using calc_x4 to
  get the boundary state.
"""->
mutable struct x4BC <: BCType
end

function (obj::x4BC)(params::ParamType2, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_x4(params, coords, t)
#  println("  u_bc = ", u_bc)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

@doc """
### AdvectionEquationMod.p0BC

  Uses the Roe solver to calculate the boundary flux using calc_p0 to
  get the boundary state
"""->

mutable struct p0BC <: BCType
end

function (obj::p0BC)(params::ParamTypes, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_p0(params, coords, t)
#  println("  u_bc = ", u_bc)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)
  return bndryflux
end

@doc """
### AdvectionEquationMod.p1BC

  Uses the Roe solver to calculate the boundary flux using calc_p1 to
  get the boundary state
"""->

mutable struct p1BC <: BCType
end

function (obj::p1BC)(params::ParamTypes, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_p1(params, coords, t)
#  println("  u_bc = ", u_bc)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)
  return bndryflux
end

@doc """
### AdvectionEquationMod.p2BC

  Uses the Roe solver to calculate the boundary flux using calc_p2 to
  get the boundary state
"""->
mutable struct p2BC <: BCType
end

function (obj::p2BC)(params::ParamTypes, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_p2(params, coords, t)
#  println("  u_bc = ", u_bc)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

@doc """
### AdvectionEquationMod.p3BC

  Uses the Roe solver to calculate the boundary flux using calc_p3 to
  get the boundary state
"""->
mutable struct p3BC <: BCType
end

function (obj::p3BC)(params::ParamTypes, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_p3(params, coords, t)
#  println("  u_bc = ", u_bc)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

@doc """
### AdvectionEquationMod.p4BC

  Uses the Roe solver to calculate the boundary flux using calc_p4 to
  get the boundary state.
"""->
mutable struct p4BC <: BCType
end

function (obj::p4BC)(params::ParamTypes, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_p4(params, coords, t)
#  println("  u_bc = ", u_bc)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

@doc """
### AdvectionEquationMod.p5BC

  Uses the Roe solver to calculate the boundary flux using calc_p5 to
  get the boundary state.
"""->

mutable struct p5BC <: BCType
end

function (obj::p5BC)(params::ParamTypes, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}
  # this is really slow: use Horner's rule!
  u_bc = calc_p5(params, coords, t)
#  println("calc_p5 @time printed above")
#  println("  u_bc = ", u_bc)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)
#  println("    RoeSolver @time printed above")

  return bndryflux
end

@doc """
### AdvectionEquationMod.exp5xplus4yplus2BC

Uses the Roe solver to calculate the boundary flux using calc_exp5xplus4yplus2
to get the boundary state.

"""->

mutable struct exp5xplus4yplus2BC <: BCType
end

function (obj::exp5xplus4yplus2BC)(params::ParamType2, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_exp5xplus4yplus2(params, coords, t)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

@doc """
### AdvectionEquationMod.exp5xplusyBC

Uses the Roe solver to calculate the boundary flux using calc_exp5xplusy to get
the boundary state.

"""->

mutable struct exp5xplusyBC <:BCType
end

function (obj::exp5xplusyBC)(params::ParamType2, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_exp5xplusy(params, coords, t)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

@doc """
### AdvectionEquationMod.exp3xplusyBC

Uses the Roe solver to calculate the boundary flux using calc_exp3xplusy to get
the boundary state.

"""->

mutable struct exp3xplusyBC <:BCType
end

function (obj::exp3xplusyBC)(params::ParamType2, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_exp3xplusy(params, coords, t)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

@doc """
### AdvectionEquationMod.exp2xplus2yBC

Uses the Roe solver to calculate the boundary flux using calc_exp2xplus2y to get
the boundary state.

"""->

mutable struct exp2xplus2yBC <:BCType
end

function (obj::exp2xplus2yBC)(params::ParamType2,  u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_exp2xplus2y(params, coords, t)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

@doc """
### AdvectionEquationMod.exp_xyBC

Uses the Roe solver to calculate the boundary flux using calc_exp_xy to get the
boundary state

"""->

mutable struct exp_xyBC <: BCType
end

function (obj::exp_xyBC)(params::ParamType2, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_exp_xy(params, coords, t)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

@doc """
### AdvectionEquationMod.xplusyBC

Uses Roe solver to calculate the boundary flux using calc_xplusy to get the
boundary state

"""->

mutable struct xplusyBC <: BCType
end

function (obj::xplusyBC)(params::ParamType2, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_xplusy(params, coords, t)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

"""
  BC for unsteadymms
"""
mutable struct unsteadymmsBC <: BCType
end

function (obj::unsteadymmsBC)(params::ParamType, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_unsteadymms(params, coords, t)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

"""
  BC for unsteadypoly
"""
mutable struct unsteadypolyBC <: BCType
end

function (obj::unsteadypolyBC)(params::ParamType, u::Tsol,
              coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol}

  u_bc = calc_unsteadypoly(params, coords, t)
  bndryflux = RoeSolver(params, u, u_bc, nrm_scaled)

  return bndryflux
end

"""
  Default BC to calculate the boundary face integral (no numerical flux
  functions)
"""
mutable struct defaultBC <: BCType
end

function (obj::defaultBC)(params::ParamType{Tdim},
              u::Tsol, coords::AbstractArray{Tmsh,1},
              nrm_scaled::AbstractArray{Tmsh,1}, t) where {Tmsh, Tsol, Tdim}

  alpha_nrm = params.alpha_x*nrm_scaled[1] + params.alpha_y*nrm_scaled[2]
  if Tdim == 3  # static analysis
    alpha_nrm += params.alpha_z*nrm_scaled[3]
  end

  return alpha_nrm*u
end

@doc """
### AdvectionEquationMod.BCDict

It stores all the possible boundary condition dictionary options. Whenever a
new boundary condition is created, it should get added to BCDict.

"""->
global const BCDict = Dict{String, BCType}(
"constantBC" => constantBC(),
"x5plusy5BC" => x5plusy5BC(),
"exp_xplusyBC" => exp_xplusyBC(),
"sinwaveBC" => sinwave_BC(),
"sinwaveyBC" => sinwavey_BC(),
"sinwavey_pertBC" => sinwavey_pertBC(),
"sinwaveamplBC" => sinwave_ampl_BC(),
"mms1BC" => mms1BC(),
"x4BC" => x4BC(),
"p0BC" => p0BC(),
"p1BC" => p1BC(),
"p2BC" => p2BC(),
"p3BC" => p3BC(),
"p4BC" => p4BC(),
"p5BC" => p5BC(),
"exp5xplus4yplus2BC" => exp5xplus4yplus2BC(),
"exp5xplusyBC" => exp5xplusyBC(),
"exp3xplusyBC" => exp3xplusyBC(),
"exp2xplus2yBC" => exp2xplus2yBC(),
"exp_xyBC" => exp_xyBC(),
"xplusyBC" => xplusyBC(),
"unsteadymmsBC" => unsteadymmsBC(),
"unsteadypolyBC" => unsteadypolyBC(),
"defaultBC" => defaultBC(),
)

@doc """
### AdvectionEquationMod.getBCFunctors

This function uses the opts dictionary to populate mesh.bndry_funcs with
the the functors

This is a high level function.

**Inputs**

*  `mesh` : Abstract mesh type
*  `sbp`  : Summation-by-parts operator
*  `eqn`  : Advection equation object
*  `opts` : Input dictionary options

**Outputs**

*  None

"""->
# use this function to populate access the needed values in BCDict
function getBCFunctors(mesh::AbstractMesh{Tmsh}, sbp::AbstractOperator,
     eqn::AdvectionData{Tsol, Tdim}, opts) where {Tmsh, Tsol, Tdim}
# populate the array mesh.bndry_funcs with the functors for the boundary condition types

  println("Entered getBCFunctors")

  for i=1:mesh.numBC
    key_i = string("BC", i, "_name")
    val = opts[key_i]
#    println("BCDict[val] = ", BCDict[val])
    mesh.bndry_funcs[i] = BCDict[val]
  end

  return nothing
end


text:
#
# Conversion between Graphs and ITensor networks
#

hasuniqueinds(args...; kwargs...) = !isempty(uniqueinds(args...; kwargs...))

function graph_dir(inds)
  dirs = dir.(inds)
  if length(dirs) == 1
    return only(dirs)
  end
  if all(==(dirs[1]), dirs)
    return dirs[1]
  end
  return ITensors.Out
end

# TODO: rename graph, dispatch on QNs to DiGraph
function Graphs.SimpleDiGraph(tn::Vector{ITensor})
  nv = length(tn)
  g = SimpleDiGraph(nv)
  for v1 in 1:nv, v2 in (v1 + 1):nv
    indsᵛ¹ᵛ² = commoninds(tn[v1], tn[v2])
    if !isempty(commoninds(tn[v1], tn[v2]))
      e = v1 => v2
      if graph_dir(indsᵛ¹ᵛ²) == ITensors.In
        e = reverse(e)
      end
      add_edge!(g, e)
    end
  end
  for v in vertices(g)
    if hasuniqueinds(tn[v], tn[all_neighbors(g, v)]...)
      # Add a self-loop
      add_edge!(g, v => v)
    end
  end
  return g
end

# TODO: rename indsgraph, dispatch on QNs to DiGraph
function MetaGraphs.MetaDiGraph(tn::Vector{ITensor})
  sg = SimpleDiGraph(tn)
  mg = MetaDiGraph(sg)
  for e in edges(mg)
    indsₑ = if is_self_loop(e)
      v = src(e)
      # For self edges, the vertex itself is included as
      # a neighbor so we must exclude it.
      uniqueinds(tn[v], tn[setdiff(all_neighbors(mg, v), v)]...)
    else
      commoninds(tn[src(e)], tn[dst(e)])
    end
    set_prop!(mg, e, :inds, indsₑ)
  end
  return mg
end

default_linkspaces() = 1
default_sitespaces() = 1

default(x, x_default) = x
default(x::Nothing, x_default) = x_default

function itensornetwork(
  g::AbstractGraph; linkspaces=default_linkspaces(), sitespaces=nothing
)
  N = nv(g)
  if !isnothing(sitespaces) && !any_self_loops(g)
    g = copy(g)
    for v in vertices(g)
      add_edge!(g, v => v)
    end
  end
  sitespaces = default(sitespaces, default_sitespaces())
  # TODO: Specialize to Index{typeof(linkspaces)}
  inds_network = [Index[] for _ in 1:N]
  for e in edges(g)
    if !is_self_loop(e)
      lₑ = Index(linkspaces; tags="l=$(src(e))↔$(dst(e))")
      push!(inds_network[src(e)], lₑ)
      push!(inds_network[dst(e)], dag(lₑ))
    else
      sₑ = Index(sitespaces; tags="s=$(src(e))")
      push!(inds_network[src(e)], sₑ)
    end
  end
  tn = Vector{ITensor}(undef, N)
  for n in 1:N
    tn[n] = ITensor(inds_network[n])
  end
  return tn
end

sites(g::Tuple{String,<:Tuple}) = g[2]
sites(g::Tuple{String,<:Tuple,<:NamedTuple}) = g[2]
sites(g::Tuple{String,Int}) = g[2]
sites(g::Tuple{String,Vararg{Int}}) = Base.tail(g)
sites(g::Tuple{String,Int,<:NamedTuple}) = g[2]

# Functionality for turning a list of gates into an ITensor
# network.
function circuit_network(gates, s::Vector{<:Index})
  s = copy(s)
  U = ITensor[]
  for g in gates
    push!(U, op(g, s))
    for n in sites(g)
      s[n] = s[n]'
    end
  end
  return U, s
end

ext: jl


text:
### Proximal gradient method
export SparseProxGradParams, fit!

mutable struct SparseProxGradParams<:AbstractParams
    stepsize::Float64 # initial stepsize
    max_iter::Int # maximum number of outer iterations
    inner_iter::Int # how many prox grad steps to take on X before moving on to Y (and vice versa)
    abs_tol::Float64 # stop if objective decrease upon one outer iteration is less than this
    min_stepsize::Float64 # use a decreasing stepsize, stop when reaches min_stepsize
end
function SparseProxGradParams(stepsize::Number=1.0; # initial stepsize
				              max_iter::Int=100, # maximum number of outer iterations
				              inner_iter::Int=1, # how many prox grad steps to take on X before moving on to Y (and vice versa)
				              abs_tol::Float64=0.00001, # stop if objective decrease upon one outer iteration is less than this
				              min_stepsize::Float64=0.01*stepsize) # stop if stepsize gets this small
    stepsize = convert(Float64, stepsize)
    return SparseProxGradParams(stepsize, max_iter, inner_iter, abs_tol, min_stepsize)
end

### FITTING
function fit!(glrm::GLRM, params::SparseProxGradParams;
			  ch::ConvergenceHistory=ConvergenceHistory("SparseProxGradGLRM"),
			  verbose=true,
			  kwargs...)
	println(params)
	### initialization
	A = glrm.A # rename these for easier local access
	losses = glrm.losses
	rx = glrm.rx
	ry = glrm.ry
	# at any time, glrm.X and glrm.Y will be the best model yet found, while
	# X and Y will be the working variables
	X = copy(glrm.X); Y = copy(glrm.Y)
	k = glrm.k

	m,n = size(A)

    # check that we didn't initialize to zero (otherwise we will never move)
    if norm(Y) == 0
    	Y = .1*randn(k,n)
    end

    # step size (will be scaled below to ensure it never exceeds 1/\|g\|_2 or so for any subproblem)
    alpha = params.stepsize
    # stopping criterion: stop when decrease in objective < tol
    tol = params.abs_tol * mapreduce(length,+,glrm.observed_features)

    # alternating updates of X and Y
    if verbose println("Fitting GLRM") end
    update_ch!(ch, 0, objective(glrm; sparse=true))
    t = time()
    steps_in_a_row = 0
    g = zeros(k)

    # cache views
    ve = [view(X,:,e) for e=1:m]
    vf = [view(Y,:,f) for f=1:n]

    for i=1:params.max_iter
# STEP 1: X update
        for inneri=1:params.inner_iter
        for e=1:m # doing this means looping over XY in row-major order, but otherwise we couldn't parallelize over Xᵢs
            rmul!(g, 0)# reset gradient to 0
            # compute gradient of L with respect to Xᵢ as follows:
            # ∇{Xᵢ}L = Σⱼ dLⱼ(XᵢYⱼ)/dXᵢ
            for f in glrm.observed_features[e]
                # but we have no function dLⱼ/dXᵢ, only dLⱼ/d(XᵢYⱼ) aka dLⱼ/du
                # by chain rule, the result is: Σⱼ (dLⱼ(XᵢYⱼ)/du * Yⱼ), where dLⱼ/du is our grad() function
                # our estimate for A[e,f] is given by dot(ve[e],vf[f])
                axpy!(grad(losses[f],dot(ve[e],vf[f]),A[e,f]), vf[f], g)
            end
            # take a proximal gradient step
            l = length(glrm.observed_features[e]) + 1
            rmul!(g, -alpha/l)
            ## gradient step: Xᵢ += -(α/l) * ∇{Xᵢ}L
            axpy!(1,g,ve[e])
            ## prox step: Xᵢ = prox_rx(Xᵢ, α/l)
            prox!(rx[e],ve[e],alpha/l)
        end
    	end
# STEP 2: Y update
        for inneri=1:params.inner_iter
        for f=1:n
            rmul!(g, 0) # reset gradient to 0
            # compute gradient of L with respect to Yⱼ as follows:
            # ∇{Yⱼ}L = Σⱼ dLⱼ(XᵢYⱼ)/dYⱼ
            for e in glrm.observed_examples[f]
                # but we have no function dLⱼ/dYⱼ, only dLⱼ/d(XᵢYⱼ) aka dLⱼ/du
                # by chain rule, the result is: Σⱼ dLⱼ(XᵢYⱼ)/du * Xᵢ, where dLⱼ/du is our grad() function
            	axpy!(grad(losses[f],dot(ve[e],vf[f]),A[e,f]), ve[e], g)
            end
            # take a proximal gradient step
            l = length(glrm.observed_examples[f]) + 1
            rmul!(g, -alpha/l)
            ## gradient step: Yⱼ += -(α/l) * ∇{Yⱼ}L
            axpy!(1,g,vf[f])
            ## prox step: Yⱼ = prox_ryⱼ(Yⱼ, α/l)
            prox!(ry[f],vf[f],alpha/l)
        end
        end
# STEP 3: Check objective
        obj = objective(glrm, X, Y; sparse=true)
        # record the best X and Y yet found
        if obj < ch.objective[end]
            t = time() - t
            update_ch!(ch, t, obj)
            copy!(glrm.X, X); copy!(glrm.Y, Y) # save new best X and Y
            alpha = alpha * 1.05
            steps_in_a_row = max(1, steps_in_a_row+1)
            t = time()
        else
            # if the objective went up, reduce the step size, and undo the step
            alpha = alpha / max(1.5, -steps_in_a_row)
            if verbose println("obj went up to $obj; reducing step size to $alpha") end
            copy!(X, glrm.X); copy!(Y, glrm.Y) # revert back to last X and Y
            steps_in_a_row = min(0, steps_in_a_row-1)
        end
# STEP 4: Check stopping criterion
        if i>10 && (steps_in_a_row > 3 && ch.objective[end-1] - obj < tol) || alpha <= params.min_stepsize
            break
        end
        if verbose && i%10==0
            println("Iteration $i: objective value = $(ch.objective[end])")
        end
    end
    t = time() - t
    update_ch!(ch, t, ch.objective[end])

    return glrm.X, glrm.Y, ch
end


text:
using ChebyshevQuantum
using Test

tol_val = 1e-8

@testset "test harmonic ω=50" begin

    ω = 50.0

    function V(x)
        return 0.5 * ω^2 * x^2
    end

    vals, vects = ChebyshevQuantum.SolveEig.solve(V=V, N=60, dosplit=false)
    #H = ChebyshevQuantum.SolveEig.setup_ham(V, 60, a=-1.0, b= 1.0);
    #vals, vects =  ChebyshevQuantum.SolveEig.solve(H);

    @test isapprox(vals[1],   (0.5+0) * ω)
    @test isapprox(vals[2],   (0.5+1) * ω)
    @test isapprox(vals[3],   (0.5+2) * ω)
    @test isapprox(vals[4],   (0.5+3) * ω)
    @test isapprox(vals[5],   (0.5+4) * ω)

end

@testset "test harmonic ω=1, " begin

    ω = 1.0

    function V(x)
        return 0.5 * ω^2 * x^2
    end

    vals, vects = ChebyshevQuantum.SolveEig.solve(V=V, N=60, a=-10.0, b= 10.0, dosplit=false)

    #    H = ChebyshevQuantum.SolveEig.setup_ham(V, 60, a=-10.0, b= 10.0);
    #    vals, vects =  ChebyshevQuantum.SolveEig.solve(H);

    @test isapprox(vals[1],   (0.5+0) * ω)
    @test isapprox(vals[2],   (0.5+1) * ω)
    @test isapprox(vals[3],   (0.5+2) * ω)
    @test isapprox(vals[4],   (0.5+3) * ω)
    @test isapprox(vals[5],   (0.5+4) * ω)

end


text:
### A Pluto.jl notebook ###
# v0.16.4

using Markdown
using InteractiveUtils

# ╔═╡ 50c6aa0a-38b8-11ec-36c5-7784eaef9e4c
begin
	using Pkg, Revise
	Pkg.activate("../../")
end

# ╔═╡ 69dd34c7-ac62-443c-8908-eba042438f07
function add_gauss_noise_fl!(img, σ)
	for i in eachindex(img)
		img[i] = img[i] + randn(eltype(img)) * eltype(img)(σ)
		# img[i] = img[i] + randn() * σ # this also works
	end
	return img
end

# ╔═╡ bb684353-9d66-4ecd-b2bc-b24aeadcf972
randn() # returns normal distribution noise

# ╔═╡ 9f6dbe12-a0f1-47a1-9b81-0f0b417b43b3
typeof(randn())

# ╔═╡ a6be7e47-6837-487d-9490-3b4689c1a3fc
begin
	arr = ones(Float32, (10, 10))
	add_gauss_noise_fl!(arr, 123213)
end

# ╔═╡ d82c4907-8de1-47f9-9959-2f562b544770
eltype(arr)

# ╔═╡ 6104149f-c39c-45d1-a9fd-50c1c8b0d71e
arr

# ╔═╡ 8092c8a3-2b09-442a-b8c0-6fce4e8e0889

# ╔═╡ 2ab1d4fc-e191-4261-9bb3-d189db987b55

# ╔═╡ Cell order:
# ╠═50c6aa0a-38b8-11ec-36c5-7784eaef9e4c
# ╠═69dd34c7-ac62-443c-8908-eba042438f07
# ╠═d82c4907-8de1-47f9-9959-2f562b544770
# ╠═bb684353-9d66-4ecd-b2bc-b24aeadcf972
# ╠═9f6dbe12-a0f1-47a1-9b81-0f0b417b43b3
# ╠═a6be7e47-6837-487d-9490-3b4689c1a3fc
# ╠═6104149f-c39c-45d1-a9fd-50c1c8b0d71e
# ╠═8092c8a3-2b09-442a-b8c0-6fce4e8e0889
# ╠═2ab1d4fc-e191-4261-9bb3-d189db987b55


text:
# This file is a part of RadiationDetectorSignals.jl, licensed under the MIT License (MIT).

# TODO: Replace by true custom types to avoid type piracy (e.g. with plotting
# recipes) on NamedTuple and TypedTables.Table (even though the types aliased
# here are quite specific, so the risk of a conflict is low).

const WaveformSamples{T<:RealQuantity} = AbstractVector{T}
const TimeAxis{T<:RealQuantity} = AbstractVector{T}

"""
    RDWaveform

Represents a radiation detector signal waveform.

Fields:

* `time`: time axis, typically a range
* `value`: sample values

Use [`ArrayOfRDWaveforms`](@ref) for arrays of `RDWaveform` that have a
compact memory layout.
"""
struct RDWaveform{
    T<:RealQuantity,U<:RealQuantity,
    TV<:TimeAxis{T},UV<:WaveformSamples{U},
}
    time::TV
    value::UV
end

export RDWaveform

RDWaveform{T,U,TV,UV}(wf::RDWaveform) where {T,U,TV,UV} = RDWaveform{T,U,TV,UV}(wf.time, wf.value)
Base.convert(::Type{RDWaveform{T,U,TV,UV}}, wf::RDWaveform) where {T,U,TV,UV} = RDWaveform{T,U,TV,UV}(wf)

# ToDo: function for waveform duration. Use IntervalSets.duration?

"""
    ArrayOfRDWaveforms = StructArray{<:RDWaveform, ...)

A `StructsArrays.StructArray` of [`RDWaveform`](@ref).

By default, uses `ArraysOfArrays.VectorOfVectors` for contiguous memory
layout.
"""
const ArrayOfRDWaveforms{
    T<:RealQuantity,U<:RealQuantity,N,
    VVT <: AbstractVector{<:AbstractVector{T}}, VVU <: AbstractVector{<:AbstractVector{U}}
} = StructArray{
    <:RDWaveform{T,U},
    N,
    NamedTuple{(:time, :value), Tuple{VVT,VVU}}
}

export ArrayOfRDWaveforms

function StructArray{RDWaveform}(
    contents::Tuple{
        AbstractArray{<:AbstractVector{<:RealQuantity},N},
        AbstractArray{<:AbstractVector{<:RealQuantity},N}
    }
) where {N}
    time, value = contents
    VT = eltype(time)
    VU = eltype(value)
    T = eltype(VT)
    U = eltype(VU)
    StructArray{RDWaveform{T,U,VT,VU}}((time, value))
end

StructArray{RDWaveform}(waveforms::AbstractVector{<:RDWaveform}) =
    StructArray{RDWaveform}((map(w -> w.time, waveforms), VectorOfVectors(map(w -> w.value, waveforms))))

Base.convert(::Type{ArrayOfRDWaveforms}, waveforms::AbstractVector{<:RDWaveform}) = StructArray{RDWaveform}(waveforms)
Base.convert(::Type{ArrayOfRDWaveforms}, waveforms::StructArray{<:RDWaveform}) = waveforms

Base.convert(::Type{StructArray{RDWaveform}}, waveforms::AbstractVector{<:RDWaveform}) = StructArray{RDWaveform}(waveforms)
Base.convert(::Type{StructArray{RDWaveform}}, waveforms::StructArray{<:RDWaveform}) = waveforms

# Specialize getindex to properly support ArraysOfArrays, preventing
# conversion to exact element type:
@inline Base.getindex(A::StructArray{<:RDWaveform}, I::Int...) =
    RDWaveform(A.time[I...], A.value[I...])

@inline ArrayOfRDWaveforms(contents) = StructArray{RDWaveform}(contents)

#=
@generated function Base.getindex(x::StructArray{T, N, NamedTuple{names, types}}, I::Int...) where {T, N, names, types}
    args = [:(getfield(cols, $i)[I...]) for i in 1:length(names)]
    return quote
        cols = fieldarrays(x)
        @boundscheck checkbounds(x, I...)
        @inbounds $(Expr(:call, :createinstance, :T, args...))
    end
end
=#

function Base.float(wf::RDWaveform)
    #...
end

#=
function StatsBase.Histogram(waveforms::AbstractVector{<:RDWaveforms})
    samples = ...

    xedge = axes(samples, 1)
    ph = Histogram((xedge, yedge), Float64, :left)
    @inbounds for evtno in axes(samples, 2)
        for i in axes(samples, 1)
            push!(ph, (i, samples[i, evtno]))
        end
    end
    ph.weights .= log10orNaN.(ph.weights)
    ph
end
=#


text:
"""
    loglik(η, x)

Compute the simple (linear) birth and death process log-likelihood where
``η = (λ, μ)^{\\prime}``, `λ` is the birth rate, `μ` is the death rate, and `x`
is the observed sample.

Define `i` as the size of the population at time 0 and `j` as the population
size at time `t`. Let ``α = (μ e^{(λ - μ) t} - μ) / (λ e^{(λ - μ) t} - μ)`` and
``β = (λ e^{(λ - μ) t} - λ) / (λ e^{(λ - μ) t} - μ)``. Transition probability
`p(j | i, t, λ, μ)` is equal to (Bailey, 1964):
``\\sum_{h = 0}^{\\min(i, j)} \\binom{i}{h} \\binom{i + j - h - 1}{i - 1}
α^{i - h} β^{j - h} (1 - α - β)^{h}``.

Suppose to observe the process at time points `t[0], t[1], ..., t[s], ..., t[S]`
with `t[S] <= T`, at which the population size is
`n[0], n[1], ..., n[s], ..., n[S]`. By the Markov property, log-likelihood is
``l(λ, μ | x) = \\log p(n[S] | n[S], T - t[S], λ, μ) +
\\sum_{s = 1}^{S} \\log p(n[s] | n[s - 1], t[s] - t[s - 1], λ, μ)``.
If we observe `M` independent processes with same parameters `(λ, μ)`, then
``l(λ, μ | x[1], ..., x[M]) = \\sum_{m = 1}^{M} l(λ, μ | x[m]).``

If the process is observed continously over time period `[0, T]`, log-likelihood
simplifies to (Darwin, 1956, Equation (24)):
``l(λ, μ | x) = B \\log λ + D \\log μ - (λ + μ) X + \\sum_{s=0}^{S - 1} n[s]``.
`B` and `D` are the total number of births and deaths observed in `[0, T]`
respectively. `X` is defined as ``n[0] (t[1] - t[0]) + n[1] (t[2] - t[1]) + ...
+ n[S - 1] (t[S] - t[S-1]) + n[S] (T - t[S])``.

# References:

Bailey, N. T. J. (1964). The elements of stochastic processes with applications
to the natural sciences. Wiley, New York, NY, USA. ISBN 0-471-04165-3.

Darwin, J. H. (1956). The behaviour of an estimator for a simple birth and death
process. Biometrika, 43(1/2), 23-31. https://doi.org/10.2307/2333575
"""
function loglik(
  η::Vector{F},
  x::ObservationContinuousTime
)::F where {
  F <: AbstractFloat
}
  x.tot_births * log(η[1]) +
  x.tot_deaths * log(η[2]) -
  (η[1] + η[2]) * x.integrated_jump +
  x.sum_log_n
end

function loglik(
  η::Vector{F},
  x::Vector{ObservationContinuousTime}
)::F where {
  F <: AbstractFloat
}
  mapreduce(y -> loglik(η, y), +, x)
end

function loglik(
  η::Vector{F},
  x::ObservationDiscreteTimeEven
) where {
  F <: AbstractFloat
}
  l = zeros(F, x.n)

  for i = 1:x.n
    itr = zip(x.state[1:(end - 1), i], x.state[2:end, i])
    l[i] = mapreduce(y -> trans_prob(y..., x.u, η), +, itr)
  end

  sum(l)
end

function loglik(
  η::Vector{F},
  x::ObservationDiscreteTimeUneven
) where {
  F <: AbstractFloat
}
  itr = zip(x.state[1:(end - 1), i], x.state[2:end, i], x.waiting_time)
  mapreduce(y -> trans_prob(y..., η), +, itr)
end

function loglik(
  η::Vector{F},
  x::Vector{ObservationDiscreteTimeUneven}
)::F where {
  F <: AbstractFloat
}
  mapreduce(y -> loglik(η, y), +, x)
end


text:
@testset "7.2.2 (d x)^m (a+b arccosh(c x))^n" begin
    (a, b, c, d, f, m, n, x, ) = @variables a b c d f m n x

    #= ::Package:: =#

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

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(b*x)^m*acosh(a*x))^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*acosh(c*x)^n=#

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [x^4*acosh(a*x), x, 6, -((8*sqrt(-1 + a*x)*sqrt(1 + a*x))/(75*a^5)) - (4*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x))/(75*a^3) - (x^4*sqrt(-1 + a*x)*sqrt(1 + a*x))/(25*a) + (1/5)*x^5*acosh(a*x)]
    @test_int [x^3*acosh(a*x), x, 5, -((3*x*sqrt(-1 + a*x)*sqrt(1 + a*x))/(32*a^3)) - (x^3*sqrt(-1 + a*x)*sqrt(1 + a*x))/(16*a) - (3*acosh(a*x))/(32*a^4) + (1/4)*x^4*acosh(a*x)]
    @test_int [x^2*acosh(a*x), x, 4, -((2*sqrt(-1 + a*x)*sqrt(1 + a*x))/(9*a^3)) - (x^2*sqrt(-1 + a*x)*sqrt(1 + a*x))/(9*a) + (1/3)*x^3*acosh(a*x)]
    @test_int [x^1*acosh(a*x), x, 3, -((x*sqrt(-1 + a*x)*sqrt(1 + a*x))/(4*a)) - acosh(a*x)/(4*a^2) + (1/2)*x^2*acosh(a*x)]
    @test_int [x^0*acosh(a*x), x, 2, -((sqrt(-1 + a*x)*sqrt(1 + a*x))/a) + x*acosh(a*x)]
    @test_int [acosh(a*x)/x^1, x, 5, (-(1/2))*acosh(a*x)^2 + acosh(a*x)*log(1 + ℯ^(2*acosh(a*x))) + (1/2)*PolyLog(2, -ℯ^(2*acosh(a*x)))]
    @test_int [acosh(a*x)/x^2, x, 3, -(acosh(a*x)/x) + a*atan(sqrt(-1 + a*x)*sqrt(1 + a*x))]
    @test_int [acosh(a*x)/x^3, x, 2, (a*sqrt(-1 + a*x)*sqrt(1 + a*x))/(2*x) - acosh(a*x)/(2*x^2)]
    @test_int [acosh(a*x)/x^4, x, 5, (a*sqrt(-1 + a*x)*sqrt(1 + a*x))/(6*x^2) - acosh(a*x)/(3*x^3) + (1/6)*a^3*atan(sqrt(-1 + a*x)*sqrt(1 + a*x))]
    @test_int [acosh(a*x)/x^5, x, 4, (a*sqrt(-1 + a*x)*sqrt(1 + a*x))/(12*x^3) + (a^3*sqrt(-1 + a*x)*sqrt(1 + a*x))/(6*x) - acosh(a*x)/(4*x^4)]
    @test_int [acosh(a*x)/x^6, x, 7, (a*sqrt(-1 + a*x)*sqrt(1 + a*x))/(20*x^4) + (3*a^3*sqrt(-1 + a*x)*sqrt(1 + a*x))/(40*x^2) - acosh(a*x)/(5*x^5) + (3/40)*a^5*atan(sqrt(-1 + a*x)*sqrt(1 + a*x))]

    @test_int [x^4*acosh(a*x)^2, x, 7, (16*x)/(75*a^4) + (8*x^3)/(225*a^2) + (2*x^5)/125 - (16*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(75*a^5) - (8*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(75*a^3) - (2*x^4*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(25*a) + (x^5*acosh(a*x)^2)/5]
    @test_int [x^3*acosh(a*x)^2, x, 6, (3*x^2)/(32*a^2) + x^4/32 - (3*x*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(16*a^3) - (x^3*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(8*a) - (3*acosh(a*x)^2)/(32*a^4) + (x^4*acosh(a*x)^2)/4]
    @test_int [x^2*acosh(a*x)^2, x, 5, (4*x)/(9*a^2) + (2*x^3)/27 - (4*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(9*a^3) - (2*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(9*a) + (x^3*acosh(a*x)^2)/3]
    @test_int [x^1*acosh(a*x)^2, x, 4, x^2/4 - (x*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(2*a) - acosh(a*x)^2/(4*a^2) + (x^2*acosh(a*x)^2)/2]
    @test_int [x^0*acosh(a*x)^2, x, 3, 2*x - (2*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/a + x*acosh(a*x)^2]
    @test_int [acosh(a*x)^2/x^1, x, 6, -acosh(a*x)^3/3 + acosh(a*x)^2*log(1 + ℯ^(2*acosh(a*x))) + acosh(a*x)*PolyLog(2, -ℯ^(2*acosh(a*x))) - PolyLog(3, -ℯ^(2*acosh(a*x)))/2]
    @test_int [acosh(a*x)^2/x^2, x, 7, -(acosh(a*x)^2/x) + 4*a*acosh(a*x)*atan(ℯ^acosh(a*x)) - (2*im)*a*PolyLog(2, (-im)*ℯ^acosh(a*x)) + (2*im)*a*PolyLog(2, im*ℯ^acosh(a*x))]
    @test_int [acosh(a*x)^2/x^3, x, 3, (a*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/x - acosh(a*x)^2/(2*x^2) - a^2*log(x)]
    @test_int [acosh(a*x)^2/x^4, x, 9, a^2/(3*x) + (a*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(3*x^2) - acosh(a*x)^2/(3*x^3) + (2*a^3*acosh(a*x)*atan(ℯ^acosh(a*x)))/3 - (im/3)*a^3*PolyLog(2, (-im)*ℯ^acosh(a*x)) + (im/3)*a^3*PolyLog(2, im*ℯ^acosh(a*x))]
    @test_int [acosh(a*x)^2/x^5, x, 5, a^2/(12*x^2) + (a*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(6*x^3) + (a^3*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(3*x) - acosh(a*x)^2/(4*x^4) - (a^4*log(x))/3]

    @test_int [x^4*acosh(a*x)^3, x, 16, -((4144*sqrt(-1 + a*x)*sqrt(1 + a*x))/(5625*a^5)) - (272*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x))/(5625*a^3) - (6*x^4*sqrt(-1 + a*x)*sqrt(1 + a*x))/(625*a) + (16*x*acosh(a*x))/(25*a^4) + (8*x^3*acosh(a*x))/(75*a^2) + (6/125)*x^5*acosh(a*x) - (8*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^2)/(25*a^5) - (4*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^2)/(25*a^3) - (3*x^4*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^2)/(25*a) + (1/5)*x^5*acosh(a*x)^3]
    @test_int [x^3*acosh(a*x)^3, x, 12, -((45*x*sqrt(-1 + a*x)*sqrt(1 + a*x))/(256*a^3)) - (3*x^3*sqrt(-1 + a*x)*sqrt(1 + a*x))/(128*a) - (45*acosh(a*x))/(256*a^4) + (9*x^2*acosh(a*x))/(32*a^2) + (3/32)*x^4*acosh(a*x) - (9*x*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^2)/(32*a^3) - (3*x^3*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^2)/(16*a) - (3*acosh(a*x)^3)/(32*a^4) + (1/4)*x^4*acosh(a*x)^3]
    @test_int [x^2*acosh(a*x)^3, x, 9, -((40*sqrt(-1 + a*x)*sqrt(1 + a*x))/(27*a^3)) - (2*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x))/(27*a) + (4*x*acosh(a*x))/(3*a^2) + (2/9)*x^3*acosh(a*x) - (2*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^2)/(3*a^3) - (x^2*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^2)/(3*a) + (1/3)*x^3*acosh(a*x)^3]
    @test_int [x^1*acosh(a*x)^3, x, 6, -((3*x*sqrt(-1 + a*x)*sqrt(1 + a*x))/(8*a)) - (3*acosh(a*x))/(8*a^2) + (3/4)*x^2*acosh(a*x) - (3*x*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^2)/(4*a) - acosh(a*x)^3/(4*a^2) + (1/2)*x^2*acosh(a*x)^3]
    @test_int [x^0*acosh(a*x)^3, x, 4, -((6*sqrt(-1 + a*x)*sqrt(1 + a*x))/a) + 6*x*acosh(a*x) - (3*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^2)/a + x*acosh(a*x)^3]
    @test_int [acosh(a*x)^3/x^1, x, 7, (-(1/4))*acosh(a*x)^4 + acosh(a*x)^3*log(1 + ℯ^(2*acosh(a*x))) + (3/2)*acosh(a*x)^2*PolyLog(2, -ℯ^(2*acosh(a*x))) - (3/2)*acosh(a*x)*PolyLog(3, -ℯ^(2*acosh(a*x))) + (3/4)*PolyLog(4, -ℯ^(2*acosh(a*x)))]
    @test_int [acosh(a*x)^3/x^2, x, 9, -(acosh(a*x)^3/x) + 6*a*acosh(a*x)^2*atan(ℯ^acosh(a*x)) - 6*im*a*acosh(a*x)*PolyLog(2, (-im)*ℯ^acosh(a*x)) + 6*im*a*acosh(a*x)*PolyLog(2, im*ℯ^acosh(a*x)) + 6*im*a*PolyLog(3, (-im)*ℯ^acosh(a*x)) - 6*im*a*PolyLog(3, im*ℯ^acosh(a*x))]
    @test_int [acosh(a*x)^3/x^3, x, 7, (3/2)*a^2*acosh(a*x)^2 + (3*a*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^2)/(2*x) - acosh(a*x)^3/(2*x^2) - 3*a^2*acosh(a*x)*log(1 + ℯ^(2*acosh(a*x))) - (3/2)*a^2*PolyLog(2, -ℯ^(2*acosh(a*x)))]
    @test_int [acosh(a*x)^3/x^4, x, 13, (a^2*acosh(a*x))/x + (a*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^2)/(2*x^2) - acosh(a*x)^3/(3*x^3) + a^3*acosh(a*x)^2*atan(ℯ^acosh(a*x)) - a^3*atan(sqrt(-1 + a*x)*sqrt(1 + a*x)) - im*a^3*acosh(a*x)*PolyLog(2, (-im)*ℯ^acosh(a*x)) + im*a^3*acosh(a*x)*PolyLog(2, im*ℯ^acosh(a*x)) + im*a^3*PolyLog(3, (-im)*ℯ^acosh(a*x)) - im*a^3*PolyLog(3, im*ℯ^acosh(a*x))]
    @test_int [acosh(a*x)^3/x^5, x, 10, -(a^3*sqrt(-1 + a*x)*sqrt(1 + a*x))/(4*x) + (a^2*acosh(a*x))/(4*x^2) + (a^4*acosh(a*x)^2)/2 + (a*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^2)/(4*x^3) + (a^3*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^2)/(2*x) - acosh(a*x)^3/(4*x^4) - a^4*acosh(a*x)*log(1 + ℯ^(2*acosh(a*x))) - (a^4*PolyLog(2, -ℯ^(2*acosh(a*x))))/2]

    @test_int [x^5*acosh(a*x)^4, x, 23, (245*x^2)/(1152*a^4) + (65*x^4)/(3456*a^2) + x^6/324 - (245*x*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(576*a^5) - (65*x^3*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(864*a^3) - (x^5*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(54*a) - (245*acosh(a*x)^2)/(1152*a^6) + (5*x^2*acosh(a*x)^2)/(16*a^4) + (5*x^4*acosh(a*x)^2)/(48*a^2) + (x^6*acosh(a*x)^2)/18 - (5*x*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^3)/(24*a^5) - (5*x^3*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^3)/(36*a^3) - (x^5*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^3)/(9*a) - (5*acosh(a*x)^4)/(96*a^6) + (x^6*acosh(a*x)^4)/6]
    @test_int [x^4*acosh(a*x)^4, x, 19, (16576*x)/(5625*a^4) + (1088*x^3)/(16875*a^2) + (24*x^5)/3125 - (16576*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(5625*a^5) - (1088*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(5625*a^3) - (24*x^4*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(625*a) + (32*x*acosh(a*x)^2)/(25*a^4) + (16*x^3*acosh(a*x)^2)/(75*a^2) + (12/125)*x^5*acosh(a*x)^2 - (32*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^3)/(75*a^5) - (16*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^3)/(75*a^3) - (4*x^4*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^3)/(25*a) + (1/5)*x^5*acosh(a*x)^4]
    @test_int [x^3*acosh(a*x)^4, x, 14, (45*x^2)/(128*a^2) + (3*x^4)/128 - (45*x*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(64*a^3) - (3*x^3*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(32*a) - (45*acosh(a*x)^2)/(128*a^4) + (9*x^2*acosh(a*x)^2)/(16*a^2) + (3/16)*x^4*acosh(a*x)^2 - (3*x*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^3)/(8*a^3) - (x^3*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^3)/(4*a) - (3*acosh(a*x)^4)/(32*a^4) + (1/4)*x^4*acosh(a*x)^4]
    @test_int [x^2*acosh(a*x)^4, x, 11, (160*x)/(27*a^2) + (8*x^3)/81 - (160*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(27*a^3) - (8*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(27*a) + (8*x*acosh(a*x)^2)/(3*a^2) + (4/9)*x^3*acosh(a*x)^2 - (8*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^3)/(9*a^3) - (4*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^3)/(9*a) + (1/3)*x^3*acosh(a*x)^4]
    @test_int [x^1*acosh(a*x)^4, x, 7, (3*x^2)/4 - (3*x*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/(2*a) - (3*acosh(a*x)^2)/(4*a^2) + (3/2)*x^2*acosh(a*x)^2 - (x*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^3)/a - acosh(a*x)^4/(4*a^2) + (1/2)*x^2*acosh(a*x)^4]
    @test_int [x^0*acosh(a*x)^4, x, 5, 24*x - (24*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x))/a + 12*x*acosh(a*x)^2 - (4*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^3)/a + x*acosh(a*x)^4]
    @test_int [acosh(a*x)^4/x^1, x, 8, (-(1/5))*acosh(a*x)^5 + acosh(a*x)^4*log(1 + ℯ^(2*acosh(a*x))) + 2*acosh(a*x)^3*PolyLog(2, -ℯ^(2*acosh(a*x))) - 3*acosh(a*x)^2*PolyLog(3, -ℯ^(2*acosh(a*x))) + 3*acosh(a*x)*PolyLog(4, -ℯ^(2*acosh(a*x))) - (3/2)*PolyLog(5, -ℯ^(2*acosh(a*x)))]
    @test_int [acosh(a*x)^4/x^2, x, 11, -(acosh(a*x)^4/x) + 8*a*acosh(a*x)^3*atan(ℯ^acosh(a*x)) - 12*im*a*acosh(a*x)^2*PolyLog(2, (-im)*ℯ^acosh(a*x)) + 12*im*a*acosh(a*x)^2*PolyLog(2, im*ℯ^acosh(a*x)) + 24*im*a*acosh(a*x)*PolyLog(3, (-im)*ℯ^acosh(a*x)) - 24*im*a*acosh(a*x)*PolyLog(3, im*ℯ^acosh(a*x)) - 24*im*a*PolyLog(4, (-im)*ℯ^acosh(a*x)) + 24*im*a*PolyLog(4, im*ℯ^acosh(a*x))]
    @test_int [acosh(a*x)^4/x^3, x, 8, 2*a^2*acosh(a*x)^3 + (2*a*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^3)/x - acosh(a*x)^4/(2*x^2) - 6*a^2*acosh(a*x)^2*log(1 + ℯ^(2*acosh(a*x))) - 6*a^2*acosh(a*x)*PolyLog(2, -ℯ^(2*acosh(a*x))) + 3*a^2*PolyLog(3, -ℯ^(2*acosh(a*x)))]
    @test_int [acosh(a*x)^4/x^4, x, 19, (2*a^2*acosh(a*x)^2)/x + (2*a*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^3)/(3*x^2) - acosh(a*x)^4/(3*x^3) - 8*a^3*acosh(a*x)*atan(ℯ^acosh(a*x)) + (4/3)*a^3*acosh(a*x)^3*atan(ℯ^acosh(a*x)) + 4*im*a^3*PolyLog(2, (-im)*ℯ^acosh(a*x)) - 2*im*a^3*acosh(a*x)^2*PolyLog(2, (-im)*ℯ^acosh(a*x)) - 4*im*a^3*PolyLog(2, im*ℯ^acosh(a*x)) + 2*im*a^3*acosh(a*x)^2*PolyLog(2, im*ℯ^acosh(a*x)) + 4*im*a^3*acosh(a*x)*PolyLog(3, (-im)*ℯ^acosh(a*x)) - 4*im*a^3*acosh(a*x)*PolyLog(3, im*ℯ^acosh(a*x)) - 4*im*a^3*PolyLog(4, (-im)*ℯ^acosh(a*x)) + 4*im*a^3*PolyLog(4, im*ℯ^acosh(a*x))]

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [x^6/acosh(a*x), x, 7, (5*SinhIntegral(acosh(a*x)))/(64*a^7) + (9*SinhIntegral(3*acosh(a*x)))/(64*a^7) + (5*SinhIntegral(5*acosh(a*x)))/(64*a^7) + SinhIntegral(7*acosh(a*x))/(64*a^7)]
    @test_int [x^5/acosh(a*x), x, 6, (5*SinhIntegral(2*acosh(a*x)))/(32*a^6) + SinhIntegral(4*acosh(a*x))/(8*a^6) + SinhIntegral(6*acosh(a*x))/(32*a^6)]
    @test_int [x^4/acosh(a*x), x, 6, SinhIntegral(acosh(a*x))/(8*a^5) + (3*SinhIntegral(3*acosh(a*x)))/(16*a^5) + SinhIntegral(5*acosh(a*x))/(16*a^5)]
    @test_int [x^3/acosh(a*x), x, 5, SinhIntegral(2*acosh(a*x))/(4*a^4) + SinhIntegral(4*acosh(a*x))/(8*a^4)]
    @test_int [x^2/acosh(a*x), x, 5, SinhIntegral(acosh(a*x))/(4*a^3) + SinhIntegral(3*acosh(a*x))/(4*a^3)]
    @test_int [x^1/acosh(a*x), x, 4, SinhIntegral(2*acosh(a*x))/(2*a^2)]
    @test_int [x^0/acosh(a*x), x, 2, SinhIntegral(acosh(a*x))/a]
    @test_int [1/(x^1*acosh(a*x)), x, 0, Unintegrable(1/(x*acosh(a*x)), x)]
    @test_int [1/(x^2*acosh(a*x)), x, 0, Unintegrable(1/(x^2*acosh(a*x)), x)]

    @test_int [x^4/acosh(a*x)^2, x, 5, -((x^4*sqrt(-1 + a*x)*sqrt(1 + a*x))/(a*acosh(a*x))) + CoshIntegral(acosh(a*x))/(8*a^5) + (9*CoshIntegral(3*acosh(a*x)))/(16*a^5) + (5*CoshIntegral(5*acosh(a*x)))/(16*a^5)]
    @test_int [x^3/acosh(a*x)^2, x, 4, -((x^3*sqrt(-1 + a*x)*sqrt(1 + a*x))/(a*acosh(a*x))) + CoshIntegral(2*acosh(a*x))/(2*a^4) + CoshIntegral(4*acosh(a*x))/(2*a^4)]
    @test_int [x^2/acosh(a*x)^2, x, 4, -((x^2*sqrt(-1 + a*x)*sqrt(1 + a*x))/(a*acosh(a*x))) + CoshIntegral(acosh(a*x))/(4*a^3) + (3*CoshIntegral(3*acosh(a*x)))/(4*a^3)]
    @test_int [x^1/acosh(a*x)^2, x, 2, -((x*sqrt(-1 + a*x)*sqrt(1 + a*x))/(a*acosh(a*x))) + CoshIntegral(2*acosh(a*x))/a^2]
    @test_int [x^0/acosh(a*x)^2, x, 3, -((sqrt(-1 + a*x)*sqrt(1 + a*x))/(a*acosh(a*x))) + CoshIntegral(acosh(a*x))/a]
    @test_int [1/(x^1*acosh(a*x)^2), x, 0, Unintegrable(1/(x*acosh(a*x)^2), x)]
    @test_int [1/(x^2*acosh(a*x)^2), x, 0, Unintegrable(1/(x^2*acosh(a*x)^2), x)]

    @test_int [x^4/acosh(a*x)^3, x, 14, -((x^4*sqrt(-1 + a*x)*sqrt(1 + a*x))/(2*a*acosh(a*x)^2)) + (2*x^3)/(a^2*acosh(a*x)) - (5*x^5)/(2*acosh(a*x)) + SinhIntegral(acosh(a*x))/(16*a^5) + (27*SinhIntegral(3*acosh(a*x)))/(32*a^5) + (25*SinhIntegral(5*acosh(a*x)))/(32*a^5)]
    @test_int [x^3/acosh(a*x)^3, x, 12, -((x^3*sqrt(-1 + a*x)*sqrt(1 + a*x))/(2*a*acosh(a*x)^2)) + (3*x^2)/(2*a^2*acosh(a*x)) - (2*x^4)/acosh(a*x) + SinhIntegral(2*acosh(a*x))/(2*a^4) + SinhIntegral(4*acosh(a*x))/a^4]
    @test_int [x^2/acosh(a*x)^3, x, 10, -((x^2*sqrt(-1 + a*x)*sqrt(1 + a*x))/(2*a*acosh(a*x)^2)) + x/(a^2*acosh(a*x)) - (3*x^3)/(2*acosh(a*x)) + SinhIntegral(acosh(a*x))/(8*a^3) + (9*SinhIntegral(3*acosh(a*x)))/(8*a^3)]
    @test_int [x^1/acosh(a*x)^3, x, 7, -((x*sqrt(-1 + a*x)*sqrt(1 + a*x))/(2*a*acosh(a*x)^2)) + 1/(2*a^2*acosh(a*x)) - x^2/acosh(a*x) + SinhIntegral(2*acosh(a*x))/a^2]
    @test_int [x^0/acosh(a*x)^3, x, 4, -((sqrt(-1 + a*x)*sqrt(1 + a*x))/(2*a*acosh(a*x)^2)) - x/(2*acosh(a*x)) + SinhIntegral(acosh(a*x))/(2*a)]
    @test_int [1/(x^1*acosh(a*x)^3), x, 0, Unintegrable(1/(x*acosh(a*x)^3), x)]
    @test_int [1/(x^2*acosh(a*x)^3), x, 0, Unintegrable(1/(x^2*acosh(a*x)^3), x)]

    @test_int [x^4/acosh(a*x)^4, x, 12, -(x^4*sqrt(-1 + a*x)*sqrt(1 + a*x))/(3*a*acosh(a*x)^3) + (2*x^3)/(3*a^2*acosh(a*x)^2) - (5*x^5)/(6*acosh(a*x)^2) + (2*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x))/(a^3*acosh(a*x)) - (25*x^4*sqrt(-1 + a*x)*sqrt(1 + a*x))/(6*a*acosh(a*x)) + CoshIntegral(acosh(a*x))/(48*a^5) + (27*CoshIntegral(3*acosh(a*x)))/(32*a^5) + (125*CoshIntegral(5*acosh(a*x)))/(96*a^5)]
    @test_int [x^3/acosh(a*x)^4, x, 9, -(x^3*sqrt(-1 + a*x)*sqrt(1 + a*x))/(3*a*acosh(a*x)^3) + x^2/(2*a^2*acosh(a*x)^2) - (2*x^4)/(3*acosh(a*x)^2) + (x*sqrt(-1 + a*x)*sqrt(1 + a*x))/(a^3*acosh(a*x)) - (8*x^3*sqrt(-1 + a*x)*sqrt(1 + a*x))/(3*a*acosh(a*x)) + CoshIntegral(2*acosh(a*x))/(3*a^4) + (4*CoshIntegral(4*acosh(a*x)))/(3*a^4)]
    @test_int [x^2/acosh(a*x)^4, x, 10, -(x^2*sqrt(-1 + a*x)*sqrt(1 + a*x))/(3*a*acosh(a*x)^3) + x/(3*a^2*acosh(a*x)^2) - x^3/(2*acosh(a*x)^2) + (sqrt(-1 + a*x)*sqrt(1 + a*x))/(3*a^3*acosh(a*x)) - (3*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x))/(2*a*acosh(a*x)) + CoshIntegral(acosh(a*x))/(24*a^3) + (9*CoshIntegral(3*acosh(a*x)))/(8*a^3)]
    @test_int [x^1/acosh(a*x)^4, x, 5, -(x*sqrt(-1 + a*x)*sqrt(1 + a*x))/(3*a*acosh(a*x)^3) + 1/(6*a^2*acosh(a*x)^2) - x^2/(3*acosh(a*x)^2) - (2*x*sqrt(-1 + a*x)*sqrt(1 + a*x))/(3*a*acosh(a*x)) + (2*CoshIntegral(2*acosh(a*x)))/(3*a^2)]
    @test_int [x^0/acosh(a*x)^4, x, 5, -(sqrt(-1 + a*x)*sqrt(1 + a*x))/(3*a*acosh(a*x)^3) - x/(6*acosh(a*x)^2) - (sqrt(-1 + a*x)*sqrt(1 + a*x))/(6*a*acosh(a*x)) + CoshIntegral(acosh(a*x))/(6*a)]
    @test_int [1/(x^1*acosh(a*x)^4), x, 0, Unintegrable(1/(x*acosh(a*x)^4), x)]
    @test_int [1/(x^2*acosh(a*x)^4), x, 0, Unintegrable(1/(x^2*acosh(a*x)^4), x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*acosh(c*x)^(n/2)=#

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [x^4*sqrt(acosh(a*x)), x, 19, (1/5)*x^5*sqrt(acosh(a*x)) - (sqrt(pi)*erf(sqrt(acosh(a*x))))/(32*a^5) - (sqrt(pi/3)*erf(sqrt(3)*sqrt(acosh(a*x))))/(64*a^5) - (sqrt(pi/5)*erf(sqrt(5)*sqrt(acosh(a*x))))/(320*a^5) - (sqrt(pi)*erfi(sqrt(acosh(a*x))))/(32*a^5) - (sqrt(pi/3)*erfi(sqrt(3)*sqrt(acosh(a*x))))/(64*a^5) - (sqrt(pi/5)*erfi(sqrt(5)*sqrt(acosh(a*x))))/(320*a^5)]
    @test_int [x^3*sqrt(acosh(a*x)), x, 14, -((3*sqrt(acosh(a*x)))/(32*a^4)) + (1/4)*x^4*sqrt(acosh(a*x)) - (sqrt(pi)*erf(2*sqrt(acosh(a*x))))/(256*a^4) - (sqrt(pi/2)*erf(sqrt(2)*sqrt(acosh(a*x))))/(32*a^4) - (sqrt(pi)*erfi(2*sqrt(acosh(a*x))))/(256*a^4) - (sqrt(pi/2)*erfi(sqrt(2)*sqrt(acosh(a*x))))/(32*a^4)]
    @test_int [x^2*sqrt(acosh(a*x)), x, 14, (1/3)*x^3*sqrt(acosh(a*x)) - (sqrt(pi)*erf(sqrt(acosh(a*x))))/(16*a^3) - (sqrt(pi/3)*erf(sqrt(3)*sqrt(acosh(a*x))))/(48*a^3) - (sqrt(pi)*erfi(sqrt(acosh(a*x))))/(16*a^3) - (sqrt(pi/3)*erfi(sqrt(3)*sqrt(acosh(a*x))))/(48*a^3)]
    @test_int [x*sqrt(acosh(a*x)), x, 9, -(sqrt(acosh(a*x))/(4*a^2)) + (1/2)*x^2*sqrt(acosh(a*x)) - (sqrt(pi/2)*erf(sqrt(2)*sqrt(acosh(a*x))))/(16*a^2) - (sqrt(pi/2)*erfi(sqrt(2)*sqrt(acosh(a*x))))/(16*a^2)]
    @test_int [sqrt(acosh(a*x)), x, 7, x*sqrt(acosh(a*x)) - (sqrt(pi)*erf(sqrt(acosh(a*x))))/(4*a) - (sqrt(pi)*erfi(sqrt(acosh(a*x))))/(4*a)]
    @test_int [sqrt(acosh(a*x))/x, x, 0, Unintegrable(sqrt(acosh(a*x))/x, x)]

    @test_int [x^4*acosh(a*x)^(3/2), x, 41, -((4*sqrt(-1 + a*x)*sqrt(1 + a*x)*sqrt(acosh(a*x)))/(25*a^5)) - (2*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x)*sqrt(acosh(a*x)))/(25*a^3) - (3*x^4*sqrt(-1 + a*x)*sqrt(1 + a*x)*sqrt(acosh(a*x)))/(50*a) + (1/5)*x^5*acosh(a*x)^(3/2) - (3*sqrt(pi)*erf(sqrt(acosh(a*x))))/(64*a^5) - (sqrt(pi/3)*erf(sqrt(3)*sqrt(acosh(a*x))))/(200*a^5) - (3*sqrt(3*pi)*erf(sqrt(3)*sqrt(acosh(a*x))))/(3200*a^5) - (3*sqrt(pi/5)*erf(sqrt(5)*sqrt(acosh(a*x))))/(3200*a^5) + (3*sqrt(pi)*erfi(sqrt(acosh(a*x))))/(64*a^5) + (sqrt(pi/3)*erfi(sqrt(3)*sqrt(acosh(a*x))))/(200*a^5) + (3*sqrt(3*pi)*erfi(sqrt(3)*sqrt(acosh(a*x))))/(3200*a^5) + (3*sqrt(pi/5)*erfi(sqrt(5)*sqrt(acosh(a*x))))/(3200*a^5)]
    @test_int [x^3*acosh(a*x)^(3/2), x, 25, (-9*x*sqrt(-1 + a*x)*sqrt(1 + a*x)*sqrt(acosh(a*x)))/(64*a^3) - (3*x^3*sqrt(-1 + a*x)*sqrt(1 + a*x)*sqrt(acosh(a*x)))/(32*a) - (3*acosh(a*x)^(3/2))/(32*a^4) + (x^4*acosh(a*x)^(3/2))/4 - (3*sqrt(pi)*erf(2*sqrt(acosh(a*x))))/(2048*a^4) - (3*sqrt(pi/2)*erf(sqrt(2)*sqrt(acosh(a*x))))/(128*a^4) + (3*sqrt(pi)*erfi(2*sqrt(acosh(a*x))))/(2048*a^4) + (3*sqrt(pi/2)*erfi(sqrt(2)*sqrt(acosh(a*x))))/(128*a^4)]
    @test_int [x^2*acosh(a*x)^(3/2), x, 22, -(sqrt(-1 + a*x)*sqrt(1 + a*x)*sqrt(acosh(a*x)))/(3*a^3) - (x^2*sqrt(-1 + a*x)*sqrt(1 + a*x)*sqrt(acosh(a*x)))/(6*a) + (x^3*acosh(a*x)^(3/2))/3 - (3*sqrt(pi)*erf(sqrt(acosh(a*x))))/(32*a^3) - (sqrt(pi/3)*erf(sqrt(3)*sqrt(acosh(a*x))))/(96*a^3) + (3*sqrt(pi)*erfi(sqrt(acosh(a*x))))/(32*a^3) + (sqrt(pi/3)*erfi(sqrt(3)*sqrt(acosh(a*x))))/(96*a^3)]
    @test_int [x*acosh(a*x)^(3/2), x, 11, (-3*x*sqrt(-1 + a*x)*sqrt(1 + a*x)*sqrt(acosh(a*x)))/(8*a) - acosh(a*x)^(3/2)/(4*a^2) + (x^2*acosh(a*x)^(3/2))/2 - (3*sqrt(pi/2)*erf(sqrt(2)*sqrt(acosh(a*x))))/(64*a^2) + (3*sqrt(pi/2)*erfi(sqrt(2)*sqrt(acosh(a*x))))/(64*a^2)]
    @test_int [acosh(a*x)^(3/2), x, 8, (-3*sqrt(-1 + a*x)*sqrt(1 + a*x)*sqrt(acosh(a*x)))/(2*a) + x*acosh(a*x)^(3/2) - (3*sqrt(pi)*erf(sqrt(acosh(a*x))))/(8*a) + (3*sqrt(pi)*erfi(sqrt(acosh(a*x))))/(8*a)]
    @test_int [acosh(a*x)^(3/2)/x, x, 0, Unintegrable(acosh(a*x)^(3/2)/x, x)]

    @test_int [x^4*acosh(a*x)^(5/2), x, 44, (2*x*sqrt(acosh(a*x)))/(5*a^4) + (x^3*sqrt(acosh(a*x)))/(15*a^2) + (3/100)*x^5*sqrt(acosh(a*x)) - (4*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^(3/2))/(15*a^5) - (2*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^(3/2))/(15*a^3) - (x^4*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^(3/2))/(10*a) + (1/5)*x^5*acosh(a*x)^(5/2) - (15*sqrt(pi)*erf(sqrt(acosh(a*x))))/(128*a^5) - (sqrt(pi/3)*erf(sqrt(3)*sqrt(acosh(a*x))))/(240*a^5) - (sqrt(3*pi)*erf(sqrt(3)*sqrt(acosh(a*x))))/(1280*a^5) - (3*sqrt(pi/5)*erf(sqrt(5)*sqrt(acosh(a*x))))/(6400*a^5) - (15*sqrt(pi)*erfi(sqrt(acosh(a*x))))/(128*a^5) - (sqrt(pi/3)*erfi(sqrt(3)*sqrt(acosh(a*x))))/(240*a^5) - (sqrt(3*pi)*erfi(sqrt(3)*sqrt(acosh(a*x))))/(1280*a^5) - (3*sqrt(pi/5)*erfi(sqrt(5)*sqrt(acosh(a*x))))/(6400*a^5)]
    @test_int [x^3*acosh(a*x)^(5/2), x, 27, -((225*sqrt(acosh(a*x)))/(2048*a^4)) + (45*x^2*sqrt(acosh(a*x)))/(256*a^2) + (15/256)*x^4*sqrt(acosh(a*x)) - (15*x*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^(3/2))/(64*a^3) - (5*x^3*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^(3/2))/(32*a) - (3*acosh(a*x)^(5/2))/(32*a^4) + (1/4)*x^4*acosh(a*x)^(5/2) - (15*sqrt(pi)*erf(2*sqrt(acosh(a*x))))/(16384*a^4) - (15*sqrt(pi/2)*erf(sqrt(2)*sqrt(acosh(a*x))))/(512*a^4) - (15*sqrt(pi)*erfi(2*sqrt(acosh(a*x))))/(16384*a^4) - (15*sqrt(pi/2)*erfi(sqrt(2)*sqrt(acosh(a*x))))/(512*a^4)]
    @test_int [x^2*acosh(a*x)^(5/2), x, 24, (5*x*sqrt(acosh(a*x)))/(6*a^2) + (5/36)*x^3*sqrt(acosh(a*x)) - (5*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^(3/2))/(9*a^3) - (5*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^(3/2))/(18*a) + (1/3)*x^3*acosh(a*x)^(5/2) - (15*sqrt(pi)*erf(sqrt(acosh(a*x))))/(64*a^3) - (5*sqrt(pi/3)*erf(sqrt(3)*sqrt(acosh(a*x))))/(576*a^3) - (15*sqrt(pi)*erfi(sqrt(acosh(a*x))))/(64*a^3) - (5*sqrt(pi/3)*erfi(sqrt(3)*sqrt(acosh(a*x))))/(576*a^3)]
    @test_int [x*acosh(a*x)^(5/2), x, 12, -((15*sqrt(acosh(a*x)))/(64*a^2)) + (15/32)*x^2*sqrt(acosh(a*x)) - (5*x*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^(3/2))/(8*a) - acosh(a*x)^(5/2)/(4*a^2) + (1/2)*x^2*acosh(a*x)^(5/2) - (15*sqrt(pi/2)*erf(sqrt(2)*sqrt(acosh(a*x))))/(256*a^2) - (15*sqrt(pi/2)*erfi(sqrt(2)*sqrt(acosh(a*x))))/(256*a^2)]
    @test_int [acosh(a*x)^(5/2), x, 9, (15/4)*x*sqrt(acosh(a*x)) - (5*sqrt(-1 + a*x)*sqrt(1 + a*x)*acosh(a*x)^(3/2))/(2*a) + x*acosh(a*x)^(5/2) - (15*sqrt(pi)*erf(sqrt(acosh(a*x))))/(16*a) - (15*sqrt(pi)*erfi(sqrt(acosh(a*x))))/(16*a)]
    @test_int [acosh(a*x)^(5/2)/x, x, 0, Unintegrable(acosh(a*x)^(5/2)/x, x)]

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [x^4/sqrt(acosh(a*x)), x, 18, -((sqrt(pi)*erf(sqrt(acosh(a*x))))/(16*a^5)) - (sqrt(3*pi)*erf(sqrt(3)*sqrt(acosh(a*x))))/(32*a^5) - (sqrt(pi/5)*erf(sqrt(5)*sqrt(acosh(a*x))))/(32*a^5) + (sqrt(pi)*erfi(sqrt(acosh(a*x))))/(16*a^5) + (sqrt(3*pi)*erfi(sqrt(3)*sqrt(acosh(a*x))))/(32*a^5) + (sqrt(pi/5)*erfi(sqrt(5)*sqrt(acosh(a*x))))/(32*a^5)]
    @test_int [x^3/sqrt(acosh(a*x)), x, 13, -((sqrt(pi)*erf(2*sqrt(acosh(a*x))))/(32*a^4)) - (sqrt(pi/2)*erf(sqrt(2)*sqrt(acosh(a*x))))/(8*a^4) + (sqrt(pi)*erfi(2*sqrt(acosh(a*x))))/(32*a^4) + (sqrt(pi/2)*erfi(sqrt(2)*sqrt(acosh(a*x))))/(8*a^4)]
    @test_int [x^2/sqrt(acosh(a*x)), x, 13, -((sqrt(pi)*erf(sqrt(acosh(a*x))))/(8*a^3)) - (sqrt(pi/3)*erf(sqrt(3)*sqrt(acosh(a*x))))/(8*a^3) + (sqrt(pi)*erfi(sqrt(acosh(a*x))))/(8*a^3) + (sqrt(pi/3)*erfi(sqrt(3)*sqrt(acosh(a*x))))/(8*a^3)]
    @test_int [x/sqrt(acosh(a*x)), x, 8, -(sqrt(pi/2)*erf(sqrt(2)*sqrt(acosh(a*x))))/(4*a^2) + (sqrt(pi/2)*erfi(sqrt(2)*sqrt(acosh(a*x))))/(4*a^2)]
    @test_int [1/sqrt(acosh(a*x)), x, 6, -((sqrt(pi)*erf(sqrt(acosh(a*x))))/(2*a)) + (sqrt(pi)*erfi(sqrt(acosh(a*x))))/(2*a)]
    @test_int [1/(x*sqrt(acosh(a*x))), x, 0, Unintegrable(1/(x*sqrt(acosh(a*x))), x)]
    @test_int [1/(x^2*sqrt(acosh(a*x))), x, 0, Unintegrable(1/(x^2*sqrt(acosh(a*x))), x)]

    @test_int [x^4/acosh(a*x)^(3/2), x, 17, -((2*x^4*sqrt(-1 + a*x)*sqrt(1 + a*x))/(a*sqrt(acosh(a*x)))) + (sqrt(pi)*erf(sqrt(acosh(a*x))))/(8*a^5) + (3*sqrt(3*pi)*erf(sqrt(3)*sqrt(acosh(a*x))))/(16*a^5) + (sqrt(5*pi)*erf(sqrt(5)*sqrt(acosh(a*x))))/(16*a^5) + (sqrt(pi)*erfi(sqrt(acosh(a*x))))/(8*a^5) + (3*sqrt(3*pi)*erfi(sqrt(3)*sqrt(acosh(a*x))))/(16*a^5) + (sqrt(5*pi)*erfi(sqrt(5)*sqrt(acosh(a*x))))/(16*a^5)]
    @test_int [x^3/acosh(a*x)^(3/2), x, 12, -((2*x^3*sqrt(-1 + a*x)*sqrt(1 + a*x))/(a*sqrt(acosh(a*x)))) + (sqrt(pi)*erf(2*sqrt(acosh(a*x))))/(4*a^4) + (sqrt(pi/2)*erf(sqrt(2)*sqrt(acosh(a*x))))/(2*a^4) + (sqrt(pi)*erfi(2*sqrt(acosh(a*x))))/(4*a^4) + (sqrt(pi/2)*erfi(sqrt(2)*sqrt(acosh(a*x))))/(2*a^4)]
    @test_int [x^2/acosh(a*x)^(3/2), x, 12, -((2*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x))/(a*sqrt(acosh(a*x)))) + (sqrt(pi)*erf(sqrt(acosh(a*x))))/(4*a^3) + (sqrt(3*pi)*erf(sqrt(3)*sqrt(acosh(a*x))))/(4*a^3) + (sqrt(pi)*erfi(sqrt(acosh(a*x))))/(4*a^3) + (sqrt(3*pi)*erfi(sqrt(3)*sqrt(acosh(a*x))))/(4*a^3)]
    @test_int [x/acosh(a*x)^(3/2), x, 6, -((2*x*sqrt(-1 + a*x)*sqrt(1 + a*x))/(a*sqrt(acosh(a*x)))) + (sqrt(pi/2)*erf(sqrt(2)*sqrt(acosh(a*x))))/a^2 + (sqrt(pi/2)*erfi(sqrt(2)*sqrt(acosh(a*x))))/a^2]
    @test_int [1/acosh(a*x)^(3/2), x, 7, -((2*sqrt(-1 + a*x)*sqrt(1 + a*x))/(a*sqrt(acosh(a*x)))) + (sqrt(pi)*erf(sqrt(acosh(a*x))))/a + (sqrt(pi)*erfi(sqrt(acosh(a*x))))/a]
    @test_int [1/(x*acosh(a*x)^(3/2)), x, 0, Unintegrable(1/(x*acosh(a*x)^(3/2)), x)]

    @test_int [x^4/acosh(a*x)^(5/2), x, 34, -((2*x^4*sqrt(-1 + a*x)*sqrt(1 + a*x))/(3*a*acosh(a*x)^(3/2))) + (16*x^3)/(3*a^2*sqrt(acosh(a*x))) - (20*x^5)/(3*sqrt(acosh(a*x))) - (sqrt(pi)*erf(sqrt(acosh(a*x))))/(12*a^5) - (3*sqrt(3*pi)*erf(sqrt(3)*sqrt(acosh(a*x))))/(8*a^5) - (5*sqrt(5*pi)*erf(sqrt(5)*sqrt(acosh(a*x))))/(24*a^5) + (sqrt(pi)*erfi(sqrt(acosh(a*x))))/(12*a^5) + (3*sqrt(3*pi)*erfi(sqrt(3)*sqrt(acosh(a*x))))/(8*a^5) + (5*sqrt(5*pi)*erfi(sqrt(5)*sqrt(acosh(a*x))))/(24*a^5)]
    @test_int [x^3/acosh(a*x)^(5/2), x, 24, -((2*x^3*sqrt(-1 + a*x)*sqrt(1 + a*x))/(3*a*acosh(a*x)^(3/2))) + (4*x^2)/(a^2*sqrt(acosh(a*x))) - (16*x^4)/(3*sqrt(acosh(a*x))) - (2*sqrt(pi)*erf(2*sqrt(acosh(a*x))))/(3*a^4) - (sqrt(2*pi)*erf(sqrt(2)*sqrt(acosh(a*x))))/(3*a^4) + (2*sqrt(pi)*erfi(2*sqrt(acosh(a*x))))/(3*a^4) + (sqrt(2*pi)*erfi(sqrt(2)*sqrt(acosh(a*x))))/(3*a^4)]
    @test_int [x^2/acosh(a*x)^(5/2), x, 22, -((2*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x))/(3*a*acosh(a*x)^(3/2))) + (8*x)/(3*a^2*sqrt(acosh(a*x))) - (4*x^3)/sqrt(acosh(a*x)) - (sqrt(pi)*erf(sqrt(acosh(a*x))))/(6*a^3) - (sqrt(3*pi)*erf(sqrt(3)*sqrt(acosh(a*x))))/(2*a^3) + (sqrt(pi)*erfi(sqrt(acosh(a*x))))/(6*a^3) + (sqrt(3*pi)*erfi(sqrt(3)*sqrt(acosh(a*x))))/(2*a^3)]
    @test_int [x/acosh(a*x)^(5/2), x, 11, -((2*x*sqrt(-1 + a*x)*sqrt(1 + a*x))/(3*a*acosh(a*x)^(3/2))) + 4/(3*a^2*sqrt(acosh(a*x))) - (8*x^2)/(3*sqrt(acosh(a*x))) - (2*sqrt(2*pi)*erf(sqrt(2)*sqrt(acosh(a*x))))/(3*a^2) + (2*sqrt(2*pi)*erfi(sqrt(2)*sqrt(acosh(a*x))))/(3*a^2)]
    @test_int [1/acosh(a*x)^(5/2), x, 8, -((2*sqrt(-1 + a*x)*sqrt(1 + a*x))/(3*a*acosh(a*x)^(3/2))) - (4*x)/(3*sqrt(acosh(a*x))) - (2*sqrt(pi)*erf(sqrt(acosh(a*x))))/(3*a) + (2*sqrt(pi)*erfi(sqrt(acosh(a*x))))/(3*a)]
    @test_int [1/(x*acosh(a*x)^(5/2)), x, 0, Unintegrable(1/(x*acosh(a*x)^(5/2)), x)]

    @test_int [x^4/acosh(a*x)^(7/2), x, 32, -((2*x^4*sqrt(-1 + a*x)*sqrt(1 + a*x))/(5*a*acosh(a*x)^(5/2))) + (16*x^3)/(15*a^2*acosh(a*x)^(3/2)) - (4*x^5)/(3*acosh(a*x)^(3/2)) + (32*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x))/(5*a^3*sqrt(acosh(a*x))) - (40*x^4*sqrt(-1 + a*x)*sqrt(1 + a*x))/(3*a*sqrt(acosh(a*x))) + (sqrt(pi)*erf(sqrt(acosh(a*x))))/(30*a^5) + (9*sqrt(3*pi)*erf(sqrt(3)*sqrt(acosh(a*x))))/(20*a^5) + (5*sqrt(5*pi)*erf(sqrt(5)*sqrt(acosh(a*x))))/(12*a^5) + (sqrt(pi)*erfi(sqrt(acosh(a*x))))/(30*a^5) + (9*sqrt(3*pi)*erfi(sqrt(3)*sqrt(acosh(a*x))))/(20*a^5) + (5*sqrt(5*pi)*erfi(sqrt(5)*sqrt(acosh(a*x))))/(12*a^5)]
    @test_int [x^3/acosh(a*x)^(7/2), x, 21, -((2*x^3*sqrt(-1 + a*x)*sqrt(1 + a*x))/(5*a*acosh(a*x)^(5/2))) + (4*x^2)/(5*a^2*acosh(a*x)^(3/2)) - (16*x^4)/(15*acosh(a*x)^(3/2)) + (16*x*sqrt(-1 + a*x)*sqrt(1 + a*x))/(5*a^3*sqrt(acosh(a*x))) - (128*x^3*sqrt(-1 + a*x)*sqrt(1 + a*x))/(15*a*sqrt(acosh(a*x))) + (16*sqrt(pi)*erf(2*sqrt(acosh(a*x))))/(15*a^4) + (4*sqrt(2*pi)*erf(sqrt(2)*sqrt(acosh(a*x))))/(15*a^4) + (16*sqrt(pi)*erfi(2*sqrt(acosh(a*x))))/(15*a^4) + (4*sqrt(2*pi)*erfi(sqrt(2)*sqrt(acosh(a*x))))/(15*a^4)]
    @test_int [x^2/acosh(a*x)^(7/2), x, 22, -((2*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x))/(5*a*acosh(a*x)^(5/2))) + (8*x)/(15*a^2*acosh(a*x)^(3/2)) - (4*x^3)/(5*acosh(a*x)^(3/2)) + (16*sqrt(-1 + a*x)*sqrt(1 + a*x))/(15*a^3*sqrt(acosh(a*x))) - (24*x^2*sqrt(-1 + a*x)*sqrt(1 + a*x))/(5*a*sqrt(acosh(a*x))) + (sqrt(pi)*erf(sqrt(acosh(a*x))))/(15*a^3) + (3*sqrt(3*pi)*erf(sqrt(3)*sqrt(acosh(a*x))))/(5*a^3) + (sqrt(pi)*erfi(sqrt(acosh(a*x))))/(15*a^3) + (3*sqrt(3*pi)*erfi(sqrt(3)*sqrt(acosh(a*x))))/(5*a^3)]
    @test_int [x/acosh(a*x)^(7/2), x, 9, -((2*x*sqrt(-1 + a*x)*sqrt(1 + a*x))/(5*a*acosh(a*x)^(5/2))) + 4/(15*a^2*acosh(a*x)^(3/2)) - (8*x^2)/(15*acosh(a*x)^(3/2)) - (32*x*sqrt(-1 + a*x)*sqrt(1 + a*x))/(15*a*sqrt(acosh(a*x))) + (8*sqrt(2*pi)*erf(sqrt(2)*sqrt(acosh(a*x))))/(15*a^2) + (8*sqrt(2*pi)*erfi(sqrt(2)*sqrt(acosh(a*x))))/(15*a^2)]
    @test_int [1/acosh(a*x)^(7/2), x, 9, -((2*sqrt(-1 + a*x)*sqrt(1 + a*x))/(5*a*acosh(a*x)^(5/2))) - (4*x)/(15*acosh(a*x)^(3/2)) - (8*sqrt(-1 + a*x)*sqrt(1 + a*x))/(15*a*sqrt(acosh(a*x))) + (4*sqrt(pi)*erf(sqrt(acosh(a*x))))/(15*a) + (4*sqrt(pi)*erfi(sqrt(acosh(a*x))))/(15*a)]
    @test_int [1/(x*acosh(a*x)^(7/2)), x, 0, Unintegrable(1/(x*acosh(a*x)^(7/2)), x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(b*x)^m*acosh(c*x)^n*with*m*symbolic=#

    @test_int [x^m*acosh(a*x)^4, x, 1, (x^(1 + m)*acosh(a*x)^4)/(1 + m) - (4*a*Unintegrable((x^(1 + m)*acosh(a*x)^3)/(sqrt(-1 + a*x)*sqrt(1 + a*x)), x))/(1 + m)]
    @test_int [x^m*acosh(a*x)^3, x, 1, (x^(1 + m)*acosh(a*x)^3)/(1 + m) - (3*a*Unintegrable((x^(1 + m)*acosh(a*x)^2)/(sqrt(-1 + a*x)*sqrt(1 + a*x)), x))/(1 + m)]
    @test_int [x^m*acosh(a*x)^2, x, 2, (x^(1 + m)*acosh(a*x)^2)/(1 + m) - (2*a*x^(2 + m)*sqrt(1 - a*x)*acosh(a*x)*HypergeometricFunctions._₂F₁(1/2, (2 + m)/2, (4 + m)/2, a^2*x^2))/((2 + 3*m + m^2)*sqrt(-1 + a*x)) - (2*a^2*x^(3 + m)*HypergeometricPFQ([1, 3/2 + m/2, 3/2 + m/2], [2 + m/2, 5/2 + m/2], a^2*x^2))/(6 + 11*m + 6*m^2 + m^3)]
    @test_int [x^m*acosh(a*x)^1, x, 4, (x^(1 + m)*acosh(a*x))/(1 + m) - (a*x^(2 + m)*sqrt(1 - a^2*x^2)*HypergeometricFunctions._₂F₁(1/2, (2 + m)/2, (4 + m)/2, a^2*x^2))/((2 + 3*m + m^2)*sqrt(-1 + a*x)*sqrt(1 + a*x))]
    @test_int [x^m/acosh(a*x)^1, x, 0, Unintegrable(x^m/acosh(a*x), x)]
    @test_int [x^m/acosh(a*x)^2, x, 0, Unintegrable(x^m/acosh(a*x)^2, x)]
    @test_int [x^m/acosh(a*x)^3, x, 0, Unintegrable(x^m/acosh(a*x)^3, x)]

    @test_int [x^m*acosh(a*x)^(3/2), x, 0, Unintegrable(x^m*acosh(a*x)^(3/2), x)]
    @test_int [x^m*sqrt(acosh(a*x)), x, 0, Unintegrable(x^m*sqrt(acosh(a*x)), x)]
    @test_int [x^m/sqrt(acosh(a*x)), x, 0, Unintegrable(x^m/sqrt(acosh(a*x)), x)]
    @test_int [x^m/acosh(a*x)^(3/2), x, 0, Unintegrable(x^m/acosh(a*x)^(3/2), x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(b*x)^m*acosh(c*x)^n*with*n*symbolic=#

    @test_int [(d*x)^m*acosh(a*x)^n, x, 0, Unintegrable((d*x)^m*acosh(a*x)^n, x)]

    @test_int [x^4*acosh(a*x)^n, x, 12, (5^(-1 - n)*acosh(a*x)^n*Gamma(1 + n, -5*acosh(a*x)))/(32*a^5*(-acosh(a*x))^n) + (acosh(a*x)^n*Gamma(1 + n, -3*acosh(a*x)))/(32*3^n*a^5*(-acosh(a*x))^n) + (acosh(a*x)^n*Gamma(1 + n, -acosh(a*x)))/(16*a^5*(-acosh(a*x))^n) + Gamma(1 + n, acosh(a*x))/(16*a^5) + Gamma(1 + n, 3*acosh(a*x))/(32*3^n*a^5) + (5^(-1 - n)*Gamma(1 + n, 5*acosh(a*x)))/(32*a^5)]
    @test_int [x^3*acosh(a*x)^n, x, 9, (acosh(a*x)^n*Gamma(1 + n, -4*acosh(a*x)))/(2^(2*(3 + n))*a^4*(-acosh(a*x))^n) + (2^(-4 - n)*acosh(a*x)^n*Gamma(1 + n, -2*acosh(a*x)))/(a^4*(-acosh(a*x))^n) + (2^(-4 - n)*Gamma(1 + n, 2*acosh(a*x)))/a^4 + Gamma(1 + n, 4*acosh(a*x))/(2^(2*(3 + n))*a^4)]
    @test_int [x^2*acosh(a*x)^n, x, 9, (3^(-1 - n)*acosh(a*x)^n*Gamma(1 + n, -3*acosh(a*x)))/(8*a^3*(-acosh(a*x))^n) + (acosh(a*x)^n*Gamma(1 + n, -acosh(a*x)))/(8*a^3*(-acosh(a*x))^n) + Gamma(1 + n, acosh(a*x))/(8*a^3) + (3^(-1 - n)*Gamma(1 + n, 3*acosh(a*x)))/(8*a^3)]
    @test_int [x*acosh(a*x)^n, x, 6, (2^(-3 - n)*acosh(a*x)^n*Gamma(1 + n, -2*acosh(a*x)))/(a^2*(-acosh(a*x))^n) + (2^(-3 - n)*Gamma(1 + n, 2*acosh(a*x)))/a^2]
    @test_int [acosh(a*x)^n, x, 4, (acosh(a*x)^n*Gamma(1 + n, -acosh(a*x)))/(2*a*(-acosh(a*x))^n) + Gamma(1 + n, acosh(a*x))/(2*a)]
    @test_int [acosh(a*x)^n/x, x, 0, Unintegrable(acosh(a*x)^n/x, x)]

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

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

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [(a + b*acosh(c*x))*x^3, x, 5, -((3*b*x*sqrt(-1 + c*x)*sqrt(1 + c*x))/(32*c^3)) - (b*x^3*sqrt(-1 + c*x)*sqrt(1 + c*x))/(16*c) - (3*b*acosh(c*x))/(32*c^4) + (1/4)*x^4*(a + b*acosh(c*x))]
    @test_int [(a + b*acosh(c*x))*x^2, x, 4, -((2*b*sqrt(-1 + c*x)*sqrt(1 + c*x))/(9*c^3)) - (b*x^2*sqrt(-1 + c*x)*sqrt(1 + c*x))/(9*c) + (1/3)*x^3*(a + b*acosh(c*x))]
    @test_int [(a + b*acosh(c*x))*x^1, x, 3, -((b*x*sqrt(-1 + c*x)*sqrt(1 + c*x))/(4*c)) - (b*acosh(c*x))/(4*c^2) + (1/2)*x^2*(a + b*acosh(c*x))]
    @test_int [(a + b*acosh(c*x))*x^0, x, 3, a*x - (b*sqrt(-1 + c*x)*sqrt(1 + c*x))/c + b*x*acosh(c*x)]
    @test_int [(a + b*acosh(c*x))/x^1, x, 5, (a + b*acosh(c*x))^2/(2*b) + (a + b*acosh(c*x))*log(1 + ℯ^(-2*acosh(c*x))) - (1/2)*b*PolyLog(2, -ℯ^(-2*acosh(c*x)))]
    @test_int [(a + b*acosh(c*x))/x^2, x, 3, -((a + b*acosh(c*x))/x) + b*c*atan(sqrt(-1 + c*x)*sqrt(1 + c*x))]
    @test_int [(a + b*acosh(c*x))/x^3, x, 2, (b*c*sqrt(-1 + c*x)*sqrt(1 + c*x))/(2*x) - (a + b*acosh(c*x))/(2*x^2)]
    @test_int [(a + b*acosh(c*x))/x^4, x, 5, (b*c*sqrt(-1 + c*x)*sqrt(1 + c*x))/(6*x^2) - (a + b*acosh(c*x))/(3*x^3) + (1/6)*b*c^3*atan(sqrt(-1 + c*x)*sqrt(1 + c*x))]
    @test_int [(a + b*acosh(c*x))/x^5, x, 4, (b*c*sqrt(-1 + c*x)*sqrt(1 + c*x))/(12*x^3) + (b*c^3*sqrt(-1 + c*x)*sqrt(1 + c*x))/(6*x) - (a + b*acosh(c*x))/(4*x^4)]

    #= ::Subsubsection:: =#
    #=n<0=#

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

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [x^2*sqrt(a + b*acosh(c*x)), x, 14, (x^3*sqrt(a + b*acosh(c*x)))/3 - (sqrt(b)*ℯ^(a/b)*sqrt(pi)*erf(sqrt(a + b*acosh(c*x))/sqrt(b)))/(16*c^3) - (sqrt(b)*ℯ^((3*a)/b)*sqrt(pi/3)*erf((sqrt(3)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(48*c^3) - (sqrt(b)*sqrt(pi)*erfi(sqrt(a + b*acosh(c*x))/sqrt(b)))/(16*c^3*ℯ^(a/b)) - (sqrt(b)*sqrt(pi/3)*erfi((sqrt(3)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(48*c^3*ℯ^((3*a)/b))]
    @test_int [x*sqrt(a + b*acosh(c*x)), x, 9, -sqrt(a + b*acosh(c*x))/(4*c^2) + (x^2*sqrt(a + b*acosh(c*x)))/2 - (sqrt(b)*ℯ^((2*a)/b)*sqrt(pi/2)*erf((sqrt(2)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(16*c^2) - (sqrt(b)*sqrt(pi/2)*erfi((sqrt(2)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(16*c^2*ℯ^((2*a)/b))]
    @test_int [sqrt(a + b*acosh(c*x)), x, 7, x*sqrt(a + b*acosh(c*x)) - (sqrt(b)*ℯ^(a/b)*sqrt(pi)*erf(sqrt(a + b*acosh(c*x))/sqrt(b)))/(4*c) - (sqrt(b)*sqrt(pi)*erfi(sqrt(a + b*acosh(c*x))/sqrt(b)))/(4*c*ℯ^(a/b))]

    @test_int [x^2*(a + b*acosh(c*x))^(3/2), x, 22, -(b*sqrt(-1 + c*x)*sqrt(1 + c*x)*sqrt(a + b*acosh(c*x)))/(3*c^3) - (b*x^2*sqrt(-1 + c*x)*sqrt(1 + c*x)*sqrt(a + b*acosh(c*x)))/(6*c) + (x^3*(a + b*acosh(c*x))^(3/2))/3 - (3*b^(3/2)*ℯ^(a/b)*sqrt(pi)*erf(sqrt(a + b*acosh(c*x))/sqrt(b)))/(32*c^3) - (b^(3/2)*ℯ^((3*a)/b)*sqrt(pi/3)*erf((sqrt(3)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(96*c^3) + (3*b^(3/2)*sqrt(pi)*erfi(sqrt(a + b*acosh(c*x))/sqrt(b)))/(32*c^3*ℯ^(a/b)) + (b^(3/2)*sqrt(pi/3)*erfi((sqrt(3)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(96*c^3*ℯ^((3*a)/b))]
    @test_int [x*(a + b*acosh(c*x))^(3/2), x, 11, (-3*b*x*sqrt(-1 + c*x)*sqrt(1 + c*x)*sqrt(a + b*acosh(c*x)))/(8*c) - (a + b*acosh(c*x))^(3/2)/(4*c^2) + (x^2*(a + b*acosh(c*x))^(3/2))/2 - (3*b^(3/2)*ℯ^((2*a)/b)*sqrt(pi/2)*erf((sqrt(2)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(64*c^2) + (3*b^(3/2)*sqrt(pi/2)*erfi((sqrt(2)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(64*c^2*ℯ^((2*a)/b))]
    @test_int [(a + b*acosh(c*x))^(3/2), x, 8, (-3*b*sqrt(-1 + c*x)*sqrt(1 + c*x)*sqrt(a + b*acosh(c*x)))/(2*c) + x*(a + b*acosh(c*x))^(3/2) - (3*b^(3/2)*ℯ^(a/b)*sqrt(pi)*erf(sqrt(a + b*acosh(c*x))/sqrt(b)))/(8*c) + (3*b^(3/2)*sqrt(pi)*erfi(sqrt(a + b*acosh(c*x))/sqrt(b)))/(8*c*ℯ^(a/b))]

    @test_int [x^2*(a + b*acosh(c*x))^(5/2), x, 24, (5*b^2*x*sqrt(a + b*acosh(c*x)))/(6*c^2) + (5*b^2*x^3*sqrt(a + b*acosh(c*x)))/36 - (5*b*sqrt(-1 + c*x)*sqrt(1 + c*x)*(a + b*acosh(c*x))^(3/2))/(9*c^3) - (5*b*x^2*sqrt(-1 + c*x)*sqrt(1 + c*x)*(a + b*acosh(c*x))^(3/2))/(18*c) + (x^3*(a + b*acosh(c*x))^(5/2))/3 - (15*b^(5/2)*ℯ^(a/b)*sqrt(pi)*erf(sqrt(a + b*acosh(c*x))/sqrt(b)))/(64*c^3) - (5*b^(5/2)*ℯ^((3*a)/b)*sqrt(pi/3)*erf((sqrt(3)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(576*c^3) - (15*b^(5/2)*sqrt(pi)*erfi(sqrt(a + b*acosh(c*x))/sqrt(b)))/(64*c^3*ℯ^(a/b)) - (5*b^(5/2)*sqrt(pi/3)*erfi((sqrt(3)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(576*c^3*ℯ^((3*a)/b))]
    @test_int [x*(a + b*acosh(c*x))^(5/2), x, 12, (-15*b^2*sqrt(a + b*acosh(c*x)))/(64*c^2) + (15*b^2*x^2*sqrt(a + b*acosh(c*x)))/32 - (5*b*x*sqrt(-1 + c*x)*sqrt(1 + c*x)*(a + b*acosh(c*x))^(3/2))/(8*c) - (a + b*acosh(c*x))^(5/2)/(4*c^2) + (x^2*(a + b*acosh(c*x))^(5/2))/2 - (15*b^(5/2)*ℯ^((2*a)/b)*sqrt(pi/2)*erf((sqrt(2)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(256*c^2) - (15*b^(5/2)*sqrt(pi/2)*erfi((sqrt(2)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(256*c^2*ℯ^((2*a)/b))]
    @test_int [(a + b*acosh(c*x))^(5/2), x, 9, (15*b^2*x*sqrt(a + b*acosh(c*x)))/4 - (5*b*sqrt(-1 + c*x)*sqrt(1 + c*x)*(a + b*acosh(c*x))^(3/2))/(2*c) + x*(a + b*acosh(c*x))^(5/2) - (15*b^(5/2)*ℯ^(a/b)*sqrt(pi)*erf(sqrt(a + b*acosh(c*x))/sqrt(b)))/(16*c) - (15*b^(5/2)*sqrt(pi)*erfi(sqrt(a + b*acosh(c*x))/sqrt(b)))/(16*c*ℯ^(a/b))]

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [x^2/sqrt(a + b*acosh(c*x)), x, 13, -(ℯ^(a/b)*sqrt(pi)*erf(sqrt(a + b*acosh(c*x))/sqrt(b)))/(8*sqrt(b)*c^3) - (ℯ^((3*a)/b)*sqrt(pi/3)*erf((sqrt(3)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(8*sqrt(b)*c^3) + (sqrt(pi)*erfi(sqrt(a + b*acosh(c*x))/sqrt(b)))/(8*sqrt(b)*c^3*ℯ^(a/b)) + (sqrt(pi/3)*erfi((sqrt(3)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(8*sqrt(b)*c^3*ℯ^((3*a)/b))]
    @test_int [x/sqrt(a + b*acosh(c*x)), x, 8, -(ℯ^((2*a)/b)*sqrt(pi/2)*erf((sqrt(2)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(4*sqrt(b)*c^2) + (sqrt(pi/2)*erfi((sqrt(2)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(4*sqrt(b)*c^2*ℯ^((2*a)/b))]
    @test_int [1/sqrt(a + b*acosh(c*x)), x, 6, -(ℯ^(a/b)*sqrt(pi)*erf(sqrt(a + b*acosh(c*x))/sqrt(b)))/(2*sqrt(b)*c) + (sqrt(pi)*erfi(sqrt(a + b*acosh(c*x))/sqrt(b)))/(2*sqrt(b)*c*ℯ^(a/b))]

    @test_int [x^2/(a + b*acosh(c*x))^(3/2), x, 12, (-2*x^2*sqrt(-1 + c*x)*sqrt(1 + c*x))/(b*c*sqrt(a + b*acosh(c*x))) + (ℯ^(a/b)*sqrt(pi)*erf(sqrt(a + b*acosh(c*x))/sqrt(b)))/(4*b^(3/2)*c^3) + (ℯ^((3*a)/b)*sqrt(3*pi)*erf((sqrt(3)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(4*b^(3/2)*c^3) + (sqrt(pi)*erfi(sqrt(a + b*acosh(c*x))/sqrt(b)))/(4*b^(3/2)*c^3*ℯ^(a/b)) + (sqrt(3*pi)*erfi((sqrt(3)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(4*b^(3/2)*c^3*ℯ^((3*a)/b))]
    @test_int [x/(a + b*acosh(c*x))^(3/2), x, 6, (-2*x*sqrt(-1 + c*x)*sqrt(1 + c*x))/(b*c*sqrt(a + b*acosh(c*x))) + (ℯ^((2*a)/b)*sqrt(pi/2)*erf((sqrt(2)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(b^(3/2)*c^2) + (sqrt(pi/2)*erfi((sqrt(2)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(b^(3/2)*c^2*ℯ^((2*a)/b))]
    @test_int [(a + b*acosh(c*x))^(-3/2), x, 7, (-2*sqrt(-1 + c*x)*sqrt(1 + c*x))/(b*c*sqrt(a + b*acosh(c*x))) + (ℯ^(a/b)*sqrt(pi)*erf(sqrt(a + b*acosh(c*x))/sqrt(b)))/(b^(3/2)*c) + (sqrt(pi)*erfi(sqrt(a + b*acosh(c*x))/sqrt(b)))/(b^(3/2)*c*ℯ^(a/b))]

    @test_int [x^2/(a + b*acosh(c*x))^(5/2), x, 22, (-2*x^2*sqrt(-1 + c*x)*sqrt(1 + c*x))/(3*b*c*(a + b*acosh(c*x))^(3/2)) + (8*x)/(3*b^2*c^2*sqrt(a + b*acosh(c*x))) - (4*x^3)/(b^2*sqrt(a + b*acosh(c*x))) - (ℯ^(a/b)*sqrt(pi)*erf(sqrt(a + b*acosh(c*x))/sqrt(b)))/(6*b^(5/2)*c^3) - (ℯ^((3*a)/b)*sqrt(3*pi)*erf((sqrt(3)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(2*b^(5/2)*c^3) + (sqrt(pi)*erfi(sqrt(a + b*acosh(c*x))/sqrt(b)))/(6*b^(5/2)*c^3*ℯ^(a/b)) + (sqrt(3*pi)*erfi((sqrt(3)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(2*b^(5/2)*c^3*ℯ^((3*a)/b))]
    @test_int [x/(a + b*acosh(c*x))^(5/2), x, 11, (-2*x*sqrt(-1 + c*x)*sqrt(1 + c*x))/(3*b*c*(a + b*acosh(c*x))^(3/2)) + 4/(3*b^2*c^2*sqrt(a + b*acosh(c*x))) - (8*x^2)/(3*b^2*sqrt(a + b*acosh(c*x))) - (2*ℯ^((2*a)/b)*sqrt(2*pi)*erf((sqrt(2)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(3*b^(5/2)*c^2) + (2*sqrt(2*pi)*erfi((sqrt(2)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(3*b^(5/2)*c^2*ℯ^((2*a)/b))]
    @test_int [(a + b*acosh(c*x))^(-5/2), x, 8, (-2*sqrt(-1 + c*x)*sqrt(1 + c*x))/(3*b*c*(a + b*acosh(c*x))^(3/2)) - (4*x)/(3*b^2*sqrt(a + b*acosh(c*x))) - (2*ℯ^(a/b)*sqrt(pi)*erf(sqrt(a + b*acosh(c*x))/sqrt(b)))/(3*b^(5/2)*c) + (2*sqrt(pi)*erfi(sqrt(a + b*acosh(c*x))/sqrt(b)))/(3*b^(5/2)*c*ℯ^(a/b))]

    @test_int [x^2/(a + b*acosh(c*x))^(7/2), x, 22, (-2*x^2*sqrt(-1 + c*x)*sqrt(1 + c*x))/(5*b*c*(a + b*acosh(c*x))^(5/2)) + (8*x)/(15*b^2*c^2*(a + b*acosh(c*x))^(3/2)) - (4*x^3)/(5*b^2*(a + b*acosh(c*x))^(3/2)) + (16*sqrt(-1 + c*x)*sqrt(1 + c*x))/(15*b^3*c^3*sqrt(a + b*acosh(c*x))) - (24*x^2*sqrt(-1 + c*x)*sqrt(1 + c*x))/(5*b^3*c*sqrt(a + b*acosh(c*x))) + (ℯ^(a/b)*sqrt(pi)*erf(sqrt(a + b*acosh(c*x))/sqrt(b)))/(15*b^(7/2)*c^3) + (3*ℯ^((3*a)/b)*sqrt(3*pi)*erf((sqrt(3)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(5*b^(7/2)*c^3) + (sqrt(pi)*erfi(sqrt(a + b*acosh(c*x))/sqrt(b)))/(15*b^(7/2)*c^3*ℯ^(a/b)) + (3*sqrt(3*pi)*erfi((sqrt(3)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(5*b^(7/2)*c^3*ℯ^((3*a)/b))]
    @test_int [x/(a + b*acosh(c*x))^(7/2), x, 9, (-2*x*sqrt(-1 + c*x)*sqrt(1 + c*x))/(5*b*c*(a + b*acosh(c*x))^(5/2)) + 4/(15*b^2*c^2*(a + b*acosh(c*x))^(3/2)) - (8*x^2)/(15*b^2*(a + b*acosh(c*x))^(3/2)) - (32*x*sqrt(-1 + c*x)*sqrt(1 + c*x))/(15*b^3*c*sqrt(a + b*acosh(c*x))) + (8*ℯ^((2*a)/b)*sqrt(2*pi)*erf((sqrt(2)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(15*b^(7/2)*c^2) + (8*sqrt(2*pi)*erfi((sqrt(2)*sqrt(a + b*acosh(c*x)))/sqrt(b)))/(15*b^(7/2)*c^2*ℯ^((2*a)/b))]
    @test_int [(a + b*acosh(c*x))^(-7/2), x, 9, (-2*sqrt(-1 + c*x)*sqrt(1 + c*x))/(5*b*c*(a + b*acosh(c*x))^(5/2)) - (4*x)/(15*b^2*(a + b*acosh(c*x))^(3/2)) - (8*sqrt(-1 + c*x)*sqrt(1 + c*x))/(15*b^3*c*sqrt(a + b*acosh(c*x))) + (4*ℯ^(a/b)*sqrt(pi)*erf(sqrt(a + b*acosh(c*x))/sqrt(b)))/(15*b^(7/2)*c) + (4*sqrt(pi)*erfi(sqrt(a + b*acosh(c*x))/sqrt(b)))/(15*b^(7/2)*c*ℯ^(a/b))]

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

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [sqrt(f*x)*(a + b*acosh(c*x))^2, x, 2, (2*(f*x)^(3/2)*(a + b*acosh(c*x))^2)/(3*f) - (8*b*c*(f*x)^(5/2)*sqrt(1 - c*x)*(a + b*acosh(c*x))*HypergeometricFunctions._₂F₁(1/2, 5/4, 9/4, c^2*x^2))/(15*f^2*sqrt(-1 + c*x)) - (16*b^2*c^2*(f*x)^(7/2)*HypergeometricPFQ([1, 7/4, 7/4], [9/4, 11/4], c^2*x^2))/(105*f^3)]

    #= ::Subsubsection:: =#
    #=n<0=#

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

    @test_int [(d*x)^m*(a + b*acosh(c*x))^2, x, 2, If(13>=8, ((d*x)^(1 + m)*(a + b*acosh(c*x))^2)/(d*(1 + m)) - (2*b*c*(d*x)^(2 + m)*sqrt(1 - c*x)*(a + b*acosh(c*x))*HypergeometricFunctions._₂F₁(1/2, (2 + m)/2, (4 + m)/2, c^2*x^2))/(d^2*(1 + m)*(2 + m)*sqrt(-1 + c*x)) - (2*b^2*c^2*(d*x)^(3 + m)*HypergeometricPFQ([1, 3/2 + m/2, 3/2 + m/2], [2 + m/2, 5/2 + m/2], c^2*x^2))/(d^3*(1 + m)*(2 + m)*(3 + m)), ((d*x)^(1 + m)*(a + b*acosh(c*x))^2)/(d*(1 + m)) - (2*b*c*(d*x)^(2 + m)*sqrt(1 - c*x)*(a + b*acosh(c*x))*HypergeometricFunctions._₂F₁(1/2, (2 + m)/2, (4 + m)/2, c^2*x^2))/(d^2*(1 + m)*(2 + m)*sqrt(-1 + c*x)) - (2*b^2*c^2*(d*x)^(3 + m)*HypergeometricPFQ([1, 3/2 + m/2, 3/2 + m/2], [2 + m/2, 5/2 + m/2], c^2*x^2))/(d^3*(3 + m)*(2 + 3*m + m^2)))]
    @test_int [(d*x)^m*(a + b*acosh(c*x))^1, x, 4, ((d*x)^(1 + m)*(a + b*acosh(c*x)))/(d*(1 + m)) - (b*c*(d*x)^(2 + m)*sqrt(1 - c^2*x^2)*HypergeometricFunctions._₂F₁(1/2, (2 + m)/2, (4 + m)/2, c^2*x^2))/(d^2*(1 + m)*(2 + m)*sqrt(-1 + c*x)*sqrt(1 + c*x))]
    @test_int [(d*x)^m/(a + b*acosh(c*x))^1, x, 0, Unintegrable((d*x)^m/(a + b*acosh(c*x)), x)]
    end

max_stars_repo_path: test/7 Inverse hyperbolic functions/7.2 Inverse hyperbolic cosine/7.2.2 (d x)^m (a+b arccosh(c x))^n.jl


text:
#################################################################
# Randomized estimators of elementary linear algebraic quantities
#################################################################

export rcond, reigmax, reigmin, rnorm, rnorms

import Base: *

"""
    randnn(el, m)
    randnn(el, m, n)

Compute randomized gaussian matrix normalized by column.

# Arguments

- `el::Type`: element type.
- `m::Int`: number of rows.
- `n::Int`: number of columns.

## Keywords

- `normalize::Bool = true`: normalize output.

# Output

**Without `n`**

- `Ω`: vector containing Gaussian random numbers of type `el`.

**With `n`**

- `Ω`: matrix of dimensions `m` x `n` containing Gaussian random numbers of
  type `el`.
"""
function randnn(el, m::Int, normalize::Bool=true)
    if el <: Real
        Ω = randn(m)
    elseif el <: Complex
        Ω = randn(m) + im*randn(m)
    else
        throw(ArgumentError("Unsupported element type: $el"))
    end
    normalize ? Ω/norm(Ω) : Ω
end
function randnn(el, m::Int, n::Int, normalize::Bool=true)
    if el <: Real
        Ω = randn(m, n)
    elseif el <: Complex
        Ω = randn(m, n) + im*randn(m, n)
    else
        throw(ArgumentError("Unsupported element type: $el"))
    end
    normalize || return Ω
    for i=1:n
        Ω[:, i] /= norm(view(Ω, :, i))
    end
    Ω
end

"""
    rnorm(A, mvps)

Compute a probabilistic upper bound on the norm of a matrix `A`.
`‖A‖ ≤ α √(2/π) maxᵢ ‖Aωᵢ‖` with probability `p=α^(-mvps)`.

# Arguments

- `A`: matrix whose norm to estimate.
- `mvps::Int`: number of matrix-vector products to compute.

## Keywords

- `p::Real=0.05`: probability of upper bound failing.

# Output

Estimate of ‖A‖.

See also [`rnorms`](@ref) for a different estimator that uses
premultiplying by both `A` and `A'`.

# References

Lemma 4.1 of Halko2011
"""
function rnorm(A, r::Int, p::Real=0.05)
    @assert 0<p≤1
    α = p^(-1.0/r)
    m, n = size(A)
    Ω = randnn(eltype(A), n, r, false)
    AΩ = A*Ω
    mx = maximum([norm(view(AΩ, :, j)) for j=1:r])
    α * √(2/π) * mx
end

"""
    rnorms(A, iters=1)

Estimate matrix norm randomly using `A'A`.

Compute a probabilistic upper bound on the norm of a matrix `A`.

	ρ = √(‖(A'A)ʲω‖/‖(A'A)ʲ⁻¹ω‖)

which is an estimate of the spectral norm of `A` produced by `iters`
steps of the power method starting with normalized `ω`, is a lower
bound on the true norm by a factor

	ρ ≤ α ‖A‖

with probability greater than `1 - p`, where `p = 4\\sqrt(n/(iters-1)) α^(-2iters)`.

# Arguments

- `A`: matrix whose norm to estimate.
- `iters::Int = 1`: mumber of power iterations to perform.

## Keywords

- `p::Real = 0.05`: probability of upper bound failing.
- `At = A'`: Transpose of `A`.

# Output

Estimate of ‖A‖.

See also [`rnorm`](@ref) for a different estimator that does not require
premultiplying by `A'`

# References

Appendix of [^Liberty2007].
"""
function rnorms(A, j::Int=1, p::Real=0.05; At = A')
    @assert 0<p≤1
     m, n = size(A)
    α = ((j-1)/n*(p/4)^2)^(-1/(4j))

    Ωold = Ω = randnn(eltype(A), n)
    for i=1:j #Power iterations
        Ω, Ωold = At*(A*Ω), Ω
    end
    ρ = √(norm(Ω)/norm(Ωold))
    α*ρ
end

"""
    rcond(A, iters=1)

Estimate matrix condition number randomly.

# Arguments

- `A`: matrix whose condition number to estimate. Must be square and
support premultiply (`A*⋅`) and solve (`A\\⋅`).
- `iters::Int = 1`: number of power iterations to run.

## Keywords

- `p::Real = 0.05`: probability that estimate fails to hold as an upper bound.

# Output

Interval `(x, y)` which contains `κ(A)` with probability `1 - p`.

# Implementation note

[^Dixon1983] originally describes this as a computation that
can be done by computing the necessary number of power iterations given p
and the desired accuracy parameter `θ=y/x`. However, these bounds were only
derived under the assumptions of exact arithmetic. Empirically, `iters≥4` has
been seen to result in incorrect results in that the computed interval does
not contain the true condition number. This implemention therefore makes `iters`
an explicitly user-controllable parameter from which to infer the accuracy
parameter and hence the interval containing `κ(A)`.
```
"""
function rcond(A, k::Int=1, p::Real=0.05)
    @assert 0<p≤1
    m, n = size(A)
    @assert m==n
    θ = (8n/(π*p^2))^(1/k)
    x = randnn(eltype(A), n)
    for i=1:k
        x = A*x
    end
    y = randnn(eltype(A), n)
    for i=1:k
        y = A\y
    end
    φ = ((x⋅x)*(y⋅y))^(1/(2k))
    (φ, θ*φ)
end

"""
    reigmax(A, iters=1)

Estimate maximal eigenvalue randomly.

# Arguments

- `A`: Matrix whose maximal eigenvalue to estimate.
Must be square and support premultiply (`A*⋅`).
- `iters::Int=1`: Number of power iterations to run. (Recommended: `iters` ≤ 3)

## Keywords

- `p::Real=0.05`: Probability that estimate fails to hold as an upper bound.

# Output

Interval `(x, y)` which contains the maximal eigenvalue of `A` with
probability `1 - p`.

# References

Corollary of Theorem 1 in [^Dixon1983]
"""
function reigmax(A, k::Int=1, p::Real=0.05)
    @assert 0<p≤1
    m, n = size(A)
    @assert m==n
    θ = (2n/(π*p^2))^(1/k)
    y = x = randnn(eltype(A), n)
    for i=1:k
        x = A*x
    end
    φ = y⋅x
    (φ, θ*φ)
end

"""
    reigmin(A, iters=1)

Estimate minimal eigenvalue randomly.

# Arguments

- `A`: Matrix whose maximal eigenvalue to estimate.
Must be square and support premultiply (`A*⋅`).
- `iters::Int=1`: Number of power iterations to run. (Recommended: `iters` ≤ 3)

## Keywords

- `p::Real=0.05`: Probability that estimate fails to hold as an upper bound.

# Output

Interval `(x, y)` which contains the maximal eigenvalue of `A` with
probability `1 - p`.

# References

Corollary of Theorem 1 in [^Dixon1983]
"""
function reigmin(A, k::Int=1, p::Real=0.05)
    @assert 0<p≤1
    m, n = size(A)
    @assert m==n
    θ = (2n/(π*p^2))^(1/k)
    y = x = randnn(eltype(A), n)
    for i=1:k
        x = A\x
    end
    φ = y⋅x
    (inv(θ*φ), inv(φ))
end

"""
    srft

A subsampled random Fourier transform.

# Fields

`l<:Integer`: number of vectors to return.

# Implements

`Base`: `*`
"""
struct srft{T<:Integer}
    l :: T
end

#(*)(A::LinearMap, B::srft) = error("method only defined to avoid ambiguity. If you need this method please open a pull request")

"""
    *

Apply a subsampled random Fourier transform to the columns of `A`.

# Arguments

- `A`: matrix to transform.
- `Ω::srft`: subsampled random Fourier transform.

# Output

- `B`: A matrix of dimensions size(A,1) x Ω.l

# References

Equation 4.6 of [^Halko2011].
"""
#Define two methods here to avoid method ambiguity with f::Function*b::Any
*(A::Function, Ω::srft) = function *(A, Ω::srft)
    m, n = size(A)
    B = A*Diagonal(exp(2π*im*rand(n))/√Ω.l)
    B = vcat([fft(A[i,:]) for i=1:m]...) #Factor of √n cancels out
    B[:, randperm(n)[1:Ω.l]]
end


text:
"""
    MOSVGP(args...; kwargs...)

Multi-Output Sparse Variational Gaussian Process

## Arguments
- `kernel::Union{Kernel,AbstractVector{<:Kernel}` : covariance function or vector of covariance functions, can be either a single kernel or a collection of kernels for multiclass and multi-outputs models
- `likelihoods::Union{AbstractLikelihood,Vector{<:Likelihood}` : Likelihood or vector of likelihoods of the model. For compatibilities, see [`Likelihood Types`](@ref likelihood_user)
- `inference` : Inference for the model, for compatibilities see the [`Compatibility Table`](@ref compat_table))
- `nLatent::Int` : Number of latent GPs
- `nInducingPoints` : number of inducing points, or collection of inducing points locations

## Keyword arguments
- `verbose::Int` : How much does the model print (0:nothing, 1:very basic, 2:medium, 3:everything)
- `optimiser` : Optimiser used for the kernel parameters. Should be an Optimiser object from the [Optimisers.jl](https://github.com/FluxML/Optimisers.jl) library. Default is `ADAM(0.001)`
- `Zoptimiser` : Optimiser used for the inducing points locations
- `Aoptimiser` : Optimiser used for the mixing parameters.
- `atfrequency::Int=1` : Choose how many variational parameters iterations are between hyperparameters optimization
- `mean=ZeroMean()` : PriorMean object, check the documentation on it [`MeanPrior`](@ref meanprior)
- `obsdim::Int=1` : Dimension of the data. 1 : X ∈ DxN, 2: X ∈ NxD
"""
mutable struct MOSVGP{
    T<:Real,
    TLikelihood,
    TInference<:AbstractInference,
    N, # Number of tasks
    Q, # Number of latent GPs
} <: AbstractGPModel{T,AbstractLikelihood,TInference,N}
    nf_per_task::NTuple{N,Int}
    f::NTuple{Q,SparseVarLatent}
    likelihood::TLikelihood
    inference::TInference
    A::Vector{Vector{Vector{T}}}
    A_opt::Any
    verbose::Int
    atfrequency::Int
    trained::Bool
end

function MOSVGP(
    kernel::Union{Kernel,AbstractVector{<:Kernel}},
    likelihoods::Union{
        AbstractVector{<:AbstractLikelihood},Tuple{Vararg{<:AbstractLikelihood}}
    },
    inference::AbstractInference,
    Zs::AbstractVector;
    verbose::Int=0,
    atfrequency::Int=1,
    mean::Union{<:Real,AbstractVector{<:Real},PriorMean}=ZeroMean(),
    optimiser=ADAM(0.01),
    Aoptimiser=ADAM(0.01),
    Zoptimiser=false,
    T::DataType=Float64,
)
    likelihoods = likelihoods isa AbstractVector ? tuple(likelihoods...) : likelihoods

    n_task = length(likelihoods)
    nf_per_task = n_latent.(likelihoods)

    inference isa AnalyticVI || error("The inference object should be of type `AnalyticVI`")
    all(implemented.(likelihoods, Ref(inference))) || error(
        "One (or more) of the likelihoods $likelihoods are not compatible or implemented with the $inference",
    )

    if mean isa Real
        mean = ConstantMean(mean)
    elseif mean isa AbstractVector{<:Real}
        mean = EmpiricalMean(mean)
    end

    if isa(optimiser, Bool)
        optimiser = optimiser ? ADAM(0.01) : nothing
    end

    if isa(Zoptimiser, Bool)
        Zoptimiser = Zoptimiser ? ADAM(0.001) : nothing
    end

    if isa(Aoptimiser, Bool)
        Aoptimiser = Aoptimiser ? ADAM(0.01) : nothing
    end

    kernel = if kernel isa Kernel
        (kernel,)
    else
        length(kernel) == n_task ||
            error("Number of kernels should be equal to the number of tasks")
        kernel
    end

    n_kernel = length(kernel)

    num_latent = length(Zs)

    latent_f = ntuple(num_latent) do i
        SparseVarLatent(T, Zs[i], kernel[mod1(i, n_kernel)], mean, optimiser, Zoptimiser)
    end

    function normalize(x)
        return x / sqrt(sum(abs2, x))
    end
    A = [[normalize(randn(T, num_latent)) for i in 1:nf_per_task[j]] for j in 1:n_task]

    return MOSVGP{T,typeof(likelihoods),typeof(inference),n_task,num_latent}(
        nf_per_task,
        latent_f,
        likelihoods,
        inference,
        A,
        Aoptimiser,
        verbose,
        atfrequency,
        false,
    )
end

function Base.show(io::IO, model::MOSVGP)
    return print(
        io,
        "Multioutput Sparse Variational Gaussian Process with the likelihoods $(likelihood(model)) infered by $(inference(model)) ",
    )
end

@traitimpl IsMultiOutput{MOSVGP}

n_output(::MOSVGP{T,L,I,N,Q}) where {T,L,I,N,Q} = Q
Zviews(m::MOSVGP) = Zview.(m.f)
objective(m::MOSVGP, state, y) = ELBO(m, state, y)


text:
using ITensors,
      Test

@testset "MPS Basics" begin

  N = 10
  sites = SiteSet(N,2)
  psi = MPS(sites)
  @test length(psi) == N
  @test length(MPS()) == 0

  str = split(sprint(show, psi), '\n')
  @test str[1] == "MPS"
  @test length(str) == length(psi) + 2

  @test siteindex(psi,2) == sites[2]
  @test hasindex(psi[3],linkindex(psi,2))
  @test hasindex(psi[3],linkindex(psi,3))

  psi[1] = ITensor(sites[1])
  @test hasindex(psi[1],sites[1])

  @testset "RandomMPS" begin
    phi = randomMPS(sites)
    @test hasindex(phi[1],sites[1])
    @test norm(phi[1])≈1.0
    @test hasindex(phi[4],sites[4])
    @test norm(phi[4])≈1.0
  end

  @testset "inner different MPS" begin
    phi = randomMPS(sites)
    psi = randomMPS(sites)
    phipsi = dag(phi[1])*psi[1]
    for j = 2:N
      phipsi *= dag(phi[j])*psi[j]
    end
    @test phipsi[] ≈ inner(phi,psi)

    badsites = SiteSet(N+1,2)
    badpsi = randomMPS(badsites)
    @test_throws DimensionMismatch inner(phi,badpsi)
  end

  @testset "inner same MPS" begin
    psi = randomMPS(sites)
    psidag = dag(psi)
    primelinks!(psidag)
    psipsi = psidag[1]*psi[1]
    for j = 2:N
      psipsi *= psidag[j]*psi[j]
    end
    @test psipsi[] ≈ inner(psi,psi)
  end

  sites = spinHalfSites(N)
  psi = MPS(sites)
  @test length(psi) == N # just make sure this works
  @test length(siteinds(psi)) == N

  psi = randomMPS(sites)
  position!(psi, N-1)
  @test ITensors.leftLim(psi) == N-2
  @test ITensors.rightLim(psi) == N
  position!(psi, 2)
  @test ITensors.leftLim(psi) == 1
  @test ITensors.rightLim(psi) == 3
end


text:
export geogebra

"""
`geogebra(G,file_name)` takes a `SimpleGraph` and writes out a
script to produce a drawing of this graph in GeoGebra.

Here is the secret sauce to make this work.

* Run `geogebra(G,file_name)` to save the script to `file_name`. By
  default, the file name is `geogebra.txt`.
* Copy the contents of `file_name` to the clipboard.
* Create a new GeoGebra document.
* In the **Input** zone at the bottom, enter `Button[]` to create a new
  button.
* Right click the button and select the **Object Properties...** menu option.
* Go to the **Scripting** tab and paste in the copied commands.
* Save and close the properties window.
* Press the newly created button.

Some properties of the vertices can be specified in this function with
named parameters as follows:

* `vertex_labels`: If set to `true`, the vertices in the drawing have
  labels. Default is `false` (no labels drawn). Note that the color of
  the labels matches the color of the vertices, so if you set the
  color to `white` the labels will be invisible.
* `vertex_color`: Use this to specify the fill color of the
  vertices. The default is `"black"`.
* `vertex_colors`: This is a dictionary mapping vertices to color
  names. Vertices are given the color specified. If, however, a vertex
  is missing from this dictionary, then `vertex_color` is used
  instead.
* `vertex_size`: Use this to specify the radius of the vertices. The
  default is `3`.
"""
function geogebra(
    G::SimpleGraph,
    file_name::String = "geogebra.txt";
    vertex_labels::Bool = false,
    vertex_color::String = "black",
    vertex_colors::Dict = Dict(),
    vertex_size::Int = 3,
)
    VV = vlist(G)
    n = NV(G)
    F = open(file_name, "w")
    xy = getxy(G)

    for i = 1:n
        v = VV[i]
        vname = string(v)
        (x, y) = xy[v]
        x = round(x, 3)
        y = round(y, 3)
        println(F, "v_{$i} = CopyFreeObject[Point[{$x,$y}]]")
        println(F, "SetPointSize[v_{$i}, $(string(vertex_size))]")

        color = ""
        try
            color = vertex_colors[v]
        catch
            color = vertex_color
        end
        println(F, "SetColor[v_{$i}, \"$color\"]")

        println(F, "SetCaption[v_{$i}, \"$vname\"]")
        println(F, "ShowLabel[v_{$i}, $(string(vertex_labels))]")
    end

    for i = 1:n-1
        u = VV[i]
        for j = i+1:n
            v = VV[j]
            if has(G, u, v)
                println(F, "e_{$i,$j} = Segment[v_{$i},v_{$j}]")
                println(F, "ShowLabel[e_{$i,$j}, false]")
            end
        end
    end
    println(F, "ShowAxes[false]")
    println(F, "ShowGrid[false]")
    println(F, "Execute[{\"Delete[button1]\"}]")
    close(F)
end


text:
@testset "Print Methods" begin
  @show @__LINE__
  selfdotq = :(
    for i ∈ eachindex(a)
      s += a[i] * a[i] + i
    end
  )
  lsselfdot = LoopVectorization.loopset(selfdotq)
  io = IOBuffer()
  println(io, LoopVectorization.operations(lsselfdot))
  s = String(take!(io))
  @test occursin("Operation[", s)
  @test occursin("s = 0", s)
  @test occursin("s = LoopVectorization.vfmadd_fast", s)

  AmulBq = :(
    for m ∈ 1:M, n ∈ 1:N
      C[m, n] = zero(eltype(B))
      for k ∈ 1:K
        C[m, n] += A[m, k] * B[k, n]
      end
    end
  )
  lsAmulB = LoopVectorization.loopset(AmulBq)
  println(io, LoopVectorization.operations(lsAmulB))
  s = String(take!(io))
  @test occursin("Operation[", s)
  @test occursin("C[m, n] = ", s)
  @test occursin(" = A[m, k]", s)
  @test occursin(" = B[k, n]", s)
  @test occursin(" = LoopVectorization.vfmadd_fast", s)
  # @test occursin(" = LoopVectorization.reduce_to_add", s)
end