text:
#==============================================================================#
# RJNN: Neural networks
#==============================================================================#

#==============================================================================#
# Notes:
# 1. I use terminology similar to that found in Bishop's 'Pattern Recognition
#  and Machine Learning' (Chapter 5). In particular, the input to a neuron's
#  activation function is called the 'activation', the resulting output is
#  called the 'activity' (this is MacKay's choice of terminology), and the
#  derivatives of the error function with respect to the activations (used
#  during backpropagation) are simply called 'errors'.
# 2. Data sets are stored as matrices in which the columns and rows correspond
#  to features and observations respectively; a data set containing a single
#  observation will be stored as a row vector. These matrices should always be
#  normalised column-wise (so that each column has a mean of 0 and a variance of
#  1) before being used in training, prediction, etc.
# 3. To complement the structure of the data matrix, the weights connecting
#  layers l and l+1 of a network are stored in a matrix of size m x n, where m
#  and n are the sizes of layers l and l+1 respectively. Activations are
#  computed via left multiplication (i.e., with the weight matrix appearing on
#  the right). The bias terms for each layer are held in row vectors; separating
#  them from the weight matrices in this way makes array handling slightly
#  simpler.
# 4. The activation function used here is f(x) = 1.7159 tanh(2x/3), as
#  recommended by LeCun ('Efficient BackProp', 1998). It is symmetric about the
#  x-axis and, when used with normalised inputs, has an output variance of
#  roughly 1 (because f(±1) ≈ ±1).
# 5. The gradient of the weight parameters is stored as part of the network
#  object; it contains one value per weight parameter (the derivative of the
#  total prediction error with respect to that parameter). Since the total error
#  is the sum of per-observation errors, the derivatives that make up the
#  gradient can be computed by summing simpler, per-observation derivatives.
#  Each of these computations (of the gradient based on a single observation)
#  involves both a feedforward and a backpropagation pass, so it makes sense to
#  store the activities computed during the forward pass for use during
#  backpropagation. Ideally we would like to compute the activities (one per
#  unit) for all of the observations at once, in a vectorised way, but storing
#  all of these values is not feasible. As a compromise, we compute and store
#  activities in batches (of observations); these activities are then used to
#  compute per-observation derivatives, which are added to the gradient vector.
#
# To do:
# 1. Apply PCA to the data in addition to standardisation.
# 2. Try different activation functions (ELU, SELU, etc.; see Keras). Consider
#  RJMCMC moves that change the activation function of a network/layer?
# 3. Only use one matrix for backpropagation errors.
# 4. Compute and store gradients with respect to bias variables as well. I think
#  this involves adding an extra row to the activity and error matrices.
#==============================================================================#

using LinearAlgebra

"An immutable neural network, containing both architecture and parameters."
struct Net{T<:AbstractMatrix{Float64}}
    # Weight-related matrices are type parameterised to allow for the use of
    # sparse (or other) matrix types.
    weights::Vector{T} # weights[l] connect l and l+1.
    biases::Vector{Matrix{Float64}} # per-unit biases, as row vectors.
    gradients::Vector{T} # per-weight error derivatives.
    activities::Vector{Matrix{Float64}} # batch feedforward activities.
    errors::Vector{Matrix{Float64}} # batch backpropagation errors.
end

function Net(weights, biases, batch_size)
    gradients = [similar(W) for W in weights]
    activities = Vector{Matrix{Float64}}(undef, length(weights) + 1)
    errors = Vector{Matrix{Float64}}(undef, length(weights) + 1)
    net = Net(weights, biases, gradients, activities, errors)
    batchsize!(net, batch_size)
    net
end

# The input layer does not count towards the depth (number of layers) of the
# network. It is handled by the width function for convenience sake, via the
# index l = 0.
inputs(net) = size(net.weights[1], 1)
outputs(net) = size(net.weights[end], 2)
depth(net) = length(net.weights)
width(net, l) = l == 0 ? inputs(net) : size(net.weights[l], 2)
batchsize(net) = size(net.activities[1], 1)

"Updates a network to operate with a given (maximum) batch size."
function batchsize!(net, batch_size)
    for l = 0:depth(net)
        w = width(net, l)
        net.activities[l+1] = Matrix{Float64}(undef, batch_size, w)
        net.errors[l+1] = Matrix{Float64}(undef, batch_size, w)
    end
end

function Base.show(io::IO, net::Net)
    print(io, "$(depth(net))-layer Net: $(inputs(net)) x ")
    join(io, (repr(width(net, l)) for l in 1:depth(net)), " x ")
end

function Base.show(io::IO, ::MIME"text/plain", net::Net)
    print(io, net, "\n")
    io = IOContext(io, :compact=>true, :limit=>true)
    for l in 1:depth(net)
        print(io, "$l: ", net.weights[l], " + ", net.biases[l], '\n')
    end
end

max_stars_repo_name: kevdur/rjnn

max_issues_repo_name: kevdur/rjnn

max_forks_repo_name: kevdur/rjnn


text:
module ThriftUtilsTests

using Thrift
using Base.Test

type _enum_TestEnum
    BOOLEAN::Int32
    INT32::Int32
    INT64::Int32
    INT96::Int32
    FLOAT::Int32
    DOUBLE::Int32
    BYTE_ARRAY::Int32
    FIXED_LEN_BYTE_ARRAY::Int32
end
const TestEnum = _enum_TestEnum(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4), Int32(5), Int32(6), Int32(7))

function test_enum()
    println("    enum...")
    @test enumstr(TestEnum, TestEnum.BOOLEAN) == "BOOLEAN"
    @test_throws ErrorException enumstr(TestEnum, Int32(11))
end

println("Testing utils functions...")
test_enum()
println("passed.")

end

ext: jl


text:
#=

include("./read_dsmacc.jl")

using blink
w = Blink.AtomShell.Window()
loadurl(w, string("file://" ,pwd(),"/scatter3d/index.html"))
@js w console.log("You have succesfully loaded index.html from Julia")

=#

len = length
normalise(x) = x/norm(x)
run() = eval(parse("""include("./read_dsmacc.jl")"""))
evlstr(x) = eval(parse(x))

using PyCall,DataFrames,RCall
unshift!(PyVector(pyimport("sys")["path"]), "")
@pyimport ncdata
filename = "dsmacc_out.nc"
data = ncdata.get(filename)

specs = names!(DataFrame(data["spec"]),[Symbol(i)for i in data["sc"]])
rates = names!(DataFrame(data["rate"]),[Symbol(i)for i in data["rc"]])

rates = rates[:,[Symbol(i) for i in filter(x->!ismatch(r"EMISS|DUMMY",x),data["rc"])]]
specs = specs[:,[Symbol(i) for i in filter(x->!ismatch(r"EMISS|DUMMY",x),data["sc"])]]

rates = rates[7:len(rates)]
rates = rates[:,[Symbol(i) for i in filter(x->sum(rates[x])>0,names(rates))]]

# reading in data and cleaning
joint = "\n"*(join(names(rates),"\n"))*"\n"

products= [split(i[4:len(i)],"+") for i in matchall(r"-->([A-z0-9+]*)",joint)]
reactants = [split(i[2:len(i)-1],"+") for i in matchall(r"\n([A-z0-9+]{1,60})([-->]{0,1})",joint)]

node_set = Set([ i[2:len(i)-1] for i in matchall(r"[\n+>-]{1}[A-z0-9]+[\n+>-]{1}",joint) ])
nodes = [i for i in node_set]
dict = Dict( [[i,x] for (x,i) in enumerate(nodes)])

flux = []
edges = []
for i in 1:length(reactants)
  rcol=[]
  for j in reactants[i]
    coeff = match( r"([\d\s\.]*)(\D[\d\D]*)", j)
    dummy = specs[Symbol(coeff[2])]
    try
      push!(rcol,parse(Float32,coeff[1]))
    catch
      push!(rcol,1.0f0)
    end

    prod = 1
    for k in rcol
      prod*=k
    end

    push!(flux, prod*rates[i])
    # make edges array
    for l in products[i]
        push!(edges,[j,l,i])
    end

  end
end

##################FOR EACH TIMESTEP##################
#using Blink
#w = Blink.Window()
sleep(1)

t=144
#using Plots
#unicodeplots();

#f4(x)=log(e,x),

fns=[f(x)=x,f1(x)=x*x, f2(x)=x*x*x, f3(x)=e^x, f5(x)=x^(1/2),f6(x)=x^(1/3),f7(x)=x^(1/e),f8(x)=abs(log(e,x)),f9(x)=abs(log(2,x))]
dcrpt = ["x","x^2","x^3", "exp^x", "x^{1/2}","x^{1/3}", "x^{1/exp}","|ln(x)|","|log_2(x)|"]

forward = []
for fn in 1:len(fns)
  bindata = []
  for t in 1:7:len(flux[1])

      links = filter(i -> flux[i][t]>0 , 1:len(flux))
      tflux = [log10(flux[i][t]) for i in links]
      weight = (tflux+abs(minimum(tflux)))
      weight = (weight/maximum(weight))*.98 +0.01

      weight = [ if q >0 fns[fn](q) else 0 end for q in weight]

      a = fit(Histogram, weight+0.1, nbins = 10, closed = :left)

      println(a.weights,[i for i in a.edges[1]])

      push!(bindata,a.weights)
  end
  push!(forward,bindata)
end

open("data.json", "w") do f
        write(f, "bardata="*JSON.json(forward)*"; timestamp ="*JSON.json(specs[1][2:len(specs[1])])*";description = "*JSON.json(dcrpt))
     end

#R"bin<- data.frame($(bindata));colnames(bin) <- NULL"

using Blink
w = Blink.AtomShell.Window()
loadurl(w, string("file://" ,pwd(),"/index.html"))
@js w console.log("You have succesfully loaded index.html from Julia")

#http://dataframesjl.readthedocs.io/en/latest/getting_started.html
#https://gist.github.com/gizmaa/7214002

max_stars_repo_name: wolfiex/RopaJL

max_issues_repo_name: wolfiex/RopaJL

max_forks_repo_name: wolfiex/RopaJL


text:
using FwiFlow
using Test
using PyCall
using PyPlot

matplotlib.use("Agg")


text:
export
    default_initializer

using ..Operators

"""
Deprecated
"""
function default_initializer(runtime::InstantRuntime, default_range_size::Int=10,
        placeholder_beliefs :: Dict = Dict())
    ensure_all!(runtime)
    for (pname, bel) in placeholder_beliefs
        node = get_node(runtime, pname)
        inst = current_instance(runtime, node)
        post_belief!(runtime, inst, bel)
    end
    order = topsort(get_initial_graph(get_network(runtime)))
    set_ranges!(runtime, Dict{Symbol, Score}(), default_range_size, 1, order, placeholder_beliefs)
end

#=
function default_initializer(runtime::DynamicRuntime{T}, default_time::T, default_range_size::Int=10) where {T}
    ensure_all!(runtime, default_time)
    net = get_network(runtime)
    ord = topsort(get_transition_graph(net))
    ranges = Dict{Symbol, Array{Any, 1}}()
    for var in ord
        inst = current_instance(runtime, var)
        sf = get_sfunc(inst)
        parranges = collect(map(p -> ranges[p.name], get_parents(net, var)))
        rng = support(sf, Tuple(parranges), default_range_size, Vector{output_type(sf)}())
        set_range!(runtime, inst, rng)
        ranges[var.name] = rng
    end
end

function default_initializer(runtime::DynamicRuntime{Int})
    default_initializer(runtime, 0)
end
=#

function  preserve_evidence(evid_i, typeOfEvidence)
    preserve_values = Vector{typeOfEvidence}()
    if isa(evid_i, Dict) # soft evidence
        push!(preserve_values, collect(keys(evid_i)))
    elseif isa(evid_i, NoiseEvidence) # noise evidence
        push!(preserve_values, evid_i.mean)
    elseif isa(evid_i, Function)# hard evidence
        # do nothing
    else # hard evidence
        push!(preserve_values, evid_i)
    end
    return preserve_values
end


text:
struct bounds
   πC_LB::Float64
   πC_UB::Float64
   πDER_LB::Float64
   πDER_UB::Float64
   s_LB::Float64
   s_UB::Float64
   τe_LB::Float64
   τe_UB::Float64
   τc_LB::Float64
   τc_UB::Float64
end

max_stars_repo_name: JipKim/IA

max_issues_repo_name: JipKim/IA

max_forks_repo_name: JipKim/IA


text:
#=
  Original Author: Anastasia Yendiki

  Copyright © 2022 The General Hospital Corporation (Boston, MA) "MGH"

  Terms and conditions for use, reproduction, distribution and contribution
  are found in the 'FreeSurfer Software License Agreement' contained
  in the file 'LICENSE' found in the FreeSurfer distribution, and here:

  https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense

  Reporting: freesurfer@nmr.mgh.harvard.edu
=#

using ColorTypes, DelimitedFiles, Statistics, ImageInTerminal, Plots

export LUT, color_lut, info, disp, show

const julia_red    = RGB(.796, .235, .200)
const julia_blue   = RGB(.251, .388, .847)
const julia_green  = RGB(.220, .596, .149)
const julia_purple = RGB(.584, .345, .698)

"Container for segmentation and tract look-up tables"
struct LUT
  id::Vector{Int}
  name::Vector{String}
  rgb::Vector{RGB}

  function LUT(infile::String)

    if !isfile(infile)
      error(infile * "is not a regular file")
    end

    # Read a look-up table from an input file
    # (assumed to have the format of FreeSurferColorLUT.txt)
    tab = readdlm(infile; comments=true, comment_char='#')

    # Label IDs
    id   = tab[:,1]

    # Label names
    name = tab[:,2]

    # Label display colors
    rgb  = RGB.(tab[:,3]/255, tab[:,4]/255, tab[:,5]/255)

    new(
      id,
      name,
      rgb
    )
  end
end

"The FreeSurfer color look-up table"
const global color_lut = LUT(pkgdir(FreeSurfer) * "/src/FreeSurferColorLUT.txt")

"""
    vol_to_rgb(vol::Array, maxint::Union{Number, Nothing}=nothing)

Convert an image array to an RGB/Gray array for display.

Determine how the image should be displayed:
- If all the values are IDs in the FreeSurfer color look-up table,
  the image is assumed to be is a segmentation map and is converted
  to RGB based on the FreeSurfer color look-up table.
- If the image has size 3 in any dimension, and the sum of squares
  of the values in that dimension is approx. 1 everywhere, the image
  is assumed to be a vector map and is converted to RGB based on
  vector orientation.
- Otherwise, the image is assumed to be a generic intensity map and
  is converted to grayscale, optionally clamping intensities above `maxint`.

Return an array of RGB/Gray values the same size as the input vol.
"""
function vol_to_rgb(vol::Array, maxint::Union{Number, Nothing}=nothing)

  if !any(isnothing.(indexin(unique(vol), color_lut.id)))
    if isempty(color_lut.id)
      error("FreeSurfer color look-up table is undefined")
    end

    # Assume the input is a segmentation map, get RGB of labels from LUT
    return color_lut.rgb[indexin(vol, color_lut.id)]
  end

  dim3 = findall(size(vol) .== 3)

  for idim in dim3
    if all(isapprox.(sum(vol.^2; dims=idim), 1) .|| all(vol.==0, dims=idim))
      # Assume the input is a vector map, get RGB based on vector orientation
      return RGB.(abs.(selectdim(vol,idim,1)),
                  abs.(selectdim(vol,idim,2)),
                  abs.(selectdim(vol,idim,3)))
    end
  end

  # Otherwise, assume the input is a generic intensity map
  if isnothing(maxint)
    return Gray.(vol / maximum(vol))
  else
    return Gray.(min.(vol, maxint) / maxint)
  end
end

"""
    info(mri::MRI)

Show basic info from the header of an `MRI` structure

"""
function info(mri::MRI)

  if !isempty(mri.fspec)
    println("Read from: " * mri.fspec)
  end
  println("Volume dimensions: " * string(collect(size(mri.vol))))
  println("Spatial resolution: " * string(Float64.(mri.volres)))
  if !isempty(mri.bval)
    println("b-values: " * string(Float64.(unique(mri.bval))))
  end
  println("Intensity range: " * string(Float64.([minimum(mri.vol),
                                                 maximum(mri.vol)])))

end

"""
    disp(mri::MRI, mrimod::Union{MRI, Nothing}=nothing)

Quick display of an `MRI` structure (an image slice and a summary of header
info) in the terminal window

# Arguments:
- `mri::MRI`: the main image to display
- `mrimod:MRI`: an optional image to modulate the main image by (e.g., an FA
  map to modulate a vector map)
"""
function disp(mri::MRI, mrimod::Union{MRI, Nothing}=nothing)

  # Find non-empty slices in z dimension
  iz = any(mri.vol .!= 0; dims=([1:2; 4:ndims(mri.vol)]))
  iz = reshape(iz, length(iz))
  iz = findall(iz)

  # Keep only the middle non-empty slice in z dimension
  iz = iz[Int(round(end/2))]

  # Find non-empty slices in x, y dimensions
  iy = any(mri.vol[:,:,iz,:] .!= 0; dims=([1; 3:ndims(mri.vol)]))
  iy = reshape(iy, length(iy))
  iy = findall(iy)

  ix = any(mri.vol[:,iy,iz,:] .!= 0; dims=(2:ndims(mri.vol)))
  ix = reshape(ix, length(ix))
  ix = findall(ix)

  # Keep only the area containing non-empty slices in x, y dimensions
  ix = ix[1]:ix[end]
  iy = iy[1]:iy[end]

  # Subsample to fit display in x, y dimensions
  if mri.ispermuted
    nsub = Int(ceil(length(iy) ./ displaysize(stdout)[2]))
  else
    nsub = Int(ceil(length(ix) ./ displaysize(stdout)[2]))
  end

  ix = ix[1:nsub:end]
  iy = iy[1:nsub:end]

  # Convert image to RGB/Gray array
  rgb = vol_to_rgb(mri.vol[ix, iy, iz, :])

  # Keep first frame only
  rgb = rgb[:, :, 1]

  # Optional intensity modulation
  if !isnothing(mrimod)
    if size(mrimod.vol)[1:3] != size(mri.vol)[1:3]
      error("Dimension mismatch between main image " *
            string(size(mri.vol)[1:3]) *
            " and modulation image " *
            string(size(mrimod.vol)[1:3]))
    end

    rgbmod = mrimod.vol[ix, iy, iz, 1] / maximum(mrimod.vol)
  end

  if hasfield(eltype(rgb), :r)		# RGB
    # Add α channel to make zero voxels transparent
    rgb = RGBA.(getfield.(rgb, :r), getfield.(rgb, :g), getfield.(rgb, :b),
                Float64.(rgb .!= RGB(0,0,0)))

    # Modulate intensity with (optional) second image
    if !isnothing(mrimod)
      rgb = RGBA.(getfield.(rgb, :r) .* rgbmod,
                  getfield.(rgb, :g) .* rgbmod,
                  getfield.(rgb, :b) .* rgbmod,
                  getfield.(rgb, :alpha))
    end
  else					# Gray
    # Add α channel to make zero voxels transparent
    rgb = GrayA.(getfield.(rgb, :val), Float64.(rgb .!= Gray(0)))

    # Modulate intensity with (optional) second image
    if !isnothing(mrimod)
      rgb = GrayA.(getfield.(rgb, :val) .* rgbmod,
                   getfield.(rgb, :alpha))
    end
  end

  # Display image
  if mri.ispermuted
    ImageInTerminal.imshow(rgb)
  else
    ImageInTerminal.imshow(permutedims(rgb, [2; 1]))
  end

  # Print image info
  println()
  info(mri)
end

"""
    view_axes(vox2ras::Matrix, plane::Char)

Given the `vox2ras` matrix of an image volume, return the axes along which the
volume has to be displayed to view the specified `plane` ('a': axial; 's': sagittal;
'c': coronal).
"""
function view_axes(vox2ras::Matrix, plane::Char)

  # Find orientation of image coordinate system
  orient = vox2ras_to_orient(vox2ras)

  # Find which axes of the volume correspond to the specified viewing plane
  if plane == 'a'			# Axial plane: A->P, R->L
    ax1 = findfirst(orient .== 'A' .|| orient .== 'P')
    (orient[ax1] == 'A') && (ax1 = -ax1)

    ax2 = findfirst(orient .== 'R' .|| orient .== 'L')
    (orient[ax2] == 'R') && (ax2 = -ax2)
  elseif plane == 's'			# Sagittal plane: S->I, P->A
    ax1 = findfirst(orient .== 'S' .|| orient .== 'I')
    (orient[ax1] == 'S') && (ax1 = -ax1)

    ax2 = findfirst(orient .== 'A' .|| orient .== 'P')
    (orient[ax2] == 'P') && (ax2 = -ax2)
  elseif plane == 'c'			# Coronal plane: S->I, R->L
    ax1 = findfirst(orient .== 'S' .|| orient .== 'I')
    (orient[ax1] == 'S') && (ax1 = -ax1)

    ax2 = findfirst(orient .== 'R' .|| orient .== 'L')
    (orient[ax2] == 'R') && (ax2 = -ax2)
  else
    error("Valid viewing planes are: " * string(('a', 's', 'c')))
  end

  return [ax1 ax2]
end

"""
    Base.show(mri::MRI; plane::Char='a', z::Union{Int64, Nothing}=nothing, t::Union{Int64, Nothing}=nothing, title::Union{String, Nothing}=nothing)

Show the `z`-th slice from the `t`-th frame of an `MRI`
structure, along the specified `plane` ('a': axial; 's': sagittal; 'c':
coronal).
"""
function Base.show(mri::MRI; plane::Char='a', z::Union{Int64, Nothing}=nothing, t::Union{Int64, Nothing}=nothing, title::Union{String, Nothing}=nothing, vec::Bool=true)

  # Find which axes of the volume correspond to the specified viewing plane
  ax = view_axes(mri.vox2ras, plane)

  (ax1, ax2) = abs.(ax)
  (flip1, flip2) = (ax .< 0)

  if mri.ispermuted
    (ax1 == 1) && (ax1 = 2)
    (ax1 == 2) && (ax1 = 1)
    (ax2 == 1) && (ax2 = 2)
    (ax2 == 2) && (ax2 = 1)
  end

  # Find the through-plane axis of the volume
  ax3 = setdiff(1:3, [ax1, ax2])[1]

  nx = size(mri.vol, ax2)
  ny = size(mri.vol, ax1)
  nz = size(mri.vol, ax3)

  # Set names and colors of axis labels
  if plane == 'a'			# Axial plane: A->P, R->L
    label1 = ["A", "P"]
    color1 = julia_green

    label2 = ["R", "L"]
    color2 = julia_red
  elseif plane == 's'			# Sagittal plane: S->I, P->A
    label1 = ["S", "I"]
    color1 = julia_blue

    label2 = ["P", "A"]
    color2 = julia_green
  else					# Coronal plane: S->I, R->L
    label1 = ["S", "I"]
    color1 = julia_blue

    label2 = ["R", "L"]
    color2 = julia_red
  end

  # Extract slice of interest
  isnothing(z) && (z = div(nz, 2))

  imslice = selectdim(mri.vol, ax3, z)

  # Extract volume of interest
  if isnothing(t)
    t = 1
    if size(imslice, 3) == 3
      imslice = imslice[:, :, 1:3]
    else
      imslice = imslice[:, :, 1]
    end
  else
    imslice = imslice[:, :, t]
  end

  # Max intensity for display (only has effect on grayscale intensity maps)
  maxint = Float32(1)

  if mri.nframes < mri.depth
    maxint = quantile(mri.vol[mri.vol .> 0], .999)
  else				# For larger volumes, only use middle slice
    imtmp = selectdim(mri.vol, ax3,  div(nz, 2))
    maxint = quantile(imtmp[imtmp .> 0], .999)
  end

  # Convert to RGB/Gray array
  rgb = vol_to_rgb(imslice, maxint)

  # Permute in-plane axes if needed
  if ax1 > ax2
    rgb = permutedims(rgb, [2;1])
  end

  # Flip in-plane axes if needed
  rgb = rgb[flip1 ? (end:-1:1) : (1:end), flip2 ? (end:-1:1) : (1:end)]

  # Display volume
  isnothing(title) && (title = basename(mri.fspec))

  p = Plots.plot(rgb, showaxis=false, ticks=[], aspect_ratio=1, title=title)

  # Add annotation for image axes
  Plots.annotate!(nx * .5, ny *.02, (label1[1], color1, 10), :top)
  Plots.annotate!(nx * .5, ny *.98, (label1[2], color1, 10), :bottom)
  Plots.annotate!(nx * .02, ny *.5, (label2[1], color2, 10), :left)
  Plots.annotate!(nx * .98, ny *.5, (label2[2], color2, 10), :right)

  # Add annotation for b-value and gradient vector
  blabel = ""

  if !isempty(mri.bval)
    blabel = "b=" * string(Int(round(mri.bval[t])))
  end

  if !isempty(mri.bvec)
    blabel = blabel * "\ng=[" *
                      string(round(mri.bvec[t,1] * 100)/100) * "," *
                      string(round(mri.bvec[t,2] * 100)/100) * "," *
                      string(round(mri.bvec[t,3] * 100)/100) * "]"
  end

  if !isempty(blabel)
    Plots.annotate!(nx * .02, ny * .02, (blabel, RGB(1,1,1), 9, :left, :top))
  end

  return p
end

max_stars_repo_name: freesurfer/jupysurf

max_issues_repo_name: freesurfer/jupysurf

max_forks_repo_name: freesurfer/jupysurf


text:
module JuliaTutorial

using Match

export Menu

 include("menu.jl")
 include("LinearAlgebra.jl")
 include("MatrixFactorizations.jl")
 include("TextMining.jl")
 include("RecSys.jl")

end


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

    #= ::Package:: =#

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

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

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

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

    @test_int [x^5*(a + b*sech(c + d*x^2)), x, 10, (a*x^6)/6 + (b*x^4*atan(ℯ^(c + d*x^2)))/d - (im*b*x^2*PolyLog(2, (-im)*ℯ^(c + d*x^2)))/d^2 + (im*b*x^2*PolyLog(2, im*ℯ^(c + d*x^2)))/d^2 + (im*b*PolyLog(3, (-im)*ℯ^(c + d*x^2)))/d^3 - (im*b*PolyLog(3, im*ℯ^(c + d*x^2)))/d^3]
    @test_int [x^4*(a + b*sech(c + d*x^2)), x, 2, (a*x^5)/5 + b*Unintegrable(x^4*sech(c + d*x^2), x)]
    @test_int [x^3*(a + b*sech(c + d*x^2)), x, 8, (a*x^4)/4 + (b*x^2*atan(ℯ^(c + d*x^2)))/d - ((im/2)*b*PolyLog(2, (-im)*ℯ^(c + d*x^2)))/d^2 + ((im/2)*b*PolyLog(2, im*ℯ^(c + d*x^2)))/d^2]
    @test_int [x^2*(a + b*sech(c + d*x^2)), x, 2, (a*x^3)/3 + b*Unintegrable(x^2*sech(c + d*x^2), x)]
    @test_int [x*(a + b*sech(c + d*x^2)), x, 4, (a*x^2)/2 + (b*atan(sinh(c + d*x^2)))/(2*d)]
    @test_int [(a + b*sech(c + d*x^2))/x, x, 2, a*log(x) + b*Unintegrable(sech(c + d*x^2)/x, x)]
    @test_int [(a + b*sech(c + d*x^2))/x^2, x, 2, -(a/x) + b*Unintegrable(sech(c + d*x^2)/x^2, x)]

    @test_int [x^5*(a + b*sech(c + d*x^2))^2, x, 15, (b^2*x^4)/(2*d) + (a^2*x^6)/6 + (2*a*b*x^4*atan(ℯ^(c + d*x^2)))/d - (b^2*x^2*log(1 + ℯ^(2*(c + d*x^2))))/d^2 - (2*im*a*b*x^2*PolyLog(2, (-im)*ℯ^(c + d*x^2)))/d^2 + (2*im*a*b*x^2*PolyLog(2, im*ℯ^(c + d*x^2)))/d^2 - (b^2*PolyLog(2, -ℯ^(2*(c + d*x^2))))/(2*d^3) + (2*im*a*b*PolyLog(3, (-im)*ℯ^(c + d*x^2)))/d^3 - (2*im*a*b*PolyLog(3, im*ℯ^(c + d*x^2)))/d^3 + (b^2*x^4*tanh(c + d*x^2))/(2*d)]
    @test_int [x^4*(a + b*sech(c + d*x^2))^2, x, 0, Unintegrable(x^4*(a + b*sech(c + d*x^2))^2, x)]
    @test_int [x^3*(a + b*sech(c + d*x^2))^2, x, 10, (a^2*x^4)/4 + (2*a*b*x^2*atan(ℯ^(c + d*x^2)))/d - (b^2*log(cosh(c + d*x^2)))/(2*d^2) - (im*a*b*PolyLog(2, (-im)*ℯ^(c + d*x^2)))/d^2 + (im*a*b*PolyLog(2, im*ℯ^(c + d*x^2)))/d^2 + (b^2*x^2*tanh(c + d*x^2))/(2*d)]
    @test_int [x^2*(a + b*sech(c + d*x^2))^2, x, 0, Unintegrable(x^2*(a + b*sech(c + d*x^2))^2, x)]
    @test_int [x*(a + b*sech(c + d*x^2))^2, x, 5, (a^2*x^2)/2 + (a*b*atan(sinh(c + d*x^2)))/d + (b^2*tanh(c + d*x^2))/(2*d)]
    @test_int [(a + b*sech(c + d*x^2))^2/x, x, 0, Unintegrable((a + b*sech(c + d*x^2))^2/x, x)]
    @test_int [(a + b*sech(c + d*x^2))^2/x^2, x, 0, Unintegrable((a + b*sech(c + d*x^2))^2/x^2, x)]

    @test_int [x*sech(a + b*x^2)^7, x, 5, (5*atan(sinh(a + b*x^2)))/(32*b) + (5*sech(a + b*x^2)*tanh(a + b*x^2))/(32*b) + (5*sech(a + b*x^2)^3*tanh(a + b*x^2))/(48*b) + (sech(a + b*x^2)^5*tanh(a + b*x^2))/(12*b)]

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

    @test_int [x^5/(a + b*sech(c + d*x^2)), x, 13, x^6/(6*a) - (b*x^4*log(1 + (a*ℯ^(c + d*x^2))/(b - sqrt(-a^2 + b^2))))/(2*a*sqrt(-a^2 + b^2)*d) + (b*x^4*log(1 + (a*ℯ^(c + d*x^2))/(b + sqrt(-a^2 + b^2))))/(2*a*sqrt(-a^2 + b^2)*d) - (b*x^2*PolyLog(2, -((a*ℯ^(c + d*x^2))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^2) + (b*x^2*PolyLog(2, -((a*ℯ^(c + d*x^2))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^2) + (b*PolyLog(3, -((a*ℯ^(c + d*x^2))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^3) - (b*PolyLog(3, -((a*ℯ^(c + d*x^2))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^3)]
    @test_int [x^4/(a + b*sech(c + d*x^2)), x, 0, Unintegrable(x^4/(a + b*sech(c + d*x^2)), x)]
    @test_int [x^3/(a + b*sech(c + d*x^2)), x, 11, x^4/(4*a) - (b*x^2*log(1 + (a*ℯ^(c + d*x^2))/(b - sqrt(-a^2 + b^2))))/(2*a*sqrt(-a^2 + b^2)*d) + (b*x^2*log(1 + (a*ℯ^(c + d*x^2))/(b + sqrt(-a^2 + b^2))))/(2*a*sqrt(-a^2 + b^2)*d) - (b*PolyLog(2, -((a*ℯ^(c + d*x^2))/(b - sqrt(-a^2 + b^2)))))/(2*a*sqrt(-a^2 + b^2)*d^2) + (b*PolyLog(2, -((a*ℯ^(c + d*x^2))/(b + sqrt(-a^2 + b^2)))))/(2*a*sqrt(-a^2 + b^2)*d^2)]
    @test_int [x^2/(a + b*sech(c + d*x^2)), x, 0, Unintegrable(x^2/(a + b*sech(c + d*x^2)), x)]
    @test_int [x/(a + b*sech(c + d*x^2)), x, 4, x^2/(2*a) - (b*atan((sqrt(a - b)*tanh((c + d*x^2)/2))/sqrt(a + b)))/(a*sqrt(a - b)*sqrt(a + b)*d)]
    @test_int [1/(x*(a + b*sech(c + d*x^2))), x, 0, Unintegrable(1/(x*(a + b*sech(c + d*x^2))), x)]
    @test_int [(a + b*sech(c + d*x^2))/x^2, x, 2, -(a/x) + b*Unintegrable(sech(c + d*x^2)/x^2, x)]

    @test_int [x^5/(a + b*sech(c + d*x^2))^2, x, 31, (b^2*x^4)/(2*a^2*(a^2 - b^2)*d) + x^6/(6*a^2) - (b^2*x^2*log(1 + (a*ℯ^(c + d*x^2))/(b - sqrt(-a^2 + b^2))))/(a^2*(a^2 - b^2)*d^2) + (b^3*x^4*log(1 + (a*ℯ^(c + d*x^2))/(b - sqrt(-a^2 + b^2))))/(2*a^2*(-a^2 + b^2)^(3/2)*d) - (b*x^4*log(1 + (a*ℯ^(c + d*x^2))/(b - sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*d) - (b^2*x^2*log(1 + (a*ℯ^(c + d*x^2))/(b + sqrt(-a^2 + b^2))))/(a^2*(a^2 - b^2)*d^2) - (b^3*x^4*log(1 + (a*ℯ^(c + d*x^2))/(b + sqrt(-a^2 + b^2))))/(2*a^2*(-a^2 + b^2)^(3/2)*d) + (b*x^4*log(1 + (a*ℯ^(c + d*x^2))/(b + sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*d) - (b^2*PolyLog(2, -((a*ℯ^(c + d*x^2))/(b - sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^3) + (b^3*x^2*PolyLog(2, -((a*ℯ^(c + d*x^2))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^2) - (2*b*x^2*PolyLog(2, -((a*ℯ^(c + d*x^2))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^2) - (b^2*PolyLog(2, -((a*ℯ^(c + d*x^2))/(b + sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^3) - (b^3*x^2*PolyLog(2, -((a*ℯ^(c + d*x^2))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^2) + (2*b*x^2*PolyLog(2, -((a*ℯ^(c + d*x^2))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^2) - (b^3*PolyLog(3, -((a*ℯ^(c + d*x^2))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^3) + (2*b*PolyLog(3, -((a*ℯ^(c + d*x^2))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^3) + (b^3*PolyLog(3, -((a*ℯ^(c + d*x^2))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^3) - (2*b*PolyLog(3, -((a*ℯ^(c + d*x^2))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^3) + (b^2*x^4*sinh(c + d*x^2))/(2*a*(a^2 - b^2)*d*(b + a*cosh(c + d*x^2)))]
    @test_int [x^4/(a + b*sech(c + d*x^2))^2, x, 0, Unintegrable(x^4/(a + b*sech(c + d*x^2))^2, x)]
    @test_int [x^3/(a + b*sech(c + d*x^2))^2, x, 22, x^4/(4*a^2) + (b^3*x^2*log(1 + (a*ℯ^(c + d*x^2))/(b - sqrt(-a^2 + b^2))))/(2*a^2*(-a^2 + b^2)^(3/2)*d) - (b*x^2*log(1 + (a*ℯ^(c + d*x^2))/(b - sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*d) - (b^3*x^2*log(1 + (a*ℯ^(c + d*x^2))/(b + sqrt(-a^2 + b^2))))/(2*a^2*(-a^2 + b^2)^(3/2)*d) + (b*x^2*log(1 + (a*ℯ^(c + d*x^2))/(b + sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*d) - (b^2*log(b + a*cosh(c + d*x^2)))/(2*a^2*(a^2 - b^2)*d^2) + (b^3*PolyLog(2, -((a*ℯ^(c + d*x^2))/(b - sqrt(-a^2 + b^2)))))/(2*a^2*(-a^2 + b^2)^(3/2)*d^2) - (b*PolyLog(2, -((a*ℯ^(c + d*x^2))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^2) - (b^3*PolyLog(2, -((a*ℯ^(c + d*x^2))/(b + sqrt(-a^2 + b^2)))))/(2*a^2*(-a^2 + b^2)^(3/2)*d^2) + (b*PolyLog(2, -((a*ℯ^(c + d*x^2))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^2) + (b^2*x^2*sinh(c + d*x^2))/(2*a*(a^2 - b^2)*d*(b + a*cosh(c + d*x^2)))]
    @test_int [x^2/(a + b*sech(c + d*x^2))^2, x, 0, Unintegrable(x^2/(a + b*sech(c + d*x^2))^2, x)]
    @test_int [x/(a + b*sech(c + d*x^2))^2, x, 6, x^2/(2*a^2) - (b*(2*a^2 - b^2)*atan((sqrt(a - b)*tanh((c + d*x^2)/2))/sqrt(a + b)))/(a^2*(a - b)^(3/2)*(a + b)^(3/2)*d) + (b^2*tanh(c + d*x^2))/(2*a*(a^2 - b^2)*d*(a + b*sech(c + d*x^2)))]
    @test_int [1/(x*(a + b*sech(c + d*x^2))^2), x, 0, Unintegrable(1/(x*(a + b*sech(c + d*x^2))^2), x)]
    @test_int [1/(x^2*(a + b*sech(c + d*x^2))^2), x, 0, Unintegrable(1/(x^2*(a + b*sech(c + d*x^2))^2), x)]
    @test_int [1/(x^3*(a + b*sech(c + d*x^2))^2), x, 0, Unintegrable(1/(x^3*(a + b*sech(c + d*x^2))^2), x)]

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

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

    @test_int [sech(1/x)^2/x^2, x, 3, -tanh(x^(-1))]

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

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

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

    @test_int [x^3*(a + b*sech(c + d*sqrt(x))), x, 20, (a*x^4)/4 + (4*b*x^(7/2)*atan(ℯ^(c + d*sqrt(x))))/d - ((14*im)*b*x^3*PolyLog(2, (-im)*ℯ^(c + d*sqrt(x))))/d^2 + ((14*im)*b*x^3*PolyLog(2, im*ℯ^(c + d*sqrt(x))))/d^2 + ((84*im)*b*x^(5/2)*PolyLog(3, (-im)*ℯ^(c + d*sqrt(x))))/d^3 - ((84*im)*b*x^(5/2)*PolyLog(3, im*ℯ^(c + d*sqrt(x))))/d^3 - ((420*im)*b*x^2*PolyLog(4, (-im)*ℯ^(c + d*sqrt(x))))/d^4 + ((420*im)*b*x^2*PolyLog(4, im*ℯ^(c + d*sqrt(x))))/d^4 + ((1680*im)*b*x^(3/2)*PolyLog(5, (-im)*ℯ^(c + d*sqrt(x))))/d^5 - ((1680*im)*b*x^(3/2)*PolyLog(5, im*ℯ^(c + d*sqrt(x))))/d^5 - ((5040*im)*b*x*PolyLog(6, (-im)*ℯ^(c + d*sqrt(x))))/d^6 + ((5040*im)*b*x*PolyLog(6, im*ℯ^(c + d*sqrt(x))))/d^6 + ((10080*im)*b*sqrt(x)*PolyLog(7, (-im)*ℯ^(c + d*sqrt(x))))/d^7 - ((10080*im)*b*sqrt(x)*PolyLog(7, im*ℯ^(c + d*sqrt(x))))/d^7 - ((10080*im)*b*PolyLog(8, (-im)*ℯ^(c + d*sqrt(x))))/d^8 + ((10080*im)*b*PolyLog(8, im*ℯ^(c + d*sqrt(x))))/d^8]
    @test_int [x^2*(a + b*sech(c + d*sqrt(x))), x, 16, (a*x^3)/3 + (4*b*x^(5/2)*atan(ℯ^(c + d*sqrt(x))))/d - ((10*im)*b*x^2*PolyLog(2, (-im)*ℯ^(c + d*sqrt(x))))/d^2 + ((10*im)*b*x^2*PolyLog(2, im*ℯ^(c + d*sqrt(x))))/d^2 + ((40*im)*b*x^(3/2)*PolyLog(3, (-im)*ℯ^(c + d*sqrt(x))))/d^3 - ((40*im)*b*x^(3/2)*PolyLog(3, im*ℯ^(c + d*sqrt(x))))/d^3 - ((120*im)*b*x*PolyLog(4, (-im)*ℯ^(c + d*sqrt(x))))/d^4 + ((120*im)*b*x*PolyLog(4, im*ℯ^(c + d*sqrt(x))))/d^4 + ((240*im)*b*sqrt(x)*PolyLog(5, (-im)*ℯ^(c + d*sqrt(x))))/d^5 - ((240*im)*b*sqrt(x)*PolyLog(5, im*ℯ^(c + d*sqrt(x))))/d^5 - ((240*im)*b*PolyLog(6, (-im)*ℯ^(c + d*sqrt(x))))/d^6 + ((240*im)*b*PolyLog(6, im*ℯ^(c + d*sqrt(x))))/d^6]
    @test_int [x*(a + b*sech(c + d*sqrt(x))), x, 12, (a*x^2)/2 + (4*b*x^(3/2)*atan(ℯ^(c + d*sqrt(x))))/d - ((6*im)*b*x*PolyLog(2, (-im)*ℯ^(c + d*sqrt(x))))/d^2 + ((6*im)*b*x*PolyLog(2, im*ℯ^(c + d*sqrt(x))))/d^2 + ((12*im)*b*sqrt(x)*PolyLog(3, (-im)*ℯ^(c + d*sqrt(x))))/d^3 - ((12*im)*b*sqrt(x)*PolyLog(3, im*ℯ^(c + d*sqrt(x))))/d^3 - ((12*im)*b*PolyLog(4, (-im)*ℯ^(c + d*sqrt(x))))/d^4 + ((12*im)*b*PolyLog(4, im*ℯ^(c + d*sqrt(x))))/d^4]
    @test_int [(a + b*sech(c + d*sqrt(x)))/x, x, 2, a*log(x) + b*Unintegrable(sech(c + d*sqrt(x))/x, x)]
    @test_int [(a + b*sech(c + d*sqrt(x)))/x^2, x, 2, -(a/x) + b*Unintegrable(sech(c + d*sqrt(x))/x^2, x)]

    @test_int [x^3*(a + b*sech(c + d*sqrt(x)))^2, x, 30, (2*b^2*x^(7/2))/d + (a^2*x^4)/4 + (8*a*b*x^(7/2)*atan(ℯ^(c + d*sqrt(x))))/d - (14*b^2*x^3*log(1 + ℯ^(2*(c + d*sqrt(x)))))/d^2 - (28*im*a*b*x^3*PolyLog(2, (-im)*ℯ^(c + d*sqrt(x))))/d^2 + (28*im*a*b*x^3*PolyLog(2, im*ℯ^(c + d*sqrt(x))))/d^2 - (42*b^2*x^(5/2)*PolyLog(2, -ℯ^(2*(c + d*sqrt(x)))))/d^3 + (168*im*a*b*x^(5/2)*PolyLog(3, (-im)*ℯ^(c + d*sqrt(x))))/d^3 - (168*im*a*b*x^(5/2)*PolyLog(3, im*ℯ^(c + d*sqrt(x))))/d^3 + (105*b^2*x^2*PolyLog(3, -ℯ^(2*(c + d*sqrt(x)))))/d^4 - (840*im*a*b*x^2*PolyLog(4, (-im)*ℯ^(c + d*sqrt(x))))/d^4 + (840*im*a*b*x^2*PolyLog(4, im*ℯ^(c + d*sqrt(x))))/d^4 - (210*b^2*x^(3/2)*PolyLog(4, -ℯ^(2*(c + d*sqrt(x)))))/d^5 + (3360*im*a*b*x^(3/2)*PolyLog(5, (-im)*ℯ^(c + d*sqrt(x))))/d^5 - (3360*im*a*b*x^(3/2)*PolyLog(5, im*ℯ^(c + d*sqrt(x))))/d^5 + (315*b^2*x*PolyLog(5, -ℯ^(2*(c + d*sqrt(x)))))/d^6 - (10080*im*a*b*x*PolyLog(6, (-im)*ℯ^(c + d*sqrt(x))))/d^6 + (10080*im*a*b*x*PolyLog(6, im*ℯ^(c + d*sqrt(x))))/d^6 - (315*b^2*sqrt(x)*PolyLog(6, -ℯ^(2*(c + d*sqrt(x)))))/d^7 + (20160*im*a*b*sqrt(x)*PolyLog(7, (-im)*ℯ^(c + d*sqrt(x))))/d^7 - (20160*im*a*b*sqrt(x)*PolyLog(7, im*ℯ^(c + d*sqrt(x))))/d^7 + (315*b^2*PolyLog(7, -ℯ^(2*(c + d*sqrt(x)))))/(2*d^8) - (20160*im*a*b*PolyLog(8, (-im)*ℯ^(c + d*sqrt(x))))/d^8 + (20160*im*a*b*PolyLog(8, im*ℯ^(c + d*sqrt(x))))/d^8 + (2*b^2*x^(7/2)*tanh(c + d*sqrt(x)))/d]
    @test_int [x^2*(a + b*sech(c + d*sqrt(x)))^2, x, 24, (2*b^2*x^(5/2))/d + (a^2*x^3)/3 + (8*a*b*x^(5/2)*atan(ℯ^(c + d*sqrt(x))))/d - (10*b^2*x^2*log(1 + ℯ^(2*(c + d*sqrt(x)))))/d^2 - (20*im*a*b*x^2*PolyLog(2, (-im)*ℯ^(c + d*sqrt(x))))/d^2 + (20*im*a*b*x^2*PolyLog(2, im*ℯ^(c + d*sqrt(x))))/d^2 - (20*b^2*x^(3/2)*PolyLog(2, -ℯ^(2*(c + d*sqrt(x)))))/d^3 + (80*im*a*b*x^(3/2)*PolyLog(3, (-im)*ℯ^(c + d*sqrt(x))))/d^3 - (80*im*a*b*x^(3/2)*PolyLog(3, im*ℯ^(c + d*sqrt(x))))/d^3 + (30*b^2*x*PolyLog(3, -ℯ^(2*(c + d*sqrt(x)))))/d^4 - (240*im*a*b*x*PolyLog(4, (-im)*ℯ^(c + d*sqrt(x))))/d^4 + (240*im*a*b*x*PolyLog(4, im*ℯ^(c + d*sqrt(x))))/d^4 - (30*b^2*sqrt(x)*PolyLog(4, -ℯ^(2*(c + d*sqrt(x)))))/d^5 + (480*im*a*b*sqrt(x)*PolyLog(5, (-im)*ℯ^(c + d*sqrt(x))))/d^5 - (480*im*a*b*sqrt(x)*PolyLog(5, im*ℯ^(c + d*sqrt(x))))/d^5 + (15*b^2*PolyLog(5, -ℯ^(2*(c + d*sqrt(x)))))/d^6 - (480*im*a*b*PolyLog(6, (-im)*ℯ^(c + d*sqrt(x))))/d^6 + (480*im*a*b*PolyLog(6, im*ℯ^(c + d*sqrt(x))))/d^6 + (2*b^2*x^(5/2)*tanh(c + d*sqrt(x)))/d]
    @test_int [x*(a + b*sech(c + d*sqrt(x)))^2, x, 18, (2*b^2*x^(3/2))/d + (a^2*x^2)/2 + (8*a*b*x^(3/2)*atan(ℯ^(c + d*sqrt(x))))/d - (6*b^2*x*log(1 + ℯ^(2*(c + d*sqrt(x)))))/d^2 - (12*im*a*b*x*PolyLog(2, (-im)*ℯ^(c + d*sqrt(x))))/d^2 + (12*im*a*b*x*PolyLog(2, im*ℯ^(c + d*sqrt(x))))/d^2 - (6*b^2*sqrt(x)*PolyLog(2, -ℯ^(2*(c + d*sqrt(x)))))/d^3 + (24*im*a*b*sqrt(x)*PolyLog(3, (-im)*ℯ^(c + d*sqrt(x))))/d^3 - (24*im*a*b*sqrt(x)*PolyLog(3, im*ℯ^(c + d*sqrt(x))))/d^3 + (3*b^2*PolyLog(3, -ℯ^(2*(c + d*sqrt(x)))))/d^4 - (24*im*a*b*PolyLog(4, (-im)*ℯ^(c + d*sqrt(x))))/d^4 + (24*im*a*b*PolyLog(4, im*ℯ^(c + d*sqrt(x))))/d^4 + (2*b^2*x^(3/2)*tanh(c + d*sqrt(x)))/d]
    @test_int [(a + b*sech(c + d*sqrt(x)))^2/x, x, 0, Unintegrable((a + b*sech(c + d*sqrt(x)))^2/x, x)]
    @test_int [(a + b*sech(c + d*sqrt(x)))^2/x^2, x, 0, Unintegrable((a + b*sech(c + d*sqrt(x)))^2/x^2, x)]

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

    @test_int [x^3/(a + b*sech(c + d*sqrt(x))), x, 23, x^4/(4*a) - (2*b*x^(7/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*d) + (2*b*x^(7/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*d) - (14*b*x^3*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^2) + (14*b*x^3*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^2) + (84*b*x^(5/2)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^3) - (84*b*x^(5/2)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^3) - (420*b*x^2*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^4) + (420*b*x^2*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^4) + (1680*b*x^(3/2)*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^5) - (1680*b*x^(3/2)*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^5) - (5040*b*x*PolyLog(6, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^6) + (5040*b*x*PolyLog(6, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^6) + (10080*b*sqrt(x)*PolyLog(7, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^7) - (10080*b*sqrt(x)*PolyLog(7, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^7) - (10080*b*PolyLog(8, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^8) + (10080*b*PolyLog(8, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^8)]
    @test_int [x^2/(a + b*sech(c + d*sqrt(x))), x, 19, x^3/(3*a) - (2*b*x^(5/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*d) + (2*b*x^(5/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*d) - (10*b*x^2*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^2) + (10*b*x^2*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^2) + (40*b*x^(3/2)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^3) - (40*b*x^(3/2)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^3) - (120*b*x*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^4) + (120*b*x*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^4) + (240*b*sqrt(x)*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^5) - (240*b*sqrt(x)*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^5) - (240*b*PolyLog(6, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^6) + (240*b*PolyLog(6, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^6)]
    @test_int [x/(a + b*sech(c + d*sqrt(x))), x, 15, x^2/(2*a) - (2*b*x^(3/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*d) + (2*b*x^(3/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*d) - (6*b*x*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^2) + (6*b*x*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^2) + (12*b*sqrt(x)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^3) - (12*b*sqrt(x)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^3) - (12*b*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^4) + (12*b*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^4)]
    @test_int [1/(x*(a + b*sech(c + d*sqrt(x)))), x, 0, Unintegrable(1/(x*(a + b*sech(c + d*sqrt(x)))), x)]
    @test_int [(a + b*sech(c + d*sqrt(x)))/x^2, x, 2, -(a/x) + b*Unintegrable(sech(c + d*sqrt(x))/x^2, x)]

    @test_int [x^3/(a + b*sech(c + d*sqrt(x)))^2, x, 61, (2*b^2*x^(7/2))/(a^2*(a^2 - b^2)*d) + x^4/(4*a^2) - (14*b^2*x^3*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a^2*(a^2 - b^2)*d^2) + (2*b^3*x^(7/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*d) - (4*b*x^(7/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*d) - (14*b^2*x^3*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a^2*(a^2 - b^2)*d^2) - (2*b^3*x^(7/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*d) + (4*b*x^(7/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*d) - (84*b^2*x^(5/2)*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^3) + (14*b^3*x^3*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^2) - (28*b*x^3*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^2) - (84*b^2*x^(5/2)*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^3) - (14*b^3*x^3*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^2) + (28*b*x^3*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^2) + (420*b^2*x^2*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^4) - (84*b^3*x^(5/2)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^3) + (168*b*x^(5/2)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^3) + (420*b^2*x^2*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^4) + (84*b^3*x^(5/2)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^3) - (168*b*x^(5/2)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^3) - (1680*b^2*x^(3/2)*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^5) + (420*b^3*x^2*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^4) - (840*b*x^2*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^4) - (1680*b^2*x^(3/2)*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^5) - (420*b^3*x^2*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^4) + (840*b*x^2*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^4) + (5040*b^2*x*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^6) - (1680*b^3*x^(3/2)*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^5) + (3360*b*x^(3/2)*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^5) + (5040*b^2*x*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^6) + (1680*b^3*x^(3/2)*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^5) - (3360*b*x^(3/2)*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^5) - (10080*b^2*sqrt(x)*PolyLog(6, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^7) + (5040*b^3*x*PolyLog(6, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^6) - (10080*b*x*PolyLog(6, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^6) - (10080*b^2*sqrt(x)*PolyLog(6, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^7) - (5040*b^3*x*PolyLog(6, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^6) + (10080*b*x*PolyLog(6, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^6) + (10080*b^2*PolyLog(7, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^8) - (10080*b^3*sqrt(x)*PolyLog(7, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^7) + (20160*b*sqrt(x)*PolyLog(7, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^7) + (10080*b^2*PolyLog(7, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^8) + (10080*b^3*sqrt(x)*PolyLog(7, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^7) - (20160*b*sqrt(x)*PolyLog(7, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^7) + (10080*b^3*PolyLog(8, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^8) - (20160*b*PolyLog(8, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^8) - (10080*b^3*PolyLog(8, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^8) + (20160*b*PolyLog(8, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^8) + (2*b^2*x^(7/2)*sinh(c + d*sqrt(x)))/(a*(a^2 - b^2)*d*(b + a*cosh(c + d*sqrt(x))))]
    @test_int [x^2/(a + b*sech(c + d*sqrt(x)))^2, x, 49, (2*b^2*x^(5/2))/(a^2*(a^2 - b^2)*d) + x^3/(3*a^2) - (10*b^2*x^2*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a^2*(a^2 - b^2)*d^2) + (2*b^3*x^(5/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*d) - (4*b*x^(5/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*d) - (10*b^2*x^2*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a^2*(a^2 - b^2)*d^2) - (2*b^3*x^(5/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*d) + (4*b*x^(5/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*d) - (40*b^2*x^(3/2)*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^3) + (10*b^3*x^2*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^2) - (20*b*x^2*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^2) - (40*b^2*x^(3/2)*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^3) - (10*b^3*x^2*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^2) + (20*b*x^2*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^2) + (120*b^2*x*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^4) - (40*b^3*x^(3/2)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^3) + (80*b*x^(3/2)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^3) + (120*b^2*x*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^4) + (40*b^3*x^(3/2)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^3) - (80*b*x^(3/2)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^3) - (240*b^2*sqrt(x)*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^5) + (120*b^3*x*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^4) - (240*b*x*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^4) - (240*b^2*sqrt(x)*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^5) - (120*b^3*x*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^4) + (240*b*x*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^4) + (240*b^2*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^6) - (240*b^3*sqrt(x)*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^5) + (480*b*sqrt(x)*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^5) + (240*b^2*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^6) + (240*b^3*sqrt(x)*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^5) - (480*b*sqrt(x)*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^5) + (240*b^3*PolyLog(6, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^6) - (480*b*PolyLog(6, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^6) - (240*b^3*PolyLog(6, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^6) + (480*b*PolyLog(6, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^6) + (2*b^2*x^(5/2)*sinh(c + d*sqrt(x)))/(a*(a^2 - b^2)*d*(b + a*cosh(c + d*sqrt(x))))]
    @test_int [x/(a + b*sech(c + d*sqrt(x)))^2, x, 37, (2*b^2*x^(3/2))/(a^2*(a^2 - b^2)*d) + x^2/(2*a^2) - (6*b^2*x*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a^2*(a^2 - b^2)*d^2) + (2*b^3*x^(3/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*d) - (4*b*x^(3/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*d) - (6*b^2*x*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a^2*(a^2 - b^2)*d^2) - (2*b^3*x^(3/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*d) + (4*b*x^(3/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*d) - (12*b^2*sqrt(x)*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^3) + (6*b^3*x*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^2) - (12*b*x*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^2) - (12*b^2*sqrt(x)*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^3) - (6*b^3*x*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^2) + (12*b*x*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^2) + (12*b^2*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^4) - (12*b^3*sqrt(x)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^3) + (24*b*sqrt(x)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^3) + (12*b^2*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^4) + (12*b^3*sqrt(x)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^3) - (24*b*sqrt(x)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^3) + (12*b^3*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^4) - (24*b*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^4) - (12*b^3*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^4) + (24*b*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^4) + (2*b^2*x^(3/2)*sinh(c + d*sqrt(x)))/(a*(a^2 - b^2)*d*(b + a*cosh(c + d*sqrt(x))))]
    @test_int [1/(x*(a + b*sech(c + d*sqrt(x)))^2), x, 0, Unintegrable(1/(x*(a + b*sech(c + d*sqrt(x)))^2), x)]
    @test_int [1/(x^2*(a + b*sech(c + d*sqrt(x)))^2), x, 0, Unintegrable(1/(x^2*(a + b*sech(c + d*sqrt(x)))^2), x)]

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

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

    @test_int [x^(3/2)*(a + b*sech(c + d*sqrt(x))), x, 14, (2*a*x^(5/2))/5 + (4*b*x^2*atan(ℯ^(c + d*sqrt(x))))/d - ((8*im)*b*x^(3/2)*PolyLog(2, (-im)*ℯ^(c + d*sqrt(x))))/d^2 + ((8*im)*b*x^(3/2)*PolyLog(2, im*ℯ^(c + d*sqrt(x))))/d^2 + ((24*im)*b*x*PolyLog(3, (-im)*ℯ^(c + d*sqrt(x))))/d^3 - ((24*im)*b*x*PolyLog(3, im*ℯ^(c + d*sqrt(x))))/d^3 - ((48*im)*b*sqrt(x)*PolyLog(4, (-im)*ℯ^(c + d*sqrt(x))))/d^4 + ((48*im)*b*sqrt(x)*PolyLog(4, im*ℯ^(c + d*sqrt(x))))/d^4 + ((48*im)*b*PolyLog(5, (-im)*ℯ^(c + d*sqrt(x))))/d^5 - ((48*im)*b*PolyLog(5, im*ℯ^(c + d*sqrt(x))))/d^5]
    @test_int [sqrt(x)*(a + b*sech(c + d*sqrt(x))), x, 10, (2*a*x^(3/2))/3 + (4*b*x*atan(ℯ^(c + d*sqrt(x))))/d - ((4*im)*b*sqrt(x)*PolyLog(2, (-im)*ℯ^(c + d*sqrt(x))))/d^2 + ((4*im)*b*sqrt(x)*PolyLog(2, im*ℯ^(c + d*sqrt(x))))/d^2 + ((4*im)*b*PolyLog(3, (-im)*ℯ^(c + d*sqrt(x))))/d^3 - ((4*im)*b*PolyLog(3, im*ℯ^(c + d*sqrt(x))))/d^3]
    @test_int [(a + b*sech(c + d*sqrt(x)))/sqrt(x), x, 4, 2*a*sqrt(x) + (2*b*atan(sinh(c + d*sqrt(x))))/d]
    @test_int [(a + b*sech(c + d*sqrt(x)))/x^(3/2), x, 2, (-2*a)/sqrt(x) + b*Unintegrable(sech(c + d*sqrt(x))/x^(3/2), x)]
    @test_int [(a + b*sech(c + d*sqrt(x)))/x^(5/2), x, 2, (-2*a)/(3*x^(3/2)) + b*Unintegrable(sech(c + d*sqrt(x))/x^(5/2), x)]

    @test_int [x^(3/2)*(a + b*sech(c + d*sqrt(x)))^2, x, 21, (2*b^2*x^2)/d + (2/5)*a^2*x^(5/2) + (8*a*b*x^2*atan(ℯ^(c + d*sqrt(x))))/d - (8*b^2*x^(3/2)*log(1 + ℯ^(2*(c + d*sqrt(x)))))/d^2 - (16*im*a*b*x^(3/2)*PolyLog(2, (-im)*ℯ^(c + d*sqrt(x))))/d^2 + (16*im*a*b*x^(3/2)*PolyLog(2, im*ℯ^(c + d*sqrt(x))))/d^2 - (12*b^2*x*PolyLog(2, -ℯ^(2*(c + d*sqrt(x)))))/d^3 + (48*im*a*b*x*PolyLog(3, (-im)*ℯ^(c + d*sqrt(x))))/d^3 - (48*im*a*b*x*PolyLog(3, im*ℯ^(c + d*sqrt(x))))/d^3 + (12*b^2*sqrt(x)*PolyLog(3, -ℯ^(2*(c + d*sqrt(x)))))/d^4 - (96*im*a*b*sqrt(x)*PolyLog(4, (-im)*ℯ^(c + d*sqrt(x))))/d^4 + (96*im*a*b*sqrt(x)*PolyLog(4, im*ℯ^(c + d*sqrt(x))))/d^4 - (6*b^2*PolyLog(4, -ℯ^(2*(c + d*sqrt(x)))))/d^5 + (96*im*a*b*PolyLog(5, (-im)*ℯ^(c + d*sqrt(x))))/d^5 - (96*im*a*b*PolyLog(5, im*ℯ^(c + d*sqrt(x))))/d^5 + (2*b^2*x^2*tanh(c + d*sqrt(x)))/d]
    @test_int [sqrt(x)*(a + b*sech(c + d*sqrt(x)))^2, x, 15, (2*b^2*x)/d + (2/3)*a^2*x^(3/2) + (8*a*b*x*atan(ℯ^(c + d*sqrt(x))))/d - (4*b^2*sqrt(x)*log(1 + ℯ^(2*(c + d*sqrt(x)))))/d^2 - (8*im*a*b*sqrt(x)*PolyLog(2, (-im)*ℯ^(c + d*sqrt(x))))/d^2 + (8*im*a*b*sqrt(x)*PolyLog(2, im*ℯ^(c + d*sqrt(x))))/d^2 - (2*b^2*PolyLog(2, -ℯ^(2*(c + d*sqrt(x)))))/d^3 + (8*im*a*b*PolyLog(3, (-im)*ℯ^(c + d*sqrt(x))))/d^3 - (8*im*a*b*PolyLog(3, im*ℯ^(c + d*sqrt(x))))/d^3 + (2*b^2*x*tanh(c + d*sqrt(x)))/d]
    @test_int [(a + b*sech(c + d*sqrt(x)))^2/sqrt(x), x, 5, 2*a^2*sqrt(x) + (4*a*b*atan(sinh(c + d*sqrt(x))))/d + (2*b^2*tanh(c + d*sqrt(x)))/d]
    @test_int [(a + b*sech(c + d*sqrt(x)))^2/x^(3/2), x, 0, Unintegrable((a + b*sech(c + d*sqrt(x)))^2/x^(3/2), x)]
    @test_int [(a + b*sech(c + d*sqrt(x)))^2/x^(5/2), x, 0, Unintegrable((a + b*sech(c + d*sqrt(x)))^2/x^(5/2), x)]

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

    @test_int [x^(3/2)/(a + b*sech(c + d*sqrt(x))), x, 17, (2*x^(5/2))/(5*a) - (2*b*x^2*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*d) + (2*b*x^2*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*d) - (8*b*x^(3/2)*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^2) + (8*b*x^(3/2)*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^2) + (24*b*x*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^3) - (24*b*x*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^3) - (48*b*sqrt(x)*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^4) + (48*b*sqrt(x)*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^4) + (48*b*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^5) - (48*b*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^5)]
    @test_int [sqrt(x)/(a + b*sech(c + d*sqrt(x))), x, 13, (2*x^(3/2))/(3*a) - (2*b*x*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*d) + (2*b*x*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*d) - (4*b*sqrt(x)*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^2) + (4*b*sqrt(x)*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^2) + (4*b*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^3) - (4*b*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^3)]
    @test_int [1/(sqrt(x)*(a + b*sech(c + d*sqrt(x)))), x, 4, (2*sqrt(x))/a - (4*b*atan((sqrt(a - b)*tanh((c + d*sqrt(x))/2))/sqrt(a + b)))/(a*sqrt(a - b)*sqrt(a + b)*d)]
    @test_int [1/(x^(3/2)*(a + b*sech(c + d*sqrt(x)))), x, 0, Unintegrable(1/(x^(3/2)*(a + b*sech(c + d*sqrt(x)))), x)]
    @test_int [1/(x^(5/2)*(a + b*sech(c + d*sqrt(x)))), x, 0, Unintegrable(1/(x^(5/2)*(a + b*sech(c + d*sqrt(x)))), x)]

    @test_int [x^(3/2)/(a + b*sech(c + d*sqrt(x)))^2, x, 43, (2*b^2*x^2)/(a^2*(a^2 - b^2)*d) + (2*x^(5/2))/(5*a^2) - (8*b^2*x^(3/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a^2*(a^2 - b^2)*d^2) + (2*b^3*x^2*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*d) - (4*b*x^2*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*d) - (8*b^2*x^(3/2)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a^2*(a^2 - b^2)*d^2) - (2*b^3*x^2*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*d) + (4*b*x^2*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*d) - (24*b^2*x*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^3) + (8*b^3*x^(3/2)*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^2) - (16*b*x^(3/2)*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^2) - (24*b^2*x*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^3) - (8*b^3*x^(3/2)*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^2) + (16*b*x^(3/2)*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^2) + (48*b^2*sqrt(x)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^4) - (24*b^3*x*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^3) + (48*b*x*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^3) + (48*b^2*sqrt(x)*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^4) + (24*b^3*x*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^3) - (48*b*x*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^3) - (48*b^2*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^5) + (48*b^3*sqrt(x)*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^4) - (96*b*sqrt(x)*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^4) - (48*b^2*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^5) - (48*b^3*sqrt(x)*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^4) + (96*b*sqrt(x)*PolyLog(4, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^4) - (48*b^3*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^5) + (96*b*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^5) + (48*b^3*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^5) - (96*b*PolyLog(5, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^5) + (2*b^2*x^2*sinh(c + d*sqrt(x)))/(a*(a^2 - b^2)*d*(b + a*cosh(c + d*sqrt(x))))]
    @test_int [sqrt(x)/(a + b*sech(c + d*sqrt(x)))^2, x, 31, (2*b^2*x)/(a^2*(a^2 - b^2)*d) + (2*x^(3/2))/(3*a^2) - (4*b^2*sqrt(x)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a^2*(a^2 - b^2)*d^2) + (2*b^3*x*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*d) - (4*b*x*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*d) - (4*b^2*sqrt(x)*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a^2*(a^2 - b^2)*d^2) - (2*b^3*x*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*d) + (4*b*x*log(1 + (a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*d) - (4*b^2*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^3) + (4*b^3*sqrt(x)*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^2) - (8*b*sqrt(x)*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^2) - (4*b^2*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^3) - (4*b^3*sqrt(x)*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^2) + (8*b*sqrt(x)*PolyLog(2, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^2) - (4*b^3*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^3) + (8*b*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^3) + (4*b^3*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^3) - (8*b*PolyLog(3, -((a*ℯ^(c + d*sqrt(x)))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^3) + (2*b^2*x*sinh(c + d*sqrt(x)))/(a*(a^2 - b^2)*d*(b + a*cosh(c + d*sqrt(x))))]
    @test_int [1/(sqrt(x)*(a + b*sech(c + d*sqrt(x)))^2), x, 6, (2*sqrt(x))/a^2 - (4*b*(2*a^2 - b^2)*atan((sqrt(a - b)*tanh((c + d*sqrt(x))/2))/sqrt(a + b)))/(a^2*(a - b)^(3/2)*(a + b)^(3/2)*d) + (2*b^2*tanh(c + d*sqrt(x)))/(a*(a^2 - b^2)*d*(a + b*sech(c + d*sqrt(x))))]
    @test_int [1/(x^(3/2)*(a + b*sech(c + d*sqrt(x)))^2), x, 0, Unintegrable(1/(x^(3/2)*(a + b*sech(c + d*sqrt(x)))^2), x)]
    @test_int [1/(x^(5/2)*(a + b*sech(c + d*sqrt(x)))^2), x, 0, Unintegrable(1/(x^(5/2)*(a + b*sech(c + d*sqrt(x)))^2), x)]

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

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

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

    @test_int [(e*x)^(-1 + n)*(a + b*sech(c + d*x^n)), x, 5, (a*(e*x)^n)/(e*n) + (b*(e*x)^n*atan(sinh(c + d*x^n)))/(d*e*n*x^n)]
    @test_int [(e*x)^(-1 + 2*n)*(a + b*sech(c + d*x^n)), x, 9, (a*(e*x)^(2*n))/(2*e*n) + (2*b*(e*x)^(2*n)*atan(ℯ^(c + d*x^n)))/(d*e*n*x^n) - (im*b*(e*x)^(2*n)*PolyLog(2, (-im)*ℯ^(c + d*x^n)))/(d^2*e*n*x^(2*n)) + (im*b*(e*x)^(2*n)*PolyLog(2, im*ℯ^(c + d*x^n)))/(d^2*e*n*x^(2*n))]
    @test_int [(e*x)^(-1 + 3*n)*(a + b*sech(c + d*x^n)), x, 11, (a*(e*x)^(3*n))/(3*e*n) + (2*b*(e*x)^(3*n)*atan(ℯ^(c + d*x^n)))/(d*e*n*x^n) - ((2*im)*b*(e*x)^(3*n)*PolyLog(2, (-im)*ℯ^(c + d*x^n)))/(d^2*e*n*x^(2*n)) + ((2*im)*b*(e*x)^(3*n)*PolyLog(2, im*ℯ^(c + d*x^n)))/(d^2*e*n*x^(2*n)) + ((2*im)*b*(e*x)^(3*n)*PolyLog(3, (-im)*ℯ^(c + d*x^n)))/(d^3*e*n*x^(3*n)) - ((2*im)*b*(e*x)^(3*n)*PolyLog(3, im*ℯ^(c + d*x^n)))/(d^3*e*n*x^(3*n))]

    @test_int [(e*x)^(-1 + n)*(a + b*sech(c + d*x^n))^2, x, 6, (a^2*(e*x)^n)/(e*n) + (2*a*b*(e*x)^n*atan(sinh(c + d*x^n)))/(d*e*n*x^n) + (b^2*(e*x)^n*tanh(c + d*x^n))/(d*e*n*x^n)]
    @test_int [(e*x)^(-1 + 2*n)*(a + b*sech(c + d*x^n))^2, x, 11, (a^2*(e*x)^(2*n))/(2*e*n) + (4*a*b*(e*x)^(2*n)*atan(ℯ^(c + d*x^n)))/(d*e*n*x^n) - (b^2*(e*x)^(2*n)*log(cosh(c + d*x^n)))/(d^2*e*n*x^(2*n)) - ((2*im)*a*b*(e*x)^(2*n)*PolyLog(2, (-im)*ℯ^(c + d*x^n)))/(d^2*e*n*x^(2*n)) + ((2*im)*a*b*(e*x)^(2*n)*PolyLog(2, im*ℯ^(c + d*x^n)))/(d^2*e*n*x^(2*n)) + (b^2*(e*x)^(2*n)*tanh(c + d*x^n))/(d*e*n*x^n)]
    @test_int [(e*x)^(-1 + 3*n)*(a + b*sech(c + d*x^n))^2, x, 16, (a^2*(e*x)^(3*n))/(3*e*n) + (b^2*(e*x)^(3*n))/(x^n*(d*e*n)) + (4*a*b*(e*x)^(3*n)*atan(ℯ^(c + d*x^n)))/(x^n*(d*e*n)) - (2*b^2*(e*x)^(3*n)*log(1 + ℯ^(2*(c + d*x^n))))/(x^(2*n)*(d^2*e*n)) - (4*im*a*b*(e*x)^(3*n)*PolyLog(2, (-im)*ℯ^(c + d*x^n)))/(x^(2*n)*(d^2*e*n)) + (4*im*a*b*(e*x)^(3*n)*PolyLog(2, im*ℯ^(c + d*x^n)))/(x^(2*n)*(d^2*e*n)) - (b^2*(e*x)^(3*n)*PolyLog(2, -ℯ^(2*(c + d*x^n))))/(x^(3*n)*(d^3*e*n)) + (4*im*a*b*(e*x)^(3*n)*PolyLog(3, (-im)*ℯ^(c + d*x^n)))/(x^(3*n)*(d^3*e*n)) - (4*im*a*b*(e*x)^(3*n)*PolyLog(3, im*ℯ^(c + d*x^n)))/(x^(3*n)*(d^3*e*n)) + (b^2*(e*x)^(3*n)*tanh(c + d*x^n))/(x^n*(d*e*n))]

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

    @test_int [(e*x)^(-1 + n)/(a + b*sech(c + d*x^n)), x, 5, (e*x)^n/(a*e*n) - (2*b*(e*x)^n*atan((sqrt(a - b)*tanh((c + d*x^n)/2))/sqrt(a + b)))/(a*sqrt(a - b)*sqrt(a + b)*d*e*n*x^n)]
    @test_int [(e*x)^(-1 + 2*n)/(a + b*sech(c + d*x^n)), x, 12, (e*x)^(2*n)/(2*a*e*n) - (b*(e*x)^(2*n)*log(1 + (a*ℯ^(c + d*x^n))/(b - sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*d*e*n*x^n) + (b*(e*x)^(2*n)*log(1 + (a*ℯ^(c + d*x^n))/(b + sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*d*e*n*x^n) - (b*(e*x)^(2*n)*PolyLog(2, -((a*ℯ^(c + d*x^n))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^2*e*n*x^(2*n)) + (b*(e*x)^(2*n)*PolyLog(2, -((a*ℯ^(c + d*x^n))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^2*e*n*x^(2*n))]
    @test_int [(e*x)^(-1 + 3*n)/(a + b*sech(c + d*x^n)), x, 14, (e*x)^(3*n)/(3*a*e*n) - (b*(e*x)^(3*n)*log(1 + (a*ℯ^(c + d*x^n))/(b - sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*d*e*n*x^n) + (b*(e*x)^(3*n)*log(1 + (a*ℯ^(c + d*x^n))/(b + sqrt(-a^2 + b^2))))/(a*sqrt(-a^2 + b^2)*d*e*n*x^n) - (2*b*(e*x)^(3*n)*PolyLog(2, -((a*ℯ^(c + d*x^n))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^2*e*n*x^(2*n)) + (2*b*(e*x)^(3*n)*PolyLog(2, -((a*ℯ^(c + d*x^n))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^2*e*n*x^(2*n)) + (2*b*(e*x)^(3*n)*PolyLog(3, -((a*ℯ^(c + d*x^n))/(b - sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^3*e*n*x^(3*n)) - (2*b*(e*x)^(3*n)*PolyLog(3, -((a*ℯ^(c + d*x^n))/(b + sqrt(-a^2 + b^2)))))/(a*sqrt(-a^2 + b^2)*d^3*e*n*x^(3*n))]

    @test_int [(e*x)^(-1 + n)/(a + b*sech(c + d*x^n))^2, x, 7, (e*x)^n/(a^2*e*n) - (2*b*(2*a^2 - b^2)*(e*x)^n*atan((sqrt(a - b)*tanh((c + d*x^n)/2))/sqrt(a + b)))/(a^2*(a - b)^(3/2)*(a + b)^(3/2)*d*e*n*x^n) + (b^2*(e*x)^n*tanh(c + d*x^n))/(a*(a^2 - b^2)*d*e*n*x^n*(a + b*sech(c + d*x^n)))]
    @test_int [(e*x)^(-1 + 2*n)/(a + b*sech(c + d*x^n))^2, x, 23, (e*x)^(2*n)/(2*a^2*e*n) + (b^3*(e*x)^(2*n)*log(1 + (a*ℯ^(c + d*x^n))/(b - sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*d*e*n*x^n) - (2*b*(e*x)^(2*n)*log(1 + (a*ℯ^(c + d*x^n))/(b - sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*d*e*n*x^n) - (b^3*(e*x)^(2*n)*log(1 + (a*ℯ^(c + d*x^n))/(b + sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*d*e*n*x^n) + (2*b*(e*x)^(2*n)*log(1 + (a*ℯ^(c + d*x^n))/(b + sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*d*e*n*x^n) - (b^2*(e*x)^(2*n)*log(b + a*cosh(c + d*x^n)))/(a^2*(a^2 - b^2)*d^2*e*n*x^(2*n)) + (b^3*(e*x)^(2*n)*PolyLog(2, -((a*ℯ^(c + d*x^n))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^2*e*n*x^(2*n)) - (2*b*(e*x)^(2*n)*PolyLog(2, -((a*ℯ^(c + d*x^n))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^2*e*n*x^(2*n)) - (b^3*(e*x)^(2*n)*PolyLog(2, -((a*ℯ^(c + d*x^n))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^2*e*n*x^(2*n)) + (2*b*(e*x)^(2*n)*PolyLog(2, -((a*ℯ^(c + d*x^n))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^2*e*n*x^(2*n)) + (b^2*(e*x)^(2*n)*sinh(c + d*x^n))/(a*(a^2 - b^2)*d*e*n*x^n*(b + a*cosh(c + d*x^n)))]
    @test_int [(e*x)^(-1 + 3*n)/(a + b*sech(c + d*x^n))^2, x, 32, (e*x)^(3*n)/(3*a^2*e*n) + (b^2*(e*x)^(3*n))/(a^2*(a^2 - b^2)*d*e*n*x^n) - (2*b^2*(e*x)^(3*n)*log(1 + (a*ℯ^(c + d*x^n))/(b - sqrt(-a^2 + b^2))))/(a^2*(a^2 - b^2)*d^2*e*n*x^(2*n)) + (b^3*(e*x)^(3*n)*log(1 + (a*ℯ^(c + d*x^n))/(b - sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*d*e*n*x^n) - (2*b*(e*x)^(3*n)*log(1 + (a*ℯ^(c + d*x^n))/(b - sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*d*e*n*x^n) - (2*b^2*(e*x)^(3*n)*log(1 + (a*ℯ^(c + d*x^n))/(b + sqrt(-a^2 + b^2))))/(a^2*(a^2 - b^2)*d^2*e*n*x^(2*n)) - (b^3*(e*x)^(3*n)*log(1 + (a*ℯ^(c + d*x^n))/(b + sqrt(-a^2 + b^2))))/(a^2*(-a^2 + b^2)^(3/2)*d*e*n*x^n) + (2*b*(e*x)^(3*n)*log(1 + (a*ℯ^(c + d*x^n))/(b + sqrt(-a^2 + b^2))))/(a^2*sqrt(-a^2 + b^2)*d*e*n*x^n) - (2*b^2*(e*x)^(3*n)*PolyLog(2, -((a*ℯ^(c + d*x^n))/(b - sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^3*e*n*x^(3*n)) + (2*b^3*(e*x)^(3*n)*PolyLog(2, -((a*ℯ^(c + d*x^n))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^2*e*n*x^(2*n)) - (4*b*(e*x)^(3*n)*PolyLog(2, -((a*ℯ^(c + d*x^n))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^2*e*n*x^(2*n)) - (2*b^2*(e*x)^(3*n)*PolyLog(2, -((a*ℯ^(c + d*x^n))/(b + sqrt(-a^2 + b^2)))))/(a^2*(a^2 - b^2)*d^3*e*n*x^(3*n)) - (2*b^3*(e*x)^(3*n)*PolyLog(2, -((a*ℯ^(c + d*x^n))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^2*e*n*x^(2*n)) + (4*b*(e*x)^(3*n)*PolyLog(2, -((a*ℯ^(c + d*x^n))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^2*e*n*x^(2*n)) - (2*b^3*(e*x)^(3*n)*PolyLog(3, -((a*ℯ^(c + d*x^n))/(b - sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^3*e*n*x^(3*n)) + (4*b*(e*x)^(3*n)*PolyLog(3, -((a*ℯ^(c + d*x^n))/(b - sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^3*e*n*x^(3*n)) + (2*b^3*(e*x)^(3*n)*PolyLog(3, -((a*ℯ^(c + d*x^n))/(b + sqrt(-a^2 + b^2)))))/(a^2*(-a^2 + b^2)^(3/2)*d^3*e*n*x^(3*n)) - (4*b*(e*x)^(3*n)*PolyLog(3, -((a*ℯ^(c + d*x^n))/(b + sqrt(-a^2 + b^2)))))/(a^2*sqrt(-a^2 + b^2)*d^3*e*n*x^(3*n)) + (b^2*(e*x)^(3*n)*sinh(c + d*x^n))/(a*(a^2 - b^2)*d*e*n*x^n*(b + a*cosh(c + d*x^n)))]
    end

max_stars_repo_path: test/6 Hyperbolic functions/6.5 Hyperbolic secant/6.5.2 (e x)^m (a+b sech(c+d x^n))^p.jl


text:
# Export DynamicLinks API
export DynamicLinks,
       getindex, length, iterate, keys, lastindex, eltype,
       handle, header

# Export Dynamic Link API
export DynamicLink,
       DynamicLinks, handle, path

# Export RPath API
export RPath,
       handle, rpaths, canonical_rpaths, find_library

# Import iteration protocol
import Base: iterate, keys, length, lastindex

"""
    DynamicLinks

This type encapsulates the list of dynamic links within an object, holding a
collection of `DynamicLink` objects.  The list of available API operations is
given below, with methods that subclasses must implement marked in emphasis:

### Creation
  - *DynamicLinks()*

### Iteration
  - *getindex()*
  - *lastindex()*
  - iterate()
  - keys()
  - eltype()

### Misc.
  - *handle()*
"""
abstract type DynamicLinks{H <: ObjectHandle} end

@mustimplement lastindex(dls::DynamicLinks)
keys(dls::DynamicLinks) = 1:length(dls)
iterate(dls::DynamicLinks, idx=1) = idx > length(dls) ? nothing : (dls[idx], idx+1)
length(dls::DynamicLinks) = lastindex(dls)
eltype(::Type{D}) where {D <: DynamicLinks} = DynamicLink
@mustimplement getindex(dls::DynamicLinks, idx)

"""
    DynamicLink

This type encapsulates the linkage of one object file to another.  The list of
available API operations is given below, with methods that subclasses must
implement marked in emphasis:

### Creation:
  - *DynamicLink()*

### Accessors:
  - *DynamicLinks()*
  - *handle()*
  - *path()*
"""
abstract type DynamicLink{H <: ObjectHandle} end

@mustimplement path(dl::DynamicLink)

"""
    RPath

This type encapsulates the search path used by an object file when looking for
a shared library.  This class enables not only looking at the path, but
querying the path for matches for given library names.  The list of available
API operations is given below, with methods that subclasses must implement
marked in emphasis:

### Creation:
  - *RPath()*

### Utility
  - *handle()*

### RPath operations
  - *rpaths()*
  - canonical_rpaths()
  - find_library()
  - lastindex()
  - iterate()
"""
abstract type RPath{H <: ObjectHandle} end

"""
    RPath(oh::ObjectHandle)

Construct an `RPath` object from the given `ObjectHandle`.
"""
@mustimplement RPath(oh::ObjectHandle)

"""
    handle(rpath::RPath)

Return the handle that this `RPath` object refers to.
"""
@mustimplement handle(rpath::RPath)

"""
    rpaths(rpath::RPath)

Return the list of paths that will be searched for shared libraries.
"""
@mustimplement rpaths(rpath::RPath)

lastindex(rpath::RPath) = lastindex(rpaths(rpath))
keys(rpath::RPath) = 1:length(rpath)
iterate(rpath::RPath, idx=1) = idx > length(rpath) ? nothing : (rpath[idx], idx+1)
length(rpath::RPath) = lastindex(rpath)
eltype(::Type{D}) where {D <: RPath} = String
getindex(rpath::RPath, idx) = rpaths(rpath)[idx]

"""
    canonical_rpaths(rpath::RPath)

Return a canonicalized list of paths that will be searched.
"""
function canonical_rpaths(rpath::RPath)
    origin = dirname(path(handle(rpath)))
    paths = rpaths(rpath)
    for idx in 1:length(paths)
        # Substitute the path of the containing handle for `$ORIGIN` and
        # `@loader_path`.  Do the same for `@executable_path` even though
        # that's technically incorrect, because we don't have a good way to
        # track the web of dependencies right now.
        paths[idx] = replace(paths[idx], "\$ORIGIN" => origin)
        paths[idx] = replace(paths[idx], "@loader_path" => origin)
        paths[idx] = replace(paths[idx], "@executable_path" => origin)
        paths[idx] = abspath(paths[idx])

        # Make sure that if it's a directory, it _always_ has a trailing path separator.
        if isdir(paths[idx]) && paths[idx][end:end] != Base.Filesystem.path_separator
            paths[idx] = paths[idx] * Base.Filesystem.path_separator
        end
    end
    return paths
end

"""
    find_library(rpath::RPath, soname::String)

Return the full path to a library, searching the given `RPath`, and then the
default library search paths.  This method takes the given `soname` and joins
it to the end of every path within the given `RPath`, returning the resultant
path if it exists, returning back the original `soname` if it doesn't.
"""
function find_library(rpath::RPath, soname::AbstractString)
    for path in canonical_rpaths(rpath)
        libpath = joinpath(path, soname)
        if isfile(libpath)
            return libpath
        end
    end
    return soname
end

### Printing
function show(io::IO, dl::DynamicLink{H}) where {H <: ObjectHandle}
    print(io, "$(format_string(H)) DynamicLink \"$(path(dl))\"")
end

show(io::IO, dls::DynamicLinks{H}) where {H <: ObjectHandle} = show_collection(io, dls, H)
show(io::IO, rp::RPath{H}) where {H <: ObjectHandle} = show_collection(io, rp, H)


text:
# This uses the `compute_forces(term, ψ, occ; kwargs...)` function defined by all terms
"""
Compute the forces of an obtained SCF solution. Returns the forces wrt. the fractional
lattice vectors. To get cartesian forces use [`compute_forces_cart`](@ref).
Returns a list of lists of forces
`[[force for atom in positions] for (element, positions) in atoms]`
which has the same structure as the `atoms` object passed to the underlying [`Model`](@ref).
"""
@timing function compute_forces(basis::PlaneWaveBasis, ψ, occ; kwargs...)
    # TODO optimize allocs here
    T = eltype(basis)
    forces = [zeros(Vec3{T}, length(positions)) for (element, positions) in basis.model.atoms]
    for term in basis.terms
        f_term = compute_forces(term, basis, ψ, occ; kwargs...)
        if !isnothing(f_term)
            forces += f_term
        end
    end
    forces
end

"""
Compute the cartesian forces of an obtained SCF solution in Hartree / Bohr.
Returns a list of lists of forces
`[[force for atom in positions] for (element, positions) in atoms]`
which has the same structure as the `atoms` object passed to the underlying [`Model`](@ref).
"""
function compute_forces_cart(basis::PlaneWaveBasis, ψ, occ; kwargs...)
    lattice = basis.model.lattice
    forces = compute_forces(basis::PlaneWaveBasis, ψ, occ; kwargs...)
    [[lattice \ f for f in forces_for_element] for forces_for_element in forces]
end

function compute_forces(scfres)
    compute_forces(scfres.basis, scfres.ψ, scfres.occupation; ρ=scfres.ρ)
end
function compute_forces_cart(scfres)
    compute_forces_cart(scfres.basis, scfres.ψ, scfres.occupation; ρ=scfres.ρ)
end


text:
using Blosc, Dates, DSP, HDF5, Logging, Printf, SeisIO, Statistics, Test
using SeisIO.FastIO, SeisIO.Nodal, SeisIO.Quake, SeisIO.RandSeis, SeisIO.SEED, SeisIO.SeisHDF, SeisIO.SUDS, SeisIO.UW

import Dates: DateTime, Hour, now
import DelimitedFiles: readdlm
import Random: rand, randperm, randstring
import SeisIO: BUF,
  KW,
  FDSN_sta_xml,
  TimeSpec,
  auto_coords,
  bad_chars,
  buf_to_double,
  checkbuf!,
  checkbuf_8!,
  code2resptyp,
  code2typ,
  datafields,
  datareq_summ,
  diff_x!,
  endtime,
  fillx_i16_le!,
  fillx_i32_be!,
  fillx_i32_le!,
  findhex,
  formats,
  get_del_ranges,
  get_http_post,
  get_http_req,
  get_views,
  ibmfloat,
  int2tstr,
  int_x!,
  mean,
  minreq,
  minreq!,
  mk_t!,
  mk_t,
  mktaper!,
  mktime,
  nx_max,
  parse_charr,
  parse_chstr,
  parse_sl,
  poly,
  polyfit,
  polyval,
  read_sacpz!,
  read_sacpz,
  read_seed_resp!,
  read_station_xml!,
  read_station_xml,
  read_sxml,
  resptyp2code,
  safe_isdir,
  safe_isfile,
  sort_segs!,
  sort_segs,
  starttime,
  sync_t,
  sμ,
  t_arr!,
  t_bounds,
  t_collapse,
  t_expand,
  t_extend,
  t_win,
  taper_seg!,
  tnote,
  trid,
  tstr2int,
  typ2code,
  w_time,
  webhdr,
  x_inds,
  xtmerge!,
  μs
import SeisIO.RandSeis: iccodes_and_units, rand_misc
import SeisIO.Quake: unsafe_convert
import SeisIO.SeisHDF:read_asdf, read_asdf!, id_match, id_to_regex
import Statistics: mean


text:
export →

## Space Operator is used to wrap other operators
# and change the domain/range space
struct SpaceOperator{O<:Operator,S<:Space,V<:Space,T} <: Operator{T}
    op::O
    domainspace::S
    rangespace::V
end

# The promote_type is needed to fix a bug in promotetimes
# not sure if its the right long term solution
SpaceOperator(o::Operator,s::Space,rs::Space) =
    SpaceOperator{typeof(o),typeof(s),typeof(rs),eltype(o)}(o,s,rs)
SpaceOperator(o,s) = SpaceOperator(o,s,s)

function convert(::Type{Operator{T}},S::SpaceOperator) where T
    if T==eltype(S)
        S
    else
        op=convert(Operator{T},S.op)
        SpaceOperator{typeof(op),typeof(S.domainspace),typeof(S.rangespace),T}(op,S.domainspace,S.rangespace)
    end
end

# Similar to wrapper, but different domain/domainspace/rangespace

@wrappergetindex SpaceOperator

# SpaceOperator can change blocks, so we need to override this
getindex(A::SpaceOperator,KR::BlockRange, JR::BlockRange) = defaultgetindex(A,KR,JR)

getindex(A::SpaceOperator,K::Block,J::Block) = A[blockrows(A,K),blockcols(A,J)]
getindex(A::SpaceOperator,K::Block,j) = A[blockrows(A,K),j]
getindex(A::SpaceOperator,k,J::Block) = A[k,blockcols(A,J)]

domain(S::SpaceOperator) = domain(domainspace(S))
domainspace(S::SpaceOperator) = S.domainspace
rangespace(S::SpaceOperator) = S.rangespace

##TODO: Do we need both max and min?
function findmindomainspace(ops::AbstractVector)::Any
    sp = UnsetSpace()

    for op in ops
        sp = union(sp,domainspace(op))
    end

    sp
end

function findmaxrangespace(ops::AbstractVector)::Any
    sp = UnsetSpace()

    for op in ops
        sp = maxspace(sp,rangespace(op))
    end

    sp
end

# The coolest definitions ever!!
# supports Derivative():Chebyshev()→Ultraspherical(1)
(:)(A::Operator,b::Space) = promotedomainspace(A,b)
→(A::Operator,b::Space) = promoterangespace(A,b)
(:)(A::UniformScaling,b::Space) = Operator(A,b)
→(A::UniformScaling,b::Space) = Operator(A,b)

promoterangespace(P::Operator,sp::Space) = promoterangespace(P,sp,rangespace(P))
promotedomainspace(P::Operator,sp::Space) = promotedomainspace(P,sp,domainspace(P))

promoterangespace(P::Operator,sp::Space,cursp::Space)::Any =
    (sp==cursp) ? P : Conversion(cursp,sp)*P
promotedomainspace(P::Operator,sp::Space,cursp::Space)::Any =
    (sp==cursp) ? P : P*Conversion(sp,cursp)

function promoterangespace(ops::AbstractVector{O}) where O<:Operator
    isempty(ops) && return ops
    k=findmaxrangespace(ops)
    #TODO: T might be incorrect
    T=mapreduce(eltype,promote_type,ops)
    Operator{T}[promoterangespace(op,k) for op in ops]
end
function promotedomainspace(ops::AbstractVector{O}) where O<:Operator
    isempty(ops) && return ops
    k=findmindomainspace(ops)
    #TODO: T might be incorrect
    T=mapreduce(eltype,promote_type,ops)
    Operator{T}[promotedomainspace(op,k) for op in ops]
end
function promotedomainspace(ops::AbstractVector{O},S::Space) where O<:Operator
    isempty(ops) && return ops
    k=conversion_type(findmindomainspace(ops),S)
    #TODO: T might be incorrect
    T=promote_type(mapreduce(eltype,promote_type,ops),prectype(S))
    Operator{T}[promotedomainspace(op,k) for op in ops]
end

####
# choosedomainspace returns a potental domainspace
# where the second argument is a target rangespace
# it defaults to the true domainspace, but if this is ambiguous
# it tries to decide a space.
###

function default_choosedomainspace(A::Operator,sp::Space)
    sp2=domainspace(A)
    isambiguous(sp2) ? sp : sp2
end

choosedomainspace(A::Operator,sp::Space) = default_choosedomainspace(A,sp)

choosedomainspace(A::Operator,f::Fun) = choosedomainspace(A,space(f))
choosedomainspace(A::Operator,f::AbstractVector{FF}) where {FF<:Fun} =
    choosedomainspace(A,Fun(f))
choosedomainspace(A::Operator,_) = choosedomainspace(A)

choosedomainspace(A) = choosedomainspace(A,UnsetSpace())

function choosedomainspace(ops::AbstractVector,spin)::Any
    sp = UnsetSpace()

    for op in ops
        sp = conversion_type(sp,choosedomainspace(op,spin))
    end

    sp
end

choosespaces(A::Operator,b) = promotedomainspace(A,choosedomainspace(A,b))

spacescompatible(A::Operator,B::Operator) =
    spacescompatible(domainspace(A),domainspace(B)) &&
    spacescompatible(rangespace(A),rangespace(B))

#It's important that domain space is promoted first as it might impact range space
promotespaces(ops::AbstractVector) = promoterangespace(promotedomainspace(ops))
function promotespaces(ops::AbstractVector,b::Fun)
    A=promotespaces(ops)
    if isa(rangespace(A),AmbiguousSpace)
        # try setting the domain space
        A=promoterangespace(promotedomainspace(ops,space(b)))
    end
    A,Fun(b,rangespace(A[end]))
end

function promotespaces(A::Operator,B::Operator)
    if spacescompatible(A,B)
        A,B
    else
        tuple(promotespaces([A,B])...)
    end
end

## algebra

ldiv_coefficients(A::SpaceOperator,b;kwds...) =
    ldiv_coefficients(A.op,b;kwds...)


text:
#### ExpandA

function _doexpand(x)
    ! isa(x,Mxpr) && return x
    n = length(x)
    args = margs(x)
    nargs = newargs(n)
    @inbounds for i in 1:n
        nargs[i] = doexpand(args[i])
    end
    return doexpand(mxpr(mhead(x),nargs))
end

doexpand(p::Mxpr{:Power}) = do_expand_power(p,base(p),expt(p))

function doexpand(prod::Mxpr{:Times})
    a = margs(prod)
    len = length(a)
    @inbounds for i in 1:len
        a[i] = doexpand(a[i])
    end
    have_sum = false
    j = 0
    @inbounds for i in 1:len # check if we have anything to do
        j += 1
        if is_Mxpr(a[i],:Plus)
            have_sum = true
            break
        end
    end
    ! have_sum && return prod
    nonsums = newargs()
    @inbounds for i in 1:j-1  # none of these are sums
        push!(nonsums,a[i])
    end
    sums = newargs()
    push!(sums,a[j]) # already know its a sum
    @inbounds for i in j+1:len   # push more sums, if there are any
        if is_Mxpr(a[i],:Plus)
            push!(sums,a[i])
        else
            push!(nonsums,a[i])
        end
    end
    sumres = length(sums) == 1 ? sums[1] : mulfacs(sums...)
    nlen = length(nonsums)
    if nlen == 0
        mxout = sumres
    elseif nlen == 1
        mxout = mulfacs(nonsums[1],sumres)
    else
        mxout = mulfacs(mxpr(:Times,nonsums),sumres)
    end
    setfixed(mxout)
    return mxout
end

doexpand(mx) = mx

do_expand_power(p,b::Mxpr{:Plus}, n::Integer) =
    length(b) != 2 ? p : do_expand_binomial(p,b[1],b[2],n)
do_expand_power(p,b,ex) = p

do_expand_binomial(p,a, b, n::Integer) = expand_binomial(a,b,n)
do_expand_binomial(a,b,n) = p

function doexpand(s::Mxpr{:Plus})
    args = margs(s)
    @inbounds for i in 1:length(args)
        args[i] = doexpand(args[i])
    end
    return s
end

## expand product of two sums
# a an b are the factors
# In test cases, this is fast. later, canonicalizing each term is slowest thing.
function mulfacs(a::Mxpr{:Plus},b::Mxpr{:Plus})
    terms = newargs(length(a)*length(b))
    i = 0
    for ax in a.args
        for bx in b.args
            i += 1
            t = flatcanon!(mxpr(:Times, ax, bx)) # TODO specialize for types of ax, bx
            mergesyms(t,ax)
            mergesyms(t,bx)
            setfixed(t)
           @inbounds  terms[i] = t
        end
    end
    mx = flatcanon!(mxpr(:Plus,terms))
    for t in terms
        mergesyms(mx,t)
    end
    setfixed(mx)
    mx
end

# Not the right way to do this. We need to expand each term, as well.
mulfacs(a,b,c) = mulfacs(mulfacs(a,b),c)
mulfacs(a,b,c,xs...) = mulfacs(mulfacs(mulfacs(a,b),c),xs...)

# Should probably write this
#function mulfacs(a::Mxpr{:Plus}, b::SJSym)
#end

function _pushfacs!(facs,mx::Mxpr{:Times})
    append!(facs,margs(mx))
end

function _pushfacs!(facs,b)
    push!(facs,b)
end

function mulfacs(a::Mxpr{:Plus},b)
    terms = newargs(length(a))
    i = 0
    for ax in a.args
        i += 1
#        if is_Mxpr(ax) && mxprtype(ax) == :Times
        if isa(ax, Mxpr{:Times})
            facs = copy(margs(ax))
            _pushfacs!(facs,b)  # I hope type inference optimizes this.
            terms[i] = mxpr(:Times,facs)
        else
            terms[i] = doeval(mmul(b, ax))
        end
        canonexpr!(terms[i])
        mergeargs(terms[i])
    end
    mx = mxpr(:Plus,terms)
    mergeargs(mx)
    mx
end

function mulfacs(a, b::Mxpr{:Plus})
    mulfacs(b,a)
end

## construct Power. Decide whether to canonicalize according to types of args
function canonpower(base,expt::Real)
    canonexpr!(mpow(base,expt))
end
function canonpower(base::SJSym,expt::Real)
    canonexpr!(mpow(base,expt))
end
function canonpower(base::SJSym,expt)
    mpow(base,expt)
end
function canonpower(base,expt)
    mpow(base,expt)
end

##

# optimize a bit for types
function _expand_binomial_aux1(a::SJSym,b::SJSym,n)
    mxpr(:Times,n,canonpower(a,(n-1)),b)
end
function _expand_binomial_aux1(a,b,n)
    flatcanon!(mxpr(:Times,n,canonpower(a,(n-1)),b))
end
function _expand_binomial_aux2(a::SJSym,b::SJSym,n)
    mxpr(:Times,n,a,canonpower(b,(n-1)))
end
function _expand_binomial_aux2(a,b,n)
    flatcanon!(mxpr(:Times,n,a,canonpower(b,(n-1))))
end

# Be careful to construct expression in canonical form.
# Lots of ways to go wrong.
# Assume a < b in canonical order
# This is the only place we are testing meta data in Mxpr giving which
# symbols it depends on.
function expand_binomial(a,b,n::T) where T<:Integer
    args = newargs(n+1)
    args[1] = canonpower(a,n)
    mergesyms(args[1],a)
    setfixed(args[1])
    args[n+1] =  canonpower(b,n)
    mergesyms(args[n+1],b)
    setfixed(args[n+1])
    if n == 2
        args[2] = flatcanon!(mxpr(:Times,2,a,b))  # we have to flatcanon
    else
        # TODO optimize for symbols a,b, as in general case below. No flatcanon.
        args[2] = _expand_binomial_aux1(a,b,n)
        args[n] = _expand_binomial_aux2(a,b,n)
        mergesyms(args[2],a)
        mergesyms(args[2],b)
        mergesyms(args[n],a)
        mergesyms(args[n],b)
        setfixed(args[2])
        setfixed(args[n])
        fac = n
        k = n
        l = one(n)
        expand_binomial_aux(k,l,n,fac,a,b,args)
    end
    mx = mxprcf(:Plus,args)
    mergesyms(mx,a)
    mergesyms(mx,b)
    apply_upvalues_to_args!(mx)  # takes some time
    setage(mx)
    setfixed(mx)
    mx
end

# Big increase in efficiency (> 10x) for both these types
# typealias ExpNoCanon Union{SJSym,Number} move to symataconstants.jl

# Expand((a+b*c)^n) is 10x slower than Expand((a+b)^n)
function _expand_mulpowers(fac,b1,e1,b2,e2)
    m1 = canonpower(b1,e1)  # adds 10-15% time
    m2 = canonpower(b2,e2)
    setfixed(m1)
    setfixed(m2)
    mergesyms(m1,b1)
    mergesyms(m2,b2)
#    return flatcanon!(mxpr(:Times, fac, m1, m2)) # flatcanon adds 10x time !, even if nothing is done
#    return flatcanon!(mmul(fac,mmul(m1,m2)))
    return doeval(mmul(fac,mmul(m1,m2)))    # flatcannon is enough because of nested call.
end

function _expand_mulpowers(fac,b1::T,e1,b2::V,e2) where {T<:ExpNoCanon, V<:ExpNoCanon}
    m1 = b1^e1
    m2 = b2^e2
    setfixed(m1)
    setfixed(m2)
    mergesyms(m1,b1)
    mergesyms(m2,b2)
    return mxpr(:Times, fac, m1, m2)
end

function _expand_mulpowers(fac,b1::T,e1,b2,e2) where T<:ExpNoCanon
    m1 = b1^e1
    m2 = canonpower(b2,e2)
    setfixed(m1)
    setfixed(m2)
    mergesyms(m1,b1)
    mergesyms(m2,b2)
    return flatcanon!(mxpr(:Times, fac, m1, m2))
end

function _expand_mulpowers(fac,b1,e1,b2::T,e2) where T<:ExpNoCanon
    m1 = canonpower(b1,e1)
    m2 = b2^e2
    setfixed(m1)
    setfixed(m2)
    mergesyms(m1,b1)
    mergesyms(m2,b2)
    return flatcanon!(mxpr(:Times, fac, m1, m2))
end

function expand_binomial_aux(k,l,n,fac,a,b,args)
        @inbounds for j in 2:n-2
            k = k - 1
            l = l + 1
            fac = mmul(fac,k)
            fac = div(fac,l)
            args[j+1] = _expand_mulpowers(fac,a,n-j,b,j)
            mergesyms(args[j+1],a)
            mergesyms(args[j+1],b)
            setfixed(args[j+1])
        end
end


text:
@testset "STDOUT" begin
    # make sure it compiles and executes
    # TODO: replace this with VT100 tests
    imshow(colorview(RGB, rand(3,2,2))); println()
    imshow(colorview(RGB, rand(3,2,2)), (2,2)); println()
    imshow256(colorview(RGB, rand(3,2,2))); println()
    imshow256(colorview(RGB, rand(3,2,2)), (2,2)); println()
    imshow24bit(colorview(RGB, rand(3,2,2))); println()
    imshow24bit(colorview(RGB, rand(3,2,2)), (2,2)); println()
end

@testset "imshow" begin
    @testset "non colorant" begin
        @test_throws ArgumentError imshow(rand(5,5))
        @test_throws ArgumentError imshow(sprand(5,5,.5))
    end
    @testset "lena" begin
        img = imresize(lena, 10, 10)
        io = IOBuffer()
        ensurecolor(imshow, io, img)
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "lena_big_imshow" res
    end
    @testset "ndarray" begin
        img = rgb_line_4d
        io = IOBuffer()
        ensurecolor(imshow, io, img)
        res = readlines(seek(io,0))
        @test_reference "ndarray_imshow" res
    end
end

@testset "imshow256" begin
    @testset "non colorant" begin
        @test_throws ArgumentError imshow256(rand(5,5))
        @test_throws ArgumentError imshow256(sprand(5,5,.5))
    end
    @testset "rgb line" begin
        io = IOBuffer()
        ensurecolor(imshow256, io, rgb_line)
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "rgbline_big_imshow256" res
        io = IOBuffer()
        ensurecolor(imshow256, io, rgb_line, (1, 45))
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "rgbline_small1_imshow256" res
        io = IOBuffer()
        ensurecolor(imshow256, io, rgb_line, (1, 19))
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "rgbline_small2_imshow256" res
    end
    @testset "lena" begin
        img = imresize(lena, 10, 10)
        io = IOBuffer()
        ensurecolor(imshow256, io, img)
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "lena_big_imshow256" res
        io = IOBuffer()
        ensurecolor(imshow256, io, img, (10, 20))
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "lena_small_imshow256" res
    end
    @testset "ndarray" begin
        img = rgb_line_4d
        io = IOBuffer()
        ensurecolor(imshow256, io, img)
        res = readlines(seek(io,0))
        @test_reference "ndarray_imshow256" res
    end
end

@testset "imshow24bit" begin
    @testset "non colorant" begin
        @test_throws ArgumentError imshow24bit(rand(5,5))
        @test_throws ArgumentError imshow24bit(sprand(5,5,.5))
    end
    @testset "rgb line" begin
        io = IOBuffer()
        ensurecolor(imshow24bit, io, rgb_line)
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "rgbline_big_imshow24bit" res
        io = IOBuffer()
        ensurecolor(imshow24bit, io, rgb_line, (1, 45))
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "rgbline_small1_imshow24bit" res
        io = IOBuffer()
        ensurecolor(imshow24bit, io, rgb_line, (1, 19))
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "rgbline_small2_imshow24bit" res
    end
    @testset "lena" begin
        img = imresize(lena, 10, 10)
        io = IOBuffer()
        ensurecolor(imshow24bit, io, img)
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "lena_big_imshow24bit" res
        io = IOBuffer()
        ensurecolor(imshow24bit, io, img, (10, 20))
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "lena_small_imshow24bit" res
    end
    @testset "ndarray" begin
        img = rgb_line_4d
        io = IOBuffer()
        ensurecolor(imshow24bit, io, img)
        res = readlines(seek(io,0))
        @test_reference "ndarray_imshow24bit" res
    end
end

@testset "imshow non1" begin
    @testset "lena" begin
        img = OffsetArray(imresize(lena, 10, 10), (-10,5))
        io = IOBuffer()
        ensurecolor(imshow, io, img)
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "lena_big_imshow" res
    end
    @testset "rotation" begin
        tfm = recenter(RotMatrix(-pi/4), center(lighthouse))
        lhr = ImageTransformations.warp(lighthouse, tfm)
        io = IOBuffer()
        ensurecolor(imshow, io, lhr)
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "lighthouse_rotated" res
    end
end


text:
# ----------------------- #
# Display-esque functions #
# ----------------------- #
function html_body(p::Plot)
    """
    <script>
        $(script_content(p))
     </script>
    """
end

function script_content(p::Plot)
    """
    gd = (function() {
      var WIDTH_IN_PERCENT_OF_PARENT = 100;
      var HEIGHT_IN_PERCENT_OF_PARENT = 100;
      var gd = Plotly.d3.select('body')
        .append('div').attr("id", "$(p.divid)")
        .style({
          width: WIDTH_IN_PERCENT_OF_PARENT + '%',
          'margin-left': (100 - WIDTH_IN_PERCENT_OF_PARENT) / 2 + '%',
          height: HEIGHT_IN_PERCENT_OF_PARENT + 'vh',
          'margin-top': (100 - HEIGHT_IN_PERCENT_OF_PARENT) / 2 + 'vh'
        })
        .node();
      var plot_json = $(json(p));
      var data = plot_json.data;
      var layout = plot_json.layout;
      Plotly.newPlot(gd, data, layout);
      window.onresize = function() {
        Plotly.Plots.resize(gd);
      };
      return gd;
    })();
    """
end

# just declare here so we can overload elsewhere
function savefig end

function savejson(p::Plot, fn::AbstractString)
    ext = split(fn, ".")[end]
    if ext == "json"
        open(f -> print(f, json(p)), fn, "w")
        return p
    else
        msg = "PlotlyBase can only save figures as JSON. For all other"
        msg *= " file types, please use PlotlyJS.jl"
        throw(ArgumentError(msg))
    end
end

# jupyterlab/nteract integration
Base.Multimedia.istextmime(::MIME"application/vnd.plotly.v1+json") = true
function Base.show(io::IO, ::MIME"application/vnd.plotly.v1+json", p::Plot)
    JSON.print(io, p)
end

function Base.show(io::IO, ::MIME"text/plain", p::Plot)
    println(io, """
    data: $(json(map(_describe, p.data), 2))
    layout: "$(_describe(p.layout))"
    """)
end

Base.show(io::IO, p::Plot) = show(io, MIME("text/plain"), p)


text:
abstract GradientBoosting

type GradientBoostingClassifier <: GradientBoosting
    n_clf::Int64
    learning_rate::Float64
    max_depth::Int64
    min_sample_split::Int64
    min_gain::Float64
    init_estimate::Union{String, Real, Vector}
    trees::Array{ClassificationTree}
end

type GradientBoostingRegressor <: GradientBoosting
    n_clf::Int64
    learning_rate::Float64
    max_depth::Int64
    min_sample_split::Int64
    min_gain::Float64
    init_estimate::Union{String, Real, Vector}
    trees::Array{RegressionTree}
end

function GradientBoostingRegressor(;
                      n_clf::Int64 = 10,
                      learning_rate::Float64 = 0.5,
                      min_sample_split::Int64 = 20,
                      min_gain::Float64 = 1e-4,
                      max_depth = 4,
                      init_estimate::Union{String, Real, Vector} = "nothing"
                      )
    trees = []
    for i = 1:n_clf
        push!(trees, RegressionTree(min_samples_split = min_sample_split,
            min_gain = min_gain, max_depth = max_depth))
    end
    return GradientBoostingRegressor(n_clf, learning_rate,
                        max_depth, min_sample_split,
                        min_gain, init_estimate, trees)
end

function GradientBoostingClassifier(;
                      n_clf::Int64 = 10,
                      learning_rate::Float64 = 0.5,
                      min_sample_split::Int64 = 20,
                      min_gain::Float64 = 1e-4,
                      max_depth = 4,
                      init_estimate::Union{String, Real, Vector} = "nothing"
                      )
    trees = []
    for i = 1:n_clf
        push!(trees, ClassificationTree(min_samples_split = min_sample_split,
            min_gain = min_gain, max_depth = max_depth))
    end
    return GradientBoostingClassifier(n_clf, learning_rate,
                        max_depth, min_sample_split,
                        min_gain, init_estimate, trees)
end

function train!(model::GradientBoosting, X::Matrix, y::Vector)
    if typeof(model) <: GradientBoostingClassifier
        y = one_hot(y)
    end
    n_sample, n_feature = size(X)
    y_pred = mean(y, 1)
    if typeof(model) <: GradientBoostingClassifier
        y_pred = repmat(y_pred,n_sample,1)
    else
        y_pred = y_pred * ones(n_sample)
    end
    count = 1
    for tree in model.trees
        residual = -(y - y_pred)
        train!(tree, X, residual)
        residual_pred = predict(tree,X)
        y_pred -= model.learning_rate * residual_pred
        #println("$(count)th tree trained!")
        count += 1
    end
end

function predict(model::GradientBoosting,
                 x::Matrix)
    count = 1
    y_pred = 0
    for tree in model.trees
        if count == 1
            y_pred = - model.learning_rate * predict(tree,x)
            count = count + 1
        else
            y_pred -= model.learning_rate * predict(tree, x)
        end
    end

    if typeof(model) <: GradientBoostingClassifier
        return sign(softmax(y_pred)-1.5)
    end
    return y_pred

end

function test_GradientBoostingRegressor()
    X_train, X_test, y_train, y_test = make_reg(n_features = 1)
    model = GradientBoostingRegressor()
    train!(model,X_train, y_train)
    predictions = predict(model,X_test)
    print("regression msea", mean_squared_error(y_test, predictions))
    PyPlot.scatter(X_test, y_test, color = "black")
    PyPlot.scatter(X_test, predictions, color = "green")
    legend(loc="upper right",fancybox="true")
end

function test_GradientBoostingClassifier()
    X_train, X_test, y_train, y_test = make_cla()
    model = GradientBoostingClassifier()
    train!(model,X_train, y_train)
    predictions = predict(model,X_test)
    println("classification accuracy: ", accuracy(y_test, predictions))
end