text:
module PolelikeCoordinate
    using NCDatasets
    using LinearAlgebra: ⋅, normalize, normalize!, norm

    include("GridFile.jl")

    """
    This function calculate the angle between the real north and the
    north defined by vertices. I call it T-north.
    Each T-grid is defined by four vertices. As shown below:
    [4] --- [3]
     |       |
     |       |
    [1] --- [2]
    T-north for this grid is defined by the average of [1]->[4] and
    [2]->[3]. Likewise, T-east is the east defined by vertices, the
    average of [1]->[2], [4]->[3]
    Positive angle α means the local T-north and T-east will align
    with true north and east if rotating true north counterclockwise
    by angle α.
    """

    d2r = π / 180.0
    r2d = 180 / π

    abstract type Grid end

    struct CurvilinearSphericalGrid <: Grid

        R     :: Float64
        Ω     :: Float64

        Nx    :: Integer
        Ny    :: Integer
        Nz    :: Integer

        ds1   :: AbstractArray{Float64, 2}
        ds2   :: AbstractArray{Float64, 2}
        ds3   :: AbstractArray{Float64, 2}
        ds4   :: AbstractArray{Float64, 2}

        α_t     :: AbstractArray{Float64, 2}
        cosα_t  :: AbstractArray{Float64, 2}
        sinα_t  :: AbstractArray{Float64, 2}

        α_uv    :: AbstractArray{Float64, 2}
        cosα_uv :: AbstractArray{Float64, 2}
        sinα_uv :: AbstractArray{Float64, 2}

        λ_T  :: AbstractArray{Float64, 3}
        λ_UV :: AbstractArray{Float64, 3}

        ϕ_T :: AbstractArray{Float64, 3}
        ϕ_UV :: AbstractArray{Float64, 3}

        z_T :: AbstractArray{Float64, 3}
        z_W :: AbstractArray{Float64, 3}
        z_UV :: AbstractArray{Float64, 3}

        Δx_T  :: AbstractArray{Float64, 3}
        Δx_U  :: AbstractArray{Float64, 3}
        Δx_V  :: AbstractArray{Float64, 3}
        Δx_W  :: AbstractArray{Float64, 3}
        Δx_UV  :: AbstractArray{Float64, 3}

        Δy_T  :: AbstractArray{Float64, 3}
        Δy_U  :: AbstractArray{Float64, 3}
        Δy_V  :: AbstractArray{Float64, 3}
        Δy_W  :: AbstractArray{Float64, 3}
        Δy_UV  :: AbstractArray{Float64, 3}

        Δz_T  :: AbstractArray{Float64, 3}
        Δz_U  :: AbstractArray{Float64, 3}
        Δz_V  :: AbstractArray{Float64, 3}
        Δz_W  :: AbstractArray{Float64, 3}
        Δz_UV  :: AbstractArray{Float64, 3}

        function CurvilinearSphericalGrid(;
            R      :: Float64,
            Ω      :: Float64,
            Nx     :: Integer,
            Ny     :: Integer,
            λ_t    :: AbstractArray{Float64, 2},  # center longitude
            ϕ_t    :: AbstractArray{Float64, 2},  # center latitude
            vs_λ   :: AbstractArray{Float64, 3},  # vertices longitude (4, Nx, Ny)
            vs_ϕ   :: AbstractArray{Float64, 3},  # vertices latitude  (4, Nx, Ny)
            z_w    :: AbstractArray{Float64, 1},
            angle_unit :: Symbol = :deg,
            sub_yrng :: Union{Colon, UnitRange} = Colon(),
        )

            Nz = length(z_w) - 1

            α_t    = zeros(Float64, Nx, Ny)
            cosα_t = zeros(Float64, Nx, Ny)
            sinα_t = zeros(Float64, Nx, Ny)

            α_uv    = zeros(Float64, Nx, Ny+1)
            cosα_uv = zeros(Float64, Nx, Ny+1)
            sinα_uv = zeros(Float64, Nx, Ny+1)

            ds1 = zeros(Float64, Nx, Ny)
            ds2 = zeros(Float64, Nx, Ny)
            ds3 = zeros(Float64, Nx, Ny)
            ds4 = zeros(Float64, Nx, Ny)

            ps = zeros(Float64, 3, 4)

            λ_t = copy(λ_t)
            ϕ_t = copy(ϕ_t)

            vs_λ = copy(vs_λ)
            vs_ϕ = copy(vs_ϕ)

            if angle_unit == :deg

                λ_t .*= d2r
                ϕ_t .*= d2r

                vs_λ .*= d2r
                vs_ϕ .*= d2r

            elseif angle_unit == :rad

                # do nothing

            else
                throw(ErrorException("Unknown `angle_unit`: " * angle_unit))
            end

            println("Compute vectors and sides")
            # First, compute the needed vectors and sides

            permuted_vs_λ = permutedims(vs_λ, [2,3,1])
            permuted_vs_ϕ = permutedims(vs_ϕ, [2,3,1])

            shp = size(permuted_vs_ϕ)
            xv = reshape( R * cos.(permuted_vs_ϕ) .* cos.(permuted_vs_λ) , 1, shp...)
            yv = reshape( R * cos.(permuted_vs_ϕ) .* sin.(permuted_vs_λ) , 1, shp...)
            zv = reshape( R * sin.(permuted_vs_ϕ)                        , 1, shp...)

            ps = vcat(xv, yv, zv)
            u1 = ps[:, :, :, 2] - ps[:, :, :, 1]
            u2 = ps[:, :, :, 3] - ps[:, :, :, 2]
            u3 = ps[:, :, :, 4] - ps[:, :, :, 3]
            u4 = ps[:, :, :, 1] - ps[:, :, :, 4]

            ds1 = mapslices(norm, u1; dims=[1,])[1, :, :]
            ds2 = mapslices(norm, u2; dims=[1,])[1, :, :]
            ds3 = mapslices(norm, u3; dims=[1,])[1, :, :]
            ds4 = mapslices(norm, u4; dims=[1,])[1, :, :]

            println("size of ds1: ", size(ds1))

            #=
            @time for j = 1:Ny
            for i = 1:Nx

                for k = 1:4
                    _λ = vs_λ[k, i, j]
                    _ϕ = vs_ϕ[k, i, j]

                    ps[1, k] = R * cos(_ϕ) * cos(_λ)
                    ps[2, k] = R * cos(_ϕ) * sin(_λ)
                    ps[3, k] = R * sin(_ϕ)
                end

                u1[:, i, j] = ps[:, 2] - ps[:, 1]
                u2[:, i, j] = ps[:, 3] - ps[:, 2]
                u3[:, i, j] = ps[:, 4] - ps[:, 3]
                u4[:, i, j] = ps[:, 1] - ps[:, 4]

                ds1[i, j] = norm(u1)
                ds2[i, j] = norm(u2)
                ds3[i, j] = norm(u3)
                ds4[i, j] = norm(u4)
            end
            end
            =#
            # Then, compute the angle α
            function computeα(
                _λ :: Float64,
                _ϕ :: Float64,
                grid_north :: Array{Float64, 1},
            )

                true_north  = [ - sin(_ϕ) * cos(_λ), - sin(_ϕ) * sin(_λ),   cos(_ϕ) ]
                true_east   = [ - sin(_λ)          ,   cos(_λ)          ,   0.0     ]
                # true_upward = [   cos(_ϕ) * cos(_λ),   cos(_ϕ) * sin(_λ),   sin(_ϕ) ]

                # α is the same "ANGLE" variable in POP2 output.
                # It is the angle between true east and grid east (true north and grid north)
                # positive means the grid is rotated counterclockwise looking top-down
                cos_α =   (grid_north ⋅ true_north) / norm(grid_north)
                sin_α = - (grid_north ⋅ true_east)  / norm(grid_north)

                # Prevent some numerical errors such that
                # cos_α is not within [-1, 1]
                if abs(cos_α) > 1.0
                    if abs(cos_α) - 1 < 1e-3   # some tolerance
                        cos_α = 1.0 * sign(cos_α)
                    else
                        throw(ErrorException("Angle calculation goes wrong. Please check."))
                    end
                end

                _α = acos(cos_α) # acos only gives number within [0, π]

                if sin_α < 0.0 # This means α is in [π, 2π]
                    _α = 2*π - _α
                end

                return _α, cos(_α), sin(_α)
            end

            println("Compute grid curved angle on T grid")
            for i=1:Nx, j=1:Ny

                α_t[i, j], cosα_t[i, j], sinα_t[i, j] = computeα(
                    λ_t[i, j],
                    ϕ_t[i, j],
                    u2[:, i, j] + (- u4[:, i, j]),
                )

            end

            println("Compute grid curved angle on UV grid")
            for i=1:Nx, j=1:Ny+1

                if j == 1
                    _λ = vs_λ[1, i, 1]
                    _ϕ = vs_ϕ[1, i, 1]
                    grid_north = - u4[:, i, 1]

                elseif j == Ny+1
                    _λ = vs_λ[4, i, Ny]
                    _ϕ = vs_ϕ[4, i, Ny]
                    grid_north = - u4[:, i, Ny]

                else
                    _λ = vs_λ[1, i, j]
                    _ϕ = vs_ϕ[1, i, j]
                    grid_north = - u4[:, i, j-1] - u4[:, i, j]
                end

                α_uv[i, j], cosα_uv[i, j], sinα_uv[i, j] = computeα(
                    _λ,
                    _ϕ,
                    grid_north,
                )

            end

            cvt13 = (arr, nx, ny) -> repeat(reshape(arr, :, 1, 1), outer=(1, nx, ny))
            cvt23 = (arr, nz) -> repeat( reshape(arr, 1, size(arr)...), outer=(nz, 1, 1) )

            λ_T  = cvt23(λ_t, Nz)
            λ_UV = zeros(Float64, Nx, Ny+1)
            λ_UV[:, 1:Ny] = vs_λ[1, :, :]
            λ_UV[:, end ] = vs_λ[4, :, end]
            λ_UV = cvt23(λ_UV, Nz)

            ϕ_T = cvt23(ϕ_t, Nz)
            ϕ_UV = zeros(Float64, Nx, Ny+1)
            ϕ_UV[:, 1:Ny] = vs_ϕ[1, :, :]
            ϕ_UV[:, end] = vs_ϕ[4, :, end]
            ϕ_UV = cvt23(ϕ_UV, Nz)

            z_t = ( z_w[1:end-1] + z_w[2:end] ) / 2.0
            z_T = cvt13(z_t, Nx, Ny)
            z_W = cvt13(z_w, Nx, Ny)
            z_UV = cvt13(z_t, Nx, Ny+1)

            println("Nz, Nx, Ny = ", (Nz, Nx, Ny))

            Δx_t = (ds1 + ds3) / 2.0
            Δx_T = cvt23(Δx_t, Nz)
            Δx_U = (Δx_T + circshift(Δx_T, (0, 1, 0)))/2.0
            Δx_v = zeros(Float64, Nx, Ny+1)
            Δx_v[:, 1:Ny] = ds1
            Δx_v[:, Ny+1] = ds4[:, end]
            Δx_V = cvt23(Δx_v, Nz)
            Δx_W = cvt23(Δx_t, Nz+1)
            Δx_uv = zeros(Float64, Nx, Ny+1)
            Δx_uv[:, 1:Ny] = (ds1 + circshift(ds1, (1, 0))) / 2.0
            Δx_uv[:, Ny+1] = (ds3[:, Ny] + circshift(ds3, (1, 0))[:, Ny]) / 2.0
            Δx_UV = cvt23(Δx_uv, Nz)

            Δy_t = (ds2 + ds4) / 2.0
            Δy_T = cvt23(Δy_t, Nz)
            Δy_U = cvt23(ds4, Nz)
            Δy_v = zeros(Float64, Nx, Ny+1)
            Δy_v[:, 2:Ny] = (Δy_t[:, 1:Ny-1] + Δy_t[:, 2:Ny]) / 2.0
            Δy_v[:, 1]    = Δy_v[:, 2]
            Δy_v[:, Ny+1] = Δy_v[:, Ny]
            Δy_V = cvt23(Δy_v, Nz)
            Δy_W = cvt23(Δy_t, Nz+1)
            Δy_uv = zeros(Float64, Nx, Ny+1)
            Δy_uv[:, 2:Ny] = (ds4[:, 1:Ny-1] + ds4[:, 2:Ny]) / 2.0
            Δy_uv[:, 1] = Δy_uv[:, 2]
            Δy_uv[:, Ny+1] = Δy_uv[:, Ny]
            Δy_UV = cvt23(Δy_uv, Nz)

            Δz_t = z_w[1:end-1] - z_w[2:end]
            Δz_T = cvt13(Δz_t, Nx, Ny)
            Δz_U = cvt13(Δz_t, Nx, Ny)
            Δz_V = cvt13(Δz_t, Nx, Ny+1)
            Δz_w = zeros(Float64, Nz+1)
            Δz_w[2:Nz] = (Δz_t[1:end-1] + Δz_t[2:end]) / 2.0
            Δz_w[1] = Δz_w[2]
            Δz_w[end] = Δz_w[end-1]
            Δz_W  = cvt13(Δz_w, Nx, Ny)
            Δz_UV = cvt13(Δz_t, Nx, Ny+1)

            f_T = 2Ω * sin.(ϕ_T)

            if sub_yrng == Colon()
                sub_yrng = 1:Ny
            end

            new_Ny = length(sub_yrng)
            sub_yrng_ext = sub_yrng[1]:sub_yrng[end]+1

            return new(
                R,
                Ω,
                Nx,
                new_Ny,
                Nz,
                ds1[:, sub_yrng],
                ds2[:, sub_yrng],
                ds3[:, sub_yrng],
                ds4[:, sub_yrng],
                α_t[:, sub_yrng],
                cosα_t[:, sub_yrng],
                sinα_t[:, sub_yrng],
                α_uv[:, sub_yrng_ext],
                cosα_uv[:, sub_yrng_ext],
                sinα_uv[:, sub_yrng_ext],
                λ_T[:, :, sub_yrng],
                λ_UV[:, :, sub_yrng_ext],
                ϕ_T[:, :, sub_yrng],
                ϕ_UV[:, :, sub_yrng_ext],
                z_T[:, :, sub_yrng],
                z_W[:, :, sub_yrng],
                z_UV[:, :, sub_yrng_ext],

                Δx_T[:, :, sub_yrng],
                Δx_U[:, :, sub_yrng],
                Δx_V[:, :, sub_yrng_ext],
                Δx_W[:, :, sub_yrng],
                Δx_UV[:, :, sub_yrng_ext],

                Δy_T[:, :, sub_yrng],
                Δy_U[:, :, sub_yrng],
                Δy_V[:, :, sub_yrng_ext],
                Δy_W[:, :, sub_yrng],
                Δy_UV[:, :, sub_yrng_ext],

                Δz_T[:, :, sub_yrng],
                Δz_U[:, :, sub_yrng],
                Δz_V[:, :, sub_yrng_ext],
                Δz_W[:, :, sub_yrng],
                Δz_UV[:, :, sub_yrng_ext],
            )

        end
    end

    function project(
        gi    :: CurvilinearSphericalGrid,
        ivf_e :: AbstractArray{Float64, 2},     # input vector field east
        ivf_n :: AbstractArray{Float64, 2};     # input vector field north
        grid      :: Symbol,
        direction :: Symbol = :Forward,
    )
        ovf_e = similar(ivf_e)
        ovf_n = similar(ivf_n)

        project!(
            gi,
            ivf_e, ivf_n,
            ovf_e, ovf_n;
            direction=direction,
            grid=grid,
        )

        return ovf_e, ovf_n
    end

    function project!(
        gi    :: CurvilinearSphericalGrid,
        ivf_e :: AbstractArray{Float64, 2},     # input vector field east
        ivf_n :: AbstractArray{Float64, 2},     # input vector field north
        ovf_e :: AbstractArray{Float64, 2},    # output vector field east
        ovf_n :: AbstractArray{Float64, 2};    # output vector field north
        grid      :: Symbol,
        direction :: Symbol = :Forward,
    )

        if grid == :T
            COSα = gi.cosα_t
            SINα = gi.sinα_t
        elseif grid == :UV
            COSα = gi.cosα_uv
            SINα = gi.sinα_uv
        else
            throw(ErrorException("Unrecognized grid: " * string(grid)))
        end
        if direction == :Forward   # from outside world onto dispalced pole grid
            @. ovf_e =   ivf_e * COSα + ivf_n * SINα
            @. ovf_n = - ivf_e * SINα + ivf_n * COSα
            #for i=1:gi.Nx, j=1:gi.Ny
            #    ovf_e[i, j] =   ivf_e[i, j] * gi.cosα[i, j] - ivf_n[i, j] * gi.sinα[i, j]
            #    ovf_n[i, j] =   ivf_e[i, j] * gi.sinα[i, j] + ivf_n[i, j] * gi.cosα[i, j]
            #end
        elseif direction == :Backward   # from displaced pole grid onto outside world
            @. ovf_e = ivf_e * COSα - ivf_n * SINα
            @. ovf_n = ivf_e * SINα + ivf_n * COSα
            #for i=1:gi.Nx, j=1:gi.Ny
            #    ovf_e[i, j] =   ivf_e[i, j] * gi.cosα[i, j] + ivf_n[i, j] * gi.sinα[i, j]
            #    ovf_n[i, j] = - ivf_e[i, j] * gi.sinα[i, j] + ivf_n[i, j] * gi.cosα[i, j]
            #end
        else
            throw(ErrorException("Unsupported direction: " * string(direction)))
        end

    end

    function genGrid(
        gf       :: GridFile,
        z_w      :: AbstractArray{Float64, 1};
        sub_yrng :: Union{Colon, UnitRange} = Colon(),
    )

        local gi

        if typeof(gf) <: CurvilinearSphericalGridFile

            gi = CurvilinearSphericalGrid(;
                R=gf.R,
                Ω=gf.Ω,
                Nx=gf.Nx,
                Ny=gf.Ny,
                λ_t=gf.xc,
                ϕ_t=gf.yc,
                vs_λ=gf.xv,
                vs_ϕ=gf.yv,
                z_w=z_w,
                angle_unit=:deg,
                sub_yrng=sub_yrng,
            )

        end

        return gi
    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
module HDFS

using Compat
using URIParser

import  Base: pwd, readdir, isdir, cd, mkdir, mv, cp, rm, rmdir,
        open, close, eof, read, write, readbytes, peek,
        readall, flush, nb_available, position, stat, filesize,
        seek, seekend, seekstart, skip,
        isequal, ==, hash

export  hdfs_connect,
        hdfs_exists, hdfs_delete,
        hdfs_flush, hdfs_get_capacity, hdfs_get_default_block_size, hdfs_get_path_info, hdfs_get_used, hdfs_is_directory,
        hdfs_cd, hdfs_pwd, hdfs_list_directory, hdfs_mkdir,
        hdfs_open, hdfs_close, hdfs_rename, hdfs_copy, hdfs_move,
        hdfs_pread, hdfs_read, hdfs_seek, hdfs_tell, hdfs_write, hdfs_available,
        hdfs_set_replication, hdfs_blocks,
        isequal, hash,
        # IO methods
        pwd, readdir, isdir, cd, mkdir, mv, cp, rm, rmdir, open, close, eof, read, write, readbytes, peek,
        readall, flush, nb_available, position, stat, filesize, seek, seekend, seekstart, skip,
        # from hdfs_types.jl
        HDFS_OBJ_FILE, HDFS_OBJ_DIR, HDFS_OBJ_INVALID,
        HdfsFS, HdfsFile, HdfsFileInfo, HdfsURL

global _debug = false
function _set_debug(d)
    global _debug
    _debug = d
end

const _libhdfs = "libhdfs"
if isless(Base.VERSION, v"0.4.0-")
dlopen(_libhdfs)
else
Libdl.dlopen(_libhdfs)
end

include("hdfs_types.jl")
include("hdfs_dfs.jl")

end

ext: jl


text:
using PartialRejectionSampling
using Random

n = 100
alphabet = ["A", "C", "G", "T"]
pattern = "ACAC"
pfs = PRS.PatternFreeString(pattern, alphabet)

rng = Random.MersenneTwister(123)
@time PRS.generate_sample_prs(rng, pfs, n)


text:
struct ODESolution{T,N,uType,uType2,DType,tType,rateType,P,A,IType} <: AbstractODESolution{T,N}
  u::uType
  u_analytic::uType2
  errors::DType
  t::tType
  k::rateType
  prob::P
  alg::A
  interp::IType
  dense::Bool
  tslocation::Int
  destats::DEStats
  retcode::Symbol
end
(sol::ODESolution)(t,deriv::Type=Val{0};idxs=nothing,continuity=:left) = sol.interp(t,idxs,deriv,sol.prob.p,continuity)
(sol::ODESolution)(v,t,deriv::Type=Val{0};idxs=nothing,continuity=:left) = sol.interp(v,t,idxs,deriv,sol.prob.p,continuity)

function build_solution(
        prob::Union{AbstractODEProblem,AbstractDDEProblem},
        alg,t,u;timeseries_errors=length(u)>2,
        dense=false,dense_errors=dense,
        calculate_error = true,
        k=[],
        interp = LinearInterpolation(t,u),
        retcode = :Default, destats=DEStats(), kwargs...)

  T = eltype(eltype(u))
  if typeof(prob.u0) <: Tuple
    N = length((size(ArrayPartition(prob.u0))..., length(u)))
  else
    N = length((size(prob.u0)..., length(u)))
  end

  if typeof(prob.f) <: Tuple
    f = prob.f[1]
  else
    f = prob.f
  end

  if has_analytic(f)
    if !(typeof(prob.u0) <: Tuple)
      u_analytic = Vector{typeof(prob.u0)}()
      errors = Dict{Symbol,real(eltype(prob.u0))}()
    else
      u_analytic = Vector{typeof(ArrayPartition(prob.u0))}()
      errors = Dict{Symbol,real(eltype(prob.u0[1]))}()
    end

    sol = ODESolution{T,N,typeof(u),typeof(u_analytic),typeof(errors),typeof(t),typeof(k),
                      typeof(prob),typeof(alg),typeof(interp)}(u,u_analytic,
                       errors,t,k,prob,alg,interp,dense,0,destats,retcode)
    if calculate_error
      calculate_solution_errors!(sol;timeseries_errors=timeseries_errors,dense_errors=dense_errors)
    end
    return sol
  else
    return ODESolution{T,N,typeof(u),Nothing,Nothing,typeof(t),typeof(k),
                       typeof(prob),typeof(alg),typeof(interp)}(u,nothing,nothing,
                       t,k,prob,alg,interp,dense,0,destats,retcode)
  end
end

function calculate_solution_errors!(sol::AbstractODESolution;fill_uanalytic=true,timeseries_errors=true,dense_errors=true)

  f = sol.prob.f

  if fill_uanalytic
    for i in 1:size(sol.u,1)
      push!(sol.u_analytic,f.analytic(sol.prob.u0,sol.prob.p,sol.t[i]))
    end
  end

  save_everystep = length(sol.u) > 2
  if !isempty(sol.u_analytic)
    sol.errors[:final] = norm(recursive_mean(abs.(sol.u[end].-sol.u_analytic[end])))

    if save_everystep && timeseries_errors
      sol.errors[:l∞] = norm(maximum(vecvecapply((x)->abs.(x),sol.u-sol.u_analytic)))
      sol.errors[:l2] = norm(sqrt(recursive_mean(vecvecapply((x)->float.(x).^2,sol.u-sol.u_analytic))))
      if sol.dense && dense_errors
        densetimes = collect(range(sol.t[1], stop=sol.t[end], length=100))
        interp_u = sol(densetimes)
        interp_analytic = VectorOfArray([f.analytic(sol.prob.u0,sol.prob.p,t) for t in densetimes])
        sol.errors[:L∞] = norm(maximum(abs.(interp_u.-interp_analytic)))
        sol.errors[:L2] = norm(sqrt(mean((interp_u.-interp_analytic).^2)))
      end
    end
  end
end

function build_solution(sol::AbstractODESolution{T,N},u_analytic,errors) where {T,N}
  ODESolution{T,N,typeof(sol.u),typeof(u_analytic),typeof(errors),typeof(sol.t),typeof(sol.k),
                     typeof(sol.prob),typeof(sol.alg),typeof(sol.interp)}(
                     sol.u,u_analytic,errors,sol.t,sol.k,sol.prob,
                     sol.alg,sol.interp,sol.dense,sol.tslocation,sol.destats,sol.retcode)
end

function solution_new_retcode(sol::AbstractODESolution{T,N},retcode) where {T,N}
  ODESolution{T,N,typeof(sol.u),typeof(sol.u_analytic),typeof(sol.errors),
                     typeof(sol.t),typeof(sol.k),
                     typeof(sol.prob),typeof(sol.alg),typeof(sol.interp)}(
                     sol.u,sol.u_analytic,sol.errors,sol.t,sol.k,sol.prob,
                     sol.alg,sol.interp,sol.dense,sol.tslocation,sol.destats,retcode)
 end

 function solution_new_tslocation(sol::AbstractODESolution{T,N},tslocation) where {T,N}
   ODESolution{T,N,typeof(sol.u),typeof(sol.u_analytic),typeof(sol.errors),
                      typeof(sol.t),typeof(sol.k),
                      typeof(sol.prob),typeof(sol.alg),typeof(sol.interp)}(
                      sol.u,sol.u_analytic,sol.errors,sol.t,sol.k,sol.prob,
                      sol.alg,sol.interp,sol.dense,tslocation,sol.destats,sol.retcode)
  end

  function solution_slice(sol::AbstractODESolution{T,N},I) where {T,N}
    ODESolution{T,N,typeof(sol.u),typeof(sol.u_analytic),typeof(sol.errors),
                       typeof(sol.t),typeof(sol.k),
                       typeof(sol.prob),typeof(sol.alg),typeof(sol.interp)}(
                       sol.u[I],
                       sol.u_analytic === nothing ? nothing : sol.u_analytic[I],
                       sol.errors,sol.t[I],
                       sol.dense ? sol.k[I] : sol.k,
                       sol.prob,
                       sol.alg,sol.interp,false,sol.tslocation,sol.destats,sol.retcode)
   end


text:
sampleFreq = 512.0
twoKp = 1.0
twoKi = 0.0
q0 = 1.0
q1 = 0.0
q2 = 0.0
q3 = 0.0
integralFBx = 0.0
integralFBy = 0.0
integralFBz = 0.0

function invsqrt(x)
	return 1/sqrt(x)
end

void AHRS_update_IMU(gx, gy, gz, ax, ay, az, q0, q1, q2, q3) {

	if(!((ax == 0.0) && (ay == 0.0) && (az == 0.0))) {

		norm = invsqrt(ax * ax + ay * ay + az * az)
		ax = ax * norm
		ay = ay * norm
		az = az * norm

		halfvx = q1 * q3 - q0 * q2
		halfvy = q0 * q1 + q2 * q3
		halfvz = q0 * q0 - 0.5 + q3 * q3

		halfex = (ay * halfvz - az * halfvy)
		halfey = (az * halfvx - ax * halfvz)
		halfez = (ax * halfvy - ay * halfvx)

		if(twoKi > 0.0) {
			integralFBx = integralFBx + twoKi * halfex * (1.0 / sampleFreq)
			integralFBy = integralFBy + twoKi * halfey * (1.0 / sampleFreq)
			integralFBz = integralFBz + twoKi * halfez * (1.0 / sampleFreq)
			gx = gx + integralFBx
			gy = gy + integralFBy
			gz = gz + integralFBz
		}
		else {
			integralFBx = 0.0
			integralFBy = 0.0
			integralFBz = 0.0
		}

		gx = gx + twoKp * halfex
		gy = gy + twoKp * halfey
		gz = gz + twoKp * halfez
	}

	gx = gx * (0.5 * (1.0 / sampleFreq))
	gy = gy * (0.5 * (1.0 / sampleFreq))
	gz = gz * (0.5 * (1.0 / sampleFreq))
	qa = q0
	qb = q1
	qc = q2
	q0 = q0 + (-qb * gx - qc * gy - q3 * gz)
	q1 = q1 + (qa * gx + qc * gz - q3 * gy)
	q2 = q2 + (qa * gy - qb * gz + q3 * gx)
	q3 = q3 + (qa * gz + qb * gy - qc * gx)

	norm = invsqrt(q0 * q0 + q1 * q1 + q2 * q2 + q3 * q3)
	q0 = q0 * norm
	q1 = q1 * norm
	q2 = q2 * norm
	q3 = q3 * norm
}

void AHRS_update(gx, gy, gz, ax, ay, az, mx, my, mz, q0, q1, q2, q3) {

	if((mx == 0.0) && (my == 0.0) && (mz == 0.0)) {
		AHRS_update_IMU(gx, gy, gz, ax, ay, az, q0, q1, q2, q3)
		return
	}

	if(!((ax == 0.0) && (ay == 0.0) && (az == 0.0))) {

		norm = invsqrt(ax * ax + ay * ay + az * az)
		ax = ax * norm
		ay = ay * norm
		az = az * norm

		norm = invsqrt(mx * mx + my * my + mz * mz)
		mx = mx * norm
		my = my * norm
		mz = mz * norm

        q0q0 = q0 * q0
        q0q1 = q0 * q1
        q0q2 = q0 * q2
        q0q3 = q0 * q3
        q1q1 = q1 * q1
        q1q2 = q1 * q2
        q1q3 = q1 * q3
        q2q2 = q2 * q2
        q2q3 = q2 * q3
        q3q3 = q3 * q3

        hx = 2.0 * (mx * (0.5 - q2q2 - q3q3) + my * (q1q2 - q0q3) + mz * (q1q3 + q0q2))
        hy = 2.0 * (mx * (q1q2 + q0q3) + my * (0.5 - q1q1 - q3q3) + mz * (q2q3 - q0q1))
        bx = sqrt(hx * hx + hy * hy)
        bz = 2.0 * (mx * (q1q3 - q0q2) + my * (q2q3 + q0q1) + mz * (0.5 - q1q1 - q2q2))

		halfvx = q1q3 - q0q2
		halfvy = q0q1 + q2q3
		halfvz = q0q0 - 0.5 + q3q3
        halfwx = bx * (0.5 - q2q2 - q3q3) + bz * (q1q3 - q0q2)
        halfwy = bx * (q1q2 - q0q3) + bz * (q0q1 + q2q3)
        halfwz = bx * (q0q2 + q1q3) + bz * (0.5 - q1q1 - q2q2)

		halfex = (ay * halfvz - az * halfvy) + (my * halfwz - mz * halfwy)
		halfey = (az * halfvx - ax * halfvz) + (mz * halfwx - mx * halfwz)
		halfez = (ax * halfvy - ay * halfvx) + (mx * halfwy - my * halfwx)

		if(twoKi > 0.0) {
			integralFBx = integralFBx + twoKi * halfex * (1.0 / sampleFreq)
			integralFBy = integralFBy + twoKi * halfey * (1.0 / sampleFreq)
			integralFBz = integralFBz + twoKi * halfez * (1.0 / sampleFreq)
			gx = gx + integralFBx
			gy = gy + integralFBy
			gz = gz + integralFBz
		}
		else {
			integralFBx = 0.0
			integralFBy = 0.0
			integralFBz = 0.0
		}

		gx = gx + twoKp * halfex
		gy = gy + twoKp * halfey
		gz = gz + twoKp * halfez
	}

	gx = gx * (0.5 * (1.0 / sampleFreq))
	gy = gy * (0.5 * (1.0 / sampleFreq))
	gz = gz * (0.5 * (1.0 / sampleFreq))
	qa = q0
	qb = q1
	qc = q2
	q0 = q0 + (-qb * gx - qc * gy - q3 * gz)
	q1 = q1 + (qa * gx + qc * gz - q3 * gy)
	q2 = q2 + (qa * gy - qb * gz + q3 * gx)
	q3 = q3 + (qa * gz + qb * gy - qc * gx)

	norm = invsqrt(q0 * q0 + q1 * q1 + q2 * q2 + q3 * q3)
	q0 = q0 * norm
	q1 = q1 * norm
	q2 = q2 * norm
	q3 = q3 * norm
}


text:
"""

```
norm_sqr(v::AbstractVector{T}) where T
```

Internal function or structure - interface may change.

# Extended help

`norm_sqr` from LinearAlgebra is not documented and is slower than this for
standard arrays.

"""
@inline function norm_sqr(v::AbstractVector{T}) where T
    n2 = zero(T)^2 # the square here is required for Units, for example
    @simd for x in v
        n2 += x^2
    end
    return n2
end

"""

```
norm(v::AbstractVector{T}) where T
```

Internal function or structure - interface may change.

# Extended help

`norm_sqr` from LinearAlgebra is not documented and is slower than this for
standard arrays. Thus we define our own `norm(x) = norm_sqr(x)`.

"""
@inline norm(v) = sqrt(norm_sqr(v))

"""

```
dot(x::AbstractVector{T1},y::AbstractVector{T2}) where {T1,T2}
```

Internal function or structure - interface may change.

# Extended help

`LinearAlgebra.dot` is slower than this for standard arrays (likely more accurate, but
that is not relevant here).

"""
@inline function dot(x::AbstractVector{T1},y::AbstractVector{T2}) where {T1,T2}
    @assert length(x) == length(y)
    d = zero(T1)*zero(T2)
    @inbounds @simd for i in eachindex(x)
        d += x[i]*y[i]
    end
    return d
end


text:
module ElectrodermalActivity

using CImGui
using CImGui.CSyntax
using CImGui.CSyntax.CStatic
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using CImGui.GLFWBackend.GLFW
using CImGui.OpenGLBackend.ModernGL
using CImGui: ImVec2, ImVec4, IM_COL32, ImU32
using CImGuiExtensions
using DataFrames
using Dates
using TimeZones
using LinearAlgebra
using DSP
using CSV # Temporarily needed for prepare_heartrate fix.
using Statistics

include("electrodermal.jl")
include("tags.jl")
include("control.jl")
include("vendors.jl")
include("duration.jl")
include("timestamps.jl")
include("timezones.jl")
include("importers.jl")
include("exporters.jl")
include("menubar.jl")
include("main.jl")
include("utility.jl")

export launch,
       prepare_heartrate,
       prepare_empatica

end # module


text:
function surface_2glvisualize(kw_args)
    result = Dict{Symbol, Any}()
    xy = []

    for (k, v) in kw_args
        if k in (:x, :y)
            push!(xy, v)
            continue
        end
        k in always_skip && continue
        if k == :colornorm
            k = :color_norm
        end
        if k == :colormap
            k = :color_map
        end
        result[k] = to_signal(v)
    end
    result[:fxaa] = true
    x, y = xy
    x_is_ur, xrange = is_unitrange(to_value(x))
    y_is_ur, yrange = is_unitrange(to_value(y))
    main = to_signal(kw_args[:z])
    if x_is_ur && y_is_ur
        result[:ranges] = (x, y)
    else
        if isa(x, AbstractMatrix) && isa(y, AbstractMatrix)
            main = to_signal.((x, y, main))
        elseif isa(x, AbstractVector) || isa(y, AbstractVector)
            xy = to_signal(lift_node(x, y) do x, y
                ngrid(x, y)
            end)
            main = (map(first, xy), map(last, xy), main)
        else
            error("surface: combination of types not supported: $(typeof(x)) $(typeof(y)) $(typeof(z))")
        end
    end
    result, main
end

function surface(scene::makie, x, y, z::AbstractMatrix{T}, attributes::Dict) where T <: AbstractFloat
    attributes[:x] = x
    attributes[:y] = y
    attributes[:z] = z
    attributes = surface_defaults(scene, attributes)
    gl_data, main = surface_2glvisualize(attributes)
    viz = visualize(main, Style(:surface), gl_data).children[]
    insert_scene!(scene, :surface, viz, attributes)
end

function surface(scene::makie, x::AbstractMatrix{T1}, y::AbstractMatrix{T2}, f::Function, attributes::Dict) where {T1, T2}
    if size(x) != size(y)
        error("x and y don't have the same size. Found: x: $(size(x)), y: $(size(y))")
    end
    z = f.(x, y)
    surface(scene, x, y, z, attributes::Dict)
end


text:
@testset "2.1 u (F^(c (a+b x)))^n" begin
    (F, a, b, c, d, e, f, g, h, m, n, x, ) = @variables F a b c d e f g h m n x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*u*(F^(c*(a+b*x)))^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*u*F^(c*(a+b*x))=#

    #= Note: The*optimal*antiderivatives*in*this*file*are*for*when*the*control*variable $UseGamma*is*False. =#

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

    @test_int [F^(c*(a + b*x))*(d + e*x)^m, x, 1, (F^(c*(a - (b*d)/e))*(d + e*x)^m*Gamma(1 + m, -((b*c*(d + e*x)*log(F))/e)))/((-((b*c*(d + e*x)*log(F))/e))^m*(b*c*log(F)))]

    @test_int [F^(c*(a + b*x))*(d + e*x)^4, x, 5, (24*e^4*F^(c*(a + b*x)))/(b^5*c^5*log(F)^5) - (24*e^3*F^(c*(a + b*x))*(d + e*x))/(b^4*c^4*log(F)^4) + (12*e^2*F^(c*(a + b*x))*(d + e*x)^2)/(b^3*c^3*log(F)^3) - (4*e*F^(c*(a + b*x))*(d + e*x)^3)/(b^2*c^2*log(F)^2) + (F^(c*(a + b*x))*(d + e*x)^4)/(b*c*log(F))]
    @test_int [F^(c*(a + b*x))*(d + e*x)^3, x, 4, -((6*e^3*F^(c*(a + b*x)))/(b^4*c^4*log(F)^4)) + (6*e^2*F^(c*(a + b*x))*(d + e*x))/(b^3*c^3*log(F)^3) - (3*e*F^(c*(a + b*x))*(d + e*x)^2)/(b^2*c^2*log(F)^2) + (F^(c*(a + b*x))*(d + e*x)^3)/(b*c*log(F))]
    @test_int [F^(c*(a + b*x))*(d + e*x)^2, x, 3, (2*e^2*F^(c*(a + b*x)))/(b^3*c^3*log(F)^3) - (2*e*F^(c*(a + b*x))*(d + e*x))/(b^2*c^2*log(F)^2) + (F^(c*(a + b*x))*(d + e*x)^2)/(b*c*log(F))]
    @test_int [F^(c*(a + b*x))*(d + e*x)^1, x, 2, -((e*F^(c*(a + b*x)))/(b^2*c^2*log(F)^2)) + (F^(c*(a + b*x))*(d + e*x))/(b*c*log(F))]
    @test_int [F^(c*(a + b*x))*(d + e*x)^0, x, 1, F^(c*(a + b*x))/(b*c*log(F))]
    @test_int [F^(c*(a + b*x))/(d + e*x)^1, x, 1, (F^(c*(a - (b*d)/e))*ExpIntegralEi((b*c*(d + e*x)*log(F))/e))/e]
    @test_int [F^(c*(a + b*x))/(d + e*x)^2, x, 2, -(F^(c*(a + b*x))/(e*(d + e*x))) + (b*c*F^(c*(a - (b*d)/e))*ExpIntegralEi((b*c*(d + e*x)*log(F))/e)*log(F))/e^2]
    @test_int [F^(c*(a + b*x))/(d + e*x)^3, x, 3, -(F^(c*(a + b*x))/(2*e*(d + e*x)^2)) - (b*c*F^(c*(a + b*x))*log(F))/(2*e^2*(d + e*x)) + (b^2*c^2*F^(c*(a - (b*d)/e))*ExpIntegralEi((b*c*(d + e*x)*log(F))/e)*log(F)^2)/(2*e^3)]
    @test_int [F^(c*(a + b*x))/(d + e*x)^4, x, 4, -(F^(c*(a + b*x))/(3*e*(d + e*x)^3)) - (b*c*F^(c*(a + b*x))*log(F))/(6*e^2*(d + e*x)^2) - (b^2*c^2*F^(c*(a + b*x))*log(F)^2)/(6*e^3*(d + e*x)) + (b^3*c^3*F^(c*(a - (b*d)/e))*ExpIntegralEi((b*c*(d + e*x)*log(F))/e)*log(F)^3)/(6*e^4)]
    @test_int [F^(c*(a + b*x))/(d + e*x)^5, x, 5, -(F^(c*(a + b*x))/(4*e*(d + e*x)^4)) - (b*c*F^(c*(a + b*x))*log(F))/(12*e^2*(d + e*x)^3) - (b^2*c^2*F^(c*(a + b*x))*log(F)^2)/(24*e^3*(d + e*x)^2) - (b^3*c^3*F^(c*(a + b*x))*log(F)^3)/(24*e^4*(d + e*x)) + (b^4*c^4*F^(c*(a - (b*d)/e))*ExpIntegralEi((b*c*(d + e*x)*log(F))/e)*log(F)^4)/(24*e^5)]

    @test_int [F^(c*(a + b*x))*expand((d + e*x)^4), x, 6, (24*e^4*F^(c*(a + b*x)))/(b^5*c^5*log(F)^5) - (24*e^3*F^(c*(a + b*x))*(d + e*x))/(b^4*c^4*log(F)^4) + (12*e^2*F^(c*(a + b*x))*(d + e*x)^2)/(b^3*c^3*log(F)^3) - (4*e*F^(c*(a + b*x))*(d + e*x)^3)/(b^2*c^2*log(F)^2) + (F^(c*(a + b*x))*(d + e*x)^4)/(b*c*log(F))]
    @test_int [F^(c*(a + b*x))*expand((d + e*x)^3), x, 5, -((6*e^3*F^(c*(a + b*x)))/(b^4*c^4*log(F)^4)) + (6*e^2*F^(c*(a + b*x))*(d + e*x))/(b^3*c^3*log(F)^3) - (3*e*F^(c*(a + b*x))*(d + e*x)^2)/(b^2*c^2*log(F)^2) + (F^(c*(a + b*x))*(d + e*x)^3)/(b*c*log(F))]
    @test_int [F^(c*(a + b*x))*expand((d + e*x)^2), x, 4, (2*e^2*F^(c*(a + b*x)))/(b^3*c^3*log(F)^3) - (2*e*F^(c*(a + b*x))*(d + e*x))/(b^2*c^2*log(F)^2) + (F^(c*(a + b*x))*(d + e*x)^2)/(b*c*log(F))]

    @test_int [F^(c*(a + b*x))/expand((d + e*x)^2), x, 3, -(F^(c*(a + b*x))/(e*(d + e*x))) + (b*c*F^(c*(a - (b*d)/e))*ExpIntegralEi((b*c*(d + e*x)*log(F))/e)*log(F))/e^2]
    @test_int [F^(c*(a + b*x))/expand((d + e*x)^3), x, 4, -(F^(c*(a + b*x))/(2*e*(d + e*x)^2)) - (b*c*F^(c*(a + b*x))*log(F))/(2*e^2*(d + e*x)) + (b^2*c^2*F^(c*(a - (b*d)/e))*ExpIntegralEi((b*c*(d + e*x)*log(F))/e)*log(F)^2)/(2*e^3)]
    @test_int [F^(c*(a + b*x))/expand((d + e*x)^4), x, 5, -(F^(c*(a + b*x))/(3*e*(d + e*x)^3)) - (b*c*F^(c*(a + b*x))*log(F))/(6*e^2*(d + e*x)^2) - (b^2*c^2*F^(c*(a + b*x))*log(F)^2)/(6*e^3*(d + e*x)) + (b^3*c^3*F^(c*(a - (b*d)/e))*ExpIntegralEi((b*c*(d + e*x)*log(F))/e)*log(F)^3)/(6*e^4)]
    @test_int [F^(c*(a + b*x))/expand((d + e*x)^5), x, 6, -(F^(c*(a + b*x))/(4*e*(d + e*x)^4)) - (b*c*F^(c*(a + b*x))*log(F))/(12*e^2*(d + e*x)^3) - (b^2*c^2*F^(c*(a + b*x))*log(F)^2)/(24*e^3*(d + e*x)^2) - (b^3*c^3*F^(c*(a + b*x))*log(F)^3)/(24*e^4*(d + e*x)) + (b^4*c^4*F^(c*(a - (b*d)/e))*ExpIntegralEi((b*c*(d + e*x)*log(F))/e)*log(F)^4)/(24*e^5)]

    @test_int [F^(c*(a + b*x))*expand((d + e*x)^n)^m, x, 2, (F^(c*(a - (b*d)/e))*((d + e*x)^n)^m*Gamma(1 + m*n, -((b*c*(d + e*x)*log(F))/e)))/((-((b*c*(d + e*x)*log(F))/e))^(m*n)*(b*c*log(F)))]

    @test_int [F^(c*(a + b*x))*expand((d + e*x)^4)^m, x, 2, (F^(c*(a - (b*d)/e))*((d + e*x)^4)^m*Gamma(1 + 4*m, -((b*c*(d + e*x)*log(F))/e)))/((-((b*c*(d + e*x)*log(F))/e))^(4*m)*(b*c*log(F)))]
    @test_int [F^(c*(a + b*x))*expand((d + e*x)^3)^m, x, 2, (F^(c*(a - (b*d)/e))*((d + e*x)^3)^m*Gamma(1 + 3*m, -((b*c*(d + e*x)*log(F))/e)))/((-((b*c*(d + e*x)*log(F))/e))^(3*m)*(b*c*log(F)))]
    @test_int [F^(c*(a + b*x))*expand((d + e*x)^2)^m, x, 2, (F^(c*(a - (b*d)/e))*((d + e*x)^2)^m*Gamma(1 + 2*m, -((b*c*(d + e*x)*log(F))/e)))/((-((b*c*(d + e*x)*log(F))/e))^(2*m)*(b*c*log(F)))]
    @test_int [F^(c*(a + b*x))*expand((d + e*x)^1)^m, x, 1, (F^(c*(a - (b*d)/e))*(d + e*x)^m*Gamma(1 + m, -((b*c*(d + e*x)*log(F))/e)))/((-((b*c*(d + e*x)*log(F))/e))^m*(b*c*log(F)))]
    @test_int [F^(c*(a + b*x))/expand((d + e*x)^1)^m, x, 1, (F^(c*(a - (b*d)/e))*Gamma(1 - m, -((b*c*(d + e*x)*log(F))/e))*(-((b*c*(d + e*x)*log(F))/e))^m)/((d + e*x)^m*(b*c*log(F)))]
    @test_int [F^(c*(a + b*x))/expand((d + e*x)^2)^m, x, 2, (F^(c*(a - (b*d)/e))*Gamma(1 - 2*m, -((b*c*(d + e*x)*log(F))/e))*(-((b*c*(d + e*x)*log(F))/e))^(2*m))/(((d + e*x)^2)^m*(b*c*log(F)))]
    @test_int [F^(c*(a + b*x))/expand((d + e*x)^3)^m, x, 2, (F^(c*(a - (b*d)/e))*Gamma(1 - 3*m, -((b*c*(d + e*x)*log(F))/e))*(-((b*c*(d + e*x)*log(F))/e))^(3*m))/(((d + e*x)^3)^m*(b*c*log(F)))]

    @test_int [F^(2 + 5*x), x, 1, F^(2 + 5*x)/(5*log(F))]
    @test_int [F^(a + b*x), x, 1, F^(a + b*x)/(b*log(F))]
    @test_int [10^(2 + 5*x), x, 1, (2^(2 + 5*x)*5^(1 + 5*x))/log(10)]

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

    @test_int [x^(7/2)*F^(a + b*x), x, 6, (105*F^a*sqrt(pi)*erfi(sqrt(b)*sqrt(x)*sqrt(log(F))))/(16*b^(9/2)*log(F)^(9/2)) - (105*F^(a + b*x)*sqrt(x))/(8*b^4*log(F)^4) + (35*F^(a + b*x)*x^(3/2))/(4*b^3*log(F)^3) - (7*F^(a + b*x)*x^(5/2))/(2*b^2*log(F)^2) + (F^(a + b*x)*x^(7/2))/(b*log(F))]
    @test_int [x^(5/2)*F^(a + b*x), x, 5, -((15*F^a*sqrt(pi)*erfi(sqrt(b)*sqrt(x)*sqrt(log(F))))/(8*b^(7/2)*log(F)^(7/2))) + (15*F^(a + b*x)*sqrt(x))/(4*b^3*log(F)^3) - (5*F^(a + b*x)*x^(3/2))/(2*b^2*log(F)^2) + (F^(a + b*x)*x^(5/2))/(b*log(F))]
    @test_int [x^(3/2)*F^(a + b*x), x, 4, (3*F^a*sqrt(pi)*erfi(sqrt(b)*sqrt(x)*sqrt(log(F))))/(4*b^(5/2)*log(F)^(5/2)) - (3*F^(a + b*x)*sqrt(x))/(2*b^2*log(F)^2) + (F^(a + b*x)*x^(3/2))/(b*log(F))]
    @test_int [x^(1/2)*F^(a + b*x), x, 3, -((F^a*sqrt(pi)*erfi(sqrt(b)*sqrt(x)*sqrt(log(F))))/(2*b^(3/2)*log(F)^(3/2))) + (F^(a + b*x)*sqrt(x))/(b*log(F))]
    @test_int [F^(a + b*x)/x^(1/2), x, 2, (F^a*sqrt(pi)*erfi(sqrt(b)*sqrt(x)*sqrt(log(F))))/(sqrt(b)*sqrt(log(F)))]
    @test_int [F^(a + b*x)/x^(3/2), x, 3, -((2*F^(a + b*x))/sqrt(x)) + 2*sqrt(b)*F^a*sqrt(pi)*erfi(sqrt(b)*sqrt(x)*sqrt(log(F)))*sqrt(log(F))]
    @test_int [F^(a + b*x)/x^(5/2), x, 4, -((2*F^(a + b*x))/(3*x^(3/2))) - (4*b*F^(a + b*x)*log(F))/(3*sqrt(x)) + (4/3)*b^(3/2)*F^a*sqrt(pi)*erfi(sqrt(b)*sqrt(x)*sqrt(log(F)))*log(F)^(3/2)]
    @test_int [F^(a + b*x)/x^(7/2), x, 5, -((2*F^(a + b*x))/(5*x^(5/2))) - (4*b*F^(a + b*x)*log(F))/(15*x^(3/2)) - (8*b^2*F^(a + b*x)*log(F)^2)/(15*sqrt(x)) + (8/15)*b^(5/2)*F^a*sqrt(pi)*erfi(sqrt(b)*sqrt(x)*sqrt(log(F)))*log(F)^(5/2)]
    @test_int [F^(a + b*x)/x^(9/2), x, 6, -((2*F^(a + b*x))/(7*x^(7/2))) - (4*b*F^(a + b*x)*log(F))/(35*x^(5/2)) - (8*b^2*F^(a + b*x)*log(F)^2)/(105*x^(3/2)) - (16*b^3*F^(a + b*x)*log(F)^3)/(105*sqrt(x)) + (16/105)*b^(7/2)*F^a*sqrt(pi)*erfi(sqrt(b)*sqrt(x)*sqrt(log(F)))*log(F)^(7/2)]

    @test_int [F^(c*(a + b*x))*(d + e*x)^(7/2), x, 6, (105*e^(7/2)*F^(c*(a - (b*d)/e))*sqrt(pi)*erfi((sqrt(b)*sqrt(c)*sqrt(d + e*x)*sqrt(log(F)))/sqrt(e)))/(16*b^(9/2)*c^(9/2)*log(F)^(9/2)) - (105*e^3*F^(c*(a + b*x))*sqrt(d + e*x))/(8*b^4*c^4*log(F)^4) + (35*e^2*F^(c*(a + b*x))*(d + e*x)^(3/2))/(4*b^3*c^3*log(F)^3) - (7*e*F^(c*(a + b*x))*(d + e*x)^(5/2))/(2*b^2*c^2*log(F)^2) + (F^(c*(a + b*x))*(d + e*x)^(7/2))/(b*c*log(F))]
    @test_int [F^(c*(a + b*x))*(d + e*x)^(5/2), x, 5, -((15*e^(5/2)*F^(c*(a - (b*d)/e))*sqrt(pi)*erfi((sqrt(b)*sqrt(c)*sqrt(d + e*x)*sqrt(log(F)))/sqrt(e)))/(8*b^(7/2)*c^(7/2)*log(F)^(7/2))) + (15*e^2*F^(c*(a + b*x))*sqrt(d + e*x))/(4*b^3*c^3*log(F)^3) - (5*e*F^(c*(a + b*x))*(d + e*x)^(3/2))/(2*b^2*c^2*log(F)^2) + (F^(c*(a + b*x))*(d + e*x)^(5/2))/(b*c*log(F))]
    @test_int [F^(c*(a + b*x))*(d + e*x)^(3/2), x, 4, (3*e^(3/2)*F^(c*(a - (b*d)/e))*sqrt(pi)*erfi((sqrt(b)*sqrt(c)*sqrt(d + e*x)*sqrt(log(F)))/sqrt(e)))/(4*b^(5/2)*c^(5/2)*log(F)^(5/2)) - (3*e*F^(c*(a + b*x))*sqrt(d + e*x))/(2*b^2*c^2*log(F)^2) + (F^(c*(a + b*x))*(d + e*x)^(3/2))/(b*c*log(F))]
    @test_int [F^(c*(a + b*x))*(d + e*x)^(1/2), x, 3, -((sqrt(e)*F^(c*(a - (b*d)/e))*sqrt(pi)*erfi((sqrt(b)*sqrt(c)*sqrt(d + e*x)*sqrt(log(F)))/sqrt(e)))/(2*b^(3/2)*c^(3/2)*log(F)^(3/2))) + (F^(c*(a + b*x))*sqrt(d + e*x))/(b*c*log(F))]
    @test_int [F^(c*(a + b*x))/(d + e*x)^(1/2), x, 2, (F^(c*(a - (b*d)/e))*sqrt(pi)*erfi((sqrt(b)*sqrt(c)*sqrt(d + e*x)*sqrt(log(F)))/sqrt(e)))/(sqrt(b)*sqrt(c)*sqrt(e)*sqrt(log(F)))]
    @test_int [F^(c*(a + b*x))/(d + e*x)^(3/2), x, 3, -((2*F^(c*(a + b*x)))/(e*sqrt(d + e*x))) + (2*sqrt(b)*sqrt(c)*F^(c*(a - (b*d)/e))*sqrt(pi)*erfi((sqrt(b)*sqrt(c)*sqrt(d + e*x)*sqrt(log(F)))/sqrt(e))*sqrt(log(F)))/e^(3/2)]
    @test_int [F^(c*(a + b*x))/(d + e*x)^(5/2), x, 4, -((2*F^(c*(a + b*x)))/(3*e*(d + e*x)^(3/2))) - (4*b*c*F^(c*(a + b*x))*log(F))/(3*e^2*sqrt(d + e*x)) + (4*b^(3/2)*c^(3/2)*F^(c*(a - (b*d)/e))*sqrt(pi)*erfi((sqrt(b)*sqrt(c)*sqrt(d + e*x)*sqrt(log(F)))/sqrt(e))*log(F)^(3/2))/(3*e^(5/2))]
    @test_int [F^(c*(a + b*x))/(d + e*x)^(7/2), x, 5, -((2*F^(c*(a + b*x)))/(5*e*(d + e*x)^(5/2))) - (4*b*c*F^(c*(a + b*x))*log(F))/(15*e^2*(d + e*x)^(3/2)) - (8*b^2*c^2*F^(c*(a + b*x))*log(F)^2)/(15*e^3*sqrt(d + e*x)) + (8*b^(5/2)*c^(5/2)*F^(c*(a - (b*d)/e))*sqrt(pi)*erfi((sqrt(b)*sqrt(c)*sqrt(d + e*x)*sqrt(log(F)))/sqrt(e))*log(F)^(5/2))/(15*e^(7/2))]
    @test_int [F^(c*(a + b*x))/(d + e*x)^(9/2), x, 6, -((2*F^(c*(a + b*x)))/(7*e*(d + e*x)^(7/2))) - (4*b*c*F^(c*(a + b*x))*log(F))/(35*e^2*(d + e*x)^(5/2)) - (8*b^2*c^2*F^(c*(a + b*x))*log(F)^2)/(105*e^3*(d + e*x)^(3/2)) - (16*b^3*c^3*F^(c*(a + b*x))*log(F)^3)/(105*e^4*sqrt(d + e*x)) + (16*b^(7/2)*c^(7/2)*F^(c*(a - (b*d)/e))*sqrt(pi)*erfi((sqrt(b)*sqrt(c)*sqrt(d + e*x)*sqrt(log(F)))/sqrt(e))*log(F)^(7/2))/(105*e^(9/2))]

    @test_int [x^(13/2)/ℯ^(b*x), x, 9, -((135135*sqrt(x))/(ℯ^(b*x)*(64*b^7))) - (45045*x^(3/2))/(ℯ^(b*x)*(32*b^6)) - (9009*x^(5/2))/(ℯ^(b*x)*(16*b^5)) - (1287*x^(7/2))/(ℯ^(b*x)*(8*b^4)) - (143*x^(9/2))/(ℯ^(b*x)*(4*b^3)) - (13*x^(11/2))/(ℯ^(b*x)*(2*b^2)) - x^(13/2)/(ℯ^(b*x)*b) + (135135*sqrt(pi)*erf(sqrt(b)*sqrt(x)))/(128*b^(15/2))]

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

    @test_int [(d + e*x)^(4/3)*F^(c*(a + b*x)), x, 1, -((e*F^(c*(a - (b*d)/e))*(d + e*x)^(1/3)*Gamma(7/3, -((b*c*(d + e*x)*log(F))/e)))/(b^2*c^2*log(F)^2*(-((b*c*(d + e*x)*log(F))/e))^(1/3)))]

    @test_int [(d + e*x)^(4/3)*(F^(c*(a + b*x)))^n, x, 2, -((e*F^(c*(a - (b*d)/e)*n - c*n*(a + b*x))*(F^(c*(a + b*x)))^n*(d + e*x)^(1/3)*Gamma(7/3, -((b*c*n*(d + e*x)*log(F))/e)))/(b^2*c^2*n^2*log(F)^2*(-((b*c*n*(d + e*x)*log(F))/e))^(1/3)))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*Px*F^(c*(a+b*x))=#

    @test_int [F^(c*(a + b*x))*(d + e*x), x, 2, -((e*F^(c*(a + b*x)))/(b^2*c^2*log(F)^2)) + (F^(c*(a + b*x))*(d + e*x))/(b*c*log(F))]
    @test_int [F^(c*(a + b*x))*(d + e*x + f*x^2), x, 8, (2*f*F^(c*(a + b*x)))/(b^3*c^3*log(F)^3) - (e*F^(c*(a + b*x)))/(b^2*c^2*log(F)^2) - (2*f*F^(c*(a + b*x))*x)/(b^2*c^2*log(F)^2) + (d*F^(c*(a + b*x)))/(b*c*log(F)) + (e*F^(c*(a + b*x))*x)/(b*c*log(F)) + (f*F^(c*(a + b*x))*x^2)/(b*c*log(F))]
    @test_int [F^(c*(a + b*x))*(d + e*x + f*x^2 + g*x^3), x, 12, -((6*F^(c*(a + b*x))*g)/(b^4*c^4*log(F)^4)) + (2*f*F^(c*(a + b*x)))/(b^3*c^3*log(F)^3) + (6*F^(c*(a + b*x))*g*x)/(b^3*c^3*log(F)^3) - (e*F^(c*(a + b*x)))/(b^2*c^2*log(F)^2) - (2*f*F^(c*(a + b*x))*x)/(b^2*c^2*log(F)^2) - (3*F^(c*(a + b*x))*g*x^2)/(b^2*c^2*log(F)^2) + (d*F^(c*(a + b*x)))/(b*c*log(F)) + (e*F^(c*(a + b*x))*x)/(b*c*log(F)) + (f*F^(c*(a + b*x))*x^2)/(b*c*log(F)) + (F^(c*(a + b*x))*g*x^3)/(b*c*log(F))]
    @test_int [F^(c*(a + b*x))*(d + e*x + f*x^2 + g*x^3 + h*x^4), x, 17, (24*F^(c*(a + b*x))*h)/(b^5*c^5*log(F)^5) - (6*F^(c*(a + b*x))*g)/(b^4*c^4*log(F)^4) - (24*F^(c*(a + b*x))*h*x)/(b^4*c^4*log(F)^4) + (2*f*F^(c*(a + b*x)))/(b^3*c^3*log(F)^3) + (6*F^(c*(a + b*x))*g*x)/(b^3*c^3*log(F)^3) + (12*F^(c*(a + b*x))*h*x^2)/(b^3*c^3*log(F)^3) - (e*F^(c*(a + b*x)))/(b^2*c^2*log(F)^2) - (2*f*F^(c*(a + b*x))*x)/(b^2*c^2*log(F)^2) - (3*F^(c*(a + b*x))*g*x^2)/(b^2*c^2*log(F)^2) - (4*F^(c*(a + b*x))*h*x^3)/(b^2*c^2*log(F)^2) + (d*F^(c*(a + b*x)))/(b*c*log(F)) + (e*F^(c*(a + b*x))*x)/(b*c*log(F)) + (f*F^(c*(a + b*x))*x^2)/(b*c*log(F)) + (F^(c*(a + b*x))*g*x^3)/(b*c*log(F)) + (F^(c*(a + b*x))*h*x^4)/(b*c*log(F))]

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

    @test_int [x^m*(a + b*x)^3/ℯ^(a + b*x), x, 6, -((a^3*x^m*Gamma(1 + m, b*x))/(ℯ^a*(b*x)^m*b)) - (3*a^2*x^m*Gamma(2 + m, b*x))/(ℯ^a*(b*x)^m*b) - (3*a*x^m*Gamma(3 + m, b*x))/(ℯ^a*(b*x)^m*b) - (x^m*Gamma(4 + m, b*x))/(ℯ^a*(b*x)^m*b)]

    @test_int [x^3*(a + b*x)^3/ℯ^(a + b*x), x, 24, -((720*ℯ^(-a - b*x))/b^4) - (360*a*ℯ^(-a - b*x))/b^4 - (72*a^2*ℯ^(-a - b*x))/b^4 - (6*a^3*ℯ^(-a - b*x))/b^4 - (720*ℯ^(-a - b*x)*x)/b^3 - (360*a*ℯ^(-a - b*x)*x)/b^3 - (72*a^2*ℯ^(-a - b*x)*x)/b^3 - (6*a^3*ℯ^(-a - b*x)*x)/b^3 - (360*ℯ^(-a - b*x)*x^2)/b^2 - (180*a*ℯ^(-a - b*x)*x^2)/b^2 - (36*a^2*ℯ^(-a - b*x)*x^2)/b^2 - (3*a^3*ℯ^(-a - b*x)*x^2)/b^2 - (120*ℯ^(-a - b*x)*x^3)/b - (60*a*ℯ^(-a - b*x)*x^3)/b - (12*a^2*ℯ^(-a - b*x)*x^3)/b - (a^3*ℯ^(-a - b*x)*x^3)/b - 30*ℯ^(-a - b*x)*x^4 - 15*a*ℯ^(-a - b*x)*x^4 - 3*a^2*ℯ^(-a - b*x)*x^4 - 6*b*ℯ^(-a - b*x)*x^5 - 3*a*b*ℯ^(-a - b*x)*x^5 - b^2*ℯ^(-a - b*x)*x^6]
    @test_int [x^2*(a + b*x)^3/ℯ^(a + b*x), x, 20, -((120*ℯ^(-a - b*x))/b^3) - (72*a*ℯ^(-a - b*x))/b^3 - (18*a^2*ℯ^(-a - b*x))/b^3 - (2*a^3*ℯ^(-a - b*x))/b^3 - (120*ℯ^(-a - b*x)*x)/b^2 - (72*a*ℯ^(-a - b*x)*x)/b^2 - (18*a^2*ℯ^(-a - b*x)*x)/b^2 - (2*a^3*ℯ^(-a - b*x)*x)/b^2 - (60*ℯ^(-a - b*x)*x^2)/b - (36*a*ℯ^(-a - b*x)*x^2)/b - (9*a^2*ℯ^(-a - b*x)*x^2)/b - (a^3*ℯ^(-a - b*x)*x^2)/b - 20*ℯ^(-a - b*x)*x^3 - 12*a*ℯ^(-a - b*x)*x^3 - 3*a^2*ℯ^(-a - b*x)*x^3 - 5*b*ℯ^(-a - b*x)*x^4 - 3*a*b*ℯ^(-a - b*x)*x^4 - b^2*ℯ^(-a - b*x)*x^5]
    @test_int [x^1*(a + b*x)^3/ℯ^(a + b*x), x, 11, -((24*ℯ^(-a - b*x))/b^2) + (6*a*ℯ^(-a - b*x))/b^2 - (24*ℯ^(-a - b*x)*(a + b*x))/b^2 + (6*a*ℯ^(-a - b*x)*(a + b*x))/b^2 - (12*ℯ^(-a - b*x)*(a + b*x)^2)/b^2 + (3*a*ℯ^(-a - b*x)*(a + b*x)^2)/b^2 - (4*ℯ^(-a - b*x)*(a + b*x)^3)/b^2 + (a*ℯ^(-a - b*x)*(a + b*x)^3)/b^2 - (ℯ^(-a - b*x)*(a + b*x)^4)/b^2]
    @test_int [x^0*(a + b*x)^3/ℯ^(a + b*x), x, 4, -((6*ℯ^(-a - b*x))/b) - (6*ℯ^(-a - b*x)*(a + b*x))/b - (3*ℯ^(-a - b*x)*(a + b*x)^2)/b - (ℯ^(-a - b*x)*(a + b*x)^3)/b]
    @test_int [(a + b*x)^3/(x^1*ℯ^(a + b*x)), x, 9, -2*ℯ^(-a - b*x) - 3*a*ℯ^(-a - b*x) - 3*a^2*ℯ^(-a - b*x) - 2*b*ℯ^(-a - b*x)*x - 3*a*b*ℯ^(-a - b*x)*x - b^2*ℯ^(-a - b*x)*x^2 + (a^3*ExpIntegralEi((-b)*x))/ℯ^a]
    @test_int [(a + b*x)^3/(x^2*ℯ^(a + b*x)), x, 8, (-b)*ℯ^(-a - b*x) - 3*a*b*ℯ^(-a - b*x) - (a^3*ℯ^(-a - b*x))/x - b^2*ℯ^(-a - b*x)*x + (3*a^2*b*ExpIntegralEi((-b)*x))/ℯ^a - (a^3*b*ExpIntegralEi((-b)*x))/ℯ^a]
    @test_int [(a + b*x)^3/(x^3*ℯ^(a + b*x)), x, 9, (-b^2)*ℯ^(-a - b*x) - (a^3*ℯ^(-a - b*x))/(2*x^2) - (3*a^2*b*ℯ^(-a - b*x))/x + (a^3*b*ℯ^(-a - b*x))/(2*x) + (3*a*b^2*ExpIntegralEi((-b)*x))/ℯ^a - (3*a^2*b^2*ExpIntegralEi((-b)*x))/ℯ^a + ((1/2)*a^3*b^2*ExpIntegralEi((-b)*x))/ℯ^a]
    @test_int [(a + b*x)^3/(x^4*ℯ^(a + b*x)), x, 12, -((a^3*ℯ^(-a - b*x))/(3*x^3)) - (3*a^2*b*ℯ^(-a - b*x))/(2*x^2) + (a^3*b*ℯ^(-a - b*x))/(6*x^2) - (3*a*b^2*ℯ^(-a - b*x))/x + (3*a^2*b^2*ℯ^(-a - b*x))/(2*x) - (a^3*b^2*ℯ^(-a - b*x))/(6*x) + (b^3*ExpIntegralEi((-b)*x))/ℯ^a - (3*a*b^3*ExpIntegralEi((-b)*x))/ℯ^a + ((3/2)*a^2*b^3*ExpIntegralEi((-b)*x))/ℯ^a - ((1/6)*a^3*b^3*ExpIntegralEi((-b)*x))/ℯ^a]

    @test_int [x^m*(e + f*x)^2*F^(a + b*(c + d*x)), x, 5, (f^2*F^(a + b*c)*x^m*Gamma(3 + m, (-b)*d*x*log(F)))/(((-b)*d*x*log(F))^m*(b^3*d^3*log(F)^3)) - (2*e*f*F^(a + b*c)*x^m*Gamma(2 + m, (-b)*d*x*log(F)))/(((-b)*d*x*log(F))^m*(b^2*d^2*log(F)^2)) + (e^2*F^(a + b*c)*x^m*Gamma(1 + m, (-b)*d*x*log(F)))/(((-b)*d*x*log(F))^m*(b*d*log(F)))]

    @test_int [x^3*(e + f*x)^2*F^(a + b*(c + d*x)), x, 17, -((120*f^2*F^(a + b*c + b*d*x))/(b^6*d^6*log(F)^6)) + (48*e*f*F^(a + b*c + b*d*x))/(b^5*d^5*log(F)^5) + (120*f^2*F^(a + b*c + b*d*x)*x)/(b^5*d^5*log(F)^5) - (6*e^2*F^(a + b*c + b*d*x))/(b^4*d^4*log(F)^4) - (48*e*f*F^(a + b*c + b*d*x)*x)/(b^4*d^4*log(F)^4) - (60*f^2*F^(a + b*c + b*d*x)*x^2)/(b^4*d^4*log(F)^4) + (6*e^2*F^(a + b*c + b*d*x)*x)/(b^3*d^3*log(F)^3) + (24*e*f*F^(a + b*c + b*d*x)*x^2)/(b^3*d^3*log(F)^3) + (20*f^2*F^(a + b*c + b*d*x)*x^3)/(b^3*d^3*log(F)^3) - (3*e^2*F^(a + b*c + b*d*x)*x^2)/(b^2*d^2*log(F)^2) - (8*e*f*F^(a + b*c + b*d*x)*x^3)/(b^2*d^2*log(F)^2) - (5*f^2*F^(a + b*c + b*d*x)*x^4)/(b^2*d^2*log(F)^2) + (e^2*F^(a + b*c + b*d*x)*x^3)/(b*d*log(F)) + (2*e*f*F^(a + b*c + b*d*x)*x^4)/(b*d*log(F)) + (f^2*F^(a + b*c + b*d*x)*x^5)/(b*d*log(F))]
    @test_int [x^2*(e + f*x)^2*F^(a + b*(c + d*x)), x, 14, (24*f^2*F^(a + b*c + b*d*x))/(b^5*d^5*log(F)^5) - (12*e*f*F^(a + b*c + b*d*x))/(b^4*d^4*log(F)^4) - (24*f^2*F^(a + b*c + b*d*x)*x)/(b^4*d^4*log(F)^4) + (2*e^2*F^(a + b*c + b*d*x))/(b^3*d^3*log(F)^3) + (12*e*f*F^(a + b*c + b*d*x)*x)/(b^3*d^3*log(F)^3) + (12*f^2*F^(a + b*c + b*d*x)*x^2)/(b^3*d^3*log(F)^3) - (2*e^2*F^(a + b*c + b*d*x)*x)/(b^2*d^2*log(F)^2) - (6*e*f*F^(a + b*c + b*d*x)*x^2)/(b^2*d^2*log(F)^2) - (4*f^2*F^(a + b*c + b*d*x)*x^3)/(b^2*d^2*log(F)^2) + (e^2*F^(a + b*c + b*d*x)*x^2)/(b*d*log(F)) + (2*e*f*F^(a + b*c + b*d*x)*x^3)/(b*d*log(F)) + (f^2*F^(a + b*c + b*d*x)*x^4)/(b*d*log(F))]
    @test_int [x^1*(e + f*x)^2*F^(a + b*(c + d*x)), x, 11, -((6*f^2*F^(a + b*c + b*d*x))/(b^4*d^4*log(F)^4)) + (4*e*f*F^(a + b*c + b*d*x))/(b^3*d^3*log(F)^3) + (6*f^2*F^(a + b*c + b*d*x)*x)/(b^3*d^3*log(F)^3) - (e^2*F^(a + b*c + b*d*x))/(b^2*d^2*log(F)^2) - (4*e*f*F^(a + b*c + b*d*x)*x)/(b^2*d^2*log(F)^2) - (3*f^2*F^(a + b*c + b*d*x)*x^2)/(b^2*d^2*log(F)^2) + (e^2*F^(a + b*c + b*d*x)*x)/(b*d*log(F)) + (2*e*f*F^(a + b*c + b*d*x)*x^2)/(b*d*log(F)) + (f^2*F^(a + b*c + b*d*x)*x^3)/(b*d*log(F))]
    @test_int [x^0*(e + f*x)^2*F^(a + b*(c + d*x)), x, 4, (2*f^2*F^(a + b*c + b*d*x))/(b^3*d^3*log(F)^3) - (2*f*F^(a + b*c + b*d*x)*(e + f*x))/(b^2*d^2*log(F)^2) + (F^(a + b*c + b*d*x)*(e + f*x)^2)/(b*d*log(F))]
    @test_int [(e + f*x)^2*F^(a + b*(c + d*x))/x^1, x, 6, e^2*F^(a + b*c)*ExpIntegralEi(b*d*x*log(F)) - (f^2*F^(a + b*c + b*d*x))/(b^2*d^2*log(F)^2) + (2*e*f*F^(a + b*c + b*d*x))/(b*d*log(F)) + (f^2*F^(a + b*c + b*d*x)*x)/(b*d*log(F))]
    @test_int [(e + f*x)^2*F^(a + b*(c + d*x))/x^2, x, 6, -((e^2*F^(a + b*c + b*d*x))/x) + 2*e*f*F^(a + b*c)*ExpIntegralEi(b*d*x*log(F)) + (f^2*F^(a + b*c + b*d*x))/(b*d*log(F)) + b*d*e^2*F^(a + b*c)*ExpIntegralEi(b*d*x*log(F))*log(F)]
    @test_int [(e + f*x)^2*F^(a + b*(c + d*x))/x^3, x, 8, -((e^2*F^(a + b*c + b*d*x))/(2*x^2)) - (2*e*f*F^(a + b*c + b*d*x))/x + f^2*F^(a + b*c)*ExpIntegralEi(b*d*x*log(F)) - (b*d*e^2*F^(a + b*c + b*d*x)*log(F))/(2*x) + 2*b*d*e*f*F^(a + b*c)*ExpIntegralEi(b*d*x*log(F))*log(F) + (1/2)*b^2*d^2*e^2*F^(a + b*c)*ExpIntegralEi(b*d*x*log(F))*log(F)^2]
    @test_int [(e + f*x)^2*F^(a + b*(c + d*x))/x^4, x, 11, -((e^2*F^(a + b*c + b*d*x))/(3*x^3)) - (e*f*F^(a + b*c + b*d*x))/x^2 - (f^2*F^(a + b*c + b*d*x))/x - (b*d*e^2*F^(a + b*c + b*d*x)*log(F))/(6*x^2) - (b*d*e*f*F^(a + b*c + b*d*x)*log(F))/x + b*d*f^2*F^(a + b*c)*ExpIntegralEi(b*d*x*log(F))*log(F) - (b^2*d^2*e^2*F^(a + b*c + b*d*x)*log(F)^2)/(6*x) + b^2*d^2*e*f*F^(a + b*c)*ExpIntegralEi(b*d*x*log(F))*log(F)^2 + (1/6)*b^3*d^3*e^2*F^(a + b*c)*ExpIntegralEi(b*d*x*log(F))*log(F)^3]
    @test_int [(e + f*x)^2*F^(a + b*(c + d*x))/x^5, x, 14, -((e^2*F^(a + b*c + b*d*x))/(4*x^4)) - (2*e*f*F^(a + b*c + b*d*x))/(3*x^3) - (f^2*F^(a + b*c + b*d*x))/(2*x^2) - (b*d*e^2*F^(a + b*c + b*d*x)*log(F))/(12*x^3) - (b*d*e*f*F^(a + b*c + b*d*x)*log(F))/(3*x^2) - (b*d*f^2*F^(a + b*c + b*d*x)*log(F))/(2*x) - (b^2*d^2*e^2*F^(a + b*c + b*d*x)*log(F)^2)/(24*x^2) - (b^2*d^2*e*f*F^(a + b*c + b*d*x)*log(F)^2)/(3*x) + (1/2)*b^2*d^2*f^2*F^(a + b*c)*ExpIntegralEi(b*d*x*log(F))*log(F)^2 - (b^3*d^3*e^2*F^(a + b*c + b*d*x)*log(F)^3)/(24*x) + (1/3)*b^3*d^3*e*f*F^(a + b*c)*ExpIntegralEi(b*d*x*log(F))*log(F)^3 + (1/24)*b^4*d^4*e^2*F^(a + b*c)*ExpIntegralEi(b*d*x*log(F))*log(F)^4]

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

    @test_int [(c + d*x)^3*(a + b*x)^4/ℯ^(a + b*x), x, 28, -((5040*d^3*ℯ^(-a - b*x))/b^4) - (2160*d^2*(b*c - a*d)*ℯ^(-a - b*x))/b^4 - (360*d*(b*c - a*d)^2*ℯ^(-a - b*x))/b^4 - (24*(b*c - a*d)^3*ℯ^(-a - b*x))/b^4 - (5040*d^3*ℯ^(-a - b*x)*(a + b*x))/b^4 - (2160*d^2*(b*c - a*d)*ℯ^(-a - b*x)*(a + b*x))/b^4 - (360*d*(b*c - a*d)^2*ℯ^(-a - b*x)*(a + b*x))/b^4 - (24*(b*c - a*d)^3*ℯ^(-a - b*x)*(a + b*x))/b^4 - (2520*d^3*ℯ^(-a - b*x)*(a + b*x)^2)/b^4 - (1080*d^2*(b*c - a*d)*ℯ^(-a - b*x)*(a + b*x)^2)/b^4 - (180*d*(b*c - a*d)^2*ℯ^(-a - b*x)*(a + b*x)^2)/b^4 - (12*(b*c - a*d)^3*ℯ^(-a - b*x)*(a + b*x)^2)/b^4 - (840*d^3*ℯ^(-a - b*x)*(a + b*x)^3)/b^4 - (360*d^2*(b*c - a*d)*ℯ^(-a - b*x)*(a + b*x)^3)/b^4 - (60*d*(b*c - a*d)^2*ℯ^(-a - b*x)*(a + b*x)^3)/b^4 - (4*(b*c - a*d)^3*ℯ^(-a - b*x)*(a + b*x)^3)/b^4 - (210*d^3*ℯ^(-a - b*x)*(a + b*x)^4)/b^4 - (90*d^2*(b*c - a*d)*ℯ^(-a - b*x)*(a + b*x)^4)/b^4 - (15*d*(b*c - a*d)^2*ℯ^(-a - b*x)*(a + b*x)^4)/b^4 - ((b*c - a*d)^3*ℯ^(-a - b*x)*(a + b*x)^4)/b^4 - (42*d^3*ℯ^(-a - b*x)*(a + b*x)^5)/b^4 - (18*d^2*(b*c - a*d)*ℯ^(-a - b*x)*(a + b*x)^5)/b^4 - (3*d*(b*c - a*d)^2*ℯ^(-a - b*x)*(a + b*x)^5)/b^4 - (7*d^3*ℯ^(-a - b*x)*(a + b*x)^6)/b^4 - (3*d^2*(b*c - a*d)*ℯ^(-a - b*x)*(a + b*x)^6)/b^4 - (d^3*ℯ^(-a - b*x)*(a + b*x)^7)/b^4]
    @test_int [(c + d*x)^2*(a + b*x)^4/ℯ^(a + b*x), x, 20, -((720*d^2*ℯ^(-a - b*x))/b^3) - (240*d*(b*c - a*d)*ℯ^(-a - b*x))/b^3 - (24*(b*c - a*d)^2*ℯ^(-a - b*x))/b^3 - (720*d^2*ℯ^(-a - b*x)*(a + b*x))/b^3 - (240*d*(b*c - a*d)*ℯ^(-a - b*x)*(a + b*x))/b^3 - (24*(b*c - a*d)^2*ℯ^(-a - b*x)*(a + b*x))/b^3 - (360*d^2*ℯ^(-a - b*x)*(a + b*x)^2)/b^3 - (120*d*(b*c - a*d)*ℯ^(-a - b*x)*(a + b*x)^2)/b^3 - (12*(b*c - a*d)^2*ℯ^(-a - b*x)*(a + b*x)^2)/b^3 - (120*d^2*ℯ^(-a - b*x)*(a + b*x)^3)/b^3 - (40*d*(b*c - a*d)*ℯ^(-a - b*x)*(a + b*x)^3)/b^3 - (4*(b*c - a*d)^2*ℯ^(-a - b*x)*(a + b*x)^3)/b^3 - (30*d^2*ℯ^(-a - b*x)*(a + b*x)^4)/b^3 - (10*d*(b*c - a*d)*ℯ^(-a - b*x)*(a + b*x)^4)/b^3 - ((b*c - a*d)^2*ℯ^(-a - b*x)*(a + b*x)^4)/b^3 - (6*d^2*ℯ^(-a - b*x)*(a + b*x)^5)/b^3 - (2*d*(b*c - a*d)*ℯ^(-a - b*x)*(a + b*x)^5)/b^3 - (d^2*ℯ^(-a - b*x)*(a + b*x)^6)/b^3]
    @test_int [(c + d*x)^1*(a + b*x)^4/ℯ^(a + b*x), x, 13, -((120*d*ℯ^(-a - b*x))/b^2) - (24*(b*c - a*d)*ℯ^(-a - b*x))/b^2 - (120*d*ℯ^(-a - b*x)*(a + b*x))/b^2 - (24*(b*c - a*d)*ℯ^(-a - b*x)*(a + b*x))/b^2 - (60*d*ℯ^(-a - b*x)*(a + b*x)^2)/b^2 - (12*(b*c - a*d)*ℯ^(-a - b*x)*(a + b*x)^2)/b^2 - (20*d*ℯ^(-a - b*x)*(a + b*x)^3)/b^2 - (4*(b*c - a*d)*ℯ^(-a - b*x)*(a + b*x)^3)/b^2 - (5*d*ℯ^(-a - b*x)*(a + b*x)^4)/b^2 - ((b*c - a*d)*ℯ^(-a - b*x)*(a + b*x)^4)/b^2 - (d*ℯ^(-a - b*x)*(a + b*x)^5)/b^2]
    @test_int [(c + d*x)^0*(a + b*x)^4/ℯ^(a + b*x), x, 5, -((24*ℯ^(-a - b*x))/b) - (24*ℯ^(-a - b*x)*(a + b*x))/b - (12*ℯ^(-a - b*x)*(a + b*x)^2)/b - (4*ℯ^(-a - b*x)*(a + b*x)^3)/b - (ℯ^(-a - b*x)*(a + b*x)^4)/b]
    @test_int [(a + b*x)^4/((c + d*x)^1*ℯ^(a + b*x)), x, 13, -((6*ℯ^(-a - b*x))/d) + (2*(b*c - a*d)*ℯ^(-a - b*x))/d^2 - ((b*c - a*d)^2*ℯ^(-a - b*x))/d^3 + ((b*c - a*d)^3*ℯ^(-a - b*x))/d^4 - (6*ℯ^(-a - b*x)*(a + b*x))/d + (2*(b*c - a*d)*ℯ^(-a - b*x)*(a + b*x))/d^2 - ((b*c - a*d)^2*ℯ^(-a - b*x)*(a + b*x))/d^3 - (3*ℯ^(-a - b*x)*(a + b*x)^2)/d + ((b*c - a*d)*ℯ^(-a - b*x)*(a + b*x)^2)/d^2 - (ℯ^(-a - b*x)*(a + b*x)^3)/d + ((b*c - a*d)^4*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/d^5]
    @test_int [(a + b*x)^4/((c + d*x)^2*ℯ^(a + b*x)), x, 11, -((2*b*ℯ^(-a - b*x))/d^2) + (4*b*(b*c - a*d)*ℯ^(-a - b*x))/d^3 - (6*b*(b*c - a*d)^2*ℯ^(-a - b*x))/d^4 - ((b*c - a*d)^4*ℯ^(-a - b*x))/(d^5*(c + d*x)) - (2*b^2*ℯ^(-a - b*x)*(c + d*x))/d^3 + (4*b^2*(b*c - a*d)*ℯ^(-a - b*x)*(c + d*x))/d^4 - (b^3*ℯ^(-a - b*x)*(c + d*x)^2)/d^4 - (4*b*(b*c - a*d)^3*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/d^5 - (b*(b*c - a*d)^4*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/d^6]
    @test_int [(a + b*x)^4/((c + d*x)^3*ℯ^(a + b*x)), x, 11, -((b^2*ℯ^(-a - b*x))/d^3) + (b^2*(3*b*c - 4*a*d)*ℯ^(-a - b*x))/d^4 - (b^3*ℯ^(-a - b*x)*x)/d^3 - ((b*c - a*d)^4*ℯ^(-a - b*x))/(2*d^5*(c + d*x)^2) + (4*b*(b*c - a*d)^3*ℯ^(-a - b*x))/(d^5*(c + d*x)) + (b*(b*c - a*d)^4*ℯ^(-a - b*x))/(2*d^6*(c + d*x)) + (6*b^2*(b*c - a*d)^2*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/d^5 + (4*b^2*(b*c - a*d)^3*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/d^6 + (b^2*(b*c - a*d)^4*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/(2*d^7)]
    @test_int [(a + b*x)^4/((c + d*x)^4*ℯ^(a + b*x)), x, 13, -((b^3*ℯ^(-a - b*x))/d^4) - ((b*c - a*d)^4*ℯ^(-a - b*x))/(3*d^5*(c + d*x)^3) + (2*b*(b*c - a*d)^3*ℯ^(-a - b*x))/(d^5*(c + d*x)^2) + (b*(b*c - a*d)^4*ℯ^(-a - b*x))/(6*d^6*(c + d*x)^2) - (6*b^2*(b*c - a*d)^2*ℯ^(-a - b*x))/(d^5*(c + d*x)) - (2*b^2*(b*c - a*d)^3*ℯ^(-a - b*x))/(d^6*(c + d*x)) - (b^2*(b*c - a*d)^4*ℯ^(-a - b*x))/(6*d^7*(c + d*x)) - (4*b^3*(b*c - a*d)*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/d^5 - (6*b^3*(b*c - a*d)^2*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/d^6 - (2*b^3*(b*c - a*d)^3*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/d^7 - (b^3*(b*c - a*d)^4*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/(6*d^8)]
    @test_int [(a + b*x)^4/((c + d*x)^5*ℯ^(a + b*x)), x, 17, -(((b*c - a*d)^4*ℯ^(-a - b*x))/(4*d^5*(c + d*x)^4)) + (4*b*(b*c - a*d)^3*ℯ^(-a - b*x))/(3*d^5*(c + d*x)^3) + (b*(b*c - a*d)^4*ℯ^(-a - b*x))/(12*d^6*(c + d*x)^3) - (3*b^2*(b*c - a*d)^2*ℯ^(-a - b*x))/(d^5*(c + d*x)^2) - (2*b^2*(b*c - a*d)^3*ℯ^(-a - b*x))/(3*d^6*(c + d*x)^2) - (b^2*(b*c - a*d)^4*ℯ^(-a - b*x))/(24*d^7*(c + d*x)^2) + (4*b^3*(b*c - a*d)*ℯ^(-a - b*x))/(d^5*(c + d*x)) + (3*b^3*(b*c - a*d)^2*ℯ^(-a - b*x))/(d^6*(c + d*x)) + (2*b^3*(b*c - a*d)^3*ℯ^(-a - b*x))/(3*d^7*(c + d*x)) + (b^3*(b*c - a*d)^4*ℯ^(-a - b*x))/(24*d^8*(c + d*x)) + (b^4*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/d^5 + (4*b^4*(b*c - a*d)*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/d^6 + (3*b^4*(b*c - a*d)^2*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/d^7 + (2*b^4*(b*c - a*d)^3*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/(3*d^8) + (b^4*(b*c - a*d)^4*ℯ^(-a + (b*c)/d)*ExpIntegralEi(-((b*(c + d*x))/d)))/(24*d^9)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*F(c*(a+b*x))*log(d*x)*(e+(f+g*x)*log(d*x))=#

    @test_int [x^m*F^(c*(a + b*x))*log(d*x)^n*(e + e*n + e*(1 + m + b*c*x*log(F))*log(d*x)), x, 1, e*F^(c*(a + b*x))*x^(1 + m)*log(d*x)^(1 + n)]

    @test_int [x^2*F^(c*(a + b*x))*log(d*x)^n*(e + e*n + e*(1 + 2 + b*c*x*log(F))*log(d*x)), x, 1, e*F^(c*(a + b*x))*x^3*log(d*x)^(1 + n)]
    @test_int [x^1*F^(c*(a + b*x))*log(d*x)^n*(e + e*n + e*(1 + 1 + b*c*x*log(F))*log(d*x)), x, 1, e*F^(c*(a + b*x))*x^2*log(d*x)^(1 + n)]
    @test_int [x^0*F^(c*(a + b*x))*log(d*x)^n*(e + e*n + e*(1 + 0 + b*c*x*log(F))*log(d*x)), x, 1, e*F^(c*(a + b*x))*x*log(d*x)^(1 + n)]
    @test_int [1/x^1*F^(c*(a + b*x))*log(d*x)^n*(e + e*n + e*(1 - 1 + b*c*x*log(F))*log(d*x)), x, 1, e*F^(c*(a + b*x))*log(d*x)^(1 + n)]
    @test_int [1/x^2*F^(c*(a + b*x))*log(d*x)^n*(e + e*n + e*(1 - 2 + b*c*x*log(F))*log(d*x)), x, 1, (e*F^(c*(a + b*x))*log(d*x)^(1 + n))/x]
    @test_int [1/x^3*F^(c*(a + b*x))*log(d*x)^n*(e + e*n + e*(1 - 3 + b*c*x*log(F))*log(d*x)), x, 1, (e*F^(c*(a + b*x))*log(d*x)^(1 + n))/x^2]

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

    @test_int [sqrt(ℯ^(a + b*x))*x^4, x, 5, (768*sqrt(ℯ^(a + b*x)))/b^5 - (384*sqrt(ℯ^(a + b*x))*x)/b^4 + (96*sqrt(ℯ^(a + b*x))*x^2)/b^3 - (16*sqrt(ℯ^(a + b*x))*x^3)/b^2 + (2*sqrt(ℯ^(a + b*x))*x^4)/b]
    @test_int [sqrt(ℯ^(a + b*x))*x^3, x, 4, -((96*sqrt(ℯ^(a + b*x)))/b^4) + (48*sqrt(ℯ^(a + b*x))*x)/b^3 - (12*sqrt(ℯ^(a + b*x))*x^2)/b^2 + (2*sqrt(ℯ^(a + b*x))*x^3)/b]
    @test_int [sqrt(ℯ^(a + b*x))*x^2, x, 3, (16*sqrt(ℯ^(a + b*x)))/b^3 - (8*sqrt(ℯ^(a + b*x))*x)/b^2 + (2*sqrt(ℯ^(a + b*x))*x^2)/b]
    @test_int [sqrt(ℯ^(a + b*x))*x^1, x, 2, -((4*sqrt(ℯ^(a + b*x)))/b^2) + (2*sqrt(ℯ^(a + b*x))*x)/b]
    @test_int [sqrt(ℯ^(a + b*x))*x^0, x, 1, (2*sqrt(ℯ^(a + b*x)))/b]
    @test_int [sqrt(ℯ^(a + b*x))/x^1, x, 2, (sqrt(ℯ^(a + b*x))*ExpIntegralEi((b*x)/2))/ℯ^((b*x)/2)]
    @test_int [sqrt(ℯ^(a + b*x))/x^2, x, 3, -(sqrt(ℯ^(a + b*x))/x) + ((1/2)*b*sqrt(ℯ^(a + b*x))*ExpIntegralEi((b*x)/2))/ℯ^((b*x)/2)]
    @test_int [sqrt(ℯ^(a + b*x))/x^3, x, 4, -(sqrt(ℯ^(a + b*x))/(2*x^2)) - (b*sqrt(ℯ^(a + b*x)))/(4*x) + ((1/8)*b^2*sqrt(ℯ^(a + b*x))*ExpIntegralEi((b*x)/2))/ℯ^((b*x)/2)]
    @test_int [sqrt(ℯ^(a + b*x))/x^4, x, 5, -(sqrt(ℯ^(a + b*x))/(3*x^3)) - (b*sqrt(ℯ^(a + b*x)))/(12*x^2) - (b^2*sqrt(ℯ^(a + b*x)))/(24*x) + ((1/48)*b^3*sqrt(ℯ^(a + b*x))*ExpIntegralEi((b*x)/2))/ℯ^((b*x)/2)]
    end

max_stars_repo_path: test/2 Exponentials/2.1 u (F^(c (a+b x)))^n.jl


text:
tok  = s -> F.find_tokens(s, F.MD_TOKENS, F.MD_1C_TOKENS)
vfn  = s -> (t = tok(s); F.validate_footnotes!(t); t)
vh   = s -> (t = vfn(s); F.validate_start_of_line!(t, Franklin.MD_HEADER_OPEN); t)
fib  = s -> (t = vh(s); F.find_indented_blocks!(t, s); t)
fib2 = s -> (t = fib(s); F.filter_lr_indent!(t, s); t)

islr(t) = t.name == :LINE_RETURN && t.ss == "\n"
istok(t, n, s) = t.name == n && t.ss == s
isind(t) = t.name == :LR_INDENT && t.ss == "\n    "

##
## FIND_TOKENS
##

@testset "P:1:find-tok" begin
    t = raw"""
    @def v = 5
    @@da
        @@db
        @@
    @@
    $A$ and \[B\] and \com{hello} etc
    """ |> tok
    @test istok(t[1], :MD_DEF_OPEN, "@def")
    @test islr(t[2])
    @test istok(t[3], :DIV_OPEN,    "@@da")
    @test islr(t[4])
    @test istok(t[5], :DIV_OPEN, "@@db")
    @test islr(t[6])
    @test istok(t[7], :DIV_CLOSE, "@@")
    @test islr(t[8])
    @test istok(t[9], :DIV_CLOSE, "@@")
    @test islr(t[10])
    @test istok(t[11], :MATH_A, "\$")
    @test istok(t[12], :MATH_A, "\$")
    @test istok(t[13], :MATH_C_OPEN, "\\[")
    @test istok(t[14], :MATH_C_CLOSE, "\\]")
    @test istok(t[15], :LX_COMMAND, "\\com")
    @test istok(t[16], :LXB_OPEN, "{")
    @test istok(t[17], :LXB_CLOSE, "}")
    @test islr(t[18])
    @test istok(t[19], :EOS, "\n")

    # complement div with `-` or `_`
    t = raw"""A @@d-a B@@ C""" |> tok
    @test istok(t[1], :DIV_OPEN, "@@d-a")
    t = raw"""A @@d-1 B@@ C""" |> tok
    @test istok(t[1], :DIV_OPEN, "@@d-1")
end

@testset "P:1:ctok" begin
    # check that tokens at EOS close properly
    # NOTE: this was avoided before by the addition of a special char
    # to denote the end of the string but we don't do that anymore.
    # see `isexactly` and `find_tokens` in the fixed pattern case.
    t = "@@d ... @@" |> tok
    @test istok(t[1], :DIV_OPEN, "@@d")
    @test istok(t[2], :DIV_CLOSE, "@@")
    @test istok(t[3], :EOS, "@")
    t = "``` ... ```" |> tok
    @test istok(t[1], :CODE_TRIPLE, "```")
    @test istok(t[2], :CODE_TRIPLE, "```")
    @test istok(t[3], :EOS, "`")
    t = "<!--...-->" |> tok
    @test istok(t[1], :COMMENT_OPEN, "<!--")
    @test istok(t[2], :COMMENT_CLOSE, "-->")
    t = "~~~...~~~" |> tok
    @test istok(t[1], :ESCAPE, "~~~")
    @test istok(t[2], :ESCAPE, "~~~")
    t = "b `j`" |> tok
    @test istok(t[1], :CODE_SINGLE, "`")
    @test istok(t[2], :CODE_SINGLE, "`")
    @test istok(t[3], :EOS, "`")
end

##
## VALIDATE_FOOTNOTE!
##

@testset "P:1:val-fn" begin
    t = raw"""
    A [^B] and
    [^B]: etc
    """ |> vfn
    @test istok(t[1], :FOOTNOTE_REF, "[^B]")
    @test islr(t[2])
    @test istok(t[3], :FOOTNOTE_DEF, "[^B]:")
end

##
## VALIDATE_HEADERS!
##

@testset "P:1:val-hd" begin
    t = raw"""
    # A
    ## B
    and # C
    """ |> vh
    @test istok(t[1], :H1_OPEN, "#")
    @test islr(t[2])
    @test istok(t[3], :H2_OPEN, "##")
    @test islr(t[4])
    @test islr(t[5])
    @test istok(t[6], :EOS, "\n")
end

##
## FIND_INDENTED_BLOCKS!
##

@testset "P:1:fib" begin
    s = raw"""
    A

        B1
        B2
        B3
    C
    @@da
        @@db,dc:dd/de
            E
        @@
    @@
    E
        F
    G
    """
    t = s |> fib
    @test islr(t[1])
    @test isind(t[2]) && t[2].lno == 3
    @test isind(t[3]) && t[3].lno == 4
    @test isind(t[4]) && t[4].lno == 5 # B3
    @test islr(t[5])
    @test islr(t[6])
    @test istok(t[7], :DIV_OPEN, "@@da")
    @test isind(t[8]) && t[8].lno == 8
    @test istok(t[9], :DIV_OPEN, "@@db,dc:dd/de")
    @test isind(t[10]) && t[10].lno == 9  # in front of E
    @test isind(t[11]) && t[11].lno == 10 # in front of @@
    @test istok(t[12], :DIV_CLOSE, "@@")
    @test islr(t[13])
    @test istok(t[14], :DIV_CLOSE, "@@")
    @test islr(t[15])
    @test isind(t[16]) && t[16].lno == 13 # F
    @test islr(t[17])
    @test islr(t[18])
    @test istok(t[19], :EOS, "\n")

    t = s |> fib2
    @test length(t) == 19

    @test isind.([t[2], t[3], t[4]]) |> all
    @test islr.([t[8], t[10], t[11], t[16]]) |> all
end


text:
using Pkg
Pkg.instantiate()
Pkg.develop(PackageSpec(path=pwd()))
Pkg.build("EzXMLPatched")


text:
module KSVM

using Reexport
using LearnBase
using LearnBase.LossFunctions
using LearnBase.Penalties
@reexport using MLKernels

using Compat
using ArrayViews
using UnicodePlots
using Base.LinAlg: BlasReal, axpy!

import StatsBase: predict, fit, fit!, coef, nobs, model_response
import MLBase: classify, labelencode, labeldecode, groupindices
import UnicodePlots: scatterplot
import Base: print, show, convert, minimum
import LearnBase: nclasses, labels, classdistribution

export

    HingeLoss,
    L1HingeLoss,
    L2HingeLoss,
    SmoothedL1HingeLoss,
    ModifiedHuberLoss,

    EpsilonInsLoss,
    L1EpsilonInsLoss,
    L2EpsilonInsLoss,

    L1Penalty,
    L2Penalty,

    DualCD,
    DualCDWithShrinking,
    Pegasos,

    LinearSVMSpec,
    SVMSpec,
      SVMRiskSpec,
        CSVM,

    PrimalSolution,
    DualSolution,

    SVM,
      PrimalSVM,
        DensePrimalSVM,
        SparsePrimalSVM,
      DualSVM,
        DenseDualSVM,
        SparseDualSVM,

    SVC,
    SVR,

    labels,

    details,
    intercept,
    isconverged,
    iterations,
    features,
    targets,
    params,
    fval,
    nsv,
    svindex,
    xsv,
    ysv,
    predmodel,
    bias,

    predict,
    classify,
    accuracy,

    svm

#include("deprecate/bridge.jl")
include("common.jl")
include("io.jl")
include("primal.jl")
include("dual.jl")
include("svm_spec.jl")
include("svm_model.jl")
include("svc_svr.jl")
include("fit.jl")
include("interface.jl")
include("linear/solver/dualcd.jl")
include("linear/solver/dualcd_shrinking.jl")
include("linear/solver/pegasos.jl")

end


text:
#%% Uniform scaling
#%% ------------------------------------------

let FT = rFFT(nᵢ=(256, 50), pᵢ=(1.0, 0.5))

	fmap = rand(Bool, Grid(FT).nxi...) |> Rmap{FT}
	gfourier = rand(Bool, Grid(FT).nki...) |> Rfourier{FT}

	Dmap = DiagOp(fmap)
	Dfourier = DiagOp(fmap)

	Id2 = 2*LinearAlgebra.I

	@inferred Id2 * fmap
	@inferred Id2 * gfourier
	@inferred Id2 \ fmap
	@inferred Id2 \ gfourier

	@inferred Dmap * Id2
	@inferred Dfourier * Id2
	@inferred Dmap + Id2
	@inferred Dfourier + Id2
	@inferred Dmap - Id2
	@inferred Dfourier - Id2

	@inferred Id2 * Dmap
	@inferred Id2 * Dfourier
	@inferred Id2 + Dmap
	@inferred Id2 + Dfourier
	@inferred Id2 - Dmap
	@inferred Id2 - Dfourier

	# @test Id2 * fmap
	# @test Id2 * gfourier
	# @test Id2 \ fmap
	# @test Id2 \ gfourier

	# @test Dmap * Id2
	# @test Dfourier * Id2
	# @test Dmap + Id2
	# @test Dfourier + Id2
	# @test Dmap - Id2
	# @test Dfourier - Id2

	# @test Id2 * Dmap
	# @test Id2 * Dfourier
	# @test Id2 + Dmap
	# @test Id2 + Dfourier
	# @test Id2 - Dmap
	# @test Id2 - Dfourier

end

#%% SField Constructors
#%% ------------------------------------------

let FT = rFFT(nᵢ=(256, 50), pᵢ=(1.0, 0.5))

	@inferred Rmap{FT}()
	@inferred Rmap{FT}(0)
	@inferred Rfourier{FT}()
	@inferred Rfourier{FT}(0)

	matx = rand(Bool, Grid(FT).nxi...)
	matk = rand(Bool, Grid(FT).nki...)

	@inferred Rmap{FT}(matx)
	@inferred Rfourier{FT}(matk)

end

#%% Test using FT for converting between Rmap and Rfourier
#%% ------------------------------------------

let FT = rFFT(nᵢ=(256, 50), pᵢ=(1.0, 0.5))

	matx = rand(Float64, Grid(FT).nxi...)
	matk = FT*rand(Float64, Grid(FT).nxi...)
	fx = matx |> Rmap{FT}
	fk = matk |> Rfourier{FT}

	@inferred FT * fx
	@inferred FT * fk
	@inferred FT \ fx
	@inferred FT * fk

	@test all((FT * fx).k .== (FT * matx))
	@test all((FT \ fx).x .== matx)

	@test all((FT * fk).k .== matk)
	@test all((FT \ fk).x .== (FT \ matk))

end

##======================================================
pᵢ  = (1.0, 0.5) # periods
nᵢ   = (256, 256) # number of samples (left endpoint included)
FT  = rFFT{Float32,nᵢ,pᵢ,length(nᵢ)}
UFT = rFFTunitary{Float32,nᵢ,pᵢ,length(nᵢ)}

FT_plan = plan(FT)
UFT_plan = plan(UFT)

grid   = Grid(FT)
λ      = wavenumber(FT)
kfulli = frequencies(FT)
xfulli = pixels(FT)

fk = FT * rand(grid.nxi...)
fx = FT \ rand(Complex{Float32}, grid.nki...)

f1 = Rmap{FT}(fx)
f2 = Rfourier{FT}(fk)
f4 = Rfourier{FT}(f1)
f3 = Rmap{FT}(f2)

@test all(f1[!] .== FT_plan * fx)
@test all(f1[:] .== fx)
@test all(f2[!] .== fk)
@test all(f2[:] .== FT_plan \ fk)

@inferred dot(f3, f1)
# real(dot(cplan(FT) * f3[:x], cplan(FT) * f1[:x])) / Grid(FT).Ωk
# real(dot(f3[:x], f1[:x])) / Grid(FT).Ωx

f1 + f2
- 2 * f2

L = DiagOp(f1)
L * f2

L1 = DiagOp(f1)
L2 = DiagOp(f2)
L3 = sqrt(L2)

@test all(L1[!] .== FT_plan * fx)
@test all(L1[:] .== fx)
@test all(L2[!] .== fk)
@test all(L2[:] .== FT_plan \ fk)

L1 * f2
L2 * L1 * f2
L2 * L1

##======================================================
pᵢ  = (1.0,)
nᵢ   = (256,)
FT  = rFFT{Float64,nᵢ,pᵢ,length(nᵢ)}

grid   = Grid(FT)
λ      = wavenumber(FT)
kfulli = frequencies(FT)
xfulli = pixels(FT)

fk = FT * rand(grid.nxi...)
fx = FT \ rand(Complex{Float64}, grid.nki...)

f1 = Rmap{FT}(fx)
f2 = Rfourier{FT}(fk)

f4 = Rfourier{FT}(f1)
f3 = Rmap{FT}(f2)

f1 + f2
- 2 * f2

L = DiagOp(f1)
L * f2

## =========================================

nᵢ  = (256, 50) # number of samples (left endpoint included)
pᵢ  = (1.0, 0.5) # periods
Δxᵢ = tuple((pn[1]/pn[2] for pn in zip(pᵢ,nᵢ))...)

FFT1 = rFFT{Float64,nᵢ,pᵢ,length(nᵢ)}
FFT2 = rFFT(nᵢ=nᵢ, pᵢ=pᵢ)
FFT3 = rFFT(nᵢ=nᵢ, Δxᵢ=Δxᵢ)
FFTu1 = rFFTunitary{Float64,nᵢ,pᵢ,length(nᵢ)}
FFTu2 = rFFTunitary(nᵢ=nᵢ, pᵢ=pᵢ)
FFTu3 = rFFTunitary(nᵢ=nᵢ, Δxᵢ=Δxᵢ)

@test FFT1 == FFT2 == FFT3
@test FFTu1 == FFTu2 == FFTu3

grid1 = Grid(FFT1)
grid2 = Grid(FFT2)
grid3 = Grid(FFT3)
gridu1 = Grid(FFTu1)
gridu2 = Grid(FFTu2)
gridu3 = Grid(FFTu3)

@test grid1.ki == grid2.ki == grid3.ki
@test gridu1.ki == gridu2.ki == gridu3.ki

## =========================================

pᵢ  = (1.0, 0.5) # periods
nᵢ   = (256, 256) # number of samples (left endpoint included)
cFT  = cFFT{Float64,nᵢ,pᵢ,length(nᵢ)}
cUFT  = cFFTunitary{Float64,nᵢ,pᵢ,length(nᵢ)}
rFT  = rFFT{Float64,nᵢ,pᵢ,length(nᵢ)}
rUFT  = rFFTunitary{Float64,nᵢ,pᵢ,length(nᵢ)}

grid = Grid(cFT)

@inferred Grid(cFT)
@inferred wavenumber(cFT)
@inferred frequencies(cFT)
@inferred pixels(cFT)

cfx = rand(Complex{Float64}, grid.nxi...)
rfx = rand(Float64, grid.nxi...)

@test grid.nxi == grid.nki

@test sum(abs2, plan(cFT)  * cfx .- cplan(rFT)  * cfx) ≈ 0.0
@test sum(abs2, plan(cUFT) * cfx .- cplan(rUFT) * cfx) ≈ 0.0

@test sum(abs2, plan(cFT)  * cfx .- cplan(cFT)  * cfx) ≈ 0.0
@test sum(abs2, plan(cUFT) * cfx .- cplan(cUFT) * cfx) ≈ 0.0

@test sum(abs2, plan(rFT)  * rfx .- rplan(rFT)  * rfx) ≈ 0.0
@test sum(abs2, plan(rUFT) * rfx .- rplan(rUFT) * rfx) ≈ 0.0

@test sum(abs2, plan(rFT)  * rfx .- rplan(cFT)  * rfx) ≈ 0.0
@test sum(abs2, plan(rUFT) * rfx .- rplan(cUFT) * rfx) ≈ 0.0

@inferred adjoint(plan(cFT))
@inferred adjoint(plan(rFT))
@inferred transpose(plan(cFT))
@inferred transpose(plan(rFT))

## =========================================

pᵢ  = (1.0, 3.5) # periods
nᵢ   = (64,64) # number of samples (left endpoint included)
FFT  = rFFT{Float64,nᵢ,pᵢ,length(nᵢ)}
UFT  = rFFTunitary{Float64,nᵢ,pᵢ,length(nᵢ)}
grid = Grid(FFT)

@inferred Grid(FFT)
@inferred wavenumber(FFT)
@inferred frequencies(FFT)
@inferred pixels(FFT)

fk = FFT * randn(nᵢ...)
fx = FFT \ randn(Complex{Float64}, grid.nki...)

f1 = Rmap{FFT}(fx)
f2 = Rfourier{FFT}(fk)
f4 = Rfourier{FFT}(f1)
f3 = Rmap{FFT}(f2)

@inferred f1 + f2
@inferred - 2 * f2

L = DiagOp(f1)
L * f2

L1 = DiagOp(f1+1)
L2 = DiagOp(f2+1)
for f ∈ (f1, f2), L ∈ (L1, L2), M ∈ (L1, L2)
	@inferred L * f
	@inferred L * M * f
	@inferred L * M \ f
	@inferred L \ M \ f
	@inferred L \ M * f
	@test mean(mean(abs.(t)) for t in fielddata((L * M * f) - (L * (M * f)))) ≈ 0.0 atol=1e-5
	@test mean(mean(abs.(t)) for t in fielddata((L * M \ f) - (inv(M) * (inv(L) * f)))) ≈ 0.0 atol=1e-5
	@test mean(mean(abs.(t)) for t in fielddata((L \ M \ f) - (inv(M) * (inv(inv(L)) * f)))) ≈ 0.0 atol=1e-5
	@test mean(mean(abs.(t)) for t in fielddata((L \ M * f) - (inv(L) * (M * f)))) ≈ 0.0 atol=1e-5
end

## =========================================

pᵢ  = (1.0,)
nᵢ   = (64,)
FFT  = rFFT{Float64,nᵢ,pᵢ,length(nᵢ)}
UFT  = rFFTunitary{Float64,nᵢ,pᵢ,length(nᵢ)}
grid = Grid(FFT)

@inferred Grid(FFT)
@inferred wavenumber(FFT)
@inferred frequencies(FFT)
@inferred pixels(FFT)

fk = FFT * rand(nᵢ...)
fx = FFT \ rand(Complex{Float64}, grid.nki...)

f1 = Rmap{FFT}(fx)
f2 = Rfourier{FFT}(fk)
f4 = Rfourier{FFT}(f1)
f3 = Rmap{FFT}(f2)

@inferred f1 + f2
@inferred - 2 * f2

L = DiagOp(f1)
L * f2

L1 = DiagOp(f1+1)
L2 = DiagOp(f2+1)
for f ∈ (f1, f2), L ∈ (L1, L2), M ∈ (L1, L2)
	@inferred L * f
	@inferred L * M * f
	@inferred L * M \ f
	@inferred L \ M \ f
	@inferred L \ M * f
	@test mean(mean(abs.(t)) for t in fielddata((L * M * f) - (L * (M * f)))) ≈ 0.0 atol=1e-5
	@test mean(mean(abs.(t)) for t in fielddata((L * M \ f) - (inv(M) * (inv(L) * f)))) ≈ 0.0 atol=1e-5
	@test mean(mean(abs.(t)) for t in fielddata((L \ M \ f) - (inv(M) * (inv(inv(L)) * f)))) ≈ 0.0 atol=1e-5
	@test mean(mean(abs.(t)) for t in fielddata((L \ M * f) - (inv(L) * (M * f)))) ≈ 0.0 atol=1e-5
end


text:
module FITSIO

export FITS,
       HDU,
       ImageHDU,
       TableHDU,
       ASCIITableHDU,
       FITSHeader,
       read_key,
       write_key,
       read_header,
       get_comment,
       set_comment!,
       copy_section

import Base: getindex,
             setindex!,
             length,
             show,
             read,
             read!,
             write,
             close,
             ndims,
             size,
             getkey,
             haskey,
             keys,
             values

# Deal with compatibility issues.
using Printf
import Base: iterate, lastindex
# Libcfitsio submodule
include("libcfitsio.jl")

using .Libcfitsio

# There are a few direct `ccall`s to libcfitsio in this module. For this, we
# need a few non-exported things from Libcfitsio: the shared library handle,
# and a helper function for raising errors. TYPE_FROM_BITPIX is awkwardly
# defined in Libcfitsio, even though it is not used there.
import .Libcfitsio: libcfitsio,
                    fits_assert_ok,
                    fits_assert_isascii,
                    TYPE_FROM_BITPIX,
                    ContiguousAbstractArray

# HDU Types
abstract type HDU end

mutable struct ImageHDU <: HDU
    fitsfile::FITSFile
    ext::Int
end

mutable struct TableHDU <: HDU
    fitsfile::FITSFile
    ext::Int
end

mutable struct ASCIITableHDU <: HDU
    fitsfile::FITSFile
    ext::Int
end

# FITS
#
# The FITS type represents a FITS file. It holds a reference to a
# FITSFile object (basically the low-level CFITSIO pointer). It also
# holds a reference to all of the previously accessed HDU
# objects. This is so that only a single HDU object is created for
# each extension in the file. It also allows a FITS object to tell
# previously created HDUs about events that happen to the file, such
# as deleting extensions. This could be done by, e.g., setting ext=-1
# in the HDU object.
"""
    FITS(filename::String, mode::String="r")

Open or create a FITS file. `mode` can be one of `"r"` (read-only),
`"r+"` (read-write) or `"w"` (write). In "write" mode, any existing
file of the same name is overwritten.

A `FITS` object is a collection of "Header-Data Units" (HDUs) and
supports the following operations:

* `f[i]`: Return the `i`-th HDU.

* `f[name]` or `f[name, ver]`: Return the HDU containing the given the
  given EXTNAME (or HDUNAME) keyword (a String), and optionally the
  given EXTVER (or HDUVER) number (an Integer).

* Iteration:
  ```
  for hdu in f
      ...
  end
  ```
"""
FITS
mutable struct FITS
    fitsfile::FITSFile
    filename::AbstractString
    mode::AbstractString
    hdus::Dict{Int, HDU}

    # Hold on to memory if backed by a julia buffer
    memhandle::FITSMemoryHandle
    hidden::Any

    function FITS(filename::AbstractString, mode::AbstractString="r")
        f = (mode == "r"                      ? fits_open_file(filename, 0)    :
             mode == "r+" && isfile(filename) ? fits_open_file(filename, 1)    :
             mode == "r+"                     ? fits_create_file(filename)     :
             mode == "w"                      ? fits_create_file("!"*filename) :
             error("invalid open mode: $mode"))

        new(f, filename, mode, Dict{Int, HDU}(), FITSMemoryHandle(), nothing)
    end

    function FITS(data::Vector{UInt8}, mode::AbstractString="r", filename = "")
        @assert mode == "r"
        f, handle = fits_open_memfile(data, 0)
        new(f, filename, mode, Dict{Int, HDU}(), handle, data)
    end
end

"""
    FITS(f::Function, args...)

Apply the function `f` to the result of `FITS(args...)` and close the
resulting file descriptor upon completion.
"""
function FITS(f::Function, args...)
    io = FITS(args...)
    try
        f(io)
    finally
        close(io)
    end
end

"""
    FITSHeader(keys::Vector{String}, values::Vector, comments::Vector{String})

An in-memory representation of the header of an HDU. It stores the
(key, value, comment) information for each 80-character "card" in a header.

Note that this structure is not linked to a FITS file in any way; it
is just a convenient structure for storing the header contents after
reading from a file. (This is similar to how an `Array` returned by
`read(f[1])` is not linked to the FITS file `f`.)  Manipulating a
`FITSHeader` will therefore have no immediate impact on any file, even
if it was created by `read_header(::HDU)`.  You can, however, write a
`FITSHeader` to a file using the `write(::FITS, ...)` methods that
append a new HDU to a file.
"""
mutable struct FITSHeader
    keys::Vector{String}
    values::Vector{Any}
    comments::Vector{String}
    map::Dict{String, Int}

    function FITSHeader(keys::Vector{String}, values::Vector,
                        comments::Vector{String})
        if ((length(keys) != length(values)) ||
            (length(keys) != length(comments)))
            error("keys, values, comments must be same length")
        end
        map = Dict{String, Int}()
        for i in 1:length(keys)
          map[keys[i]] = i
        end
        new(keys, Vector{Any}(values), comments, map)
    end
end

include("fits.jl")  # FITS methods
include("header.jl")  # FITSHeader methods
include("image.jl")  # ImageHDU methods
include("table.jl")  # TableHDU & ASCIITableHDU methods

"""
    libcfitsio_version() -> VersionNumber

Return the version of the underlying CFITSIO library

# Example

```julia
julia> FITSIO.libcfitsio_version()
v"3.37.0"
```

"""
function libcfitsio_version(version=fits_get_version())
    # fits_get_version returns a float. e.g., 3.341f0. We parse that
    # into a proper version number. E.g., 3.341 -> v"3.34.1"
    v = Int(round(1000 * version))
    x = div(v, 1000)
    y = div(rem(v, 1000), 10)
    z = rem(v, 10)
    VersionNumber(x, y, z)
end

end # module


text:
## @lc code=start
using LeetCode

## Counting sort
function count_difference(nums::Vector{Int}, k::Int)::Int
    res, mi, ma = 0, minimum(nums), maximum(nums)
    table = Dict{Int, Int}(i => 0 for i in mi:ma)
    for i in nums
        table[i] += 1
    end
    for i in mi:(ma - k)
        res += (table[i] * table[i+k])
    end
    res
end
## @lc code=end

## @lc test=start
@testset "2006.count-number-of-pairs-with-absolute-difference-k.jl" begin
    @test count_difference([1,2,2,1], 1) == 4
    @test count_difference([1,3], 3) == 0
    @test count_difference([3,2,1,5,4], 2) == 3
end
## @lc test=end

## @lc info=start
# Date: 2022-02-09
# link: [solution to 2006](https://leetcode-cn.com/problems/count-number-of-pairs-with-absolute-difference-k/solution/2006-chai-de-jue-dui-zhi-wei-k-de-shu-du-q6k9/)
## @lc info=end


text:
# This script demonstrates how to use MERAKit.jl to find an approximation to the ground
# state of the critical Ising model.
using Random
using LinearAlgebra
using TensorKit
using MERAKit
# DemoTools holds some utility functions needed by this script, for creating Hamiltonians
# and storing and reading MERAs to/from disk. It exports nothing, so all calls need to be
# qualified as DemoTools.whatever.
include("demo_tools.jl")
using .DemoTools

function main()
    # Fix the random seed for the sake of reproducibility. MERA algorithms sometimes get
    # stuck in local minima, so you should usually try a few different seeds for the same
    # optimisation. This one is going to do well though.
    Random.seed!(1)

    # Create the local Hamiltonian term `h` for the critical Ising model. Use Z2 symmetric
    # tensors, and block two sites together, so that the local state space dimension is 4.
    # The blocking simply obviates the need for a trivial, unitary layer at the bottom of
    # the MERA.
    h = DemoTools.ising_hamiltonian(; symmetry = :Z2, block_size = 2)
    exact_energy = -4/pi
    V_physical = space(h, 1)

    # Choose the type of MERA to use.
    meratype = TernaryMERA
    layers = 3
    # It's typically beneficial to gradually increase the bond dimension during the
    # optimisation. We do this here using the following steps: Each element of
    # `V_virtual_steps` is a vector space for the virtual indices of the MERA, and we run
    # the optimisation for each of them, using as the starting point the result of the
    # previous step.
    # Note that for Z2 symmetric tensors we have to specify how the bond dimension is
    # divided between the symmetry sectors. It is a priori not obvious what the optimal
    # choice is, and although here it happens to be such that both sectors get the same
    # dimension, this isn't always the case. In general, finding the optimal proportions
    # requires manually trying different options.
    # Note also that we could use different vector spaces for the different indices in the
    # MERA, both the ones between layers and the ones within layers (between isometries and
    # disentanglers, making the disentanglers isometric instead of unitary). This would
    # probably be better for accuracy/speed, but using the same `V_virtual` everywhere works
    # well too, so for simplicity we stick to that.
    V_virtual_steps = (
        Z2Space(0 => 2, 1 => 2),
        Z2Space(0 => 3, 1 => 3),
        Z2Space(0 => 4, 1 => 4),
    )

    # Specify the optimisation method. We use here Riemannian LBFGS as described in
    # https://arxiv.org/abs/2007.03638
    method = :lbfgs
    verbosity = 2
    gradient_delta = 1e-7
    # In practice the above threshold for the gradient norm, which determines convergence,
    # will not be reached. Instead, we set the maximum number of iterations to be done. We
    # specify this separately for the different bond dimensions in `V_virtual_steps`.
    # Choosing these numbers is more art than science, but many different choices would work
    # roughly equally well.
    maxiter_steps = (1000, 500, 500)

    # Once the optimisation is done, the resulting MERA will be written to a file at this
    # path.
    path = "./demo_result.jld2"

    # To keep track of how the optimisation is progressing, we will every once in a while
    # print out the energy error, the entropies of the density matrices at different scales
    # in the MERA, and the scaling dimensions. For this we use the `finalize!` keyword
    # argument of `minimize_expectation`, which can be a function that is run at the end of
    # every iteration. `finalize!` could be used to do plenty of things, including modifying
    # the current MERA `m` during the optimisation, but we simply use it to print things out
    # every `checkpoint_frequency` iterations. See the documentation of MERAKit.jl for more
    # details on `finalize!`.
    checkpoint_frequency = 100
    function finalize!(m, energy, g, repnum)
        repnum % checkpoint_frequency != 0 && (return m, energy, g)
        rhoees = densitymatrix_entropies(m)
        scaldims = scalingdimensions(m)
        @info("Energy error: $(energy - exact_energy)")
        @info("Density matrix entropies: $rhoees")
        @info("Scaling dimensions: $scaldims")
        return m, energy, g
    end

    # The main loop, that for every vector space in `V_virtual_steps` optimises the MERA.
    # For the first one, we use a randomly initialised MERA as the starting point, for the
    # latter ones we expand the bond dimension by padding tensors with zeros.
    local m
    for i in 1:length(V_virtual_steps)
        V_virtual = V_virtual_steps[i]
        if i == 1
            @info("Creating a random MERA with bond dimension $(dim(V_virtual)).")
            # The lowest index is of the physical bond dimension, the others are virtual.
            Vs = (V_physical, ntuple(x -> V_virtual, layers-1)...)
            m = random_MERA(meratype, Float64, Vs)
        else
            @info("Expanding to bond dimension $(dim(V_virtual)).")
            for j in 2:layers
                m = expand_bonddim(m, j-1, V_virtual)
                m = expand_internal_bonddim(m, j, V_virtual)
            end
        end
        # Many more optimisation parameters can be specified, see the documentation on
        # `minimize_expectation`.
        pars = (
            gradient_delta = gradient_delta,
            maxiter = maxiter_steps[i],
            method = method,
            verbosity = verbosity
        )
        m = minimize_expectation(m, h, pars; finalize! = finalize!)
    end

    # Write the resulting MERA to disc. It can then be read using
    # `DemoTools.load_mera(path)`.
    DemoTools.store_mera(path, m)
    return m
end

# Run the actual script. The whole thing is in a function mainly to avoid polluting the
# global namespace.
main()


text:
module TestContrasts

using Base.Test
using DataFrames

d = DataFrame(x = @pdata( [:a, :b, :c, :a, :a, :b] ))

mf = ModelFrame(Formula(nothing, :x), d)

# Dummy coded contrasts by default:
@test ModelMatrix(mf).m == [1  0  0
                            1  1  0
                            1  0  1
                            1  0  0
                            1  0  0
                            1  1  0]
@test coefnames(mf) == ["(Intercept)"; "x: b"; "x: c"]

mmm = ModelMatrix(mf).m
setcontrasts!(mf, x = DummyCoding())
@test ModelMatrix(mf).m == mmm

setcontrasts!(mf, x = EffectsCoding())
@test ModelMatrix(mf).m == [1 -1 -1
                            1  1  0
                            1  0  1
                            1 -1 -1
                            1 -1 -1
                            1  1  0]
@test coefnames(mf) == ["(Intercept)"; "x: b"; "x: c"]

# change base level of contrast
setcontrasts!(mf, x = EffectsCoding(base = :b))
@test ModelMatrix(mf).m == [1  1  0
                            1 -1 -1
                            1  0  1
                            1  1  0
                            1  1  0
                            1 -1 -1]
@test coefnames(mf) == ["(Intercept)"; "x: a"; "x: c"]

# change levels of contrast
setcontrasts!(mf, x = EffectsCoding(levels = [:c, :b, :a]))
@test ModelMatrix(mf).m == [1  0  1
                            1  1  0
                            1 -1 -1
                            1  0  1
                            1  0  1
                            1  1  0]
@test coefnames(mf) == ["(Intercept)"; "x: b"; "x: a"]

# change levels and base level of contrast
setcontrasts!(mf, x = EffectsCoding(levels = [:c, :b, :a], base = :a))
@test ModelMatrix(mf).m == [1 -1 -1
                            1  0  1
                            1  1  0
                            1 -1 -1
                            1 -1 -1
                            1  0  1]
@test coefnames(mf) == ["(Intercept)"; "x: c"; "x: b"]

# Helmert coded contrasts
setcontrasts!(mf, x = HelmertCoding())
@test ModelMatrix(mf).m == [1 -1 -1
                            1  1 -1
                            1  0  2
                            1 -1 -1
                            1 -1 -1
                            1  1 -1]
@test coefnames(mf) == ["(Intercept)"; "x: b"; "x: c"]

# Mismatching types of data and contrasts levels throws an error:
@test_throws ArgumentError setcontrasts!(mf, x = EffectsCoding(levels = ["a", "b", "c"]))

# Missing data is handled gracefully, dropping columns when a level is lost
d[3, :x] = NA
mf_missing = ModelFrame(Formula(nothing, :x), d, contrasts = Dict(:x => EffectsCoding()))
@test ModelMatrix(mf_missing).m == [1 -1
                                    1  1
                                    1 -1
                                    1 -1
                                    1  1]
@test coefnames(mf_missing) == ["(Intercept)"; "x: b"]

# Things that are bad to do:
# Applying contrasts that only have a subset of data levels:
@test_throws ArgumentError setcontrasts!(mf, x = EffectsCoding(levels = [:a, :b]))
# Applying contrasts that expect levels not found in data:
@test_throws ArgumentError setcontrasts!(mf, x = EffectsCoding(levels = [:a, :b, :c, :d]))
# Asking for base level that's not found in data
@test_throws ArgumentError setcontrasts!(mf, x = EffectsCoding(base = :e))

# Manually specified contrasts
contrasts = [0  1
             -1 -.5
             1  -.5]
setcontrasts!(mf, x = ContrastsCoding(contrasts))
@test ModelMatrix(mf).m == [1  0  1
                            1 -1 -.5
                            1  1 -.5
                            1  0  1
                            1  0  1
                            1 -1 -.5]

# throw argument error if number of levels mismatches
@test_throws ArgumentError setcontrasts!(mf, x = ContrastsCoding(contrasts[1:2, :]))
@test_throws ArgumentError setcontrasts!(mf, x = ContrastsCoding(hcat(contrasts, contrasts)))

# contrasts types must be instaniated
@test_throws ArgumentError setcontrasts!(mf, x = DummyCoding)

end


text:
# Check Output Component
function cpntlib2D()
    CpntLib = ("Vz","Vx","Tzz","Txx","P","Txz","Wy");
    println("Output Component Library: $CpntLib")
end

function CheckCpnt2D(
    rec::Union{Void,receiver},
    Slices::Union{Void,UnitRange{Int64},StepRange{Int64,Int64},Array{Int64,1},Int64},
    OptCpnt::Union{Void,Array{String}})
    if !(rec == nothing && Slices == nothing)
        if OptCpnt == nothing
            error("Please provide output component")
        else
            #CpntLib : Component Library
            CpntLib = ("Vz","Vx","Tzz","Txx","P","Txz","Wy","V","Tii");

            nOptCpnt = length(OptCpnt);
            # Check if Output component is in Component Library
            bool_rslt = map(x->x in CpntLib, OptCpnt);
            num_rslt = zeros(nOptCpnt);
            for i in 1:nOptCpnt
                num_rslt[i] = bool_rslt[i]*1
            end
            if countnz(num_rslt)!=nOptCpnt
                for i in 1:nOptCpnt
                    if bool_rslt[i] == false
                        println("$(OptCpnt[i]) is not a valid output component.")
                    end
                end
                error("Check output component library by typing 'cpntlib2D()'.")
            else println("Output Components Check Pass!")
            end
        end
    end
end

function message(
    rec::Union{Void,receiver},
    RecDataPath::Union{String,Void},
    Slices::Union{Void,UnitRange{Int64},StepRange{Int64,Int64},Array{Int64,1},Int64},
    WFDataPath::Union{String,Void},
    OptCpnt::Union{Void,Array{String}})

    if rec == nothing
        println("Warning: No recordings are output")
    else
        if RecDataPath == nothing
            error("Please provide recording output path")
        else
            println("Receiver locations are successfully initiated")
        end
    end
    if Slices == nothing
        println("Warning: No wavefields are output")
    else
        if WFDataPath == nothing
            error("Please provide wavefields output path")
        else
            println("Wavefield slices are successfully initiated")
        end
    end
    CheckCpnt2D(rec,Slices,OptCpnt)
end

function message(
    Slices::Union{Void,UnitRange{Int64},StepRange{Int64,Int64},Array{Int64,1},Int64},
    WFDataPath::Union{String,Void})

    if Slices == nothing
        error("Please give time slices for model order reduction")
    else
        if WFDataPath == nothing
            error("Please provide wavefields output path")
        else
            println("Wavefield slices are successfully initiated")
        end
    end
end


text:
@testset "Sparse matrices" begin

  @testset "Hermite normal form" begin
    for i in 1:10
      r = 200
      c = 100
      A = matrix(FlintZZ, rand([0,0,0,0,0,0,0,0,0,0,1], r, c))
      As = SMat(A)
      @test sub(hnf(A), 1:c, 1:c) == fmpz_mat(hnf(As))
    end
  end
end


text:
##### ----- Power Law ----- #####
# Calculate probability density for a residence time τ given α of power law,
# minimum and maximum residence times τ₋ and τ₊
function E_powerLaw(τ, (q, α, τ₋, τ₊, k)) # Note - destructured arguments so that all pdfs have τ as first argument and p as second
    if(τ < τ₋)
        0
    elseif(τ > τ₊)
        0
    else
        ((τ^(-α) - τ₊^(-α))*(-α+1.0)/(τ₊^(-α+1.0)-τ₋^(-α+1.0)-τ₊^(-α)*(τ₊ - τ₋)*(1.0-α)))
    end
    # else
    #     (α-1)*τ₋^(α-1)*τ^(-α)
    # end
end

# Calculate rate of hyporheic exchange [T^-1] given a ratio of hyporheic volume to channel volume (V_frac)
# and parameters describing the pdf of the residence time distribution
function qCalc_powerLaw(α, τ₋, τ₊, V_frac)
    V_frac/quadgk(τ -> (τ₊^(-α+1) - τ^(-α+1)- (τ₊^(-α)*(τ₊-τ)*(-α+1)))/(τ₊^(-α+1) - τ₋^(-α+1)- (τ₊^(-α)*(τ₊-τ₋)*(-α+1))), τ₋, τ₊)[1]
    # V_frac/quadgk(τ -> (-τ^(-α+1))/(- τ₋^(-α+1)), τ₋, Inf)[1]
end

##### ----- Exponential ----- #####


text:
function assemble_system!(u, problem::FDGZProblem{TBC, TBC}) where {TBC<:DirichletBC}

    problem.F[1] =problem.ϵ * (u[2] - 2*u[1] + problem.left_bc.u_bc)/problem.Δx^2 +  u[1] * (1-u[1]^2);

    for i in 2:problem.n-1
        problem.F[i]= problem.ϵ * (u[i+1] - 2*u[i] +u[i-1])/problem.Δx^2 +u[i] * (1-u[i]^2);
    end
    problem.F[end] = problem.ϵ * (problem.right_bc.u_bc - 2*u[end] + u[end-1])/problem.Δx^2 +  u[end] * (1-u[end]^2);

    problem
end


text:
# Autogenerated wrapper script for ROCmLLVM_jll for x86_64-linux-musl-cxx03
export clang, libclang, libclang_cpp, liblto, llc, lld, llvm_config, llvm_mca, opt

using Zlib_jll
JLLWrappers.@generate_wrapper_header("ROCmLLVM")
JLLWrappers.@declare_library_product(libclang, "libclang.so.12git")
JLLWrappers.@declare_library_product(libclang_cpp, "libclang-cpp.so.12git")
JLLWrappers.@declare_library_product(liblto, "libLTO.so.12git")
JLLWrappers.@declare_executable_product(clang)
JLLWrappers.@declare_executable_product(llc)
JLLWrappers.@declare_executable_product(lld)
JLLWrappers.@declare_executable_product(llvm_config)
JLLWrappers.@declare_executable_product(llvm_mca)
JLLWrappers.@declare_executable_product(opt)
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_library_product(
        libclang,
        "lib/libclang.so",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libclang_cpp,
        "lib/libclang-cpp.so",
        nothing,
    )

    JLLWrappers.@init_library_product(
        liblto,
        "lib/libLTO.so",
        nothing,
    )

    JLLWrappers.@init_executable_product(
        clang,
        "bin/clang",
    )

    JLLWrappers.@init_executable_product(
        llc,
        "tools/llc",
    )

    JLLWrappers.@init_executable_product(
        lld,
        "bin/lld",
    )

    JLLWrappers.@init_executable_product(
        llvm_config,
        "tools/llvm-config",
    )

    JLLWrappers.@init_executable_product(
        llvm_mca,
        "tools/llvm-mca",
    )

    JLLWrappers.@init_executable_product(
        opt,
        "tools/opt",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Load Julia packages (libraries)

using Pkg, DrWatson

#@quickactivate "SR2StanPluto"
using StanSample
using StatisticalRethinking

# Define the Stan language model

stan8_5s = "
data{
  int N;
  vector[N] y;
}
parameters{
  real sigma;
  real a1;
  real a2;
}
model{
  real mu;
  a1 ~ normal(0, 10);
  a2 ~ normal(0, 10);
  mu = a1 + a2;
  y ~ normal( mu , sigma );
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

m8_5s = SampleModel("m8.5s", stan8_5s);

# Input data for cmdstan

m8_5_data = Dict("N" => 100, "y" => rand(Normal(0, 1), 100));

# Sample using cmdstan

rc8_5s = stan_sample(m8_5s, data=m8_5_data);

rethinking = "
       mean   sd   5.5% 94.5% n_eff Rhat
a1    -0.08 7.15 -11.34 11.25  1680    1
a2    -0.05 7.15 -11.37 11.19  1682    1
sigma  0.90 0.07   0.81  1.02  2186    1
";

# Describe the draws
if success(rc8_5s)
  part8_5s = read_samples(m8_5s, :particles)
  part8_5s |> display
end