text:
using SparseArrays
using KernelAbstractions
using KernelAbstractions.Extras
using CUDA
using Random
using Adapt
using LinearAlgebra
using Printf

function D2_split_naive(idata,odata,Nx,Ny,h,::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
    tidx = threadIdx().x
    tidy = threadIdx().y

    i = (blockIdx().x - 1) * TILE_DIM1 + tidx
    j = (blockIdx().y - 1) * TILE_DIM2 + tidy

    if 1 <= i <= Nx && 1 <= j <= Ny
        odata[i,j] = 0 # maybe do this on local memory
        # odata[i,j] = 0 # maybe do this on local memory
    end
    sync_threads()

    if 2 <= i <= Nx-1 && 2 <= j <= Ny - 1
        @inbounds   odata[i,j] = (idata[i-1,j] + idata[i+1,j] + idata[i,j-1] + idata[i,j+1] - 4*idata[i,j])
        # odata[i,j] = (idata[i-1,j] + idata[i+1,j] + idata[i,j-1] + idata[i,j+1] - 4*idata[i,j])
    end

    nothing
end

function D2_split_naive_v2(idata,odata,Nx,Ny,h,::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
    tidx = threadIdx().x
    tidy = threadIdx().y

    i = (blockIdx().x - 1) * TILE_DIM1 + tidx
    j = (blockIdx().y - 1) * TILE_DIM2 + tidy

    if i <= Nx && j <= Ny
        if i == 1 || i == Nx || j == 1 || j == Ny
            @inbounds odata[i,j] = 0
        else
            @inbounds  odata[i,j] = (idata[i-1,j] + idata[i+1,j] + idata[i,j-1] + idata[i,j+1] - 4*idata[i,j])
       end
    #    @inbounds odata[i,j] = idata[i,j]
    end

    nothing
end

function copy_naive!(b, a)
    N = size(a, 1)

    # which thread are we in our block
    tidx = threadIdx().x
    tidy = threadIdx().y

    # which block of threads are we in
    bidx = blockIdx().x
    bidy = blockIdx().y

    # What is the size of the thread block
    dimx = blockDim().x
    dimy = blockDim().y

    # what index am I in the global thread space
    i = tidx + dimx * (bidx - 1)
    j = tidy + dimy * (bidy - 1)

    if i <= N && j <= N
        # aval = a[i, j]
        # b[i, j] = aval
        @inbounds b[i, j] = a[i, j]
    end

    return nothing
end

function D2_split(idata,odata,Nx,Ny,h,::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
    tidx = threadIdx().x
    tidy = threadIdx().y

    i = (blockIdx().x - 1) * TILE_DIM1 + tidx
    j = (blockIdx().y - 1) * TILE_DIM2 + tidy

    si = tidx + 1
    sj = tidy + 1

    HALO_WIDTH = 2 # For second order derivative

    tile = @cuStaticSharedMem(eltype(idata), (TILE_DIM1+2, TILE_DIM2+2))

    global_index = (i-1)*Ny+j
    if 0 <= i <= Nx && 1 <= j <= Ny
        @inbounds  odata[i,j] = 0 # maybe do this on local memory
    end

    if 1 <= i <= Nx && 1 <= j <= Ny
        @inbounds  tile[si,sj] = idata[i,j]
    end
    sync_threads()

    # if tidx == 1 && 1 <= j <= Ny
    #     tile[si-1,sj] = tile[Nx+si-2,sj]
    # #     tile[si+Nx,sj] = tile[si,sj]
    # end

    # sync_threads()

    # if 1 <= i <= Nx && j <= 1
    #     tile[si,sj-1] = tile[si,Ny+sj-2]
    #     tile[si,sj+Ny] = tile[si,sj]
    # end

    sync_threads()

    if 2 <= i <= Nx -1 && 2 <= j <= Ny - 1 &&  1 <= si <= TILE_DIM1 + 1 && 1 <= sj <= TILE_DIM2 + 1
        @inbounds odata[i,j] = (tile[si-1,sj] + tile[si+1,sj] + tile[si,sj-1] + tile[si,sj+1] - 4*tile[si,sj])
    end

    nothing
end

struct boundary_containers
    CPU_W
    CPU_W_T
    CPU_E
    CPU_E_T
    CPU_N
    CPU_S
    CPU_OUT_W
    CPU_OUT_E
    CPU_OUT_N
    CPU_OUT_S
end

function matrix_free_A(idata,odata)
    # odata .= 0
    Nx,Ny = size(idata)
    h = 1/(Nx-1)
    TILE_DIM_1 = 16
    TILE_DIM_2 = 16
    griddim = (div(Nx,TILE_DIM_1) + 1, div(Ny,TILE_DIM_2) + 1)
	blockdim = (TILE_DIM_1,TILE_DIM_2)

    # CPU_W = zeros(Nx,3)
    # CPU_W_T = zeros(3,Nx)
    # CPU_E = zeros(Nx,3)
    # CPU_E_T = zeros(3,Nx)
    # CPU_N = zeros(3,Ny)
    # CPU_S = zeros(3,Ny)

    CPU_W = Array{Float64,2}(undef,Nx,3)
    CPU_W_T = Array{Float64,2}(undef,3,Nx)
    CPU_E = Array{Float64,2}(undef,Nx,3)
    CPU_E_T = Array{Float64,2}(undef,3,Nx)
    CPU_N = Array{Float64,2}(undef,3,Ny)
    CPU_S = Array{Float64,2}(undef,3,Ny)

    # CPU_OUT_W = Array{Float64,2}(undef,Nx,3)
    # CPU_OUT_W_T = Array{Float64,2}(undef,3,Nx)
    # CPU_OUT_E = Array{Float64,2}(undef,Nx,3)
    # CPU_OUT_E_T = Array{Float64,2}(undef,3,Nx)
    # CPU_OUT_N = Array{Float64,2}(undef,3,Ny)
    # CPU_OUT_S = Array{Float64,2}(undef,3,Ny)

    copyto!(CPU_W,view(idata,:,1:3))
    copyto!(CPU_E,view(idata,:,Ny-2:Ny))
    copyto!(CPU_N,view(idata,1:3,:))
    copyto!(CPU_S,view(idata,Nx-2:Nx,:))

    @cuda threads=blockdim blocks=griddim D2_split_naive(idata,odata,Nx,Ny,h,Val(TILE_DIM_1), Val(TILE_DIM_2))

    CPU_OUT_W = zeros(Nx,3)
    CPU_OUT_W_T = zeros(3,Nx)
    CPU_OUT_E = zeros(Nx,3)
    CPU_OUT_E_T = zeros(3,Nx)
    CPU_OUT_N = zeros(3,Ny)
    CPU_OUT_S = zeros(3,Ny)

    # CPU_OUT_W = Array{Float64,2}(undef,Nx,3)
    # CPU_OUT_W_T = Array{Float64,2}(undef,3,Nx)
    # CPU_OUT_E = Array{Float64,2}(undef,Nx,3)
    # CPU_OUT_E_T = Array{Float64,2}(undef,3,Nx)
    # CPU_OUT_N = Array{Float64,2}(undef,3,Ny)
    # CPU_OUT_S = Array{Float64,2}(undef,3,Ny)

    ## CPU Calculation

    CPU_W_T .= CPU_W'
    CPU_E_T .= CPU_E'
    # CPU_OUT_W_T = Array{Float64,2}(undef,3,Ny)
    # CPU_OUT_E_T = Array{Float64,2}(undef,3,Ny)

    alpha1 = alpha2 = -13/h
    alpha3 = alpha4 = -1
    beta = 1

    i = 1
    # Threads.@threads for j in 2:Ny-1
    # @inbounds for j in 2:Ny-1
    for j in 2:Ny-1
        @inbounds CPU_OUT_N[1,j] += (CPU_N[1,j] - 2*CPU_N[2,j] + CPU_N[3,j] + CPU_N[1,j-1] - 2* CPU_N[1,j] + CPU_N[1,j+1] + 2 * alpha3 * (1.5 * CPU_N[1,j] - 2*CPU_N[2,j] + 0.5*CPU_N[3,j])) / 2
    end
    synchronize()

    i = Nx
    # # Threads.@threads for j in 2:Ny-1
    # @inbounds for j in 2:Ny-1
    for j in 2:Ny-1
        @inbounds CPU_OUT_S[3,j] += (CPU_S[3,j] - 2*CPU_S[2,j] + CPU_S[1,j] + CPU_S[3,j-1] - 2* CPU_S[3,j] + CPU_S[3,j+1] + 2 * alpha4 * (1.5 * CPU_S[3,j] - 2*CPU_S[2,j] + 0.5*CPU_S[1,j])) / 2
    end
    # synchronize()

    j = 1

    # @inbounds for i in 2:Nx-1
    for i in 2:Nx-1
        @inbounds CPU_OUT_W_T[1,i] += (CPU_W_T[1,i-1] - 2*CPU_W_T[1,i] + CPU_W_T[1,i+1] + CPU_W_T[1,i] - 2*CPU_W_T[2,i] + CPU_W_T[3,i]) / 2
        @inbounds CPU_OUT_W_T[1,i] += (2 * beta * (1.5 * CPU_W_T[1,i]) + 2 * alpha2 * CPU_W_T[1,i] * h) / 2 # should be using alpha1 maybe
        @inbounds CPU_OUT_W_T[2,i] += (2 * beta * (-1 * CPU_W_T[1,i]))
        @inbounds CPU_OUT_W_T[3,i] += (0.5 * beta * CPU_W_T[1,i])
    end

    # # odata[:,1] .+= odata_W_T[1,:]

    j = Ny
    # @inbounds for i in 2:Nx-1
    for i in 2:Nx-1
        @inbounds CPU_OUT_E_T[3,i] += (CPU_E_T[3,i-1] - 2*CPU_E_T[3,i] + CPU_E_T[3,i+1] + CPU_E_T[3,i] - 2*CPU_E_T[2,i] + CPU_E_T[1,i]) / 2
        @inbounds CPU_OUT_E_T[3,i] += (2 * beta * (1.5 * CPU_E_T[3,i]) + 2 * alpha1 * CPU_E_T[3,i] * h) / 2
        @inbounds CPU_OUT_E_T[2,i] += (2 * beta * (-1 * CPU_E_T[3,i]))
        @inbounds CPU_OUT_E_T[1,i] += (0.5 * beta * CPU_E_T[3,i])
    end

    (i,j) = (1,1)

    @inbounds CPU_OUT_N[1,j] += (CPU_W[i,j] - 2*CPU_W[i+1,j] + CPU_W[i+2,j] + CPU_W[i,j] - 2*CPU_W[i,j+1] + CPU_W[i,j+2]) / 4 # D2

    @inbounds CPU_OUT_N[1,j] += 2 * alpha3 * (( 1.5* CPU_W[i,j] - 2*CPU_W[i+1,j] + 0.5*CPU_W[i+2,j])) / 4 # Neumann

    @inbounds CPU_OUT_N[1,j] += (2 * beta * (1.5 * CPU_W[i,j]) + 2 * alpha1 * (CPU_W[i,j]) * h) / 4 # Dirichlet
    @inbounds CPU_OUT_N[1,j+1] += (2 * beta * (-1 * CPU_W[i,j])) / 2 # Dirichlet
    @inbounds CPU_OUT_N[1,j+2] += (0.5 * beta * (CPU_W[i,j])) / 2 # Dirichlet

    (i,j) = (1,Ny)
    @inbounds CPU_OUT_N[1,j] += (CPU_E[i,3] - 2*CPU_E[i+1,3] + CPU_E[i+2,3] + CPU_E[i,3] - 2*CPU_E[i,2] + CPU_E[i,1]) / 4 # D2

    @inbounds CPU_OUT_N[1,j] += 2 * alpha3 * (1.5 * CPU_E[i,3] - 2*CPU_E[i+1,3] + 0.5 * CPU_E[i+2,3]) / 4 # Neumann
    @inbounds CPU_OUT_N[1,j] += (2 * beta * (1.5 * CPU_E[i,3]) + 2 * alpha2 * (CPU_E[i,3]) * h) / 4 # Dirichlet
    @inbounds CPU_OUT_N[1,j-1] += (2 * beta * (-1 * CPU_E[i,3])) / 2 # Dirichlet
    @inbounds CPU_OUT_N[1,j-2] += (0.5 * beta * (CPU_E[i,3])) / 2 # Dirichlet

    (i,j) = (Nx,1)
    @inbounds CPU_OUT_S[3,j] += (CPU_W[i,j] - 2*CPU_W[i-1,j] + CPU_W[i-2,j] + CPU_W[i,j] - 2*CPU_W[i,j+1] + CPU_W[i,j+2]) / 4 # D2

    @inbounds CPU_OUT_S[3,j] += 2 * alpha4 * (( 1.5* CPU_W[i,j] - 2*CPU_W[i-1,j] + 0.5*CPU_W[i-2,j])) / 4 # Neumann
    @inbounds CPU_OUT_S[3,j] += (2 * beta * (1.5 * CPU_W[i,j]) + 2 * alpha1 * (CPU_W[i,j]) * h) / 4 # Dirichlet
    @inbounds CPU_OUT_S[3,j+1] += (2 * beta * (-1 * CPU_W[i,j])) / 2 # Dirichlet
    @inbounds CPU_OUT_S[3,j+2] += (0.5 * beta * (CPU_W[i,j])) / 2 # Dirichlet

    (i,j) = (Nx,Ny)
    @inbounds CPU_OUT_S[3,j] += (CPU_E[Nx,3] - 2*CPU_E[Nx-1,3] + CPU_E[Nx-2,3] + CPU_E[Nx,3] - 2*CPU_E[Nx,2] + CPU_E[Nx,1]) / 4 # D2

    @inbounds CPU_OUT_S[3,j] += 2 * alpha4 * (1.5 * CPU_E[Nx,3] - 2*CPU_E[Nx-1,3] + 0.5 * CPU_E[Nx-2,3]) / 4 # Neumann
    @inbounds CPU_OUT_S[3,j] += (2 * beta * (1.5 * CPU_E[Nx,3]) + 2 * alpha2 * (CPU_E[Nx,3]) * h) / 4 # Dirichlet
    @inbounds CPU_OUT_S[3,j-1] += (2 * beta * (-1 * CPU_E[Nx,3])) / 2 # Dirichlet
    @inbounds CPU_OUT_S[3,j-2] += (0.5 * beta * (CPU_E[Nx,3])) / 2 # Dirichlet

    ## End CPU calculation

    @inbounds CPU_OUT_E .= CPU_OUT_E_T'
    @inbounds CPU_OUT_W .= CPU_OUT_W_T'

    synchronize()

    # @show CPU_OUT_N
    # @show CPU_OUT_S
    # @show CPU_OUT_W
    # @show CPU_OUT_E

    # show(stdout, "text/plain", odata)
    # println()

    # show(stdout, "text/plain", CPU_OUT_W)
    # println()
    # show(stdout, "text/plain", CPU_OUT_E)
    # println()
    # show(stdout, "text/plain", CPU_OUT_N)
    # println()
    # show(stdout, "text/plain", CPU_OUT_S)
    # println()

    ## Copy W & E boundary
    copyto!(view(odata,1:Nx,1:3),view(odata,1:Nx,1:3) + CuArray(CPU_OUT_W))
    copyto!(view(odata,1:Nx,Ny-2:Ny),view(odata,1:Nx,Ny-2:Ny) + CuArray(CPU_OUT_E))

    # Copy N & S boundary
    copyto!(view(odata,1,1:Ny),view(odata,1,1:Ny) + CuArray(CPU_OUT_N[1,:]))
    copyto!(view(odata,Nx,1:Ny),view(odata,Nx,1:Ny) + CuArray(CPU_OUT_S[end,:]))

    # show(stdout, "text/plain", odata)
    # println()
    synchronize()
    nothing
end

function SBP_N!(idata,odata,alpha1,alpha2,alpha3,beta,h)
    (Nx,Ny) = size(idata)
    tidx = threadIdx().x
    idx = (blockIdx().x - 1) * blockDim().x + tidx
    if 4 <= idx <= Ny - 3
        odata[1,idx] = (idata[1,idx] - 2*idata[2,idx] + idata[3,idx] + idata[1,idx-1] - 2*idata[1,idx] + idata[1,idx+1] + 2 * alpha3 * (1.5 * idata[1,idx] - 2*idata[2,idx] + 0.5*idata[3,idx])) / 2
    end
    # sync_threads()
    if idx == 1
        odata[1,idx] = (idata[1,idx] - 2*idata[1,idx+1] + idata[1,idx+2] + idata[1,idx] - 2*idata[2,idx] + idata[3,idx]
                        + 2 * alpha3 * (( 1.5* idata[1,idx] - 2*idata[2,idx] + 0.5*idata[3,idx])) ) / 4
                        # + 2 * beta * (1.5 * idata[1,idx]) + 2 * alpha1 * (idata[1,idx]) * h

        odata[1,idx+1] =  (idata[1,idx+1] - 2*idata[2,idx+1] + idata[3,idx+1] + idata[1,idx] - 2*idata[1,idx+1] + idata[1,idx+2] + 2 * alpha3 * (1.5 * idata[1,idx+1] - 2*idata[2,idx+1] + 0.5*idata[3,idx+1])) / 2 # Dirichlet
                        # (2 * beta * (-1 * idata[1,idx])) / 2 +
        odata[1,idx+2] = (idata[1,idx+2] - 2*idata[2,idx+2] + idata[3,idx+2] + idata[1,idx+1] - 2*idata[1,idx+2] + idata[1,idx+3] + 2 * alpha3 * (1.5 * idata[1,idx+2] - 2*idata[2,idx+2] + 0.5*idata[3,idx+2])) / 2# Dirichlet
                        # (0.5 * beta * (idata[1,idx])) / 2 +
    end
    # sync_threads()
    if idx == Ny
        odata[1,idx] = (idata[1,idx] - 2*idata[1,idx-1] + idata[1,idx-2] + idata[1,idx] - 2*idata[2,idx] + idata[3,idx]
                        + 2 * alpha3 * (( 1.5* idata[1,idx] - 2*idata[2,idx] + 0.5*idata[3,idx]))  ) / 4
                        # + 2 * beta * (1.5 * idata[1,idx]) + 2 * alpha2 * (idata[1,idx]) * h

        odata[1,idx-1] =  (idata[1,idx-1] - 2*idata[2,idx-1] + idata[3,idx-1] + idata[1,idx-2] - 2*idata[1,idx-1] + idata[1,idx] + 2 * alpha3 * (1.5 * idata[1,idx-1] - 2*idata[2,idx-1] + 0.5*idata[3,idx-1])) / 2 # Dirichlet
                        # (2 * beta * (-1 * idata[1,idx])) / 2 +
        odata[1,idx-2] =  (idata[1,idx-2] - 2*idata[2,idx-2] + idata[3,idx-2] + idata[1,idx-3] - 2*idata[1,idx-2] + idata[1,idx-1] + 2 * alpha3 * (1.5 * idata[1,idx-2] - 2*idata[2,idx-2] + 0.5*idata[3,idx-2])) / 2# Dirichlet
                        # (0.5 * beta * (idata[1,idx])) / 2 +
    end
    # sync_threads()
    nothing
end

function SBP_S!(idata,odata,alpha1,alpha2,alpha4,beta,h)
    (Nx,Ny) = size(idata)
    tidx = threadIdx().x
    idx = (blockIdx().x - 1) * blockDim().x + tidx
    if 4 <= idx <= Ny - 3
        odata[end,idx] = (idata[3,idx] - 2*idata[2,idx] + idata[1,idx] + idata[3,idx-1] - 2*idata[3,idx] + idata[3,idx+1] + 2 * alpha4 * (1.5 * idata[3,idx] - 2*idata[2,idx] + 0.5*idata[1,idx])) / 2
    end
    # sync_threads()
    if idx == 1
        odata[end,idx] = (idata[3,idx] - 2*idata[3,idx+1] + idata[3,idx+2] + idata[3,idx] - 2*idata[2,idx] + idata[1,idx]
                        + 2 * alpha4 * (( 1.5* idata[3,idx] - 2*idata[2,idx] + 0.5*idata[1,idx]))  ) / 4
                        # + 2 * beta * (1.5 * idata[3,idx]) + 2 * alpha1 * (idata[3,idx]) * h

        odata[end,idx+1] = (idata[3,idx+1] - 2*idata[2,idx+1] + idata[1,idx+1] + idata[3,idx] - 2*idata[3,idx+1] + idata[3,idx+2] + 2 * alpha4 * (1.5 * idata[3,idx+1] - 2*idata[2,idx+1] + 0.5*idata[1,idx+1])) / 2  # Dirichlet
                        #   (2 * beta * (-1 * idata[3,idx])) / 2
        odata[end,idx+2] = (idata[3,idx+2] - 2*idata[2,idx+2] + idata[1,idx+2] + idata[3,idx+1] - 2*idata[3,idx+2] + idata[3,idx+3] + 2 * alpha4 * (1.5 * idata[3,idx+2] - 2*idata[2,idx+2] + 0.5*idata[1,idx+2])) / 2# Dirichlet
                        #  (0.5 * beta * (idata[3,idx])) / 2 +
    end
    # sync_threads()
    if idx == Ny
        odata[end,idx] = (idata[3,idx] - 2*idata[3,idx-1] + idata[3,idx-2] + idata[3,idx] - 2*idata[2,idx] + idata[1,idx]
                        + 2 * alpha4 * (( 1.5* idata[3,idx] - 2*idata[2,idx] + 0.5*idata[1,idx]))  ) / 4
                        # + 2 * beta * (1.5 * idata[3,idx]) + 2 * alpha2 * (idata[3,idx]) * h

        odata[end,idx-1] = (idata[3,idx-1] - 2*idata[2,idx-1] + idata[1,idx-1] + idata[3,idx-2] - 2*idata[3,idx-1] + idata[3,idx] + 2 * alpha4 * (1.5 * idata[3,idx-1] - 2*idata[2,idx-1] + 0.5*idata[1,idx-1])) / 2 # Dirichlet
                        # (2 * beta * (-1 * idata[3,idx])) / 2 +
        odata[end,idx-2] =  (idata[3,idx-2] - 2*idata[2,idx-2] + idata[1,idx-2] + idata[3,idx-3] - 2*idata[3,idx-2] + idata[3,idx-1] + 2 * alpha4 * (1.5 * idata[3,idx-2] - 2*idata[2,idx-2] + 0.5*idata[1,idx-2])) / 2# Dirichlet
                        # (0.5 * beta * (idata[3,idx])) / 2 +
    end
    # sync_threads()
    nothing
end

function SBP_W!(idata,odata,alpha1,beta,h)
    (Nx,Ny) = size(idata)
    tidx = threadIdx().x
    idx = (blockIdx().x - 1) * blockDim().x + tidx
    if 2 <= idx <= Nx-1
        odata[idx,1] = (idata[idx-1,1] - 2*idata[idx,1] + idata[idx+1,1] + idata[idx,1] - 2*idata[idx,2] + idata[idx,3] + 2 * beta * (1.5 * idata[idx,1]) + 2 * alpha1 * idata[idx,1] * h) / 2
        # (2 * beta * (1.5 * CPU_W_T[1,i]) + 2 * alpha2 * CPU_W_T[1,i] * h) / 2
        odata[idx,2] = (2 * beta * (-1 * idata[idx,1]))
        odata[idx,3] = (0.5 * beta * idata[idx,1])
    end
    sync_threads()
    if idx == 1 || idx == Nx
        odata[idx,1] =  (2 * beta * (1.5 * idata[idx,1]) + 2 * alpha1 * (idata[idx,1]) * h )/ 4
        odata[idx,2] = (2 * beta * (-1 * idata[idx,1])) / 2
        odata[idx,3] = (0.5 * beta * idata[idx,1]) / 2
    end
    nothing
end

function SBP_E!(idata,odata,alpha2,beta,h)
    (Nx,Ny) = size(idata)
    tidx = threadIdx().x
    idx = (blockIdx().x - 1) * blockDim().x + tidx
    if 2 <= idx <= Nx-1
        odata[idx,3] = (idata[idx-1,3] - 2*idata[idx,3] + idata[idx+1,3] + idata[idx,3] - 2*idata[idx,2] + idata[idx,1] + 2 * beta * (1.5 * idata[idx,3]) + 2 * alpha2 * idata[idx,3] * h) / 2
        # (2 * beta * (1.5 * CPU_W_T[1,i]) + 2 * alpha2 * CPU_W_T[1,i] * h) / 2
        odata[idx,2] = (2 * beta * (-1 * idata[idx,3]))
        odata[idx,1] = (0.5 * beta * idata[idx,3])
    end
    sync_threads()
    if idx == 1 || idx == Nx
        odata[idx,3] = ( 2 * beta * (1.5 * idata[idx,3]) + 2 * alpha2 * idata[idx,3] * h ) / 4
        odata[idx,2] = (2 * beta * (-1 * idata[idx,3])) / 2
        odata[idx,1] = (0.5 * beta * idata[idx,3]) / 2
        # odata[idx,3] = 0
     end
    nothing
end

function show_nonzero_index(A)
    (M,N) = size(A)
    for i in 1:M
        for j in 1:N
            if !(isapprox(A[i,j],0,atol=1e-8))
                print((i,j))
            end
        end
    end
    println()
end

function matrix_free_A_full_GPU(idata,odata)
    Nx,Ny = size(idata)
    h = 1/(Nx-1)
    TILE_DIM_1 = 16
    TILE_DIM_2 = 16
    griddim_2d = (div(Nx,TILE_DIM_1) + 1, div(Ny,TILE_DIM_2) + 1)
	blockdim_2d = (TILE_DIM_1,TILE_DIM_2)

    alpha1 = alpha2 = -13/h
    alpha3 = alpha4 = -1
    beta = 1
    @cuda threads=blockdim_2d blocks=griddim_2d D2_split_naive_v2(idata,odata,Nx,Ny,h,Val(TILE_DIM_1), Val(TILE_DIM_2))
    tile_dim_1d = 256
    griddim_1d = cld(Nx,tile_dim_1d)
    GPU_OUT_W = CuArray{Float64,2}(undef,Nx,3)
    GPU_OUT_E = CuArray{Float64,2}(undef,Nx,3)
    GPU_OUT_N = CuArray{Float64,2}(undef,1,Ny)
    GPU_OUT_S = CuArray{Float64,2}(undef,1,Ny)

    @cuda threads=tile_dim_1d blocks=griddim_1d SBP_N!(view(idata,1:3,:),GPU_OUT_N,alpha1,alpha2,alpha3,beta,h)
    # println("finished N")
    @cuda threads=tile_dim_1d blocks=griddim_1d SBP_S!(view(idata,Nx-2:Nx,:),GPU_OUT_S,alpha1,alpha2,alpha4,beta,h)
    # println("finished S")
    # @cuda threads=tile_dim_1d blocks=griddim_1d SBP_W!_v2(view(idata,:,1:3),GPU_OUT_W,alpha1,beta,h)
    @cuda threads=tile_dim_1d blocks=griddim_1d SBP_W!(view(idata,:,1:3),GPU_OUT_W,alpha1,beta,h)

    # println("finished W")
    # println("Finished launching kernel for N,S,N\n")
    # @cuda threads=tile_dim_1d blocks=griddim_1d SBP_E!_v2(view(idata,:,Ny-2:Ny),GPU_OUT_E,alpha2,beta,h)
    @cuda threads=tile_dim_1d blocks=griddim_1d SBP_E!(view(idata,:,Ny-2:Ny),GPU_OUT_E,alpha2,beta,h)
    # println("finished E")

     # Copy W & E boundary
    copyto!(view(odata,1:Nx,1:3), view(odata,1:Nx,1:3) + view(GPU_OUT_W,1:Nx,1:3))
    copyto!(view(odata,1:Nx,Ny-2:Ny),view(odata,1:Nx,Ny-2:Ny) + view(GPU_OUT_E,1:Nx,1:3))

    # Copy N & S boundary
    copyto!(view(odata,1,1:Ny), view(odata,1,1:Ny) + view(GPU_OUT_N,1,:))
    copyto!(view(odata,Nx,1:Ny), view(odata,Nx,1:Ny) + view(GPU_OUT_S,1,:))

    nothing
end

function test_matrix_free_A(level;TILE_DIM_1=16,TILE_DIM_2=16)
    Nx = Ny = 2^level + 1
    memsize = sizeof(Float64) * Nx * Ny / 1024^3
    h = 1/(Nx-1)
    println("")
    println("Starting Test")
    println("2D Domain Size: $Nx by $Ny")
    @printf("memsize %f GB\n",memsize)
    Random.seed!(0)
    idata = CuArray(randn(Nx,Ny))
    odata = similar(idata)
    odata_v2 = similar(idata)
    odata_full_GPU = similar(idata)

    TILE_DIM_1 = TILE_DIM_2 = 16

    griddim = (cld(Nx,TILE_DIM_1),cld(Ny,TILE_DIM_2))
    blockdim = (TILE_DIM_1,TILE_DIM_2)
    @cuda threads=blockdim blocks=griddim D2_split_naive_v2(idata,odata_v2,Nx,Ny,h,Val(TILE_DIM_1), Val(TILE_DIM_2))
    @cuda threads=blockdim blocks=griddim D2_split_naive(idata,odata,Nx,Ny,h,Val(TILE_DIM_1), Val(TILE_DIM_2))
    @show norm(odata-odata_v2)
    @cuda threads=blockdim blocks=griddim copy_naive!(odata,idata)
    # @cuda threads=blockdim blocks=griddim copy_naive!(d_b,d_a)
    synchronize()
    # boundary_data = boundary_containers(zeros(Nx,3),zeros(3,Nx),zeros(Nx,3),zeros(3,Nx),zeros(3,Ny),zeros(3,Ny),zeros(Nx,3),zeros(Nx,3),zeros(3,Ny),zeros(3,Ny))
    matrix_free_A(idata,odata)
    matrix_free_A_full_GPU(idata,odata_full_GPU)
    # @show odata - odata_full_GPU
    @show norm(odata-odata_full_GPU)

    iter_times = 1000
    @show iter_times

    # Evaluating only D2_naive
    t_D2_naive = @elapsed begin
        for _ in 1:iter_times
            @cuda threads=blockdim blocks=griddim D2_split_naive(idata,odata,Nx,Ny,h,Val(TILE_DIM_1), Val(TILE_DIM_2))
        end
        synchronize()
    end

    t_D2_naive = t_D2_naive * 1000 / iter_times
    @show t_D2_naive
    # End evaluating D2_naive

     # Evaluating only D2_naive_v2
     t_D2_naive_v2 = @elapsed begin
        for _ in 1:iter_times
            @cuda threads=blockdim blocks=griddim D2_split_naive_v2(idata,odata,Nx,Ny,h,Val(TILE_DIM_1), Val(TILE_DIM_2))
        end
        synchronize()
    end

    t_D2_naive_v2 = t_D2_naive_v2 * 1000 / iter_times
    @show t_D2_naive_v2
    # End evaluating D2_naive_v2

    t_A = @elapsed begin
        for _ in 1:iter_times
            matrix_free_A(idata,odata)
        end
        synchronize()
    end
    t_A = t_A * 1000 / iter_times
    @show t_A

    t_A_full_GPU = @elapsed begin
        for _ in 1:iter_times
            matrix_free_A_full_GPU(idata,odata)
        end
        synchronize()
    end
    t_A_full_GPU = t_A_full_GPU * 1000 / iter_times
    @show t_A_full_GPU

    # through_put = 2 * memsize / (t_D2_naive_v2/1000)
    through_put = 2 * memsize / (t_A_full_GPU/1000)
    @printf("through_put:      %f GiB / s\n", through_put)

end

function CG_GPU(b_reshaped_GPU,x_GPU)
    (Nx,Ny) = size(b_reshaped_GPU);
    odata = CuArray(zeros(Nx,Ny))
    matrix_free_A(x_GPU,odata);
    r_GPU = b_reshaped_GPU - odata;
    p_GPU = copy(r_GPU);
    # rsold_GPU = sum(r_GPU .* r_GPU)
    rsold_GPU = norm(r_GPU)^2
    Ap_GPU = CuArray(zeros(Nx,Ny))
    num_iter_steps = 0
    machine_eps = sqrt(eps(real(eltype(b_reshaped_GPU))))
    rel_tol = machine_eps * max(sqrt(rsold_GPU),1)
    # @show rel_tol
    norms = [sqrt(rsold_GPU)]
    # @show rsold_GPU
    for i in 1:Nx*Ny
        num_iter_steps += 1
        matrix_free_A(p_GPU,Ap_GPU);
        # sum(p_GPU .* Ap_GPU)
        alpha_GPU = rsold_GPU / (sum(p_GPU .* Ap_GPU))
        r_GPU .-= alpha_GPU .* Ap_GPU
        x_GPU .+= alpha_GPU .* p_GPU
        rsnew_GPU = norm(r_GPU)^2
        if sqrt(rsnew_GPU) < rel_tol
            break
        end
        ratio = (rsnew_GPU/rsold_GPU)
        p_GPU .= r_GPU .+ ratio .* p_GPU;
        rsold_GPU = rsnew_GPU
        append!(norms,sqrt(rsold_GPU))
        # if i < 20
        #     @show rsold_GPU
        # end
    end
    @show
    # @show num_iter_steps
    (num_iter_steps,norms[end])
end

function CG_GPU_v2(b_reshaped_GPU,x_GPU,odata,r_GPU,p_GPU,Ap_GPU)
    # preallocation version
    (Nx,Ny) = size(b_reshaped_GPU);
    matrix_free_A(x_GPU,odata);
    r_GPU .= b_reshaped_GPU .- odata;
    p_GPU .= r_GPU
    rsold_GPU = sum(r_GPU .* r_GPU)
    num_iter_steps = 0
    machine_eps = sqrt(eps(real(eltype(b_reshaped_GPU))))
    rel_tol = machine_eps * sqrt(rsold_GPU)
    norms = [sqrt(rsold_GPU)]
    for i in 1:Nx*Ny
        num_iter_steps += 1
        matrix_free_A(p_GPU,Ap_GPU);
        alpha_GPU = rsold_GPU / (sum(p_GPU .* Ap_GPU))
        r_GPU .-= alpha_GPU .* Ap_GPU
        x_GPU .+= alpha_GPU .* p_GPU
        # CUDA.CUBLAS.axpy!()
        rsnew_GPU = sum(r_GPU .* r_GPU)
        if sqrt(rsnew_GPU) < rel_tol
            break
        end
        p_GPU .= r_GPU .+ (rsnew_GPU/rsold_GPU) .* p_GPU;
        rsold_GPU = rsnew_GPU
        append!(norms,sqrt(rsold_GPU))
        # if i < 20
        #     @show rsold_GPU
        # end
    end
    # @show num_iter_steps
    (num_iter_steps,rel_tol,rsold_GPU)
end

function CG_full_GPU(b_reshaped_GPU,x_GPU)
    (Nx,Ny) = size(b_reshaped_GPU);
    odata = CuArray(zeros(Nx,Ny))
    # matrix_free_A_full_GPU(x_GPU,odata);
    matrix_free_A_full_GPU(x_GPU,odata)
    r_GPU = b_reshaped_GPU - odata;
    p_GPU = copy(r_GPU);
    rsold_GPU = sum(r_GPU .* r_GPU)
    # Ap_GPU = CUDA.zeros(Nx,Ny);
    Ap_GPU = CuArray(zeros(Nx,Ny))
    num_iter_steps = 0
    machine_eps = sqrt(eps(real(eltype(b_reshaped_GPU))))
    # rel_tol = machine_eps * max(sqrt(rsold_GPU),1)
    rel_tol = machine_eps
    norms = [sqrt(rsold_GPU)]
    for i in 1:Nx*Ny
    # for i in 1:20
        num_iter_steps += 1
        matrix_free_A_full_GPU(p_GPU,Ap_GPU);
        alpha_GPU = rsold_GPU / (sum(p_GPU .* Ap_GPU))
        # x_GPU = x_GPU + alpha_GPU * p_GPU;
        # r_GPU = r_GPU - alpha_GPU * Ap_GPU;
        r_GPU .-= alpha_GPU .* Ap_GPU
        x_GPU .+= alpha_GPU .* p_GPU
        # CUDA.CUBLAS.axpy!()
        rsnew_GPU = sum(r_GPU .* r_GPU)
        if sqrt(rsnew_GPU) < rel_tol
            break
        end
        p_GPU .= r_GPU .+ (rsnew_GPU/rsold_GPU) .* p_GPU;
        rsold_GPU = rsnew_GPU
        append!(norms,sqrt(rsold_GPU))
        # if i < 20
        #     @show rsold_GPU
        # end
    end
    # @show num_iter_steps
    (num_iter_steps,rel_tol,norms[end])
end

function CG_CPU(A,b,x)
    r = b - A * x;
    p = r;
    rsold = r' * r
    # Ap = spzeros(length(b))
    Ap = similar(b);

    num_iter_steps = 0
    # @show rsold
    for _ = 1:length(b)
    # for _ in 1:20
        num_iter_steps += 1
        mul!(Ap,A,p);
        alpha = rsold / (p' * Ap)
        x .= x .+ alpha * p;
        r .= r .- alpha * Ap;
        rsnew = r' * r
        if sqrt(rsnew) < sqrt(eps(real(eltype(b))))
              break
        end
        p = r + (rsnew / rsold) * p;
        rsold = rsnew
        # @show rsold
    end
    # @show num_iter_steps
    num_iter_steps
end


text:
"""
    intensity(ts, events, background, kappa, kernel)
Evaluate the intensity over a timegrid for some observed events with given Hawkes parameters.

# Arguments
* `ts::Array{<:Number,1}`: time grid to evalue the intensity over
* `events::Array{<:Number,1}`: Events of the process
* `background`: Background rate
* `kappa::Float64`: Kappa value
* `kernel::Function`: Kernel function

# Notes
* `kappa` must be between 0 and 1

"""
function intensity(ts::Array{<:Number, 1}, events::Array{<:Number, 1}, background, kappa::Float64, kernel::Function)
  output = Array{Number, 1}(undef, length(ts))
  for (i, timestamp) in enumerate(ts)
    output[i] = intensity(timestamp, events, background, kappa, kernel)
  end
  return output
end

function intensity(ts::Array{<:Number, 1}, events::Array{<:Number, 1}, background, kappa::Float64, kernel::Distributions.Distribution)
  kernel_pdf, kernel_cdf = generate_kernel_functions(kernel)
  intensity(ts, events, background, kappa, kernel_pdf)
end

function intensity(ts::Number, events::Array{<:Number, 1},
                   background::Number, kappa::Float64, kernel::Function)
  bgContrib = background
  bgContrib + kernel_contribution(ts, events, kappa, kernel)
end

function intensity(ts::Number, events::Array{<:Number, 1},
                   background::Function, kappa::Float64, kernel::Function)
  bgContrib = background(ts)
  bgContrib + kernel_contribution(ts, events, kappa, kernel)
end

function kernel_contribution(ts::Number, events::Array{<:Number, 1}, kappa::Float64, kernel::Function)
  evalEvents = events[events .< ts]
  s::Number = zero(ts)
  if !isempty(evalEvents)
    kernEval::Array{Number, 1} = kernel.(ts .- evalEvents)
    s = kappa * sum(kernEval)
  end
  return s
end

ext: jl


text:
using ABCDGraphGenerator
using Random

@info "Usage: julia graph_sampler.jl networkfile communityfile degreefile communitysizesfile mu|xi fraction isCL islocal [seed]"
@info "Example: julia graph_sampler.jl network.dat community.dat degrees.dat community_sizes.dat xi 0.2 true true 42"

networkfile = ARGS[1]
communityfile = ARGS[2]
degreefile = ARGS[3]
communitysizesfile = ARGS[4]
muxi = ARGS[5]
fraction = parse(Float64, ARGS[6])
isCL = parse(Bool, ARGS[7])
islocal = parse(Bool, ARGS[8])

muxi in ["mu","xi"] || throw(ArgumentError("only mu or xi names are allowed for"))
μ, ξ = nothing, nothing
if muxi == "mu"
    μ = fraction
else
    ξ = fraction
end
length(ARGS) == 9 &&  Random.seed!(parse(Int, ARGS[9]))

p = ABCDGraphGenerator.ABCDParams(parse.(Int, readlines(degreefile)),
                                  parse.(Int, readlines(communitysizesfile)),
                                  μ, ξ, isCL, islocal)

edges, clusters = ABCDGraphGenerator.gen_graph(p)

open(networkfile, "w") do io
    for (a, b) in sort!(collect(edges))
        println(io, a, "\t", b)
    end
end

open(communityfile, "w") do io
    for (i, c) in enumerate(clusters)
        println(io, i, "\t", c)
    end
end


text:
struct Problem{fType,S}

    f::fType  # Objective function

    LB::Vector{S}
    UB::Vector{S}

end

# If LB and UB are not specified, use default of +/- 1000
function Problem(objFunc, numVars)

    # Check that nDims > 0
    numVars > 0 || throw(ArgumentError("N must be greater than zero."))

    # Get types
    fType = typeof(objFunc)

    # Initialize lower and upper bounds
    LB = Vector{Int}(undef, numVars)
    UB = Vector{Int}(undef, numVars)
    @inbounds for i in 1:numVars
        LB[i] = -1000
        UB[i] = 1000
    end

    return Problem{fType,Int}(objFunc, LB, UB)
end

function Problem{S}(objFunc, LB, UB) where {S}

    # Check that LB and UB are the same length
    N = length(LB)
    if length(UB) != N
        throw(ArgumentError("Lower bounds and upper bound vectors must be the same length."))
    end

    # Check that LB[i] < UB[i] ∀ i
    @inbounds for i in 1:N
        LB[i] < UB[i] || throw(ArgumentError("LB[i] must be greater than UB[i] for all i."))
    end

    # Get f type
    fType = typeof(objFunc)

    # Initialize lower and upper bounds
    sLB = Vector{S}(LB)
    sUB = Vector{S}(UB)

    return Problem{fType,S}(objFunc, sLB, sUB)
end

function Problem(objFunc, LB::AbstractArray{S}, UB::AbstractArray{S}) where {S}
    return Problem{S,N}(objFunc, LB, UB)
end


text:
function loadInput(inputFile::AbstractString)::Vector{Any}
    raw = readchomp(inputFile)

    str = "[" * replace(replace(raw, "[" => "Any["), "\n" => ",\n") * "]"

    return include_string(Main, str)
end

function explodeSnail!(root, num, depth, indices=[])
    if depth >= 4 && isa(num[1], Number) && isa(num[2], Number)
        # Has to be a pair if both are numbers
        # Explode for depth >= 4

        pairValues  = [num[1], num[2]]

        # Get the containing pair
        pair = root
        for i in 1:length(indices)-1
            pair = pair[indices[i]]
        end

        # Replace with 0
        pair[indices[end]] = 0

        # Add left and right
        for (i, j) in [(1, 2), (2, 1)]
            last = findlast(x -> x==j, indices)

            if last !== nothing

                # Get the containing pair
                pair = root

                for i in 1:last-1
                    pair = pair[indices[i]]
                end

                if isa(pair[i], Number)
                    pair[i] += pairValues[i]
                else
                    pair = pair[i]
                    while true
                        if isa(pair[j], Number)
                            pair[j] += pairValues[i]
                            break
                        end
                        pair = pair[j]
                    end
                end
            end

        end

        return true
    end

    if !isa(num[1], Number)
        push!(indices, 1)
        if explodeSnail!(root, num[1], depth + 1, indices)
            return true
        end
    end

    if !isa(num[2], Number)
        push!(indices, 2)
        if explodeSnail!(root, num[2], depth + 1, indices)
            return true
        end
    end

    if !isempty(indices)
        pop!(indices)
    end

    return false
end

function explodeSnail!(root, depth=0)
    explodeSnail!(root, root, depth)
end

function splitSnail!(num)

    if isa(num[1], Number)
        if num[1] >= 10
            left, right = round(num[1]/2, RoundDown), round(num[1]/2, RoundUp)
            num[1] = Any[left, right]
            return true
        end
    end

    for i in 1:2
        if !isa(num[i], Number)
            if splitSnail!(num[i])
                return true
            end
        end
    end

    if isa(num[2], Number)
        if num[2] >= 10
            left, right = round(num[2]/2, RoundDown), round(num[2]/2, RoundUp)
            num[2] = Any[left, right]
            return true
        end
    end

    return false
end

function reduceSnail(num)
    while true
        explodeSnail!(num) && continue
        splitSnail!(num) && continue
        break
    end

    return num
end

function addSnails(left, right)
    return reduceSnail([left, right])
end

function magnitude(num)
    if isa(num, Number)
        return num
    else
        return 3*magnitude(num[1]) + 2*magnitude(num[2])
    end
end

#! Part 1

snailNumbers = loadInput("input.txt")
sum = snailNumbers[1]

for i in 2:length(snailNumbers)
    sum = addSnails(sum, snailNumbers[i])
end

magnitude(sum)

#! Part 2

snailNumbers = loadInput("input.txt")
maximum([magnitude(addSnails(deepcopy(snailNumbers[i]), deepcopy(snailNumbers[j]))) for i in 1:length(snailNumbers) for j in 1:length(snailNumbers) if i!=j])

max_stars_repo_path: Day 18/main.jl

max_issues_repo_path: Day 18/main.jl

max_forks_repo_path: Day 18/main.jl


text:
@testset "adt" begin
    @testset "adt List" begin
        # 1. define List
        @data List{T} begin
            Nil{T}
            Cons{T}(head :: T, tail :: List{T})
        end
        # 2. define interfaces
        len(xs::List{T}) where T = @match xs begin
            Nil{T}() => 0
            Cons{T}(_, tail) => 1 + len(tail)
        end

        @test len(Nil{Any}()) == 0
        xs = Cons(3,Cons(2, Cons(1, Nil{Int}())))
        @test len(xs) == 3
    end
    @testset "adt Arith" begin
        # 1. define Arith
        @data Arith begin
            Number(v :: Int)
            Minus(fst :: Arith, snd :: Arith)
            Add(fst :: Arith, snd :: Arith)
            Mult(fst :: Arith, snd :: Arith)
            Divide(fst :: Arith, snd :: Arith)
        end
        # 2. define interfaces
        function eval_arith(arith :: Arith)
            @match arith begin
                Number(v)       => v
                Add(fst, snd) => eval_arith(fst) + eval_arith(snd)
                Minus(fst, snd) => eval_arith(fst) - eval_arith(snd)
                Mult(fst, snd)   => eval_arith(fst) * eval_arith(snd)
                Divide(fst, snd) => eval_arith(fst) / eval_arith(snd)
            end
        end

        @test eval_arith(
            Add(Number(1),
                Minus(Number(2),
                    Divide(Number(20),
                            Mult(Number(2),
                                Number(5)))))) == 1
    end

end

# for type definition not allowed in the local scope, define it at top level.
abstract type A end

@testset "@case" begin
    @case C{T}(a :: Int, b::T)
    @case D(a, b)

    @case E <: A

    @test E <: A
    @test fieldnames(D) == (:a, :b)
    @test_throws MethodError C(3.0, :abc)

end


text:
# RSA protocol (RFC-3447) implementation in Julia
# Copyright (C): Jaka Smrekar (vinctux) <vinctux@outlook.com>, 2016
#
# PyCrypto (https://github.com/dlitz/pycrypto) as Pythonic reference

using Krypto
using SHA

abstract RSAKey

type RSAPubKey <: RSAKey
    n::BigInt             # RSA modulus
    e::BigInt             # RSA public component
    o::Array{UInt8, 1}    # Octet array representation of RSA key (n, e)
    hex::AbstractString   # Hex representation of RSA key (n, e)
end
RSAPubKey() = RSAPubKey(0, 65537, Array{UInt8, 1}(), "")

type RSAPrivKey <: RSAKey
    n::BigInt             # RSA modulus
    d::BigInt             # RSA private component
    o::Array{UInt8, 1}    # Octet array representation of RSA key (n, e)
    hex::AbstractString   # Hex representation of RSA key (n, e)
end
RSAPrivKey() = RSAPrivKey(0, 65537, Array{UInt8, 1}(), "")

# RFC-3447/4.1
function IP2OS{T<:Integer, L<:Integer}(x::T, xLen::L)
    #if x >= 256^BigInt(xLen) error("Integer too large.") end
    X = Array{UInt8, 1}()
    while x > 256 append!(X, x % 256); x = div(x, 256) end
    while length(X) < xLen append!(X, 0o0) end
    return reverse(X)
end
IP2OS{T<:Integer}(x::T) = IP2OS(x, div(bitsize(x), 8))

# RFC-3447/4.2
function OS2IP(X::Array{UInt8, 1})
    x = BigInt(0)
    for i in 1:length(X) x += X[i] * 256^BigInt(length(X) - i) end
    return x
end

function OS2STR(X::Array{UInt8, 1})
    r = ""
    for i in 1:length(X) r *= string(Char('\0' + d)) end
    return r
end

# RFC-3447/5.1.1
function RSAEP(pubkey::RSAPubKey, m::Integer)
    if m > pubkey.n - 1 || m < 0 error("Message representative out of range.") end
    return powermod(m, pubkey.e, pubkey.n)
end

# RFC-3447/5.1.2
function RSADP(privkey::RSAPrivKey, m::Integer)
    if m > privkey.n - 1 || m < 0 error("Message representative out of range.") end
    return powermod(m, privkey.d, privkey.n)
end

# RFC-3447/5.2.1
function RSASP1(privkey::RSAPrivKey, m::Integer)
    if m > privkey.n - 1 || m < 0 error("Message representative out of range.") end
    return powermod(m, privkey.d, privkey.n)
end

# RFC-3447/5.2.2
function RSAVP1(pubkey::RSAPubKey, m::Integer)
    if m > pubkey.n - 1 || m < 0 error("Signature representative out of range.") end
    return powermod(m, pubkey.e, pubkey.n)
end

# FIXME: RSAES_PKCS1_V1_5_ENCRYPT not working properly.
# RFC-3447/7.2.1
function RSAES_PKCS1_V1_5_ENCRYPT(pubkey::RSAPubKey, M::Array{UInt8, 1})
    k, mLen = length(IP2OS(pubkey.n)), length(M)
    mLen > k - 11 && error("Message too long.")
    PS = [csrand(0o1:0o255) for i in 1:(k - mLen - 3)]
    EM = vcat([0o0, 0o2], PS, [0o0], M)
    m, c, C = OS2IP(EM), RSAEP(pubkey, m), IP2OS(c, k)
    return C
end

# FIXME: RSAES_PKCS1_V1_5_DECRYPT not working properly.
# RFC-3447/7.2.2
function RSAES_PKCS1_V1_5_DECRYPT(privkey::RSAPrivKey, C::Array{UInt8, 1})
    k = length(IP2OS(privkey.n))
    length(C) != k && error("Decryption error.")
    EM = IP2OS(RSADP(privkey, OS2IP(C)), k)
    (EM[1] != 0o0 || EM[2] != 0o2) && error("Decryption error.")
    i, PS = 3, Array{UInt8, 1}()
    while EM[i] != 0o0 vcat(PS, EM[i]); i += 1 end
    i < 8 && error("Decryption error.")
    M = EM[(i + 1):end]
    return M
end

# FIXME: RSASSA_PKCS1_V1_5_SIGN not working properly.
# RFC-3447/8.2.1
function RSASSA_PKCS1_V1_5_SIGN(privkey::RSAPrivKey, M::Array{UInt8, 1})
    EM = EMSA_PKCS1_V1_5_ENCODE(M, length(IP2OS(privkey.n)))
    return IP2OS(RSASP1(privkey, OS2IP(EM)))
end

# FIXME: RSASSA_PKCS1_V1_5_VERIFY not working properly.
# RFC-3447/8.2.2
function RSASSA_PKCS1_V1_5_VERIFY(pubkey::RSAPubKey, M::Array{UInt8, 1}, S::Array{UInt8, 1})
    k = length(IP2OS(pubkey.n))
    length(S) != k && error("Invalid signature.")
    EM1, EM2 = IP2OS(RSAVP1(pubkey, OS2IP(S)), k), EMSA_PKCS1_V1_5_ENCODE(M, k)
    return EM1 == EM2
end

# TODO: Use DER/BER encoding standard before concating.
# RFC-3447/9.1
EMSA_PKCS1_V1_5_ENCODE{C<:SHA.SHA_CTX}(M::Array{UInt8, 1}, k::Integer, H::C = SHA.SHA3_CTX) =
    vcat([0o0, 0o1], [0o377 for i in 1:(k - length(M) - 3)], [0o0], update!(H, M))

function RSAKeyGen{T<:Integer}(b::T, e::T = 65537, newe::Bool = false, ERR::Int8 = Int8(75))
    if (b - 1024) % 256 != 0 error("Number of bits should be in form b = 1024 + 256x.") end
    p = q = BigInt(1)
    while bitsize(BigInt(p * q)) < b
        p = BigInt(csprime(div(b, 2), ERR))
        q = BigInt(csprime(div(b, 2), ERR))
    end
    N = BigInt(p * q)
    phi = BigInt(N - (p + q - 1))
    if newe
        e = BigInt(csrand(1:phi))
        while gcd(e, phi) != 1 e = csrand(1:phi) end
    end
    PUB  = RSAPubKey(N, e, Array{UInt8, 1}(), "")
    PRIV = RSAPrivKey(N, invmod(e, phi), Array{UInt8, 1}(), "")
    return PUB, PRIV
end

## SHORTHANDS ##
RSAEncrypt(pubkey::RSAPubKey,   M::Array{UInt8, 1}) = RSAES_PKCS1_V1_5_ENCRYPT(pubkey, M)
RSADecrypt(privkey::RSAPrivKey, C::Array{UInt8, 1}) = RSAES_PKCS1_V1_5_DECRYPT(privkey, C)
RSASign(privkey::RSAPrivKey,    M::Array{UInt8, 1}) = RSASSA_PKCS1_V1_5_SIGN(privkey, M)
RSAVerify(pubkey::RSAPubKey,    M::Array{UInt8, 1}, S::Array{UInt8, 1}) = RSASSA_PKCS1_V1_5_VERIFY(pubkey, M, S)

encrypt(::Type{Krypto.RSA}, K::RSAPubKey,  M::Array{UInt8, 1}) = RSAEncrypt(K, M)
decrypt(::Type{Krypto.RSA}, K::RSAPrivKey, C::Array{UInt8, 1}) = RSADecrypt(K, C)
sign(::Type{Krypto.RSA},    K::RSAPrivKey, M::Array{UInt8, 1}) = RSASign(K, M)
verify(::Type{Krypto.RSA},  K::RSAPubKey,  C::Array{UInt8, 1}) = RSAVerify(K, C)


text:
require("bm3d")

using Images, ImageView
using NumericExtensions
using BM3D

function getPsnr(img::Matrix{Float64},imgClean::Matrix{Float64})
	return 20.0*log10(255.0*sqrt(length(imgClean)/sumsqdiff(imgClean,img)))
end

# Convert test image to grayscale
img = imread("lena.tiff")
tmp = mean(convert(Array{Float64},img.data),1)
tmp = reshape(tmp,512,512)'
img = convert(Image,tmp)
img.properties["limits"] = (0.0,255.0)

# Noisy image
imgn = deepcopy(img)
sig = 15.0
imgn.data += sig*randn(512,512)

# Perform BM3D denoising
@time imgBasic = bm3d_thr(imgn.data,sig)
@time imgOut = bm3d_wie(imgn.data,imgBasic,sig)

imgBasic = convert(Image,imgBasic)
imgBasic.properties["limits"] = (0.0,255.0)

imgOut = convert(Image,imgOut)
imgOut.properties["limits"] = (0.0,255.0)

errNoisy = getPsnr(imgn.data,img.data)
errBasic = getPsnr(img.data,imgBasic.data)
errFinal = getPsnr(img.data,imgOut.data)

println("PSNR noisy: ", errNoisy, " dB")
println("PSNR basic: ", errBasic, " dB")
println("PSNR denoised: ", errFinal, " dB")

#display([imgn imgBasic])


text:
## Compound objects

export Line
function Line(p₁::RealVector, p₂::RealVector, r::Real; infty=10^5)
    v = normalize(p₂-p₁)
    return Cylinder(p₁-infty*v,p₂+infty*v,r)
end

export HalfLine
function HalfLine(p₁::RealVector, p₂::RealVector, r::Real; infty=10^5)
    v = normalize(p₂-p₁)
    return Cylinder(p₁,p₂+infty*v,r)
end

export Arrow
function Arrow(end1::RealVector, end2::RealVector, r::Real)
    n = normalize(end2-end1)
    end3 = end2-6*r*n
    return csgUnion(Cylinder(end1,end3,r),Cone(end3,end2,2*r))
end

function Blocks³(p₁::RealVector, p₂::RealVector, p₃::RealVector, thickness::Float64)
    center = Circumcenter(p₁,p₂,p₃)
    e₃ = NormalVector(p₁,p₂,p₃)
    R = norm(p₁-center)
    e₁ = normalize(p₁-center)
    e₂ = cross(e₃,e₁)
    box = Box([R+thickness,R+thickness,thickness],-[R+thickness,0,thickness])
    box1 = AffineTransform(box,hcat(e₁,e₂,e₃),center)
    e₁ = -normalize(p₃-center)
    e₂ = cross(e₃,e₁)
    box2 = AffineTransform(box,hcat(e₁,e₂,e₃),center)
    if det(hcat(p₁-center,p₃-center,e₃)) > 0
        return csgIntersection(box1,box2)
    else
        return csgMerge(box1,box2)
    end
end

export Arc
function Arc(p₁::RealVector, p₂::RealVector, p₃::RealVector, r::Float64; ε=10^(-4))
    α₁,α₂,α₃ = Angles(p₁,p₂,p₃)
    if (α₁>π-ε)
        return csgUnion(HalfLine(p₁,p₁+(p₂-p₃),r),HalfLine(p₃,p₃-(p₂-p₃),r))
    elseif (α₂>π-ε)
        return Cylinder(p₁,p₃,r)
    elseif (α₃>π-ε)
        return csgUnion(HalfLine(p₁,p₁+(p₁-p₂),r),HalfLine(p₃,p₃-(p₁-p₂),r))
    else
        return csgBound(csgIntersection(Torus(p₁,p₂,p₃,r),Blocks³(p₁,p₂,p₃,2r)),Torus(p₁,p₂,p₃,r))
    end
end

export RoundedBox
function RoundedBox(vertex1::RealVector, vertex2::RealVector, radius::Real)
    if radius == 0
        return Box(vertex1, vertex2)
    elseif radius < 0
        error("radius must be non-negative")
    end
    x=vertex1[1], vertex2[1]
    y=vertex1[2], vertex2[2]
    z=vertex1[3], vertex2[3]
    x_min, x_max=minimum(x), maximum(x)
    y_min, y_max=minimum(y), maximum(y)
    z_min, z_max=minimum(z), maximum(z)

    v_min=[x_min, y_min, z_min] .+ radius
    v_max=[x_max, y_max, z_max] .- radius

    x_r=v_min[1], v_max[1]
    y_r=v_min[2], v_max[2]
    z_r=v_min[3], v_max[3]

    spheres=[Sphere([x_r[i], y_r[j], z_r[k]], radius) for i in 1:2, j in 1:2, k in 1:2]

    cylinders_x=[Cylinder([x_r[1], y_r[j], z_r[k]], [x_r[2], y_r[j], z_r[k]], radius) for j in 1:2, k in 1:2]
    cylinders_y=[Cylinder([x_r[i], y_r[1], z_r[k]], [x_r[i], y_r[2], z_r[k]], radius) for i in 1:2, k in 1:2]
    cylinders_z=[Cylinder([x_r[i], y_r[j], z_r[1]], [x_r[i], y_r[j], z_r[2]], radius) for i in 1:2, j in 1:2]

    box_x=Box([x_min, y_r[1], z_r[1]], [x_max, y_r[2], z_r[2]])
    box_y=Box([x_r[1], y_min, z_r[1]], [x_r[2], y_max, z_r[2]])
    box_z=Box([x_r[1], y_r[1], z_min], [x_r[2], y_r[2], z_max])

    return csgMerge(spheres..., cylinders_x..., cylinders_y..., cylinders_z..., box_x, box_y, box_z)
end

max_stars_repo_name: hyrodium/JuliRay

max_issues_repo_name: hyrodium/JuliRay

max_forks_repo_name: hyrodium/JuliRay


text:
@testset "fglm finite field" begin

    R, (x, y) = PolynomialRing(GF(2^31-1), ["x","y"], ordering=:degrevlex)

    noon = Groebner.change_ordering(Groebner.noonn(2, ground=GF(2^31-1)), :degrevlex)
    gb = Groebner.groebner(noon, ordering=:degrevlex)
    x1, x2 = gens(parent(first(Groebner.fglm(gb))))
    @test Groebner.fglm(gb) == [x2^5 + 1952257860*x2^4 + 1288490186*x2^3 + 2*x2^2 + 1028839893*x2 + 644245093,
                                x1 + x2^4 + 1952257860*x2^3 + 644245093*x2^2 + x2 + 1952257860]

    gb = Groebner.groebner(noon, ordering=:deglex)
    x1, x2 = gens(parent(first(Groebner.fglm(gb))))
    @test Groebner.fglm(gb) == [x2^5 + 1952257860*x2^4 + 1288490186*x2^3 + 2*x2^2 + 1028839893*x2 + 644245093,
                                x1 + x2^4 + 1952257860*x2^3 + 644245093*x2^2 + x2 + 1952257860]

    gb = Groebner.groebner(noon, ordering=:lex)
    x1, x2 = gens(parent(first(Groebner.fglm(gb))))
    @test Groebner.fglm(gb) == [x2^5 + 1952257860*x2^4 + 1288490186*x2^3 + 2*x2^2 + 1028839893*x2 + 644245093,
                                x1 + x2^4 + 1952257860*x2^3 + 644245093*x2^2 + x2 + 1952257860]

end

@testset "fglm rationals" begin
    R, (x, y, z) = PolynomialRing(QQ, ["x","y","z"], ordering=:degrevlex)

    gb = [  z^2 - (9//490)*y - 201//980*z + 13//980,
            y*z - 4//35*y + 2//35*z - 1//35,
            y^2 - 1//10*y - 6//5*z + 1//10,
            x^2 - 4//5*y + 2//5*z - 1//5]

    x, y, z = gens(parent(first(Groebner.fglm(gb))))
    @test Groebner.fglm(gb) == [
                                z^3 - 313//980*z^2 + 37//980*z - 1//490,
                                y - 490//9*z^2 + 67//6*z - 13//18,
                                x^2 - 392//9*z^2 + 28//3*z - 7//9,
                                ]

end


text:
@doc raw"""
    quasi_Newton(M, F, gradF, x)

Perform a quasi Newton iteration for `F` on the manifold `M` starting
in the point `x` using a retraction ``R`` and a vector transport ``T``

The ``k``th iteration consists of
1. Compute the search direction ``η_k = -\mathcal{B}_k [\operatorname{grad}f (x_k)]`` or solve ``\mathcal{H}_k [η_k] = -\operatorname{grad}f (x_k)]``.
2. Determine a suitable stepsize ``α_k`` along the curve ``\gamma(α) = R_{x_k}(α η_k)`` e.g. by using [`WolfePowellLineseach`](@ref).
3. Compute ``x_{k+1} = R_{x_k}(α_k η_k)``.
4. Define ``s_k = T_{x_k, α_k η_k}(α_k η_k)`` and ``y_k = \operatorname{grad}f(x_{k+1}) - T_{x_k, α_k η_k}(\operatorname{grad}f(x_k))``.
5. Compute the new approximate Hessian ``H_{k+1}`` or its inverse ``B_k``.

# Input
* `M` – a manifold ``\mathcal{M}``.
* `F` – a cost function ``F : \mathcal{M} →ℝ`` to minimize.
* `gradF`– the gradient ``\operatorname{grad}F : \mathcal{M} →T_x\mathcal M`` of ``F``.
* `x` – an initial value ``x ∈ \mathcal{M}``.

    stopping_criterion::StoppingCriterion=StopWhenAny(
        StopAfterIteration(max(1000, memory_size)), StopWhenGradientNormLess(10^(-6))
    ),
    return_options=false,

# Optional
* `basis` – (`DefaultOrthonormalBasis()`) basis within the tangent space(s) to represent the
  Hessian (inverse).
* `cautious_update` – (`false`) – whether or not to use
  a [`QuasiNewtonCautiousDirectionUpdate`](@ref)
* `cautious_function` – (`(x) -> x*10^(-4)`) – a monotone increasing function that is zero
  at 0 and strictly increasing at 0 for the cautious update.
* `direction_update` – ([`InverseBFGS`](@ref)`()`) the update rule to use.
* `evaluation` – ([`AllocatingEvaluation`](@ref)) specify whether the gradient works by
   allocation (default) form `gradF(M, x)` or [`MutatingEvaluation`](@ref) in place, i.e.
   is of the form `gradF!(M, X, x)`.
* `initial_operator` – (`Matrix{Float64}(I,n,n)`) initial matrix to use die the
  approximation, where `n=manifold_dimension(M)`, see also `scale_initial_operator`.
* `memory_size` – (`20`) limited memory, number of ``s_k, y_k`` to store. Set to a negative
  value to use a full memory representation
* `retraction_method` – (`ExponentialRetraction()`) a retraction method to use, by default
  the exponential map.
* `scale_initial_operator` - (`true`) scale initial operator with
  ``\frac{⟨s_k,y_k⟩_{x_k}}{\lVert y_k\rVert_{x_k}}`` in the computation
* `step_size` – ([`WolfePowellLineseach`](@ref)`(retraction_method, vector_transport_method)`)
  specify a [`Stepsize`](@ref).
* `stopping_criterion` - (`StopWhenAny(StopAfterIteration(max(1000, memory_size)), StopWhenGradientNormLess(10^(-6))`)
  specify a [`StoppingCriterion`](@ref)
* `vector_transport_method` – (`ParallelTransport()`) a vector transport to use.
* `return_options` – (`false`) – specify whether to return just the result `x` (default) or the complete [`Options`](@ref), e.g. to access recorded values. if activated, the extended result, i.e. the

# Output
* `x_opt` – the resulting (approximately critical) point of the quasi–Newton method
OR
* `options` – the options returned by the solver (see `return_options`)
"""
function quasi_Newton(
    M::AbstractManifold, F::Function, gradF::G, x::P; kwargs...
) where {P,G}
    x_res = allocate(x)
    copyto!(M, x_res, x)
    return quasi_Newton!(M, F, gradF, x_res; kwargs...)
end
@doc raw"""
    quasi_Newton!(M, F, gradF, x; options...)

Perform a quasi Newton iteration for `F` on the manifold `M` starting
in the point `x` using a retraction ``R`` and a vector transport ``T``.

# Input
* `M` – a manifold ``\mathcal{M}``.
* `F` – a cost function ``F: \mathcal{M} →ℝ`` to minimize.
* `gradF`– the gradient ``\operatorname{grad}F : \mathcal{M} → T_x\mathcal M`` of ``F``.
* `x` – an initial value ``x ∈ \mathcal{M}``.

For all optional parameters, see [`quasi_Newton`](@ref).
"""
function quasi_Newton!(
    M::AbstractManifold,
    F::Function,
    gradF::G,
    x::P;
    cautious_update::Bool=false,
    cautious_function::Function=x -> x * 10^(-4),
    retraction_method::AbstractRetractionMethod=ExponentialRetraction(),
    vector_transport_method::AbstractVectorTransportMethod=ParallelTransport(),
    basis::AbstractBasis=DefaultOrthonormalBasis(),
    direction_update::AbstractQuasiNewtonUpdateRule=InverseBFGS(),
    evaluation::AbstractEvaluationType=AllocatingEvaluation(),
    memory_size::Int=20,
    initial_operator::AbstractMatrix=Matrix{Float64}(
        I, manifold_dimension(M), manifold_dimension(M)
    ),
    scale_initial_operator::Bool=true,
    step_size::Stepsize=WolfePowellLineseach(retraction_method, vector_transport_method),
    stopping_criterion::StoppingCriterion=StopWhenAny(
        StopAfterIteration(max(1000, memory_size)), StopWhenGradientNormLess(10^(-6))
    ),
    return_options=false,
    kwargs...,
) where {P,G}
    if memory_size >= 0
        local_dir_upd = QuasiNewtonLimitedMemoryDirectionUpdate(
            direction_update,
            zero_vector(M, x),
            memory_size;
            scale=scale_initial_operator,
            vector_transport_method=vector_transport_method,
        )
    else
        local_dir_upd = QuasiNewtonMatrixDirectionUpdate(
            direction_update,
            basis,
            initial_operator;
            scale=scale_initial_operator,
            vector_transport_method=vector_transport_method,
        )
    end
    if cautious_update == true
        local_dir_upd = QuasiNewtonCautiousDirectionUpdate(
            local_dir_upd; θ=cautious_function
        )
    end

    p = GradientProblem(M, F, gradF; evaluation=evaluation)
    o = QuasiNewtonOptions(
        x,
        get_gradient(p, x),
        local_dir_upd,
        stopping_criterion,
        step_size;
        retraction_method=retraction_method,
        vector_transport_method=vector_transport_method,
    )

    o = decorate_options(o; kwargs...)
    resultO = solve(p, o)

    if return_options
        return resultO
    else
        return get_solver_result(resultO)
    end
end

function initialize_solver!(::GradientProblem, ::QuasiNewtonOptions) end

function step_solver!(p::GradientProblem, o::QuasiNewtonOptions, iter)
    o.gradient = get_gradient(p, o.x)
    η = o.direction_update(p, o)
    α = o.stepsize(p, o, iter, η)
    x_old = deepcopy(o.x)
    retract!(p.M, o.x, o.x, α * η, o.retraction_method)
    β = locking_condition_scale(
        p.M, o.direction_update, x_old, α * η, o.x, o.vector_transport_method
    )
    vector_transport_to!(
        p.M, o.sk, x_old, α * η, o.x, get_update_vector_transport(o.direction_update)
    )
    vector_transport_to!(
        p.M,
        o.gradient,
        x_old,
        o.gradient,
        o.x,
        get_update_vector_transport(o.direction_update),
    )
    o.yk = get_gradient(p, o.x) / β - o.gradient
    update_hessian!(o.direction_update, p, o, x_old, iter)
    return o
end

function locking_condition_scale(
    M::AbstractManifold, ::AbstractQuasiNewtonDirectionUpdate, x_old, v, x, vt
)
    return norm(M, x_old, v) / norm(M, x, vector_transport_to(M, x_old, v, x, vt))
end

@doc raw"""
    update_hessian!(d, p, o, x_old, iter)

update the hessian wihtin the [`QuasiNewtonOptions`](@ref) `o` given a [`Problem`](@ref) `p`
as well as the an [`AbstractQuasiNewtonDirectionUpdate`](@ref) `d` and the last iterate `x_old`.
Note that the current (`iter`th) iterate is already stored in `o.x`.

See also [`AbstractQuasiNewtonUpdateRule`](@ref) for the different rules that are available
within `d`.
"""
update_hessian!(d::AbstractQuasiNewtonDirectionUpdate, ::Any, ::Any, ::Any, ::Any)

function update_hessian!(
    d::QuasiNewtonMatrixDirectionUpdate{InverseBFGS}, p, o, x_old, iter
)
    update_basis!(d.basis, p.M, x_old, o.x, d.vector_transport_method)
    yk_c = get_coordinates(p.M, o.x, o.yk, d.basis)
    sk_c = get_coordinates(p.M, o.x, o.sk, d.basis)
    skyk_c = inner(p.M, o.x, o.sk, o.yk)
    if iter == 1 && d.scale == true
        d.matrix = skyk_c / inner(p.M, o.x, o.yk, o.yk) * d.matrix
    end
    d.matrix =
        (I - sk_c * yk_c' / skyk_c) * d.matrix * (I - yk_c * sk_c' / skyk_c) +
        sk_c * sk_c' / skyk_c
    return d
end

# BFGS update
function update_hessian!(d::QuasiNewtonMatrixDirectionUpdate{BFGS}, p, o, x_old, iter)
    update_basis!(d.basis, p.M, x_old, o.x, d.vector_transport_method)
    yk_c = get_coordinates(p.M, o.x, o.yk, d.basis)
    sk_c = get_coordinates(p.M, o.x, o.sk, d.basis)
    skyk_c = inner(p.M, o.x, o.sk, o.yk)
    if iter == 1 && d.scale == true
        d.matrix = inner(p.M, o.x, o.yk, o.yk) / skyk_c * d.matrix
    end
    d.matrix =
        d.matrix + yk_c * yk_c' / skyk_c -
        d.matrix * sk_c * sk_c' * d.matrix / dot(sk_c, d.matrix * sk_c)
    return d
end

# Inverse DFP update
function update_hessian!(d::QuasiNewtonMatrixDirectionUpdate{InverseDFP}, p, o, x_old, iter)
    update_basis!(d.basis, p.M, x_old, o.x, d.vector_transport_method)
    yk_c = get_coordinates(p.M, o.x, o.yk, d.basis)
    sk_c = get_coordinates(p.M, o.x, o.sk, d.basis)
    skyk_c = inner(p.M, o.x, o.sk, o.yk)
    if iter == 1 && d.scale == true
        d.matrix = inner(p.M, o.x, o.sk, o.sk) / skyk_c * d.matrix
    end
    d.matrix =
        d.matrix + sk_c * sk_c' / skyk_c -
        d.matrix * yk_c * yk_c' * d.matrix / dot(yk_c, d.matrix * yk_c)
    return d
end

# DFP update
function update_hessian!(d::QuasiNewtonMatrixDirectionUpdate{DFP}, p, o, x_old, iter)
    update_basis!(d.basis, p.M, x_old, o.x, d.vector_transport_method)
    yk_c = get_coordinates(p.M, o.x, o.yk, d.basis)
    sk_c = get_coordinates(p.M, o.x, o.sk, d.basis)
    skyk_c = inner(p.M, o.x, o.sk, o.yk)
    if iter == 1 && d.scale == true
        d.matrix = skyk_c / inner(p.M, o.x, o.sk, o.sk) * d.matrix
    end
    d.matrix =
        (I - yk_c * sk_c' / skyk_c) * d.matrix * (I - sk_c * yk_c' / skyk_c) +
        yk_c * yk_c' / skyk_c
    return d
end

# Inverse SR-1 update
function update_hessian!(
    d::QuasiNewtonMatrixDirectionUpdate{InverseSR1}, p, o, x_old, ::Int
)
    update_basis!(d.basis, p.M, x_old, o.x, d.vector_transport_method)
    yk_c = get_coordinates(p.M, o.x, o.yk, d.basis)
    sk_c = get_coordinates(p.M, o.x, o.sk, d.basis)

    # computing the new matrix which represents the approximating operator in the next iteration
    srvec = sk_c - d.matrix * yk_c
    if d.update.r < 0 || abs(dot(srvec, yk_c)) >= d.update.r * norm(srvec) * norm(yk_c)
        d.matrix = d.matrix + srvec * srvec' / (srvec' * yk_c)
    end
    return d
end

# SR-1 update
function update_hessian!(d::QuasiNewtonMatrixDirectionUpdate{SR1}, p, o, x_old, ::Int)
    update_basis!(d.basis, p.M, x_old, o.x, d.vector_transport_method)
    yk_c = get_coordinates(p.M, o.x, o.yk, d.basis)
    sk_c = get_coordinates(p.M, o.x, o.sk, d.basis)

    # computing the new matrix which represents the approximating operator in the next iteration
    srvec = yk_c - d.matrix * sk_c
    if d.update.r < 0 || abs(dot(srvec, sk_c)) >= d.update.r * norm(srvec) * norm(sk_c)
        d.matrix = d.matrix + srvec * srvec' / (srvec' * sk_c)
    end
    return d
end

# Inverse Broyden update
function update_hessian!(
    d::QuasiNewtonMatrixDirectionUpdate{InverseBroyden}, p, o, x_old, ::Int
)
    update_basis!(d.basis, p.M, x_old, o.x, d.vector_transport_method)
    yk_c = get_coordinates(p.M, o.x, o.yk, d.basis)
    sk_c = get_coordinates(p.M, o.x, o.sk, d.basis)
    skyk_c = inner(p.M, o.x, o.sk, o.yk)
    ykBkyk_c = yk_c' * d.matrix * yk_c

    φ = update_broyden_factor!(d, sk_c, yk_c, skyk_c, ykBkyk_c, d.update.update_rule)
    # computing the new matrix which represents the approximating operator in the next iteration
    d.matrix =
        d.matrix - (d.matrix * yk_c * yk_c' * d.matrix) / ykBkyk_c +
        (sk_c * sk_c') / skyk_c +
        φ * ykBkyk_c * (sk_c / skyk_c - (d.matrix * yk_c) / ykBkyk_c) *
        (sk_c / skyk_c - (d.matrix * yk_c) / ykBkyk_c)'
    return d
end

# Broyden update
function update_hessian!(d::QuasiNewtonMatrixDirectionUpdate{Broyden}, p, o, x_old, ::Int)
    update_basis!(d.basis, p.M, x_old, o.x, d.vector_transport_method)
    yk_c = get_coordinates(p.M, o.x, o.yk, d.basis)
    sk_c = get_coordinates(p.M, o.x, o.sk, d.basis)
    skyk_c = inner(p.M, o.x, o.sk, o.yk)
    skHksk_c = sk_c' * d.matrix * sk_c

    φ = update_broyden_factor!(d, sk_c, yk_c, skyk_c, skHksk_c, d.update.update_rule)
    # computing the new matrix which represents the approximating operator in the next iteration
    d.matrix =
        d.matrix - (d.matrix * sk_c * sk_c' * d.matrix) / skHksk_c +
        (yk_c * yk_c') / skyk_c +
        φ * skHksk_c * (yk_c / skyk_c - (d.matrix * sk_c) / skHksk_c) *
        (yk_c / skyk_c - (d.matrix * sk_c) / skHksk_c)'
    return d
end

function update_broyden_factor!(d, sk_c, yk_c, skyk_c, skHksk_c, s::Symbol)
    return update_broyden_factor!(d, sk_c, yk_c, skyk_c, skHksk_c, Val(s))
end

function update_broyden_factor!(d, ::Any, ::Any, ::Any, ::Any, ::Val{:constant})
    return d.update.φ
end

function update_broyden_factor!(d, ::Any, yk_c, skyk_c, skHksk_c, ::Val{:Davidon})
    yk_c_c = d.matrix \ yk_c
    ykyk_c_c = yk_c' * yk_c_c
    if skyk_c <= 2 * (skHksk_c * ykyk_c_c) / (skHksk_c + ykyk_c_c)
        return d.update.φ =
            (skyk_c * (ykyk_c_c - skyk_c)) / (skHksk_c * ykyk_c_c - skyk_c^2)
    else
        return d.update.φ = skyk_c / (skyk_c - skHksk_c)
    end
end

function update_broyden_factor!(d, sk_c, ::Any, skyk_c, ykBkyk_c, ::Val{:InverseDavidon})
    sk_c_c = d.matrix \ sk_c
    sksk_c_c = sk_c' * sk_c_c
    if skyk_c <= 2 * (ykBkyk_c * sksk_c_c) / (ykBkyk_c + sksk_c_c)
        return d.update.φ =
            (skyk_c * (sksk_c_c - skyk_c)) / (ykBkyk_c * sksk_c_c - skyk_c^2)
    else
        return d.update.φ = skyk_c / (skyk_c - ykBkyk_c)
    end
end

function update_basis!(
    b::AbstractBasis, ::AbstractManifold, ::P, ::P, ::AbstractVectorTransportMethod
) where {P}
    return b
end

function update_basis!(
    b::CachedBasis, M::AbstractManifold, x::P, y::P, m::AbstractVectorTransportMethod
) where {P}
    # transport all basis tangent vectors in the tangent space of the next iterate
    for v in b.data
        vector_transport_to!(M, v, x, v, y, m)
    end
    return b
end

# Cautious update
function update_hessian!(
    d::QuasiNewtonCautiousDirectionUpdate{U}, p, o, x_old, iter
) where {U<:AbstractQuasiNewtonDirectionUpdate}
    # computing the bound used in the decission rule
    bound = d.θ(norm(p.M, o.x, o.gradient))
    sk_normsq = norm(p.M, o.x, o.sk)^2
    if sk_normsq != 0 && (inner(p.M, o.x, o.sk, o.yk) / sk_normsq) >= bound
        update_hessian!(d.update, p, o, x_old, iter)
    end
    return d
end

# Limited-memory update
function update_hessian!(
    d::QuasiNewtonLimitedMemoryDirectionUpdate{U}, p, o, x_old, ::Int
) where {U<:InverseBFGS}
    (capacity(d.memory_s) == 0) && return d
    # only transport the first if it does not get overwritten at the end
    start = length(d.memory_s) == capacity(d.memory_s) ? 2 : 1

    for i in start:length(d.memory_s)
        # transport all stored tangent vectors in the tangent space of the next iterate
        vector_transport_to!(
            p.M, d.memory_s[i], x_old, d.memory_s[i], o.x, d.vector_transport_method
        )
        vector_transport_to!(
            p.M, d.memory_y[i], x_old, d.memory_y[i], o.x, d.vector_transport_method
        )
    end

    # add newest
    push!(d.memory_s, o.sk)
    push!(d.memory_y, o.yk)
    return d
end

# all Cautious Limited Memory
function update_hessian!(
    d::QuasiNewtonCautiousDirectionUpdate{QuasiNewtonLimitedMemoryDirectionUpdate{NT,T,VT}},
    p,
    o,
    x_old,
    iter,
) where {NT<:AbstractQuasiNewtonUpdateRule,T,VT<:AbstractVectorTransportMethod}
    # computing the bound used in the decission rule
    bound = d.θ(norm(p.M, x_old, get_gradient(p, x_old)))
    sk_normsq = norm(p.M, o.x, o.sk)^2

    # if the decission rule is fulfilled, the new sk and yk are added
    if sk_normsq != 0 && (inner(p.M, o.x, o.sk, o.yk) / sk_normsq) >= bound
        update_hessian!(d.update, p, o, x_old, iter)
    else
        # the stored vectores are just transported to the new tangent space, sk and yk are not added
        for i in 1:length(d.update.memory_s)
            vector_transport_to!(
                p.M,
                d.update.memory_s[i],
                x_old,
                d.update.memory_s[i],
                o.x,
                d.update.vector_transport_method,
            )
            vector_transport_to!(
                p.M,
                d.update.memory_y[i],
                x_old,
                d.update.memory_y[i],
                o.x,
                d.update.vector_transport_method,
            )
        end
    end
    return d
end

get_solver_result(o::QuasiNewtonOptions) = o.x


text:
module Results
    #
    # This module models a complete set of calculation results for all the members of a single household (tax due, benefits entitled to, marginal rate, etc.)
    # It also contains code to map from the Dictionary to Array representation of incomes.
    #
    # FIXME this badly needs rewritten to mirror the structure of the Household in ModelHousehold.jl
    #
    using Parameters: @with_kw
    using DataFrames
    using StaticArrays

    using ScottishTaxBenefitModel

    using .Definitions

    using .GeneralTaxComponents: RateBands

    using .ModelHousehold:
        BenefitUnit,
        BenefitUnits,
        Household,
        Person,
        Pid_Array,
        get_benefit_units

    using .STBIncomes
    using .RunSettings
    using .Utils:
        to_md_table

    export
        BenefitUnitResult,
        HouseholdResult,
        HousingResult,
        IndividualResult,
        ITResult,
        LMTCanApplyFor,
        LMTIncomes,
        LMTResults,
        LocalTaxes,
        NIResult,

        aggregate_tax,
        aggregate!,
        get_indiv_result,
        get_net_income,
        gross_pension_contributions,
        has_any,
        init_benefit_unit_result,
        init_household_result,
        map_incomes,
        to_string,
        total,
        tozero!

    @with_kw mutable struct BenefitCapResults{RT<:Real}
        cap :: RT = zero(RT)
        cap_benefits :: RT = zero(RT)
        reduction :: RT = zero(RT)
    end

    @with_kw mutable struct UCResults{RT<:Real}
        basic_conditions_satisfied :: Bool = false
        disqualified_on_capital :: Bool = false
        work_allowance :: RT = zero(RT)
        earnings_before_allowances :: RT = zero(RT)
        earned_income :: RT = zero(RT)
        other_income :: RT = zero(RT)
        tariff_income :: RT = zero(RT)
        standard_allowance  :: RT = zero(RT)
        child_element ::  RT = zero(RT)
        limited_capacity_for_work_activity_element ::  RT = zero(RT)
        carer_element ::  RT = zero(RT)
        childcare_costs :: RT = zero(RT)
        housing_element :: RT = zero(RT)
        other_housing_costs :: RT = zero(RT)
        assets :: RT = zero(RT)
        maximum  :: RT = zero(RT)
        total_income :: RT = zero(RT)
        transitional_protection :: RT = zero(RT)
        ctr :: RT = zero(RT)
        recipient :: BigInt = -1
    end

    @with_kw mutable struct LMTIncomes{RT<:Real}
        gross_earnings :: RT = zero(RT)
        net_earnings   :: RT = zero(RT)
        other_income   :: RT = zero(RT)
        total_income   :: RT = zero(RT)
        disregard :: RT = zero(RT)
        childcare :: RT = zero(RT)
        capital :: RT = zero(RT)
        tariff_income :: RT = zero(RT)
        disqualified_on_capital :: Bool = false
    end

    #
    @with_kw mutable struct LMTCanApplyFor
        ctc  :: Bool = false
        esa :: Bool = false
        hb  :: Bool = false
        is :: Bool = false
        jsa :: Bool = false
        pc  :: Bool = false
        sc  :: Bool = false
        wtc  :: Bool = false
        ctr :: Bool = false
    end

    @with_kw mutable struct LMTResults{RT<:Real}
        #ctc :: RT = zero(RT) XX
        #esa :: RT = zero(RT) x
        #hb  :: RT = zero(RT) x
        #is :: RT = zero(RT) x
        #jsa :: RT = zero(RT) x
        #pc  :: RT = zero(RT) x
        #sc   :: RT = zero(RT) x
        #wtc  :: RT = zero(RT)
        #ctr  :: RT = zero(RT) x

        ndds :: RT = zero(RT)
        mig  :: RT = zero(RT)
        mig_recipient :: BigInt = -1

        # total_benefits :: RT = zero(RT) #  hb and ctr

        # FIXME better name than MIG here
        # FIXME rename premia => premium everywhere here
        mig_premia :: RT = zero(RT)
        mig_allowances :: RT = zero(RT)
        mig_incomes = LMTIncomes{RT}()
        sc_incomes = LMTIncomes{RT}()

        ctc_elements :: RT = zero(RT)
        wtc_income  :: RT = zero(RT)
        wtc_elements :: RT = zero(RT)
        wtc_ctc_threshold :: RT = zero(RT)
        wtc_ctc_tapered_excess :: RT = zero(RT)
        wtc_recipient :: BigInt = -1
        cost_of_childcare :: RT = zero(RT)

        hb_passported :: Bool = false

        # TODO FIXME
        other_housing_costs :: RT = zero(RT)

        hb_premia :: RT = zero(RT)
        hb_allowances :: RT = zero(RT)
        hb_incomes = LMTIncomes{RT}()
        hb_eligible_rent :: RT = zero(RT)
        hb_recipient :: BigInt = -1

        ctr :: RT = zero(RT)
        ctr_passported :: Bool = false
        ctr_premia :: RT = zero(RT)
        ctr_allowances :: RT = zero(RT)
        ctr_incomes = LMTIncomes{RT}()
        ctr_eligible_amount :: RT = zero(RT)
        ctr_recipient :: BigInt = -1

        pc_premia :: RT = zero(RT)
        pc_allowances :: RT = zero(RT)
        pc_incomes = LMTIncomes{RT}()
        pc_recipient :: BigInt = -1

        premia = LMTPremiaSet()
        can_apply_for = LMTCanApplyFor()
    end

    @with_kw mutable struct NIResult{RT<:Real}
        above_lower_earnings_limit :: Bool = false
        # total_ni :: RT = 0.0
        class_1_primary    :: RT = 0.0
        class_1_secondary  :: RT = 0.0
        class_2   :: RT = 0.0
        class_3   :: RT = 0.0
        class_4   :: RT = 0.0
        assumed_gross_wage :: RT = 0.0
    end

    @with_kw mutable struct ITResult{RT<:Real}
        # total_tax :: RT = 0.0
        taxable_income :: RT = 0.0
        adjusted_net_income :: RT = 0.0
        total_income :: RT = 0.0
        allowance   :: RT = 0.0

        non_savings_tax :: RT = 0.0
        non_savings_band :: Integer = 0
        non_savings_income :: RT = 0.0
        non_savings_taxable :: RT = 0.0

        savings_tax :: RT = 0.0
        savings_band :: Integer = 0
        savings_income :: RT = 0.0
        savings_taxable :: RT = 0.0

        dividends_tax :: RT = 0.0
        dividend_band :: Integer = 0
        dividends_income :: RT = 0.0
        dividends_taxable :: RT = 0.0

        unused_allowance :: RT = 0.0
        mca :: RT = 0.0
        transferred_allowance :: RT = 0.0
        pension_eligible_for_relief :: RT = 0.0
        pension_relief_at_source :: RT = 0.0

        non_savings_thresholds :: RateBands = zeros(RT,0)
        savings_thresholds  :: RateBands = zeros(RT,0)
        dividend_thresholds :: RateBands = zeros(RT,0)

        savings_rates  :: RateBands = zeros(RT,0)
        dividend_rates :: RateBands = zeros(RT,0)

        personal_savings_allowance :: RT = 0.0
    end

    @with_kw mutable struct IndividualResult{RT<:Real}

       net_income :: RT = zero(RT)
       ni = NIResult{RT}()
       it = ITResult{RT}()
       metr :: RT = -12345.0
       replacement_rate :: RT = zero(RT)
       income = STBIncomes.make_a( RT );
    end

    function to_string( ir :: IndividualResult, depth=2 )::String
        s = to_md_table( ir, exclude=[:income], depth=depth )
        s *= #*repeat("#",depth)*"Incomes";
        s *= inctostr( ir.income )
        return s
    end

    @with_kw mutable struct BenefitUnitResult{RT<:Real} # FIXME make this a dynamic thing to match get_benefit_units
        income = STBIncomes.make_a( RT )
        net_income    :: RT = zero(RT)
        eq_net_income :: RT = zero(RT)
        legacy_mtbens = LMTResults{RT}()
        uc = UCResults{RT}()
        bencap = BenefitCapResults{RT}()
        other_benefits  :: RT = zero(RT)
        pers = Dict{BigInt,IndividualResult{RT}}() # FIXME name change to `people`
        adults = Pid_Array()
        route :: LegacyOrUC = legacy_bens
    end

    function to_string( br :: BenefitUnitResult, depth=1 )::String
        s = to_md_table( br, exclude=[:pers,:adults,:income], depth=depth )
        s *= "#### Benefit Unit Incomes "
        s *= inctostr(br.income)
        for pid in sort(collect(keys(br.pers)))
            s *= "### Individual Result; Person $(pid)"
            s *= to_string( br.pers[pid] )
        end
        return s
    end

    function has_any( bur :: BenefitUnitResult, things ... ) :: Bool
        for (pid,pers) in bur.pers
            if any_positive( pers.income, IncomesSet(things))
                return true
            end
        end
        return false
    end

    @with_kw mutable struct LocalTaxes{RT<:Real}
        # council_tax :: RT = zero(RT)
        # this isn't really used at present since LOCAL_TAXES
        # are in the incomes list, but we'll keep it around for modelling other types of local
        # tax in future.
        thing_just_to_make_with_kw_work :: RT = -99
    end

    function calc_net_income(incs::AbstractArray{T})::T where T
        # HB/CT treated seperately at household level
        # println( "incs="*inctostr( incs ))
        # println( ALL_INCOMES_EXCEPT_HOUSING_BENEFITS )
        # println( DIRECT_TAXES_AND_DEDUCTIONS )
        n :: T = isum(incs,
            ALL_INCOMES, # _EXCEPT_HOUSING_BENEFITS,
            deducted=DIRECT_TAXES_AND_DEDUCTIONS )
        # println( "n=$n")
        return n
    end

    @with_kw mutable struct HousingResult{RT<:Real}
        allowed_rooms :: Int = 0
        excess_rooms :: Int = 0
        allowed_rent :: RT = zero(RT) # FIXME this name
        gross_rent :: RT = zero(RT)
        # might not just be gross-allowed if we model other deductions
        rooms_rent_reduction :: RT = zero(RT) #
     end

    @with_kw mutable struct HouseholdResult{RT<:Real}
        income = STBIncomes.make_a( RT )

        bhc_net_income :: RT = zero(RT)
        eq_bhc_net_income :: RT = zero(RT)
        ahc_net_income :: RT = zero(RT)
        eq_ahc_net_income :: RT = zero(RT)

        net_housing_costs :: RT = zero(RT)
        housing = HousingResult{RT}()
        # FIXME note this is at the household level, which makes local income taxes, etc. akward. OK for now.
        local_tax = LocalTaxes{RT}()
        #
        # FIXME make this `people` to match the ModelHousehold
        # and adapt the get_benefit_units function
        #
        bus = Vector{BenefitUnitResult{RT}}(undef,0)
    end

    function get_net_income( hres :: HouseholdResult; target :: TargetBCIncomes ) :: Real
        if target == ahc_hh
            return hres.ahc_net_income
        elseif target == bhc_hh
            return hres.bhc_net_income
        elseif target == total_bens
            return isum(hres.income, BENEFITS )
        elseif target == total_taxes
            return isum(hres.income, INCOME_TAXES)
        end
    end

    function to_string( hr :: HouseholdResult, depth=0 )::String
        s = "# Results for Household"
        s *= to_md_table( hr, exclude=[:bus,:income], depth=depth )
        s *= "## HH Income"
        s *= inctostr( hr.income )
        s *= "\n\n"
        for bn in eachindex(hr.bus)
            s *= repeat("#",depth+1)*" Benefit Unit Result $(bn)"
            s *= to_string( hr.bus[bn], depth+1)
        end
        return s
    end

    """
    Test the individual incomes in hhr for some values
    """
    function has_any( hhr :: HouseholdResult, things ... ) :: Bool
        for bu in hhr.bus
            if has_any( bu, things ... )
                return true
            end
        end
        return false
    end

    function total( bur :: BenefitUnitResult{T}, which :: Incomes ) ::T where T
        t = zero(T)
        for (pid,pers) in bur.pers
            t += pers.income[which]
        end
        return t
    end

    function tozero!( bur :: BenefitUnitResult, which... )
        for (pid,pers) in bur.pers
            # println( pers.income )
            for w in which
                pers.income[w] = 0
            end
        end
    end

    function gross_pension_contributions( pres :: IndividualResult ) :: Real
        p = pres.income[PENSION_CONTRIBUTIONS_EMPLOYEE] +
            pres.it.pension_relief_at_source
        return max(0.0, p )
    end

    function gross_pension_contributions( bures :: BenefitUnitResult ) :: Real
        p = 0.0
        for pid in bures.adults
            p += gross_pension_contributions( bu.pers[pid])
        end
        return p;
    end

    function total( hhr :: HouseholdResult{T}, which :: Incomes ) :: T where T
        t = zero(T)
        for bu in hhr.bus
            t += total( bu, which )
        end
        return t
    end

    function add_to!( ni :: NIResult, ni2 :: NIResult )
        # ni.above_lower_earnings_limit += ni2.above_lower_earnings_limit
        # ni.total_ni += ni2.total_ni
        ni.class_1_primary    += ni2.class_1_primary
        ni.class_1_secondary  += ni2.class_1_secondary
        ni.class_2   += ni2.class_2
        ni.class_3   += ni2.class_3
        ni.class_4   += ni2.class_4
        ni.assumed_gross_wage += ni2.assumed_gross_wage
    end

    function add_to!( it :: ITResult, it2 :: ITResult )
        # it.total_tax += it2.total_tax
        it.taxable_income += it2.taxable_income
        it.adjusted_net_income += it2.adjusted_net_income
        it.total_income += it2.total_income
        it.allowance   += it2.allowance

        it.non_savings_tax += it2.non_savings_tax
        it.non_savings_band += it2.non_savings_band
        it.non_savings_income += it2.non_savings_income
        it.non_savings_taxable += it2.non_savings_taxable

        it.savings_tax += it2.savings_tax
        it.savings_band += it2.savings_band
        it.savings_income += it2.savings_income
        it.savings_taxable += it2.savings_taxable

        it.dividends_tax += it2.dividends_tax
        it.dividend_band += it2.dividend_band
        it.dividends_income += it2.dividends_income
        it.dividends_taxable += it2.dividends_taxable

        it.unused_allowance += it2.unused_allowance
        it.mca += it2.mca
        it.transferred_allowance += it2.transferred_allowance
        it.pension_eligible_for_relief += it2.pension_eligible_for_relief
        it.pension_relief_at_source += it2.pension_relief_at_source

        it.personal_savings_allowance += it2.personal_savings_allowance
    end

    """
    This is used to get aggregate taxable income for means-tested benefits;
    we don't really need bu/hh aggregated data in the results record other than total tax
    """
    function aggregate_tax( bu :: BenefitUnitResult{T}; include_children :: Bool = true ) :: Tuple where T
        pids = include_children ? keys( bu.pers ) : bu.adults
        it = ITResult{T}()
        ni = NIResult{T}()
        for pid in pids
            add_to!( it, bu.pers[pid].it )
            add_to!( ni, bu.pers[pid].ni )
        end
        return (it,ni)
    end

    function aggregate!( bures :: BenefitUnitResult{T} ) where T
        bures.income .= zero(T)
        for (pid,pers) in bures.pers
            bures.income .+= pers.income
        end
        bures.net_income = calc_net_income( bures.income )
    end

    function aggregate!( hh :: Household{T}, hres :: HouseholdResult{T} ) where T
        hres.income .= zero(T)
        for bu in hres.bus
            aggregate!( bu )
            hres.income .+= bu.income
        end
        #
        # FIXME it's possible we need BU net incomes as well
        #
        hres.bhc_net_income = calc_net_income( hres.income )
        # -
        #     hres.bus[1].uc.housing_element
        #    hres.income[HOUSING_BENEFIT] -
        #    hres.income[COUNCIL_TAX_BENEFIT]

        #
        # count pension relief at source as net income (but not as negative income tax)
        #
        for bu in hres.bus
            for (pid,pers) in bu.pers
                hres.bhc_net_income += pers.it.pension_relief_at_source
            end
        end
        if hres.bhc_net_income <= 0
            println("zero bhc_net_income for seq=$(hh.sequence); hid=$(hh.hid) year=$(hh.interview_year) ")
            println("income")
            println( inctostr( hres.income ))
            println( "bhc_net_income $(hres.bhc_net_income) HOUSING_BENEFIT=$(hres.income[HOUSING_BENEFIT]) COUNCIL_TAX_BENEFIT=$(hres.income[COUNCIL_TAX_BENEFIT]) hres.income[COUNCIL_TAX_BENEFIT]=$(hres.income[COUNCIL_TAX_BENEFIT])")
        end
        # FIXME this doesn't work with Universal Credit
        hres.net_housing_costs = hh.gross_rent +
            hres.income[LOCAL_TAXES] +
            hh.mortgage_payment +
            hh.other_housing_charges +
            hh.water_and_sewerage -
            hres.income[HOUSING_BENEFIT] -
            hres.income[COUNCIL_TAX_BENEFIT] - 0
            # hres.bus[1].uc.housing_element

        hres.ahc_net_income = hres.bhc_net_income -
            hh.gross_rent -
            hres.income[LOCAL_TAXES] -
            hh.mortgage_payment -
            hh.other_housing_charges -
            hh.water_and_sewerage

        # hres.net_housing_costs
        hres.eq_bhc_net_income = hres.bhc_net_income/hh.equivalence_scales.oecd_bhc
        hres.eq_ahc_net_income = hres.ahc_net_income/hh.equivalence_scales.oecd_ahc
    end

   function init_benefit_unit_result( T::Type, bu :: BenefitUnit ) :: BenefitUnitResult{T}
        # FIXME add an explicit head/spouse pid here
        bur = BenefitUnitResult{T}()
        bur.adults = bu.adults
        for pid in keys( bu.people )
            bur.pers[pid] = IndividualResult{T}()
            bur.pers[pid].income = map_incomes( bu.people[pid])

        end
        return bur
    end

    # create results that mirror some
    # allocation of people to benefit units
    # FIXME remove the type and use where RT
    function init_household_result( hh :: Household{T} ) :: HouseholdResult{T} where T
        bus = get_benefit_units(hh)
        hr = HouseholdResult{T}()
        for bu in bus
            push!( hr.bus, init_benefit_unit_result( T, bu ))
        end
        return hr
    end

    function get_indiv_result( hres :: HouseholdResult, pid :: BigInt  ) :: IndividualResult
        for bu in hres.bus
            if haskey( bu.pers, pid )
                return bu.pers[pid]
            end
        end
        @assert false "no pid $pid in household result"
        ## deliberate fail here
    end

    function map_incomes( pers :: Person{T}; include_calculated :: Bool=false ) :: MVector{INC_ARRAY_SIZE,T} where T
        out = MVector{INC_ARRAY_SIZE,T}( zeros(T,INC_ARRAY_SIZE ))
        incd = pers.income
        # println( "incd = $(incd)")
        if haskey(incd, Definitions.wages )
            out[WAGES] = incd[Definitions.wages]
        end
        if haskey(incd, Definitions.self_employment_income )
            out[SELF_EMPLOYMENT_INCOME] = incd[Definitions.self_employment_income]
        end
        if haskey(incd, Definitions.odd_jobs )
            out[ODD_JOBS] = incd[Definitions.odd_jobs]
        end
        if haskey(incd, Definitions.private_pensions )
            out[PRIVATE_PENSIONS] = incd[Definitions.private_pensions]
        end
        if haskey(incd, Definitions.national_savings )
            out[NATIONAL_SAVINGS] = incd[Definitions.national_savings]
        end
        if haskey(incd, Definitions.bank_interest )
            out[BANK_INTEREST] = incd[Definitions.bank_interest]
        end
        if haskey(incd, Definitions.stocks_shares )
            out[STOCKS_SHARES] = incd[Definitions.stocks_shares]
        end
        if haskey(incd, Definitions.individual_savings_account )
            out[INDIVIDUAL_SAVINGS_ACCOUNT] = incd[Definitions.individual_savings_account]
        end
        if haskey(incd, Definitions.property )
            out[PROPERTY] = incd[Definitions.property]
        end
        if haskey(incd, Definitions.royalties )
            out[ROYALTIES] = incd[Definitions.royalties]
        end
        if haskey(incd, Definitions.bonds_and_gilts )
            out[BONDS_AND_GILTS] = incd[Definitions.bonds_and_gilts]
        end
        if haskey(incd, Definitions.other_investment_income )
            out[OTHER_INVESTMENT_INCOME] = incd[Definitions.other_investment_income]
        end
        if haskey(incd, Definitions.other_income )
            out[OTHER_INCOME] = incd[Definitions.other_income]
        end
        if haskey(incd, Definitions.alimony_and_child_support_received )
            out[ALIMONY_AND_CHILD_SUPPORT_RECEIVED] = incd[Definitions.alimony_and_child_support_received]
        end
        if haskey(incd, Definitions.private_sickness_scheme_benefits )
            out[PRIVATE_SICKNESS_SCHEME_BENEFITS] = incd[Definitions.private_sickness_scheme_benefits]
        end
        if haskey(incd, Definitions.accident_insurance_scheme_benefits )
            out[ACCIDENT_INSURANCE_SCHEME_BENEFITS] = incd[Definitions.accident_insurance_scheme_benefits]
        end
        if haskey(incd, Definitions.hospital_savings_scheme_benefits )
            out[HOSPITAL_SAVINGS_SCHEME_BENEFITS] = incd[Definitions.hospital_savings_scheme_benefits]
        end
        if haskey(incd, Definitions.unemployment_or_redundancy_insurance )
            out[UNEMPLOYMENT_OR_REDUNDANCY_INSURANCE] = incd[Definitions.unemployment_or_redundancy_insurance]
        end
        if haskey(incd, Definitions.permanent_health_insurance )
            out[PERMANENT_HEALTH_INSURANCE] = incd[Definitions.permanent_health_insurance]
        end
        if haskey(incd, Definitions.any_other_sickness_insurance )
            out[ANY_OTHER_SICKNESS_INSURANCE] = incd[Definitions.any_other_sickness_insurance]
        end
        if haskey(incd, Definitions.critical_illness_cover )
            out[CRITICAL_ILLNESS_COVER] = incd[Definitions.critical_illness_cover]
        end
        if haskey(incd, Definitions.trade_union_sick_or_strike_pay )
            out[TRADE_UNION_SICK_OR_STRIKE_PAY] = incd[Definitions.trade_union_sick_or_strike_pay]
        end
        if haskey(incd, Definitions.health_insurance )
            out[HEALTH_INSURANCE] = incd[Definitions.health_insurance]
        end
        if haskey(incd, Definitions.alimony_and_child_support_paid )
            out[ALIMONY_AND_CHILD_SUPPORT_PAID] = incd[Definitions.alimony_and_child_support_paid]
        end
        if haskey(incd, Definitions.trade_unions_etc )
            out[TRADE_UNIONS_ETC] = incd[Definitions.trade_unions_etc]
        end
        if haskey(incd, Definitions.friendly_societies )
            out[FRIENDLY_SOCIETIES] = incd[Definitions.friendly_societies]
        end
        if haskey(incd, Definitions.work_expenses )
            out[WORK_EXPENSES] = incd[Definitions.work_expenses]
        end
        if haskey(incd, Definitions.avcs )
            out[AVCS] = incd[Definitions.avcs]
        end
        if haskey(incd, Definitions.other_deductions )
            out[OTHER_DEDUCTIONS] = incd[Definitions.other_deductions]
        end
        if haskey(incd, Definitions.loan_repayments )
            out[LOAN_REPAYMENTS] = incd[Definitions.loan_repayments]
        end
        if haskey(incd, Definitions.pension_contributions_employee )
            out[PENSION_CONTRIBUTIONS_EMPLOYEE] = incd[Definitions.pension_contributions_employee]
        end
        if haskey(incd, Definitions.pension_contributions_employer )
            out[PENSION_CONTRIBUTIONS_EMPLOYER] = incd[Definitions.pension_contributions_employer]
        end

        ### passed through benefits
        if haskey(incd, Definitions.other_benefits )
            out[OTHER_BENEFITS] = incd[Definitions.other_benefits]
        end
        if haskey(incd, Definitions.student_grants )
            out[STUDENT_GRANTS] = incd[Definitions.student_grants]
        end
        if haskey(incd, Definitions.student_loans )
            out[STUDENT_LOANS] = incd[Definitions.student_loans]
        end
        if haskey(incd, Definitions.free_school_meals )
            out[FREE_SCHOOL_MEALS] = incd[Definitions.free_school_meals]
        end
        if haskey(incd, Definitions.social_fund_loan_repayment_from_is_or_pc)
            out[SOCIAL_FUND_LOAN_REPAYMENT] = incd[Definitions.social_fund_loan_repayment_from_is_or_pc]
        end
        if haskey(incd, Definitions.social_fund_loan_repayment_from_jsa_or_esa)
            out[SOCIAL_FUND_LOAN_REPAYMENT] += incd[Definitions.social_fund_loan_repayment_from_jsa_or_esa]
        end
        if haskey(incd, Definitions.student_loan_repayments )
            out[STUDENT_LOAN_REPAYMENTS] = incd[Definitions.student_loan_repayments]
        end
        if haskey(incd, Definitions.armed_forces_compensation_scheme )
            out[ARMED_FORCES_COMPENSATION_SCHEME] = incd[Definitions.armed_forces_compensation_scheme]
        end
        if haskey(incd, Definitions.war_widows_or_widowers_pension )
            out[WAR_WIDOWS_PENSION] = incd[Definitions.war_widows_or_widowers_pension]
        end
        if haskey(incd, Definitions.severe_disability_allowance )
            out[SEVERE_DISABILITY_ALLOWANCE] = incd[Definitions.severe_disability_allowance]
        end
        if haskey(incd, Definitions.attendance_allowance )
            out[ATTENDANCE_ALLOWANCE] = incd[Definitions.attendance_allowance]
        end
        if haskey(incd, Definitions.foster_care_payments )
            out[FOSTER_CARE_PAYMENTS] = incd[Definitions.foster_care_payments]
        end
        if haskey(incd, Definitions.maternity_grant_from_social_fund )
            out[MATERNITY_GRANT] = incd[Definitions.maternity_grant_from_social_fund]
        end
        if haskey(incd, Definitions.funeral_grant_from_social_fund )
            out[FUNERAL_GRANT] = incd[Definitions.funeral_grant_from_social_fund]
        end
        if haskey(incd, Definitions.any_other_ni_or_state_benefit )
            out[ANY_OTHER_NI_OR_STATE_BENEFIT] = incd[Definitions.any_other_ni_or_state_benefit]
        end
        if haskey(incd, Definitions.friendly_society_benefits )
            out[FRIENDLY_SOCIETY_BENEFITS] = incd[Definitions.friendly_society_benefits]
        end
        if haskey(incd, Definitions.government_training_allowances )
            out[GOVERNMENT_TRAINING_ALLOWANCES] = incd[Definitions.government_training_allowances]
        end

        #
        #
        #

        if include_calculated
            if haskey(incd, Definitions.income_tax )
                out[INCOME_TAX] = incd[Definitions.income_tax]
            end
            if haskey(incd, Definitions.national_insurance )
                out[NATIONAL_INSURANCE] = incd[Definitions.national_insurance]
            end
            if haskey(incd, Definitions.local_taxes )
                out[LOCAL_TAXES] = incd[Definitions.local_taxes]
            end

            if haskey(incd, Definitions.care_insurance )
                out[CARE_INSURANCE] = incd[Definitions.care_insurance]
            end
            if haskey(incd, Definitions.child_benefit )
                out[CHILD_BENEFIT] = incd[Definitions.child_benefit]
            end
            if haskey(incd, Definitions.state_pension )
                out[STATE_PENSION] = incd[Definitions.state_pension]
            end
            if haskey(incd, Definitions.bereavement_allowance_or_widowed_parents_allowance_or_bereavement )
                out[BEREAVEMENT_ALLOWANCE] = incd[Definitions.bereavement_allowance_or_widowed_parents_allowance_or_bereavement]
            end
            if haskey(incd, Definitions.carers_allowance )
                out[CARERS_ALLOWANCE] = incd[Definitions.carers_allowance]
            end
            if haskey(incd, Definitions.industrial_injury_disablement_benefit )
                out[INDUSTRIAL_INJURY_BENEFIT] = incd[Definitions.industrial_injury_disablement_benefit]
            end
            if haskey(incd, Definitions.incapacity_benefit )
                out[INCAPACITY_BENEFIT] = incd[Definitions.incapacity_benefit]
            end
            if haskey(incd, Definitions.personal_independence_payment_daily_living )
                out[PERSONAL_INDEPENDENCE_PAYMENT_DAILY_LIVING] = incd[Definitions.personal_independence_payment_daily_living]
            end
            if haskey(incd, Definitions.personal_independence_payment_mobility )
                out[PERSONAL_INDEPENDENCE_PAYMENT_MOBILITY] = incd[Definitions.personal_independence_payment_mobility]
            end
            if haskey(incd, Definitions.dlaself_care )
                out[DLA_SELF_CARE] = incd[Definitions.dlaself_care]
            end
            if haskey(incd, Definitions.dlamobility )
                out[DLA_MOBILITY] = incd[Definitions.dlamobility]
            end
            if haskey(incd, Definitions.education_allowances )
                out[EDUCATION_ALLOWANCES] = incd[Definitions.education_allowances]
            end
            if haskey(incd, Definitions.maternity_allowance )
                out[MATERNITY_ALLOWANCE] = incd[Definitions.maternity_allowance]
            end
            if haskey(incd, Definitions.jobseekers_allowance )
                if pers.jsa_type == contributory_jsa
                    out[CONTRIB_JOBSEEKERS_ALLOWANCE] = incd[Definitions.jobseekers_allowance]
                elseif pers.jsa_type == income_related_jsa
                    out[NON_CONTRIB_JOBSEEKERS_ALLOWANCE] = incd[Definitions.jobseekers_allowance]
                elseif pers.jsa_type == both_jsa
                    out[NON_CONTRIB_JOBSEEKERS_ALLOWANCE] = incd[Definitions.jobseekers_allowance]/2
                    out[CONTRIB_JOBSEEKERS_ALLOWANCE] = incd[Definitions.jobseekers_allowance]/2
                else
                    @assert false "jsa is positive but jsa_type unset"
                end
            end
            if haskey(incd, Definitions.widows_payment )
                out[WIDOWS_PAYMENT] = incd[Definitions.widows_payment]
            end
            if haskey(incd, Definitions.winter_fuel_payments )
                out[WINTER_FUEL_PAYMENTS] = incd[Definitions.winter_fuel_payments]
            end
            if haskey(incd, Definitions.working_tax_credit )
                out[WORKING_TAX_CREDIT] = incd[Definitions.working_tax_credit]
            end
            if haskey(incd, Definitions.child_tax_credit )
                out[CHILD_TAX_CREDIT] = incd[Definitions.child_tax_credit]
            end
            if haskey(incd, Definitions.employment_and_support_allowance )
                if pers.esa_type == contributory_jsa
                    out[CONTRIB_EMPLOYMENT_AND_SUPPORT_ALLOWANCE] = incd[Definitions.employment_and_support_allowance]
                elseif pers.esa_type == income_related_jsa
                    out[NON_EMPLOYMENT_AND_SUPPORT_ALLOWANCE] = incd[Definitions.employment_and_support_allowance]
                elseif pers.esa_type == both_jsa
                    out[NON_CONTRIB_EMPLOYMENT_AND_SUPPORT_ALLOWANCE] = incd[Definitions.employment_and_support_allowance]/2
                    out[CONTRIB_EMPLOYMENT_AND_SUPPORT_ALLOWANCE] = incd[Definitions.employment_and_support_allowance]/2
                else
                    @assert false "jsa is positive but jsa_type unset"
                end
            end
            if haskey(incd, Definitions.income_support )
                out[INCOME_SUPPORT] = incd[Definitions.income_support]
            end
            if haskey(incd, Definitions.pension_credit )
                out[PENSION_CREDIT] = incd[Definitions.pension_credit]
            end
            # merged with pension credit in the frs, I think
            # if haskey(incd, Definitions.savings_credit )
            #    out[SAVINGS_CREDIT] = incd[Definitions.savings_credit]
            # end
            if haskey(incd, Definitions.housing_benefit )
                out[HOUSING_BENEFIT] = incd[Definitions.housing_benefit]
            end

            if haskey(incd, Definitions.extended_hb )
                out[HOUSING_BENEFIT] += incd[Definitions.extended_hb]
            end

            if haskey(incd, Definitions.working_tax_credit_lump_sum )
                out[WORKING_TAX_CREDIT] += incd[Definitions.working_tax_credit_lump_sum]
            end
            if haskey(incd, Definitions.child_tax_credit_lump_sum )
                out[CHILD_TAX_CREDIT] += incd[Definitions.child_tax_credit_lump_sum]
            end

            if haskey(incd, Definitions.universal_credit )
                out[UNIVERSAL_CREDIT] = incd[Definitions.universal_credit]
            end
            if haskey(incd, Definitions.guardians_allowance )
                out[GUARDIANS_ALLOWANCE] = incd[Definitions.guardians_allowance]
            end
            # not in the income list
            # if haskey(incd, Definitions.council_tax_rebate )
            #     out[COUNCIL_TAX_REBATE] = incd[Definitions.council_tax_rebate]
            # end
        end # include calculated
        return out
    end
end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
function get_flow_rate(obj::ThrustPropulsionModel, arg0::FieldSpacecraftState, arg1::Vector{RealFieldElement})
    return jcall(obj, "getFlowRate", RealFieldElement, (FieldSpacecraftState, Vector{RealFieldElement}), arg0, arg1)
end

function get_flow_rate(obj::ThrustPropulsionModel, arg0::SpacecraftState)
    return jcall(obj, "getFlowRate", jdouble, (SpacecraftState,), arg0)
end

function get_flow_rate(obj::ThrustPropulsionModel, arg0::SpacecraftState, arg1::Vector{jdouble})
    return jcall(obj, "getFlowRate", jdouble, (SpacecraftState, Vector{jdouble}), arg0, arg1)
end

function get_thrust_vector(obj::ThrustPropulsionModel, arg0::FieldSpacecraftState, arg1::Vector{RealFieldElement})
    return jcall(obj, "getThrustVector", FieldVector3D, (FieldSpacecraftState, Vector{RealFieldElement}), arg0, arg1)
end

function get_thrust_vector(obj::ThrustPropulsionModel, arg0::SpacecraftState)
    return jcall(obj, "getThrustVector", Vector3D, (SpacecraftState,), arg0)
end

function get_thrust_vector(obj::ThrustPropulsionModel, arg0::SpacecraftState, arg1::Vector{jdouble})
    return jcall(obj, "getThrustVector", Vector3D, (SpacecraftState, Vector{jdouble}), arg0, arg1)
end