text:
using QuantumControlBase
using QuantumPropagators
using Documenter
using Pkg

DocMeta.setdocmeta!(
    QuantumControlBase,
    :DocTestSetup,
    :(using QuantumControlBase);
    recursive=true
)

PROJECT_TOML = Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))
VERSION = PROJECT_TOML["version"]
NAME = PROJECT_TOML["name"]
AUTHORS = join(PROJECT_TOML["authors"], ", ") * " and contributors"
GITHUB = "https://github.com/JuliaQuantumControl/QuantumControlBase.jl"

println("Starting makedocs")

makedocs(;
    authors=AUTHORS,
    sitename="QuantumControlBase.jl",
    modules=[QuantumControlBase],
    format=Documenter.HTML(;
        prettyurls=true,
        canonical="https://juliaquantumcontrol.github.io/QuantumControlBase.jl",
        assets=String[],
        footer="[$NAME.jl]($GITHUB) v$VERSION docs powered by [Documenter.jl](https://github.com/JuliaDocs/Documenter.jl)."
    ),
    pages=[
        "Home" => "index.md",
        hide("QuantumPropagators" => "quantumpropagators.md"),
        "History" => "history.md",
    ]
)

println("Finished makedocs")

deploydocs(; repo="github.com/JuliaQuantumControl/QuantumControlBase.jl")


text:
using LinearAlgebra, Plots, Statistics, TaylorSeries, Printf

path1=realpath(dirname(@__FILE__)*"/..")
include(string(path1,"/src/TaylorBigF.jl"))

# 3.2 Function Approximation

p=2000;setprecision(p); setrounding(BigFloat, RoundUp);
precision(BigFloat)==p
L=BigFloat("1"); n=BigInt(201); iN=convert(Int64,n); dx=BigFloat("2")*L/(n-BigFloat("1"));
x=convert(Array{BigFloat},-L:dx:L)
length(x)==n
# x=-L:dx:L
# x=convert(Array{BigFloat},(rand(n).-1/2)/10)
I=zeros(BigFloat, n, n);for i=1:n I[i,i]=BigFloat("1") end

V=TaylorBigF.calc_Vandermonde(x,n,n)

function func(x) sin(x) end

f=func.(x)
a=V\f
# plot(a)
# inv_V_comp = inv(V)
# a=inv_V_comp*f

# inv_V_analyt,det_V_analyt=TaylorBigF.calc_Vandermonde_inv_formula(x,n);
# a=inv_V_analyt*f

# inv_V_comp = V\I
# a=inv_V_comp*f

# a=convert(Array{BigFloat},convert(Array{Float64},V)\convert(Array{Float64},f))

xi=(x.+dx/BigFloat("2"))[1:end-1];
Vi=TaylorBigF.calc_Vandermonde(xi,n-1,n);

fi_comp=Vi*a;
fi_analyt=func.(xi)
diff_f=maximum(abs.(f-V*a));@sprintf("diff_f=%1.3e",diff_f)
diff_fi=maximum(abs.(fi_analyt-fi_comp));@sprintf("diff_fi=%1.3e",diff_fi)

# l1=30;fF=convert(Array{Float64},fi_comp);xF=convert(Array{Float64},x);xiF=convert(Array{Float64},xi);
# ticksY = round.(range(minimum(fF[1:l1]),stop=maximum(fF[1:l1]),length = 5),digits=2);
# ticksX = round.(range(minimum(xiF[1:l1]),stop=maximum(xiF[1:l1]),length = 5),digits=2);
# scatter(xF[1:l1],f[1:l1],color=:black,markersize=3,xlabel="x",
# ylabel="f(x)",yticks=ticksY,xticks=ticksX,label="f(x) analytical",
# legend=:bottomright,markershape=:rtriangle)
# scatter!(xiF[1:l1],fi_comp[1:l1],color=:black,markersize=3,
# label=string("f(x) numerical p=",p),markershape=:circle)# circle,rect,cross

# savefig(string("f__.pdf"))

ext: jl


text:
module KernelGradients

import KernelAbstractions: Kernel, CPUCTX, CPUCtx
import Cassette
import Enzyme
using Requires

@inline function Cassette.overdub(::CPUCtx, ::typeof(Enzyme.autodiff_deferred), f, annotation::Enzyme.Annotation, args...)
    f′ = (args...) -> (Base.@_inline_meta; Cassette.overdub(CPUCTX, f, args...))
    Enzyme.autodiff_deferred(f′, annotation, args...)
end

@inline function Cassette.overdub(::CPUCtx, ::typeof(Enzyme.autodiff_deferred), f, args...)
    f′ = (args...) -> (Base.@_inline_meta; Cassette.overdub(CPUCTX, f, args...))
    Enzyme.autodiff_deferred(f′, args...)
end

function Enzyme.autodiff(kernel::Kernel{<:Any, <:Any, <:Any, Fun}) where Fun
    function df(ctx, args...)
        Enzyme.autodiff_deferred(kernel.f, Enzyme.Const, ctx, args...)
    end
    similar(kernel, df)
end

function __init__()
    @require CUDAKernels="72cfdca4-0801-4ab0-bf6a-d52aa10adc57" include("cuda.jl")
    @require ROCKernels="7eb9e9f0-4bd3-4c4c-8bef-26bd9629d9b9" include("roc.jl")
end

end # module


text:
using OpenTelemetrySDK
using OpenTelemetryAPI
using Test
using Base.Threads

@testset "OpenTelemetrySDK.jl" begin
    include("trace.jl")
end


text:
function preprocessST(H, g, params::TparamsST, n1, n2)
    ζ = params.ζ

    return  PDataST(H, g, ζ, true)
end


text:
@doc raw"""
    postproc(model::AshkinTeller, param::Parameter, obs::MCObservableSet)

# Observables to be calculated
In the following, $m$ is total magnetization per site and $\epsilon$ is total energy per site.

- `"Binder Ratio"`
    - $R := \frac{\left \langle m^4 \right \rangle}{\left \langle m^2 \right\rangle^2}$
- `"Susceptibility"`
    - $\chi := \frac{N}{T}\left(\left\langle m^2\right\rangle\right)$
- `"Connected Susceptibility"`
    - $\frac{N}{T}\left(\left\langle m^2\right\rangle - \left\langle |m| \right\rangle^2\right)$
- `"Binder Ratio sigma"`
    - Binder ratio with respct to ``\sigma``
- `"Susceptibility sigma"`
    - Susceptibility with respct to ``\sigma``
- `"Connected Susceptibility sigma"`
    - Connected susceptibility with respct to ``\sigma``
- `"Binder Ratio tau"`
    - Binder ratio with respct to ``\tau``
- `"Susceptibility tau"`
    - Susceptibility with respct to ``\tau``
- `"Connected Susceptibility tau"`
    - Connected susceptibility with respct to ``\tau``
- `"Specific Heat"`
    - $\frac{N}{T^2}\left(\left\langle \epsilon^2\right\rangle - \left\langle \epsilon \right\rangle^2\right)$
"""
function postproc(model::AshkinTeller, param::Parameter, obs::MCObservableSet)
    nsites = numsites(model)
    T = convert(Float64, param["T"])
    beta = 1.0/T

    jk = jackknife(obs)

    jk["Binder Ratio"] = jk["|Magnetization|^4"] / (jk["|Magnetization|^2"]^2)
    jk["Susceptibility"] = (nsites*beta)*jk["|Magnetization|^2"]
    jk["Connected Susceptibility"] = (nsites*beta)*(jk["|Magnetization|^2"] - jk["|Magnetization|"]^2)

    jk["Binder Ratio sigma"] = jk["Magnetization sigma^4"] / (jk["Magnetization sigma^2"]^2)
    jk["Susceptibility sigma"] = (nsites*beta)*jk["Magnetization sigma^2"]
    jk["Connected Susceptibility sigma"] = (nsites*beta)*(jk["Magnetization sigma^2"] - jk["|Magnetization sigma|"]^2)

    jk["Binder Ratio tau"] = jk["Magnetization tau^4"] / (jk["Magnetization tau^2"]^2)
    jk["Susceptibility tau"] = (nsites*beta)*jk["Magnetization tau^2"]
    jk["Connected Susceptibility tau"] = (nsites*beta)*(jk["Magnetization tau^2"] - jk["|Magnetization tau|"]^2)

    jk["Specific Heat"] = (nsites*beta*beta)*(jk["Energy^2"] - jk["Energy"]^2)
    return jk
end


text:
"""
    VRTCreate(int,
              int) -> VRTDatasetH
"""
function vrtcreate(arg1::Integer, arg2::Integer)
    checknull(ccall((:VRTCreate, libgdal), Ptr{VRTDatasetH}, (Cint, Cint), arg1, arg2))
end

"""
    VRTFlushCache(VRTDatasetH) -> void
"""
function vrtflushcache(arg1::Ref{VRTDatasetH})
    ccall((:VRTFlushCache, libgdal), Void, (Ptr{Void},), arg1)
end

"""
    VRTSerializeToXML(VRTDatasetH,
                      const char *) -> CPLXMLNode *
"""
function vrtserializetoxml(arg1::Ref{VRTDatasetH}, arg2)
    ccall((:VRTSerializeToXML, libgdal), Ptr{CPLXMLNode}, (Ptr{Void}, Cstring), arg1, arg2)
end

"""
    VRTAddBand(VRTDatasetH,
               GDALDataType,
               char **) -> int
"""
function vrtaddband(arg1::Ref{VRTDatasetH}, arg2::GDALDataType, arg3)
    ccall((:VRTAddBand, libgdal), Cint, (Ptr{Void}, GDALDataType, Ptr{Cstring}), arg1, arg2, arg3)
end

"""
    VRTAddSource(VRTSourcedRasterBandH,
                 VRTSourceH) -> CPLErr
"""
function vrtaddsource(arg1::Ref{VRTSourcedRasterBandH}, arg2::Ref{VRTSourceH})
    ccall((:VRTAddSource, libgdal), CPLErr, (Ptr{Void}, Ptr{Void}), arg1, arg2)
end

"""
    VRTAddSimpleSource(VRTSourcedRasterBandH,
                       GDALRasterBandH,
                       int,
                       int,
                       int,
                       int,
                       int,
                       int,
                       int,
                       int,
                       const char *,
                       double) -> CPLErr
"""
function vrtaddsimplesource(arg1::Ref{VRTSourcedRasterBandH}, arg2::Ref{GDALRasterBandH}, arg3::Integer, arg4::Integer, arg5::Integer, arg6::Integer, arg7::Integer, arg8::Integer, arg9::Integer, arg10::Integer, arg11, arg12::Real)
    ccall((:VRTAddSimpleSource, libgdal), CPLErr, (Ptr{Void}, Ptr{Void}, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cstring, Cdouble), arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9, arg10, arg11, arg12)
end

"""
    VRTAddComplexSource(VRTSourcedRasterBandH,
                        GDALRasterBandH,
                        int,
                        int,
                        int,
                        int,
                        int,
                        int,
                        int,
                        int,
                        double,
                        double,
                        double) -> CPLErr
"""
function vrtaddcomplexsource(arg1::Ref{VRTSourcedRasterBandH}, arg2::Ref{GDALRasterBandH}, arg3::Integer, arg4::Integer, arg5::Integer, arg6::Integer, arg7::Integer, arg8::Integer, arg9::Integer, arg10::Integer, arg11::Real, arg12::Real, arg13::Real)
    ccall((:VRTAddComplexSource, libgdal), CPLErr, (Ptr{Void}, Ptr{Void}, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cdouble, Cdouble, Cdouble), arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9, arg10, arg11, arg12, arg13)
end

"""
    VRTAddFuncSource(VRTSourcedRasterBandH,
                     VRTImageReadFunc,
                     void *,
                     double) -> CPLErr
"""
function vrtaddfuncsource(arg1::Ref{VRTSourcedRasterBandH}, arg2::Ref{VRTImageReadFunc}, arg3, arg4::Real)
    ccall((:VRTAddFuncSource, libgdal), CPLErr, (Ptr{Void}, Ptr{Void}, Ptr{Void}, Cdouble), arg1, arg2, arg3, arg4)
end


text:
"""動力学関連
"""
module Dynamics

using LinearAlgebra

"""慣性行列"""
function M!(q::Vector{Float64}, out::Matrix{Float64})
    ccall(
        (:M, "o/soft_robot/derived/ikko_dake/eqs/c_so/M.so"),
        Cvoid,
        (Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
        q[1], q[2], q[3], out
    )
end

"""慣性行列"""
function M(q::Vector{Float64})
    Z = Matrix{Float64}(undef, 3, 3)
    M!(q, Z)
    return Z
end

"""慣性行列の逆行列"""
function invM!(q::Vector{Float64}, out::Matrix{Float64})
    ccall(
        (:invM, "o/soft_robot/derived/ikko_dake/eqs/c_so/invM.so"),
        Cvoid,
        (Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
        q[1], q[2], q[3], out
    )
end

"""慣性行列の逆行列"""
function invM(q::Vector{Float64})
    Z = Matrix{Float64}(undef, 3, 3)
    invM!(q, Z)
    return Z
end

function C!(q::Vector{Float64}, q_dot::Vector{Float64}, out::Matrix{Float64})
    ccall(
        (:C, "o/soft_robot/derived/ikko_dake/eqs/c_so/C.so"),
        Cvoid,
        (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
        q[1], q[2], q[3], q_dot[1], q_dot[2], q_dot[3], out
    )
end

"""コリオリ＋遠心力行列"""
function C(q::Vector{Float64}, q_dot::Vector{Float64})
    Z = Matrix{Float64}(undef, 3, 3)

    C!(q, q_dot, Z)

    return Z
end

"""Gのラッパー"""
function G!(q::Vector{Float64}, out::Vector{Float64})
    ccall(
        (:G, "o/soft_robot/derived/ikko_dake/eqs/c_so/G.so"),
        Cvoid,
        (Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
        q[1], q[2], q[3], out
    )
end

"""重力行列"""
function G(q::Vector{Float64})
    Z = Vector{Float64}(undef, 3)

    G!(q, Z)

    return Z
end

const K = diagm([1700.0, 1700.0, 1700.0])  # 剛性行列の真値
const D = diagm([110.0, 110.0, 110.0])  # 減衰行列の真値

const uncertain_K = diagm([1020.0, 1020.0, 1020.0])  # 不確かな剛性行列
const uncertain_D = diagm([77.0, 77.0, 77.0])  # 不確かな減衰行列

const αh = 23.705
const βh = 1.7267
const γh = -42.593

function h_dot(h::T, q::T, q_dot::T) where T
    q * (αh - (βh * sign(q_dot*h) + γh)*abs(h))
end

"""ヒステリシスベクトルの更新式"""
function H_dot(H::Vector{T}, q::Vector{T}, q_dot::Vector{T}) where T
    [
        h_dot(H[1], q[1], q_dot[1]),
        h_dot(H[2], q[2], q_dot[2]),
        h_dot(H[3], q[3], q_dot[3]),
    ]
end

"""加速度"""
function calc_q_dot_dot(
    τ::Vector{T}, q::Vector{T}, q_dot::Vector{T},
    H::Vector{T},
    ) where T
    inv(M(q)) * (τ .- (C(q, q_dot) .+ D)*q_dot .- K*q .- G(q) .- H) |> vec
end

"""加速度

outに加速度ベクトルを書き出す
"""
function q_dot_dot!(
    τ::Vector{Float64}, q::Vector{Float64}, q_dot::Vector{Float64},
    H::Vector{Float64}, out::Vector{Float64}
    )
    ccall(
        (:q_dot_dot, "o/soft_robot/derived/ikko_dake/eqs/c_so/q_dot_dot.so"),
        Cvoid,
        (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
        H[1], H[2], H[3], q[1], q_dot[1], q[2], q_dot[2], q[3], q_dot[3], τ[1], τ[2], τ[3], out
    )
end

end

# using LinearAlgebra
# using .Dynamics

# m = Dynamics.invM([0.0, 0.01, 0.0])
# m2 = Dynamics.invM([0.0, 0.02, 0.0])
# println((m .- m2) ./ m)
# # e,u=eigen(m)
# # e


text:
using ArgParse
s = ArgParseSettings()
@add_arg_table s begin
    "--part2"
    help = "Do part 2"
    action = :store_true
    "--test"
    help = "Run test data"
    action = :store_true
end
parsed_args = parse_args(ARGS, s)

test_suffix = parsed_args["test"] ? "_test" : ""
input = readlines("input/day_17$test_suffix.txt")
txt = split(input[1][14:end], ", ")

function get_min_max(rng)
    mn = [parse(Int16, p) for p in split(rng[3:end],"..")]
    (min = mn[1], max = mn[2])
end

bd = map( get_min_max, txt)
bounds = (x = bd[1], y = bd[2])

# target area: x=20..30, y=-10..-5
function hitsTarget(vel)
    pos = (x=0, y=0)
    ymax = 0
    while (pos.x < bounds.x.max && pos.y >= bounds.y.min)
        pos = (x= (pos.x + vel.x), y= (pos.y + vel.y))
        vel = (x= vel.x > 0 ? vel.x - 1 : 0, y= vel.y -1)
        ymax = max(ymax, pos.y)
        if pos.x >= bounds.x.min && pos.y >= bounds.y.min &&
            pos.x <= bounds.x.max && pos.y <= bounds.y.max
            return (true, false, ymax)
        end
    end
   (false, (pos.x >= bounds.x.max) || vel.x==0, ymax)
end

function part1()
    # 2 * bounds.min.x = (nn+n) >
    minx = Int16(floor(sqrt(2 * bounds.x.min)))
    maxx = bounds.x.max
    ymax = 0
    for vx in minx:maxx
        vy  = 0
        while true
            (hit, over, ym) = hitsTarget((x=vx, y=vy))
            if hit && ym > ymax
                ymax = ym
            end
            if vy > 1000
                break
            end
            vy +=1
        end
    end
    ymax
end

function part2()
    minx = Int16(floor(sqrt(2 * bounds.x.min)))
    maxx = bounds.x.max
    hits = 0
    for vx in minx:maxx
        vy  = bounds.y.min
        while true
            (hit, _, _) = hitsTarget((x=vx, y=vy))
            if hit
                hits += 1
            end
            if vy > 1000
                break
            end
            vy +=1
        end
    end
    hits
end

p = parsed_args["part2"] ? part2 : part1
answer = p()
println("The answer is: $answer")


text:
"""
Resource generic data structure
https://github.com/frictionlessdata/datapackage-jl#resource
"""

mutable struct Resource
    descriptor::Dict
    strict::Bool
    name::String
    path::String
    profile::String
    dialect::Dict
    schema::Schema
    errors::Array{PackageError}

    function Resource(d::Dict ; strict::Bool=false)
        schema = haskey(d, "schema") ?
            Schema(d["schema"], strict) : Schema()
        name = haskey(d, "name") ?
            d["name"] : nothing
        path = haskey(d, "path") ?
            d["path"] : nothing
        profile = haskey(d, "profile") ?
            d["profile"] : nothing
        dialect = haskey(d, "dialect") ?
            d["dialect"] : Dict()

        new(d, strict, name, path, profile, dialect, schema, [])
    end

    function Resource(path::String ; strict::Bool=false, name::String=nothing)
        name = isempty(name) ? path.split('/')[-1] : name
        new(
            Dict(), strict, name, path, "tabular-data-resource", Dict(), Schema(), []
        )
    end
end

function fields_to_descriptor(s::Schema)
    f_dict = []
    for f in s.fields
        push!(f_dict, Dict(
            "name" => f.name,
            "type" => f.typed
        ))
    end
    Dict("fields" => f_dict)
end

function get_table(r::Resource)
    if is_empty(r.schema)
        s = Schema()
        t = Table(r.path)
        tr = read(t, cast=false)
        infer(s, tr, t.headers)
        s.errors = []
        s.descriptor = fields_to_descriptor(s)
        validate(s, r.strict)
        r.schema = t.schema = s
        t
    else
        Table(r.path, r.schema)
    end
end

function read(r::Resource)
    if r.profile == "tabular-data-resource"
        t = get_table(r)
        read(t, cast=false)
    else
        throw(ErrorException("Not supported"))
    end
end


text:
using Printf

# Convert some native Julia data types to MathML
const _INT_FMT=String("%i")
const _FLT_FMT=String("%.2f")

_num_format(::Type{<:Integer}) = _INT_FMT
_num_format(::Type{<:Rational}) = _INT_FMT
_num_format(::Type{<:AbstractFloat}) = _FLT_FMT
_num_format(::Type{Complex{T}}) where T <: Integer = _INT_FMT
_num_format(::Type{Complex{T}}) where T <: Rational = _INT_FMT
_num_format(::Type{Complex{T}}) where T <: AbstractFloat = _FLT_FMT
_num_format(x) = _FLT_FMT

function num_to_expr(x::Integer; fmt::String=_INT_FMT)
    if x >= 0
        return NUM(@eval @sprintf($fmt, $x))
    else
        return Group(0, OPID('-'), NUM(@eval @sprintf($fmt, $(-x))))
    end
end

function num_to_expr(x::AbstractFloat; fmt::String=_FLT_FMT)
    if x >= 0
        return NUM(@eval @sprintf($fmt, $x))
    else
        return Group(0, OPID('-'), NUM(@eval @sprintf($fmt, $(-x))))
    end
end

function num_to_expr(x::Rational{T}; fmt::String=_num_format(T)) where T
    x == 0 && return NUM("0")
    r = Int(sign(x))
    if (x < 0) x = -x end
    if x.den != 1
        if r > 0
            return Frac(num_to_expr(x.num; fmt), num_to_expr(x.den; fmt))
        else
            return Group(0, OPID('-'), Frac(num_to_expr(x.num; fmt), num_to_expr(x.den; fmt)))
        end
    else
        if r > 0
            return num_to_expr(x.num; fmt)
        else
            return Group(0, OPID('-'), num_to_expr(x.num; fmt))
        end
    end
end

function num_to_expr(x::Complex{T}; fmt::String=_num_format(T)) where T
    return Group(0, num_to_expr(x.re; fmt), OPID('+'), num_to_expr(x.im; fmt), ID("𝑖"))
end

function vector_to_expr(
    x::Vector{T};
    mtype="pmatrix",
    fmt::String=_num_format(T)
) where T<:Number

    N = length(x)
    v = CompoExpr[]

    for (i, a) in enumerate(x)
        push!(v, num_to_expr(a; fmt))
        i < N && push!(v, Newline())
    end
    return Environment(mtype, v)
end

function matrix_to_expr(
    x::Matrix{T};
    mtype="pmatrix",
    fmt::String=_num_format(T)
) where T<:Number

    sz, N = size(x), length(x)
    v = CompoExpr[]

    for i in 1:sz[1], j in 1:sz[2]
        push!(v, num_to_expr(x[i,j]; fmt))
        if j < sz[2]
            push!(v, Ampersand())
        elseif i < sz[1]
            push!(v, Newline())
        end
    end
    return Environment(mtype, v)
end

# for LinearAlgebra.Adjoint type matrix
matrix_to_expr(x::AbstractMatrix) = matrix_to_expr(Matrix{eltype(x)}(x))

# Convert some native Julia data types to TeX/LaTeX
num_to_tex(x::Integer; fmt::String=_INT_FMT) = @eval @sprintf($fmt, $x)
num_to_tex(x::AbstractFloat; fmt::String=_FLT_FMT) = @eval @sprintf($fmt, $x)

function num_to_tex(x::Rational{T}; fmt::String=_num_format(T)) where T
    x == 0 && return "0"
    r = Int(sign(x))
    if (x < 0) x = -x end
    if x.den != 1
        if r > 0
            return "\\frac{$(num_to_tex(x.num; fmt))}{$(num_to_tex(x.den; fmt))}"
        else
            return "-\\frac{$(num_to_tex(x.num; fmt))}{$(num_to_tex(x.den; fmt))}"
        end
    else
        if r > 0
            return "$(num_to_tex(x.num; fmt))"
        else
            return "-$(num_to_tex(x.num; fmt))"
        end
    end
end

function num_to_tex(x::Complex{T}; fmt::String=_num_format(T)) where T
    return "$(num_to_tex(x.re; fmt))+$(num_to_tex(x.im; fmt))i"
end

function vector_to_tex(
    x::Vector{T};
    mtype::String="pmatrix",
    fmt::String=_num_format(T)
) where T<:Number

    N = length(x)
    s = "\\begin{$(mtype)}"

    for (i, a) in enumerate(x)
        s *= num_to_tex(a; fmt)
        if i < N
            s *= "\\\\"
        end
    end
    return s*"\\end{$(mtype)}"
end

function matrix_to_tex(
    x::Matrix{T};
    mtype::String="pmatrix",
    fmt::String=_num_format(T)
) where T<:Number

    sz, N = size(x), length(x)
    s = "\\begin{$(mtype)}"

    for i in 1:sz[1], j in 1:sz[2]
        s *= num_to_tex(x[i,j]; fmt)
        if j < sz[2]
            s *= "&"
        elseif i < sz[1]
            s *= "\\\\"
        end
    end
    return s*"\\end{$(mtype)}"
end
# for LinearAlgebra.Adjoint type matrix
matrix_to_tex(x::AbstractMatrix) = matrix_to_tex(Matrix{eltype(x)}(x))


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

    #= ::Package:: =#

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

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*x^m*(a+b*F^(e+f*x))^p=#

    #= ::Subsection:: =#
    #=p>0=#

    #= ::Subsection::Closed:: =#
    #=p<0=#

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

    @test_int [1/(a + b*ℯ^(c - d*x)), x, 4, x/a + log(a + b*ℯ^(c - d*x))/(a*d)]
    @test_int [1/(a + b*ℯ^(-c - d*x)), x, 4, x/a + log(a + b*ℯ^(-c - d*x))/(a*d)]

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

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

    @test_int [x^3/(a + b*ℯ^(c + d*x))^3, x, 26, (3*x^2)/(2*a^3*d^2) - (3*x^2)/(2*a^2*d^2*(a + b*ℯ^(c + d*x))) - (3*x^3)/(2*a^3*d) + x^3/(2*a*d*(a + b*ℯ^(c + d*x))^2) + x^3/(a^2*d*(a + b*ℯ^(c + d*x))) + x^4/(4*a^3) - (3*x*log(1 + (b*ℯ^(c + d*x))/a))/(a^3*d^3) + (9*x^2*log(1 + (b*ℯ^(c + d*x))/a))/(2*a^3*d^2) - (x^3*log(1 + (b*ℯ^(c + d*x))/a))/(a^3*d) - (3*PolyLog(2, -((b*ℯ^(c + d*x))/a)))/(a^3*d^4) + (9*x*PolyLog(2, -((b*ℯ^(c + d*x))/a)))/(a^3*d^3) - (3*x^2*PolyLog(2, -((b*ℯ^(c + d*x))/a)))/(a^3*d^2) - (9*PolyLog(3, -((b*ℯ^(c + d*x))/a)))/(a^3*d^4) + (6*x*PolyLog(3, -((b*ℯ^(c + d*x))/a)))/(a^3*d^3) - (6*PolyLog(4, -((b*ℯ^(c + d*x))/a)))/(a^3*d^4)]
    @test_int [x^2/(a + b*ℯ^(c + d*x))^3, x, 23, x/(a^3*d^2) - x/(a^2*d^2*(a + b*ℯ^(c + d*x))) - (3*x^2)/(2*a^3*d) + x^2/(2*a*d*(a + b*ℯ^(c + d*x))^2) + x^2/(a^2*d*(a + b*ℯ^(c + d*x))) + x^3/(3*a^3) - log(a + b*ℯ^(c + d*x))/(a^3*d^3) + (3*x*log(1 + (b*ℯ^(c + d*x))/a))/(a^3*d^2) - (x^2*log(1 + (b*ℯ^(c + d*x))/a))/(a^3*d) + (3*PolyLog(2, -((b*ℯ^(c + d*x))/a)))/(a^3*d^3) - (2*x*PolyLog(2, -((b*ℯ^(c + d*x))/a)))/(a^3*d^2) + (2*PolyLog(3, -((b*ℯ^(c + d*x))/a)))/(a^3*d^3)]
    @test_int [x^1/(a + b*ℯ^(c + d*x))^3, x, 15, -(1/(2*a^2*d^2*(a + b*ℯ^(c + d*x)))) - (3*x)/(2*a^3*d) + x/(2*a*d*(a + b*ℯ^(c + d*x))^2) + x/(a^2*d*(a + b*ℯ^(c + d*x))) + x^2/(2*a^3) + (3*log(a + b*ℯ^(c + d*x)))/(2*a^3*d^2) - (x*log(1 + (b*ℯ^(c + d*x))/a))/(a^3*d) - PolyLog(2, -((b*ℯ^(c + d*x))/a))/(a^3*d^2)]
    @test_int [x^0/(a + b*ℯ^(c + d*x))^3, x, 3, 1/(2*a*d*(a + b*ℯ^(c + d*x))^2) + 1/(a^2*d*(a + b*ℯ^(c + d*x))) + x/a^3 - log(a + b*ℯ^(c + d*x))/(a^3*d)]
    @test_int [1/(x^1*(a + b*ℯ^(c + d*x))^3), x, 0, Unintegrable(1/((a + b*ℯ^(c + d*x))^3*x), x)]
    @test_int [1/(x^2*(a + b*ℯ^(c + d*x))^3), x, 0, Unintegrable(1/((a + b*ℯ^(c + d*x))^3*x^2), x)]

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

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

    #= ::Subsection::Closed:: =#
    #=p>0=#

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

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

    @test_int [(c + d*x)^3*(a + b*(F^(g*(e + f*x)))^n)^3, x, 14, (a^3*(c + d*x)^4)/(4*d) - (18*a^2*b*d^3*(F^(e*g + f*g*x))^n)/(f^4*g^4*n^4*log(F)^4) - (9*a*b^2*d^3*(F^(e*g + f*g*x))^(2*n))/(8*f^4*g^4*n^4*log(F)^4) - (2*b^3*d^3*(F^(e*g + f*g*x))^(3*n))/(27*f^4*g^4*n^4*log(F)^4) + (18*a^2*b*d^2*(F^(e*g + f*g*x))^n*(c + d*x))/(f^3*g^3*n^3*log(F)^3) + (9*a*b^2*d^2*(F^(e*g + f*g*x))^(2*n)*(c + d*x))/(4*f^3*g^3*n^3*log(F)^3) + (2*b^3*d^2*(F^(e*g + f*g*x))^(3*n)*(c + d*x))/(9*f^3*g^3*n^3*log(F)^3) - (9*a^2*b*d*(F^(e*g + f*g*x))^n*(c + d*x)^2)/(f^2*g^2*n^2*log(F)^2) - (9*a*b^2*d*(F^(e*g + f*g*x))^(2*n)*(c + d*x)^2)/(4*f^2*g^2*n^2*log(F)^2) - (b^3*d*(F^(e*g + f*g*x))^(3*n)*(c + d*x)^2)/(3*f^2*g^2*n^2*log(F)^2) + (3*a^2*b*(F^(e*g + f*g*x))^n*(c + d*x)^3)/(f*g*n*log(F)) + (3*a*b^2*(F^(e*g + f*g*x))^(2*n)*(c + d*x)^3)/(2*f*g*n*log(F)) + (b^3*(F^(e*g + f*g*x))^(3*n)*(c + d*x)^3)/(3*f*g*n*log(F))]
    @test_int [(c + d*x)^2*(a + b*(F^(g*(e + f*x)))^n)^3, x, 11, (a^3*(c + d*x)^3)/(3*d) + (6*a^2*b*d^2*(F^(e*g + f*g*x))^n)/(f^3*g^3*n^3*log(F)^3) + (3*a*b^2*d^2*(F^(e*g + f*g*x))^(2*n))/(4*f^3*g^3*n^3*log(F)^3) + (2*b^3*d^2*(F^(e*g + f*g*x))^(3*n))/(27*f^3*g^3*n^3*log(F)^3) - (6*a^2*b*d*(F^(e*g + f*g*x))^n*(c + d*x))/(f^2*g^2*n^2*log(F)^2) - (3*a*b^2*d*(F^(e*g + f*g*x))^(2*n)*(c + d*x))/(2*f^2*g^2*n^2*log(F)^2) - (2*b^3*d*(F^(e*g + f*g*x))^(3*n)*(c + d*x))/(9*f^2*g^2*n^2*log(F)^2) + (3*a^2*b*(F^(e*g + f*g*x))^n*(c + d*x)^2)/(f*g*n*log(F)) + (3*a*b^2*(F^(e*g + f*g*x))^(2*n)*(c + d*x)^2)/(2*f*g*n*log(F)) + (b^3*(F^(e*g + f*g*x))^(3*n)*(c + d*x)^2)/(3*f*g*n*log(F))]
    @test_int [(c + d*x)^1*(a + b*(F^(g*(e + f*x)))^n)^3, x, 8, (a^3*(c + d*x)^2)/(2*d) - (3*a^2*b*d*(F^(e*g + f*g*x))^n)/(f^2*g^2*n^2*log(F)^2) - (3*a*b^2*d*(F^(e*g + f*g*x))^(2*n))/(4*f^2*g^2*n^2*log(F)^2) - (b^3*d*(F^(e*g + f*g*x))^(3*n))/(9*f^2*g^2*n^2*log(F)^2) + (3*a^2*b*(F^(e*g + f*g*x))^n*(c + d*x))/(f*g*n*log(F)) + (3*a*b^2*(F^(e*g + f*g*x))^(2*n)*(c + d*x))/(2*f*g*n*log(F)) + (b^3*(F^(e*g + f*g*x))^(3*n)*(c + d*x))/(3*f*g*n*log(F))]
    @test_int [(c + d*x)^0*(a + b*(F^(g*(e + f*x)))^n)^3, x, 4, a^3*x + (3*a^2*b*(F^(g*(e + f*x)))^n)/(f*g*n*log(F)) + (3*a*b^2*(F^(g*(e + f*x)))^(2*n))/(2*f*g*n*log(F)) + (b^3*(F^(g*(e + f*x)))^(3*n))/(3*f*g*n*log(F))]
    @test_int [(a + b*(F^(g*(e + f*x)))^n)^3/(c + d*x)^1, x, 8, (3*a^2*b*F^((e - (c*f)/d)*g*n - g*n*(e + f*x))*(F^(e*g + f*g*x))^n*ExpIntegralEi((f*g*n*(c + d*x)*log(F))/d))/d + (3*a*b^2*F^(2*(e - (c*f)/d)*g*n - 2*g*n*(e + f*x))*(F^(e*g + f*g*x))^(2*n)*ExpIntegralEi((2*f*g*n*(c + d*x)*log(F))/d))/d + (b^3*F^(3*(e - (c*f)/d)*g*n - 3*g*n*(e + f*x))*(F^(e*g + f*g*x))^(3*n)*ExpIntegralEi((3*f*g*n*(c + d*x)*log(F))/d))/d + (a^3*log(c + d*x))/d]
    @test_int [(a + b*(F^(g*(e + f*x)))^n)^3/(c + d*x)^2, x, 11, -(a^3/(d*(c + d*x))) - (3*a^2*b*(F^(e*g + f*g*x))^n)/(d*(c + d*x)) - (3*a*b^2*(F^(e*g + f*g*x))^(2*n))/(d*(c + d*x)) - (b^3*(F^(e*g + f*g*x))^(3*n))/(d*(c + d*x)) + (3*a^2*b*f*F^((e - (c*f)/d)*g*n - g*n*(e + f*x))*(F^(e*g + f*g*x))^n*g*n*ExpIntegralEi((f*g*n*(c + d*x)*log(F))/d)*log(F))/d^2 + (6*a*b^2*f*F^(2*(e - (c*f)/d)*g*n - 2*g*n*(e + f*x))*(F^(e*g + f*g*x))^(2*n)*g*n*ExpIntegralEi((2*f*g*n*(c + d*x)*log(F))/d)*log(F))/d^2 + (3*b^3*f*F^(3*(e - (c*f)/d)*g*n - 3*g*n*(e + f*x))*(F^(e*g + f*g*x))^(3*n)*g*n*ExpIntegralEi((3*f*g*n*(c + d*x)*log(F))/d)*log(F))/d^2]
    @test_int [(a + b*(F^(g*(e + f*x)))^n)^3/(c + d*x)^3, x, 14, -(a^3/(2*d*(c + d*x)^2)) - (3*a^2*b*(F^(e*g + f*g*x))^n)/(2*d*(c + d*x)^2) - (3*a*b^2*(F^(e*g + f*g*x))^(2*n))/(2*d*(c + d*x)^2) - (b^3*(F^(e*g + f*g*x))^(3*n))/(2*d*(c + d*x)^2) - (3*a^2*b*f*(F^(e*g + f*g*x))^n*g*n*log(F))/(2*d^2*(c + d*x)) - (3*a*b^2*f*(F^(e*g + f*g*x))^(2*n)*g*n*log(F))/(d^2*(c + d*x)) - (3*b^3*f*(F^(e*g + f*g*x))^(3*n)*g*n*log(F))/(2*d^2*(c + d*x)) + (3*a^2*b*f^2*F^((e - (c*f)/d)*g*n - g*n*(e + f*x))*(F^(e*g + f*g*x))^n*g^2*n^2*ExpIntegralEi((f*g*n*(c + d*x)*log(F))/d)*log(F)^2)/(2*d^3) + (6*a*b^2*f^2*F^(2*(e - (c*f)/d)*g*n - 2*g*n*(e + f*x))*(F^(e*g + f*g*x))^(2*n)*g^2*n^2*ExpIntegralEi((2*f*g*n*(c + d*x)*log(F))/d)*log(F)^2)/d^3 + (9*b^3*f^2*F^(3*(e - (c*f)/d)*g*n - 3*g*n*(e + f*x))*(F^(e*g + f*g*x))^(3*n)*g^2*n^2*ExpIntegralEi((3*f*g*n*(c + d*x)*log(F))/d)*log(F)^2)/(2*d^3)]

    #= ::Subsection::Closed:: =#
    #=p<0=#

    @test_int [(c + d*x)^3/(a + b*(F^(g*(e + f*x)))^n), x, 6, (c + d*x)^4/(4*a*d) - ((c + d*x)^3*log(1 + (b*(F^(g*(e + f*x)))^n)/a))/(a*f*g*n*log(F)) - (3*d*(c + d*x)^2*PolyLog(2, -((b*(F^(g*(e + f*x)))^n)/a)))/(a*f^2*g^2*n^2*log(F)^2) + (6*d^2*(c + d*x)*PolyLog(3, -((b*(F^(g*(e + f*x)))^n)/a)))/(a*f^3*g^3*n^3*log(F)^3) - (6*d^3*PolyLog(4, -((b*(F^(g*(e + f*x)))^n)/a)))/(a*f^4*g^4*n^4*log(F)^4)]
    @test_int [(c + d*x)^2/(a + b*(F^(g*(e + f*x)))^n), x, 5, (c + d*x)^3/(3*a*d) - ((c + d*x)^2*log(1 + (b*(F^(g*(e + f*x)))^n)/a))/(a*f*g*n*log(F)) - (2*d*(c + d*x)*PolyLog(2, -((b*(F^(g*(e + f*x)))^n)/a)))/(a*f^2*g^2*n^2*log(F)^2) + (2*d^2*PolyLog(3, -((b*(F^(g*(e + f*x)))^n)/a)))/(a*f^3*g^3*n^3*log(F)^3)]
    @test_int [(c + d*x)^1/(a + b*(F^(g*(e + f*x)))^n), x, 4, (c + d*x)^2/(2*a*d) - ((c + d*x)*log(1 + (b*(F^(g*(e + f*x)))^n)/a))/(a*f*g*n*log(F)) - (d*PolyLog(2, -((b*(F^(g*(e + f*x)))^n)/a)))/(a*f^2*g^2*n^2*log(F)^2)]
    @test_int [(c + d*x)^0/(a + b*(F^(g*(e + f*x)))^n), x, 5, x/a - log(a + b*(F^(g*(e + f*x)))^n)/(a*f*g*n*log(F))]
    @test_int [1/((c + d*x)^1*(a + b*(F^(g*(e + f*x)))^n)), x, 1, Unintegrable(1/((a + b*(F^(e*g + f*g*x))^n)*(c + d*x)), x)]
    @test_int [1/((c + d*x)^2*(a + b*(F^(g*(e + f*x)))^n)), x, 1, Unintegrable(1/((a + b*(F^(e*g + f*g*x))^n)*(c + d*x)^2), x)]

    @test_int [(c + d*x)^3/(a + b*(F^(g*(e + f*x)))^n)^2, x, 13, (c + d*x)^4/(4*a^2*d) - (c + d*x)^3/(a^2*f*g*n*log(F)) + (c + d*x)^3/(a*f*(a + b*(F^(g*(e + f*x)))^n)*g*n*log(F)) + (3*d*(c + d*x)^2*log(1 + (b*(F^(g*(e + f*x)))^n)/a))/(a^2*f^2*g^2*n^2*log(F)^2) - ((c + d*x)^3*log(1 + (b*(F^(g*(e + f*x)))^n)/a))/(a^2*f*g*n*log(F)) + (6*d^2*(c + d*x)*PolyLog(2, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^2*f^3*g^3*n^3*log(F)^3) - (3*d*(c + d*x)^2*PolyLog(2, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^2*f^2*g^2*n^2*log(F)^2) - (6*d^3*PolyLog(3, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^2*f^4*g^4*n^4*log(F)^4) + (6*d^2*(c + d*x)*PolyLog(3, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^2*f^3*g^3*n^3*log(F)^3) - (6*d^3*PolyLog(4, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^2*f^4*g^4*n^4*log(F)^4)]
    @test_int [(c + d*x)^2/(a + b*(F^(g*(e + f*x)))^n)^2, x, 11, (c + d*x)^3/(3*a^2*d) - (c + d*x)^2/(a^2*f*g*n*log(F)) + (c + d*x)^2/(a*f*(a + b*(F^(g*(e + f*x)))^n)*g*n*log(F)) + (2*d*(c + d*x)*log(1 + (b*(F^(g*(e + f*x)))^n)/a))/(a^2*f^2*g^2*n^2*log(F)^2) - ((c + d*x)^2*log(1 + (b*(F^(g*(e + f*x)))^n)/a))/(a^2*f*g*n*log(F)) + (2*d^2*PolyLog(2, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^2*f^3*g^3*n^3*log(F)^3) - (2*d*(c + d*x)*PolyLog(2, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^2*f^2*g^2*n^2*log(F)^2) + (2*d^2*PolyLog(3, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^2*f^3*g^3*n^3*log(F)^3)]
    @test_int [(c + d*x)^1/(a + b*(F^(g*(e + f*x)))^n)^2, x, 11, (c + d*x)^2/(2*a^2*d) - (d*x)/(a^2*f*g*n*log(F)) + (c + d*x)/(a*f*(a + b*(F^(g*(e + f*x)))^n)*g*n*log(F)) + (d*log(a + b*(F^(g*(e + f*x)))^n))/(a^2*f^2*g^2*n^2*log(F)^2) - ((c + d*x)*log(1 + (b*(F^(g*(e + f*x)))^n)/a))/(a^2*f*g*n*log(F)) - (d*PolyLog(2, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^2*f^2*g^2*n^2*log(F)^2)]
    @test_int [(c + d*x)^0/(a + b*(F^(g*(e + f*x)))^n)^2, x, 4, x/a^2 + 1/(a*f*(a + b*(F^(g*(e + f*x)))^n)*g*n*log(F)) - log(a + b*(F^(g*(e + f*x)))^n)/(a^2*f*g*n*log(F))]
    @test_int [1/((c + d*x)^1*(a + b*(F^(g*(e + f*x)))^n)^2), x, 1, Unintegrable(1/((a + b*(F^(e*g + f*g*x))^n)^2*(c + d*x)), x)]
    @test_int [1/((c + d*x)^2*(a + b*(F^(g*(e + f*x)))^n)^2), x, 1, Unintegrable(1/((a + b*(F^(e*g + f*g*x))^n)^2*(c + d*x)^2), x)]

    @test_int [(c + d*x)^3/(a + b*(F^(g*(e + f*x)))^n)^3, x, 26, (c + d*x)^4/(4*a^3*d) + (3*d*(c + d*x)^2)/(2*a^3*f^2*g^2*n^2*log(F)^2) - (3*d*(c + d*x)^2)/(2*a^2*f^2*(a + b*(F^(g*(e + f*x)))^n)*g^2*n^2*log(F)^2) - (3*(c + d*x)^3)/(2*a^3*f*g*n*log(F)) + (c + d*x)^3/(2*a*f*(a + b*(F^(g*(e + f*x)))^n)^2*g*n*log(F)) + (c + d*x)^3/(a^2*f*(a + b*(F^(g*(e + f*x)))^n)*g*n*log(F)) - (3*d^2*(c + d*x)*log(1 + (b*(F^(g*(e + f*x)))^n)/a))/(a^3*f^3*g^3*n^3*log(F)^3) + (9*d*(c + d*x)^2*log(1 + (b*(F^(g*(e + f*x)))^n)/a))/(2*a^3*f^2*g^2*n^2*log(F)^2) - ((c + d*x)^3*log(1 + (b*(F^(g*(e + f*x)))^n)/a))/(a^3*f*g*n*log(F)) - (3*d^3*PolyLog(2, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^3*f^4*g^4*n^4*log(F)^4) + (9*d^2*(c + d*x)*PolyLog(2, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^3*f^3*g^3*n^3*log(F)^3) - (3*d*(c + d*x)^2*PolyLog(2, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^3*f^2*g^2*n^2*log(F)^2) - (9*d^3*PolyLog(3, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^3*f^4*g^4*n^4*log(F)^4) + (6*d^2*(c + d*x)*PolyLog(3, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^3*f^3*g^3*n^3*log(F)^3) - (6*d^3*PolyLog(4, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^3*f^4*g^4*n^4*log(F)^4)]
    @test_int [(c + d*x)^2/(a + b*(F^(g*(e + f*x)))^n)^3, x, 24, (c + d*x)^3/(3*a^3*d) + (d^2*x)/(a^3*f^2*g^2*n^2*log(F)^2) - (d*(c + d*x))/(a^2*f^2*(a + b*(F^(g*(e + f*x)))^n)*g^2*n^2*log(F)^2) - (3*(c + d*x)^2)/(2*a^3*f*g*n*log(F)) + (c + d*x)^2/(2*a*f*(a + b*(F^(g*(e + f*x)))^n)^2*g*n*log(F)) + (c + d*x)^2/(a^2*f*(a + b*(F^(g*(e + f*x)))^n)*g*n*log(F)) - (d^2*log(a + b*(F^(g*(e + f*x)))^n))/(a^3*f^3*g^3*n^3*log(F)^3) + (3*d*(c + d*x)*log(1 + (b*(F^(g*(e + f*x)))^n)/a))/(a^3*f^2*g^2*n^2*log(F)^2) - ((c + d*x)^2*log(1 + (b*(F^(g*(e + f*x)))^n)/a))/(a^3*f*g*n*log(F)) + (3*d^2*PolyLog(2, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^3*f^3*g^3*n^3*log(F)^3) - (2*d*(c + d*x)*PolyLog(2, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^3*f^2*g^2*n^2*log(F)^2) + (2*d^2*PolyLog(3, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^3*f^3*g^3*n^3*log(F)^3)]
    @test_int [(c + d*x)^1/(a + b*(F^(g*(e + f*x)))^n)^3, x, 17, (c + d*x)^2/(2*a^3*d) - d/(2*a^2*f^2*(a + b*(F^(g*(e + f*x)))^n)*g^2*n^2*log(F)^2) - (3*d*x)/(2*a^3*f*g*n*log(F)) + (c + d*x)/(2*a*f*(a + b*(F^(g*(e + f*x)))^n)^2*g*n*log(F)) + (c + d*x)/(a^2*f*(a + b*(F^(g*(e + f*x)))^n)*g*n*log(F)) + (3*d*log(a + b*(F^(g*(e + f*x)))^n))/(2*a^3*f^2*g^2*n^2*log(F)^2) - ((c + d*x)*log(1 + (b*(F^(g*(e + f*x)))^n)/a))/(a^3*f*g*n*log(F)) - (d*PolyLog(2, -((b*(F^(g*(e + f*x)))^n)/a)))/(a^3*f^2*g^2*n^2*log(F)^2)]
    @test_int [(c + d*x)^0/(a + b*(F^(g*(e + f*x)))^n)^3, x, 4, x/a^3 + 1/(2*a*f*(a + b*(F^(g*(e + f*x)))^n)^2*g*n*log(F)) + 1/(a^2*f*(a + b*(F^(g*(e + f*x)))^n)*g*n*log(F)) - log(a + b*(F^(g*(e + f*x)))^n)/(a^3*f*g*n*log(F))]
    @test_int [1/((c + d*x)^1*(a + b*(F^(g*(e + f*x)))^n)^3), x, 1, Unintegrable(1/((a + b*(F^(e*g + f*g*x))^n)^3*(c + d*x)), x)]
    @test_int [1/((c + d*x)^2*(a + b*(F^(g*(e + f*x)))^n)^3), x, 1, Unintegrable(1/((a + b*(F^(e*g + f*g*x))^n)^3*(c + d*x)^2), x)]

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

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

    @test_int [sqrt(c + d*x)*(a + b*ℯ^x), x, 5, b*ℯ^x*sqrt(c + d*x) + (2*a*(c + d*x)^(3/2))/(3*d) - ((1/2)*b*sqrt(d)*sqrt(pi)*erfi(sqrt(c + d*x)/sqrt(d)))/ℯ^(c/d)]

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

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

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

    @test_int [sqrt(c + d*x)/(a + b*ℯ^x), x, 0, Unintegrable(sqrt(c + d*x)/(a + b*ℯ^x), x)]

    @test_int [sqrt(c + d*x)/(a + b*ℯ^x)^2, x, 0, Unintegrable(sqrt(c + d*x)/(a + b*ℯ^x)^2, x)]

    @test_int [sqrt(c + d*x)/(a + b*ℯ^x)^3, x, 0, Unintegrable(sqrt(c + d*x)/(a + b*ℯ^x)^3, x)]

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

    @test_int [(c + d*x)^m*(a + b*(F^(g*(e + f*x)))^n)^3, x, 8, (a^3*(c + d*x)^(1 + m))/(d*(1 + m)) + (3^(-1 - m)*b^3*F^(3*(e - (c*f)/d)*g*n - 3*g*n*(e + f*x))*(F^(e*g + f*g*x))^(3*n)*(c + d*x)^m*Gamma(1 + m, -((3*f*g*n*(c + d*x)*log(F))/d)))/((-((f*g*n*(c + d*x)*log(F))/d))^m*(f*g*n*log(F))) + (3*2^(-1 - m)*a*b^2*F^(2*(e - (c*f)/d)*g*n - 2*g*n*(e + f*x))*(F^(e*g + f*g*x))^(2*n)*(c + d*x)^m*Gamma(1 + m, -((2*f*g*n*(c + d*x)*log(F))/d)))/((-((f*g*n*(c + d*x)*log(F))/d))^m*(f*g*n*log(F))) + (3*a^2*b*F^((e - (c*f)/d)*g*n - g*n*(e + f*x))*(F^(e*g + f*g*x))^n*(c + d*x)^m*Gamma(1 + m, -((f*g*n*(c + d*x)*log(F))/d)))/((-((f*g*n*(c + d*x)*log(F))/d))^m*(f*g*n*log(F)))]
    @test_int [(c + d*x)^m*(a + b*(F^(g*(e + f*x)))^n)^2, x, 6, (a^2*(c + d*x)^(1 + m))/(d*(1 + m)) + (2^(-1 - m)*b^2*F^(2*(e - (c*f)/d)*g*n - 2*g*n*(e + f*x))*(F^(e*g + f*g*x))^(2*n)*(c + d*x)^m*Gamma(1 + m, -((2*f*g*n*(c + d*x)*log(F))/d)))/((-((f*g*n*(c + d*x)*log(F))/d))^m*(f*g*n*log(F))) + (2*a*b*F^((e - (c*f)/d)*g*n - g*n*(e + f*x))*(F^(e*g + f*g*x))^n*(c + d*x)^m*Gamma(1 + m, -((f*g*n*(c + d*x)*log(F))/d)))/((-((f*g*n*(c + d*x)*log(F))/d))^m*(f*g*n*log(F)))]
    @test_int [(c + d*x)^m*(a + b*(F^(g*(e + f*x)))^n)^1, x, 4, (a*(c + d*x)^(1 + m))/(d*(1 + m)) + (b*F^((e - (c*f)/d)*g*n - g*n*(e + f*x))*(F^(e*g + f*g*x))^n*(c + d*x)^m*Gamma(1 + m, -((f*g*n*(c + d*x)*log(F))/d)))/((-((f*g*n*(c + d*x)*log(F))/d))^m*(f*g*n*log(F)))]
    @test_int [(c + d*x)^m/(a + b*(F^(g*(e + f*x)))^n)^1, x, 1, Unintegrable((c + d*x)^m/(a + b*(F^(e*g + f*g*x))^n), x)]
    @test_int [(c + d*x)^m/(a + b*(F^(g*(e + f*x)))^n)^2, x, 1, Unintegrable((c + d*x)^m/(a + b*(F^(e*g + f*g*x))^n)^2, x)]

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

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

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

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

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

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

    @test_int [x^3*F^(c + d*x)/(a + b*F^(c + d*x))^3, x, 12, -((3*x^2)/(2*a^2*b*d^2*log(F)^2)) + (3*x^2)/(2*a*b*d^2*(a + b*F^(c + d*x))*log(F)^2) + x^3/(2*a^2*b*d*log(F)) - x^3/(2*b*d*(a + b*F^(c + d*x))^2*log(F)) + (3*x*log(1 + (b*F^(c + d*x))/a))/(a^2*b*d^3*log(F)^3) - (3*x^2*log(1 + (b*F^(c + d*x))/a))/(2*a^2*b*d^2*log(F)^2) + (3*PolyLog(2, -((b*F^(c + d*x))/a)))/(a^2*b*d^4*log(F)^4) - (3*x*PolyLog(2, -((b*F^(c + d*x))/a)))/(a^2*b*d^3*log(F)^3) + (3*PolyLog(3, -((b*F^(c + d*x))/a)))/(a^2*b*d^4*log(F)^4)]
    @test_int [x^2*F^(c + d*x)/(a + b*F^(c + d*x))^3, x, 11, -(x/(a^2*b*d^2*log(F)^2)) + x/(a*b*d^2*(a + b*F^(c + d*x))*log(F)^2) + x^2/(2*a^2*b*d*log(F)) - x^2/(2*b*d*(a + b*F^(c + d*x))^2*log(F)) + log(a + b*F^(c + d*x))/(a^2*b*d^3*log(F)^3) - (x*log(1 + (b*F^(c + d*x))/a))/(a^2*b*d^2*log(F)^2) - PolyLog(2, -((b*F^(c + d*x))/a))/(a^2*b*d^3*log(F)^3)]
    @test_int [x^1*F^(c + d*x)/(a + b*F^(c + d*x))^3, x, 4, 1/(2*a*b*d^2*(a + b*F^(c + d*x))*log(F)^2) + x/(2*a^2*b*d*log(F)) - x/(2*b*d*(a + b*F^(c + d*x))^2*log(F)) - log(a + b*F^(c + d*x))/(2*a^2*b*d^2*log(F)^2)]
    @test_int [x^0*F^(c + d*x)/(a + b*F^(c + d*x))^3, x, 2, -(1/(2*b*d*(a + b*F^(c + d*x))^2*log(F)))]
    @test_int [1/x^1*F^(c + d*x)/(a + b*F^(c + d*x))^3, x, 1, -(1/(2*b*d*(a + b*F^(c + d*x))^2*x*log(F))) - Unintegrable(1/((a + b*F^(c + d*x))^2*x^2), x)/(2*b*d*log(F))]
    @test_int [1/x^2*F^(c + d*x)/(a + b*F^(c + d*x))^3, x, 1, -(1/(2*b*d*(a + b*F^(c + d*x))^2*x^2*log(F))) - Unintegrable(1/((a + b*F^(c + d*x))^2*x^3), x)/(b*d*log(F))]
    end

max_stars_repo_path: test/2 Exponentials/2.2 (c+d x)^m (F^(g (e+f x)))^n (a+b (F^(g (e+f x)))^n)^p.jl


text:
using PkgBenchmark: benchmarkpkg, writeresults, export_markdown

mkpath("tmp")
results = benchmarkpkg("LyapunovExponents")
writeresults("tmp/benchmarks_current.json", results)
export_markdown("tmp/benchmarks_current.md", results)
showall(results)
println()


text:
module ProgrammingProjects

using HTTP
using Gumbo
using Cascadia

include("webcrawler.jl")
# using .webcrawler: PAGE_URL
# using .webcrawler: LINKS
using .webcrawler: fetchpage
using .webcrawler: extractlinks
using .webcrawler: extracttitle
using .webcrawler: extractimage
using .webcrawler: articleinfo
using .webcrawler: articledom

# export PAGE_URL
# export LINKS
export fetchpage
export extractlinks
export extracttitle
export extractimage
export articleinfo
export articledom

end


text:
export GitReflog, GitReflogEntry,
       set_target, set_symbolic_target, resolve,
       rename, target, symbolic_target, name,
       is_valid_ref, has_reflog, peel, log!

function GitReference(ptr::Ptr{Void})
    ty = ccall((:git_reference_type, libgit2), Cint, (Ptr{Void},), ptr)
    RType = ty == 1 ? Oid : Sym
    ref = GitReference{RType}(ptr)
    finalizer(ref, free!)
    return ref
end

Base.show{T}(io::IO, ref::GitReference{T}) = begin
    write(io, "GitReference{$T}(\"$(name(ref))\")")
end

is_valid_ref(ref::String) = bool(ccall((:git_reference_is_valid_name, libgit2), Cint, (Ptr{Uint8},), ref))

function set_symbolic_target(r::GitReference, target::String;
                             sig::MaybeSignature=nothing,
                             logmsg::MaybeString=nothing)
    ref_ptr = Ptr{Void}[0]
    @check ccall((:git_reference_set_symbolic_target, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Uint8}, Ptr{SignatureStruct}, Ptr{Uint8}),
                 ref_ptr, r, target,
                 sig != nothing ? sig : C_NULL,
                 msg != nothing ? msg : C_NULL)
    return GitReference(ref_ptr[1])
end

function set_target(r::GitReference, id::Oid;
                    sig::MaybeSignature=nothing,
                    logmsg::MaybeString=nothing)
    ref_ptr = Ptr{Void}[0]
    @check ccall((:git_reference_set_target, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Oid}, Ptr{SignatureStruct}, Ptr{Uint8}),
                 ref_ptr, r, &id,
                 sig != nothing ? sig : C_NULL,
                 logmsg != nothing ? logmsg : C_NULL)
    return GitReference(ref_ptr[1])
end

function resolve(r::GitReference)
    ref_ptr = Ptr{Void}[0]
    @check ccall((:git_reference_resolve, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}), ref_ptr, r)
    return GitReference(ref_ptr[1])
end

function rename(r::GitReference, name::String;
                force::Bool=false,
                sig::MaybeSignature=nothing,
                logmsg::MaybeString=nothing)
    ref_ptr = Ptr{Void}[0]
    @check ccall((:git_reference_rename, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Uint8}, Cint, Ptr{SignatureStruct}, Ptr{Uint8}),
                 ref_ptr, r, name, force,
                 sig != nothing ? sig : C_NULL,
                 logmsg != nothing ? logmsg : C_NULL)
    return GitReference(ref_ptr[1])
end

function target(r::GitReference)
    oid_ptr = ccall((:git_reference_target, libgit2), Ptr{Oid}, (Ptr{Void},), r)
    if oid_ptr == C_NULL
        return nothing
    end
    return Oid(oid_ptr)
end

function symbolic_target(r::GitReference)
    sym_ptr = ccall((:git_reference_symbolic_target, libgit2), Ptr{Uint8}, (Ptr{Void},), r)
    if sym_ptr == C_NULL
        return utf8("")
    end
    return utf8(bytestring(sym_ptr))
end

name(r::GitReference) = utf8(bytestring(ccall((:git_reference_name, libgit2), Ptr{Uint8}, (Ptr{Void},), r)))

function peel{T}(r::GitReference{T})
    obj_ptr = Ptr{Void}[0]
    err = ccall((:git_reference_peel, libgit2), Cint,
                (Ptr{Ptr{Void}}, Ptr{Void}, Cint), obj_ptr, r, GitConst.OBJ_ANY)
    if err == GitErrorConst.ENOTFOUND
        return nothing
    elseif err != GitErrorConst.GIT_OK
        if obj_ptr[1] != C_NULL
            ccall((:git_object_free, libgit2), Void, (Ptr{Void},), obj_ptr[1])
        end
        throw(LibGitError(err))
    end
    if is(T, Oid)
        oid = Oid(ccall((:git_object_id, libgit2), Ptr{Oid}, (Ptr{Void},), obj_ptr[1]))
        tid = Oid(ccall((:git_reference_target, libgit2), Ptr{Oid}, (Ptr{Void},), r))
        if oid == tid
            ccall((:git_object_free, libgit2), Void, (Ptr{Void},), obj_ptr[1])
            return nothing
        end
    end
    id = Oid(ccall((:git_object_id, libgit2), Ptr{Oid}, (Ptr{Void},), obj_ptr[1]))
    ccall((:git_object_free, libgit2), Void, (Ptr{Void},), obj_ptr[1])
    return id
end

type GitReflog
    ptr::Ptr{Void}

    function GitReflog(ptr::Ptr{Void})
        @assert ptr != C_NULL
        rl = new(ptr)
        finalizer(fl, free!)
        return rl
    end
end

free!(r::GitReflog) = begin
    if r.ptr != C_NULL
        ccall((:git_reflog_free, libgit2), Void, (Ptr{Void},), r.ptr)
        r.ptr = C_NULL
    end
end

Base.convert(::Type{Ptr{Void}}, r::GitReflog) = r.ptr

type GitReflogEntry
    id_old::Oid
    id_new::Oid
    committer::Signature
    message::UTF8String
end

function new_reflog_entry(entry_ptr::Ptr{Void})
    id_old  = Oid(ccall((:git_reflog_entry_id_old, libgit2), Ptr{Oid}, (Ptr{Void},), entry_ptr))
    id_new  = Oid(ccall((:git_reflog_entry_id_new, libgit2), Ptr{Oid}, (Ptr{Void},), entry_ptr))
    sig_ptr = ccall((:git_reflog_entry_committer, libgit2), Ptr{SignatureStruct}, (Ptr{Void},), entry_ptr)
    msg_ptr = ccall((:git_reflog_entry_message, libgit2), Ptr{Uint8}, (Ptr{Void},), entry_ptr)
    msg = msg_ptr == C_NULL ? utf8("") : utf8(bytestring(msg_ptr))
    sig = Signature(sig_ptr)
    return GitReflogEntry(id_old, id_new, sig, msg)
end

GitReflog(r::GitReference) = begin
    reflog_ptr = Ptr{Void}[0]
    owner_ptr = ccall((:git_reference_owner, libgit2), Ptr{Void}, (Ptr{Void},), r)
    @check ccall((:git_reflog_read, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Uint8}), reflog_ptr, owner_ptr, name(r))
    refcount = ccall((:git_reflog_entrycount, libgit2), Csize_t, (Ptr{Void},), reflog_ptr[1])
    entries = GitReflogEntry[]
    for i in 0:refcount-1
        entry_ptr = ccall((:git_reflog_entry_byindex, libgit2), Ptr{Void},
                          (Ptr{Void}, Csize_t), reflog_ptr[1], refcount - i - 1)
        @assert entry_ptr != C_NULL
        push!(entries, new_reflog_entry(entry_ptr))
    end
    ccall((:git_reflog_free, libgit2), Void, (Ptr{Void},), reflog_ptr[1])
    return entries
end

has_reflog(r::GitReference) = bool(ccall((:git_reference_has_log, libgit2), Cint, (Ptr{Void},), r.ptr))

function log!(r::GitReference, msg::String="", committer::MaybeSignature=nothing)
    reflog_ptr = Ptr{Void}[0]
    repo_ptr = ccall((:git_reference_owner, libgit2), Ptr{Void}, (Ptr{Void},), r)
    name_ptr = ccall((:git_reference_name, libgit2), Ptr{Uint8}, (Ptr{Void},), r)
    @check ccall((:git_reflog_read, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Uint8}), reflog_ptr, repo_ptr, name_ptr)
    local sig_ptr::Ptr{SignatureStruct}
    if committer == nothing
        ptr = Ptr{SignatureStruct}[0]
        @check ccall((:git_signature_default, libgit2), Cint,
                     (Ptr{Ptr{SignatureStruct}}, Ptr{Void}), ptr, repo_ptr)
        sig_ptr = ptr[1]
    else
        sig_ptr = convert(Ptr{SignatureStruct}, committer)
    end
    target_ptr = ccall((:git_reference_target, libgit2), Ptr{Void}, (Ptr{Void},), r)
    err = ccall((:git_reflog_append, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{SignatureStruct}, Ptr{Uint8}),
                 reflog_ptr[1], target_ptr, sig_ptr, msg)
    if err == GitErrorConst.GIT_OK
        err = ccall((:git_reflog_write, libgit2), Cint, (Ptr{Void},), reflog_ptr[1])
    end
    @assert reflog_ptr[1] != C_NULL
    ccall((:git_reflog_free, libgit2), Void, (Ptr{Void},), reflog_ptr[1])
    if err != GitErrorConst.GIT_OK
        throw(LibGitError(err))
    end
    return nothing
end

git_reftype{T}(r::GitReference{T}) = begin
    T <: Sym && return GitConst.REF_SYMBOLIC
    T <: Oid && return GitConst.REF_OID
    throw(ArgumentError("Unknown reference type $T"))
end


text:
cd("INSERT_PATH")

using Distributions
using PolynomialRoots
using DataFrames

include("toolBox.jl")
include("quadDNRR.jl")

#########################################################
# Parameter Declaration
#########################################################

ss = 1000

# same \eta, different \sigma
# w1, lam1 = 0.6, 4.8e04 # f_{2,1},/ seed: 13
# w1, lam1 = 0.7, 3.5e04 # f_{2,2}, dseed: 116, sseed = 122, sp = 3e-06
# w1, lam1 = 0.8, 2.5e04 # f_{2,3}, dseed: 13, sseed = 1236, sp = 4.5e-06
w1, lam1 = 0.9, 1.3e04 # f_{2,4}, dseed: 13, sseed = 1214, sp = 4.5e-06

# different \eta, same \sigma
# w1, lam1 = 0.6, 2e04 # f_{2,1}, seed: 8152
# w1, lam1 = 0.7, 2e04 # f_{2,2}, seed: 1134
# w1, lam1 = 0.8, 2e04 # f_{2,3}, seed: 89
# w1, lam1 = 0.9, 2e04 # f_{2,4}, seed: 58

lam2 = 1e-02 * lam1 # noise mixture 3 {dynamical noise}
# w1, lam1, lam2 = 0.3, 2.5e03, 2.5e03 # noise mixture 3 {dynamical noise}

# # map parameters {for sampler}
# # theta, x0, x00 = [1.4, 0., 0.3, 0., -1., 0.], 0.5, 1.0
# theta = [1.31, 0., 0.23, 0., -1., 0.]
theta, x0, x00 = [1.38, 0., 0.27, 0., -1., 0.], 0.5, 0.5

# x0, x00 =  0.1, 0.1 # initial conditions

# map parameters {for data generation}
# thetamap = [1.4, -1., 0.3] # pars1 (standard)
# thetamap = [1., -1.31, 0.23] # pars2 (SJH)
# thetamap = [1.38, -1., 0.27] # pars3 (Kantz)

# corresponding deterministic orbit
xdet = zeros(ss)
xdet[1] = henon(theta, x0, x00)
xdet[2] = henon(theta, xdet[1], x0)
for i in 3:ss
  xdet[i] = henon(theta, xdet[i-1], xdet[i-2])
end

# orbit contaminated with dynamical noise (latent in State Space)
# first argument is  type of noise: 0 is parametric, 1-3 nonparametric
dataSeed = 13
data = genData(3, ss, theta, x0, x00, w1, lam1, lam2, dataSeed)
eta = sqrt(w1*1/lam1+(1-w1)*1/lam2)/std(data)

# read data set
path = "/home/kkaloudis/Documents/Mathematics/PhD/Julia/Prediction"
x1 = readdlm(string(path,"/lynx.txt"))
x=zeros(114)
for i in 1:114
	x[i] = x1[i]
end

#########################################################
# GSB Reconstruction - Denoising
#########################################################

# Dynamical Noise Only

# include("quadDPDeno.jl")
# path = "/home/kkaloudis/Documents/DNRR paper codes/full quadratic/Results"
# data = readdlm(string(path,"/data1.txt"))
# xdet = readdlm(string(path,"/xdet1.txt"))

k, gibbsIter, burnIn, dela, delb, wdist, gEps1, gEps2, thLow, thUp, zLow, zUp, papr, pbpr, thin, samplerSeed =
  0, 1000000, 1000000, 1e05, 1e-03, 1e02, 1e-03, 1e-03, -10.0, 10.0, -20.0, 20.0, 0.5, 0.5, 100, 76169;

filename = "/Results"
savelocation = string(pwd(), filename, "/seed$samplerSeed")
mkpath(savelocation)
writedlm(string(savelocation,"/xdet.txt"), xdet)

## DNRR application
@time thetas = quadDNRR(x, k, gibbsIter, burnIn, dela, delb, wdist, gEps1, gEps2, thLow,
  thUp, zLow, zUp, papr, pbpr, thin, samplerSeed);

thetas = 0;
gc()

# ################################################################
# ################################################################

max_stars_repo_name: KosKal/DNRR

max_issues_repo_name: KosKal/DNRR

max_forks_repo_name: KosKal/DNRR


text:
#Read/write EasyPlot plots to/from HDF5
#-------------------------------------------------------------------------------

#==Useful constants
===============================================================================#

#==Axis read/write functions
===============================================================================#
#(Axis is a bit strange to write out)

function _write_typed(grp::HDF5.Group, name::String, v::Axis{T}) where T
	grp[name] = String(T)
	_write_datatype_attr(grp[name], Axis)
end

function __read(::Type{Axis}, ds::HDF5.Dataset)
	scale = Symbol(HDF5.read(ds))
	return Axis(scale)
end

#==Plot read/write functions
===============================================================================#

function _writeplot(grp::HDF5.Group, plot::Plot)
	_write_typed(grp, "title", plot.title)
	_write_typed(grp, "xlabel", plot.xlabel)
	_write_typed(grp, "legend", plot.legend)
	_write_typed(grp, "xaxis", plot.xaxis)
	_write_typed(grp, "xext", plot.xext)

	grpslist = creategrp(grp, "ystriplist")
	_write_length_attr(grpslist, plot.ystriplist)
	for (i, strip) in enumerate(plot.ystriplist)
		_write_typed(grpslist, "strip$i", strip)
	end

	grpwlist = creategrp(grp, "wfrmlist")
	_write_length_attr(grpwlist, plot.wfrmlist)
	for (i, wfrm) in enumerate(plot.wfrmlist)
		_write_typed(grpwlist, "wfrm$i", wfrm)
	end
	return
end

function _readplot(grp::HDF5.Group)
	plot = Plot(title=_read_typed(grp, "title"))
	plot.ystriplist = []
	plot.xlabel = _read_typed(grp, "xlabel")
	plot.legend = _read_typed(grp, "legend")
	plot.xaxis = _read_typed(grp, "xaxis")
	plot.xext = _read_typed(grp, "xext")

	grpslist = opengrp(grp, "ystriplist")
	_len = _read_length_attr(Vector, grpslist)
	for i in 1:_len
		strip = _read_typed(grpslist, "strip$i")
		push!(plot.ystriplist, strip)
	end

	grpwlist = opengrp(grp, "wfrmlist")
	_len = _read_length_attr(Vector, grpwlist)
	for i in 1:_len
		wfrm = _read_typed(grpwlist, "wfrm$i")
		push!(plot, wfrm)
	end

	return plot
end

#==PlotCollection read/write functions
===============================================================================#
function _writepcoll(grp::HDF5.Group, pcoll::PlotCollection)
	_write_typed(grp, "title", pcoll.title)
	_write_typed(grp, "ncolumns", pcoll.ncolumns)

	grpplist = creategrp(grp, "plotlist")
	_write_length_attr(grpplist, pcoll.plotlist)
	for (i, plot) in enumerate(pcoll.plotlist)
		grpplot = creategrp(grpplist, "plot$i")
		_writeplot(grpplot, plot)
	end
	return
end

function _readpcoll(grp::HDF5.Group)
	pcoll = PlotCollection(title=_read_typed(grp, "title"))
	pcoll.ncolumns = _read_typed(grp, "ncolumns")

	grpplist = opengrp(grp, "plotlist")
	_len = _read_length_attr(Vector, grpplist)
	for i in 1:_len
		grpplot = opengrp(grpplist, "plot$i")
		plot = _readplot(grpplot)
		push!(pcoll, plot)
	end

	return pcoll
end

#==User-level interface:
===============================================================================#

function Base.write(w::EasyDataWriter, pcoll::PlotCollection, name::String)
	grp = creategrp(w, "$hdf5plotcollroot/$name")
	return _writepcoll(grp, pcoll)
end
function readplot(r::EasyDataReader, name::String) #Avoids specifying ::Type{PlotCollection}
	grp = opengrp(r, "$hdf5plotcollroot/$name")
	return _readpcoll(grp)
end
Base.read(::Type{PlotCollection}, r::EasyDataReader, name::String) = readplot(r, name)

#Interface to read/write by filepath:
#-------------------------------------------------------------------------------
function writeplot(filepath::String, pcoll::PlotCollection; name::String="_unnamed")
	openwriter(filepath) do w
		write(w, pcoll, name)
	end
end

function readplot(filepath::String; name::String="_unnamed")
	openreader(filepath) do r
		readplot(r, name)
	end
end

#Last Line


text:
# ---
# title: 19. Remove Nth Node From End of List
# id: problem19
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Linked List, Two Pointers
# link: <https://leetcode.com/problems/remove-nth-node-from-end-of-list/description/>
# hidden: true
# ---
#
# Given the `head` of a linked list, remove the `nth` node from the end of the
# list and return its head.
#
# **Follow up:**  Could you do this in one pass?
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/03/remove_ex1.jpg)
#
#
#
#     Input: head = [1,2,3,4,5], n = 2
#     Output: [1,2,3,5]
#
#
# **Example 2:**
#
#
#
#     Input: head = [1], n = 1
#     Output: []
#
#
# **Example 3:**
#
#
#
#     Input: head = [1,2], n = 1
#     Output: [1]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the list is `sz`.
#   * `1 <= sz <= 30`
#   * `0 <= Node.val <= 100`
#   * `1 <= n <= sz`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Provides FFT functionality, regridding,

module visibilities

using AbstractFFTs
using FFTW
using HDF5
using ..image
using ..gridding
using ..constants

import Base.conj! # extend this for DataVis
import Base.- # extend this for FullModelVis

export DataVis, ModelVis, RawModelVis, FullModelVis, fillModelVis, ResidVis, ModelVisRotate, DataVisReal, ModelVisReal
export get_qq
export plan_interpolate, interpolate_uv
export transform, rfftfreq, fftfreq, phase_shift, phase_shift!, max_baseline, get_nyquist_pixel
export lnprob, lnprob_iter, lnprob_real
export -

using Test

# Stores the data visibilities for a single channel
mutable struct DataVis
    lam::Float64 # [μm] Wavelength (in microns) corresponding to channel
    uu::Vector{Float64} # [kλ] Vectors of the u, v locations in kilolambda, shape (nchan, nvis)
    vv::Vector{Float64} # [kλ]
    VV::Vector{ComplexF64} # [Jy] Complex visibilities
    invsig::Vector{Float64} # 1/sigma [1/Jy]
    # NOTE that there is no flags field
end

"
    DataVis(fname::AbstractString, flagged::Bool=false)

Read all the visibilities from an HDF5 string and return them as an array of DataVis objects
`flagged` means whether we should be loading the flagged points or not."
function DataVis(fname::AbstractString, flagged::Bool=false)
    fid = h5open(fname, "r")

    freqs = read(fid["freqs"]) # [Hz]
    # Convert from Hz to wavelengths in μm
    lams = cc ./ freqs * 1e4 # [μm]

    uu = read(fid, "uu")
    vv = read(fid, "vv")
    real = read(fid, "real")
    imag = read(fid, "imag")
    VV = real + imag .* im # Complex visibility
    weight = read(fid, "weight") # [1/Jy^2]

    if !flagged
        flag = convert(Array{Bool}, read(fid, "flag"))
    end

    close(fid)

    nlam = length(lams)

    # Do we want to include the flagged visibility points? In nearly all cases, the answer
    # should be no. The one exception is when we want to export model visibilities in `write_model.jl`

    # Do the flagging channel by channel

    # Return an array of DataVis
    out = Array{DataVis}(undef, nlam)
    if !flagged
        for i=1:nlam
            ch_flag = flag[:,i]
            out[i] = DataVis(lams[i], uu[ch_flag,i], vv[ch_flag,i], VV[ch_flag, i], sqrt.(weight[ch_flag, i]))
        end
    else
        for i=1:nlam
            out[i] = DataVis(lams[i], uu[:,i], vv[:,i], VV[:, i], sqrt.(weight[:, i]))
        end
    end
    return out
end

"
    DataVis(fname::AbstractString, index::Int, flagged::Bool=false)

Read just one channel of visibilities from the HDF5 file."
function DataVis(fname::AbstractString, index::Int, flagged::Bool=false)
    fid = h5open(fname, "r")
    # the indexing and `vec` are necessary here because HDF5 doesn't naturally
    # squeeze trailing dimensions of length-1

    freq = read(fid["freqs"])[index][1] # [Hz]
    # Convert from Hz to wavelengths in μm
    lam = cc / freq * 1e4 # [μm]

    len = length(fid["uu"][:, index])

    if !flagged
        flag = convert(Array{Bool}, read(fid["flag"]))
        ch_flag = flag[:,index]
    else
        ch_flag = ones(Bool, len) # Keep all visibilities, regardless of what flag says
    end

    uu = vec(fid["uu"][:,index])[ch_flag]
    vv = vec(fid["vv"][:,index])[ch_flag]
    real = vec(fid["real"][:,index])[ch_flag]
    imag = vec(fid["imag"][:,index])[ch_flag]

    VV = real + imag .* im # Complex visibility

    invsig = vec(sqrt.(fid["weight"][:,index]))[ch_flag] # [1/Jy^2]

    close(fid)

    #Return a DataVis object
    return DataVis(lam, uu, vv, VV, invsig)
end

"
    DataVis(fname::AbstractString, indices::Vector{Int}, flagged::Bool=false)

Read just a subset of channels from the HDF5 file and return an array of DataVis."
function DataVis(fname::AbstractString, indices::Vector{Int}, flagged::Bool=false)
    nchan = length(indices)
    out = Array{DataVis}(undef, nchan)
    for i=1:nchan
        out[i] = DataVis(fname, indices[i], flagged)
    end
    return out
end

"
    DataVis(fname::AbstractString, mask::BitArray, flagged::Bool=false)

Read just a subset of channels from the HDF5 file and return an array of DataVis."
function DataVis(fname::AbstractString, mask::BitArray{1}, flagged::Bool=false)
    # load the array and mask it
    return DataVis(fname, flagged)[mask]
end

struct DataVisReal
    lam::Float64 # [μm] Wavelength (in microns) corresponding to channel
    uu::Vector{Float64} # [kλ] Vectors of the u, v locations in kilolambda, shape (nchan, nvis)
    vv::Vector{Float64} # [kλ]
    VVr::Vector{Float64} # [Jy] Complex visibilities
    VVi::Vector{Float64} # [Jy] Complex visibilities
    invsig::Vector{Float64} # 1/sigma [1/Jy]
    # NOTE that there is no flags field
end

"
    DataVisReal(dvis::DataVis)

Convert a `DataVis` object, with visibilities stored as `Comlpex128`, into an object with visibilities stored with real and complex values stored separately as `Float64`."
function DataVisReal(dvis::DataVis)
    return DataVisReal(dvis.lam, dvis.uu, dvis.vv, real(dvis.VV), imag(dvis.VV), dvis.invsig)
end

"
    conj!(dv::DataVis)

Import the complex conjugate function from Base, and extend it to work
on a DataVis.
I think this is necessary because the SMA and ALMA baseline conventions are swapped from
what I'm using in the NRAO Synthesis Summer School textbook."
function conj!(dv::DataVis)
    conj!(dv.VV)
end

"
    conj!(dvarr::Array{DataVis, 1})

Apply ``conj!`` to an entire DataVis array."
function conj!(dvarr::Array{DataVis, 1})
    for dset in dvarr
        conj!(dset) # Swap UV convention
    end
end

"
    get_qq(dv::DataVis)

Given a DataVis, convert the Cartesian spatial frequency coordinates ``(u,v)`` into a radial spatial coordinate ``q``."
function get_qq(dv::DataVis)
  sqrt.(dv.uu.^2 .+ dv.vv.^2)
end

"
    write(dv::DataVis, fname::AbstractString)

Take in a visibility data set and then write it to the HDF5 file
The HDF5 file actually expects multi-channel data, so instead we will need to
store all of this information with arrays of shape (..., 1) [an extra trailing]
dimension of 1."
function write(dv::DataVis, fname::AbstractString)
    nvis = length(dv.uu)
    VV = reshape(dv.VV, (nvis, 1))
    fid = h5open(fname, "w")

    fid["freqs"] = [cc/(1e-4 * dv.lam)] # expects 1D array [Hz]
    fid["uu"] = reshape(dv.uu, (nvis, 1)) # expects 2D array
    fid["vv"] = reshape(dv.vv, (nvis, 1)) # expects 2D array
    fid["real"] = real(VV)
    fid["imag"] = imag(VV)
    fid["weight"] = reshape(dv.invsig, (nvis, 1)).^2 # expects 2D array

    close(fid)
end

"
    write(dvarr::Array{DataVis, 1}, fname::AbstractString)

Write an array of DataVis to a single HDF5 file."
function write(dvarr::Array{DataVis, 1}, fname::AbstractString)
    nvis = length(dvarr[1].VV)
    nlam = length(dvarr)
    fid = h5open(fname, "w")
    # hcat here stacks the individual channel data sets into a big block
    # of shape (nvis, nlam)

    fid["freqs"] = Float64[cc / (1e-4 * dv.lam) for dv in dvarr] # [Hz]
    fid["uu"] = hcat([dv.uu for dv in dvarr]...)
    fid["vv"] = hcat([dv.vv for dv in dvarr]...)
    fid["real"] = hcat([real(dv.VV) for dv in dvarr]...)
    fid["imag"] = hcat([imag(dv.VV) for dv in dvarr]...)
    fid["weight"] = hcat([dv.invsig.^2 for dv in dvarr]...)
    close(fid)

end

"
    copy_flags(source::AbstractString, dest::AbstractString)

Copy the flags from one dataset to another."
function copy_flags(source::AbstractString, dest::AbstractString)
    println("Copying flags from $source to $dest")

    fid_source = h5open(source, "r") # read only
    fid_dest = h5open(dest, "r+") # append mode

    @assert isapprox(read(fid_source, "uu"), read(fid_dest, "uu")) "UU spacings between datasets do not match, make sure you have loaded $source with flagged=true."
    @assert isapprox(read(fid_source, "vv"), read(fid_dest, "vv")) "VV spacings between dataset do not match, make sure you have loaded $source with flagged=true."
    @assert isapprox(read(fid_source, "weight"), read(fid_dest, "weight")) "Weights between dataset do not match, make sure you have loaded $source with flagged=true."

    fid_dest["flag"] = read(fid_source, "flag")

    close(fid_source)
    close(fid_dest)
end

"
    max_baseline(dvarr)

Determine the maximum uu or vv baseline contained in the dataset, so we know at what resolution we will need to synthesize the images.

returned in kilolambda.
"
function max_baseline(dvarr::Array{DataVis, 1})
    max = 0.0

    for dv in dvarr
        max_uu = maximum(dv.uu)
        if max_uu > max
            max = max_uu
        end

        max_vv = maximum(dv.vv)
        if max_vv > max
            max = max_vv
        end
    end

    return max
end

"
    get_nyquist_pixel(max_base::Float64, angular_width::Float64)

Determine how many pixels the image needs at a given distance in order to satisfy the Nyquist sampling theorem.

`max_base` is in kilolambda.
`angular_width` is in radians."
function get_nyquist_pixel(max_base::Float64, angular_width::Float64)

    nyquist_factor = 4 # normally 2, but we want to be extra sure.
    #Calculate the maximum dRA and dDEC from max_base
    dRA_max = 1/(nyquist_factor * max_base * 1e3) # [radians]

    npix = 64
    dRA = angular_width/npix

    while dRA > dRA_max
        npix *= 2
        dRA = angular_width/npix
    end

    return npix
end

"
    RawModelVis

Type to store the product of FFT'ing a single channel of a SkyImage."
mutable struct RawModelVis
    lam::Float64 # [μm] Wavelength (in microns) corresponding to channel
    uu::Vector{Float64} # [kλ] Vectors of the u, v locations
    vv::Vector{Float64} # [kλ]
    VV::Matrix{ComplexF64} # Output from rfft
end

# Taking the complex conjugate to make a full grid over all visibility space.
mutable struct FullModelVis
    lam::Float64 # [μm] Wavelength (in microns) corresponding to channel
    uu::Vector{Float64} # [kλ] Vectors of the u, v locations
    vv::Vector{Float64} # [kλ]
    VV::Matrix{ComplexF64} # Output from rfft
end

"
    -(vis1::FullModelVis, vis2::FullModelVis)

Subtract two `FullModelVis`'s.
"
function -(vis1::FullModelVis, vis2::FullModelVis)
    # @assert vis1.lam == vis2.lam "Visibilities must have the same wavelengths."
    @assert vis1.uu == vis2.uu "Visibilities must have same uu sampling."
    @assert vis1.vv == vis2.vv "Visibilities must have same vv sampling."
    VV = vis1.VV - vis2.VV
    return FullModelVis(vis1.lam, vis1.uu, vis1.vv, VV)
end

# Produced by gridding a RawModelVis to match the data
mutable struct ModelVis
    dvis::DataVis # A reference to the matching dataset
    VV::Vector{ComplexF64} # vector of complex visibilities directly
    # corresponding to the u, v locations in DataVis
end

"
    ModelVis(dvis::DataVis, fmvis::FullModelVis)

Given a DataSet and a FullModelVis, go through and interpolate at the
u,v locations of the DataVis."
function ModelVis(dvis::DataVis, fmvis::FullModelVis)
    nvis = length(dvis.VV)
    VV = Array{ComplexF64}(undef, nvis)
    for i=1:nvis
        VV[i] = interpolate_uv(dvis.uu[i], dvis.vv[i], fmvis)
    end

    return ModelVis(dvis, VV)
end

"
    ModelVisRotate(dvis::DataVis, fmvis::FullModelVis, PA::Real)

Given a DataSet, FullModelVis, and a desired position angle rotation, use the Fourier rotation theorem to sample the image at the rotated baselines and produce a rotated (sampled) model. A positive `PA` value (in degrees) means that the new model is rotated ``PA`` number of degrees counter-clockwise in the image plane (towards East, from North)."
function ModelVisRotate(dvis::DataVis, fmvis::FullModelVis, PA::Real)

    # Convert from degrees to radians
    PA = PA * deg

    nvis = length(dvis.VV)

    VV = Array{ComplexF64}(undef, nvis)
    for i=1:nvis
        # Rotate points according to PA
        uuprime = cos(PA) * dvis.uu[i] + sin(PA) * dvis.vv[i]
        vvprime = -sin(PA) * dvis.uu[i] + cos(PA) * dvis.vv[i]

        # Interpolate at rotated points
        VV[i] = interpolate_uv(uuprime, vvprime, fmvis)
    end

    return ModelVis(dvis, VV)
end

"
    ModelVis2DataVis(mvis::ModelVis)

Collapse a ModelVis into a DataVis in order to write to disk."
function ModelVis2DataVis(mvis::ModelVis)
    DV = mvis.dvis
    return DataVis(DV.lam, DV.uu, DV.vv, mvis.VV, DV.invsig)
end

function conj!(mv::ModelVis)
    conj!(mv.VV)
end

"
    ResidVis(dvarr::Array{DataVis, 1}, mvarr)

Return a model visibility file that contains the residuals."
function ResidVis(dvarr::Array{DataVis, 1}, mvarr)
    nchan = length(dvarr)
    @assert length(mvarr) == nchan # make sure data and model are the same length.
    rvarr = Array{DataVis}(undef, nchan)
    for i=1:nchan
        dvis = dvarr[i]
        mvis = mvarr[i]
        VV = dvis.VV - mvis.VV
        rvarr[i] = DataVis(dvis.lam, dvis.uu, dvis.vv, VV, dvis.invsig)
    end
    return rvarr
end

mutable struct ModelVisReal
    dvis::DataVis # A reference to the matching dataset
    VVr::Vector{Float64} # vector of complex visibilities directly
    VVi::Vector{Float64} # vector of complex visibilities directly
    # corresponding to the u, v locations in DataVis
end

"
    ModelVisReal(mvis::ModelVis)

Convert a `ModelVis` object, with visibilities stored as `ComplexF64`, into a `ModelVisReal`, with real and imaginary visibility values stored separately as `Float64`."
function ModelVisReal(mvis::ModelVis)
    return ModelVisReal(mvis.dvis, real(mvis.VV), imag(mvis.VV))
end

"
    lnprob(dvis::DataVis, mvis::DataVis)

Lnprob for computing the difference between two datasets. Does not check whether the baselines of the two datasets are the same. Use this function at your own discretion."
function lnprob(dvis::DataVis, mvis::DataVis)
    return -0.5 * sum(abs2, dvis.invsig .* (dvis.VV - mvis.VV)) # Basic chi2
end

"
    lnprob(dvis::DataVis, mvis::ModelVis)

Likelihood function used to calculate the likelihood of the data visibilities given a current model of them. Asserts that the model visibilities were sampled at the same baselines as the data."
function lnprob(dvis::DataVis, mvis::ModelVis)
    @assert dvis === mvis.dvis # Using the wrong ModelVis, otherwise!
    return -0.5 * sum(abs2, dvis.invsig .* (dvis.VV - mvis.VV)) # Basic chi2
end

"
    lnprob(dvis::DataVis, mvis::ModelVis, mu_RA::Real, mu_DEC::Real)

Calculate the lnlikelihood using a single loop to do the phase shifting and summation. No extra memory needed for re-copying the model.
"
function lnprob(dvis::DataVis, mvis::ModelVis, mu_RA::Float64, mu_DEC::Float64)
    lnp::Float64 = 0.0
    prefac::ComplexF64 = -2pi * 1.0im * 1e3 * arcsec
    for i=1:length(dvis.VV)
        lnp += abs2(dvis.invsig[i] * (dvis.VV[i] - mvis.VV[i] * exp(prefac * (dvis.uu[i] * mu_RA + dvis.vv[i] * mu_DEC))))
    end
    return -0.5 * lnp
end

"
    lnprob_real(dvis::DataVis, mvis::ModelVis, mu_RA::Float64, mu_DEC::Float64)

Calculate the lnlikelihood using a single loop to do the phase shifting and summation. No extra memory needed for re-copying the model. Uses Euler formula to do the calculation, rather than exponentiating a complex number.
"
function lnprob_real(dreal::DataVisReal, mreal::ModelVisReal, mu)
    lnp::Float64 = 0.0
    val::Float64 = 0.0
    sval::Float64 = 0.0 # sin of val
    cval::Float64 = 0.0 # cos of val

    # dreal and mreal both have the real and imaginary terms stored separately as Float64 vectors.
    prefac::Float64 = -2pi * 1e3 * arcsec

    for i=1:length(dreal.VVr)
        val = prefac * (mu[1] * dreal.uu[i] + mu[2] * dreal.vv[i])
        sval = sin(val)
        cval = cos(val)

        lnp += ((dreal.VVr[i] - cval * mreal.VVr[i] + sval * mreal.VVi[i]) * dreal.invsig[i])^2 + ((dreal.VVi[i] - sval * mreal.VVr[i] - cval * mreal.VVi[i]) * dreal.invsig[i])^2
    end
    return -0.5 * lnp
end

function lnprob_iter(dvis::DataVis, mvis::ModelVis)
    lnp = 0.0
    for i=1:length(dvis.VV)
        lnp += abs2(dvis.invsig[i] * (dvis.VV[i] - mvis.VV[i]))
    end
    return -0.5 * lnp
end

"
    chi2(dvis::DataVis, mvis::ModelVis)

Calculate the ``\\chi^2`` value of the current model. This function is not used in inference (see `lnprob`), but it is sometimes a useful number to report."
function chi2(dvis::DataVis, mvis::ModelVis)
    @assert dvis === mvis.dvis # Using the wrong ModelVis, otherwise!
    return sum(abs2, dvis.invsig .* (dvis.VV - mvis.VV)) # Basic chi2
end

"
    phase_shift!(mvis::ModelVis, mu_RA, mu_DEC)

Given a new model centroid in the image plane (in arcseconds), shift the model
visibilities by corresponding amount."
function phase_shift!(mvis::ModelVis, mu_RA, mu_DEC)

    mu = Float64[mu_RA, mu_DEC] * arcsec # [radians]

    nvis = length(mvis.VV)
    # Go through each visibility and apply the phase shift
    for i=1:nvis
        # Convert from [kλ] to [λ]
        R = Float64[mvis.dvis.uu[i], mvis.dvis.vv[i]] * 1e3 #[λ]
        # Not actually in polar phase form
        shift = exp(-2pi * 1.0im * (R' * mu)[1])
        mvis.VV[i] = mvis.VV[i] * shift
    end
end

function phase_shift(mvis::ModelVis, mu_RA, mu_DEC)
    mv = deepcopy(mvis)
    mv.dvis = mvis.dvis
    phase_shift!(mv, mu_RA, mu_DEC)
    return mv
end

# Shift the (0,0) data phase center to these new RA, DEC centroids. The mu_RA, mu_DEC should be
# the same sign as those in the modelVis phase shift routine.
# This is only necessary if the data is so off centered relative to the RADMC image size
# that you are running into problems with tapering cutting off a majority of the image.
# function phase_shift!(dvis::DataVis, mu_RA, mu_DEC)
#
#     mu = -Float64[mu_RA, mu_DEC] * arcsec # [radians]
#
#     nvis = length(dvis.VV)
#     # Go through each visibility and apply the phase shift
#     for i=1:nvis
#         # Convert from [kλ] to [λ]
#         R = Float64[dvis.dvis.uu[i], dvis.dvis.vv[i]] * 1e3 #[λ]
#         # Not actually in polar phase form
#         shift = exp(-2pi * 1.0im * (R' * mu)[1])
#         dvis.VV[i] = dvis.VV[i] * shift
#     end
# end

"
    phase_shift!(fvis::FullModelVis, mu_RA, mu_DEC)

Given a new model centroid in the image plane (in arcseconds), shift the model
visibilities by corresponding amount."
function phase_shift!(fvis::FullModelVis, mu_RA, mu_DEC)

    mu = Float64[mu_RA, mu_DEC] * arcsec # [radians]

    nu = length(fvis.uu)
    nv = length(fvis.vv)
    # Go through each visibility and apply the phase shift
    for i=1:nu
        for j=1:nv
        # Convert from [kλ] to [λ]
        R = Float64[fvis.uu[i], fvis.vv[j]] * 1e3 #[λ]
        # Not actually in polar phase form
        shift = exp(-2pi * 1.0im * (R' * mu)[1])
        fvis.VV[j,i] = fvis.VV[j,i] * shift
        end
    end
end

"
    transform(img::SkyImage, index::Int=1)

Transform the SkyImage produced by RADMC using FFT."
function transform(img::SkyImage, index::Int=1)

    # By default, select the first channel of any spectral hypercube. This
    # routine can only transform one channel at a time.
    data = img.data[:, :, index]

    lam = img.lams[index]

    # convert ra and dec in [arcsec] to radians, and then take the sin to
    # convert to ll, mm
    # ll = sin(img.ra * arcsec)
    # mm = sin(img.dec * arcsec)

    # Sin is not necessary, since we are using the small angle approximation
    ll = img.ra * arcsec
    mm = img.dec * arcsec

    # number of elements in each array
    nl = length(ll)
    nm = length(mm)

    # find the spacing between the elements
    dl = abs(ll[2] - ll[1]) # [radians]
    dm = abs(mm[2] - mm[1]) # [radians]

    # println("Transform using dl ", dl)

    # determine uv plane coordinates in kλ
    uu = fftshift(fftfreq(nl, dl)) * 1e-3 # [kλ]
    vv = fftshift(fftfreq(nm, dm)) * 1e-3 # [kλ]

    # properly pack the data for input using fftshift to move the component at
    # RA=0,DEC=0 to the first array element: data[1,1]
    # For RADMC images, we'll be moving the RA=~0.5, DEC=~0.5 element
    # (1/2 pixel away), and so there will need to be a corresponding phase
    # shift to correct.

    # We also want to normalize the result by the input array spacings, so that
    # they are directly comparable with the analytic transforms
    # (Numerical Recipes ed. 3, Press, Eqn 12.1.6)
    out = dl * dm * fftshift(fft(fftshift(data)))
    return FullModelVis(lam, uu, vv, out)
end

"
    fillModelVis(vis::RawModelVis)

This function is designed to copy the partial arrays in RawModelVis into a
full image for easy plotting. This means that the u axis can remain the same
but we'll need to make the complex conjugate of the v axis."
function fillModelVis(vis::RawModelVis)

    # The full image will stretch from -(n/2 - 1) to (n/2 -1) and
    # have one less point in the u,v directions

    #Swap the +u and -u quadrants so that we have u increasing from -n/2 to 0 to n/2 - 1
    VV = fftshift(vis.VV, 2)

    # Create the top chunk (v >= 0), by flipping the array such that
    # u=n/2 is in the first row of the array and u=0 on the last row.
    # then trim off the first row and column, so that the array stretches
    # 0 <= v <= n/2 - 1  and -(n/2 - 1) <= u <= (n/2 - 1)
    top = flipdim(VV, 1)[2:end, 2:end]

    # Create the bottom chunk (v < 0) by using the original orientation,
    # and trimming off rows v = 0 and v = -n/2, and column u = n/2.
    # we use the property that the FFT of a real image is Hermitian, and thus use
    # the complex conjugate property to fill in the missing quadrants.
    bottom = conj(flipdim(VV, 2))[2:end-1, 1:end-1]

    # make the coordinates
    uu = fftshift(vis.uu)[2:end] # clip off u=-n/2 frequency

    vv_top = flipdim(vis.vv, 1)[2:end] # clip off v=n/2 frequency
    vv_bottom = -1 .* vis.vv[2:end-1] # clip off v=0 and v=-n/2 frequency
    vv = vcat(vv_top, vv_bottom)

    return FullModelVis(vis.lam, uu, vv, vcat(top, bottom))
end

"
    plan_interpolate(dvis::DataVis, uu::Vector{Float64}, vv::Vector{Float64})

Return a function that is used to interpolate the visibilities, in the
spirit of interpolate_uv but *much* faster.
Closures save time and money!

The uu and vv vectors are the visibility coordinates for the dense visibility *model*.

The point is that if the distance to the source and the size of the sythesized image
are not changing, then we will always be interpolating from a dense Visibility grid that
has the exact same U and V spacings, meaning that the weighting terms used to evaluate the
interpolation for a specific visibility can be cached.

So, using a closure, those weighting terms are calculated once and then cached for further use."
function plan_interpolate(dvis::DataVis, uu::Vector{Float64}, vv::Vector{Float64})

    nvis = length(dvis.VV)
    uinds = Array{UnitRange{Int64}}(undef, nvis)
    vinds = Array{UnitRange{Int64}}(undef, nvis)
    uws = Array{Float64}(undef, 6, nvis) #stored along columns
    vws = Array{Float64}(undef, 6, nvis)

    for i=1:nvis
        u = dvis.uu[i]
        v = dvis.vv[i]
        iu0 = argmin(abs.(u .- uu))
        iv0 = argmin(abs.(v .- vv))

        # now find the relative distance to this nearest grid point (not absolute)
        u0 = u - uu[iu0]
        v0 = v - vv[iv0]

        # determine the uu and vv distance for 3 grid points (could be later taken out)
        du = abs(uu[4] - uu[1])
        dv = abs(vv[4] - vv[1])

        # Calculate the appropriate u and v indexes for the 6 nearest pixels
        # (3 on either side)

        # Are u0 and v0 to the left or the right of the index?
        # we want to index three to the left, three to the right

        # First check that we are still in bounds of the array
        # Check to make sure that at least three grid points exist in all directions
        lenu = length(uu)
        lenv = length(vv)
        @assert iu0 >= 4 "The model visibility does not have the resolution to cover the largest baselines in the dataset."
        @assert iv0 >= 4 "The model visibility does not have the resolution to cover the largest baselines in the dataset."
        @assert lenu - iu0 >= 4
        @assert lenv - iv0 >= 4

        if u0 >= 0.0
            # To the right of the index
            uind = iu0-2:iu0+3
        else
            # To the left of the index
            uind = iu0-3:iu0+2
        end

        if v0 >= 0.0
            # To the right of the index
            vind = iv0-2:iv0+3
        else
            # To the left of the index
            vind = iv0-3:iv0+2
        end

        # Store these slices so we can later index into the visibility array
        uinds[i] = uind
        vinds[i] = vind

        # 3. Calculate the weights corresponding to these 6 nearest pixels (gcffun)
        etau = (uu[uind] .- u)/du
        etav = (vv[vind] .- v)/dv

        uw = gcffun(etau)
        vw = gcffun(etav)

        # 4. Normalization such that it has an area of 1. Divide by w later.
        uw = uw/sum(uw)
        vw = vw/sum(vw)

        uws[:,i] = uw
        vws[:,i] = vw
    end

    # tol = 1e-5 .* ones(Float64, length(uu))
    tol = 1e-5

    # This function inherits all of the variables just defined in this scope (uu, vv)
    function interpolate(data::DataVis, fmvis::FullModelVis)
        # Assert that we calculated the same UU and VV spacings for the FT'ed image, otherwise we did something wrong!
        # The 1e-5 addition is to prevent an undetermined error from the uu = 0.0 point.
        @assert all(abs.((uu .- fmvis.uu) ./ (uu .+ 1e-5)) .< tol)
        @assert all(abs.((vv .- fmvis.vv) ./ (vv .+ 1e-5)) .< tol)

        # output array
        Vmodel = Array{ComplexF64}(undef, nvis)

        for i=1:nvis

            Vdata = fmvis.VV[vinds[i], uinds[i]] # Array is packed like the image

            # Loop over all 36 grid indices and sum to find the interpolation.
            cum::ComplexF64 = 0.0 + 0.0im
            for k=1:6
                for l=1:6
                    cum += uws[k,i] * vws[l,i] * Vdata[l,k] # Array is packed like the image
                end
            end

            Vmodel[i] = cum

        end

        return ModelVis(data, Vmodel)
    end

    return interpolate
end

# called ModGrid in gridding.c (KR code) and in Model.for (MIRIAD)
# Uses spheroidal wave functions to interpolate a model to a (u,v) coordinate.
# u,v are in [kλ]

"
    interpolate_uv(u::Float64, v::Float64, vis::FullModelVis)

Interpolates a dense grid of visibilities (e.g., from FFT of an image) to a specfic (u,v) point using spheroidal functions in a band-limited manner designed to reduce aliasing.
"
function interpolate_uv(u::Float64, v::Float64, vis::FullModelVis)

    # Note that vis.uu goes from positive to negative (East-West)
    # and vis.vv goes from negative to positive (North-South)

    # 1. Find the nearest gridpoint in the FFT'd image.
    iu0 = argmin(abs.(u .- vis.uu))
    iv0 = argmin(abs.(v .- vis.vv))

    # now find the relative distance from (u,v) to this nearest grid point (not absolute)
    u0 = u - vis.uu[iu0]
    v0 = v - vis.vv[iv0]

    # determine the uu and vv distance for 3 grid points (could be later taken out)
    du = abs.(vis.uu[4] - vis.uu[1])
    dv = abs.(vis.vv[4] - vis.vv[1])

    # 2. Calculate the appropriate u and v indexes for the 6 nearest pixels
    # (3 on either side)

    # Are u0 and v0 to the left or the right of the index?
    # we want to index three to the left, three to the right

    # First check that our (u,v) point still exists within the appropriate margins of the
    # dense visibility array
    # This is to make sure that at least three grid points exist in all directions
    # If this fails, this means that the synthesized image is too large compared to the sampled visibilities (meaning that the dense FFT grid is too small).
    # The max u,v sampled is 2/dRA or 2/dDec. This means that if dRA or dDEC is too large, then this
    # will fail
    lenu = length(vis.uu)
    lenv = length(vis.vv)
    @assert iu0 >= 4
    @assert iv0 >= 4
    @assert lenu - iu0 >= 4
    @assert lenv - iv0 >= 4

    if u0 >= 0.0
        # To the right of the index
        uind = iu0-2:iu0+3
    else
        # To the left of the index
        uind = iu0-3:iu0+2
    end

    if v0 >= 0.0
        # To the right of the index
        vind = iv0-2:iv0+3
    else
        # To the left of the index
        vind = iv0-3:iv0+2
    end

    etau = (vis.uu[uind] .- u)/du
    etav = (vis.vv[vind] .- v)/dv
    VV = vis.VV[vind, uind] # Array is packed like the image

    # 3. Calculate the weights corresponding to these 6 nearest pixels (gcffun)
    uw = gcffun(etau)
    vw = gcffun(etav)

    # 4. Normalization such that it has an area of 1. Divide by w later.
    w = sum(uw) * sum(vw)

    # 5. Loop over all 36 grid indices and sum to find the interpolation.
    cumulative::ComplexF64 = 0.0 + 0.0im
    for i=1:6
        for j=1:6
            cumulative += uw[i] * vw[j] * VV[j,i] # Array is packed like the image
        end
    end

    cumulative = cumulative/w

    return cumulative
end

"
    rfftfreq(n::Int, d::Float64)

Return the frequencies corresponding to the output of the real FFT. After numpy.fft.rfftfreq
    f = [0, 1, ...,     n/2-1,     n/2] / (d*n)   if n is even
    f = [0, 1, ..., (n-1)/2-1, (n-1)/2] / (d*n)   if n is odd."
function rfftfreq(n::Int, d::Float64)
    # n even
    if n % 2 == 0
        #Array contains n/2 + 1 elements from [0, n/2]
        return linspace(0, n/2, int(n/2) + 1) / (d * n)

    # n odd
    else
        #Array contains n/2 elements from [0, (n-1)/2]
        return linspace(0, (n - 1)/2, int(n/2)) / (d * n)
    end
end

"
    fftfreq(n::Int, d::Float64)

After numpy.fft.fftfreq

    f = [0, 1, ...,   n/2-1,     -n/2, ..., -1] / (d*n)   if n is even
    f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (d*n)   if n is odd."
function fftfreq(n::Int, d::Float64)
    val = 1.0/(n * d)
    results = Array{Float64}(undef, n)
    N = floor(Int, (n  - 1)/2) + 1

    p1 = Float64[i for i=0:(N-1)]
    results[1:N] = p1

    p2 = Float64[i for i=-(floor(Int, n/2)):-1]
    results[N+1:end] = p2

    return results * val
end

end