text:
const libfasttransforms = find_library("libfasttransforms", [joinpath(dirname(@__DIR__), "deps")])

if libfasttransforms ≡ nothing || length(libfasttransforms) == 0
    error("FastTransforms is not properly installed. Please run Pkg.build(\"FastTransforms\") ",
          "and restart Julia.")
end

function ft_fftw_plan_with_nthreads(n::Integer)
    ccall((:ft_fftw_plan_with_nthreads, libfasttransforms), Cvoid, (Cint, ), n)
end

function __init__()
    ccall((:ft_fftw_init_threads, libfasttransforms), Cint, ())
    ft_fftw_plan_with_nthreads(ceil(Int, Sys.CPU_THREADS/2))
end

"""
    mpfr_t <: AbstractFloat

A Julia struct that exactly matches `mpfr_t`.
"""
struct mpfr_t <: AbstractFloat
    prec::Clong
    sign::Cint
    exp::Clong
    d::Ptr{Limb}
end

mpfr_t(x::BigFloat) = mpfr_t(x.prec, x.sign, x.exp, x.d)

function BigFloat(x::mpfr_t)
    nb = ccall((:mpfr_custom_get_size,:libmpfr), Csize_t, (Clong,), precision(BigFloat))
    nb = (nb + Core.sizeof(Limb) - 1) ÷ Core.sizeof(Limb) # align to number of Limb allocations required for this
    str = unsafe_string(Ptr{UInt8}(x.d), nb * Core.sizeof(Limb))
    _BigFloat(x.prec, x.sign, x.exp, str)
end

set_num_threads(n::Integer) = ccall((:ft_set_num_threads, libfasttransforms), Cvoid, (Cint, ), n)

const LEG2CHEB              = 0
const CHEB2LEG              = 1
const ULTRA2ULTRA           = 2
const JAC2JAC               = 3
const LAG2LAG               = 4
const JAC2ULTRA             = 5
const ULTRA2JAC             = 6
const JAC2CHEB              = 7
const CHEB2JAC              = 8
const ULTRA2CHEB            = 9
const CHEB2ULTRA           = 10
const SPHERE               = 11
const SPHEREV              = 12
const DISK                 = 13
const TRIANGLE             = 14
const TETRAHEDRON          = 15
const SPHERESYNTHESIS      = 16
const SPHEREANALYSIS       = 17
const SPHEREVSYNTHESIS     = 18
const SPHEREVANALYSIS      = 19
const DISKSYNTHESIS        = 20
const DISKANALYSIS         = 21
const TRIANGLESYNTHESIS    = 22
const TRIANGLEANALYSIS     = 23
const TETRAHEDRONSYNTHESIS = 24
const TETRAHEDRONANALYSIS  = 25

let k2s = Dict(LEG2CHEB             => "Legendre--Chebyshev",
               CHEB2LEG             => "Chebyshev--Legendre",
               ULTRA2ULTRA          => "ultraspherical--ultraspherical",
               JAC2JAC              => "Jacobi--Jacobi",
               LAG2LAG              => "Laguerre--Laguerre",
               JAC2ULTRA            => "Jacobi--ultraspherical",
               ULTRA2JAC            => "ultraspherical--Jacobi",
               JAC2CHEB             => "Jacobi--Chebyshev",
               CHEB2JAC             => "Chebyshev--Jacobi",
               ULTRA2CHEB           => "ultraspherical--Chebyshev",
               CHEB2ULTRA           => "Chebyshev--ultraspherical",
               SPHERE               => "Spherical harmonic--Fourier",
               SPHEREV              => "Spherical vector field--Fourier",
               DISK                 => "Zernike--Chebyshev×Fourier",
               TRIANGLE             => "Proriol--Chebyshev²",
               TETRAHEDRON          => "Proriol--Chebyshev³",
               SPHERESYNTHESIS      => "FFTW Fourier synthesis on the sphere",
               SPHEREANALYSIS       => "FFTW Fourier analysis on the sphere",
               SPHEREVSYNTHESIS     => "FFTW Fourier synthesis on the sphere (vector field)",
               SPHEREVANALYSIS      => "FFTW Fourier analysis on the sphere (vector field)",
               DISKSYNTHESIS        => "FFTW Chebyshev×Fourier synthesis on the disk",
               DISKANALYSIS         => "FFTW Chebyshev×Fourier analysis on the disk",
               TRIANGLESYNTHESIS    => "FFTW Chebyshev synthesis on the triangle",
               TRIANGLEANALYSIS     => "FFTW Chebyshev analysis on the triangle",
               TETRAHEDRONSYNTHESIS => "FFTW Chebyshev synthesis on the tetrahedron",
               TETRAHEDRONANALYSIS  => "FFTW Chebyshev analysis on the tetrahedron")
    global kind2string
    kind2string(k::Integer) = k2s[Int(k)]
end

struct ft_plan_struct end

mutable struct FTPlan{T, N, K}
    plan::Ptr{ft_plan_struct}
    n::Int
    l::Int
    m::Int
    function FTPlan{T, N, K}(plan::Ptr{ft_plan_struct}, n::Int) where {T, N, K}
        p = new(plan, n)
        finalizer(destroy_plan, p)
        p
    end
    function FTPlan{T, N, K}(plan::Ptr{ft_plan_struct}, n::Int, m::Int) where {T, N, K}
        p = new(plan, n, -1, m)
        finalizer(destroy_plan, p)
        p
    end
    function FTPlan{T, N, K}(plan::Ptr{ft_plan_struct}, n::Int, l::Int, m::Int) where {T, N, K}
        p = new(plan, n, l, m)
        finalizer(destroy_plan, p)
        p
    end
end

eltype(p::FTPlan{T}) where {T} = T
ndims(p::FTPlan{T, N}) where {T, N} = N
show(io::IO, p::FTPlan{T, 1, K}) where {T, K} = print(io, "FastTransforms ", kind2string(K), " plan for $(p.n)-element array of ", T)
show(io::IO, p::FTPlan{T, 2, SPHERE}) where T = print(io, "FastTransforms ", kind2string(SPHERE), " plan for $(p.n)×$(2p.n-1)-element array of ", T)
show(io::IO, p::FTPlan{T, 2, SPHEREV}) where T = print(io, "FastTransforms ", kind2string(SPHEREV), " plan for $(p.n)×$(2p.n-1)-element array of ", T)
show(io::IO, p::FTPlan{T, 2, DISK}) where T = print(io, "FastTransforms ", kind2string(DISK), " plan for $(p.n)×$(4p.n-3)-element array of ", T)
show(io::IO, p::FTPlan{T, 2, TRIANGLE}) where T = print(io, "FastTransforms ", kind2string(TRIANGLE), " plan for $(p.n)×$(p.n)-element array of ", T)
show(io::IO, p::FTPlan{T, 3, TETRAHEDRON}) where T = print(io, "FastTransforms ", kind2string(TETRAHEDRON), " plan for $(p.n)×$(p.n)×$(p.n)-element array of ", T)
show(io::IO, p::FTPlan{T, 2, K}) where {T, K} = print(io, "FastTransforms plan for ", kind2string(K), " for $(p.n)×$(p.m)-element array of ", T)
show(io::IO, p::FTPlan{T, 3, K}) where {T, K} = print(io, "FastTransforms plan for ", kind2string(K), " for $(p.n)×$(p.l)×$(p.m)-element array of ", T)

function checksize(p::FTPlan{T}, x::Array{T}) where T
    if p.n != size(x, 1)
        throw(DimensionMismatch("FTPlan has dimensions $(p.n) × $(p.n), x has leading dimension $(size(x, 1))"))
    end
end

for K in (SPHERE, SPHEREV, DISK)
    @eval function checksize(p::FTPlan{T, 2, $K}, x::Matrix{T}) where T
        if p.n != size(x, 1)
            throw(DimensionMismatch("FTPlan has dimensions $(p.n) × $(p.n), x has leading dimension $(size(x, 1))"))
        end
        if iseven(size(x, 2))
            throw(DimensionMismatch("This FTPlan only operates on arrays with an odd number of columns."))
        end
    end
end

unsafe_convert(::Type{Ptr{ft_plan_struct}}, p::FTPlan) = p.plan
unsafe_convert(::Type{Ptr{mpfr_t}}, p::FTPlan) = unsafe_convert(Ptr{mpfr_t}, p.plan)

destroy_plan(p::FTPlan{Float32, 1}) = ccall((:ft_destroy_tb_eigen_FMMf, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
destroy_plan(p::FTPlan{Float64, 1}) = ccall((:ft_destroy_tb_eigen_FMM, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
destroy_plan(p::FTPlan{BigFloat, 1}) = ccall((:ft_mpfr_destroy_plan, libfasttransforms), Cvoid, (Ptr{mpfr_t}, Cint), p, p.n)
destroy_plan(p::FTPlan{Float64, 2}) = ccall((:ft_destroy_harmonic_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
destroy_plan(p::FTPlan{Float64, 3}) = ccall((:ft_destroy_tetrahedral_harmonic_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
destroy_plan(p::FTPlan{Float64, 2, SPHERESYNTHESIS}) = ccall((:ft_destroy_sphere_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
destroy_plan(p::FTPlan{Float64, 2, SPHEREANALYSIS}) = ccall((:ft_destroy_sphere_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
destroy_plan(p::FTPlan{Float64, 2, SPHEREVSYNTHESIS}) = ccall((:ft_destroy_sphere_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
destroy_plan(p::FTPlan{Float64, 2, SPHEREVANALYSIS}) = ccall((:ft_destroy_sphere_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
destroy_plan(p::FTPlan{Float64, 2, DISKSYNTHESIS}) = ccall((:ft_destroy_disk_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
destroy_plan(p::FTPlan{Float64, 2, DISKANALYSIS}) = ccall((:ft_destroy_disk_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
destroy_plan(p::FTPlan{Float64, 2, TRIANGLESYNTHESIS}) = ccall((:ft_destroy_triangle_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
destroy_plan(p::FTPlan{Float64, 2, TRIANGLEANALYSIS}) = ccall((:ft_destroy_triangle_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
destroy_plan(p::FTPlan{Float64, 3, TETRAHEDRONSYNTHESIS}) = ccall((:ft_destroy_tetrahedron_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
destroy_plan(p::FTPlan{Float64, 3, TETRAHEDRONANALYSIS}) = ccall((:ft_destroy_tetrahedron_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)

struct AdjointFTPlan{T, S}
    parent::S
end

AdjointFTPlan(p::FTPlan) = AdjointFTPlan{eltype(p), typeof(p)}(p)

adjoint(p::FTPlan) = AdjointFTPlan(p)
adjoint(p::AdjointFTPlan) = p.parent

eltype(p::AdjointFTPlan{T, S}) where {T, S} = T
ndims(p::AdjointFTPlan{T, S}) where {T, S} = ndims(p.parent)
function show(io::IO, p::AdjointFTPlan{T, S}) where {T, S}
    print(io, "Adjoint ")
    show(io, p.parent)
end

checksize(p::AdjointFTPlan, x) = checksize(p.parent, x)

unsafe_convert(::Type{Ptr{ft_plan_struct}}, p::AdjointFTPlan{T, FTPlan{T, N, K}}) where {T, N, K} = unsafe_convert(Ptr{ft_plan_struct}, p.parent)
unsafe_convert(::Type{Ptr{mpfr_t}}, p::AdjointFTPlan{T, FTPlan{T, N, K}}) where {T, N, K} = unsafe_convert(Ptr{mpfr_t}, p.parent)

struct TransposeFTPlan{T, S}
    parent::S
end

TransposeFTPlan(p::FTPlan) = TransposeFTPlan{eltype(p), typeof(p)}(p)

transpose(p::FTPlan) = TransposeFTPlan(p)
transpose(p::TransposeFTPlan) = p.parent

eltype(p::TransposeFTPlan{T, S}) where {T, S} = T
ndims(p::TransposeFTPlan{T, S}) where {T, S} = ndims(p.parent)
function show(io::IO, p::TransposeFTPlan{T, S}) where {T, S}
    print(io, "Transpose ")
    show(io, p.parent)
end

checksize(p::TransposeFTPlan, x) = checksize(p.parent, x)

unsafe_convert(::Type{Ptr{ft_plan_struct}}, p::TransposeFTPlan{T, FTPlan{T, N, K}}) where {T, N, K} = unsafe_convert(Ptr{ft_plan_struct}, p.parent)
unsafe_convert(::Type{Ptr{mpfr_t}}, p::TransposeFTPlan{T, FTPlan{T, N, K}}) where {T, N, K} = unsafe_convert(Ptr{mpfr_t}, p.parent)

for f in (:leg2cheb, :cheb2leg, :ultra2ultra, :jac2jac,
          :lag2lag, :jac2ultra, :ultra2jac, :jac2cheb,
          :cheb2jac, :ultra2cheb, :cheb2ultra,
          :sph2fourier, :sphv2fourier, :disk2cxf, :tri2cheb, :tet2cheb)
    plan_f = Symbol("plan_", f)
    @eval begin
        $plan_f(x::AbstractArray{T}, y...; z...) where T = $plan_f(T, size(x, 1), y...; z...)
        $plan_f(::Type{Complex{T}}, y...; z...) where T <: Real = $plan_f(T, y...; z...)
        $f(x::AbstractArray, y...; z...) = $plan_f(x, y...; z...)*x
    end
end

for (f, plan_f) in ((:fourier2sph, :plan_sph2fourier), (:fourier2sphv, :plan_sphv2fourier),
                    (:cxf2disk2, :plan_disk2cxf), (:cheb2tri, :plan_tri2cheb),
                    (:cheb2tet, :plan_tet2cheb))
    @eval begin
        $f(x::AbstractArray, y...; z...) = $plan_f(x, y...; z...)\x
    end
end

function plan_leg2cheb(::Type{Float32}, n::Integer; normleg::Bool=false, normcheb::Bool=false)
    plan = ccall((:ft_plan_legendre_to_chebyshevf, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint), normleg, normcheb, n)
    return FTPlan{Float32, 1, LEG2CHEB}(plan, n)
end

function plan_cheb2leg(::Type{Float32}, n::Integer; normcheb::Bool=false, normleg::Bool=false)
    plan = ccall((:ft_plan_chebyshev_to_legendref, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint), normcheb, normleg, n)
    return FTPlan{Float32, 1, CHEB2LEG}(plan, n)
end

function plan_ultra2ultra(::Type{Float32}, n::Integer, λ, μ; norm1::Bool=false, norm2::Bool=false)
    plan = ccall((:ft_plan_ultraspherical_to_ultrasphericalf, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Float32, Float32), norm1, norm2, n, λ, μ)
    return FTPlan{Float32, 1, ULTRA2ULTRA}(plan, n)
end

function plan_jac2jac(::Type{Float32}, n::Integer, α, β, γ, δ; norm1::Bool=false, norm2::Bool=false)
    plan = ccall((:ft_plan_jacobi_to_jacobif, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Float32, Float32, Float32, Float32), norm1, norm2, n, α, β, γ, δ)
    return FTPlan{Float32, 1, JAC2JAC}(plan, n)
end

function plan_lag2lag(::Type{Float32}, n::Integer, α, β; norm1::Bool=false, norm2::Bool=false)
    plan = ccall((:ft_plan_laguerre_to_laguerref, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Float32, Float32), norm1, norm2, n, α, β)
    return FTPlan{Float32, 1, LAG2LAG}(plan, n)
end

function plan_jac2ultra(::Type{Float32}, n::Integer, α, β, λ; normjac::Bool=false, normultra::Bool=false)
    plan = ccall((:ft_plan_jacobi_to_ultrasphericalf, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Float32, Float32, Float32), normjac, normultra, n, α, β, λ)
    return FTPlan{Float32, 1, JAC2ULTRA}(plan, n)
end

function plan_ultra2jac(::Type{Float32}, n::Integer, λ, α, β; normultra::Bool=false, normjac::Bool=false)
    plan = ccall((:ft_plan_ultraspherical_to_jacobif, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Float32, Float32, Float32), normultra, normjac, n, λ, α, β)
    return FTPlan{Float32, 1, ULTRA2JAC}(plan, n)
end

function plan_jac2cheb(::Type{Float32}, n::Integer, α, β; normjac::Bool=false, normcheb::Bool=false)
    plan = ccall((:ft_plan_jacobi_to_chebyshevf, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Float32, Float32), normjac, normcheb, n, α, β)
    return FTPlan{Float32, 1, JAC2CHEB}(plan, n)
end

function plan_cheb2jac(::Type{Float32}, n::Integer, α, β; normcheb::Bool=false, normjac::Bool=false)
    plan = ccall((:ft_plan_chebyshev_to_jacobif, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Float32, Float32), normcheb, normjac, n, α, β)
    return FTPlan{Float32, 1, CHEB2JAC}(plan, n)
end

function plan_ultra2cheb(::Type{Float32}, n::Integer, λ; normultra::Bool=false, normcheb::Bool=false)
    plan = ccall((:ft_plan_ultraspherical_to_chebyshevf, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Float32), normultra, normcheb, n, λ)
    return FTPlan{Float32, 1, ULTRA2CHEB}(plan, n)
end

function plan_cheb2ultra(::Type{Float32}, n::Integer, λ; normcheb::Bool=false, normultra::Bool=false)
    plan = ccall((:ft_plan_chebyshev_to_ultrasphericalf, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Float32), normcheb, normultra, n, λ)
    return FTPlan{Float32, 1, CHEB2ULTRA}(plan, n)
end

function plan_leg2cheb(::Type{Float64}, n::Integer; normleg::Bool=false, normcheb::Bool=false)
    plan = ccall((:ft_plan_legendre_to_chebyshev, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint), normleg, normcheb, n)
    return FTPlan{Float64, 1, LEG2CHEB}(plan, n)
end

function plan_cheb2leg(::Type{Float64}, n::Integer; normcheb::Bool=false, normleg::Bool=false)
    plan = ccall((:ft_plan_chebyshev_to_legendre, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint), normcheb, normleg, n)
    return FTPlan{Float64, 1, CHEB2LEG}(plan, n)
end

function plan_ultra2ultra(::Type{Float64}, n::Integer, λ, μ; norm1::Bool=false, norm2::Bool=false)
    plan = ccall((:ft_plan_ultraspherical_to_ultraspherical, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Float64, Float64), norm1, norm2, n, λ, μ)
    return FTPlan{Float64, 1, ULTRA2ULTRA}(plan, n)
end

function plan_jac2jac(::Type{Float64}, n::Integer, α, β, γ, δ; norm1::Bool=false, norm2::Bool=false)
    plan = ccall((:ft_plan_jacobi_to_jacobi, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Float64, Float64, Float64, Float64), norm1, norm2, n, α, β, γ, δ)
    return FTPlan{Float64, 1, JAC2JAC}(plan, n)
end

function plan_lag2lag(::Type{Float64}, n::Integer, α, β; norm1::Bool=false, norm2::Bool=false)
    plan = ccall((:ft_plan_laguerre_to_laguerre, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Float64, Float64), norm1, norm2, n, α, β)
    return FTPlan{Float64, 1, LAG2LAG}(plan, n)
end

function plan_jac2ultra(::Type{Float64}, n::Integer, α, β, λ; normjac::Bool=false, normultra::Bool=false)
    plan = ccall((:ft_plan_jacobi_to_ultraspherical, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Float64, Float64, Float64), normjac, normultra, n, α, β, λ)
    return FTPlan{Float64, 1, JAC2ULTRA}(plan, n)
end

function plan_ultra2jac(::Type{Float64}, n::Integer, λ, α, β; normultra::Bool=false, normjac::Bool=false)
    plan = ccall((:ft_plan_ultraspherical_to_jacobi, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Float64, Float64, Float64), normultra, normjac, n, λ, α, β)
    return FTPlan{Float64, 1, ULTRA2JAC}(plan, n)
end

function plan_jac2cheb(::Type{Float64}, n::Integer, α, β; normjac::Bool=false, normcheb::Bool=false)
    plan = ccall((:ft_plan_jacobi_to_chebyshev, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Float64, Float64), normjac, normcheb, n, α, β)
    return FTPlan{Float64, 1, JAC2CHEB}(plan, n)
end

function plan_cheb2jac(::Type{Float64}, n::Integer, α, β; normcheb::Bool=false, normjac::Bool=false)
    plan = ccall((:ft_plan_chebyshev_to_jacobi, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Float64, Float64), normcheb, normjac, n, α, β)
    return FTPlan{Float64, 1, CHEB2JAC}(plan, n)
end

function plan_ultra2cheb(::Type{Float64}, n::Integer, λ; normultra::Bool=false, normcheb::Bool=false)
    plan = ccall((:ft_plan_ultraspherical_to_chebyshev, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Float64), normultra, normcheb, n, λ)
    return FTPlan{Float64, 1, ULTRA2CHEB}(plan, n)
end

function plan_cheb2ultra(::Type{Float64}, n::Integer, λ; normcheb::Bool=false, normultra::Bool=false)
    plan = ccall((:ft_plan_chebyshev_to_ultraspherical, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Float64), normcheb, normultra, n, λ)
    return FTPlan{Float64, 1, CHEB2ULTRA}(plan, n)
end

function plan_leg2cheb(::Type{BigFloat}, n::Integer; normleg::Bool=false, normcheb::Bool=false)
    plan = ccall((:ft_mpfr_plan_legendre_to_chebyshev, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Clong, Int32), normleg, normcheb, n, precision(BigFloat), Base.MPFR.ROUNDING_MODE[])
    return FTPlan{BigFloat, 1, LEG2CHEB}(plan, n)
end

function plan_cheb2leg(::Type{BigFloat}, n::Integer; normcheb::Bool=false, normleg::Bool=false)
    plan = ccall((:ft_mpfr_plan_chebyshev_to_legendre, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Clong, Int32), normcheb, normleg, n, precision(BigFloat), Base.MPFR.ROUNDING_MODE[])
    return FTPlan{BigFloat, 1, CHEB2LEG}(plan, n)
end

function plan_ultra2ultra(::Type{BigFloat}, n::Integer, λ, μ; norm1::Bool=false, norm2::Bool=false)
    plan = ccall((:ft_mpfr_plan_ultraspherical_to_ultraspherical, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Ref{BigFloat}, Ref{BigFloat}, Clong, Int32), norm1, norm2, n, λ, μ, precision(BigFloat), Base.MPFR.ROUNDING_MODE[])
    return FTPlan{BigFloat, 1, ULTRA2ULTRA}(plan, n)
end

function plan_jac2jac(::Type{BigFloat}, n::Integer, α, β, γ, δ; norm1::Bool=false, norm2::Bool=false)
    plan = ccall((:ft_mpfr_plan_jacobi_to_jacobi, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Clong, Int32), norm1, norm2, n, α, β, γ, δ, precision(BigFloat), Base.MPFR.ROUNDING_MODE[])
    return FTPlan{BigFloat, 1, JAC2JAC}(plan, n)
end

function plan_lag2lag(::Type{BigFloat}, n::Integer, α, β; norm1::Bool=false, norm2::Bool=false)
    plan = ccall((:ft_mpfr_plan_laguerre_to_laguerre, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Ref{BigFloat}, Ref{BigFloat}, Clong, Int32), norm1, norm2, n, α, β, precision(BigFloat), Base.MPFR.ROUNDING_MODE[])
    return FTPlan{BigFloat, 1, LAG2LAG}(plan, n)
end

function plan_jac2ultra(::Type{BigFloat}, n::Integer, α, β, λ; normjac::Bool=false, normultra::Bool=false)
    plan = ccall((:ft_mpfr_plan_jacobi_to_ultraspherical, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Clong, Int32), normjac, normultra, n, α, β, λ, precision(BigFloat), Base.MPFR.ROUNDING_MODE[])
    return FTPlan{BigFloat, 1, JAC2ULTRA}(plan, n)
end

function plan_ultra2jac(::Type{BigFloat}, n::Integer, λ, α, β; normultra::Bool=false, normjac::Bool=false)
    plan = ccall((:ft_mpfr_plan_ultraspherical_to_jacobi, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Clong, Int32), normultra, normjac, n, λ, α, β, precision(BigFloat), Base.MPFR.ROUNDING_MODE[])
    return FTPlan{BigFloat, 1, ULTRA2JAC}(plan, n)
end

function plan_jac2cheb(::Type{BigFloat}, n::Integer, α, β; normjac::Bool=false, normcheb::Bool=false)
    plan = ccall((:ft_mpfr_plan_jacobi_to_chebyshev, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Ref{BigFloat}, Ref{BigFloat}, Clong, Int32), normjac, normcheb, n, α, β, precision(BigFloat), Base.MPFR.ROUNDING_MODE[])
    return FTPlan{BigFloat, 1, JAC2CHEB}(plan, n)
end

function plan_cheb2jac(::Type{BigFloat}, n::Integer, α, β; normcheb::Bool=false, normjac::Bool=false)
    plan = ccall((:ft_mpfr_plan_chebyshev_to_jacobi, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Ref{BigFloat}, Ref{BigFloat}, Clong, Int32), normcheb, normjac, n, α, β, precision(BigFloat), Base.MPFR.ROUNDING_MODE[])
    return FTPlan{BigFloat, 1, CHEB2JAC}(plan, n)
end

function plan_ultra2cheb(::Type{BigFloat}, n::Integer, λ; normultra::Bool=false, normcheb::Bool=false)
    plan = ccall((:ft_mpfr_plan_ultraspherical_to_chebyshev, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Ref{BigFloat}, Clong, Int32), normultra, normcheb, n, λ, precision(BigFloat), Base.MPFR.ROUNDING_MODE[])
    return FTPlan{BigFloat, 1, ULTRA2CHEB}(plan, n)
end

function plan_cheb2ultra(::Type{BigFloat}, n::Integer, λ; normcheb::Bool=false, normultra::Bool=false)
    plan = ccall((:ft_mpfr_plan_chebyshev_to_ultraspherical, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Ref{BigFloat}, Clong, Int32), normcheb, normultra, n, λ, precision(BigFloat), Base.MPFR.ROUNDING_MODE[])
    return FTPlan{BigFloat, 1, CHEB2ULTRA}(plan, n)
end

function plan_sph2fourier(::Type{Float64}, n::Integer)
    plan = ccall((:ft_plan_sph2fourier, libfasttransforms), Ptr{ft_plan_struct}, (Cint, ), n)
    return FTPlan{Float64, 2, SPHERE}(plan, n)
end

function plan_sphv2fourier(::Type{Float64}, n::Integer)
    plan = ccall((:ft_plan_sph2fourier, libfasttransforms), Ptr{ft_plan_struct}, (Cint, ), n)
    return FTPlan{Float64, 2, SPHEREV}(plan, n)
end

function plan_disk2cxf(::Type{Float64}, n::Integer)
    plan = ccall((:ft_plan_disk2cxf, libfasttransforms), Ptr{ft_plan_struct}, (Cint, ), n)
    return FTPlan{Float64, 2, DISK}(plan, n)
end

function plan_tri2cheb(::Type{Float64}, n::Integer, α, β, γ)
    plan = ccall((:ft_plan_tri2cheb, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Float64, Float64, Float64), n, α, β, γ)
    return FTPlan{Float64, 2, TRIANGLE}(plan, n)
end

function plan_tet2cheb(::Type{Float64}, n::Integer, α, β, γ, δ)
    plan = ccall((:ft_plan_tet2cheb, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Float64, Float64, Float64, Float64), n, α, β, γ, δ)
    return FTPlan{Float64, 3, TETRAHEDRON}(plan, n)
end

for (fJ, fC, fE, K) in ((:plan_sph_synthesis, :ft_plan_sph_synthesis, :ft_execute_sph_synthesis, SPHERESYNTHESIS),
                    (:plan_sph_analysis, :ft_plan_sph_analysis, :ft_execute_sph_analysis, SPHEREANALYSIS),
                    (:plan_sphv_synthesis, :ft_plan_sphv_synthesis, :ft_execute_sphv_synthesis, SPHEREVSYNTHESIS),
                    (:plan_sphv_analysis, :ft_plan_sphv_analysis, :ft_execute_sphv_analysis, SPHEREVANALYSIS),
                    (:plan_disk_synthesis, :ft_plan_disk_synthesis, :ft_execute_disk_synthesis, DISKSYNTHESIS),
                    (:plan_disk_analysis, :ft_plan_disk_analysis, :ft_execute_disk_analysis, DISKANALYSIS),
                    (:plan_tri_synthesis, :ft_plan_tri_synthesis, :ft_execute_tri_synthesis, TRIANGLESYNTHESIS),
                    (:plan_tri_analysis, :ft_plan_tri_analysis, :ft_execute_tri_analysis, TRIANGLEANALYSIS))
    @eval begin
        $fJ(x::Matrix{T}) where T = $fJ(T, size(x, 1), size(x, 2))
        $fJ(::Type{Complex{T}}, x...) where T <: Real = $fJ(T, x...)
        function $fJ(::Type{Float64}, n::Integer, m::Integer)
            plan = ccall(($(string(fC)), libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint), n, m)
            return FTPlan{Float64, 2, $K}(plan, n, m)
        end
        function lmul!(p::FTPlan{Float64, 2, $K}, x::Matrix{Float64})
            if p.n != size(x, 1) || p.m != size(x, 2)
                throw(DimensionMismatch("FTPlan has dimensions $(p.n) × $(p.m), x has dimensions $(size(x, 1)) × $(size(x, 2))"))
            end
            ccall(($(string(fE)), libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, Ptr{Float64}, Cint, Cint), p, x, size(x, 1), size(x, 2))
            return x
        end
    end
end

plan_tet_synthesis(x::Array{T, 3}) where T = plan_tet_synthesis(T, size(x, 1), size(x, 2), size(x, 3))
plan_tet_synthesis(::Type{Complex{T}}, x...) where T <: Real = plan_tet_synthesis(T, x...)

function plan_tet_synthesis(::Type{Float64}, n::Integer, l::Integer, m::Integer)
    plan = ccall((:ft_plan_tet_synthesis, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint), n, l, m)
    return FTPlan{Float64, 3, TETRAHEDRONSYNTHESIS}(plan, n, l, m)
end

function lmul!(p::FTPlan{Float64, 3, TETRAHEDRONSYNTHESIS}, x::Array{Float64, 3})
    if p.n != size(x, 1) || p.l != size(x, 2) || p.m != size(x, 3)
        throw(DimensionMismatch("FTPlan has dimensions $(p.n) × $(p.l) × $(p.m), x has dimensions $(size(x, 1)) × $(size(x, 2)) × $(size(x, 3))"))
    end
    ccall((:ft_execute_tet_synthesis, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, Ptr{Float64}, Cint, Cint, Cint), p, x, size(x, 1), size(x, 2), size(x, 3))
    return x
end

plan_tet_analysis(x::Array{T, 3}) where T = plan_tet_analysis(T, size(x, 1), size(x, 2), size(x, 3))
plan_tet_analysis(::Type{Complex{T}}, x...) where T <: Real = plan_tet_analysis(T, x...)

function plan_tet_analysis(::Type{Float64}, n::Integer, l::Integer, m::Integer)
    plan = ccall((:ft_plan_tet_analysis, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint), n, l, m)
    return FTPlan{Float64, 3, TETRAHEDRONANALYSIS}(plan, n, l, m)
end

function lmul!(p::FTPlan{Float64, 3, TETRAHEDRONANALYSIS}, x::Array{Float64, 3})
    if p.n != size(x, 1) || p.l != size(x, 2) || p.m != size(x, 3)
        throw(DimensionMismatch("FTPlan has dimensions $(p.n) × $(p.l) × $(p.m), x has dimensions $(size(x, 1)) × $(size(x, 2)) × $(size(x, 3))"))
    end
    ccall((:ft_execute_tet_analysis, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, Ptr{Float64}, Cint, Cint, Cint), p, x, size(x, 1), size(x, 2), size(x, 3))
    return x
end

*(p::FTPlan{T}, x::Array{T}) where T = lmul!(p, deepcopy(x))
*(p::AdjointFTPlan{T}, x::Array{T}) where T = lmul!(p, deepcopy(x))
*(p::TransposeFTPlan{T}, x::Array{T}) where T = lmul!(p, deepcopy(x))
\(p::FTPlan{T}, x::Array{T}) where T = ldiv!(p, deepcopy(x))
\(p::AdjointFTPlan{T}, x::Array{T}) where T = ldiv!(p, deepcopy(x))
\(p::TransposeFTPlan{T}, x::Array{T}) where T = ldiv!(p, deepcopy(x))

for (fJ, fC, elty) in ((:lmul!, :ft_bfmvf, :Float32),
                       (:ldiv!, :ft_bfsvf, :Float32),
                       (:lmul!, :ft_bfmv , :Float64),
                       (:ldiv!, :ft_bfsv , :Float64))
    @eval begin
        function $fJ(p::FTPlan{$elty, 1}, x::Vector{$elty})
            checksize(p, x)
            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Ptr{ft_plan_struct}, Ptr{$elty}), 'N', p, x)
            return x
        end
        function $fJ(p::AdjointFTPlan{$elty, FTPlan{$elty, 1, K}}, x::Vector{$elty}) where K
            checksize(p, x)
            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Ptr{ft_plan_struct}, Ptr{$elty}), 'T', p, x)
            return x
        end
        function $fJ(p::TransposeFTPlan{$elty, FTPlan{$elty, 1, K}}, x::Vector{$elty}) where K
            checksize(p, x)
            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Ptr{ft_plan_struct}, Ptr{$elty}), 'T', p, x)
            return x
        end
    end
end

for (fJ, fC) in ((:lmul!, :ft_mpfr_trmv),
                 (:ldiv!, :ft_mpfr_trsv))
    @eval begin
        function $fJ(p::FTPlan{BigFloat, 1}, x::Vector{BigFloat})
            checksize(p, x)
            xt = deepcopy.(x)
            xc = mpfr_t.(xt)
            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Ptr{mpfr_t}, Cint, Ptr{mpfr_t}, Int32), 'N', p.n, p, p.n, xc, Base.MPFR.ROUNDING_MODE[])
            x .= BigFloat.(xc)
            return x
        end
        function $fJ(p::AdjointFTPlan{BigFloat, FTPlan{BigFloat, 1, K}}, x::Vector{BigFloat}) where K
            checksize(p, x)
            xt = deepcopy.(x)
            xc = mpfr_t.(xt)
            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Ptr{mpfr_t}, Cint, Ptr{mpfr_t}, Int32), 'T', p.parent.n, p, p.parent.n, xc, Base.MPFR.ROUNDING_MODE[])
            x .= BigFloat.(xc)
            return x
        end
        function $fJ(p::TransposeFTPlan{BigFloat, FTPlan{BigFloat, 1, K}}, x::Vector{BigFloat}) where K
            checksize(p, x)
            xt = deepcopy.(x)
            xc = mpfr_t.(xt)
            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Ptr{mpfr_t}, Cint, Ptr{mpfr_t}, Int32), 'T', p.parent.n, p, p.parent.n, xc, Base.MPFR.ROUNDING_MODE[])
            x .= BigFloat.(xc)
            return x
        end
    end
end

for (fJ, fC, elty) in ((:lmul!, :ft_bfmmf, :Float32),
                       (:ldiv!, :ft_bfsmf, :Float32),
                       (:lmul!, :ft_bfmm , :Float64),
                       (:ldiv!, :ft_bfsm , :Float64))
    @eval begin
        function $fJ(p::FTPlan{$elty, 1}, x::Matrix{$elty})
            checksize(p, x)
            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Ptr{ft_plan_struct}, Ptr{$elty}, Cint, Cint), 'N', p, x, size(x, 1), size(x, 2))
            return x
        end
        function $fJ(p::AdjointFTPlan{$elty, FTPlan{$elty, 1, K}}, x::Matrix{$elty}) where K
            checksize(p, x)
            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Ptr{ft_plan_struct}, Ptr{$elty}, Cint, Cint), 'T', p, x, size(x, 1), size(x, 2))
            return x
        end
        function $fJ(p::TransposeFTPlan{$elty, FTPlan{$elty, 1, K}}, x::Matrix{$elty}) where K
            checksize(p, x)
            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Ptr{ft_plan_struct}, Ptr{$elty}, Cint, Cint), 'T', p, x, size(x, 1), size(x, 2))
            return x
        end
    end
end

for (fJ, fC) in ((:lmul!, :ft_mpfr_trmm),
                 (:ldiv!, :ft_mpfr_trsm))
    @eval begin
        function $fJ(p::FTPlan{BigFloat, 1}, x::Matrix{BigFloat})
            checksize(p, x)
            xt = deepcopy.(x)
            xc = mpfr_t.(xt)
            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Ptr{mpfr_t}, Cint, Ptr{mpfr_t}, Cint, Cint, Int32), 'N', p.n, p, p.n, xc, size(x, 1), size(x, 2), Base.MPFR.ROUNDING_MODE[])
            x .= BigFloat.(xc)
            return x
        end
        function $fJ(p::AdjointFTPlan{BigFloat, FTPlan{BigFloat, 1, K}}, x::Matrix{BigFloat}) where K
            checksize(p, x)
            xt = deepcopy.(x)
            xc = mpfr_t.(xt)
            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Ptr{mpfr_t}, Cint, Ptr{mpfr_t}, Cint, Cint, Int32), 'T', p.parent.n, p, p.parent.n, xc, size(x, 1), size(x, 2), Base.MPFR.ROUNDING_MODE[])
            x .= BigFloat.(xc)
            return x
        end
        function $fJ(p::TransposeFTPlan{BigFloat, FTPlan{BigFloat, 1, K}}, x::Matrix{BigFloat}) where K
            checksize(p, x)
            xt = deepcopy.(x)
            xc = mpfr_t.(xt)
            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Ptr{mpfr_t}, Cint, Ptr{mpfr_t}, Cint, Cint, Int32), 'T', p.parent.n, p, p.parent.n, xc, size(x, 1), size(x, 2), Base.MPFR.ROUNDING_MODE[])
            x .= BigFloat.(xc)
            return x
        end
    end
end

for (fJ, fC, K) in ((:lmul!, :ft_execute_sph2fourier, SPHERE),
                    (:ldiv!, :ft_execute_fourier2sph, SPHERE),
                    (:lmul!, :ft_execute_sphv2fourier, SPHEREV),
                    (:ldiv!, :ft_execute_fourier2sphv, SPHEREV),
                    (:lmul!, :ft_execute_disk2cxf, DISK),
                    (:ldiv!, :ft_execute_cxf2disk, DISK),
                    (:lmul!, :ft_execute_tri2cheb, TRIANGLE),
                    (:ldiv!, :ft_execute_cheb2tri, TRIANGLE))
    @eval begin
        function $fJ(p::FTPlan{Float64, 2, $K}, x::Matrix{Float64})
            checksize(p, x)
            ccall(($(string(fC)), libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, Ptr{Float64}, Cint, Cint), p, x, size(x, 1), size(x, 2))
            return x
        end
    end
end

function lmul!(p::FTPlan{Float64, 3, TETRAHEDRON}, x::Array{Float64, 3})
    checksize(p, x)
    ccall((:ft_execute_tet2cheb, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, Ptr{Float64}, Cint, Cint, Cint), p, x, size(x, 1), size(x, 2), size(x, 3))
    return x
end

function ldiv!(p::FTPlan{Float64, 3, TETRAHEDRON}, x::Array{Float64, 3})
    checksize(p, x)
    ccall((:ft_execute_cheb2tet, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, Ptr{Float64}, Cint, Cint, Cint), p, x, size(x, 1), size(x, 2), size(x, 3))
    return x
end

*(p::FTPlan{T}, x::Array{Complex{T}}) where T = lmul!(p, deepcopy(x))
*(p::AdjointFTPlan{T}, x::Array{Complex{T}}) where T = lmul!(p, deepcopy(x))
*(p::TransposeFTPlan{T}, x::Array{Complex{T}}) where T = lmul!(p, deepcopy(x))
\(p::FTPlan{T}, x::Array{Complex{T}}) where T = ldiv!(p, deepcopy(x))
\(p::AdjointFTPlan{T}, x::Array{Complex{T}}) where T = ldiv!(p, deepcopy(x))
\(p::TransposeFTPlan{T}, x::Array{Complex{T}}) where T = ldiv!(p, deepcopy(x))

for fJ in (:lmul!, :ldiv!)
    @eval begin
        function $fJ(p::FTPlan{T}, x::AbstractArray{Complex{T}}) where T
            x .= complex.($fJ(p, real(x)), $fJ(p, imag(x)))
            return x
        end
        function $fJ(p::AdjointFTPlan{T, FTPlan{T, N, K}}, x::AbstractArray{Complex{T}}) where {T, N, K}
            x .= complex.($fJ(p, real(x)), $fJ(p, imag(x)))
            return x
        end
        function $fJ(p::TransposeFTPlan{T, FTPlan{T, N, K}}, x::AbstractArray{Complex{T}}) where {T, N, K}
            x .= complex.($fJ(p, real(x)), $fJ(p, imag(x)))
            return x
        end
    end
end


text:
"""
A = JopShift(spc, b)
A = JopShift(b)

A(x) = x + b

where `spc::JetSpace` is the domain/range of `A`, and `b::Array` or `b::Number`
is the affine translation or shift.  If `b<:Number`, then shift is constant
and one must specify `spc`.

# Examples
```julia
A = JopShift(spc, b)
```
where `spc::JetSpace` is the domain and range of `A`, and `b<:Array` or
`b<:Number`

```julia
A = JopShift(b)
```
where `b<:AbstractArray`.  The domain and range of `A` are determined by the
size and type of `b`.

```julia
A = JopShift([1.0, 2.0, 3.0])
m = ones(domain(A))
d = A*m # d = [2.0 ; 3.0 ; 4.0]
```
"""
JopShift(spc::JetAbstractSpace, shift) = JopNl(f! = JopShift_f!, df! = JopShift_df!, dom = spc, rng = spc, s = (shift=shift,))
JopShift(shift::AbstractArray{T}) where {T} = JopShift( JetSpace(T, size(shift)),shift)
export JopShift

JopShift_f!(d, m; shift, kwargs...) = d .=  m .+ shift
JopShift_df!(δd, δm; mₒ, shift, kwargs...) = δd .=  δm #Derivative is identity

ext: jl


text:
# Test SortedVector
v = SortedVector(collect([1.; 10.; 10:15.]))
A = AxisArray(reshape(1:16, 8, 2), v, [:a, :b])
@test A[ClosedInterval(8.,12.), :] == A[2:5, :]
@test A[1. .. 1., :] == A[1:1, :]
@test A[10. .. 10., :] == A[2:3, :]

# Test SortedVector with a hierarchical index (indexed using Tuples)
srand(1234)
data = reshape(1.:40., 20, 2)
v = collect(zip([:a, :b, :c][rand(1:3,20)], [:x,:y][rand(1:2,20)], [:x,:y][rand(1:2,20)]))
idx = sortperm(v)
A = AxisArray(data[idx,:], SortedVector(v[idx]), [:a, :b])
@test A[:b, :] == A[5:12, :]
@test A[[:a,:c], :] == A[[1:4;13:end], :]
@test A[(:a,:y), :] == A[2:4, :]
@test A[(:c,:y,:y), :] == A[16:end, :]
@test A[ClosedInterval(:a,:b), :] == A[1:12, :]
@test A[ClosedInterval((:a,:x),(:b,:x)), :] ==  A[1:9, :]
@test A[[ClosedInterval((:a,:x),(:b,:x)),:c], :] ==  A[[1:9;13:end], :]


text:
using Pkg: @pkg_str
pkg"activate ."

#---

using Cassette
using JuliaInterpreter: JuliaInterpreter, BreakpointSignature, BreakpointRef
using MagneticReadHead
using Debugger
using Dates

using JuliaInterpreter: Frame, enter_call, get_return
@inline function run_interpretted(f, args...; kwargs...)
    fr = enter_call(f, args...)
    cmd = :c
    while true
        ret = JuliaInterpreter.debug_command(JuliaInterpreter.finish_and_return!, fr, cmd)
        ret == nothing && break
        fr, pc = ret
    end
    return get_return(fr)
end

Cassette.@context MixModeDebugCtx

function setup_overdubs(breakpoints=JuliaInterpreter._breakpoints)
    for bp in breakpoints
        create_overdub(bp)
    end
end

function create_overdub(bp)
    ft = MagneticReadHead.functiontypeof(bp.f)
    #TODO: select based on all the method sig info, not just the function
    @eval @inline function Cassette.overdub(::MixModeDebugCtx, f::$ft, args...)
        run_interpretted(f, args...)
    end
end

function run_mixedmode(f, args)
    Cassette.recurse(MixModeDebugCtx(), f, args)
end
##############

using Base: gc_num, time_ns, GC_Diff, gc_alloc_count, time_print, show_unquoted
macro time_trial(expr)
    run = quote
        GC.gc()
        local stats = gc_num()
        local elapsedtime = time_ns()
        local val = $(esc(expr))
        elapsedtime = time_ns() - elapsedtime
        local diff = GC_Diff(gc_num(), stats)
        time_print(elapsedtime, diff.allocd, diff.total_time,
                   gc_alloc_count(diff))
        println()
    end

    show_op = :(println($(sprint(show_unquoted, expr))))

    return quote
        $show_op
        print(" - 1st Run: ")
        $run
        print(" - 2nd Run: ")
        $run
        print(" - 3nd Run: ")
        $run
        println()
    end
end

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

# Setup:
function resetup()
    # Reinclude so that it recompiles
    include("../test/demo_funs.jl")
    JuliaInterpreter.remove()
    @breakpoint exp(1.0)
    setup_overdubs()  # for MixedModeDebugger
end

const x = rand(1_000, 500)
# Timing
resetup()
@time_trial winter(x)
resetup()
@time_trial run_interpretted(winter, x)
resetup()
@time_trial run_mixedmode(winter, x)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Rclone"
version = v"1.53.2"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/rclone/rclone.git", "30dae752c9258f6aca9521a205696648fde4a98e")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd rclone/

make

# install manually as `make install` doesn't include $exeext
install -d ${bindir}
install -t ${bindir} ${GOPATH}/bin/rclone${exeext}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("rclone", :rclone)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; compilers = [:go, :c])

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
"""
Draw a line between two points and add an arrowhead at the end. The arrowhead length will be
the length of the side of the arrow's head, and the arrowhead angle is the angle between the
sloping side of the arrowhead and the arrow's shaft.

    arrow(startpoint::Point, endpoint::Point; linewidth=1.0, arrowheadlength=10, arrowheadangle=pi/8)

Arrows don't use the current linewidth setting (`setline()`); you can specify the linewidth.

It doesn't need stroking/filling, the shaft is `stroke()`d and the head `fill()`ed with the
current color.

"""
function arrow(startpoint::Point, endpoint::Point; linewidth=1.0, arrowheadlength=10, arrowheadangle=pi/8)
  gsave()
  setlinejoin("butt")
  setline(linewidth)
  shaftlength = norm(startpoint, endpoint)
  shaftangle = atan2(startpoint.y - endpoint.y, startpoint.x - endpoint.x)
  arrowheadtopsideangle = shaftangle + arrowheadangle
  # shorten the length so that lines
  # stop before we get to the arrow
  # thus wide shafts won't stick out through the head of the arrow.
  max_undershoot = shaftlength - ((linewidth/2) / tan(arrowheadangle))
  true_arrowheadlength = arrowheadlength * cos(arrowheadangle)
  if true_arrowheadlength < max_undershoot
    ratio = (shaftlength - true_arrowheadlength)/shaftlength
  else
    ratio = max_undershoot/shaftlength
  end
  tox = startpoint.x + (endpoint.x - startpoint.x) * ratio
  toy = startpoint.y + (endpoint.y - startpoint.y) * ratio

  #  in case we ever have to do  start arrows:
  #    tox=endpoint.x
  #    toy=endpoint.y
  #    fromx = startpoint.x + (endpoint.x  - startpoint.x) * (1-ratio)
  #    fromy = startpoint.y + (endpoint.y - startpoint.y) *  (1-ratio)

  fromx = startpoint.x
  fromy = startpoint.y

  # draw the shaft of the arrow
  newpath()
  move(fromx,fromy)
  line(tox,toy)
  stroke()

  # draw the arrowhead
  # this is a bit boring, isn't it
  topx = endpoint.x + cos(arrowheadtopsideangle) * arrowheadlength
  topy = endpoint.y + sin(arrowheadtopsideangle) * arrowheadlength
  arrowheadbottomsideangle = shaftangle - arrowheadangle
  botx = endpoint.x + cos(arrowheadbottomsideangle) * arrowheadlength
  boty = endpoint.y + sin(arrowheadbottomsideangle) * arrowheadlength
  poly([Point(topx,topy), endpoint, Point(botx,boty)], :fill)
  grestore()
end

"""
Draw a curved arrow, an arc centered at `centerpos` starting at `startangle` and ending at
`endangle` with an arrowhead at the end. Angles are measured clockwise from the positive
x-axis.

Arrows don't use the current linewidth setting (`setline()`); you can specify the linewidth.

    arrow(centerpos::Point, radius, startangle, endangle; linewidth=1.0, arrowheadlength=10, arrowheadangle=pi/8)

"""
function arrow(centerpos::Point, radius, startangle, endangle; linewidth=1.0, arrowheadlength=10, arrowheadangle=pi/8)
  gsave()
  setlinejoin("butt")
  setline(linewidth)
  if startangle > endangle
      endangle += 2pi
  end
  # don't bother with them if theyre too small
  if isapprox(startangle, endangle, rtol = 0.1)
    return
  end
  # shorten the length so that lines
  # stop before we get to the arrow
  # thus wide shafts won't stick out through the head of the arrow.
  startpoint = Point(radius * cos(startangle), radius * sin(startangle))
  endpoint   = Point(radius * cos(endangle), radius * sin(endangle))
  arclength = radius * mod2pi(endangle - startangle)
  max_undershoot = arclength - ((linewidth/2) / tan(arrowheadangle))
  true_arrowheadlength = arrowheadlength * cos(arrowheadangle)
  if true_arrowheadlength < max_undershoot
    ratio = (arclength - true_arrowheadlength)/arclength
  else
    ratio = max_undershoot/arclength
  end
  newarclength = arclength * ratio
  newendangle = (newarclength/radius) + startangle
  translate(centerpos)
  # draw the arrow
  newpath()
  move(radius * cos(startangle), radius * sin(startangle))
  arc(0, 0, radius, startangle, newendangle, :stroke)
  closepath()
  # draw head
  # rotation of head should be based on end of shaft, not end of arrow
  newendpoint =              Point(radius * cos(newendangle), radius * sin(newendangle))
  shaftangle =               mod2pi(-pi/2 + atan2(0 - newendpoint.y, 0 - newendpoint.x))
  #
  arrowheadoutersideangle = shaftangle + pi - arrowheadangle
  arrowheadinnersideangle    = shaftangle + pi + arrowheadangle
  topx =                     endpoint.x + cos(arrowheadinnersideangle) * arrowheadlength
  topy =                     endpoint.y + sin(arrowheadinnersideangle) * arrowheadlength
  botx =                     endpoint.x + cos(arrowheadoutersideangle) * arrowheadlength
  boty =                     endpoint.y + sin(arrowheadoutersideangle) * arrowheadlength
  poly([Point(topx,topy), Point(endpoint.x, endpoint.y), Point(botx,boty)], :fill)
  grestore()
end


text:
## Cauchy

# pseudocauchy does not normalize at ∞
pseudostieltjes(S::Space{<:Arc},f,z) = stieltjes(setcanonicaldomain(S),f,mobius(S,z))
pseudohilbert(S::Space{<:Arc},f,z) = hilbert(setdomain(S,ChebyshevInterval()),f,mobius(S,z))

stieltjes(S::Space{<:Arc},f,z) =
    stieltjes(setcanonicaldomain(S),f,mobius(S,z))-stieltjes(setcanonicaldomain(S),f,mobius(S,Inf))
function hilbert(S::Space{<:Arc},f,z)
    y = mobius(S,z)
    abs(imag(y)) < 10E-15 || throw(ArgumentError())
    hilbert(setcanonicaldomain(S),f,real(y))+(1/π)*stieltjes(setcanonicaldomain(S),f,mobius(S,Inf))
end

function PseudoHilbert(sp::JacobiWeight{LS,RR},k::Integer) where {LS,RR<:Arc}
    @assert k==1
    csp=setcanonicaldomain(sp)
    H=Hilbert(csp)
    PseudoHilbertWrapper(SpaceOperator(H,setdomain(domainspace(H),domain(sp)),setdomain(rangespace(H),domain(sp))),k)
end

function Hilbert(sp::JacobiWeight{LS,RR},k::Integer) where {LS,RR<:Arc}
    @assert k==1
    csp=setcanonicaldomain(sp)
    St = Stieltjes(csp,mobius(sp,Inf))
    H=Hilbert(csp)+(1/π)*SpaceOperator(St,csp,ConstantSpace(eltype(St)))
    HilbertWrapper(SpaceOperator(H,setdomain(domainspace(H),domain(sp)),setdomain(rangespace(H),domain(sp))),k)
end

## stieltjesintegral

function stieltjesintegral(sp::Space{<:Arc},w,z)
    g=Fun(setcanonicaldomain(sp),w)*fromcanonicalD(sp)
    stieltjesintegral(g,mobius(sp,z))-
        stieltjesintegral(g,mobius(sp,Inf))+
        sum(Fun(sp,w))*log(z-fromcanonical(sp,Inf))
end

function linesumstieltjesintegral(sp::Space{<:Arc},w,z)
    g=Fun(setcanonicaldomain(sp),w)*abs(fromcanonicalD(sp))
    stieltjesintegral(g,mobius(sp,z))-
        stieltjesintegral(g,mobius(sp,Inf))+
        linesum(Fun(sp,w))*log(z-fromcanonical(sp,Inf))
end

function logkernel(sp::Space{<:Arc},w,z)
    g=Fun(setcanonicaldomain(sp),w)*abs(fromcanonicalD(sp))
    logkernel(g,mobius(sp,z))-
        logkernel(g,mobius(sp,Inf))+
        linesum(Fun(sp,w))*logabs(z-fromcanonical(sp,Inf))/π
end

function SingularIntegral(S::JacobiWeight{JW,RR},k::Integer) where {JW,RR<:Arc}
    d=domain(S)
    if k==0
        tol=1E-15
        # the mapped logkernel
        csp=setcanonicaldomain(S)
        Σ=SingularIntegral(csp,0)
        M=Multiplication(abs(fromcanonicalD(d,Fun(identity,csp))),csp)

        z∞=mobius(d,Inf)
        cnst=Vector{Float64}()
        for j = 1:10000
            push!(cnst,logkernel(Fun(csp,[zeros(j-1);1.]),z∞))
            if j ≥ 3 && norm(cnst[end-2:end]) < tol
                break
            end
        end
        L∞=FiniteOperator(cnst',csp,ConstantSpace(eltype(cnst)))

        x=Fun(identity,S)
        SpaceOperator((Σ-L∞)*M,S,setdomain(rangespace(Σ),d)) +
            (logabs(x-fromcanonical(d,Inf))/π)*DefiniteLineIntegral(S)
    else
        # multiply by abs(M')/M' to change to dz to ds
        Mp=fromcanonicalD(d)
        Hilbert(S,k)[setdomain(abs(Mp)/Mp,d)]
    end
end


text:
using Diana
using Test

@testset "Diana" begin
include("cliente.jl");
include("validationast.jl");
include("parser.jl");
include("lexer.jl");
include("execute.jl");

end;


text:
function get_axis_label(i::Integer)
    if 0 <= i <= 9
        return " $(i)"
    elseif 10 <= i <= 99
        return "$(i)"
    else
        return "  "
    end
end

function visualize(io::IO, image::AbstractMatrix{Bool})
    height, width = size(image)

    print(io, "  ")
    for j in 1:width
        print(io, get_axis_label(j))
    end
    println(io)

    for i in 1:height
        print(io, get_axis_label(i))

        for j in 1:width
            if image[i, j]
                print(io, "██")
            else
                if iseven(i + j)
                    print(io, "░░")
                else
                    print(io, "▒▒")
                end
            end
        end
        println(io)
    end

    return nothing
end

visualize(image::AbstractMatrix) = visualize(stdout, image)


text:
#################### Posterior Statistics ####################
"""
    autocor(chn;
        lags=[1, 5, 10, 50],
        demean=true,
        relative=true
        showall=false,
        append_chains=true,
        sections=[:parameters])

Compute the autocorrelation of each parameter for the chain. Setting `append_chains=false` will return a vector of dataframes containing the autocorrelations for each chain.
"""
function autocor(chn::AbstractChains;
        lags::Vector=[1, 5, 10, 50],
        demean::Bool=true,
        relative::Bool=true,
        showall=false,
        append_chains = false,
        sections::Union{Symbol, Vector{Symbol}}=Symbol[:parameters])
    funs = Function[]
    func_names = String[]
    for i in lags
        push!(funs, x -> autocor(x, [i], demean=demean)[1])
        push!(func_names, "lag $i")
    end
    return summarize(chn, funs...;
        func_names = func_names,
        showall = showall,
        append_chains = append_chains,
        name = "Autocorrelation")
end

"""
    cor(chn; showall=false, append_chains=true, sections=[:parameters])

Compute the Pearson correlation matrix for the chain. Setting `append_chains=false` will
return a vector of dataframes containing a correlation matrix for each chain.
"""
function cor(chn::AbstractChains;
        showall=false,
        append_chains=true,
        sections::Union{Symbol, Vector{Symbol}}=Symbol[:parameters])
    df = DataFrame(chn, sections, showall=showall, append_chains=append_chains)

    # Generate
    cormat = if append_chains
        cor(convert(Matrix, df[:, 1:end]))
    else
        [cor(convert(Matrix, i[:, 1:end])) for i in df]
    end
    nms = append_chains ? names(df) : names(df[1])
    columns = if append_chains
        [nms, [cormat[:, i] for i in 1:size(cormat, 2)]...]
    else
        [[nms, [cm[:, i] for i in 1:size(cm, 2)]...] for cm in cormat]
    end
    colnames = vcat([:parameters], nms...)
    df_summary = if append_chains
        ChainDataFrame("Correlation", DataFrame(columns, colnames))
    else
        [ChainDataFrame("Correlation", DataFrame(c, colnames)) for c in columns]
    end
    return df_summary
end

"""
    changerate(chn;
        append_chains=true,
        showall=false,
        sections=[:parameters])

Computes the change rate for the chain. Setting `append_chains=false` will
return a vector of dataframes containing the change rates for each chain.
"""
function changerate(chn::AbstractChains;
    append_chains=true,
    showall=false,
    sections::Union{Symbol, Vector{Symbol}}=Symbol[:parameters])
    # Check for missing values.
    @assert !any(ismissing.(chn.value)) "Change rate comp. doesn't support missing values."

    df = DataFrame(chn, append_chains=true, showall=showall)
    n, p = size(df[1])
    m = length(chains(chn))

    r = zeros(Float64, p, 1, 1)
    r_mv = 0.0
    delta = Array{Bool}(undef, p)
    for k in 1:m
        dfk = df[k]
        prev = convert(Vector, dfk[1,:])
        for i in 2:n
            for j in 1:p
                x = dfk[i, j]
                dx = x != prev[j]
                r[j] += dx
                delta[j] = dx
                prev[j] = x
            end
            r_mv += any(delta)
        end
    end
    vals = round.([r..., r_mv] / (m * (n - 1)), digits = 3)
    rownames = push!(names(df[1]), :multivariate)
    return DataFrame(parameters = rownames, change_rate = vals, name="Change Rate")
end

describe(c::AbstractChains; args...) = describe(stdout, c; args...)

"""
    describe(io::IO,
        c::AbstractChains;
        q = [0.025, 0.25, 0.5, 0.75, 0.975],
        etype=:bm,
        showall=false,
        sections::Union{Symbol, Vector{Symbol}}=Symbol[:parameters],
        args...)

Prints the summary statistics and quantiles for the chain.
"""
function describe(io::IO,
                  c::AbstractChains;
                  q = [0.025, 0.25, 0.5, 0.75, 0.975],
                  etype=:bm,
                  showall::Bool=false,
                  sections::Union{Symbol, Vector{Symbol}}=Symbol[:parameters],
                  args...
                 )
    dfs = [summarystats(c, showall=showall, sections=sections, etype=etype, args...),
           quantile(c, showall=showall, sections=sections, q=q)]
    return dfs
end

function _hpd(x::Vector{T}; alpha::Real=0.05) where {T<:Real}
    n = length(x)
    m = max(1, ceil(Int, alpha * n))

    y = sort(x)
    a = y[1:m]
    b = y[(n - m + 1):n]
    _, i = findmin(b - a)

    return [a[i], b[i]]
end

function hpd(chn::AbstractChains; alpha::Real=0.05,
        append_chains=true,
        showall=false,
        sections::Union{Symbol, Vector{Symbol}}=Symbol[:parameters])
    labels = [:upper, :lower]
    u(x) = _hpd(x, alpha=alpha)[1]
    l(x) = _hpd(x, alpha=alpha)[2]
    return summarize(chn, u, l; func_names = labels, showall=showall, name="HPD")
end

"""
    quantile(chn;
        q=[0.025, 0.25, 0.5, 0.75, 0.975],
        append_chains=true,
        showall=false,
        sections=[:parameters])

Computes the quantiles for each parameter in the chain. Setting `append_chains=false` will
return a vector of dataframes containing the quantiles for each chain.
"""
function quantile(chn::AbstractChains;
        q::Vector=[0.025, 0.25, 0.5, 0.75, 0.975],
        append_chains=true,
        showall=false,
        sections::Union{Symbol, Vector{Symbol}}=Symbol[:parameters])
    # compute quantiles
    funs = Function[]
    func_names = String[]
    for i in q
        push!(funs, x -> quantile(cskip(x), i))
        push!(func_names, "$(string(100*i))%")
    end
    return summarize(chn, funs...;
        func_names=func_names,
        showall=showall,
        sections=sections,
        name = "Quantiles")
end

"""
	ess(chn::AbstractChains;
		showall=false,
		sections::Union{Symbol, Vector{Symbol}}=Symbol[:parameters],
		maxlag = 250)

Compute a chain's number of effective samples. More information can be found in the Gelman et al. (2014) book "Bayesian Data Analysis", or in [this article](https://arxiv.org/abs/1903.08008).
"""
function ess(chn::AbstractChains;
    showall=false,
    sections::Union{Symbol, Vector{Symbol}}=Symbol[:parameters],
    maxlag = 250
    )
	param = showall ? names(chn) : names(chn, sections)
	n_chain_orig = size(chn, 3)

	# Split the chains.
	parameter_vec = Vector(undef, length(param))
	midpoint = Int32(floor(size(chn, 1) / 2))
	for i in 1:length(param)
		parameter_vec[i] = []
		for j in 1:n_chain_orig
            c1 = vec(cskip(chn[1:midpoint, param[i], j].value.data))
            c2 = vec(cskip(chn[midpoint+1:end, param[i], j].value.data))
            push!(parameter_vec[i], c1, c2)
		end
	end

    # Misc allocations.
    m = n_chain_orig * 2
    n = min(length.(parameter_vec[1])...)
    maxlag = min(maxlag, n-1)
    lags = collect(0:maxlag)

    # Preallocate B, W, varhat, and Rhat vectors for each param.
    B = Vector(undef, length(param))
    W = Vector(undef, length(param))
    varhat = Vector(undef, length(param))
    Rhat = Vector(undef, length(param))

    # calculate B, W, varhat, and Rhat for each param.
    for i in 1:length(param)
		draw = parameter_vec[i]
		p = param[i]
        allchain = mean(vcat([d for d in draw]...))
        eachchain = [mean(draw[j]) for j in 1:m]
        s = [sum((draw[j] .- eachchain[j]).^2) / (n-1) for j in 1:m]
        B[i] = (n / (m - 1)) * sum((eachchain .- allchain).^2)
        W[i] = sum(s) / m
        varhat[i] = (n-1)/n * W[i] + (B[i] / n)
        Rhat[i] = sqrt(varhat[i] / W[i])
    end

	V = Vector(undef, length(param))
    ρ = Vector(undef, length(param))
    for p in eachindex(V)
        V[p] = Vector(undef, length(lags))
		ρ[p] = Vector(undef, length(lags))
    end

    # Calculate ρ
    c_autocor = Vector(undef, length(param))
    for i in 1:length(param)
        c_autocor[i] = [0.0]
    end

    for t in eachindex(lags)
        lag = lags[t]
        range1 = lag+1:n
        range2 = 1:(n-lag)
        for i in 1:length(param)
            draw = parameter_vec[i]
			p = param[i]
            z = [draw[j][range1] .- draw[j][range2] for j in 1:m]
			z = sum([zi .^ 2 for zi in z])
            V[i][t] = 1 / (m * (n-lag)) * sum(z)
            autocors = [_autocorrelation(draw[j], lag) for j in 1:m]
			ρ[i][t] = 1 - V[i][t] / (2 * varhat[i])
        end
    end

	# Find first odd positive integer where ρ[p][T+1] + ρ[p][T+2] is negative
    P = Vector(undef, length(param))
    ess = Vector(undef, length(param))
	for i in 1:length(param)
        big_P = 0.0
		ρ_val = Float64.(ρ[i])

        # Big P.
        P[i] = Float64[ρ_val[1]]
        k = tprime = 1
        for tprime in 1:Int(floor((length(lags)/2 - 1)))
            sumvals = ρ_val[2*tprime] + ρ_val[2*tprime+1]
            if sumvals < 0
                break
            else
                push!(P[i], sumvals)
                k = tprime
            end
        end

        # Create monotone.
        P_monotone = [min(P[i][t], P[i][1:t]...) for t in 1:length(P[i])]

        ess[i] = (n*m) / (-1 + 2*sum(P_monotone))
	end

    df = DataFrame(parameters = Symbol.(param), ess = ess, r_hat = Rhat)
	return ChainDataFrame("ESS", df)
end

# this function is sourced from https://github.com/tpapp/MCMCDiagnostics.jl/blob/master/src/MCMCDiagnostics.jl
function _autocorrelation(x::AbstractVector, k::Integer, v = var(x))
    x1 = @view(x[1:(end-k)])
    x2 = @view(x[(1+k):end])
    V = sum((x1 .- x2).^2) / length(x1)
    1 - V / (2*v)
end

"""
    summarystats(chn;
        append_chains=true,
        showall=false,
        sections=[:parameters],
        args...)

Computes the mean, standard deviation, naive standard error, Monte Carlo standard error, and effective sample size for each parameter in the chain. Setting `append_chains=false` will return a vector of dataframes containing the summary statistics for each chain. `args...` is passed to the `msce` function.
"""
function summarystats(chn::MCMCChains.AbstractChains;
        append_chains::Bool=true,
        showall::Bool=false,
        sections::Union{Symbol, Vector{Symbol}}=Symbol[:parameters],
        etype=:bm, args...
    )

    # Make some functions.
    sem(x) = sqrt(var(x) / length(x))
    df_mcse(x) = length(x) < 200 ?
        missing :
        mcse(cskip(x), etype, args...)

    # Store everything.
    funs = [mean∘cskip, std∘cskip, sem∘cskip, df_mcse]
    func_names = [:mean, :std, :naive_se, :mcse]

    # Caluclate ESS separately.
    ess_df = ess(chn, sections=sections, showall=showall).df

    # Summarize.
    summary_df = summarize(chn, funs...;
        sections=sections,
        func_names=func_names,
        showall=showall,
        name="Summary Statistics",
        additional_df = ess_df)

    return summary_df
end


text:
# Now that you know how to find low-risk paths in the cave, you can try to find
# your way out.
#
# The entire cave is actually five times larger in both dimensions than you
# thought; the area you originally scanned is just one tile in a 5x5 tile area
# that forms the full map. Your original map tile repeats to the right and
# downward; each time the tile repeats to the right or downward, all of its
# risk levels are 1 higher than the tile immediately up or left of it. However,
# risk levels above 9 wrap back around to 1. So, if your original map had some
# position with a risk level of 8, then that same position on each of the 25
# total tiles would be as follows:
#
# 8 9 1 2 3
# 9 1 2 3 4
# 1 2 3 4 5
# 2 3 4 5 6
# 3 4 5 6 7
#
# Each single digit above corresponds to the example position with a value of 8
# on the top-left tile. Because the full map is actually five times larger in
# both dimensions, that position appears a total of 25 times, once in each
# duplicated tile, with the values shown above.
#
# Here is the full five-times-as-large version of the first example above, with
# the original map in the top left corner highlighted:
#
# 11637517422274862853338597396444961841755517295286
# 13813736722492484783351359589446246169155735727126
# 21365113283247622439435873354154698446526571955763
# 36949315694715142671582625378269373648937148475914
# 74634171118574528222968563933317967414442817852555
# 13191281372421239248353234135946434524615754563572
# 13599124212461123532357223464346833457545794456865
# 31254216394236532741534764385264587549637569865174
# 12931385212314249632342535174345364628545647573965
# 23119445813422155692453326671356443778246755488935
# 22748628533385973964449618417555172952866628316397
# 24924847833513595894462461691557357271266846838237
# 32476224394358733541546984465265719557637682166874
# 47151426715826253782693736489371484759148259586125
# 85745282229685639333179674144428178525553928963666
# 24212392483532341359464345246157545635726865674683
# 24611235323572234643468334575457944568656815567976
# 42365327415347643852645875496375698651748671976285
# 23142496323425351743453646285456475739656758684176
# 34221556924533266713564437782467554889357866599146
# 33859739644496184175551729528666283163977739427418
# 35135958944624616915573572712668468382377957949348
# 43587335415469844652657195576376821668748793277985
# 58262537826937364893714847591482595861259361697236
# 96856393331796741444281785255539289636664139174777
# 35323413594643452461575456357268656746837976785794
# 35722346434683345754579445686568155679767926678187
# 53476438526458754963756986517486719762859782187396
# 34253517434536462854564757396567586841767869795287
# 45332667135644377824675548893578665991468977611257
# 44961841755517295286662831639777394274188841538529
# 46246169155735727126684683823779579493488168151459
# 54698446526571955763768216687487932779859814388196
# 69373648937148475914825958612593616972361472718347
# 17967414442817852555392896366641391747775241285888
# 46434524615754563572686567468379767857948187896815
# 46833457545794456865681556797679266781878137789298
# 64587549637569865174867197628597821873961893298417
# 45364628545647573965675868417678697952878971816398
# 56443778246755488935786659914689776112579188722368
# 55172952866628316397773942741888415385299952649631
# 57357271266846838237795794934881681514599279262561
# 65719557637682166874879327798598143881961925499217
# 71484759148259586125936169723614727183472583829458
# 28178525553928963666413917477752412858886352396999
# 57545635726865674683797678579481878968159298917926
# 57944568656815567976792667818781377892989248891319
# 75698651748671976285978218739618932984172914319528
# 56475739656758684176786979528789718163989182927419
# 67554889357866599146897761125791887223681299833479
#
# Equipped with the full map, you can now find a path from the top left corner
# to the bottom right corner with the lowest total risk:
#
# 11637517422274862853338597396444961841755517295286
# 13813736722492484783351359589446246169155735727126
# 21365113283247622439435873354154698446526571955763
# 36949315694715142671582625378269373648937148475914
# 74634171118574528222968563933317967414442817852555
# 13191281372421239248353234135946434524615754563572
# 13599124212461123532357223464346833457545794456865
# 31254216394236532741534764385264587549637569865174
# 12931385212314249632342535174345364628545647573965
# 23119445813422155692453326671356443778246755488935
# 22748628533385973964449618417555172952866628316397
# 24924847833513595894462461691557357271266846838237
# 32476224394358733541546984465265719557637682166874
# 47151426715826253782693736489371484759148259586125
# 85745282229685639333179674144428178525553928963666
# 24212392483532341359464345246157545635726865674683
# 24611235323572234643468334575457944568656815567976
# 42365327415347643852645875496375698651748671976285
# 23142496323425351743453646285456475739656758684176
# 34221556924533266713564437782467554889357866599146
# 33859739644496184175551729528666283163977739427418
# 35135958944624616915573572712668468382377957949348
# 43587335415469844652657195576376821668748793277985
# 58262537826937364893714847591482595861259361697236
# 96856393331796741444281785255539289636664139174777
# 35323413594643452461575456357268656746837976785794
# 35722346434683345754579445686568155679767926678187
# 53476438526458754963756986517486719762859782187396
# 34253517434536462854564757396567586841767869795287
# 45332667135644377824675548893578665991468977611257
# 44961841755517295286662831639777394274188841538529
# 46246169155735727126684683823779579493488168151459
# 54698446526571955763768216687487932779859814388196
# 69373648937148475914825958612593616972361472718347
# 17967414442817852555392896366641391747775241285888
# 46434524615754563572686567468379767857948187896815
# 46833457545794456865681556797679266781878137789298
# 64587549637569865174867197628597821873961893298417
# 45364628545647573965675868417678697952878971816398
# 56443778246755488935786659914689776112579188722368
# 55172952866628316397773942741888415385299952649631
# 57357271266846838237795794934881681514599279262561
# 65719557637682166874879327798598143881961925499217
# 71484759148259586125936169723614727183472583829458
# 28178525553928963666413917477752412858886352396999
# 57545635726865674683797678579481878968159298917926
# 57944568656815567976792667818781377892989248891319
# 75698651748671976285978218739618932984172914319528
# 56475739656758684176786979528789718163989182927419
# 67554889357866599146897761125791887223681299833479
#
# The total risk of this path is 315 (the starting position is still never
# entered, so its risk is not counted).
#
# Using the full map, what is the lowest total risk of any path from the top
# left to the bottom right?

# load risk map
miniriskmap = hcat(map(l -> parse.(Int8, split(l, "")), eachline("input.txt"))...)
(miniwidth, miniheight) = size(miniriskmap)

# build whole map
(width, height) = (5*miniwidth, 5*miniheight)
riskmap = Array{Int8}(undef, width, height)
for i in 0:4
  riskmap[(i*miniwidth+1):((i+1)*miniwidth), 1:miniheight] = ((miniriskmap .+ (i-1)) .% 9) .+ 1
end
for i in 1:4
  riskmap[1:width, (i*miniheight+1):((i+1)*miniheight)] = ((riskmap[1:width, 1:miniheight] .+ (i-1)) .% 9) .+ 1
end

# some data structures to track things
cost = zeros(Int, size(riskmap))
queue = [(1,1)]
function update_cost(x, y, c)
  c += riskmap[x,y]
  if cost[x,y] == 0
    # first time setting a cost on this node - add to queue
    push!(queue, (x,y))
  elseif cost[x,y] <= c
    # we found a cheaper route to this node already
    return
  end
  cost[x,y] = c
end

# Dijkstra's
visited = falses(size(riskmap))
while !visited[width, height]
  # find the next node with the smallest tentative cost
  sort!(queue; by=((x,y),) -> cost[x,y], rev=true, alg=InsertionSort)
  (x,y) = pop!(queue)
  visited[x,y] = true

  # update cost of neighbors
  c = cost[x,y]
  x > 1 && !visited[x-1, y] && update_cost(x-1, y, c)
  y > 1 && !visited[x, y-1] && update_cost(x, y-1, c)
  x < width && !visited[x+1, y] && update_cost(x+1, y, c)
  y < height && !visited[x, y+1] && update_cost(x, y+1, c)
end

println("Result: ", cost[width, height])


text:
using SymbolicWedderburn
using PermutationGroups

using SparseArrays

using DynamicPolynomials
using SumOfSquares
using SCS

include(joinpath(@__DIR__, "action_polynomials.jl"))

OPTIMIZER = optimizer_with_attributes(
    SCS.Optimizer,
    "acceleration_lookback" => 0,
    "max_iters" => 10_000,
    "eps" => 2e-6,
    "linear_solver" => SCS.DirectSolver,
)

@polyvar x y z

motzkin = x^4 * y^2 + y^4 * x^2 - 3 * x^2 * y^2 + 1
g = (x^2 + y^2 + 1)
monomial_basis = monomials([x, y], 0:7)

ts, st = let f = motzkin, basis = monomial_basis, m = SOSModel(OPTIMIZER)
    @variable m t >=0
    # @objective m Max t
    @variable m sos SOSPoly(basis)
    @constraint m g*f - t == sos
    optimize!(m)
    @info (m,) termination_status(m) objective_value(m) solve_time(m)
    termination_status(m), solve_time(m)
end

ts_sa, st_sa = let f = motzkin,
    basis = monomial_basis,
    G = PermGroup(perm"(1,2)")

    sa_basis, symmetry_adaptation_time, = @timed let
        ssimple_basis = SymbolicWedderburn.symmetry_adapted_basis(
            Float64,
            G,
            VariablePermutation(),
            exponents.(basis),
            semisimple=true
        )
        sp = sparse.(ssimple_basis)
        droptol!.(sp, 1e-12)
    end

    let m = SOSModel(OPTIMIZER), basis = basis, R = sa_basis
        @variable m t >= 0
        # @objective m Max t

        soses = @variable m [r in R] SOSPoly(FixedPolynomialBasis(r * basis))
        @constraint m g*f - t == sum(soses)

        optimize!(m)
        @info m, termination_status(m) objective_value(m) solve_time(m) symmetry_adaptation_time
        termination_status(m), solve_time(m)
    end
end

@assert ts == ts_sa == SumOfSquares.MOI.OPTIMAL
@assert st_sa < st


text:
function step_script(s::mjSim, robot)
    elapsed_time = round(Int, s.d.time * 1000)  # elapsed time in milliseconds (non-paused simulation)

    # check every 100 milliseconds for another action to take
    if !s.paused && elapsed_time % 100 == 0 && elapsed_time > 0
        #println(elapsed_time, ": ", elapsed_time, "\tframecount: ", round(Int, s.framecount))

        # Start trotting 2 seconds after starting the simulation
        if elapsed_time == 2000 && robot.command.horizontal_velocity[1] > 0
            toggle_trot(robot)
            println("Beginning march with velocity", robot.command.horizontal_velocity)
        end
    end
end