text:
"""
    half_mass_radius(
        pos::AbstractMatrix{<:Number},
        mass::Union{AbstractVector{<:Number},Nothing},
        [r²::AbstractVector{<:Number}];
        rmax::Union{Number,Nothing}=nothing,
    )

Returns the half-mass radius of the particles within radius `rmax`, or of all particles if `rmax`
is `nothing`.
"""
function half_mass_radius(
    pos::AbstractMatrix{<:Number},
    mass::AbstractVector{<:Number},
    r²::AbstractVector{<:Number};
    rmax::Union{Number,Nothing}=nothing,
)
    # check dimensions
    @assert size(pos, 1) == 3
    @assert length(mass) == size(pos, 2) == length(r²)

    if isnothing(rmax)
        mask = (:) # all particles
    else
        mask = r² .≤ rmax^2
    end

    mass_half = @views 1 // 2 * sum(mass[mask])

    # mass sorted by distance to center
    ind = @views sortperm(r²[mask])
    mass_sorted = mass[mask][ind]

    # traverse masses outwardly until half of the total mass is reached
    mass_cumul = zero(eltype(mass))
    i = 1
    @inbounds while mass_cumul < mass_half
        mass_cumul += mass_sorted[i]
        i += 1
    end

    r²_masked = @views r²[mask]
    return 1 // 2 * (sqrt(r²_masked[ind[i - 1]]) + sqrt(r²_masked[ind[i - 2]]))
end

function half_mass_radius(
    pos::AbstractMatrix{<:Number},
    mass::AbstractVector{<:Number};
    rmax::Union{Number,Nothing}=nothing,
)
    r² = r²_sphere(pos)
    half_mass_radius(pos, mass, r²; rmax)
end

"""
    half_mass_radius(pos::AbstractMatrix{<:Number}; rmax::Union{Number,Nothing}=nothing)

Returns the half-mass radius of the particles within radius `rmax` assuming equal masses of all
particles, or of all particles if `rmax` is `nothing`.
"""
function half_mass_radius(pos::AbstractMatrix{<:Number}; rmax::Union{Number,Nothing}=nothing)
    # check dimensions
    @assert size(pos, 1) == 3

    r² = r²_sphere(pos)

    if isnothing(rmax)
        return quantile(sqrt.(r²), 0.5)
    else
        mask = r² .≤ rmax^2
        return @views quantile(sqrt.(r²[mask]), 0.5)
    end
end

half_mass_radius(pos::AbstractMatrix{<:Number}, mass::Nothing; kwargs...) = half_mass_radius(pos, kwargs...)

"""
    half_mass_radius(g::AbstractGalaxy, ptype::Symbol; kwargs...)
    half_mass_radius(p::Particles; kwargs...)

Returns the particles' half-mass radius, depending on if a mass vector is available or not.
If the mass is not available, equal particle masses are assumed.
"""
function half_mass_radius(p::Particles; kwargs...)
    if haskey(p, :mass) && p.mass isa AbstractVector
        half_mass_radius(p.pos, p.mass; kwargs...)
    else
        half_mass_radius(p.pos; kwargs...)
    end
end

half_mass_radius(g::AbstractGalaxy, ptype::Symbol; kwargs...) = half_mass_radius(g[ptype]; kwargs...)

"""
    half_mass_radius_2D(
        pos::AbstractMatrix{<:Number},
        [mass::Union{AbstractVector{<:Number},Nothing}];
        rmax::Union{Number,Nothing}=nothing,
        perspective::Symbol=:edgeon,
    )

Returns the 2D half-mass radius of the particles within radius `rmax` when viewed from a given
`perspective` (`:edgeon`, `:sideon`, `:faceon`), or of all particles if `rmax` is `nothing`.
"""
function half_mass_radius_2D(
    pos::AbstractMatrix{<:Number},
    mass::AbstractVector{<:Number};
    rmax::Union{Number,Nothing}=nothing,
    perspective::Symbol=:edgeon,
)
    # check dimensions
    @assert size(pos, 1) == 3

    dims = get_dims(perspective)
    r² = @views r²_circle(pos[dims, :])

    return half_mass_radius(pos, mass, r²; rmax)
end

function half_mass_radius_2D(
    pos::AbstractMatrix{<:Number};
    rmax::Union{Number,Nothing}=nothing,
    perspective::Symbol=:edgeon,
)
    # check dimensions
    @assert size(pos, 1) == 3

    dims = get_dims(perspective)
    r² = @views r²_circle(pos[dims, :])

    if isnothing(rmax)
        return quantile(sqrt.(r²), 0.5)
    else
        mask = r² .≤ rmax^2
        return @views quantile(sqrt.(r²[mask]), 0.5)
    end
end

function half_mass_radius_2D(pos::AbstractMatrix{<:Number}, mass::Nothing; kwargs...)
    half_mass_radius_2D(pos; kwargs...)
end

"""
    half_mass_radius_2D(g::AbstractGalaxy, ptype::Symbol; kwargs...)
    half_mass_radius_2D(p::Particles; kwargs...)

Returns the particles' 2D half-mass radius, depending on if a mass vector is available or not.
If the mass is not available, equal particle masses are assumed.
"""
function half_mass_radius_2D(p::Particles; kwargs...)
    if haskey(p, :mass) && p.mass isa AbstractVector
        half_mass_radius_2D(p.pos, p.mass; kwargs...)
    else
        half_mass_radius_2D(p.pos; kwargs...)
    end
end

half_mass_radius_2D(g::AbstractGalaxy, ptype::Symbol; kwargs...) = half_mass_radius_2D(g[ptype]; kwargs...)


text:
using SpecialFunctions

abstract type Volume end

#General type of a volume decomposed into spherical shells
abstract type SphericalVolume <: Volume end

"""A volume described by shell-wise spherical harmonics expansions"""
type SphericalHarmonicsVolume <: SphericalVolume
    #spherical harmonics coefficients
    coeff::Array{Array{Complex{Float64}}}

    #bandlimit
    LMAX::Int64

    #radial parameters
    KMAX::Int64
    rmax::Float64

    #radial interpolation
    radial_interp::Bool

    function SphericalHarmonicsVolume(coeff, LMAX::Int64, KMAX::Int64, rmax::Float64, radial_interp::Bool=true)
        new(coeff, LMAX, KMAX, rmax, radial_interp)
    end

    function SphericalHarmonicsVolume(LMAX::Int64, KMAX::Int64, rmax::Float64, radial_interp::Bool=true)
        coeff = [zeros(Complex{Float64}, num_coeff(LMAX)) for k=1:KMAX]
        new(coeff, LMAX, KMAX, rmax, radial_interp)
    end
end

"""A volume described by uniform samples of spherical shells"""
type SurfaceVolume <: SphericalVolume
    #spherical harmonics coefficients
    surf::Array{Array{Complex{Float64}}}

    #bandlimit
    LMAX::Int64

    #radial parameters
    KMAX::Int64
    rmax::Float64

    #radial interpolation
    radial_interp::Bool

    function SurfaceVolume(surf, LMAX::Int64, KMAX::Int64, rmax::Float64, radial_interp::Bool=true)
        new(surf, LMAX, KMAX, rmax, radial_interp)
    end

    function SurfaceVolume(LMAX::Int64, KMAX::Int64, rmax::Float64, radial_interp::Bool=true)
        surf = [zeros(Complex{Float64}, num_val(LMAX)) for k=1:KMAX]
        new(surf, LMAX, KMAX, rmax, radial_interp)
    end
end

#############################################################################
# General Spherical Harmonics functions
#############################################################################

"Wigner 3j symbol from the GSL library"
function wigner(l1::Int64, m1::Int64, l2::Int64, m2::Int64, l3::Int64, m3::Int64)
    return ccall((:gsl_sf_coupling_3j, "libgsl"), Float64, (Int32, Int32, Int32, Int32, Int32, Int32), 2*l1,2*l2,2*l3, 2*m1, 2*m2, 2*m3)
end
#Tested with https://www.wolframalpha.com/input/?i=ThreeJSymbol%5B%7B8,8%7D,+%7B8,-7%7D,+%7B8,-1%7D%5D

"Wrapper for the normalized Legendre polynomials"
function sphPlm(l::Int64,m::Int64,theta::Float64)
    if m>= 0
        return ccall((:gsl_sf_legendre_sphPlm, "libgsl"), Float64, (Int32, Int32, Float64), l,m,cos(theta))
    else
        ma = abs(m)
        return result = (-1.0)^ma * ccall((:gsl_sf_legendre_sphPlm, "libgsl"), Float64, (Int32, Int32, Float64), l,ma,cos(theta))
    end
end

"Wrapper for the assocaited legendre polynomials"
function Plm(l::Int64,m::Int64,theta::Float64)
    if m>= 0
        return ccall((:gsl_sf_legendre_Plm, "libgsl"), Float64, (Int32, Int32, Float64), l,m,cos(theta))
    else
        ma = abs(m)
        return result = (-1.0)^ma*Float64(factorial(BigInt(l-ma))/factorial(BigInt(l+ma))) * ccall((:gsl_sf_legendre_Plm, "libgsl"), Float64, (Int32, Int32, Float64), l,ma,cos(theta))
    end
end

"Spherical harmonics functions"
function Ylm(l::Int64,m::Int64,phi::Float64,theta::Float64)
    return sphPlm(l,m,theta) * exp(1im*m*phi)
end

#############################################################################
# Spherical Harmonics Expansion functions
#############################################################################

function get_size(LMAX::Int64) return 2*LMAX end
function num_val(LMAX::Int64) return (2*LMAX)^2 end
function num_coeff(LMAX::Int64) return LMAX^2 end

function dr(volume::SphericalVolume) return volume.rmax/volume.KMAX end
function dq(volume::SphericalVolume) return volume.rmax/volume.KMAX end
function qmax(volume::SphericalVolume) return pi/dr(volume) end
function qmax(KMAX::Int64, rmax::Float64) return pi/(rmax/KMAX) end

function +(a::SurfaceVolume, b::SurfaceVolume) return SurfaceVolume(a.surf + b.surf, a.LMAX, a.KMAX, a.rmax) end
function -(a::SurfaceVolume, b::SurfaceVolume) return SurfaceVolume(a.surf - b.surf, a.LMAX, a.KMAX, a.rmax) end
function *(a::SurfaceVolume, b::Number) return SurfaceVolume(a.surf * b, a.LMAX, a.KMAX, a.rmax) end
function *(a::Number, b::SurfaceVolume) return SurfaceVolume(b.surf * a, b.LMAX, b.KMAX, b.rmax) end
function /(a::SurfaceVolume, b::Number) return SurfaceVolume(a.surf / b, a.LMAX, a.KMAX, a.rmax) end
function sum(a::SurfaceVolume) return Base.sum(map(sum,a.surf)) end
function sum(f::Union{Function, Type}, a::SurfaceVolume) return Base.sum(f, map((x)->Base.sum(f,x),a.surf)) end
function length(a::SurfaceVolume) return length(a.surf) end
function maximum(a::SurfaceVolume) return Base.maximum(Base.map(maximum, real(a.surf))) end

"""Averages the structure of a range of results for a given run."""
function average(list::Array{SphericalHarmonicsVolume})
    return getSphericalHarmonicsVolume(reduce(+, map(getSurfaceVolume, list)) / length(list))
end

"""Disables the radial interpolation of the spherical harmonics volume"""
function disable_radial_interpolation(volume::SphericalHarmonicsVolume)
    volume.radial_interp = false
end

"""calculate the coefficients from surface"""
function transform(surface::Vector{Complex{Float64}}, LMAX::Int64)
    rcoeff = zeros(Float64, num_coeff(LMAX))
    icoeff = zeros(Float64, num_coeff(LMAX))

    ccall( (:transform_sample, "libsh.so"), Ptr{Void},
    (Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64}, Int32, Int32),
    real(surface), imag(surface), rcoeff, icoeff, LMAX, 1)

    return complex.(rcoeff,icoeff)
end

"""Calculates surface from coefficients"""
function backtransform(coeff::Vector{Complex{Float64}}, LMAX::Int64)
    rdata = zeros(Float64, num_val(LMAX))
    idata = zeros(Float64, num_val(LMAX))

    ccall((:transform_sample, "libsh.so"), Ptr{Void},
    (Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64,}, Int32, Int32),
    real(coeff), imag(coeff), rdata, idata, LMAX, 0)

    return complex.(rdata, idata)
end

"""SphericalHarmonicsVolume to SurfaceVolume"""
function getSurfaceVolume(volume::SphericalHarmonicsVolume)
    surf = map((x) -> backtransform(x, volume.LMAX), volume.coeff)
    return SurfaceVolume(surf, volume.LMAX, volume.KMAX, volume.rmax, volume.radial_interp)
end

"""SurfaceVolume to SphericalHarmonicsVolume"""
function getSphericalHarmonicsVolume(volume::SurfaceVolume)
    coeff = map((x) -> transform(x, volume.LMAX), volume.surf)
    return SphericalHarmonicsVolume(coeff, volume.LMAX, volume.KMAX, volume.rmax, volume.radial_interp)
end

"""Helper function to get index in linear-array for (l,m)-doublet"""
function seanindex(l::Int64,m::Int64,LMAX::Int64)
    bigL = LMAX - 1
    if m >= 0
        return round(Int64, ( m * ( bigL + 1 ) - ( ( m * (m - 1) ) /2 ) + ( l - m ) ) +1 )
    else
        return round(Int64,  ( ( ( bigL * ( bigL + 3 ) ) /2 ) + 1 + ( ( bigL + m ) * ( bigL + m + 1 ) / 2 ) + ( l - abs( m ) ) ) + 1 )
    end
end

"""Helper function: gets specific spherical harmonics coefficient!"""
function getc(volume::SphericalHarmonicsVolume,k::Int64,l::Int64,m::Int64)
    return volume.coeff[k][seanindex(l,m,volume.LMAX)]
end
"""Helper function: sets specific spherical harmonics coefficient"""
function setc(volume::SphericalHarmonicsVolume,k::Int64,l::Int64,m::Int64,val::Complex{Float64})
    volume.coeff[k][seanindex(l,m,volume.LMAX)] = val
end

"""Helper function: sets specific spherical harmonics coefficient"""
function setc(volume::SphericalHarmonicsVolume,k::Int64,l::Int64,m::Int64,val::Float64)
    setc(volume, k, l, m, Complex(val))
end

"""Gets the whole 2*l+1 dimensional coefficients vector"""
function cvec_get(volume::SphericalHarmonicsVolume,k,l)
    v = Array{Complex{Float64}}(2*l+1)
    for m=-l:l
        v[m+l+1] = getc(volume,k,l,m)
    end
    return v
end
"""Sets the whole 2*l+1 dimensional coefficient vector"""
function cvec_set(volume::SphericalHarmonicsVolume,k, l, v)
    for m=-l:l
        setc(volume,k, l, m, v[m+l+1])
    end
end

"""helper function that translates an index to two angles, according to the size
returns PHI, THETA"""
function i_to_ang(index, LMAX)
    size = get_size(LMAX)
    iphi = mod(index-1,size)
    itheta = floor((index-1)/size)
    phi = iphi / (size) * 2.0*pi
    theta = pi*(2*itheta + 1)/(2*size)
    return (phi,theta)
end

"""Calculates index in linear surface array for given angles"""
function ang_to_i(phi::Float64,theta::Float64, LMAX)
    size = get_size(LMAX)
    phi = mod(phi, 2.0*pi)
    theta = mod(theta, pi)
    ip = round(Int64, phi/(2.0*pi)*size)
    it = round(Int64, ((theta*2*size)/pi - 1)/ 2)
    i = clamp(it*size+ip + 1,1,size^2)
    return Int64(i)
end

"""From a surface function, calculate the spherical harmonics coefficients up to LMAX"""
function getSurface(func, LMAX::Int64)
    surface = Array{Complex{Float64}}(num_val(LMAX))
    for i = 1:num_val(LMAX)
        phi,theta = i_to_ang(i,LMAX)
        surface[i] = func(phi,theta)
    end
    return surface
end

"""Calculates the SH coefficients for a density"""
function getDensity(atomlist, LMAX::Int64, KMAX::Int64, rmax::Float64)
    surflist = Array{Complex{Float64},1}[]
    dr = rmax/KMAX
    for k = 1:KMAX
        r = k*dr
        push!(surflist, getSurface(densityShell(atomlist, k*dr), LMAX))
    end
    surfVolume = SurfaceVolume(surflist, LMAX, KMAX, rmax)

    return getSphericalHarmonicsVolume(surfVolume)
end

"""wraps the density as a function only of the angle with fixed radius"""
function densityShell(atomlist, r)
    return function(phi,theta) return getDensity_cartesian(atomlist, sphericalToCartesian(r,phi,theta)) end
end

"""Calculates the Fourier shell"""
function fourierShell(atomlist, k)
    return function(phi,theta) return getFourier_cartesian(atomlist,sphericalToCartesian(k,phi,theta)) end
end
"""Gets intensity shell"""
function intensityShell(atomlist, k)
    return function(phi,theta) return getIntensity_cartesian(atomlist,sphericalToCartesian(k,phi,theta)) end
end

# """Calculates the spherical harmonics coefficients from a cubic representation"""
# function cubic_to_spherical(cube)
#     cubic_shell = function(r)
#       return function(phi, theta)
#         x = clamp(round(Int64,sphericalToCartesian(r,phi,theta)) + round(Int64,[(cubesize-1)/2 + 1, (cubesize-1)/2 + 1, (cubesize-1)/2 + 1] ), 1, cubesize)
#         return cube[x[1],x[2],x[3]]
#       end
#     end
#     coeff = []
#     for i = 1:KMAX
#         r = i
#         c = get_coefficients(cubic_shell(r))
#         push!(coeff,c)
#     end
#     return coeff
# end

"""Deletes upper terms of a spherical harmonics representation
deleteTerms(coeff,kuct,L, LMAX)"""
function deleteTerms(volume::SphericalHarmonicsVolume, K::Int64, L::Int64)
    newvol = deepcopy(volume)
    for k = 1:newvol.KMAX
        for l = 0:newvol.LMAX-1
            for m = -l:l
                if k > K || l > L
                    setc(newvol,k,l,m,0.0)
                end
            end
        end
    end
    return newvol
end

# Hankel transform

"""Numerical integration of \$f \\cdot r dr\$ """
function frdr(f::Vector{Float64}, r::Vector{Float64})
    l = length(r)
    res = Array{Float64}(l)

    @inbounds res[1] = sqrt(0.5) * (r[1] + r[2])^2
    for j = 2:l-1
        @inbounds res[j] = (r[j+1] - r[j-1]) * (r[j-1] + 2*r[j] + r[j+1])
    end
    @inbounds res[l] = 4*r[l]^2 - (r[l-1] + r[l])^2

    return pi/4 * (res.^(1.5) .*f)
end

"""Calculates the Hankel transform
Source:
# *http://www.mathworks.de/matlabcentral/fileexchange/13371-hankeltransform
# *http://documents.epfl.ch/users/l/le/leuteneg/www/MATLABToolbox/HankelTransform.html"""
function hat!(res::Vector{Float64}, h::Vector{Float64}, I::Matrix{Float64}, forward::Bool=true)
    l = length(h)
    if forward
        r = collect(1.0:l)
        res[:] = I*frdr(h, r)
    else
        k = pi/l * collect(1.0:l)
        res[:] = transpose(I)*frdr(h, k)/(2.0*pi)^2
    end
end

"""Precalcualtes Imatrix for numerical integration of Hankel transform"""
function precalcImatrices(LMAX::Int64, KMAX::Int64)

    Imatrix = function(len, n)
        r = Float64[1:len;]
        k = pi/len * Float64[1:len;]
        return SpecialFunctions.besselj.(n, k*transpose(r))
    end

    #Calculate the matrices, if not already available
    if !isdefined(:Ilist) || length(Ilist) != LMAX || size(Ilist[1])[1] != KMAX
        global Ilist = []
        for l = 0:LMAX-1
            push!(Ilist, Imatrix(KMAX, l+0.5))
        end
    end

end

"""On a set of spherical harmonics coefficients, this performs a hankel transform with a given max value of rho (rmax or qmax)"""
function sphBesselTransform(volume::SphericalHarmonicsVolume, forward::Bool)
    newvol = deepcopy(volume)

    #Precalculate the basis
    precalcImatrices(volume.LMAX, volume.KMAX)
    samp_r, samp_c = zeros(Float64, volume.KMAX), zeros(Float64, volume.KMAX)

    for l = 0:volume.LMAX-1
        ipart = Ilist[l+1]
        for m = -l:l
            rho = Complex{Float64}[getc(volume,k,l,m) for k=1:volume.KMAX]
            hat!(samp_r, real(rho), ipart, forward)
            hat!(samp_c, imag(rho), ipart, forward)

            for k=1:volume.KMAX
                fac = forward ? Complex{Float64}((1im)^l *sqrt(1/k)) : Complex{Float64}((-1im)^l *sqrt(1/k))
                setc(newvol, k, l, m, fac*complex(samp_r[k],samp_c[k]))
            end
        end
    end
    newvol.rmax = qmax(volume) #this also gives rmax for backtransform
    return newvol
end

"""Takes density SH and performs hankel transform to get amplitude SH"""
function forward(volume::SphericalHarmonicsVolume)
    return sphBesselTransform(volume, true)
end

"""Takes density SH and performs hankel transform to get amplitude SH"""
function backward(volume::SphericalHarmonicsVolume)
    return sphBesselTransform(volume, false)
end

"Converts a density surface array into its Fourier counterpart"
function forward(volume::SurfaceVolume)
    volume = getSphericalHarmonicsVolume(volume)
    volume = forward(volume)
    return getSurfaceVolume(volume)
end

"Converts a fourier surface array into its real space counterpart"
function backward(volume::SurfaceVolume)
    volume = getSphericalHarmonicsVolume(volume)
    volume = backward(volume)
    return getSurfaceVolume(volume)
end

"""Indirectly calculate intensity coefficients by squaring on surface"""
function absoluteSquare(volume::SphericalHarmonicsVolume)
    newvol = deepcopy(volume)
    newvol.coeff = map((x) -> transform( complex(abs.(backtransform(x, volume.LMAX)).^2), volume.LMAX), volume.coeff)
    return newvol
end

"Calculates the Intensity Shell Correlation for a set of coefficients"
function shell_correlation(volume1::SphericalHarmonicsVolume, volume2::SphericalHarmonicsVolume, K_range::UnitRange{Int64}=1:0, cor_method=cor_FSC)
    surf1 = getSurfaceVolume(volume1)
    surf2 = getSurfaceVolume(volume2)
    @assert length(surf1) == length(surf2)
    K_range = (length(K_range) == 0 ? (1:length(surf1)) : K_range)

    return Float64[cor_method(surf1.surf[k], surf2.surf[k]) for k in K_range]
end

function shell_correlation_ISC(volume1::SphericalHarmonicsVolume, volume2::SphericalHarmonicsVolume, K_range::UnitRange{Int64}=1:0)
    shell_correlation(volume1, volume2, K_range, cor_ISC)
end

function shell_correlation_FSC(volume1::SphericalHarmonicsVolume, volume2::SphericalHarmonicsVolume, K_range::UnitRange{Int64}=1:0)
    shell_correlation(volume1, volume2, K_range, cor_FSC)
end

function ISC(density::SphericalHarmonicsVolume, intensity::SphericalHarmonicsVolume, K_range::UnitRange{Int64}=1:0)
    return shell_correlation_ISC(intensity, absoluteSquare(forward(density)), K_range)
end

function FSC(density::SphericalHarmonicsVolume, fourier::SphericalHarmonicsVolume, K_range::UnitRange{Int64}=1:0)
    return shell_correlation_FSC(fourier, forward(density), K_range)
end

"Calculates the similarity between two sets of coefficients via Fourier Shell Correlation"
function similarity(volume1::SphericalHarmonicsVolume, volume2::SphericalHarmonicsVolume, K::Int64)
    @assert K>0
    return sum(shell_correlation_ISC(volume1, volume2, 1:K))/K
end

"Calculates the similarity between two sets of coefficients via Fourier Shell Correlation"
function similarity(volume1::SphericalHarmonicsVolume, volume2::SphericalHarmonicsVolume, K_range::UnitRange{Int64})
    @assert length(K_range)>0
    return sum(shell_correlation_ISC(volume1, volume2, K_range))/length(K_range)
end

"Calculates the similarity between two sets of coefficients via Fourier Shell Correlation"
function similarity(volume1::SphericalHarmonicsVolume, volume2::SphericalHarmonicsVolume)
    return similarity(volume1, volume2, volume1.KMAX)
end

#Complex to real and real to complex

"""Calculates the U matrix for order l that translates complex to real coefficients"""
function Umat(l)
    l = float(l)
    Gamma = (m)-> m>0 ? 1 : 0

    Delta = (a,b)-> a==b ? 1 : 0

    U = function(u,m)
        return ( Delta(m,0)*Delta(u,0) + 1/sqrt(2)*( Gamma(u)*Delta(m,u) + Gamma(-u)*1im*(-1)^m*Delta(m,u)+
        Gamma(-u)*(-1im)*Delta(m,-u)+Gamma(u)*(-1)^m*Delta(m,-u)))
    end

    return Complex{Float64}[ U(x,y) for x = -l:l, y= -l:l]
end

"""input array of [k,l] entries with [-m ... m] spherical harmonics coefficients in a vector"""
function comp_to_real(volume::SphericalHarmonicsVolume)
    newvol = deepcopy(volume)
    for l = 1:volume.LMAX-1
        mat = Umat(l)
        for k = 1:volume.KMAX
            cvec_set(newvol, k,l, mat*cvec_get(volume, k, l))
        end
    end
    return newvol
end

"""Calculates real spherical harmonics coefficients from complex"""
function real_to_comp(volume::SphericalHarmonicsVolume)
    newvol = deepcopy(volume)
    for l = 1:volume.LMAX-1
        mat = ctranspose(Umat(l))
        for k = 1:volume.KMAX
            cvec_set(newvol,k,l, mat*cvec_get(volume, k, l))
        end
    end
    return newvol
end

"""Given a set of surfaces, calculates the cartesian value via interpolation"""
function getVolumeInterpolated(volume::SurfaceVolume, vec::Vector{Float64})
    q,p,t = cartesianToSpherical(vec)
    size = get_size(volume.LMAX)
    dp = 2.0*pi / size
    dt = pi / size
    surf = volume.surf
    k = q/dr(volume)

    k0 = Int64(clamp(round(Int64, round(k)), 1, volume.KMAX))
    k1 = k0

    if volume.radial_interp
        k0 = Int64(clamp(round(Int64, floor(k)), 1, volume.KMAX))
        k1 = Int64(clamp(round(Int64, ceil(k)), 1, volume.KMAX))
    end

    #phi_k = 2*pi*k/2B,
    p0 = floor(p/dp)*dp
    p1 = ceil(p/dp)*dp

    #theta_j = (2j + 1)*pi/4B
    # t0 = clamp(floor(t/dt-0.5)*dt,0,pi)
    # t1 = clamp(ceil(t/dt-0.5)*dt,0,pi)
    t0 = max(floor( ((t*2*size)/pi - 1)/ 2) *dt + dt/2, 0)
    t1 = min(ceil( ((t*2*size)/pi - 1)/ 2) *dt+ dt/2, pi)

    a = k-k0
    b = (p-p0)/dp
    c = (t-t0)/dt

    V000 = (1-a) * (1-b) * (1-c)
    V100 = a *(1-b) *(1-c)
    V010 = (1-a)* b* (1-c)
    V001 = (1-a)* (1-b)* c
    V101 = a*(1-b)* c
    V011 = (1-a)* b* c
    V110 = a* b* (1-c)
    V111 = a* b* c

    return      V000 * surf[k0][ang_to_i(p0,t0,volume.LMAX)] +
    V100 * surf[k1][ang_to_i(p0,t0,volume.LMAX)] +
    V010 * surf[k0][ang_to_i(p1,t0,volume.LMAX)] +
    V001 * surf[k0][ang_to_i(p0,t1,volume.LMAX)] +
    V101 * surf[k1][ang_to_i(p0,t1,volume.LMAX)] +
    V011 * surf[k0][ang_to_i(p1,t1,volume.LMAX)] +
    V110 * surf[k1][ang_to_i(p1,t0,volume.LMAX)] +
    V111 * surf[k1][ang_to_i(p1,t1,volume.LMAX)]
end

"""Calculates a cube from spherical shells"""
function getCube(volume::SurfaceVolume)
    cubesize = 2*volume.KMAX+1
    getCube(volume, cubesize)
end

"""Calculates a cube from spherical harmonics coefficients"""
function getCube(volume::SphericalHarmonicsVolume)
    return getCube(getSurfaceVolume(volume))
end

"""Calculates a cube from spherical shells"""
function getCube(volume::SurfaceVolume, cubesize::Int64)
    r = linspace(-volume.rmax, volume.rmax, cubesize)
    return CubeVolume(Complex{Float64}[getVolumeInterpolated(volume,[x,y,z]) for x=r, y=r, z=r], cubesize, volume.rmax)
end

"""Calculates a cube from spherical harmonics coefficients"""
function getCube(volume::SphericalHarmonicsVolume, cubesize::Int64)
    return getCube(getSurfaceVolume(volume), cubesize)
end

"""Saving the cube calculated from a spherical harmonics volume"""
function saveCube(volume::SphericalHarmonicsVolume, filename)
    saveCube(getCube(volume), filename)
end

##########################################################
# rotations
##########################################################

"""Big number factorial function"""
function myfactorial(x)
    return factorial(BigInt(x))
end

"""working calculation of the matrix d elements taken from:
Source: unknown, maybe Ludger?"""
function d_work(l::Int64, m1::Int64, m2::Int64, beta::Float64)
    t1 = BigInt(l+m1)
    t2 = BigInt(l-m2)
    t3 = BigInt(m2 - m1)
    t4 = BigInt(2*l + m1 - m2)
    vmin = (-t3 < 0) ? 0 : -t3
    vmax = (t2 < t1) ? t2 : t1

    factor = (myfactorial(l+m1)*myfactorial(l-m1)*myfactorial(l+m2)*myfactorial(l-m2))^(1/2)
    val = 0

    for s = vmin:vmax
        val += cos(beta/2)^(2*l+m1-m2-2*s) * sin(beta/2)^(m2-m1+s*2) * (-1.0)^(s+m2-m1)/(myfactorial(l+m1-s)*myfactorial(s)*myfactorial(m2-m1+s)*myfactorial(l-m2-s))
    end
    return Complex{Float64}(val*factor)
end

function Dlms(l::Int64, m::Int64, s::Int64, theta::Float64, phi::Float64, gamma::Float64)
    return exp(-1im*s*gamma-1im*m*phi) * d_work(l, m, s, theta)
end

function rotationMatrix(theta::Float64, phi::Float64, gamma::Float64, l::Int64)
    return Complex{Float64}[ Dlms(l, m, s, theta, phi, gamma) for m=-l:l,s=-l:l]
end

"""Rotates a volume via rotating the coefficients"""
function rotateStructure(volume::SphericalHarmonicsVolume, theta::Float64, phi::Float64, gamma::Float64, K::Int64, lrange::Range)
    newvol = deepcopy(volume)

    for l in lrange
        rot = rotationMatrix(theta, phi, gamma, l)
        for k = 1:K
            cvec_set(newvol,k,l, rot*cvec_get(volume, k,l) )
        end
    end
    return newvol
end
"""Overloaded version of the rotateStruture function"""
function rotateStructure(volume::SphericalHarmonicsVolume, theta::Float64, phi::Float64, gamma::Float64, K::Int64, L::Int64)
    return rotateStructure(volume, theta, phi, gamma, K, 1:L)
end

"""Calculates the function at k,phi,theta for a set of spherical harmonic coefficient"""
function get_value(volume::SphericalHarmonicsVolume,k, phi, theta)
    val = 0
    for l = 0:volume.KMAX-1
        for m=-l:l
            val += getc(volume,k,l,m)*(Ylm(l,m,phi, theta))
        end
    end
    return val
end

function negativityCheck(volume::SphericalHarmonicsVolume, range::UnitRange{Int64})
    volume = getSurfaceVolume(volume)
    mapreduce((x) -> max.(sum(abs, min.(real(x),0))/sum(abs,x),0), +, volume.surf[range])
end

ext: jl


text:
struct Simulation{ET <: Equation, IT <: Integrator, ST <: Solution}
    equation::ET
    integrator::IT
    solution::ST
    ncycle::Int
    run_id::String
    filename::String
end

function Simulation(equ::ET, int::IT, sol::ST, run_id::String, filename::String; ncycle=1) where {ET,IT,ST}
    Simulation{ET,IT,TT}(equ, int, sol, ncycle, run_id, filename)
end

function Simulation(equ::Equation, int::Integrator, Δt, run_id, filename, ntime, nsave; ncycle=1)
    Simulation(equ, int, Solution(equ, Δt, ntime, nsave), run_id, filename; ncycle=ncycle)
end

function Simulation(equ::Equation, tableau::AbstractTableau, Δt, run_id, filename, ntime, nsave; ncycle=1)
    Simulation(equ, Integrator(equ, tableau, Δt), Δt, run_id, filename, ntime, nsave; ncycle=ncycle)
end

function Simulation(equ::Equation, integrator, tableau::AbstractTableau, Δt, run_id, filename, ntime, nsave; ncycle=1)
    Simulation(equ, integrator(equ, tableau, Δt), Δt, run_id, filename, ntime, nsave; ncycle=ncycle)
end

function run!(sim::Simulation)

    println("Running ", sim.run_id, "...")

    h5 = create_hdf5(sim.solution, sim.filename)

    try
        @showprogress 5 for c in 1:sim.ncycle
            integrate!(sim.integrator, sim.solution)
            write_to_hdf5(sim.solution, h5)
            reset!(sim.solution)
        end
    catch ex
        if isa(ex, DomainError)
            warn("DOMAIN ERROR")
        elseif isa(ex, ErrorException)
            warn("Simulation exited early.")
            warn(ex.msg)
        else
            throw(ex)
        end
    end

    close(h5)

    return sim.solution
end


text:
@inline function parent_cell(p::NTuple{N,A}, dxi::NTuple{N,B}) where {N,A,B}
    return ntuple(i -> Int16(p[i] ÷ dxi[i]) + 1, Val(N))
end

# dimension-agnostic fully unrolled euclidean distance
@generated function distance(a::NTuple{N,T}, b::NTuple{N,T}) where {N,T}
    ex = zero(T)
    @inbounds for i in 1:N
        ex = :((a[$i] - b[$i]) * (a[$i] - b[$i]) + $ex)
    end
    return :(√($ex))
end

# check whether particle is inside the grid (includes boundary)
@inline function isinside(px::Real, py::Real, x, y)
    xmin, xmax = extrema(x)
    ymin, ymax = extrema(y)
    @assert (xmin ≤ px ≤ xmax) && (ymin ≤ py ≤ ymax)
end

@inline function isinside(px::Real, py::Real, pz::Real, x, y, z)
    xmin, xmax = extrema(x)
    ymin, ymax = extrema(y)
    zmin, zmax = extrema(z)
    @assert (xmin ≤ px ≤ xmax) && (ymin ≤ py ≤ ymax) && (zmin ≤ pz ≤ zmax)
end

@inline isinside(p::NTuple{2,T1}, x::NTuple{2,T2}) where {T1,T2} = isinside(p[1], p[2], x[1], x[2])

@inline function isinside(p::NTuple{3,T1}, x::NTuple{3,T2}) where {T1,T2}
    return isinside(p[1], p[2], p[3], x[1], x[2], x[3])
end

# normalize coordinates
@inline function normalize_coordinates(
    p::NTuple{N,A}, xi::NTuple{N,B}, dxi::NTuple{N,C}, idx::NTuple{N,D}
) where {N,A,B,C,D}
    return ntuple(i -> (p[i] - xi[i][idx[i]]) * (1 / dxi[i]), Val(N))
end

# compute grid size
function grid_size(x::NTuple{N,T}) where {T,N}
    return ntuple(i -> x[i][2] - x[i][1], Val(N))
end

# Get field F at the corners of a given cell
@inline function field_corners(F::AbstractArray{T,2}, idx::NTuple{2,Integer}) where {T}
    idx_x, idx_y = idx
    return (
        F[idx_x, idx_y], F[idx_x + 1, idx_y], F[idx_x, idx_y + 1], F[idx_x + 1, idx_y + 1]
    )
end

@inline function field_corners(F::AbstractArray{T,3}, idx::NTuple{3,Integer}) where {T}
    idx_x, idx_y, idx_z = idx
    return (
        F[idx_x, idx_y, idx_z],   # v000
        F[idx_x + 1, idx_y, idx_z],   # v100
        F[idx_x, idx_y, idx_z + 1], # v001
        F[idx_x + 1, idx_y, idx_z + 1], # v101
        F[idx_x, idx_y + 1, idx_z],   # v010
        F[idx_x + 1, idx_y + 1, idx_z],   # v110
        F[idx_x, idx_y + 1, idx_z + 1], # v011
        F[idx_x + 1, idx_y + 1, idx_z + 1], # v111
    )
end

@inline particle2tuple(p::NTuple{N,AbstractArray}, ix) where {N} = ntuple(i -> p[i][ix], Val(N))

# 2D random particle generator for regular grids
function random_particles(nxcell, x, y, dx, dy, nx, ny)
    # number of cells
    ncells = (nx - 1) * (ny - 1)
    # allocate particle coordinate arrays
    px, py = zeros(nxcell * ncells), zeros(nxcell * ncells)
    Threads.@threads for i in 1:(nx - 1)
        @inbounds for j in 1:(ny - 1)
            # lowermost-left corner of the cell
            x0, y0 = x[i], y[j]
            # cell index
            cell = i + (nx - 1) * (j - 1)
            for l in 1:nxcell
                px[(cell - 1) * nxcell + l] = rand() * dx + x0
                py[(cell - 1) * nxcell + l] = rand() * dy + y0
            end
        end
    end

    return px, py
end

# 3D random particle generator for regular grids
function random_particles(nxcell, x, y, z, dx, dy, dz, nx, ny, nz)
    # number of cells
    ncells = (nx - 1) * (ny - 1) * (nz - 1)
    # allocate particle coordinate arrays
    px, py, pz = zeros(nxcell * ncells), zeros(nxcell * ncells), zeros(nxcell * ncells)
    Threads.@threads for i in 1:(nx - 1)
        @inbounds for j in 1:(ny - 1), k in 1:(nz - 1)
            # lowermost-left corner of the cell
            x0, y0, z0 = x[i], y[j], z[k]
            # cell index
            cell = i + (nx - 1) * (j - 1) + (nx - 1) * (ny - 1) * (k - 1)
            for l in 1:nxcell
                px[(cell - 1) * nxcell + l] = rand() * dx + x0
                py[(cell - 1) * nxcell + l] = rand() * dy + y0
                pz[(cell - 1) * nxcell + l] = rand() * dz + z0
            end
        end
    end

    return px, py, pz
end


text:
PyArg_CheckNumArgsEq(name::String, args::PyPtr, n::Integer) = begin
    m = isnull(args) ? 0 : PyTuple_Size(args)
    if m == n
        return 0
    else
        PyErr_SetString(PyExc_TypeError(), "$name() takes $n arguments ($m given)")
        return -1
    end
end

PyArg_CheckNumArgsGe(name::String, args::PyPtr, n::Integer) = begin
    m = isnull(args) ? 0 : PyTuple_Size(args)
    if m ≥ n
        return 0
    else
        PyErr_SetString(PyExc_TypeError(), "$name() takes at least $n arguments ($m given)")
        return -1
    end
end

PyArg_CheckNumArgsLe(name::String, args::PyPtr, n::Integer) = begin
    m = isnull(args) ? 0 : PyTuple_Size(args)
    if m ≤ n
        return 0
    else
        PyErr_SetString(PyExc_TypeError(), "$name() takes at most $n arguments ($m given)")
        return -1
    end
end

PyArg_CheckNumArgsBetween(name::String, args::PyPtr, n0::Integer, n1::Integer) = begin
    m = isnull(args) ? 0 : PyTuple_Size(args)
    if n0 ≤ m ≤ n1
        return 0
    else
        PyErr_SetString(PyExc_TypeError(), "$name() takes $n0 to $n1 arguments ($m given)")
        return -1
    end
end

PyArg_CheckNoKwargs(name::String, kwargs::PyPtr) = begin
    if isnull(kwargs) || PyObject_Length(kwargs) == 0
        return 0
    else
        argnames = PyIterable_Collect(kwargs, String)
        isempty(argnames) && PyErr_IsSet() && return -1
        PyErr_SetString(
            PyExc_TypeError(),
            "$name() got unexpected keyword arguments: $(join(["'$n'" for n in argnames], ", "))",
        )
        return -1
    end
end

struct NODEFAULT end

PyArg_Find(args::PyPtr, kwargs::PyPtr, i::Union{Int,Nothing}, k::Union{String,Nothing}) =
    if i !== nothing && !isnull(args) && 0 ≤ i < PyTuple_Size(args)
        return PyTuple_GetItem(args, i)
    elseif k !== nothing &&
           !isnull(kwargs) &&
           (ro = PyDict_GetItemString(kwargs, k)) != PyNULL
        return ro
    else
        return PyNULL
    end

"""
    PyArg_GetArg(T, funcname, [args, i], [kwargs, k], [default])

Attempt to find and convert the specified argument to a `T`. Return 0 on success, in which case the result can be retrieved with `takeresult(T)`. Return -1 on failure, with a Python error set.

- `funcname::String` is the name of the function this is used in, for constructing error messages.
- `args::PyPtr` is a tuple of arguments, and `i::Int` an index.
- `kwargs::PyPtr` is a dict of keyword arguments, and `k::String` a key.
- `default` specifies a default value if the argument is not found. NOTE: It need not be a `T`, so when retrieving the result use `takeresult(Union{T,typeof(default)})`.
"""
PyArg_GetArg(
    ::Type{T},
    name::String,
    args::PyPtr,
    i::Union{Int,Nothing},
    kwargs::PyPtr,
    k::Union{String,Nothing},
    d = NODEFAULT(),
) where {T} = begin
    ro = PyArg_Find(args, kwargs, i, k)
    if isnull(ro)
        if d !== NODEFAULT()
            putresult(d)
            return 0
        elseif k !== nothing
            PyErr_SetString(PyExc_TypeError(), "$name() did not get required argument '$k'")
            return -1
        elseif i !== nothing && i ≥ 0
            PyErr_SetString(
                PyExc_TypeError(),
                "$name() takes at least $(i+1) arguments (got $(isnull(args) ? 0 : PyTuple_Size(args)))",
            )
            return -1
        else
            error("impossible to satisfy this argument")
        end
    end
    r = PyObject_TryConvert(ro, T)
    if r == -1
        return -1
    elseif r == 0
        PyErr_SetString(
            PyExc_TypeError(),
            "Argument $(k !== nothing ? "'$k'" : i !== nothing ? "$i" : error("impossible")) to $name() must be convertible to a Julia '$T'",
        )
        return -1
    else
        return 0
    end
end
PyArg_GetArg(::Type{T}, name::String, args::PyPtr, i::Int, d = NODEFAULT()) where {T} =
    PyArg_GetArg(T, name, args, i, PyNULL, nothing, d)
PyArg_GetArg(::Type{T}, name::String, kwargs::PyPtr, k::String, d = NODEFAULT()) where {T} =
    PyArg_GetArg(T, name, PyNULL, nothing, kwargs, k, d)


text:
module Lorenz63
export lorenz_63

using Parameters: @with_kw, @unpack
using ..ExampleBase: LEDemo, ContinuousExample

@with_kw struct Lorenz63Param
    σ::Float64 = 10.0
    ρ::Float64 = 28.0
    β::Float64 = 8/3
end

function phase_dynamics!(du, u, p, t)
    @unpack σ, ρ, β = p
    du[1] = σ * (u[2] - u[1])
    du[2] = u[1] * (ρ - u[3]) - u[2]
    du[3] = u[1] * u[2] - β*u[3]
end

@inline @views function tangent_dynamics!(du, u, p, t)
    @unpack σ, ρ, β = p
    phase_dynamics!(du[:, 1], u[:, 1], p, t)
    du[1, 2:end] .= σ .* (u[2, 2:end] .- u[1, 2:end])
    du[2, 2:end] .=
        u[1, 2:end] .* (ρ .- u[3, 1]) .-
        u[1, 1] .* u[3, 2:end] .-
        u[2, 2:end]
    du[3, 2:end] .=
        u[1, 2:end] .* u[2, 1] .+ u[1, 1] .* u[2, 2:end] .-
        β .* u[3, 2:end]
end

"""
Return a [`LEDemo`](@ref) for the Lorenz system.

* <https://en.wikipedia.org/wiki/Lorenz_system>
* <http://sprott.physics.wisc.edu/chaos/comchaos.htm>
* E. N. Lorenz, J. Atmos. Sci. 20, 130-141 (1963)
"""
function lorenz_63(;
        u0=[0.1, 0.1, 0.1],
        t_renorm=1.0,
        t_attr=100000,
        atol=1e-2, rtol=1e-2,
        terminator_options = [:atol => atol, :rtol => rtol],
        kwargs...)
    LEDemo(ContinuousExample(
        "Lorenz (1963)",
        phase_dynamics!, u0, t_renorm, Lorenz63Param(),
        tangent_dynamics!,
        t_attr,
        [0.9056, 0, -14.5723],  # known_exponents
        atol, rtol,
        terminator_options,
    ); kwargs...)
end

end


text:
n = "Kevin"
a = 36


text:
# helper convertion functions
Base.convert(::Type{ImVec2}, x::NTuple{2}) = ImVec2(x...)
Base.convert(::Type{ImVec2}, x::Vector) = ImVec2(x...)
Base.convert(::Type{ImVec4}, x::NTuple{4}) = ImVec4(x...)
Base.convert(::Type{ImVec4}, x::Vector) = ImVec4(x...)
Base.convert(::Type{ImVec4}, x::ImU32) = ColorConvertU32ToFloat4(x)
Base.convert(::Type{ImU32}, x::ImVec4) = ColorConvertFloat4ToU32(x)

# TODO: put these in CEnum.jl
Base.:~(x::Cenum{UInt32}) = ~UInt32(x)
Base.:(:)(a::T, b::Cenum) where {T<:Integer} = a:T(b)
Base.:(:)(a::Cenum, b::T) where {T<:Integer} = T(a):b

function ShowFlags(::Type{T}) where {T<:Cenum}
    io = IOBuffer()
    s = "```\n"
    for (n,v) in name_value_pairs(T)
        s *= string(n)*" $v \n"
    end
    s*"\n```"
end
GetFlags(::Type{T}) where {T<:Cenum} = name_value_pairs(T) |> collect

# simple unsafe destruction helper
function UnsafeGetPtr(x::Ptr{T}, name::Symbol) where {T}
    offset = x
    type = T
    flag = false
    for i = 1:fieldcount(T)
        name == fieldname(T, i) || continue
        flag = true
        type = fieldtype(T, i)
        offset += fieldoffset(T, i)
        break
    end
    flag || throw(ArgumentError("$T has no field named $name."))
    return Ptr{type}(offset)
end

function Get(x::Ptr{T}, name::Symbol) where {T}
    GC.@preserve x begin
        value = unsafe_load(UnsafeGetPtr(x, name))
    end
    return value
end

function Set(x::Ptr{T}, name::Symbol, value::S) where {T,S}
    GC.@preserve x value begin
        unsafe_store!(UnsafeGetPtr(x, name), value)
    end
    return value
end

# emulate ImGui::GetIO().xxx
function Base.getproperty(io::Ptr{ImGuiIO}, x::Symbol)
    x == :ConfigFlags && return ImGuiIO_Get_ConfigFlags(io)
    x == :BackendFlags && return ImGuiIO_Get_BackendFlags(io)
    x == :DisplaySize && return ImGuiIO_Get_DisplaySize(io)
    x == :DeltaTime && return ImGuiIO_Get_DeltaTime(io)
    x == :IniSavingRate && return ImGuiIO_Get_IniSavingRate(io)
    x == :IniFilename && return ImGuiIO_Get_IniFilename(io)
    x == :LogFilename && return ImGuiIO_Get_LogFilename(io)
    x == :MouseDoubleClickTime && return ImGuiIO_MouseDoubleClickTime(io)
    x == :MouseDoubleClickMaxDist && return ImGuiIO_Get_MouseDoubleClickMaxDist(io)
    x == :MouseDragThreshold && return ImGuiIO_MouseDragThreshold(io)
    x == :KeyRepeatDelay && return ImGuiIO_Get_KeyRepeatDelay(io)
    x == :KeyRepeatRate && return ImGuiIO_Get_KeyRepeatRate(io)
    x == :UserData && return ImGuiIO_Get_UserData(io)
    x == :Fonts && return ImGuiIO_Get_Fonts(io)
    x == :FontGlobalScale && return ImGuiIO_Get_FontGlobalScale(io)
    x == :FontAllowUserScaling && return ImGuiIO_Get_FontAllowUserScaling(io)
    x == :FontDefault && return ImGuiIO_Get_FontDefault(io)
    x == :DisplayFramebufferScale && return ImGuiIO_Get_DisplayFramebufferScale(io)
    x == :MouseDrawCursor && return ImGuiIO_Get_MouseDrawCursor(io)
    x == :ConfigMacOSXBehaviors && return ImGuiIO_Get_ConfigMacOSXBehaviors(io)
    x == :ConfigInputTextCursorBlink && return ImGuiIO_Get_ConfigInputTextCursorBlink(io)
    x == :ConfigWindowsResizeFromEdges && return ImGuiIO_Get_ConfigWindowsResizeFromEdges(io)
    x == :ConfigWindowsMoveFromTitleBarOnly && return ImGuiIO_Get_ConfigWindowsMoveFromTitleBarOnly(io)
    x == :BackendPlatformName && return ImGuiIO_Get_BackendPlatformName(io)
    x == :BackendRendererName && return ImGuiIO_Get_BackendRendererName(io)
    x == :ImeWindowHandle && return ImGuiIO_Get_ImeWindowHandle(io)
    x == :MousePos && return ImGuiIO_Get_MousePos(io)
    x == :MouseWheel && return ImGuiIO_Get_MouseWheel(io)
    x == :MouseWheelH && return ImGuiIO_Get_MouseWheelH(io)
    x == :KeyCtrl && return ImGuiIO_Get_KeyCtrl(io)
    x == :KeyShift && return ImGuiIO_Get_KeyShift(io)
    x == :KeyAlt && return ImGuiIO_Get_KeyAlt(io)
    x == :KeySuper && return ImGuiIO_Get_KeySuper(io)
    x == :WantCaptureMouse && return ImGuiIO_Get_WantCaptureMouse(io)
    x == :WantCaptureKeyboard && return ImGuiIO_Get_WantCaptureKeyboard(io)
    x == :WantTextInput && return ImGuiIO_Get_WantTextInput(io)
    x == :WantSetMousePos && return ImGuiIO_Get_WantSetMousePos(io)
    x == :WantSaveIniSettings && return ImGuiIO_Get_WantSaveIniSettings(io)
    x == :NavActive && return ImGuiIO_Get_NavActive(io)
    x == :NavVisible && return ImGuiIO_Get_NavVisible(io)
    x == :Framerate && return ImGuiIO_Get_Framerate(io)
    x == :MetricsRenderVertices && return ImGuiIO_Get_MetricsRenderVertices(io)
    x == :MetricsRenderIndices && return ImGuiIO_Get_MetricsRenderIndices(io)
    x == :MetricsRenderWindows && return ImGuiIO_Get_MetricsRenderWindows(io)
    x == :MetricsActiveWindows && return ImGuiIO_Get_MetricsActiveWindows(io)
    x == :MetricsActiveAllocations && return ImGuiIO_Get_MetricsActiveAllocations(io)
    x == :MouseDelta && return ImGuiIO_Get_MouseDelta(io)
    x == :MousePosPrev && return ImGuiIO_Get_MousePosPrev(io)
    throw(ArgumentError("field $x is not supported to be used like this, please use `Get_$x` instead."))
end

function Base.setproperty!(io::Ptr{ImGuiIO}, x::Symbol, v)
    x == :ConfigFlags && return ImGuiIO_Set_ConfigFlags(io, v)
    x == :BackendFlags && return ImGuiIO_Set_BackendFlags(io, v)
    x == :DisplaySize && return ImGuiIO_Set_DisplaySize(io, v)
    x == :DeltaTime && return ImGuiIO_Set_DeltaTime(io, v)
    x == :IniSavingRate && return ImGuiIO_Set_IniSavingRate(io, v)
    x == :IniFilename && return ImGuiIO_Set_IniFilename(io, v)
    x == :LogFilename && return ImGuiIO_Set_LogFilename(io, v)
    x == :MouseDoubleClickTime && return ImGuiIO_MouseDoubleClickTime(io, v)
    x == :MouseDoubleClickMaxDist && return ImGuiIO_Set_MouseDoubleClickMaxDist(io, v)
    x == :MouseDragThreshold && return ImGuiIO_MouseDragThreshold(io, v)
    x == :KeyRepeatDelay && return ImGuiIO_Set_KeyRepeatDelay(io, v)
    x == :KeyRepeatRate && return ImGuiIO_Set_KeyRepeatRate(io, v)
    x == :UserData && return ImGuiIO_Set_UserData(io, v)
    x == :Fonts && return ImGuiIO_Set_Fonts(io, v)
    x == :FontGlobalScale && return ImGuiIO_Set_FontGlobalScale(io, v)
    x == :FontAllowUserScaling && return ImGuiIO_Set_FontAllowUserScaling(io, v)
    x == :FontDefault && return ImGuiIO_Set_FontDefault(io, v)
    x == :DisplayFramebufferScale && return ImGuiIO_Set_DisplayFramebufferScale(io, v)
    x == :MouseDrawCursor && return ImGuiIO_Set_MouseDrawCursor(io, v)
    x == :ConfigMacOSXBehaviors && return ImGuiIO_Set_ConfigMacOSXBehaviors(io, v)
    x == :ConfigInputTextCursorBlink && return ImGuiIO_Set_ConfigInputTextCursorBlink(io, v)
    x == :ConfigWindowsResizeFromEdges && return ImGuiIO_Set_ConfigWindowsResizeFromEdges(io, v)
    x == :ConfigWindowsMoveFromTitleBarOnly && return ImGuiIO_Set_ConfigWindowsMoveFromTitleBarOnly(io, v)
    x == :BackendPlatformName && return ImGuiIO_Set_BackendPlatformName(io, v)
    x == :BackendRendererName && return ImGuiIO_Set_BackendRendererName(io, v)
    x == :GetClipboardTextFn && return ImGuiIO_Set_GetClipboardTextFn(io, v)
    x == :SetClipboardTextFn && return ImGuiIO_Set_SetClipboardTextFn(io, v)
    x == :ClipboardUserData && return ImGuiIO_Set_ClipboardUserData(io, v)
    x == :ImeWindowHandle && return ImGuiIO_Set_ImeWindowHandle(io, v)
    x == :MousePos && return ImGuiIO_Set_MousePos(io, v)
    x == :MouseWheel && return ImGuiIO_Set_MouseWheel(io, v)
    x == :MouseWheelH && return ImGuiIO_Set_MouseWheelH(io, v)
    x == :KeyCtrl && return ImGuiIO_Set_KeyCtrl(io, v)
    x == :KeyShift && return ImGuiIO_Set_KeyShift(io, v)
    x == :KeyAlt && return ImGuiIO_Set_KeyAlt(io, v)
    x == :KeySuper && return ImGuiIO_Set_KeySuper(io, v)
    throw(ArgumentError("field $x is not supported to be used like this!"))
end

# emulate ImGui::GetStyle().xxx
function Base.getproperty(s::Ptr{ImGuiStyle}, x::Symbol)
    x == :Alpha && return ImGuiStyle_Get_Alpha(s)
    x == :WindowPadding && return ImGuiStyle_Get_WindowPadding(s)
    x == :WindowRounding && return ImGuiStyle_Get_WindowRounding(s)
    x == :WindowBorderSize && return ImGuiStyle_Get_WindowBorderSize(s)
    x == :WindowMinSize && return ImGuiStyle_Get_WindowMinSize(s)
    x == :WindowTitleAlign && return ImGuiStyle_Get_WindowTitleAlign(s)
    x == :ChildRounding && return ImGuiStyle_Get_ChildRounding(s)
    x == :ChildBorderSize && return ImGuiStyle_Get_ChildBorderSize(s)
    x == :PopupRounding && return ImGuiStyle_Get_PopupRounding(s)
    x == :PopupBorderSize && return ImGuiStyle_Get_PopupBorderSize(s)
    x == :FramePadding && return ImGuiStyle_Get_FramePadding(s)
    x == :FrameRounding && return ImGuiStyle_Get_FrameRounding(s)
    x == :FrameBorderSize && return ImGuiStyle_Get_FrameBorderSize(s)
    x == :ItemSpacing && return ImGuiStyle_Get_ItemSpacing(s)
    x == :ItemInnerSpacing && return ImGuiStyle_Get_ItemInnerSpacing(s)
    x == :TouchExtraPadding && return ImGuiStyle_Get_TouchExtraPadding(s)
    x == :IndentSpacing && return ImGuiStyle_Get_IndentSpacing(s)
    x == :ColumnsMinSpacing && return ImGuiStyle_Get_ColumnsMinSpacing(s)
    x == :ScrollbarSize && return ImGuiStyle_Get_ScrollbarSize(s)
    x == :ScrollbarRounding && return ImGuiStyle_Get_ScrollbarRounding(s)
    x == :GrabMinSize && return ImGuiStyle_Get_GrabMinSize(s)
    x == :GrabRounding && return ImGuiStyle_Get_GrabRounding(s)
    x == :TabRounding && return ImGuiStyle_Get_TabRounding(s)
    x == :TabBorderSize && return ImGuiStyle_Get_TabBorderSize(s)
    x == :ButtonTextAlign && return ImGuiStyle_Get_ButtonTextAlign(s)
    x == :SelectableTextAlign && return ImGuiStyle_Get_SelectableTextAlign(s)
    x == :DisplayWindowPadding && return ImGuiStyle_Get_DisplayWindowPadding(s)
    x == :DisplaySafeAreaPadding && return ImGuiStyle_Get_DisplaySafeAreaPadding(s)
    x == :MouseCursorScale && return ImGuiStyle_Get_MouseCursorScale(s)
    x == :AntiAliasedLines && return ImGuiStyle_Get_AntiAliasedLines(s)
    x == :AntiAliasedFill && return ImGuiStyle_Get_AntiAliasedFill(s)
    x == :CurveTessellationTol && return ImGuiStyle_Get_CurveTessellationTol(s)
    throw(ArgumentError("field $x is not supported to be used like this, please use `Get_$x` instead."))
end

function Base.setproperty!(s::Ptr{ImGuiStyle}, x::Symbol, v)
    x == :Alpha && return ImGuiStyle_Get_Alpha(s, v)
    x == :WindowPadding && return ImGuiStyle_Get_WindowPadding(s, v)
    x == :WindowRounding && return ImGuiStyle_Get_WindowRounding(s, v)
    x == :WindowBorderSize && return ImGuiStyle_Get_WindowBorderSize(s, v)
    x == :WindowMinSize && return ImGuiStyle_Get_WindowMinSize(s, v)
    x == :WindowTitleAlign && return ImGuiStyle_Get_WindowTitleAlign(s, v)
    x == :ChildRounding && return ImGuiStyle_Get_ChildRounding(s, v)
    x == :ChildBorderSize && return ImGuiStyle_Get_ChildBorderSize(s, v)
    x == :PopupRounding && return ImGuiStyle_Get_PopupRounding(s, v)
    x == :PopupBorderSize && return ImGuiStyle_Get_PopupBorderSize(s, v)
    x == :FramePadding && return ImGuiStyle_Get_FramePadding(s, v)
    x == :FrameRounding && return ImGuiStyle_Get_FrameRounding(s, v)
    x == :FrameBorderSize && return ImGuiStyle_Get_FrameBorderSize(s, v)
    x == :ItemSpacing && return ImGuiStyle_Get_ItemSpacing(s, v)
    x == :ItemInnerSpacing && return ImGuiStyle_Get_ItemInnerSpacing(s, v)
    x == :TouchExtraPadding && return ImGuiStyle_Get_TouchExtraPadding(s, v)
    x == :IndentSpacing && return ImGuiStyle_Get_IndentSpacing(s, v)
    x == :ColumnsMinSpacing && return ImGuiStyle_Get_ColumnsMinSpacing(s, v)
    x == :ScrollbarSize && return ImGuiStyle_Get_ScrollbarSize(s, v)
    x == :ScrollbarRounding && return ImGuiStyle_Get_ScrollbarRounding(s, v)
    x == :GrabMinSize && return ImGuiStyle_Get_GrabMinSize(s, v)
    x == :GrabRounding && return ImGuiStyle_Get_GrabRounding(s, v)
    x == :TabRounding && return ImGuiStyle_Get_TabRounding(s, v)
    x == :TabBorderSize && return ImGuiStyle_Get_TabBorderSize(s, v)
    x == :ButtonTextAlign && return ImGuiStyle_Get_ButtonTextAlign(s, v)
    x == :SelectableTextAlign && return ImGuiStyle_Get_SelectableTextAlign(s, v)
    x == :DisplayWindowPadding && return ImGuiStyle_Get_DisplayWindowPadding(s, v)
    x == :DisplaySafeAreaPadding && return ImGuiStyle_Get_DisplaySafeAreaPadding(s, v)
    x == :MouseCursorScale && return ImGuiStyle_Get_MouseCursorScale(s, v)
    x == :AntiAliasedLines && return ImGuiStyle_Get_AntiAliasedLines(s, v)
    x == :AntiAliasedFill && return ImGuiStyle_Get_AntiAliasedFill(s, v)
    x == :CurveTessellationTol && return ImGuiStyle_Get_CurveTessellationTol(s, v)
    throw(ArgumentError("field $x is not supported to be used like this!"))
end

# emulate draw_data->xxx
function Base.getproperty(data::Ptr{ImDrawData}, x::Symbol)
    x == :Valid && return ImDrawData_Get_Valid(data)
    x == :CmdListsCount && return ImDrawData_Get_CmdListsCount(data)
    x == :TotalIdxCount && return ImDrawData_Get_TotalIdxCount(data)
    x == :TotalVtxCount && return ImDrawData_Get_TotalVtxCount(data)
    x == :DisplayPos && return ImDrawData_Get_DisplayPos(data)
    x == :DisplaySize && return ImDrawData_Get_DisplaySize(data)
    x == :FramebufferScale && return ImDrawData_Get_FramebufferScale(data)
    throw(ArgumentError("field $x is not supported to be used like this, please use `Get_$x` instead."))
end

# emulate cmd.xxx
function Base.getproperty(cmd::Ptr{ImDrawCmd}, x::Symbol)
    x == :ElemCount && return ImDrawCmd_Get_ElemCount(cmd)
    x == :ClipRect && return ImDrawCmd_Get_ClipRect(cmd)
    x == :TextureId && return ImDrawCmd_Get_TextureId(cmd)
    x == :VtxOffset && return ImDrawCmd_Get_VtxOffset(cmd)
    x == :IdxOffset && return ImDrawCmd_Get_IdxOffset(cmd)
    x == :UserCallback && return ImDrawCmd_Get_UserCallback(cmd)
    x == :UserCallbackData && return ImDrawCmd_Get_UserCallbackData(cmd)
    throw(ArgumentError("field $x is not supported to be used like this, please use `Get_$x` instead."))
end

# emulate io.Fonts->.xxx
function Base.getproperty(f::Ptr{ImFontAtlas}, x::Symbol)
    x == :TexID && return ImFontAtlas_Get_TexID(f)
    x == :TexWidth && return ImFontAtlas_Get_TexWidth(f)
    x == :TexHeight && return ImFontAtlas_Get_TexHeight(f)
    throw(ArgumentError("field $x is not supported to be used like this, please use `Get_$x` instead."))
end


text:
#############################
# Wrappers for easy calling
#############################
linesearch(prob::Problem{L}, z::Vector{T}, v::Vector{T}) where {T<:Real, L} =
    linesearch(L, prob.y, z, v)

linesearch(prob::Problem{L, B, CS}, z::Vector{T}, ω::B) where {T<:Real, L, B, CS} =
    linesearch(L, prob.y, z, predictor(prob, ω))

linesearch(prob::Problem{L, B}, ωs::Vector{B}, ω::B) where {L, B} =
    linesearch(L, prob.y, prediction(prob, ωs), ω)

#############################
# For various Losses
#############################
linesearch(::Type{L}, args...) where {L<:AbLoss} = error("Not Implemented.")

linesearch(::Type{MSELoss}, y, z, v) = ((y-z)'*v)/(v'*v)

function linesearch(::Type{LogitLoss},
                    y::Vector{T},
                    z::Vector{T},
                    v::Vector{T},
                    ) where T<:Real
    ρ = zero(T)
    for i in 1:10
        s = sigmoid.(z+ρ*v)
        grad = v'*(s-y)
        hess = v'*(v.*s.*(1-s))
        dec = grad/hess
        ρ -= dec
        @debug "$i) $dec $ρ"
        grad*dec < 1e-6 && break         # Newton Decrement
    end
    ρ
end

function weighted_median(data::Vector{T1}, weights::Vector{T})::T1 where {T1,T<:Real}
    ordered = sortperm(data)
    cumordwts = cumsum(weights[ordered])
    data[ordered][findfirst(cumordwts .≥ cumordwts[end]/2)]
end

linesearch(::Type{MADLoss}, y, z, v) = weighted_median((y-z)./v, abs.(v))


text:
"""
    elimGaussWithPivoting(A,B)

adalah fungsi untuk menyelesaikan SPL `Ax=B` dengan `A` adalah matriks sembarang
dan `B` adalah vektor ruas kanan menggunakan metode eliminasi Gauss dengan pivoting.

# Examples
```jl
julia> A = [ 2 5 4 2 3
            4 1 0.4 6 5
            0.9 4 8 7 10
            5 7 0.5 2 7
            2 6 8 9 4];

julia> B = [26; 24.4; 32.8; 30.5; 53];

julia> Xn = elimGaussWithPivoting(A,B)
5×1 Array{Float64,2}:
  1.9999999999999993
  3.000000000000001
  1.0
  3.0
 -1.0000000000000007
```
return solusi SPL `X`

"""
function elimGaussWithPivoting(A,b)
    # Hitung ukuran matriks dan inisiasi solusi X
    n = length(b);
    X = zeros(n,1);
    # Buatlah matriks Gandeng [A|b]
    Aug = [A b];
    # Lakukan operasi baris dasar terhadap matriks gandeng
    for p = 1:n-1
        # Lakukan operasi pindah baris untuk menentukan nilai pivot
        val,j = findmax(abs.(Aug[p:n,p]));
        C = Aug[p,:];
        Aug[p,:] = Aug[j+p-1,:];
        Aug[j+p-1,:] = C;
        # Jika pivot bernilai nol, maka gagal.
        if Aug[p,p] == 0;
            error("Matriks tidak punya solusi tunggal")
        end
        # Lakukan eliminasi menggunakan operasi Eij(k)
        for i = p+1:n
            k = Aug[i,p]/Aug[p,p];
            Aug[i,:] = Aug[i,:] - k*Aug[p,:];
        end
    end
    # Pisahkan matriks gandeng, kemudian substitusi mundur.
    A = Aug[:,1:n];
    b = Aug[:,1+n];
    X = backsub(A,b);
end


text:
module JungleHelperFirefly

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/Firefly" Firefly(x::Integer, y::Integer, number::Integer=1, sprite::String="")

const placements = Ahorn.PlacementDict(
    "Firefly (Jungle Helper)" => Ahorn.EntityPlacement(
        Firefly
    )
)

sprite = "JungleHelper/Firefly/firefly00"

function Ahorn.selection(entity::Firefly)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Firefly, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 0)

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
@recipe function f(sol::AbstractFVSolution; tidx = size(gettimes(sol),1), vars=nothing)
    seriestype  :=  :path
    xguide --> "x"
    yguide --> "u"
    labels = String[]
    nvars = getnvars(sol)
    if nvars == 1
      push!(labels,"u")
    else
      for i in 1:nvars
        push!(labels,"u$i")
      end
    end
    yvector = getvalues(sol)[tidx]
    if vars != nothing && nvars > 1
      yvector = getvalues(sol)[tidx][vars,:]
      labels = labels[vars]
    end
    if typeof(labels) <: String
      label --> labels
    else
      label --> reshape(labels,1,length(labels))
    end
    cell_centers(getmesh(sol)), yvector
end

# @recipe function f(sol::DGSolution; tidx = size(sol.t,1), vars=nothing)
#     xguide --> "x"
#     yguide --> "u"
#     labels = String[]
#     for i in 1:size(sol.u[tidx],2)
#       push!(labels,"u$i")
#     end
#     yvector = sol.u[tidx]
#     if vars != nothing
#       yvector = sol.u[tidx][:,vars]
#       labels = labels[vars]
#     end
#     if typeof(labels) <: String
#       label --> labels
#     else
#       label --> reshape(labels,1,length(labels))
#     end
#     @series begin
#       seriestype  :=  :path
#       sol.nodes, yvector
#     end
# end

# @recipe function f(ooc::FVOOCTable)
#     seriestype  :=  :path
#     xguide --> "M"
#     yguide --> "error"
#     xscale --> :log
#     yscale --> :log
#     lab --> ooc.alg_name
#     ooc.h, ooc.errors
# end


text:
macro ArfFloat(x)
    convert(ArfFloat, string(:($x)))
end
macro ArfFloat(p,x)
    convert(ArfFloat{:($p)}, string(:($x)))
end

macro ArbFloat(x)
    quote
       if isa($x, Irrational)
           convert(ArbFloat, $x)
       else
            convert(ArbFloat, string($x))
       end
    end
end

macro ArbFloat(p,x)
    quote
       if isa($x, Irrational)
           convert(ArbFloat{$p}, $x)
       else
           convert(ArbFloat{$p}, string($x))
       end
    end
end

convert(::Type{T}, x::T) where {T <: ArfFloat} = x
convert(::Type{T}, x::T) where {T <: ArbFloat} = x
convert(::Type{ArfFloat{P}}, x::ArfFloat{P}) where {P} = x
convert(::Type{ArbFloat{P}}, x::ArbFloat{P}) where {P} = x

#=
function convert{Q}(::Type{ArfFloat}, x::ArfFloat{Q})
   P = precision(ArfFloat)
   z = initializer(ArfFloat{P})
   ccall(@libarb(arf_set_round), Void, (Ref{ArfFloat{P}}, Ref{ArfFloat{Q}}, Clong), z, x, Clong(P))
   return z
end
=#
function convert(::Type{ArfFloat{P}}, x::ArfFloat{Q}) where {P,Q}
   z = initializer(ArfFloat{P})
   ccall(@libarb(arf_set_round), Cvoid, (Ref{ArfFloat{P}}, Ref{ArfFloat{Q}}, Clong), z, x, Clong(P))
   return z
end
#=
function convert{Q}(::Type{ArbFloat}, x::ArbFloat{Q})
   P = precision(ArbFloat)
   z = initializer(ArbFloat{P})
   ccall(@libarb(arb_set_round), Void, (Ref{ArbFloat{P}}, Ref{ArbFloat{Q}}, Clong), z, x, Clong(P))
   return z
end
=#
function convert(::Type{ArbFloat{P}}, x::ArbFloat{Q}) where {P,Q}
   z = initializer(ArbFloat{P})
   ccall(@libarb(arb_set_round), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{Q}}, Clong), z, x, Clong(P))
   return z
end

function convert(::Type{ArbFloat{P}}, x::ArfFloat{P}) where {P}
   z = initializer(ArbFloat{P})
   ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbFloat{P}}, Ref{ArfFloat{P}}), z, x)
   return z
end

function convert(::Type{ArfFloat{P}}, x::ArbFloat{P}) where {P}
    z = initializer(ArfFloat{P})
    z.exponentOf2  = x.exponentOf2
    z.nwords_sign  = x.nwords_sign
    z.significand1 = x.significand1
    z.significand2 = x.significand2
    return z
end

function convert(::Type{ArbFloat{P}}, x::ArfFloat{Q}) where {P,Q}
    y = convert(ArfFloat{P}, x)
    z = convert(ArbFloat{P}, y)
    return z
end
#=
function convert{Q}(::Type{ArbFloat}, x::ArfFloat{Q})
    P = precision(ArbFloat)
    y = convert(ArfFloat{P}, x)
    z = convert(ArbFloat{P}, y)
    return z
end
=#
function convert(::Type{ArfFloat{P}}, x::ArbFloat{Q}) where {P,Q}
    y = convert(ArbFloat{P}, x)
    z = convert(ArfFloat{P}, y)
    return z
end
function convert(::Type{ArfFloat}, x::ArbFloat{Q}) where {Q}
    P = precision(ArfFloat)
    y = convert(ArbFloat{P}, x)
    z = convert(ArfFloat{P}, y)
    return z
end

# convert ArbFloat with other types

function convert(::Type{T}, x::UInt64) where {T <: ArbFloat}
    P = precision(T)
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_set_ui), Cvoid, (Ref{ArbFloat}, UInt64), z, x)
    return z
end
function convert(::Type{T}, x::Int64) where {T <: ArbFloat}
    P = precision(T)
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_set_si), Cvoid, (Ref{ArbFloat}, Int64), z, x)
    return z
end

function convert(::Type{T}, x::Float64) where {T <: ArbFloat}
    P = precision(T)
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_set_d), Cvoid, (Ref{ArbFloat}, Float64), z, x)
    return z
end

function convert(::Type{T}, x::String) where {T <: ArbFloat}
    P = precision(T)
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_set_str), Cvoid, (Ref{ArbFloat}, Ptr{UInt8}, Int), z, x, P)
    return z
end

convert(::Type{T}, x::UInt32) where {T <: ArbFloat}  = convert(T, x%UInt64)
convert(::Type{T}, x::UInt16) where {T <: ArbFloat}  = convert(T, x%UInt64)
convert(::Type{T}, x::UInt8) where {T <: ArbFloat}   = convert(T, x%UInt64)
convert(::Type{T}, x::UInt128) where {T <: ArbFloat} = convert(T, string(x))

convert(::Type{T}, x::Int32) where {T <: ArbFloat}  = convert(T, x%Int64)
convert(::Type{T}, x::Int16) where {T <: ArbFloat}  = convert(T, x%Int64)
convert(::Type{T}, x::Int8) where {T <: ArbFloat}   = convert(T, x%Int64)
convert(::Type{T}, x::Int128) where {T <: ArbFloat} = convert(T, string(x))

convert(::Type{T}, x::Float32) where {T <: ArbFloat} = convert(T, convert(Float64,x))
convert(::Type{T}, x::Float16) where {T <: ArbFloat} = convert(T, convert(Float64,x))

function convert(::Type{Float64}, x::T) where {T <: ArbFloat}
    ptr2mid = ptr_to_midpoint(x)
    fl = ccall(@libarb(arf_get_d), Float64, (Ref{ArfFloat}, Int), ptr2mid, 4) # round nearest
    return fl
end

function convert(::Type{Float32}, x::T) where {T <: ArbFloat}
    return convert(Float32, convert(Float64, x))
end
function convert(::Type{Float16}, x::T) where {T <: ArbFloat}
    return convert(Float16, convert(Float64, x))
end

for I in (:UInt64, :UInt128)
  @eval begin
    function convert(::Type{$I}, x::T) where {T <: ArbFloat}
        if isinteger(x)
           if notnegative(x)
               return convert($I, convert(BigInt,x))
           else
               throw( DomainError() )
           end
        else
           throw( InexactError() )
        end
    end
  end
end

convert(::Type{UInt32}, x::T) where {T <: ArbFloat} = convert(UInt32, convert(UInt64,x))
convert(::Type{UInt16}, x::T) where {T <: ArbFloat} = convert(UInt16, convert(UInt64,x))
convert(::Type{UInt8}, x::T) where {T <: ArbFloat} = convert(UInt8, convert(UInt64,x))

for I in (:Int64, :Int128)
  @eval begin
    function convert(::Type{$I}, x::T) where {T <: ArbFloat}
        if isinteger(x)
           return convert($I, convert(BigInt,x))
        else
           throw( InexactError() )
        end
    end
  end
end

convert(::Type{Int32}, x::T) where {T <: ArbFloat} = convert(Int32, convert(Int64,x))
convert(::Type{Int16}, x::T) where {T <: ArbFloat} = convert(Int16, convert(Int64,x))
convert(::Type{Int8}, x::T) where {T <: ArbFloat} = convert(Int8, convert(Int64,x))

function parse(::Type{T}, x::String) where {T <: ArbFloat}
    return T(x)
end

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

convert(::Type{BigInt}, x::String) = parse(BigInt,x)
convert(::Type{BigFloat}, x::String) = parse(BigFloat,x)

function convert(::Type{T}, x::BigFloat) where {T <: ArbFloat}
     P = precision(T)+24
     x = round(x, digits=P, base=2)
     s = string(x)
     z = T(s)
     return z
end

#=
function convert{T<:ArbFloat}(::Type{BigFloat}, x::T)
     s = string(midpoint(x))
     return parse(BigFloat, s)
end
=#
function convert(::Type{BigFloat}, x::T) where {T <: ArbFloat}
    ptr2mid = ptr_to_midpoint(x)
    bf = zero(BigFloat)
    rounddir = ccall(@libarb(arf_get_mpfr), Int, (Ref{BigFloat}, Ref{ArfFloat}, Int), bf, ptr2mid, 4) # round nearest
    return bf
end

function convert(::Type{Rational{I}}, x::ArbFloat{P}) where {I <: Integer,P}
    bf = convert(BigFloat, x)
    return convert(Rational{I}, bf)
end

for T in (:Integer, :Signed)
  @eval begin
    function convert(::Type{$T}, x::ArbFloat{P}) where {P}
        y = trunc(x)
        try
           return convert(Int64, x)
        catch
           try
              return convert(Int128, x)
           catch
              DomainError()
           end
        end
    end
  end
end

for F in (:BigInt, :Rational)
  @eval begin
    function convert(::Type{T}, x::$F) where {T <: ArbFloat}
        P = precision(T)
        B = precision(BigFloat)
        if B < P+24
            return convert(ArbFloat{P}, string(x))
        else
            return convert(ArbFloat{P}, convert(BigFloat, x))
        end
    end
    function convert(::Type{ArbFloat{P}}, x::$F) where {P}
        B = precision(BigFloat)
        if B < P+24
            return convert(ArbFloat{P}, string(x))
        else
            return convert(ArbFloat{P}, convert(BigFloat, x))
        end
    end
  end
end

function convert(::Type{T}, x::Irrational{Sym}) where {T <: ArbFloat,Sym}
    P = precision(T)
    bf_precision = precision(BigFloat)
    setprecision(BigFloat, precision(T)+24)
    bf_x = convert(BigFloat, x)
    z = convert(ArbFloat{P}, bf_x)
    setprecision(BigFloat, bf_precision)
    return z
end

function convert(::Type{BigInt}, x::ArbFloat{P}) where {P}
   z = trunc(convert(BigFloat, x))
   return convert(BigInt, z)
end

# Promotion
for T in (:Int128, :Int64, :Int32, :Int16, :Float64, :Float32, :Float16,
          :(Rational{Int64}), :(Rational{Int32}), :(Rational{Int16}),
          :String)
  @eval promote_rule(::Type{ArbFloat{P}}, ::Type{$T}) where {P} = ArbFloat
end

float(x::ArbFloat{P}) where {P} = x

promote_rule(::Type{ArbFloat{P}}, ::Type{BigFloat}) where {P} = ArbFloat
promote_rule(::Type{ArbFloat{P}}, ::Type{BigInt}) where {P} = ArbFloat
promote_rule(::Type{ArbFloat{P}}, ::Type{Rational{BigInt}}) where {P} = Rational{BigInt}

promote_rule(::Type{ArbFloat{P}}, ::Type{ArbFloat{Q}}) where {P,Q} =
    ifelse(P>Q, ArbFloat{P}, ArbFloat{Q})

# convert a vector of ArbFloats to another numeric type
for T in (:BigFloat, :Float64, :Float32, :BigInt, :Int128, :Int64, :Int32, :Rational)
    @eval ($T)(x::Array{ArbFloat{P},1}) where {P} = ($T).(x)
end
# convert a numeric vector into ArbFloats
for T in (:Float64, :Float32, :BigInt, :Int128, :Int64, :Int32, :Rational)
    @eval begin
        function ArbFloat(x::Array{($T),1})
             P=precision(ArbFloat)
             return map(v->ArbFloat{P}(v), x)
        end
    end
end


text:
@testset "Penalty Parameter Interface" begin
    struct Fake <: PH.AbstractPenaltyParameter end
    fake = Fake()
    xhid = PH.XhatID(PH.NodeID(0), PH.Index(0))
    phd = fake_phdata(3, 4; add_leaves=false)
    @test_throws PH.UnimplementedError PH.get_penalty_value(fake)
    @test_throws PH.UnimplementedError PH.get_penalty_value(fake, xhid)
    @test_throws PH.UnimplementedError PH.is_initial_value_dependent(Fake)
    @test_throws PH.UnimplementedError PH.is_subproblem_dependent(Fake)
    @test_throws PH.UnimplementedError PH.is_variable_dependent(Fake)
    @test_throws PH.UnimplementedError PH.penalty_map(fake)
    @test_throws PH.UnimplementedError PH.process_penalty_initial_value(fake, phd)
    @test_throws PH.UnimplementedError PH.process_penalty_subproblem(fake,
                                                                     phd,
                                                                     PH.scid(0),
                                                                     Dict{PH.VariableID,
                                                                          Float64}()
                                                                  )
end

@testset "Scalar Penalty" begin
    @test PH.is_initial_value_dependent(PH.ScalarPenaltyParameter) == false
    @test PH.is_subproblem_dependent(PH.ScalarPenaltyParameter) == false
    @test PH.is_variable_dependent(PH.ScalarPenaltyParameter) == false

    rflt = 10.0*rand()
    r = PH.ScalarPenaltyParameter(rflt)
    @test PH.get_penalty_value(r) == rflt
    @test PH.get_penalty_value(r, PH.XhatID(PH.NodeID(0), PH.Index(0))) == rflt
end

@testset "Proportional Penalty" begin
    @test PH.is_initial_value_dependent(PH.ProportionalPenaltyParameter) == false
    @test PH.is_subproblem_dependent(PH.ProportionalPenaltyParameter) == true
    @test PH.is_variable_dependent(PH.ProportionalPenaltyParameter) == true

    N = 2
    st = PH.two_stage_tree(N)
    var_map = Dict{PH.ScenarioID,Dict{PH.VariableID,PH.VariableInfo}}()
    coef_map = Dict{Int,Dict{PH.VariableID,Float64}}()
    for s in 0:(N-1)
        js = two_stage_model(PH.scid(s))
        svm = PH.report_variable_info(js, st)
        var_map[PH.scid(s)] = svm
        (br_vars, lf_vars) = PH._split_variables(st, collect(keys(svm)))
        coef_map[s] = PH.report_penalty_info(js,
                                             br_vars,
                                             PH.ProportionalPenaltyParameter
                                             )
    end

    c_rho = 10.0 * (rand() + 1e-12) # Avoid the possibility of c_rho == 0.0
    phd = PH.PHData(PH.ProportionalPenaltyParameter(c_rho),
                    st,
                    Dict(1 => Set([PH.scid(0), PH.scid(1)])),
                    var_map,
                    TimerOutputs.TimerOutput()
                    )

    for s in 0:(N-1)
        PH.process_penalty_subproblem(phd.r, phd, PH.scid(s), coef_map[s])
    end

    for (xhid, rho) in pairs(phd.r.penalties)
        @test st.tree_map[xhid.node].stage == PH.stid(1)
        @test isapprox(rho, c_rho) || isapprox(rho, phd.r.constant)
    end
end

@testset "SEP Penalty" begin
    @test PH.is_initial_value_dependent(PH.SEPPenaltyParameter) == true
    @test PH.is_subproblem_dependent(PH.SEPPenaltyParameter) == true
    @test PH.is_variable_dependent(PH.SEPPenaltyParameter) == true

    N = 2
    st = PH.two_stage_tree(N)
    var_map = Dict{PH.ScenarioID,Dict{PH.VariableID,PH.VariableInfo}}()
    var_vals = Dict{Int,Dict{PH.VariableID,Float64}}()
    coef_map = Dict{Int,Dict{PH.VariableID,Float64}}()
    for s in 0:(N-1)
        js = two_stage_model(PH.scid(s))
        svm = PH.report_variable_info(js, st)
        var_map[PH.scid(s)] = svm
        (br_vars, lf_vars) = PH._split_variables(st, collect(keys(svm)))
        PH.solve_subproblem(js)
        var_vals[s] = PH.report_values(js, collect(keys(svm)))
        coef_map[s] = PH.report_penalty_info(js,
                                             br_vars,
                                             PH.SEPPenaltyParameter
                                             )
    end

    phd = PH.PHData(PH.SEPPenaltyParameter(),
                    st,
                    Dict(1 => Set([PH.scid(0), PH.scid(1)])),
                    var_map,
                    TimerOutputs.TimerOutput()
                    )

    for s in 0:(N-1)
        PH._copy_values(phd.scenario_map[PH.scid(s)].branch_vars, var_vals[s])
        PH.process_penalty_subproblem(phd.r, phd, PH.scid(s), coef_map[s])
    end

    PH.compute_and_save_xhat(phd)
    PH.process_penalty_initial_value(phd.r, phd)

    value = 0.5 * (10.0 + 3.0)
    value = 0.5 * (10.0 - value) + 0.5 * (value - 3.0)
    rho_val = 1.0 / value
    for (xhid, rho) in pairs(phd.r.penalties)
        @test st.tree_map[xhid.node].stage == PH.stid(1)
        @test isapprox(rho, rho_val) || isapprox(rho, phd.r.default)
    end

    if isdefined(Main, :Xpress)
        N = 8
        st = PH.two_stage_tree(N)
        var_map = Dict{PH.ScenarioID,Dict{PH.VariableID,PH.VariableInfo}}()
        var_vals = Dict{Int,Dict{PH.VariableID,Float64}}()
        coef_map = Dict{Int,Dict{PH.VariableID,Float64}}()
        for s in 0:(N-1)
            js = two_stage_int(PH.scid(s))
            svm = PH.report_variable_info(js, st)
            var_map[PH.scid(s)] = svm
            (br_vars, lf_vars) = PH._split_variables(st, collect(keys(svm)))
            PH.solve_subproblem(js)
            var_vals[s] = PH.report_values(js, collect(keys(svm)))
            coef_map[s] = PH.report_penalty_info(js,
                                                 br_vars,
                                                 PH.SEPPenaltyParameter
                                                 )
        end

        phd = PH.PHData(PH.SEPPenaltyParameter(),
                        st,
                        Dict(1 => Set([PH.scid(s) for s in 0:(N-1)])),
                        var_map,
                        TimerOutputs.TimerOutput()
                        )

        for s in 0:(N-1)
            PH._copy_values(phd.scenario_map[PH.scid(s)].branch_vars, var_vals[s])
            PH.process_penalty_subproblem(phd.r, phd, PH.scid(s), coef_map[s])
        end

        PH.compute_and_save_xhat(phd)
        PH.process_penalty_initial_value(phd.r, phd)

        value = (10 - 4 + 1)
        rho_val = 2.0 / value
        for (xhid, rho) in pairs(phd.r.penalties)
            @test st.tree_map[xhid.node].stage == PH.stid(1)
            @test isapprox(rho, rho_val)
        end
    else
        println("Skipping SEP integer setting tests.")
    end
end

@testset "Reset penalty parameters" begin
    for r in (PH.ProportionalPenaltyParameter(1.0), PH.SEPPenaltyParameter(1.0))
        @test isempty(r.penalties) == true
        r.penalties[PH.XhatID(PH.NodeID(0), PH.Index(0))] = 1.0
        @test isempty(r.penalties) == false
        PH.reset_penalty_map(r)
        @test isempty(r.penalties) == true
    end
end