text:
export iauNumat
"""
**  Form the matrix of nutation.
**
**  This function is part of the International Astronomical Union's
**  SOFA (Standards Of Fundamental Astronomy) software collection.
**
**  Status:  support function.
**
**  Given:
**     epsa        double         mean obliquity of date (Note 1)
**     dpsi,deps   double         nutation (Note 2)
**
**  Returned:
**     rmatn       double[3][3]   nutation matrix (Note 3)
**
**  Notes:
**
**
**  1) The supplied mean obliquity epsa, must be consistent with the
**     precession-nutation models from which dpsi and deps were obtained.
**
**  2) The caller is responsible for providing the nutation components;
**     they are in longitude and obliquity, in radians and are with
**     respect to the equinox and ecliptic of date.
**
**  3) The matrix operates in the sense V(true) = rmatn * V(mean),
**     where the p-vector V(true) is with respect to the true
**     equatorial triad of date and the p-vector V(mean) is with
**     respect to the mean equatorial triad of date.
**
**  Called:
**     iauIr        initialize r-matrix to identity
**     iauRx        rotate around X-axis
**     iauRz        rotate around Z-axis
**
**  Reference:
**
**     Explanatory Supplement to the Astronomical Almanac,
**     P. Kenneth Seidelmann (ed), University Science Books (1992),
**     Section 3.222-3 (p114).
**
**  This revision:  2013 June 18
**
**  SOFA release 2018-01-30
**
**  Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauNumat(double epsa, double dpsi, double deps, double rmatn[3][3])
function iauNumat(epsa::Real, dpsi::Real, deps::Real)
    # Allocate return values
    rmatn = zeros(Float64, 3, 3)

    ccall((:iauNumat, libsofa_c), Cvoid,
         (Cdouble, Cdouble, Cdouble, Ptr{Cdouble},),
         convert(Float64, epsa), convert(Float64, dpsi), convert(Float64, deps),
         rmatn)

    return SMatrix{3,3}(rmatn')
 end


text:
using ReTest
include("GaPLACTests.jl")
GaPLACTests.runtests()

using GaPLAC
GaPLAC.runtests()

ext: jl

max_stars_repo_name: biobakery/GaPLAC

max_issues_repo_name: biobakery/GaPLAC

max_forks_repo_name: biobakery/GaPLAC


text:
function longer(a, b)
    try
        b[endof(a)]
    catch
        return true
    end
    return false
end

function printlongest(io::IO)
    lines = longest = ""
    while !eof(io)
        line = readline(io)
        if longer(line, longest)
            longest = lines = line
        elseif !longer(longest, line)
            lines *= "\n" * line
        end
    end
    println(lines)
end
printlongest(str::String) = printlongest(IOBuffer(str))

printlongest("a\nbb\nccc\nddd\nee\nf\nggg")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function pysparse_to_julia(A)
    A = sp.csc_matrix(A)
    m, n = A.shape
    colPtr = Int[i+1 for i in PyArray(A."indptr")]
    rowVal = Int[i+1 for i in PyArray(A."indices")]
    nzVal = A.data
    SparseMatrixCSC{eltype(nzVal),Int}(m, n, colPtr, rowVal, nzVal)
end

function qobj_to_jl(A,getρ=false)
    if A."type" == "oper"
        return getρ ? A.full() : pysparse_to_julia(A.data)
    elseif A."type" == "ket"
        ψ = reshape(A.full(),:)
        return getρ ? ψ * ψ' : ψ
    else
        error("Unkown type $(A."type") of Qobj!")
    end
end

function myqload(file,silent=true)
    # redirect python stdout to avoid "Loaded tuple object." message from qload
    if silent
        tmp = nothing # because of scoping rules, tmp has to be defined before assignment
        @pywith pyimport("contextlib").redirect_stdout(nothing) begin
            tmp = qt.qload(file)
        end
        tmp
    else
        qt.qload(file)
    end
end

function load_input_from_qutip(filename)
    data = myqload(filename)
    if length(data) == 4
        H_q, c_ops_q, ρ0_q, ts = data
        saveddata = :probabilities
    elseif length(data) == 5
        H_q, c_ops_q, ρ0_q, ts, e_ops_q = data
        saveddata = e_ops_q isa String ? Symbol(e_ops_q) : qobj_to_jl.(e_ops_q)
    else
        raise(ValueError("length(data) must be 4 or 5, got length(data) = $(length(data))."))
    end
    H = qobj_to_jl(H_q)
    J = qobj_to_jl.(c_ops_q)
    ρ0 = qobj_to_jl(ρ0_q,true)
    (H,J,ρ0,ts), saveddata
end


text:
module GR_Spherical

using DifferentialEquations
using BoundaryValueDiffEq
using OrdinaryDiffEq
using Fun1d
using DataFrames
using CSV
using Plots
using Roots

# using Random

struct GBSSN_Variables{S,T} <: AbstractArray{T,2}
    α::GridFun{S,T}
    A::GridFun{S,T}
    βr::GridFun{S,T}
    Br::GridFun{S,T}
    χ::GridFun{S,T}
    γtrr::GridFun{S,T}
    γtθθ::GridFun{S,T}
    Arr::GridFun{S,T}
    K::GridFun{S,T}
    Γr::GridFun{S,T}
    𝜙::GridFun{S,T}
    K𝜙::GridFun{S,T}
end

cont(x::GBSSN_Variables) = (x.α, x.A, x.βr, x.Br, x.χ, x.γtrr, x.γtθθ, x.Arr, x.K, x.Γr, x.𝜙, x.K𝜙)
numvar = 12

# Iteration
Base.IteratorSize(::Type{<:GBSSN_Variables}) = Iterators.HasShape{2}()
Base.eltype(::Type{GBSSN_Variables{S,T}}) where {S,T} = T
Base.isempty(x::GBSSN_Variables) = isempty(x.α)
function Base.iterate(x::GBSSN_Variables, state...)
    return iterate(Iterators.flatten(cont(x)), state...)
end
Base.size(x::GBSSN_Variables) = (length(x.α), numvar)
Base.size(x::GBSSN_Variables, d) = size(x)[d]

# Indexing
function lin2cart(x::GBSSN_Variables, i::Number)
    n = length(x.α)
    return (i - 1) % n + 1, (i - 1) ÷ n + 1
end
Base.firstindex(x::GBSSN_Variables) = error("not implemented")
Base.getindex(x::GBSSN_Variables, i) = getindex(x, i.I...)
Base.getindex(x::GBSSN_Variables, i::Number) = getindex(x, lin2cart(x, i)...)
Base.getindex(x::GBSSN_Variables, i, j) = getindex(cont(x)[j], i)
Base.lastindex(x::GBSSN_Variables) = error("not implemented")
Base.setindex!(x::GBSSN_Variables, v, i) = setindex!(x, v, i.I...)
Base.setindex!(x::GBSSN_Variables, v, i::Number) = setindex!(x, v, lin2cart(x, i))
Base.setindex!(x::GBSSN_Variables, v, i, j) = setindex!(cont(x)[j], v, i)

# Abstract Array
Base.IndexStyle(::GBSSN_Variables) = IndexCartesian()
Base.similar(x::GBSSN_Variables) = GBSSN_Variables(map(similar, cont(x))...)
function Base.similar(x::GBSSN_Variables, ::Type{T}) where {T}
    return GBSSN_Variables(map(y -> similar(y,T), cont(x))...)
end
Base.similar(x::GBSSN_Variables, ::Dims) = similar(x)
Base.similar(x::GBSSN_Variables, ::Dims, ::Type{T}) where {T} = similar(x, T)

# Broadcasting
Base.BroadcastStyle(::Type{<:GBSSN_Variables}) = Broadcast.ArrayStyle{GBSSN_Variables}()
function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{GBSSN_Variables}},
                      ::Type{T}) where {T}
    x = find_GBSSN_Variables(bc)
    return similar(x, T)
end
find_GBSSN_Variables(bc::Base.Broadcast.Broadcasted) = find_GBSSN_Variables(bc.args)
find_GBSSN_Variables(args::Tuple) = find_GBSSN_Variables(find_GBSSN_Variables(args[1]), Base.tail(args))
find_GBSSN_Variables(x) = x
find_GBSSN_Variables(::Tuple{}) = nothing
find_GBSSN_Variables(a::GBSSN_Variables, rest) = a
find_GBSSN_Variables(::Any, rest) = find_GBSSN_Variables(rest)

# Others
function Base.map(fun, x::GBSSN_Variables, ys::GBSSN_Variables...)
    return GBSSN_Variables(
        map(fun, x.α, (y.α for y in ys)...),
        map(fun, x.A, (y.A for y in ys)...),
        map(fun, x.βr, (y.βr for y in ys)...),
        map(fun, x.Br, (y.Br for y in ys)...),
        map(fun, x.χ, (y.χ for y in ys)...),
        map(fun, x.γtrr, (y.γtrr for y in ys)...),
        map(fun, x.γtθθ, (y.γtθθ for y in ys)...),
        map(fun, x.Arr, (y.Arr for y in ys)...),
        map(fun, x.K, (y.K for y in ys)...),
        map(fun, x.Γr, (y.Γr for y in ys)...),
        map(fun, x.𝜙, (y.𝜙 for y in ys)...),
        map(fun, x.K𝜙, (y.K𝜙 for y in ys)...)
        )
end

# function Base.rand(rng::AbstractRNG, ::Random.SamplerType{GBSSN_Variables{T}}) where {T}
#     return GBSSN_Variables{T}(rand(rng, T), rand(rng, T))
# end

Base.zero(::Type{<:GBSSN_Variables}) = error("not implemented")
Base.zero(x::GBSSN_Variables) = GBSSN_Variables(map(zero, cont(x))...)

Base.:+(x::GBSSN_Variables) = map(+, x)
Base.:-(x::GBSSN_Variables) = map(-, x)

Base.:+(x::GBSSN_Variables, y::GBSSN_Variables) = map(+, x, y)
Base.:-(x::GBSSN_Variables, y::GBSSN_Variables) = map(-, x, y)

Base.:*(x::GBSSN_Variables, a::Number) = map(b -> b * a, x)
Base.:*(a::Number, x::GBSSN_Variables) = map(b -> a * b, x)
Base.:/(x::GBSSN_Variables, a::Number) = map(b -> b / a, x)
Base.:\(a::Number, x::GBSSN_Variables) = map(b -> a \ b, x)

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

function printlogo()
    # Just a fancy ASCII logo for the program
    println(
"\n",
"           /\\\\\\\\\\\\\\\\\\       /\\\\\\\\\\\\\\\\\\\\\\\\\n",
" ________/\\\\\\\\\\\\\\\\\\\\\\\\\\\\___/ /\\\\\\\\/ / /\\\\\\\\_____________________________\n",
"  ______//\\\\\\\\\\/ / / / /___\\/ /\\\\\\\\/_/ / /\\\\\\\\___________________________\n",
"   _____//\\\\\\\\/_/_/_/_/_____\\/ /\\\\\\\\ \\/ / /\\\\\\\\___________________________\n",
"    ____/ /\\\\\\\\____/\\\\\\\\\\\\___\\/ /\\\\\\\\\\\\\\\\\\\\\\\\\\\\_______in___________________\n",
"     ___\\/ /\\\\\\\\__/ /  /\\\\\\___\\/ /\\\\\\\\\\\\\\\\\\\\\\\\/________Spherical____________\n",
"      ___\\/ /\\\\\\\\_\\/__/ /\\\\\\___\\/ /\\\\\\\\ / /\\\\\\\\__________Symmetry____________\n",
"       ___\\/ /\\\\\\\\____\\/ /\\\\\\___\\/ /\\\\\\\\_/ //\\\\\\\\_____________________________\n",
"        ___\\/ /\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\___\\/ /\\\\\\\\\\// //\\\\\\\\____________________________\n",
"         ___\\/ / /\\\\\\\\\\\\\\\\\\  /____\\/ /\\\\\\\\_\\// //\\\\\\\\_____by Erik Schnetter_____\n",
"          ___\\/_/ / / / / /_/______\\/ /  /___\\// /  /_______and Conner Dailey____\n",
"           _____\\/_/_/_/_/__________\\/__/______\\/__/______________________________\n"
)

end

function setup(::Type{S}, rspan, points::Int) where {S}

    # Specify the radial domain and the
    # radial step amount and make a Grid object
    rmin,rmax = rspan
    domain = Domain{S}(rmin, rmax)
    grid = Grid(domain, points)
    return grid

end

function dissipation(f::GridFun{S,T}) where {S,T}

    ############################################
    # Calculates the numerical dissipation terms
    #
    # Finite differencing fails to model high
    # frequency modes in the system, and these
    # modes can lead to instabilities. Adding
    # numerical dissipation terms damps these
    # high frequency modes by subtracting off
    # a high order derivative from each dynamical
    # variable in the system.
    ############################################

    dx = spacing(f.grid)
    n = f.grid.ncells + 4
    dvalues = Array{T}(undef, n)

    dvalues[1:4] .= 0.

    # dvalues[3] = -((-48*f.values[3] + 96*f.values[4] - 48*f.values[5])/17)
    #
    # dvalues[4] = -((96*f.values[3] - 240*f.values[4] + 192*f.values[5]
    #     - 48*f.values[6])/59)
    #
    # dvalues[5] = -((-48*f.values[3] + 192*f.values[4] - 288*f.values[5]
    #     + 192*f.values[6] - 48*f.values[7])/43)
    #
    # dvalues[6] = -((-48*f.values[4] + 192*f.values[5] - 288*f.values[6]
    #     + 192*f.values[7] - 48*f.values[8])/49)

    for i in 5:(n - 4)
        dvalues[i] = (f.values[i-2] - 4*f.values[i-1] + 6*f.values[i]
        - 4*f.values[i+1] + f.values[i+2])
    end
    dvalues[(n-3):n] .= 0.

    return GridFun(f.grid, dvalues)
end

function init(::Type{T}, grid::Grid, param) where {T}

    ############################################
    # Specifies the Initial Conditions
    ############################################

    n = grid.ncells + 4
    domain = grid.domain
    initgrid = grid
    drt = spacing(grid)
    r = param[4]
    drdrt = param[5]
    d2rdrt = param[6]
    m = param[7]
    rtspan = param[8]

    num = 0

    # Initial conditions for Schwarzschild metric (Ker-Schild Coordinates)

    # Mass (no real reason not to use 1 here)
    #M = 1

    fα(M,rt) = real((1+2*M/(r(rt))+0im)^(-1/2))
    fA(rt) = 0.
    fβr(M,rt) = (2*M/r(rt))*fα(M,rt)^2
    fBr(rt) = 0.
    fχ(rt) = 1.
    fγtrr(M,rt) = 1+2*M/r(rt)
    fγtθθ(rt) = r(rt)^2
    fArr(M,∂M,rt) = (4/3)*(r(rt)*(M+r(rt))*∂M-M*(3*M+2*r(rt)))/real(((r(rt)^5)*(r(rt)+2*M)+0im)^(1/2))
    fK(M,∂M,rt) = (2*M*(3*M+r(rt))+2*r(rt)*∂M*(M+r(rt)))/real((r(rt)*(r(rt)+2*M)+0im)^(3/2))
    fΓr(M,∂M,rt) = (r(rt)*∂M-2*r(rt)-5*M)/(r(rt)+2*M)^2

    # fArr(M,rt) = -(4*M/3)*(3*M+2*r(rt))/real(((r(rt)^5)*(r(rt)+2*M)+0im)^(1/2))
    # fK(M,rt) = (2*M)*(3*M+r(rt))/real((r(rt)*(r(rt)+2*M)+0im)^(3/2))
    # fΓr(M,rt) = -(2*r(rt)+5*M)/(r(rt)+2*M)^2

    r0 = 10.
    σr = 0.5
    #Amp = 1.
    Amp = 0*0.1

    f𝜙(rt) = Amp*(1/r(rt))*exp(-(1/2)*((r(rt)-r0)/σr)^2)
    f∂𝜙(rt) = Amp*exp(-(1/2)*((r(rt)-r0)/σr)^2)*(r(rt)*r0-r(rt)^2-σr^2)/(r(rt)^2*σr^2)
    #f∂t𝜙(rt) = Amp*exp(-(1/2)*((r(rt)-r0)/σr)^2)*(r0-r(rt))/(r(rt)*σr^2)
    fK𝜙(rt) = 0.

    f∂χ(rt) = 0.
    f∂γtrr(M,∂M,rt) = 2*(r(rt)*∂M-M)/(r(rt)^2)
    f∂γtθθ(rt) = 2*r(rt)
    # f∂Arr(M,rt) = (4*M/3)*(15*M^2+15*M*r(rt)+4*r(rt)^2)/real(((r(rt)^7)*((r(rt)+2*M)^3)+0im)^(1/2))
    # f∂K(M,rt) = -2*M*(9*M^2+10*M*r(rt)+2*r(rt)^2)/real((r(rt)*(r(rt)+2*M)+0im)^(5/2))
    # f∂Γr(M,rt) = 2*(r(rt)+3*M)/(r(rt)+2*M)^3
    #
    f∂2γtθθ(rt) = 2.
    f∂2χ(rt) = 0.

    # fgtt(M,rt) = -(1/fα(M,rt)^2)
    # fgtr(M,rt) = fβr(M,rt)/fα(M,rt)^2
    # fgrr(M,χ,rt) = χ/fγtrr(M,rt) - (fβr(M,rt)/fα(M,rt))^2
    #
    # # Lagrangian Density for scalar field
    #
    # f𝓛(M,χ,rt) = ((1/2)*fgtt(M,rt)*f∂t𝜙(rt)^2 + (1/2)*fgrr(M,χ,rt)*f∂𝜙(rt)^2
    #  + fgtr(M,rt)*f∂𝜙(rt)*f∂t𝜙(rt) - (1/2)*(m^2)*f𝜙(rt)^2)
    #
    # # Stresss Energy components (contravariant indices)
    #
    # fTtt(M,χ,rt) = (fgtt(M,rt)*f∂t𝜙(rt) + fgtr(M,rt)*f∂𝜙(rt)).^2 - fgtt(M,rt)*f𝓛(M,χ,rt)
    # fTtr(M,χ,rt) = ((fgtt(M,rt)*f∂t𝜙(rt) + fgtr(M,rt)*f∂𝜙(rt))*(fgtr(M,rt)*f∂t𝜙(rt)
    #  + fgrr(M,χ,rt)*f∂𝜙(rt)) - fgtr(M,rt)*f𝓛(M,χ,rt))

    #fK𝜙(rt) = -(f∂t𝜙(rt) - fβr(rt)*f∂𝜙(rt))/(2*fα(M,rt))

    fρ(M,rt) = (2*fK𝜙(rt)^2 + (1/2)*(fχ(rt)/fγtrr(M,rt))*f∂𝜙(rt)^2
        + (1/2)*m^2*f𝜙(rt)^2)

    fSr(rt) = 2*fK𝜙(rt)*f∂𝜙(rt)

    #f∂χ(rt,(χ, ∂χ, M)) = ∂χ

    # function f∂2χ(rt,(χ, X, Arr))
    #  -(1/2)*fγtrr(rt)*(-(3/2)*(Arr/fγtrr(rt))^2 + (2/3)*fK(rt)^2
    #  - (5/2)*((X^2)/χ)/fγtrr(rt)
    #  + 2*χ/fγtθθ(rt) - 2*χ*(f∂2γtθθ(rt)/fγtθθ(rt))/fγtrr(rt)
    #  + 2*X*(f∂γtθθ(rt)/fγtθθ(rt))/fγtrr(rt)
    #  + χ*(f∂γtrr(rt)/(fγtrr(rt)^2))*(f∂γtθθ(rt)/fγtθθ(rt))
    #  - X*f∂γtrr(rt)/(fγtrr(rt)^2)
    #  + (1/2)*χ*((f∂γtθθ(rt)/fγtθθ(rt))^2)/fγtrr(rt) - 16*pi*fρ(rt,χ))
    # end
    #
    # function f∂Arr(rt,(χ, X, Arr))
    #  -fγtrr(rt)*(-(2/3)*f∂K(rt) - (3/2)*Arr*(X/χ)/fγtrr(rt)
    #  + (3/2)*Arr*(f∂γtθθ(rt)/fγtθθ(rt))/fγtrr(rt) - Arr*f∂γtrr(rt)/(fγtrr(rt)^2)
    #  - 8*pi*fγtrr(rt)*fSr(rt,χ)/χ)
    # end

    f∂M(M,rt) = 4*pi*(r(rt)^2)*fρ(M,rt)

    function f𝓗(M,∂M,rt)
         (-(3/2)*(fArr(M,∂M,rt)/fγtrr(M,rt))^2 + (2/3)*fK(M,∂M,rt)^2
         - (5/2)*((f∂χ(rt)^2)/fχ(rt))/fγtrr(M,rt) + 2*f∂2χ(rt)/fγtrr(M,rt)
         + 2*fχ(rt)/fγtθθ(rt) - 2*fχ(rt)*(f∂2γtθθ(rt)/fγtθθ(rt))/fγtrr(M,rt)
         + 2*f∂χ(rt)*(f∂γtθθ(rt)/fγtθθ(rt))/fγtrr(M,rt)
         + fχ(rt)*(f∂γtrr(M,∂M,rt)/(fγtrr(M,rt)^2))*(f∂γtθθ(rt)/fγtθθ(rt))
         - f∂χ(rt)*f∂γtrr(M,∂M,rt)/(fγtrr(M,rt)^2)
         + (1/2)*fχ(rt)*((f∂γtθθ(rt)/fγtθθ(rt))^2)/fγtrr(M,rt) - 16*pi*fρ(M,rt))
    end

    # function f𝓜r(M,∂M,rt)
    #  (f∂Arr(M,∂M,rt)/fγtrr(rt) - (2/3)*f∂K(rt) - (3/2)*fArr(rt)*(f∂χ(rt)/fχ(rt))/fγtrr(rt)
    #  + (3/2)*fArr(rt)*(f∂γtθθ(rt)/fγtθθ(rt))/fγtrr(rt) - fArr(rt)*f∂γtrr(rt)/(fγtrr(rt)^2)
    #  - 8*pi*fγtrr(rt)*fSr(rt,fχ(rt))/fχ(rt))
    # end

    # grid = Grid(domain,Int((2^num)*(n-5)+1))
    # n = grid.ncells + 4

    # ∂χ = sample(T, grid, rt->0)
    # ∂γtrr = sample(T, grid, f∂γtrr)
    # ∂γtθθ = sample(T, grid, f∂γtθθ)
    # ∂K = sample(T, grid, f∂K)
    #∂Arr = sample(T, grid, rt->(4*M/3)*(15*M^2+15*M*r(rt)+4*r(rt)^2)/real(((r(rt)^7)*((r(rt)+2*M)^3)+0im)^(1/2)))

    # ∂2γtθθ = sample(T, grid, f∂2γtθθ)
    #
    # order = 4
    #
    # rr = sample(T, grid, param[4])
    # drdrt = sample(Float64, grid, param[5])
    #
    # ∂rt𝜙 = deriv(𝜙,order,1)
    # ∂𝜙 = ∂rt𝜙./drdrt

    # fArrreg(M,∂M,rt) = real((r(rt)+ 0im)^(5/2))*fArr(M,∂M,rt)
    # fKreg(M,∂M,rt) = real((r(rt)+ 0im)^(3/2))*fK(M,∂M,rt)

    fαreg(M,rt) = 2*M
    fγtrrreg(M,rt) = 2*M
    fArrreg(M,∂M,rt) = real((r(rt)+ 0im)^(5/2))*fArr(M,∂M,rt)
    fKreg(M,∂M,rt) = real((r(rt)+ 0im)^(3/2))*fK(M,∂M,rt)

    # Constraint Equations

    rtspan = (rtspan[1], rtspan[2])
    #rtspan = (rtspan[2], 0.5)

    function constraintSystem(M, param, rt)
        f∂M(M,rt)*drdrt(rt)
    end

    # function boundaryCondition!(residual, M, param, rt)
    #     residual = M[1] - 1. #inner boundary condition
    # end

    atol = 1e-15

    BVP = ODEProblem(constraintSystem, 1., rtspan, param)
    M = solve(BVP, Tsit5(), abstol=atol, dt=drt, adaptive=false)

    ∂M(rt) = f∂M(M(rt),rt)

    # M(rt) = 1.
    # ∂M(rt) = 0.

    α = sample(T, grid, rt -> fα(M(rt),rt) )
    A = sample(T, grid, fA)
    βr = sample(T, grid, rt -> fβr(M(rt),rt) )
    Br = sample(T, grid, fBr)
    χ = sample(T, grid, fχ)
    γtrr = sample(T, grid, rt -> fγtrr(M(rt),rt) )
    γtθθreg = sample(T, grid, rt -> 0)
    Arr = sample(T, grid, rt -> fArr(M(rt),∂M(rt),rt) )
    K = sample(T, grid, rt -> fK(M(rt),∂M(rt),rt) )
    Γr = sample(T, grid, rt -> fΓr(M(rt),∂M(rt),rt))
    𝜙 = sample(T, grid, f𝜙)
    K𝜙 = sample(T, grid, fK𝜙)

    # ∂𝜙 = sample(T, grid, f∂𝜙)
    # ∂t𝜙 = sample(T, grid, f∂t𝜙)
    #
    # K𝜙 = -(∂t𝜙 - βr.*∂𝜙)./(2*α)
    #

    state = GBSSN_Variables(α, A, βr, Br, χ, γtrr, γtθθreg, Arr, K, Γr, 𝜙, K𝜙)

    #
    #cons = constraints(T,state,param)

    #println(f∂M(14.,(1., 0., 1.)))
    #println(f∂2χ(14.,(1., 0., 1.)))

    #rr = sample(T, grid, param[4])
    # 𝓗 = sample(T, grid, f𝓗)
    # 𝓜r = sample(T, grid, f𝓜r)
    #
    # ∂2χ = sample(T, grid, rt -> f∂Arr(rt,(1, 0, fArr(rt)))-f∂Arr(rt))
    #
    #println(χ[1:10].-1)
    # println(χ[n-10:n].-1)
    #
    # temp1 = sample(T, grid, rt -> solution(rt)[3]-0*fArr(rt))
    #temp = sample(T, grid, rt -> f𝓗(M(rt),∂M(rt),rt))
    #temp = sample(T, grid, rt -> M(rt))

    #plot(rr[5:n-2],temp[5:n-2])
    #plot!(rr[5:n-2],temp1[5:n-2])

    #plot(solution, vars=(0,1))
    #println(temp[1:10])
    #println(temp[n-15:n-5])

    #plot(rr[5:n-4],cons[1][5:n-4])

    return state

end

function rhs(state::GBSSN_Variables, param, t)

    ############################################
    # Caculates the right hand ride of the
    # evolved variables
    #
    # This is the main meat of the program.
    # This function contains all of the boundary
    # conditions, coordinate conversions,
    # spatial derivative calculations,
    # evolution equations, and numerical
    # dissipation. Each time the Julia DiffEq
    # Solver moves one time step, it calls
    # this function to calculate the new
    # values of the evolved variables.
    ############################################

    # Unpack the Variables

    α = state.α
    A = state.A
    βr = state.βr
    Br = state.Br
    χ = state.χ
    γtrr = state.γtrr
    γtθθreg = state.γtθθ
    Arr = state.Arr
    K = state.K
    Γr = state.Γr
    𝜙 = state.𝜙
    K𝜙 =state.K𝜙

    drt = spacing(α.grid)
    n = α.grid.ncells + 4

    m = param[7]

    # Boundary Conditions

    # These inner boundary conditions are necessary for stable
    # evolution for the specified gauge condition and do not
    # specify anything physical about the system.

    # Eulerian/Lagrangian condition (0/1)
    v = param[3]

    # if v == 0 # Eulerian Condition
    #     βr[2] = (17*βr[3] + 9*βr[4] - 5*βr[5] + βr[6])/22
    # elseif v == 1 # Lagrangian Condition
    #     γtθθreg[2] = ((-315*γtθθreg[3] + 210*γtθθreg[4]
    #     - 126*γtθθreg[5] + 45*γtθθreg[6] - 7*γtθθreg[7])/63)
    #     #βr[2] = (-315*βr[3] + 210*βr[4] - 126*βr[5] + 45*βr[6] - 7*βr[7])/63
    # end

    # Spatial Derivatives (finite differences) with respect to coordinate rt

    # Accuarcy order, 2 for 2nd order, 4 for 4th order
    order = 4

    # First derivatives
    ∂rtα = deriv(α,order,1)
    ∂rtA = deriv(A,order,1)
    ∂rtβr = deriv(βr,order,-1)
    ∂rtBr = deriv(Br,order,-1)
    ∂rtχ = deriv(χ,order,1)
    ∂rtγtrr = deriv(γtrr,order,1)
    ∂rtγtθθreg = deriv(γtθθreg,order,1)
    ∂rtArr = deriv(Arr,order,1)
    ∂rtK = deriv(K,order,1)
    ∂rtΓr = deriv(Γr,order,-1)
    ∂rt𝜙 = deriv(𝜙,order,1)
    ∂rtK𝜙 = deriv(K𝜙,order,-1)

    # Second derivatives

    # ∂2rtα = deriv2(α,order,1)
    # ∂2rtβr = deriv2(βr,order,-1)
    # ∂2rtχ = deriv2(χ,order,1)
    # ∂2rtγtrr = deriv2(γtrr,order,1)
    # ∂2rtγtθθreg = deriv2(γtθθreg,order,1)
    # ∂2rt𝜙 = deriv2(𝜙,order,1)

    ∂2rtα = deriv(∂rtα,order,1)
    ∂2rtβr = deriv(∂rtβr,order,-1)
    ∂2rtχ = deriv(∂rtχ,order,1)
    ∂2rtγtrr = deriv(∂rtγtrr,order,1)
    ∂2rtγtθθreg = deriv(∂rtγtθθreg,order,1)
    ∂2rt𝜙 = deriv(∂rt𝜙,order,1)

    # Coordinate transformations from computational rt coordinate
    # to physical r coordinate

    r = sample(Float64, α.grid, param[4])
    drdrt = sample(Float64, α.grid, param[5])
    d2rdrt = sample(Float64, α.grid, param[6])

    ∂α = ∂rtα./drdrt
    ∂A = ∂rtA./drdrt
    ∂βr = ∂rtβr./drdrt
    ∂Br = ∂rtBr./drdrt
    ∂χ = ∂rtχ./drdrt
    ∂γtrr = ∂rtγtrr./drdrt
    ∂γtθθreg = ∂rtγtθθreg./drdrt
    ∂Arr = ∂rtArr./drdrt
    ∂K = ∂rtK./drdrt
    ∂Γr = ∂rtΓr./drdrt
    ∂𝜙 = ∂rt𝜙./drdrt
    ∂K𝜙 = ∂rtK𝜙./drdrt

    ∂2α = (∂2rtα - d2rdrt.*∂α)./(drdrt.^2)
    ∂2βr = (∂2rtβr - d2rdrt.*∂βr)./(drdrt.^2)
    ∂2χ = (∂2rtχ - d2rdrt.*∂χ)./(drdrt.^2)
    ∂2γtrr = (∂2rtγtrr - d2rdrt.*∂γtrr)./(drdrt.^2)
    ∂2γtθθreg = (∂2rtγtθθreg - d2rdrt.*∂γtθθreg)./(drdrt.^2)
    ∂2𝜙 = (∂2rt𝜙 - d2rdrt.*∂𝜙)./(drdrt.^2)

    # r[1:2] .= 0.
    # βr[1:2] .= 0.
    # γtθθreg[1:2] .= 0.
    # Γr[1:2] .= 0.

    # Conversions from regularized variables to canonical variables

    γtθθ = (r.^2).*(γtθθreg .+ 1)
    ∂γtθθ = (r.^2).*∂γtθθreg + (2*r).*(γtθθreg .+ 1)
    ∂2γtθθ = (r.^2).*∂2γtθθreg + (4*r).*∂γtθθreg + 2*(γtθθreg .+ 1)

    # γtθθ[1:2] .= 0.
    # ∂γtθθ[1:2] .= 0.
    # ∂2γtθθ[1:2] .= 0.

    # K = real((r .+ 0im).^(-3/2)).*Kreg
    # ∂K = real((r .+ 0im).^(-3/2)).*∂Kreg - (3/2)*real((r .+ 0im).^(-5/2)).*Kreg
    #
    # Arr = real((r .+ 0im).^(-5/2)).*Arrreg
    # ∂Arr = real((r .+ 0im).^(-5/2)).*∂Arrreg - (5/2)*real((r .+ 0im).^(-7/2)).*Arrreg

    # Γr = -(2 ./r).*Γreg
    # ∂Γr = -(2 ./r).*∂Γreg + (2 ./(r.^2)).*Γreg

    # Γr = Γreg
    # ∂Γr = ∂Γreg

    # Gauge Conditions

    # Coordinate drift parameter.
    # Positive values lead to continued evolution
    # zero gives eventual steady state
    η = 0

    #Superscript condition...1 is a plus, 0 is a minus
    a = 0

    #Subscript condition...1 is a plus, 0 is a minus
    b = 0

    # Zero condition, 1 includes shift, 0 for vanishing shift
    # if this is 0, a and b don't matter as long as the
    # initial shift is zero
    c = 1

    #########################################################
    # Evolution Equations
    #
    # This is the full suite of evolution equations
    # for GR in spherical symmetry in the BSSN framework.
    # I have tried to keep them looking as close to their
    # mathematically written counterpart as possible.
    #
    # They are written in the order they appear in the
    # reference (arXiv:0705.3845v2) except for the ∂tBr
    # equation since it contains a ∂tΓr term.
    #
    #########################################################

    ∂tα = a*βr.*∂α - 2*α.*A

    ∂tβr = c*((3/4)*Br + b*βr.*∂βr)

    ∂tχ = ((2/3)*K.*α.*χ - (1/3)*v*βr.*χ.*∂γtrr./γtrr - (2/3)*v*βr.*χ.*∂γtθθ./γtθθ
     - (2/3)*v*χ.*∂βr + βr.*∂χ)

    ∂tγtrr = (-2*Arr.*α - (1/3)*v*βr.*∂γtrr + βr.*∂γtrr
     - (2/3)*v*γtrr.*βr.*∂γtθθ./γtθθ + 2*γtrr.*∂βr - (2/3)*v*γtrr.*∂βr)

    ∂tγtθθ = (Arr.*γtθθ.*α./γtrr - (1/3)*v*γtθθ.*βr.*∂γtrr./γtrr - (2/3)*v*βr.*∂γtθθ
     + βr.*∂γtθθ - (2/3)*v*γtθθ.*∂βr)

    ∂tArr = (-2*α.*(Arr.^2)./γtrr + K.*α.*Arr - (1/3)*v*βr.*Arr.*∂γtrr./γtrr
     - (2/3)*v*βr.*Arr.*∂γtθθ./γtθθ - (2/3)*v*Arr.*∂βr + 2*Arr.*∂βr
     + (2/3)*α.*χ.*(∂γtrr./γtrr).^2 - (1/3)*α.*χ.*(∂γtθθ./γtθθ).^2
     - (1/6)*α.*(∂χ.^2)./χ - (2/3)*α.*χ.*γtrr./γtθθ + βr.*∂Arr
     + (2/3)*α.*χ.*γtrr.*∂Γr - (1/2)*α.*χ.*(∂γtrr./γtrr).*(∂γtθθ./γtθθ)
     + (1/3)*χ.*∂γtrr.*∂α./γtrr + (1/3)*χ.*∂α.*∂γtθθ./γtθθ - (1/6)*α.*∂γtrr.*∂χ./γtrr
     - (1/6)*α.*∂γtθθ.*∂χ./γtθθ - (2/3)*∂α.*∂χ - (1/3)*α.*χ.*∂2γtrr./γtrr
     + (1/3)*α.*χ.*∂2γtθθ./γtθθ - (2/3)*χ.*∂2α + (1/3)*α.*∂2χ)

    ∂tK = ((3/2)*α.*(Arr./γtrr).^2 + (1/3)*α.*K.^2 + βr.*∂K
     + (1/2)*χ.*∂γtrr.*∂α./(γtrr.^2) - χ.*∂α.*(∂γtθθ./γtθθ)./γtrr
     + (1/2)*∂α.*∂χ./γtrr - χ.*∂2α./γtrr)

    ∂tΓr = (-v*βr.*((∂γtθθ./γtθθ).^2)./γtrr + α.*Arr.*(∂γtθθ./γtθθ)./(γtrr.^2)
     - (1/3)*v*∂βr.*(∂γtθθ./γtθθ)./γtrr + ∂βr.*(∂γtθθ./γtθθ)./γtrr
     + βr.*∂Γr + α.*Arr.*∂γtrr./(γtrr.^3)
     - (4/3)*α.*∂K./γtrr - 2*Arr.*∂α./(γtrr.^2) + (1/2)*v*∂βr.*∂γtrr./(γtrr.^2)
     - (1/2)*∂βr.*∂γtrr./(γtrr.^2) - 3*α.*Arr.*(∂χ./χ)./(γtrr.^2)
     + (1/6)*v*βr.*∂2γtrr./(γtrr.^2) + (1/3)*v*βr.*(∂2γtθθ./γtθθ)./γtrr
     + (1/3)*v*∂2βr./γtrr + ∂2βr./γtrr)

    ∂tA = ∂tK

    ∂tBr = c*(∂tΓr + b*βr.*∂Br - b*βr.*∂Γr - η*Br)

    #########################################################
    # Source Terms and Source Evolution
    #
    # This currently includes the addition of source terms
    # to GR that come from a Klein-Gordon scalar field
    #
    #########################################################

    # Klein-Gordon System

    ∂t𝜙 = βr.*∂𝜙 - 2*α.*K𝜙
    ∂tK𝜙 = (βr.*∂K𝜙 + α.*K.*K𝜙 - (1/2)*α.*χ.*∂2𝜙./γtrr
        + (1/4)*α.*χ.*∂γtrr.*∂𝜙./γtrr.^2 - (1/4)*α.*∂χ.*∂𝜙./γtrr
        - (1/2)*χ.*∂α.*∂𝜙./γtrr - (1/2)*χ.*∂γtθθ.*∂𝜙./(γtrr.*γtθθ)
        + (1/2)*∂χ.*∂𝜙./(γtrr) + (1/2)*m^2*𝜙)

    # Inverse metric (contravariant indices)

    # gtt = -(1 ./α.^2)
    # gtr = βr./α.^2
    # grr = χ./γtrr - (βr./α).^2
    # gθθ = χ./γtθθ
    #
    # # Lagrangian Density for scalar field
    #
    # 𝓛 = (1/2)*gtt.*∂t𝜙.^2 + (1/2)*grr.*∂𝜙.^2 + gtr.*∂𝜙.*∂t𝜙 - (1/2)*(m^2)*𝜙.^2
    #
    # # Stresss Energy components (contravariant indices)
    #
    # Ttt = (gtt.*∂t𝜙 + gtr.*∂𝜙).^2 - gtt.*𝓛
    # Trr = (gtr.*∂t𝜙 + grr.*∂𝜙).^2 - grr.*𝓛
    # Ttr = (gtt.*∂t𝜙 + gtr.*∂𝜙).*(gtr.*∂t𝜙 + grr.*∂𝜙) - gtr.*𝓛
    # Tθθ = -gθθ.*𝓛

    # Source Terms to GR
    # Sr here is a contravariant vector component
    # Srr here is a covariant tensor component

    # ρ = (α.^2).*Ttt
    # Sr = α.*Ttr
    # Srr = ((γtrr.^2)./(χ.^2)).*Trr
    # S = (γtrr.*Trr + 2*γtθθ.*Tθθ)./χ

    ρ = 2*K𝜙.^2 + (1/2)*(χ./γtrr).*∂𝜙.^2 + (1/2)*m^2*𝜙.^2
    Sr = 2*K𝜙.*∂𝜙
    S = 6*K𝜙.^2 - (1/2)*(χ./γtrr).*∂𝜙.^2 - (3/2)*m^2*𝜙.^2
    Srr = (γtrr./χ).*(2*K𝜙.^2 + (1/2)*(χ./γtrr).*∂𝜙.^2 - (1/2)*m^2*𝜙.^2)

    ∂tArr .+= -8*pi*α.*(χ.*Srr - (1/3)*S.*γtrr)
    ∂tK .+= 4*pi*α.*(ρ + S)
    ∂tΓr .+= -16*pi*α.*Sr./γtrr

    # Convert back to regularized variables

    # ∂tΓreg = -(r/2).*∂tΓr
    # ∂tγtθθreg = (1 ./r.^2).*∂tγtθθ

    ∂tγtθθreg = (1 ./r.^2).*∂tγtθθ
    # ∂tArrreg = real((r .+ 0im).^(5/2)).*∂tArr
    # ∂tKreg = real((r .+ 0im).^(3/2)).*∂tK

    # Numerical Dissipation terms

    ∂4α = dissipation(α)
    ∂4A = dissipation(A)
    ∂4βr = dissipation(βr)
    ∂4Br = dissipation(Br)
    ∂4χ = dissipation(χ)
    ∂4γtrr = dissipation(γtrr)
    ∂4γtθθreg = dissipation(γtθθreg)
    ∂4Arr = dissipation(Arr)
    ∂4K = dissipation(K)
    ∂4Γr = dissipation(Γr)
    ∂4𝜙 = dissipation(𝜙)
    ∂4K𝜙 = dissipation(K𝜙)

    #sign = -1 seems the best
    sign = -1
    σ = 0.2

    # ∂tα .+= (1/(2^6))*sign*σ*(drt^5)*∂6α
    # ∂tβr .+= (1/(2^6))*sign*σ*(drt^5)*∂6βr
    # ∂tBr .+= (1/(2^6))*sign*σ*(drt^5)*∂6Br
    # ∂tχ .+= (1/(2^6))*sign*σ*(drt^5)*∂6χ
    # ∂tγtrr .+= (1/(2^6))*sign*σ*(drt^5)*∂6γtrr
    # ∂tγtθθ .+= (1/(2^6))*sign*σ*(drt^5)*∂6γtθθ
    # ∂tArr .+= (1/(2^6))*sign*σ*(drt^5)*∂6Arr
    # ∂tK .+= (1/(2^6))*sign*σ*(drt^5)*∂6K
    # ∂tΓreg .+= (1/(2^6))*sign*σ*(drt^5)*∂6Γreg

    ∂tα .+= (1/16)*sign*σ*∂4α
    ∂tA .+= (1/16)*sign*σ*∂4A
    ∂tβr .+= (1/16)*sign*σ*∂4βr
    ∂tBr .+= (1/16)*sign*σ*∂4Br
    ∂tχ .+= (1/16)*sign*σ*∂4χ
    ∂tγtrr .+= (1/16)*sign*σ*∂4γtrr
    ∂tγtθθreg .+= (1/16)*sign*σ*∂4γtθθreg
    ∂tArr .+= (1/16)*sign*σ*∂4Arr
    ∂tK .+= (1/16)*sign*σ*∂4K
    ∂tΓr .+= (1/16)*sign*σ*∂4Γr
    ∂t𝜙 .+= (1/16)*sign*σ*∂4𝜙
    ∂tK𝜙 .+= (1/16)*sign*σ*∂4K𝜙

    # Inner temporal boundary Conditions

    ∂tα[1:2] .= 0
    ∂tA[1:2] .= 0
    ∂tβr[1:2] .= 0
    ∂tBr[1:2] .= 0
    ∂tχ[1:2] .= 0
    ∂tγtrr[1:2] .= 0
    ∂tγtθθreg[1:2] .= 0
    ∂tArr[1:2] .= 0
    ∂tK[1:2] .= 0
    ∂tΓr[1:2] .= 0
    ∂t𝜙[1:2] .= 0
    ∂tK𝜙[1:2] .= 0

    # ∂tα[1] = ((-25*∂tα[1]+48*∂tα[2]-36*∂tα[3]+16*∂tα[4]-3*∂tα[5])/(12*drt))
    # ∂tA[1] = ((-25*∂tα[1]+48*∂tα[2]-36*∂tα[3]+16*∂tα[4]-3*∂tα[5])/(12*drt))
    # ∂tβr[1] = ((-25*∂tα[1]+48*∂tα[2]-36*∂tα[3]+16*∂tα[4]-3*∂tα[5])/(12*drt))
    # ∂tBr[1] = ((-25*∂tα[1]+48*∂tα[2]-36*∂tα[3]+16*∂tα[4]-3*∂tα[5])/(12*drt))
    # ∂tχ[1] = ((-25*∂tα[1]+48*∂tα[2]-36*∂tα[3]+16*∂tα[4]-3*∂tα[5])/(12*drt))
    # ∂tγtrr[1] = ((-25*∂tα[1]+48*∂tα[2]-36*∂tα[3]+16*∂tα[4]-3*∂tα[5])/(12*drt))
    # ∂tγtθθreg[1] = ((-25*∂tα[1]+48*∂tα[2]-36*∂tα[3]+16*∂tα[4]-3*∂tα[5])/(12*drt))
    # ∂tArr[1] = ((-25*∂tα[1]+48*∂tα[2]-36*∂tα[3]+16*∂tα[4]-3*∂tα[5])/(12*drt))
    # ∂tK[1] = ((-25*∂tα[1]+48*∂tα[2]-36*∂tα[3]+16*∂tα[4]-3*∂tα[5])/(12*drt))
    # ∂tΓr[1] = ((-25*∂tα[1]+48*∂tα[2]-36*∂tα[3]+16*∂tα[4]-3*∂tα[5])/(12*drt))
    # ∂t𝜙[1] = ((-25*∂tα[1]+48*∂tα[2]-36*∂tα[3]+16*∂tα[4]-3*∂tα[5])/(12*drt))
    # ∂tK𝜙[1] = ((-25*∂tα[1]+48*∂tα[2]-36*∂tα[3]+16*∂tα[4]-3*∂tα[5])/(12*drt))
    #
    #
    #
    # dvalues[1] = ((-25*f.values[1] + 48*f.values[2] - 36*f.values[3]
    #     +16*f.values[4]-3*f.values[5])/(12*dx))
    #
    # dvalues[2] = ((-3*f.values[1] - 10*f.values[2] + 18*f.values[3]
    #     - 6*f.values[4] + f.values[5])/(12*dx))

    # Outer temporal boundary conditions

    ∂tα[(n-1):n] .= 0.
    ∂tA[(n-1):n] .= 0.
    ∂tβr[(n-1):n] .= 0.
    ∂tBr[(n-1):n] .= 0.
    ∂tχ[(n-1):n] .= 0.
    ∂tγtrr[(n-1):n] .= 0.
    ∂tγtθθreg[(n-1):n] .= 0.
    ∂tArr[(n-1):n] .= 0.
    ∂tK[(n-1):n] .= 0.
    ∂tΓr[(n-1):n] .= 0.
    ∂t𝜙[(n-1):n] .= 0.
    ∂tK𝜙[(n-1):n] .= 0.

    # In case you want to freeze all GR variables

    # ∂tα[2:(n-1)] .= 0.
    # ∂tβr[2:(n-1)] .= 0.
    # ∂tBr[2:(n-1)] .= 0.
    # ∂tχ[2:(n-1)] .= 0.
    # ∂tγtrr[2:(n-1)] .= 0.
    # ∂tγtθθreg[2:(n-1)] .= 0.
    # ∂tArr[2:(n-1)] .= 0.
    # ∂tK[2:(n-1)] .= 0.
    # ∂tΓreg[2:(n-1)] .= 0.

    #Values at infinity
    # α0 = 1
    # βr0 = 0
    # Br0 = 0
    # χ0 = 1
    # grr0 = 1
    # gθθreg0 = 0
    # Arr0 = 0
    # K0 = 0
    # Γreg0 = 1
    #
    # ∂tα[(n-1):n] .= (α[(n-1):n] .- α0)./r[(n-1):n] - ∂α[(n-1):n] + hα./(r[(n-1):n].^w)
    # ∂tβr[(n-1):n] .= (βr[(n-1):n] .- βr0)./r[(n-1):n] - ∂βr[(n-1):n] + hβr./(r[(n-1):n].^w)
    # ∂tBr[(n-1):n] .= (Br[(n-1):n] .- Br0)./r[(n-1):n] - ∂Br[(n-1):n] + hBr./(r[(n-1):n].^w)
    # ∂tχ[(n-1):n] .= (χ[(n-1):n] .- χ0)./r[(n-1):n] - ∂χ[(n-1):n] + hχ./(r[(n-1):n].^w)
    # ∂tgrr[(n-1):n] .= (grr[(n-1):n] .- grr0)./r[(n-1):n] - ∂grr[(n-1):n] + hgrr./(r[(n-1):n].^w)
    # ∂tgθθreg[(n-1):n] .= (gθθreg[(n-1):n] .- gθθreg0)./r[(n-1):n] - ∂gθθreg[(n-1):n] + hgθθreg./(r[(n-1):n].^w)
    # ∂tArr[(n-1):n] .= (Arr[(n-1):n] .- Arr0)./r[(n-1):n] - ∂Arr[(n-1):n] + hArr./(r[(n-1):n].^w)
    # ∂tK[(n-1):n] .= (K[(n-1):n] .- K0)./r[(n-1):n] - ∂K[(n-1):n] + hK./(r[(n-1):n].^w)
    # ∂tΓreg[(n-1):n] .= (Γreg[(n-1):n] .- Γreg0)./r[(n-1):n] - ∂Γreg[(n-1):n] + hΓreg./(r[(n-1):n].^w)

    return GBSSN_Variables(∂tα,∂tA,∂tβr,∂tBr,∂tχ,∂tγtrr,∂tγtθθreg,∂tArr,∂tK,∂tΓr,∂t𝜙,∂tK𝜙)

end

function constraints(T,state::GBSSN_Variables,param)

    ############################################
    # Caculates the constraints of the system
    #
    # These are outputed at every saved time step.
    # They should limit to zero as the spatial
    # resolution is increased. They are important to
    # monitor to make sure the physics of the system
    # is being properly modeled.
    ############################################

    # Unpack Variables

    α = state.α
    βr = state.βr
    χ = state.χ
    γtrr = state.γtrr
    γtθθreg = state.γtθθ
    Arr = state.Arr
    K = state.K
    Γr = state.Γr
    𝜙 = state.𝜙
    K𝜙 = state.K𝜙

    m = param[7]

    # Gauge conditions

    v = param[3]

    # if v == 1 # Lagrangian Condition
    #     γtθθreg[2] = ((-315*γtθθreg[3] + 210*γtθθreg[4] - 126*γtθθreg[5]
    #     + 45*γtθθreg[6] - 7*γtθθreg[7])/63)
    # end

    # Spatial Derivatives

    order = 4

    # First derivatives
    ∂rtχ = deriv(χ,order,1)
    ∂rtγtrr = deriv(γtrr,order,1)
    ∂rtγtθθreg = deriv(γtθθreg,order,1)
    ∂rtArr = deriv(Arr,order,1)
    ∂rtK = deriv(K,order,1)
    ∂rt𝜙 = deriv(𝜙,order,1)

    # Second derivatives
    ∂2rtχ = deriv2(χ,order,1)
    ∂2rtγtθθreg = deriv2(γtθθreg,order,1)

    # Coordinate transformations from computational rt coordinate
    # to physical r coordinate

    r = sample(Float64, χ.grid, param[4])
    drdrt = sample(Float64, χ.grid, param[5])
    d2rdrt = sample(Float64, χ.grid, param[6])

    ∂χ = ∂rtχ./drdrt
    ∂γtrr = ∂rtγtrr./drdrt
    ∂γtθθreg = ∂rtγtθθreg./drdrt
    ∂Arr = ∂rtArr./drdrt
    ∂K = ∂rtK./drdrt
    ∂𝜙 = ∂rt𝜙./drdrt

    ∂2χ = (∂2rtχ - d2rdrt.*∂χ)./(drdrt.^2)
    ∂2γtθθreg = (∂2rtγtθθreg - d2rdrt.*∂γtθθreg)./(drdrt.^2)

    #∂t𝜙 = βr.*∂𝜙 - 2*α.*K𝜙

    # Conversions from regularized variables to canonical variables

    γtθθ = (r.^2).*(γtθθreg .+ 1)
    ∂γtθθ = (r.^2).*∂γtθθreg + (2*r).*(γtθθreg .+ 1)
    ∂2γtθθ = (r.^2).*∂2γtθθreg + (4*r).*∂γtθθreg + 2*(γtθθreg .+ 1)

    # K = real((r .+ 0im).^(-3/2)).*Kreg
    # ∂K = real((r .+ 0im).^(-3/2)).*∂Kreg - (3/2)*real((r .+ 0im).^(-5/2)).*Kreg
    #
    # Arr = real((r .+ 0im).^(-5/2)).*Arrreg
    # ∂Arr = real((r .+ 0im).^(-5/2)).*∂Arrreg - (5/2)*real((r .+ 0im).^(-7/2)).*Arrreg

    # Γr = -(2 ./r).*Γreg

    # Inverse metric (contravariant indices)

    # gtt = -(1 ./α.^2)
    # gtr = βr./α.^2
    # grr = χ./γtrr - (βr./α).^2
    # gθθ = χ./γtθθ
    #
    # # Lagrangian Density for scalar field
    #
    # 𝓛 = (1/2)*gtt.*∂t𝜙.^2 + (1/2)*grr.*∂𝜙.^2 + gtr.*∂𝜙.*∂t𝜙 - (1/2)*(m^2)*𝜙.^2
    #
    # # Stresss Energy components (contravariant indices)
    #
    # Ttt = (gtt.*∂t𝜙 + gtr.*∂𝜙).^2 - gtt.*𝓛
    # Ttr = (gtt.*∂t𝜙 + gtr.*∂𝜙).*(gtr.*∂t𝜙 + grr.*∂𝜙) - gtr.*𝓛
    #
    # # Source Terms to the constraints
    #
    # ρ = (α.^2).*Ttt
    # Sr = α.*Ttr

    ρ = 2*K𝜙.^2 + (1/2)*(χ./γtrr).*∂𝜙.^2 + (1/2)*m^2*𝜙.^2

    Sr = 2*K𝜙.*∂𝜙

    # Constraint Equations

    𝓗 = (-(3/2)*(Arr./γtrr).^2 + (2/3)*K.^2 - (5/2)*((∂χ.^2)./χ)./γtrr
     + 2*∂2χ./γtrr + 2*χ./γtθθ - 2*χ.*(∂2γtθθ./γtθθ)./γtrr + 2*∂χ.*(∂γtθθ./γtθθ)./γtrr
     + χ.*(∂γtrr./(γtrr.^2)).*(∂γtθθ./γtθθ) - ∂χ.*∂γtrr./(γtrr.^2)
     + (1/2)*χ.*((∂γtθθ./γtθθ).^2)./γtrr - 16*pi*ρ)

    𝓜r = (∂Arr./γtrr - (2/3)*∂K - (3/2)*Arr.*(∂χ./χ)./γtrr
     + (3/2)*Arr.*(∂γtθθ./γtθθ)./γtrr - Arr.*∂γtrr./(γtrr.^2)
     - 8*pi*Sr)

    𝓖r = -(1/2)*∂γtrr./(γtrr.^2) + Γr + (∂γtθθ./γtθθ)./γtrr

    𝓗[1:2] .= 0.
    𝓜r[1:2] .= 0.
    𝓖r[1:2] .= 0.

    return (𝓗, 𝓜r, 𝓖r, ρ)

end

function horizon(T,state::GBSSN_Variables,param)

    ############################################
    # Caculates the apparent horizon
    #
    # Where the function crosses zero is the
    # apparent horizon of the black hole.
    ############################################

    v = param[3]

    # Unpack Variables

    χ = state.χ
    γtrr = state.γtrr
    γtθθreg = state.γtθθ
    Arr = state.Arr
    K = state.K

    # Gauge condition

    # if v == 1 # Lagrangian Condition
    #     γtθθreg[2] = ((-315*γtθθreg[3] + 210*γtθθreg[4] - 126*γtθθreg[5]
    #     + 45*γtθθreg[6] - 7*γtθθreg[7])/63)
    # end

    r = sample(T, χ.grid, param[4])
    drdrt = sample(T, χ.grid, param[5])

    # Conversions from regularized variables to canonical variables

    γtθθ = (r.^2).*(γtθθreg .+ 1)

    # K = real((r .+ 0im).^(-3/2)).*Kreg
    #
    # Arr = real((r .+ 0im).^(-5/2)).*Arrreg

    # Intermediate calculations

    Kθθ = ((1/3)*γtθθ.*K - (1/2)*Arr.*γtθθ./γtrr)./χ

    grr =  γtrr./χ

    gθθ =  γtθθ./χ

    # Spatial Derivatives

    ∂rtgθθ = deriv(gθθ,4,1)

    # Coordinate transformations from computational rt coordinate
    # to physical r coordinate

    ∂gθθ = ∂rtgθθ./drdrt

    # Apparent horizon function

    Θ = (∂gθθ./gθθ)./real((grr .+ 0im).^(1/2)) - 2*Kθθ./gθθ

    # cross = GridFun(χ.grid, sign.(Θ))

    return Θ

end

# function crossings(Θ::GridFun)
#
#     GridFun(χ.grid, sign.(Θ))
#
#     for
#
# end

function custom_progress_message(dt,state,param,t)

    ###############################################
    # Outputs status numbers while the program runs
    ###############################################

    if param[1]==param[2]
        println("")
        println("| # | Time Step | Time | max α'(t) | max χ'(t) | max γtrr'(t) | max γtθθ'(t) | max Arr'(t) | max K'(t) | max Γr'(t) |")
        println("|___|___________|______|___________|___________|______________|______________|_____________|___________|____________|")
        println("")
    end

    derivstate = rhs(state,param,t)

    println("  ",
    rpad(string(param[1]),6," "),
    rpad(string(round(dt,digits=3)),10," "),
    rpad(string(round(t,digits=3)),10," "),
    rpad(string(round(maximum(abs.(derivstate.α)),digits=3)),12," "),
    rpad(string(round(maximum(abs.(derivstate.χ)),digits=3)),12," "),
    rpad(string(round(maximum(abs.(derivstate.γtrr)),digits=3)),14," "),
    rpad(string(round(maximum(abs.(derivstate.γtθθ)),digits=3)),14," "),
    rpad(string(round(maximum(abs.(derivstate.Arr)),digits=3)),12," "),
    # rpad(string(round(maximum(abs.(derivstate.𝜙)),digits=3)),12," "),
    # rpad(string(round(maximum(abs.(derivstate.K𝜙)),digits=3)),14," ")
    rpad(string(round(maximum(abs.(derivstate.K)),digits=3)),12," "),
    rpad(string(round(maximum(abs.(derivstate.Γr)),digits=3)),14," ")
    )

    #PrettyTables.jl

    param[1] += param[2]

end

function solution_saver(T,grid,sol,param,folder)

    ###############################################
    # Saves all of the variables in nice CSV files
    # in the choosen data folder directory
    ###############################################

    vars = (["α","A","βr","Br","χ","γtrr","γtθθ","Arr","K","Γr","𝜙","K𝜙",
    "∂tα","∂tA","∂tβr","∂tBr","∂tχ","∂tγtrr","∂tγtθθ","∂tArr","∂tK","∂tΓr","∂t𝜙","∂tK𝜙",
    "H","Mr","Gr","ρ","appHorizon"])
    varlen = length(vars)
    #mkdir(string("data\\",folder))
    tlen = size(sol)[3]
    rlen = grid.ncells + 4
    loc = sample(T, grid, param[4])
    #loc[1] =
    cons = Array{GridFun,2}(undef,tlen,4)
    derivs = Array{GBSSN_Variables,1}(undef,tlen)
    apphorizon = Array{GridFun,1}(undef,tlen)

    for i in 1:tlen
        derivs[i] = rhs(sol[i],param,0)
        cons[i,1:4] .= constraints(T,sol[i],param)
        apphorizon[i] = horizon(T,sol[i],param)
    end

    array = Array{T,2}(undef,tlen+1,rlen+1)

    array[1,1] = 0
    array[1,2:end] .= loc

    for j = 1:varlen
        if j < 13
            for i = 2:tlen+1
                array[i,1] = sol.t[i-1]
                array[i,2:end] .= sol[:,j,i-1]
            end
        elseif j < 25
            for i = 2:tlen+1
                array[i,1] = sol.t[i-1]
                array[i,2:end] .= derivs[i-1][:,j-12]
            end
        elseif j < 29
            for i = 2:tlen+1
                array[i,1] = sol.t[i-1]
                array[i,2:end] .= cons[i-1,j-24]
            end
        elseif j == varlen
            for i = 2:tlen+1
                array[i,1] = sol.t[i-1]
                array[i,2:end] .= apphorizon[i-1]
            end
        end

        CSV.write(
        string("data/",folder,"/",vars[j],".csv"),
        DataFrame(array, :auto),
        header=false
        )

    end

end

function main(points)

    ###############################################
    # Main Program
    #
    # Calls each of the above functions to run a
    # simulation. Sets up the numerical grid,
    # sets the gauge conditions, sets up a new
    # computational rt coordinate that makes the
    # physical r coordinate step larger near the
    # outer boundary, sets the initial conditions,
    # and finally runs the numerical DiffEq
    # package to run the time integration.
    #
    # All data is saved in the folder specified to
    # the solution_saver, each in their own CSV
    # file.
    ###############################################

    T = Float64
    rspan = T[1,210]
    rtspan = T[1,21]

    grid = setup(T, rtspan, points)
    drt = spacing(grid)
    dt = drt/4

    tspan = T[0,25]
    v = 1

    m = 0

    f(b) = b*tan(rtspan[2]/b)-rspan[2]

    scale = find_zero(f, 0.64*rtspan[2])

    r(rt) = scale*tan(rt/scale)
    drdrt(rt) = sec(rt/scale)^2
    d2rdrt(rt) = (2/scale)*(sec(rt/scale)^2)*tan(rt/scale)

    # r(rt) = rt
    # drdrt(rt) = 1
    # d2rdrt(rt) = 0

    atol = eps(T)^(T(3) / 4)

    #alg = KuttaPRK2p5()
    alg = RK4()

    #printlogo()

    printtimes = 1
    custom_progress_step = round(Int, printtimes/dt)
    step_iterator = custom_progress_step
    param = [step_iterator, custom_progress_step, v, r, drdrt, d2rdrt, m, rtspan]
    println("Defining Initial State...")
    state = init(T, grid, param)::GBSSN_Variables
    println("Defining Problem...")
    prob = ODEProblem(rhs, state, tspan, param)
    println("Starting Solution...")

    #init(T, grid, param)

    sol = solve(
        prob, alg,
        abstol = atol,
        dt = drt/4,
        adaptive = false,
        saveat = 1,
        progress = true,
        progress_steps=custom_progress_step,
        progress_message=custom_progress_message
    )

    solution_saver(T,grid,sol,param,"ScalarTests_old")

end

end


text:
mutable struct TileSelection
    fg::Bool
    tiles::Array{Char, 2}
    selection::Rectangle
    startX::Number
    startY::Number
    offsetX::Number
    offsetY::Number
end

TileSelection(fg::Bool, tiles::Array{Char, 2}, selection::Rectangle) = TileSelection(fg, tiles, selection, selection.x, selection.y, 0, 0)

selectableLayers = ["fgTiles", "bgTiles", "entities", "triggers", "fgDecals", "bgDecals"]
selectionTargets = Dict{String, Function}(
    "entities" => room -> room.entities,
    "triggers" => room -> room.triggers,
    "bgDecals" => room -> room.bgDecals,
    "fgDecals" => room -> room.fgDecals
)

function getSelection(trigger::Maple.Trigger, room::Maple.Room, node::Number=0)
    return triggerSelection(trigger, room, node)
end

function getSelection(decal::Maple.Decal, room::Maple.Room, node::Number=0)
    return decalSelection(decal)
end

function getSelection(entity::Maple.Entity, room::Maple.Room, node::Number=0)
    return entitySelection(entity, room, node)
end

function getSelection(target::TileSelection, room::Maple.Room, node::Number=0)
    return target.selection
end

function position(target::TileSelection)
    return target.startX, target.startY
end

# TODO - Use mouse position and check if its closer to the center as well
# Area is "good enough" for now
function bestSelection(set::Set{Tuple{String, Rectangle, Any, Number}})
    best = nothing
    bestVal = typemax(Int)

    for selection in set
        layer, rect, target, node = selection
        area = rect.w * rect.h

        if area < bestVal
            best = selection
            bestVal = area
        end
    end

    return best
end

function getSelected(room::Room, name::String, selection::Rectangle)
    res = Set{Tuple{String, Rectangle, Any, Number}}()

    # Rectangular based selection - Triggers, Entities, Decals
    if haskey(selectionTargets, name)
        targets = selectionTargets[name](room)

        for target in targets
            rect = getSelection(target, room)

            if isa(rect, Rectangle)
                if checkCollision(selection, rect)
                    push!(res, (name, rect, target, 0))
                end

            elseif isa(rect, Array{Rectangle, 1})
                for (i, r) in enumerate(rect)
                    if checkCollision(selection, r)
                        # The first rect is the main entity itself, followed by the nodes
                        push!(res, (name, r, target, i - 1))
                    end
                end
            end
        end

    # Tile based selections
    elseif name == "fgTiles" || name == "bgTiles"
        fg = name == "fgTiles"
        tiles = fg ? room.fgTiles.data : room.bgTiles.data

        tx, ty = floor(Int, selection.x / 8), floor(Int, selection.y / 8)
        tw, th = ceil(Int, selection.w / 8), ceil(Int, selection.h / 8)

        gx, gy = tx * 8, ty * 8
        gw, gh = tw * 8, th * 8
        gridSelection = Rectangle(gx, gy, gw, gh)

        drawingTiles = fill('0', (th + 2, tw + 2))
        drawingTiles[2:end - 1, 2:end - 1] = get(tiles, (ty + 1:ty + th, tx + 1:tx + tw), '0')

        target = TileSelection(
            fg,
            drawingTiles,
            gridSelection
        )

        push!(res, (name, gridSelection, target, 0))
    end

    return res
end

getSelected(room::Room, layer::Layer, selection::Rectangle) = getSelected(room, layerName(layer), selection)

# Currently no support for tiles
function getFilteredSelections(selections::Set{Tuple{String, Rectangle, Any, Number}}, rect::Rectangle, room::Room, name::String, func::Function)
    res = Set{Tuple{String, Rectangle, Any, Number}}()

    best = hasSelectionAt(selections, rect, room)
    if best == false
        return res
    end

    # Rectangular based selection - Triggers, Entities, Decals
    if haskey(selectionTargets, name)
        targets = selectionTargets[name](room)

        for target in targets
            rect = getSelection(target, room)

            if isa(rect, Rectangle)
                if func(best, target)
                    push!(res, (name, rect, target, 0))
                end

            elseif isa(rect, Array{Rectangle, 1})
                for (i, r) in enumerate(rect)
                    if func(best, target)
                        # The first rect is the main entity itself, followed by the nodes
                        push!(res, (name, r, target, i - 1))
                    end
                end
            end
        end
    end

    return res
end

function getFilteredSelections(selections::Set{Tuple{String, Rectangle, Any, Number}}, rect::Rectangle, room::Room, layer::Layer, func::Function)
    return getFilteredSelections(selections, rect, room, layerName(layer), func)
end

# Return a function that can compare the base target with all the other targets
# Returning true if they are "same-ish"
function getFilterFunction(name::String, strict::Bool=false)
    if name == "entities" || name == "triggers"
        # Strict mode is all flags the same, excluding x, y and nodes
        # Normal mode is just name

        if strict
            return function(base::Union{Entity, Trigger}, target::Union{Entity, Trigger})
                # Don't count nodes as a key
                if (length(base.data) - haskey(base.data, "nodes")) != (length(target.data) - haskey(target.data, "nodes"))
                    return false
                end

                if base.name != target.name
                    return false
                end

                for (k, v) in base.data
                    if k == "x" || k == "y" || k == "nodes"
                        continue
                    end

                    if get(target.data, k, nothing) != v
                        return false
                    end
                end

                return true
            end

        else
            return function(base::Union{Entity, Trigger}, target::Union{Entity, Trigger})
                return base.name == target.name
            end
        end

    elseif name == "bgDecals" || name == "fgDecals"
        # Strict mode is same scales as well as texture
        # Normal mode is just texture

        if strict
            return function(base::Decal, target::Decal)
                return base.texture == target.texture && base.scaleX == target.scaleX && base.scaleY == target.scaleY
            end

        else
            return function(base::Decal, target::Decal)
                return base.texture == target.texture
            end
        end
    end

    return function(base, target)
        return false
    end
end

getFilterFunction(layer::Layer, strict::Bool=false) = getFilterFunction(layerName(layer), strict)

function hasSelectionAt(selections::Set{Tuple{String, Rectangle, Any, Number}}, rect::Rectangle, room::Maple.Room)
    for (layer, box, target, node) in selections
        targetRect = getSelection(target, room)
        if isa(targetRect, Rectangle) && checkCollision(rect, targetRect) && node == 0
            return target

        elseif isa(targetRect, Array{Rectangle, 1})
            for (i, r) in enumerate(targetRect)
                if checkCollision(rect, r) && node == i - 1
                    return target
                end
            end
        end
    end

    return false
end

function updateSelections!(selections::Set{Tuple{String, Rectangle, Any, Number}}, room::Room, name::String, rect::Rectangle; retain::Bool=false, best::Bool=false, mass::Bool=false)
    # Return selections that are no longer selected
    unselected = Set{Tuple{String, Rectangle, Any, Number}}()
    newlySelected = Set{Tuple{String, Rectangle, Any, Number}}()

    # Holding shift keeps the last selection as well
    if !retain
        unselected = deepcopy(selections)
        empty!(selections)
    end

    validSelections = Set{Tuple{String, Rectangle, Any, Number}}()
    if name == "all"
        for n in selectableLayers
            if mass
                func = getFilterFunction(n, best)
                union!(validSelections, getFilteredSelections(retain ? selections : unselected, rect, room, n, func))

            else
                union!(validSelections, getSelected(room, n, rect))
            end
        end

    else
        if mass
            func = getFilterFunction(name, best)
            union!(validSelections, getFilteredSelections(retain ? selections : unselected, rect, room, name, func))

        else
            union!(validSelections, getSelected(room, name, rect))
        end
    end

    if best && !mass
        target = bestSelection(validSelections)
        if target !== nothing
            push!(selections, target)
            push!(newlySelected, target)
        end

    else
        union!(selections, validSelections)
        union!(newlySelected, validSelections)
    end

    return unselected, newlySelected
end

function fixSelections(room::Maple.Room, selections::Set{Tuple{String, Rectangle, Any, Number}})
    res = Set{Tuple{String, Rectangle, Any, Number}}()

    for (layer, box, target, node) in selections
        if haskey(selectionTargets, layer)
            targets = selectionTargets[layer](room)

            index = findfirst(isequal(target), targets)
            if index !== nothing
                push!(res, (layer, box, targets[index], node))
            end

        else
            push!(res, (layer, box, target, node))
        end
    end

    return res
end

updateSelections!(selections::Set{Tuple{String, Rectangle, Any, Number}}, room::Room, layer::Layer, rect::Rectangle) = updateSelections!(selections, room, layerName(layer), rect)

# TODO - Verify that scales are respected properly
function getSpriteRectangle(texture::String, x::Number, y::Number; jx::Number=0.5, jy::Number=0.5, sx::Number=1, sy::Number=1)
    sprite = getTextureSprite(texture)

    width, height = sprite.width, sprite.height
    realWidth, realHeight = sprite.realWidth, sprite.realHeight

    drawX, drawY = floor(Int, x - (realWidth * jx + sprite.offsetX) * sx), floor(Int, y - (realHeight * jy + sprite.offsetY) * sy)
    drawX += sx < 0 ? width * sx : 0
    drawY += sy < 0 ? height * sy : 0

    return Rectangle(drawX, drawY, sprite.width * abs(sx), sprite.height * abs(sy))
end

function getEntityRectangle(entity::Maple.Entity)
    x, y = position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return Rectangle(x, y, width, height)
end


text:
module NonMeansTestedBenefits
#
# This module provides routines for benefits which are (mostly) not means-tested benefits,
# both current and historic. See [HistoricBenefits.jl] for the imputation routines used here.
#
    using Base: has_tight_type
    using Dates
    using Dates: Date, now, TimeType, Year
    using ScottishTaxBenefitModel

    using .Definitions
    using .STBIncomes

    using .Utils: nearest

    using .ModelHousehold:
        BenefitUnit,
        Household,
        Person,
        get_benefit_units,
        has_income,
        num_children

    using .BenefitGenerosity: change_status

    using .Intermediate:
        has_limited_capactity_for_work_activity,
        has_limited_capactity_for_work,
        make_recipient

    using .STBParameters:
        AgeLimits,
        AttendanceAllowance,
        BereavementSupport,
        CarersAllowance,
        ChildBenefit,
        ContributoryESA,
        DisabilityLivingAllowance,
        HoursLimits,
        JobSeekersAllowance,
        MaternityAllowance,
        NonMeansTestedSys,
        PersonalIndependencePayment,
        RetirementPension,
        WidowsPensions,
        reached_state_pension_age

    using .Results:
        BenefitUnitResult,
        HouseholdResult,
        IndividualResult,
        LMTIncomes

    export
        calc_child_benefit!,
        calc_dla,
        calc_pip,
        calc_post_tax_non_means_tested!,
        calc_pre_tax_non_means_tested!,
        calc_state_pension,
        calc_widows_benefits

    """
    Child Benefit - this has to be done *after* income tax, so we have
    income tax `total_income` for each adult in the `BenefitUnitResult` struct.
    """
    function calc_child_benefit!(
        bures :: BenefitUnitResult,
        bu    :: BenefitUnit,
        cb    :: ChildBenefit{T} ) :: T where T
        c = zero(T)
        #= FIXME something to check the exact type of children -
        for now, assume the BU allocation has got this right.
        for pid in bu.children
            #
        end
        =#
        if cb.abolished
            return c
        end
        nc = num_children(bu)
        if nc == 0
            return 0.0
        end
        c += cb.first_child
        if nc > 1
            c += (nc-1)*cb.other_children
        end
        recipient = make_recipient( bu, CHILD_BENEFIT )
        # guardian's allowance
        bures.pers[recipient].income[GUARDIANS_ALLOWANCE] = 0.0

        # println( bu.people[bu.head].relationships)

        if c > 0 # fixme not quite right - qualify for CB for each child
            for cp in bu.children
                # this checks if anyone in the BU has a parent-like
                # relationship to the child
                is_guardian = true
                for pid in bu.adults
                    if ! (bu.people[pid].relationships[cp] in [Grand_parent, Other_relative, Other_non_relative])
                        is_guardian = false
                        break;
                    end
                end
                if is_guardian
                    # FIXME we could also guess approx number of kids
                    # qualifying for GA from receipt in HistoricBenefits || has_income( bu, guardians_allowance )
                    bures.pers[recipient].income[GUARDIANS_ALLOWANCE] +=
                        cb.guardians_allowance
                end
            end
        end # guardian's allowance
        # cb but not guardians withdrawn with high incomes
        # high income thing - highest of the BU's *individual* income; see cpag ch.27
        max_inc = zero(T)
        for pid in bu.adults
            max_inc = max( max_inc, bures.pers[pid].it.total_income )
            if bu.people[pid].sex == Female
                recipient = pid
            end
        end
        if max_inc > cb.high_income_thresh
            # this should really be done in steps £1 for every £100, but since
            # everything is weekly here we'll just multiply
            withdrawn = cb.withdrawal*(max_inc - cb.high_income_thresh)
            c = max(0.0, c-withdrawn)
        end
        bures.pers[recipient].income[CHILD_BENEFIT] = c
    end

    function calc_widows_benefits(
        pers :: Person{T},
        has_kids :: Bool,
        bp :: BereavementSupport{T},
        wp :: WidowsPensions{T}) :: T where T
        wid = zero(T)
        if wp.abolished || bp.abolished
            return wid
        end
        #
        # We don't know when someone was widowed
        # so we rely on this. FIXME Obviously this gets progressively worse as time goes on and
        # we keep using the old years of FRS.
        #
        # new-style: payable for 18 months at a flat rate
        # CHECK 2829 2017 why is that bereavment_support?
        if pers.bereavement_type in [bereavement_allowance,bereavement_support]
            # payable for 18 months so we'll allocate
            # 2/3rds of the weekly value of the lump-sum
            if has_kids
                wid = bp.higher + bp.lump_sum_higher*2/3
            else
                wid = bp.lower + bp.lump_sum_lower*2/3
            end
        elseif pers.bereavement_type == widowed_parents # At least 3 years ago so no
            # need to worry about lump-sums; just scale
            # the standard rate by the ratio of
            # their receipt to the standard rate at the time
            # of interview. CHECK hhld 477 2015 for an example
            # of someone with the lump sum but no standard benefit
            # so no ratio; for us, just assign the standard
            # rate in those cases.
            #
            # println( "pers.benefit_ratios $(pers.benefit_ratios)")
            if haskey( pers.benefit_ratios, bereavement_allowance_or_widowed_parents_allowance_or_bereavement )
                wid = pers.benefit_ratios[bereavement_allowance_or_widowed_parents_allowance_or_bereavement]*
                    wp.standard_rate
            else
                wid = wp.standard_rate
            end
        else
            # check we're not missing anyone
            @assert ! haskey( pers.benefit_ratios, bereavement_allowance_or_widowed_parents_allowance_or_bereavement)
        end
        return wid
    end

    """
    PIP and DLA (below) rely on all the types being sorted out earlier
    either in some kind of probit or by being inferred from receipts (see [HistoricBenefits.jl] for
    what we have so far).
    """
    function calc_pip(
        pers :: Person{T},
        pip  :: PersonalIndependencePayment{T}) :: Tuple{T,T} where T
        pl = zero(T)
        pm = zero(T)
        daily_type = pers.pip_daily_living_type
        mob_type = pers.pip_mobility_type
        if pip.abolished
            return (pl, pm )
        end
        daily_type = change_status(
            candidates=pip.dl_candidates,
            pid=pers.pid,
            change=pip.extra_people,
            choices=[standard_pip,enhanced_pip],
            current_value=daily_type,
            disqual_value=no_pip )
        mob_type = change_status(
            candidates=pip.mobility_candidates,
            pid=pers.pid,
            choices=[standard_pip,enhanced_pip],
            change=pip.extra_people,
            current_value=mob_type,
            disqual_value=no_pip )

        # fixme consistent names dl->daily living etc.
        if daily_type == standard_pip
            pl = pip.dl_standard
        elseif daily_type == enhanced_pip
            pl = pip.dl_enhanced
        end
        if mob_type == standard_pip
            pm = pip.mobility_standard
        elseif mob_type == enhanced_pip
            pm = pip.mobility_enhanced
        end
        return (pl, pm )
    end # pip calc

    function calc_attendance_allowance(
        pers :: Person{T},
        aa  :: AttendanceAllowance{T},
         ) :: T where T
        a =zero(T)
        if aa.abolished
            return a
        end
        at = pers.attendance_allowance_type
        at = change_status(
            candidates=aa.candidates,
            pid=pers.pid,
            choices=[high,low],
            change=aa.extra_people,
            current_value=at,
            disqual_value=missing_lmh  )
        if at == missing_lmh
            a = zero(T)
        elseif at == high
            a = aa.higher
        else
            a = aa.lower;
        end
        return a
    end

    function calc_dla(
        pers :: Person{T},
        dla  :: DisabilityLivingAllowance{T} ) :: Tuple{T,T} where T
        dc = zero(T)
        dm = zero(T)
        if dla.abolished
            return (dc,dm)
        end
        dla_s = pers.dla_self_care_type
        dla_m = pers.dla_mobility_type
        # FIXME we use the same list for both mob and self
        # I think because of small sample size (kids only)
        dla_s = change_status(
            candidates=dla.candidates,
            pid=pers.pid,
            choices=[high,mid,low],
            change=dla.extra_people,
            current_value=dla_s,
            disqual_value=missing_lmh  )
        dla_m = change_status(
            candidates=dla.candidates,
            pid=pers.pid,
            choices=[high,mid,low],
            change=dla.extra_people,
            current_value=dla_m,
            disqual_value=missing_lmh  )
        # FIXME make all these names constisent (mid/middle,care->self_care etc.)
        if dla_s == high
            dc = dla.care_high
        elseif dla_s == mid
            dc = dla.care_middle
        elseif dla_s == low
            dc = dla.care_low
        end
        if dla_m == high
            dm = dla.mob_high
        elseif dla_m in (low,mid)
            dm = dla.mob_low
        end
        # println( "setting DLA as $dc, $dm")
        return (dc,dm);
    end # dla calc

    function calc_state_pension(
        pers :: Person{T},
        rp :: RetirementPension{T},
        age_limits :: AgeLimits ) :: T where T
        pen = zero(T)
        if rp.abolished
            return pen
        end
        if reached_state_pension_age(
            age_limits,
            pers.age,
            pers.sex )
            # so, was someone over pension age *before* the new state pension?
            # if so, use the proportion of the `class_a` they seemed to be on at
            # the time, times the current `class_a`.
            if reached_state_pension_age(
                # old style
                age_limits,
                pers.age,
                pers.sex,
                age_limits.savings_credit_to_new_state_pension )
                # println( "pers.benefit_ratios $(pers.benefit_ratios)")
                if ! haskey( pers.benefit_ratios, state_pension )
                    ratio =  pers.age < 70 ? 0.0 : 1.0 # kinda crude deferrment
                else
                    ratio = pers.benefit_ratios[state_pension]
                end
                pen = ratio*rp.cat_a
            else
                #
                # new style
                pen = rp.new_state_pension
            end
        end
        return pen
    end

    """
    FIXME Model this properly
    """
    function calc_esa(
        pers :: Person{T},
        esa  :: ContributoryESA{T}) :: T where T
        e = zero(T)
        if esa.abolished
            return e
        end
        if pers.esa_type == contributory_jsa
            if pers.age < 25
                # FIXME not quite right since
                # could be past assessment stage;
                # maybe check time out of work?
                e = esa.assessment_u25
            else
                e = esa.main
            end
            if has_limited_capactity_for_work_activity( pers )
                e += esa.support
            end
        end
        return e
    end

    function calc_maternity_allowance(
        pers :: Person{T},
        ma :: MaternityAllowance ) :: T where T
        m = zero(T)
        if ma.abolished
            return m
        end
        # fixme the design means you should never have to check the incomes dict here
        if has_income( pers, maternity_allowance )
            m = ma.rate
        end
        return m
    end

    function calc_carers_allowance(
        pers :: Person{T},
        pres :: IndividualResult{T},
        carers :: CarersAllowance{T}) :: T where T
        c = zero(T)
        if carers.abolished
            return c
        end
        earnings :: T = isum(
            pres.income,
            carers.earnings;
            deducted=carers.deductions )
        # println( "earnings=$earnings carers.gainful_employment_min=$(carers.gainful_employment_min)")
        if pers.hours_of_care_given >= carers.hours &&
            earnings < carers.gainful_employment_min
            c = carers.allowance
        end
        return c
    end

    function calc_jsa(
        pers :: Person{T},
        jsa :: JobSeekersAllowance{T},
        hrs  :: HoursLimits ) :: T where T
        j = zero(T)
        if jsa.abolished
            return j
        end
        if pers.jsa_type == contributory_jsa &&
           pers.usual_hours_worked <= hrs.med &&
           (! has_limited_capactity_for_work( pers ))
            j = pers.age < 25 ? jsa.u25 : jsa.o24
        end
        return j
    end

    """
    Household level calculations for all nmt benefits that don't require knowlege
    of income tax/ni liabilties - not necessarily taxable benefits, just things
    that don't require any kind of net income calculation and so can be done
    before IT/NI.
    """
    function calc_pre_tax_non_means_tested!(
        hhres :: HouseholdResult,
        hh    :: Household,
        sys   :: NonMeansTestedSys,
        hours_limits   :: HoursLimits,
        age_limits :: AgeLimits )
        ## maybe add a benefit unit allocator
        bus = get_benefit_units( hh )
        buno = 1
        for bu in bus
            has_children = size( bu.children )[1] > 0
            bures = hhres.bus[buno]
            for adno in bu.adults
                pers = bu.people[adno]
                pres = bures.pers[adno]
                pres.income[STATE_PENSION] = calc_state_pension(
                    pers,
                    sys.pensions, age_limits );
                pres.income[WIDOWS_PAYMENT] = calc_widows_benefits(
                    pers, has_children, sys.bereavement, sys.widows_pension )
                #
                # FIXME
                # pip/dla can only be claimed
                # if ! reached_state_pension_age( age_limits, pers.age, pers.sex )
                # but claims can run on indefinitely and for now we're just using
                # receipts, so ignore any upper age limits until we model these fully.
                #
                pres.income[sys.dla.care_slot],
                pres.income[sys.dla.mob_slot] = calc_dla( pers, sys.dla );
                pres.income[sys.pip.care_slot],
                pres.income[sys.pip.mob_slot] = calc_pip( pers, sys.pip )

                #
                # .. conversely, this age limit seems safe: a 62 yo female recieving
                # in the data should be disallowed now the pension age has increased.
                #
                if reached_state_pension_age( age_limits, pers.age, pers.sex )
                    pres.income[sys.attendance_allowance.slot] = calc_attendance_allowance( pers, sys.attendance_allowance )
                else
                    pres.income[CONTRIB_EMPLOYMENT_AND_SUPPORT_ALLOWANCE] = calc_esa( pers, sys.esa )
                    pres.income[CONTRIB_JOBSEEKERS_ALLOWANCE] = calc_jsa( pers, sys.jsa, hours_limits )
                    pres.income[MATERNITY_ALLOWANCE] = calc_maternity_allowance( pers, sys.maternity )
                end
                # NON-overlapping rules p1178 go here
            end # ad loop
            buno += 1
        end # bu loop
    end # calc_non_means_tested

    """
    NMT Benefits that require knowlege of tax and NI liabilities and so have to be done
    after a tax calculation - not necessarily tax free as such (CB higher charge). Kinda-sorta
    means-tested bens, I suppose.
    """
    function calc_post_tax_non_means_tested!(
        hhres :: HouseholdResult,
        hh    :: Household,
        sys   :: NonMeansTestedSys,
        age_limits :: AgeLimits )
        ## maybe add a benefit unit allocator
        bus = get_benefit_units( hh )
        buno = 1
        for bu in bus
            has_children = size( bu.children )[1] > 0
            bures = hhres.bus[buno]
            calc_child_benefit!(
                bures,
                bu,
                sys.child_benefit )
            for adno in bu.adults
                pers = bu.people[adno]
                pres = bures.pers[adno]
                pres.income[sys.carers.slot] =
                    calc_carers_allowance( pers, pres, sys.carers )

                if hh.region == Scotland
                    if pres.income[sys.carers.slot] > 0
                        pres.income[SCOTTISH_CARERS_SUPPLEMENT] = sys.carers.scottish_supplement
                    end
                end
                # NON-overlapping rules p1178 go here
            end # ad loop
            buno += 1
        end # bu loop
    end # calc_non_means_tested

end # package non-means-tested

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
using ..ChemistryFeaturization.AbstractType: AbstractCodec

"""
    OneHotOneCold(encode_f, decode_f, nbins, logspaced)

AbstractCodec type which uses a dummy variable (as defined in statistical literature), i.e., which employs
one-hot encoding and a one-cold decoding scheme.
"""
struct OneHotOneCold <: AbstractCodec
    encode_f::Function
    decode_f::Function
    nbins::Integer
    logspaced::Bool
end

"A flexible version of Flux.onehot that can handle both categorical and continuous-valued encoding."
function build_onehot_vec(val, bins, categorical)
    local bin_index, onehot_vec
    if categorical
        onehot_vec = [0.0 for i = 1:length(bins)]
        bin_index = findfirst(isequal(val), bins)
    else
        @assert eltype(bins) <: Number "Your bins aren't numbers...are you sure you didn't mean for this feature to be categorical?"
        @assert bins[1] <= val <= bins[end] "The value $val is outside the range of bins $bins"
        onehot_vec = [0.0 for i = 1:(length(bins)-1)]
        bin_index = searchsorted(bins, val).stop
        if bin_index == length(bins) # got the max value
            bin_index = bin_index - 1
        elseif isapprox(val, bins[1]) # sometimes will get 0 if this doesn't get checked
            bin_index = 1
        end
    end
    onehot_vec[bin_index] = 1.0
    return onehot_vec
end


text:
using OMEinsum: _unique

@testset "utils" begin
    @test _unique(Int,(1,2,3,3,)) == [1,2,3]
    @test asarray(3) isa Array
    @test asarray(3, randn(3,3)) isa Array
    @test asarray(randn(3,3)) isa Array
    @test asscalar(fill(3,())) === 3
    @test asscalar(3) === 3
end

@testset "allunique" begin
    @test allunique(())
    @test allunique(('i',))
    @test allunique((1,2))
    @test !allunique((1,2,1))
end

@testset "tensorpermute" begin
    a = randn(100, 100)
    @test OMEinsum.tensorpermute(a, [1,2]) == a
    @test OMEinsum.tensorpermute(a, (2,1)) == transpose(a)
end

@testset "align_types" begin
    a = randn(100, 100)
    b = randn(ComplexF64, 100, 100)
    res = OMEinsum.align_eltypes(a, b)
    @test res[2] === b
    @test res == (ComplexF64.(a), b)
end

@testset "batched gemm" begin
    A = randn(10, 10, 10)
    B = randn(10, 10, 10)
    for C1 in ['N', 'T']
        for C2 in ['N', 'T']
            A_ = Array{Any}(A)
            B_ = Array{Any}(B)
            @test OMEinsum._batched_gemm(C1, C2, A, B) ≈ OMEinsum._batched_gemm(C1, C2, A_, B_)
        end
    end
end


text:
@testset "SIFTS Mappings" begin

    sifts_file = joinpath(pwd(), "data", "2vqc.xml.gz")

    @testset "parse" begin

        dbs = [(dbUniProt, "P20220"), (dbPfam, "PF09645"), (dbNCBI, "244589"),
               (dbPDB, "2vqc"), (dbSCOP, "153426"), (dbPDBe, "2vqc")]
        for (to,toid) in dbs, (from, fromid) in dbs
            map = siftsmapping(sifts_file, from, fromid, to, toid)
            for (k,v) in map
                to == from && @test k == v
            end
        end
    end

    @testset "missings = false" begin

        map = siftsmapping(sifts_file,dbPDBe,"2vqc",dbPDB,"2vqc",chain="A",missings=false)
        @test_throws KeyError map["9"]  # Missing
        @test_throws KeyError map["80"] # Missing
        @test_throws KeyError map["1"]  # Missing
        @test map["10"] == "4"
        @test map["79"] == "73"
    end

    @testset "missings = true" begin

        map = siftsmapping(sifts_file, dbPDBe, "2vqc", dbPDB, "2vqc", chain="A")
        @test map["9"]  == "3"   # Missing
        @test map["80"] == "74"  # Missing
        @test map["1"]  == "-5"  # Missing # Negative Resnum
        @test map["10"] == "4"
        @test map["79"] == "73"
    end

    @testset "Insert codes" begin
    # 1SSX : Residues with insert codes: 15A 15B

        _1ssx_file = joinpath(pwd(), "data", "1ssx.xml.gz")

        map = siftsmapping(_1ssx_file, dbPDBe, "1ssx", dbPDB, "1ssx", chain="A")
        residue_A = map["1"]
        residue_B = map["2"]
        residue_C = map["3"]
        @test residue_A == "15A"
        @test residue_B == "15B"
        @test residue_C == "16"

        mapII = siftsmapping(_1ssx_file, dbPDB, "1ssx", dbUniProt, "P00778", chain="A")
        @test mapII["15A"] == "200"
        @test mapII["15B"] == "201"
        @test mapII["16"]  == "202"
    end

    @testset "Multiple InterProt annotations" begin
    # 1CBN : Multiple InterProt annotations, the last is used.
    # Identical PDBe ResNum for Residue 22:
    #
    #         <residue dbSource="PDBe" dbCoordSys="PDBe" dbResNum="22" dbResName="SER">
    #           <crossRefDb dbSource="PDB" dbCoordSys="PDBresnum" dbAccessionId="1cbn" dbResNum="22" dbResName="SER" dbChainId="A"/>
    #           <crossRefDb dbSource="UniProt" dbCoordSys="UniProt" dbAccessionId="P01542" dbResNum="22" dbResName="P"/>
    #           ....
    #         <residue dbSource="PDBe" dbCoordSys="PDBe" dbResNum="22" dbResName="PRO">
    #           <crossRefDb dbSource="PDB" dbCoordSys="PDBresnum" dbAccessionId="1cbn" dbResNum="22" dbResName="PRO" dbChainId="A"/>
    #           <crossRefDb dbSource="UniProt" dbCoordSys="UniProt" dbAccessionId="P01542" dbResNum="22" dbResName="P"/>
    #           ...

        _1cbn_file = joinpath(pwd(), "data", "1cbn.xml.gz")

        map = siftsmapping(_1cbn_file, dbPDBe, "1cbn", dbInterPro, "IPR001010", chain="A")
        @test_throws KeyError map["1"] # Without InterPro
        @test map["2"] == "2" # Same ResNum for different InterPros

        mapII = read(_1cbn_file, SIFTSXML, chain="A")
        @test length(mapII[1].InterPro) == 0 # Without InterPro
        @test length(mapII[2].InterPro) == 4 # Same ResNum for different InterPros
    end

    @testset "None" begin
    # 4CPA : It has "None"
    # <residue dbSource="PDBe" dbCoordSys="PDBe" dbResNum="2" dbResName="GLX">
    #   <crossRefDb dbSource="PDB" dbCoordSys="PDBresnum" dbAccessionId="4cpa" dbResNum="2" dbResName="GLX" dbChainId="J"/>
    #   <crossRefDb dbSource="SCOP" dbCoordSys="PDBresnum" dbAccessionId="118683" dbResNum="2" dbResName="GLX" dbChainId="J"/>
    #   <crossRefDb dbSource="InterPro" dbCoordSys="UniProt" dbAccessionId="IPR011052" dbResNum="None" dbResName="None" dbEvidence="SSF57027"/>
    #   <crossRefDb dbSource="InterPro" dbCoordSys="UniProt" dbAccessionId="IPR021142" dbResNum="None" dbResName="None" dbEvidence="PD884054"/>
    #   <crossRefDb dbSource="InterPro" dbCoordSys="UniProt" dbAccessionId="IPR004231" dbResNum="None" dbResName="None" dbEvidence="PF02977"/>
    #   <residueDetail dbSource="PDBe" property="codeSecondaryStructure">T</residueDetail>
    #   <residueDetail dbSource="PDBe" property="nameSecondaryStructure">loop</residueDetail>
    # </residue>

        _4cpa_file = joinpath(pwd(), "data", "4cpa.xml.gz")

        map = read(_4cpa_file, SIFTSXML, chain="J")
        res = filter(r -> r.PDBe.number == "2", map)[1]
        @test length(res.InterPro) == 3
        for i in 1:3
            @test res.InterPro[i].name == ""
            @test res.InterPro[i].number == ""
        end
    end

    @testset "NMR" begin
    # 1AS5 : NMR

        _1as5_file = joinpath(pwd(), "data", "1as5.xml.gz")

        map = siftsmapping(_1as5_file, dbPDBe,"1as5", dbUniProt, "P56529", chain="A")
                          # missings=true : NMR there are not missing residues
        @test map["23"] == "73"
        @test_throws KeyError map["24"] # Without UniProt

        mapII = siftsmapping(_1as5_file, dbPDBe,"1as5", dbPDB, "1as5", chain="A")
        @test mapII["24"] == "24"
    end

    @testset "Inserted residues lack insertion code" begin
    # 1DPO : Inserted residues lack insertion letters
    # Single unnamed chain in 1DPO contains insertions at postions 184 (Gly, Phe),
    # 188 (Gly, Lys), and 221 (Ala, Leu) but no insertion letters.

        _1dpo_file = joinpath(pwd(), "data", "1dpo.xml.gz")
        map = siftsmapping(_1dpo_file, dbPDBe, "1dpo", dbPDB, "1dpo", chain="A")
        # Unnamed chain is "A" in SIFTS
        @test map["164"] == "184"
        @test map["165"] == "184A" # Has insertion code in SIFTS
        @test map["169"] == "188"
        @test map["170"] == "188A" # Has insertion code in SIFTS
        @test map["198"] == "221"
        @test map["199"] == "221A" # Has insertion code in SIFTS
    end

    @testset "Insertion block" begin
    # 1IGY : Insertions have more than one copy of the same amino acid in a single
    # insertion block. For example, chain B in 1IGY contains a block of four residues
    # inserted at sequence position 82. The block contains Leu-Ser-Ser-Leu.

        _1igy_file = joinpath(pwd(), "data", "1igy.xml.gz")
        map = siftsmapping(_1igy_file, dbPDBe, "1igy", dbCATH, "2.60.40.10", chain="B")
        @test map["82"] == "82"
        @test map["83"] == "82A"
        @test map["84"] == "82B"
        @test map["85"] == "82C"
    end

    @testset "1HAG" begin
    # 1HAG : Chain E begins with 1H, 1G, 1F, ... 1A, then 1 (in reverse alphabetic order)

        _1hag_file = joinpath(pwd(), "data", "1hag.xml.gz")
        map = siftsmapping(_1hag_file, dbPDBe, "1hag", dbPDB, "1hag", chain="E")
        @test map["1"] == "1H"
        @test map["2"] == "1G"
        @test map["3"] == "1F"
        @test map["4"] == "1E"
        @test map["5"] == "1D"
        @test map["6"] == "1C"
        @test map["7"] == "1B"
        @test map["8"] == "1A"
        @test map["9"] == "1"
    end
end

@testset "1NSA" begin
# 1NSA : Contains a single (unnamed) protein chain with sequence 7A-95A that continues 4-308.

    _1nsa_file = joinpath(pwd(), "data", "1nsa.xml.gz")
    mapping = read(_1nsa_file, SIFTSXML)

    @testset "findall & read" begin

        four = findall(db -> db.id == "1nsa" && db.number == "4" , mapping, dbPDB)[1]
        @test findall(db -> db.id == "1nsa" && db.number == "95A", mapping, dbPDB)[1]+1 == four
        @test mapping[findall(db -> db.id == "1nsa" && db.number == "95A",mapping,dbPDB)][1].PDB.number == "95A"
    end

    @testset "capture fields" begin

        cap = map(mapping) do res
            number = get(res, dbPDB, :number, "")
            if get(res, dbPDB, :id, "") == "1nsa" && number == "95A"
                number
            else
                missing
            end
        end
        @test cap[ [!ismissing(x) for x in cap] ][1] == "95A"
    end
end

@testset "find & filter" begin
# 1IAO : Contains in chain B (in this order) 1S, 323P-334P, 6-94, 94A, 95-188, 1T, 2T

    mapp = read(joinpath(pwd(), "data", "1iao.xml.gz"), SIFTSXML)

    @test filter(db -> db.id == "1iao" && db.number == "1S" && db.chain == "B", mapp, dbPDB)[1].PDBe.number == "1"
    i = findall(db -> db.id == "1iao" && db.number == "1S" && db.chain == "B", mapp, dbPDB)[1]
    res = mapp[i+2]
    @test get(res,dbPDB,:id,"") == "1iao" &&  get(res,dbPDB,:chain,"") == "B" &&
          get(res,dbPDB,:number,"") == "323P"
end

@testset "MIToS 1.0 error" begin

    sf = joinpath(pwd(), "data", "4gcr.xml.gz")

    mapping = siftsmapping(sf, dbPfam, "PF00030", dbPDB, "4gcr")
    @test mapping["3"] == "2"
end

@testset "download" begin

    pdb = "2vqc"
    mapping = read(joinpath(pwd(), "data", "$(pdb).xml.gz"), SIFTSXML)

    @test_throws AssertionError downloadsifts(pdb, filename="bad_name.txt")
    @test_throws ErrorException downloadsifts("2vqc_A")
    # FTP causes problem with Travis-CI
    # filename = downloadsifts(pdb, filename=tempname()*".xml.gz")
    # @test length(read(filename, SIFTSXML)) == length(mapping)
end

@testset "Ensembl" begin
# 18GS has multiple Ensembl annotations for each residue
# 18GS also has EC and GO annotations

    mapping = read(joinpath(pwd(), "data", "18gs.xml.gz"), SIFTSXML)

    last_res = mapping[end]
    @test length(last_res.Ensembl) == 2
    @test last_res.Ensembl[1].id == last_res.Ensembl[2].id
    @test last_res.Ensembl[2].transcript == "ENST00000642444"
    @test last_res.Ensembl[1].transcript == "ENST00000398606"
    @test last_res.Ensembl[2].transcript == "ENST00000642444"
    @test last_res.Ensembl[1].translation == "ENSP00000381607"
    @test last_res.Ensembl[2].translation == "ENSP00000493538"
    @test last_res.Ensembl[1].exon == "ENSE00001921020"
    @test last_res.Ensembl[2].exon == "ENSE00003822108"
end


text:
leftgroupjoin(left, right) = leftgroupjoin(identity, identity, left, right)
leftgroupjoin(lkey::Callable, rkey::Callable, left, right) = leftgroupjoin(lkey, rkey, merge, left, right)
leftgroupjoin(lkey::Callable, rkey::Callable, f::Callable, left, right) = leftgroupjoin(lkey, rkey, f, isequal, left, right)

"""
    leftgroupjoin(lkey, rkey, f, comparison, left, right)

Creates a collection if groups labelled by `lkey(l)` where each group contains elements
`f(l, r)` which satisfy `comparison(lkey(l), rkey(r))`. If there are no matches, the group
is still created (with an empty collection).

This operation shares some similarities with an SQL left outer join.

### Example

```jldoctest
julia> leftgroupjoin(iseven, iseven, tuple, ==, [1,2,3,4], [0,1,2])
Dictionary{Bool,Array{Tuple{Int64,Int64},1}} with 2 entries:
  false │ Tuple{Int64,Int64}[(1, 1), (3, 1)]
  true  │ Tuple{Int64,Int64}[(2, 0), (2, 2), (4, 0), (4, 2)]
```
"""
function leftgroupjoin(lkey::Callable, rkey::Callable, f::Callable, comparison::Callable, left, right)
    # The O(length(left)*length(right)) generic method when nothing about `comparison` is known

    # TODO Do this inference-free, like comprehensions...
    T = Core.Compiler.return_type(f, Tuple{eltype(left), eltype(right)})
    K = Core.Compiler.return_type(lkey, Tuple{eltype(left)})
    V = Vector{T}
    out = Dictionary{K, V}()
    for a ∈ left
        key = lkey(a)
        group = get!(V, out, key)
        for b ∈ right
            if comparison(key, rkey(b))
                push!(group, f(a, b))
            end
        end
    end
    return out
end

function leftgroupjoin(lkey::Callable, rkey::Callable, f::Callable, ::typeof(isequal), left, right)
    # isequal heralds a hash-based approach, roughly O(length(left) * log(length(right)))

    # TODO Do this inference-free, like comprehensions...
    T = Core.Compiler.return_type(f, Tuple{eltype(left), eltype(right)})
    K = Core.Compiler.return_type(rkey, Tuple{eltype(right)})
    V = eltype(right)
    dict = Dictionary{K,Vector{V}}() # maybe a different stategy if right is unique
    for b ∈ right
        key = rkey(b)
        push!(get!(Vector{V}, dict, key), b)
    end

    K2 = Core.Compiler.return_type(lkey, Tuple{eltype(left)})
    out = Dictionary{K2, Vector{T}}()
    for a ∈ left
        key = lkey(a)
        token = gettoken(out, key)
        group = get!(Vector{T}, out, key)
        (has_index, token) = gettoken(dict, key)
        if has_index
            for b ∈ gettokenvalue(dict, token)
                push!(group, f(a, b))
            end
        end
    end
    return out
end


text:
using Test
using DFTK
using LinearMaps
import DFTK: ortho_qr, pack_ψ, unpack_ψ, reinterpret_real, reinterpret_complex
import DFTK: proj_tangent, proj_tangent!
import DFTK: apply_K, apply_Ω, filled_occupation, select_occupied_orbitals
import DFTK: precondprep!, FunctionPreconditioner

include("testcases.jl")

if mpi_nprocs() == 1  # Distributed implementation not yet available

    function eigen_ΩplusK(basis::PlaneWaveBasis{T}, ψ, numval) where T
        # check that there are no virtual orbitals
        model = basis.model
        filled_occ = filled_occupation(model)
        n_spin = model.n_spin_components
        n_bands = div(model.n_electrons, n_spin * filled_occ, RoundUp)
        @assert n_bands == size(ψ[1], 2)

        occupation = [filled_occ * ones(T, n_bands) for kpt = basis.kpoints]

        pack(ψ) = reinterpret_real(pack_ψ(ψ))
        unpack(x) = unpack_ψ(reinterpret_complex(x), size.(ψ))

        # compute quantites at the point which define the tangent space
        ρ = compute_density(basis, ψ, occupation)
        _, H = energy_hamiltonian(basis, ψ, occupation; ρ=ρ)

        # preconditioner
        Pks = [PreconditionerTPA(basis, kpt) for kpt in basis.kpoints]
        for (ik, ψk) in enumerate(ψ)
            precondprep!(Pks[ik], ψk)
        end
        function f_ldiv!(x, y)
            for n in 1:size(y, 2)
                δψ = unpack(y[:, n])
                proj_tangent!(δψ, ψ)
                Pδψ = [ Pks[ik] \ δψk for (ik, δψk) in enumerate(δψ)]
                proj_tangent!(Pδψ, ψ)
                x[:, n] .= pack(Pδψ)
            end
            x
        end

        # random starting point on the tangent space to avoid eigenvalue 0
        x0 = map(1:numval) do n
            xn = proj_tangent([randn(Complex{eltype(basis)}, length(G_vectors(kpt)), n_bands)
                               for kpt in basis.kpoints], ψ)
            pack(xn)
        end
        x0 = hcat(x0...)

        # Rayleigh-coefficients
        Λ = [ψk'Hψk for (ψk, Hψk) in zip(ψ, H * ψ)]

        # mapping of the linear system on the tangent space
        function ΩpK(x)
            δψ = unpack(x)
            Kδψ = apply_K(basis, δψ, ψ, ρ, occupation)
            Ωδψ = apply_Ω(δψ, ψ, H, Λ)
            pack(Ωδψ + Kδψ)
        end
        J = LinearMap{T}(ΩpK, size(x0, 1))

        # compute smallest eigenvalue of Ω with LOBPCG
        lobpcg_hyper(J, x0, prec=FunctionPreconditioner(f_ldiv!), tol=1e-7)
    end

    @testset "Compute eigenvalues" begin
        @testset "Compute smallest eigenvalue of Ω" begin
            Ecut = 5
            tol = 1e-12
            numval = 3  # number of eigenvalues we want to compute

            Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.psp))
            model = model_atomic(silicon.lattice, [Si => silicon.positions])
            basis = PlaneWaveBasis(model; Ecut, kgrid=[1,1,1])

            scfres = self_consistent_field(basis; tol=tol)
            ψ = select_occupied_orbitals(basis, scfres.ψ)

            res = eigen_ΩplusK(basis, ψ, numval)
            gap = scfres.eigenvalues[1][5] - scfres.eigenvalues[1][4]

            # in the linear case, the smallest eigenvalue of Ω is the gap
            @test abs(res.λ[1] - gap) < 1e-5
            @test res.λ[1] > 1e-3
        end

        @testset "Compute smallest eigenvalue of Ω+K" begin
            Ecut = 5
            tol = 1e-12
            numval = 3  # number of eigenvalues we want to compute

            Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.psp))
            model = model_LDA(silicon.lattice, [Si => silicon.positions])
            basis = PlaneWaveBasis(model; Ecut, kgrid=[1,1,1])

            scfres = self_consistent_field(basis; tol=tol)
            ψ = select_occupied_orbitals(basis, scfres.ψ)

            res = eigen_ΩplusK(basis, ψ, numval)

            @test res.λ[1] > 1e-3
        end
    end
end


text:
using AlphaVantage, Test

@testset "http" begin
  @testset "client instantiation" begin
      c = AlphaVantageClient()
      @test c.key == "demo"

      c = AlphaVantageClient(;key="123")
      @test c.key == "123"

      c = AlphaVantageClient("xyz")
      @test c.key == "xyz"
  end

  @testset "direct get" begin
      c = AlphaVantageClient()
      p = Dict(
          "function" => "TIME_SERIES_MONTHLY_ADJUSTED",
          "symbol" => "IBM",
          "apikey" => "demo"
      )
      resp = get(c, p)
      @test typeof(resp) == Vector{UInt8}
      @test length(resp) > 0
  end
end

@testset "handlers" begin
  c = AlphaVantageClient()

  @testset "time series" begin
    ts = get(c, TIME_SERIES_MONTHLY(), Dict("symbol" => "AAPL"))
    @test typeof(ts) == DataFrame
    @test nrow(ts) > 0
  end

  @testset "fundamentals" begin

  end

  @testset "currency" begin

  end

  @testset "economic" begin

  end

  @testset "technical" begin

  end
end


text:
devices = Dict{Symbol, DeviceModel}(:Generators => DeviceModel(PSY.ThermalStandard, ThermalDispatch),
                                    :Loads =>  DeviceModel(PSY.PowerLoad, StaticPowerLoad))
branches = Dict{Symbol, DeviceModel}(:L => DeviceModel(PSY.Line, StaticLineUnbounded))
services = Dict{Symbol, ServiceModel}()

@testset "Operation set ref models" begin
    template = OperationsTemplate(CopperPlatePowerModel, devices, branches, services);
    op_problem = OperationsProblem(TestOpProblem, template, c_sys5)
    set_transmission_model!(op_problem, DCPLLPowerModel)
    @test op_problem.template.transmission == DCPLLPowerModel

    template = OperationsTemplate(CopperPlatePowerModel, devices, branches, services);
    op_problem = OperationsProblem(TestOpProblem, template, c_sys5)
    new_devices = Dict{Symbol, DeviceModel}(:Generators => DeviceModel(PSY.ThermalStandard, ThermalBasicUnitCommitment),
                                            :Loads =>  DeviceModel(PSY.PowerLoad, StaticPowerLoad))
    set_devices_template!(op_problem, new_devices)
    @test op_problem.template.devices[:Generators].formulation == ThermalBasicUnitCommitment
    jump_model = op_problem.canonical.JuMPmodel
    @test ((JuMP.VariableRef, MOI.ZeroOne) in JuMP.list_of_constraint_types(jump_model)) == true

    template = OperationsTemplate(DCPPowerModel, devices, branches, services);
    op_problem = OperationsProblem(TestOpProblem, template, c_sys5)
    new_branches = Dict{Symbol, DeviceModel}(:L => DeviceModel(PSY.Line, StaticLine))
    set_branches_template!(op_problem, new_branches)
    @test op_problem.template.branches[:L].formulation == StaticLine
end

@testset "Operation set models" begin
    template = OperationsTemplate(CopperPlatePowerModel, devices, branches, services);
    op_problem = OperationsProblem(TestOpProblem, template, c_sys5)
    set_device_model!(op_problem, :Generators, DeviceModel(PSY.ThermalStandard, ThermalBasicUnitCommitment))
    @test op_problem.template.devices[:Generators].formulation == ThermalBasicUnitCommitment
    jump_model = op_problem.canonical.JuMPmodel
    @test ((JuMP.VariableRef, MOI.ZeroOne) in JuMP.list_of_constraint_types(jump_model)) == true

    template = OperationsTemplate(DCPPowerModel, devices, branches, services);
    op_problem = OperationsProblem(TestOpProblem, template, c_sys5)
    set_branch_model!(op_problem, :L, DeviceModel(PSY.Line, StaticLine))
    @test op_problem.template.branches[:L].formulation == StaticLine
end