text:
"""
Constants for 27pt 3D Helmholtz matrix-vector product
"""
struct helm3d_params{T}
    w3a::T
    wm2::T
    wm3::T
    wm4::T
    cNNN::T
    cx::T
    cy::T
    cz::T
    xy_coef::T
    xz_coef::T
    yz_coef::T
    px_lo::Array{Complex{T},1}
    px_hi::Array{Complex{T},1}
    px::Array{Complex{T},1}
    py_lo::Array{Complex{T},1}
    py_hi::Array{Complex{T},1}
    py::Array{Complex{T},1}
    pz_lo::Array{Complex{T},1}
    pz_hi::Array{Complex{T},1}
    pz::Array{Complex{T},1}
end

function helm3d_operto_mvp_forw_impl!(wn, x, y, n, params, deriv_mode)
    M, N, P = 1, 2, 3
    zero_x = complex(0.0,0.0)
    zero_w = zero(eltype(wn))
    nx,ny,nz = n
    is_deriv_mode = deriv_mode ? complex(0.0,0.0) : complex(1.0,0.0)
    num_blocks = Threads.nthreads()
    block_size = ceil(Int64,nz/num_blocks)
    Threads.@threads for i=1:num_blocks
        z_idx = (i-1)*block_size+1:min(nz,i*block_size)
        wn_window = zeros(eltype(wn),3,3,3)
        coef = zeros(Complex{Float64},3,3,3)
        x_window = zeros(Complex{Float64},3,3,3)

        @inbounds for k = z_idx
            cxz = (params.xz_coef*params.pz[k])
            cxzlo = (params.xz_coef*params.pz_lo[k])
            cxzhi = (params.xz_coef*params.pz_hi[k])
            cyz = (params.yz_coef*params.pz[k])
            cyzlo = (params.yz_coef*params.pz_lo[k])
            cyzhi = (params.yz_coef*params.pz_hi[k])

            for j = 1:ny
                c1 = params.xy_coef*params.py[j] + cxz
                c2_lo = params.cy*params.py_lo[j] + cyz
                c2_hi = params.cy*params.py_hi[j] + cyz
                c3_lo = params.yz_coef*params.py[j] + params.cz*params.pz_lo[k]
                c3_hi = params.yz_coef*params.py[j] + params.cz*params.pz_hi[k]

                dyzLL = - params.yz_coef*params.py_lo[j] - params.yz_coef*params.pz_lo[k]
                dyzLH = - params.yz_coef*params.py_lo[j] - params.yz_coef*params.pz_hi[k]
                dyzHL = - params.yz_coef*params.py_hi[j] - params.yz_coef*params.pz_lo[k]
                dyzHH = - params.yz_coef*params.py_hi[j] - params.yz_coef*params.pz_hi[k]

                eyzL = params.w3a*params.py[j] - params.xz_coef*params.pz_lo[k]
                eyzH = params.w3a*params.py[j] - params.xz_coef*params.pz_hi[k]
                fyzL = - params.xy_coef*params.py_lo[j] + params.w3a*params.pz[k]
                fyzH = - params.xy_coef*params.py_hi[j] + params.w3a*params.pz[k]
                gyzLL = params.py_lo[j] + params.pz_lo[k]
                gyzLH = params.py_lo[j] + params.pz_hi[k]
                gyzHL = params.py_hi[j] + params.pz_lo[k]
                gyzHH = params.py_hi[j] + params.pz_hi[k]

                for i = 1:nx
                    # Load wn_window + x_window
                    for kk=1:3
                        load_z = k+kk-2 > 0 && k+kk-2<=nz;
                        for jj=1:3
                            load_y = j+jj-2 > 0 && j+jj-2<=ny;
                            for ii=1:3
                                load_x = i+ii-2 > 0 && i+ii-2<=nx;
                                if load_x & load_y & load_z
                                    x_window[ii,jj,kk] = x[i+ii-2,j+jj-2,k+kk-2]
                                    wn_window[ii,jj,kk] = wn[i+ii-2,j+jj-2,k+kk-2]
                                else
                                    x_window[ii,jj,kk] = zero_x
                                    wn_window[ii,jj,kk] = zero_w
                                end
                            end
                        end
                    end

                    coef[M,N,N] = is_deriv_mode*(params.cx*params.px_lo[i] + c1) - params.wm2 * wn_window[M,N,N]
                    coef[P,N,N] = is_deriv_mode*(params.cx*params.px_hi[i] + c1) - params.wm2 * wn_window[P,N,N]
                    coef[N,M,N] = is_deriv_mode*(params.xy_coef*params.px[i] + c2_lo) - params.wm2 * wn_window[N,M,N]
                    coef[N,P,N] = is_deriv_mode*(params.xy_coef*params.px[i] + c2_hi) - params.wm2 * wn_window[N,P,N]
                    coef[N,N,M] = is_deriv_mode*(params.xz_coef*params.px[i] + c3_lo) - params.wm2 * wn_window[N,N,M]
                    coef[N,N,P] = is_deriv_mode*(params.xz_coef*params.px[i] + c3_hi) - params.wm2 * wn_window[N,N,P]

                    coef[N,M,M] = is_deriv_mode*(params.w3a*params.px[i] + dyzLL) - params.wm3 * wn_window[N,M,M]
                    coef[N,M,P] = is_deriv_mode*(params.w3a*params.px[i] + dyzLH) - params.wm3 * wn_window[N,M,P]
                    coef[N,P,M] = is_deriv_mode*(params.w3a*params.px[i] + dyzHL) - params.wm3 * wn_window[N,P,M]
                    coef[N,P,P] = is_deriv_mode*(params.w3a*params.px[i] + dyzHH) - params.wm3 * wn_window[N,P,P]
                    coef[M,N,M] = is_deriv_mode*(-params.xz_coef*params.px_lo[i] + eyzL) - params.wm3 * wn_window[M,N,M]
                    coef[M,N,P] = is_deriv_mode*(-params.xz_coef*params.px_lo[i] + eyzH) - params.wm3 * wn_window[M,N,P]
                    coef[P,N,M] = is_deriv_mode*(-params.xz_coef*params.px_hi[i] + eyzL) - params.wm3 * wn_window[P,N,M]
                    coef[P,N,P] = is_deriv_mode*(-params.xz_coef*params.px_hi[i] + eyzH) - params.wm3 * wn_window[P,N,P]
                    coef[M,M,N] = is_deriv_mode*(-params.xy_coef*params.px_lo[i] + fyzL) - params.wm3 * wn_window[M,M,N]
                    coef[M,P,N] = is_deriv_mode*(-params.xy_coef*params.px_lo[i] + fyzH) - params.wm3 * wn_window[M,P,N]
                    coef[P,M,N] = is_deriv_mode*(-params.xy_coef*params.px_hi[i] + fyzL) - params.wm3 * wn_window[P,M,N]
                    coef[P,P,N] = is_deriv_mode*(-params.xy_coef*params.px_hi[i] + fyzH) - params.wm3 * wn_window[P,P,N]
                    coef[M,M,M] = is_deriv_mode*(-params.w3a*(params.px_lo[i] + gyzLL)) - params.wm4 * wn_window[M,M,M]
                    coef[M,M,P] = is_deriv_mode*(-params.w3a*(params.px_lo[i] + gyzLH)) - params.wm4 * wn_window[M,M,P]
                    coef[M,P,M] = is_deriv_mode*(-params.w3a*(params.px_lo[i] + gyzHL)) - params.wm4 * wn_window[M,P,M]
                    coef[M,P,P] = is_deriv_mode*(-params.w3a*(params.px_lo[i] + gyzHH)) - params.wm4 * wn_window[M,P,P]
                    coef[P,M,M] = is_deriv_mode*(-params.w3a*(params.px_hi[i] + gyzLL)) - params.wm4 * wn_window[P,M,M]
                    coef[P,M,P] = is_deriv_mode*(-params.w3a*(params.px_hi[i] + gyzLH)) - params.wm4 * wn_window[P,M,P]
                    coef[P,P,M] = is_deriv_mode*(-params.w3a*(params.px_hi[i] + gyzHL)) - params.wm4 * wn_window[P,P,M]
                    coef[P,P,P] = is_deriv_mode*(-params.w3a*(params.px_hi[i] + gyzHH)) - params.wm4 * wn_window[P,P,P]
                    coef[N,N,N] = is_deriv_mode*(-params.cx*params.px[i] - params.cy*params.py[j] - params.cz*params.pz[k]) + params.cNNN*wn_window[N,N,N]

                    t = complex(0.0,0.0);
                    for kk=1:3
                        for jj=1:3
                            for ii=1:3
                                t += coef[ii,jj,kk]*x_window[ii,jj,kk];
                            end
                        end
                    end
                    y[i,j,k] = t;
                end
            end
        end
    end
end

function helm3d_operto_mvp_adj_impl!(wn, x, y, n, params, deriv_mode)
    M, N, P = 1, 2, 3
    zero_x = complex(0.0,0.0)
    zero_w = zero(eltype(wn))
    nx,ny,nz = n
    is_deriv_mode = deriv_mode ? complex(0.0,0.0) : complex(1.0,0.0)
    num_blocks = Threads.nthreads()
    block_size = ceil(Int64,nz/num_blocks)
    Threads.@threads for i=1:num_blocks
        z_idx = (i-1)*block_size+1:min(nz,i*block_size)
        coef = zeros(Complex{Float64},3,3,3)
        x_window = zeros(Complex{Float64},3,3,3)

        # (i+1,j+1,k+1) is the current index for the current point
        # with respect to the pml functions
        @inbounds for k = z_idx
            kM,kN,kP = k,k+1,k+2
            for j=1:ny
                jM,jN,jP = j,j+1,j+2
                for i=1:nx
                    wNNN = wn[i,j,k]
                    for kk=1:3
                        load_z = k+kk-2 > 0 && k+kk-2<=nz;
                        for jj=1:3
                            load_y = j+jj-2 > 0 && j+jj-2<=ny;
                            for ii=1:3
                                load_x = i+ii-2 > 0 && i+ii-2<=nx;
                                if load_x & load_y & load_z
                                    x_window[ii,jj,kk] = x[i+ii-2,j+jj-2,k+kk-2]
                                else
                                    x_window[ii,jj,kk] = zero_x
                                end
                            end
                        end
                    end

                    iM,iN,iP = i,i+1,i+2
                    coef[M,M,M] = -params.wm4*wNNN + is_deriv_mode*(-params.w3a*(params.px_hi[iM] + params.py_hi[jM] + params.pz_hi[kM]))
                    coef[N,M,M] = -params.wm3*wNNN + is_deriv_mode*((-params.yz_coef) * (params.pz_hi[kM] + params.py_hi[jM]) + params.w3a*params.px[iN])
                    coef[P,M,M] = -params.wm4*wNNN + is_deriv_mode*(-params.w3a*( params.px_lo[iP] + params.py_hi[jM] + params.pz_hi[kM]))
                    coef[M,N,M] = -params.wm3*wNNN + is_deriv_mode*(-params.xz_coef * (params.px_hi[iM] + params.pz_hi[kM]) + params.w3a*params.py[jN])
                    coef[N,N,M] = -params.wm2*wNNN + is_deriv_mode*(params.cz*params.pz_hi[kM] + params.yz_coef*params.py[jN] + params.xz_coef*params.px[iN])
                    coef[P,N,M] = -params.wm3*wNNN + is_deriv_mode*(-params.xz_coef * (params.pz_hi[kM] + params.px_lo[iP]) + params.w3a*params.py[jN])
                    coef[M,P,M] = -params.wm4*wNNN + is_deriv_mode*(-params.w3a*( params.px_hi[iM] + params.py_lo[jP] + params.pz_hi[kM] ));
                    coef[N,P,M] = -params.wm3*wNNN + is_deriv_mode*(-params.yz_coef * (params.py_lo[jP] + params.pz_hi[kM]) + params.w3a*params.px[iN])
                    coef[P,P,M] = -params.wm4*wNNN + is_deriv_mode*(-params.w3a*( params.px_lo[iP] + params.py_lo[jP] + params.pz_hi[kM]))
                    coef[M,M,N] = -params.wm3*wNNN + is_deriv_mode*(-params.xy_coef * (params.px_hi[iM] + params.py_hi[jM]) + params.w3a*params.pz[kN])
                    coef[N,M,N] = -params.wm2*wNNN + is_deriv_mode*(params.cy*params.py_hi[jM] + params.yz_coef*params.pz[kN] + params.xy_coef*params.px[iN])
                    coef[P,M,N] = -params.wm3*wNNN + is_deriv_mode*(-params.xy_coef * (params.px_lo[iP] + params.py_hi[jM]) + params.w3a*params.pz[kN]	)
                    coef[M,N,N] = -params.wm2*wNNN + is_deriv_mode*(params.cx*params.px_hi[iM] + params.xz_coef*params.pz[kN] + params.xy_coef*params.py[jN])
                    coef[N,N,N] = params.cNNN*wNNN - is_deriv_mode*( params.cx*params.px[iN] + params.cy*params.py[jN] + params.cz*params.pz[kN])
                    coef[P,N,N] = -params.wm2*wNNN + is_deriv_mode*(params.cx*params.px_lo[iP] + params.xz_coef*params.pz[kN] + params.xy_coef*params.py[jN])
                    coef[M,P,N] = -params.wm3*wNNN + is_deriv_mode*(-params.xy_coef * (params.px_hi[iM] + params.py_lo[jP]) + params.w3a*params.pz[kN])
                    coef[N,P,N] = -params.wm2*wNNN + is_deriv_mode*(params.cy*params.py_lo[jP] + params.yz_coef*params.pz[kN] + params.xy_coef*params.px[iN])
                    coef[P,P,N] = -params.wm3*wNNN + is_deriv_mode*(-params.xy_coef * (params.px_lo[iP] + params.py_lo[jP]) + params.w3a*params.pz[kN])
                    coef[M,M,P] = -params.wm4*wNNN + is_deriv_mode*(-params.w3a*( params.px_hi[iM] + params.py_hi[jM] + params.pz_lo[kP] ))
                    coef[N,M,P] = -params.wm3*wNNN + is_deriv_mode*(-params.yz_coef * (params.py_hi[jM] + params.pz_lo[kP]) + params.w3a*params.px[iN])
                    coef[P,M,P] = -params.wm4*wNNN + is_deriv_mode*(-params.w3a*( params.px_lo[iP] + params.py_hi[jM] + params.pz_lo[kP]))
                    coef[M,N,P] = -params.wm3*wNNN + is_deriv_mode*(-params.xz_coef * (params.pz_lo[kP] + params.px_hi[iM]) + params.w3a*params.py[jN])
                    coef[N,N,P] = -params.wm2*wNNN + is_deriv_mode*(params.cz*params.pz_lo[kP] + params.yz_coef*params.py[jN] + params.xz_coef*params.px[iN])
                    coef[P,N,P] = -params.wm3*wNNN + is_deriv_mode*(-params.xz_coef * (params.pz_lo[kP] + params.px_lo[iP]) + params.w3a*params.py[jN])
                    coef[M,P,P] = -params.wm4*wNNN + is_deriv_mode*(-params.w3a*( params.px_hi[iM] + params.py_lo[jP] + params.pz_lo[kP]))
                    coef[N,P,P] = -params.wm3*wNNN + is_deriv_mode*(-params.yz_coef * (params.pz_lo[kP] + params.py_lo[jP]) + params.w3a*params.px[iN])
                    coef[P,P,P] = -params.wm4*wNNN + is_deriv_mode*(-params.w3a*( params.px_lo[iP] + params.py_lo[jP] + params.pz_lo[kP]))

                    t = complex(0.0,0.0);
                    for kk=1:3
                        for jj=1:3
                            for ii=1:3
                                t += conj(coef[ii,jj,kk])*x_window[ii,jj,kk];
                            end
                        end
                    end
                    y[i,j,k] = t;
                end
            end
        end
    end
end

function helm3d_operto_mvp_mt(wn,Δ,n,freq,npml,x; forw_mode::Bool=true, deriv_mode::Bool=false)
    M,N,P = 1,2,3

    nx,ny,nz = n[1],n[2],n[3]
    hx,hy,hz = Δ[1]^2,Δ[2]^2,Δ[3]^2
    hyz  = hy + hz
    hxy = hy + hx
    hxz = hx + hz
    hxyz = hx + hy + hz

    npxlo,npxhi,npylo,npyhi,npzlo,npzhi = npml[1,1],npml[2,1],npml[1,2],npml[2,2],npml[1,3],npml[2,3]

    y = zeros(Complex{Float64},nx,ny,nz)

    # Create PML functions
    px_lo,px_hi = pml_func(nx,npxlo,npxhi)
    if ~forw_mode
        prepend!(px_lo,[px_lo[1]])
        append!(px_lo,[px_lo[end]])
        prepend!(px_hi,[px_hi[1]])
        append!(px_hi,[px_hi[end]])
    end
    px = px_lo+px_hi
    py_lo,py_hi = pml_func(ny,npylo,npyhi)
    if ~forw_mode
        prepend!(py_lo,[py_lo[1]])
        append!(py_lo,[py_lo[end]])
        prepend!(py_hi,[py_hi[1]])
        append!(py_hi,[py_hi[end]])
    end

    py = py_lo+py_hi
    pz_lo,pz_hi = pml_func(nz,npzlo,npzhi)
    if ~forw_mode
        prepend!(pz_lo,[pz_lo[1]])
        append!(pz_lo,[pz_lo[end]])
        prepend!(pz_hi,[pz_hi[1]])
        append!(pz_hi,[pz_hi[end]])
    end

    pz = pz_lo+pz_hi

    # Weights
    w1  = 1.8395262e-5
    w2  = 0.29669233333333334
    w3  = 0.02747615
    wm1 = 0.49649658
    wm2 = 0.07516874999999999
    wm3 = 0.004373916666666667
    wm4 = 5.690375e-7
    w3a = 2*(w3*3/(4*hxyz))

    cx = - (w1/hx + w2/hx + w2/hxz + w2/hxy + 4*w3a)
    cy = - (w1/hy + w2/hy + w2/hyz + w2/hxy + 4*w3a)
    cz = - (w1/hz + w2/hz + w2/hxz + w2/hyz + 4*w3a)
    cNNN = - (w1 + 3*w2 + (8*w3a*hxyz)/3 + wm1-1)

    xz_coef = w2/(2*hxz)
    xy_coef = w2/(2*hxy)
    yz_coef = w2/(2*hyz)
    params = helm3d_params{Float64}(w3a,wm2,wm3,wm4,cNNN,cx,cy,cz,xy_coef,xz_coef,yz_coef,px_lo,px_hi,px,py_lo,py_hi,py,pz_lo,pz_hi,pz)
    if forw_mode
        helm3d_operto_mvp_forw_impl!(wn, x, y, n, params, deriv_mode)
    else
        helm3d_operto_mvp_adj_impl!(wn, x, y, n, params, deriv_mode)
    end
    return vec(y)
end

function helm3d_operto_mvp(wn::Union{AbstractArray{F,3},AbstractArray{Complex{F},3}},Δ::AbstractArray{F,1},n::AbstractArray{I,1},freq::Union{F,Complex{F}},npml::AbstractArray{I,2},x::AbstractArray{Complex{F},3};forw_mode::Bool=true,deriv_mode::Bool=false) where {F<:Real,I<:Integer}

    M,N,P = 1,2,3

    nx,ny,nz = n[1],n[2],n[3]
    hx,hy,hz = Δ[1]^2,Δ[2]^2,Δ[3]^2
    hyz  = hy + hz
    hxy = hy + hx
    hxz = hx + hz
    hxyz = hx + hy + hz

    npxlo,npxhi,npylo,npyhi,npzlo,npzhi = npml[1,1],npml[2,1],npml[1,2],npml[2,2],npml[1,3],npml[2,3]

    y = zeros(Complex{F},tuple(n...))

    # Create PML functions
    px_lo,px_hi = pml_func(nx,npxlo,npxhi)
    if ~forw_mode
        prepend!(px_lo,[px_lo[1]])
        append!(px_lo,[px_lo[end]])
        prepend!(px_hi,[px_hi[1]])
        append!(px_hi,[px_hi[end]])
    end
    px = px_lo+px_hi
    py_lo,py_hi = pml_func(ny,npylo,npyhi)
    if ~forw_mode
        prepend!(py_lo,[py_lo[1]])
        append!(py_lo,[py_lo[end]])
        prepend!(py_hi,[py_hi[1]])
        append!(py_hi,[py_hi[end]])
    end

    py = py_lo+py_hi
    pz_lo,pz_hi = pml_func(nz,npzlo,npzhi)
    if ~forw_mode
        prepend!(pz_lo,[pz_lo[1]])
        append!(pz_lo,[pz_lo[end]])
        prepend!(pz_hi,[pz_hi[1]])
        append!(pz_hi,[pz_hi[end]])
    end

    pz = pz_lo+pz_hi

    # Weights
    w1  = 1.8395262e-5
    w2  = 0.29669233333333334
    w3  = 0.02747615
    wm1 = 0.49649658
    wm2 = 0.07516874999999999
    wm3 = 0.004373916666666667
    wm4 = 5.690375e-7
    w3a = 2*(w3*3/(4*hxyz))::F

    cx = - (w1/hx + w2/hx + w2/hxz + w2/hxy + 4*w3a)
    cy = - (w1/hy + w2/hy + w2/hyz + w2/hxy + 4*w3a)
    cz = - (w1/hz + w2/hz + w2/hxz + w2/hyz + 4*w3a)
    cNNN = - (w1 + 3*w2 + (8*w3a*hxyz)/3 + wm1-1)

    xz_coef = w2/(2*hxz)
    xy_coef = w2/(2*hxy)
    yz_coef = w2/(2*hyz)
    wn_window = zeros(F,3,3,3)
    coef = zeros(Complex{F},3,3,3)
    x_window = zeros(eltype(x),3,3,3)
    zero_x = complex(0.0,0.0)
    zero_w = zero(eltype(wn))
    is_deriv_mode = deriv_mode ? complex(0.0,0.0) : complex(1.0,0.0)
    if forw_mode
        for k = 1:nz
            cxz = xz_coef*pz[k]
	        cxzlo = xz_coef*pz_lo[k]
	        cxzhi = xz_coef*pz_hi[k]
	        cyz = yz_coef*pz[k]
	        cyzlo = yz_coef*pz_lo[k]
	        cyzhi = yz_coef*pz_hi[k]
            for j = 1:ny
                c1 = xy_coef*py[j] + cxz
		        c2_lo = cy*py_lo[j] + cyz
		        c2_hi = cy*py_hi[j] + cyz
		        c3_lo = yz_coef*py[j] + cz*pz_lo[k]
		        c3_hi = yz_coef*py[j] + cz*pz_hi[k]

                dyzLL = - yz_coef*py_lo[j] - yz_coef*pz_lo[k]
                dyzLH = - yz_coef*py_lo[j] - yz_coef*pz_hi[k]
                dyzHL = - yz_coef*py_hi[j] - yz_coef*pz_lo[k]
                dyzHH = - yz_coef*py_hi[j] - yz_coef*pz_hi[k]

                eyzL = w3a*py[j] - xz_coef*pz_lo[k]
                eyzH = w3a*py[j] - xz_coef*pz_hi[k]
                fyzL = - xy_coef*py_lo[j] + w3a*pz[k]
                fyzH = - xy_coef*py_hi[j] + w3a*pz[k]
                gyzLL = py_lo[j] + pz_lo[k]
                gyzLH = py_lo[j] + pz_hi[k]
                gyzHL = py_hi[j] + pz_lo[k]
                gyzHH = py_hi[j] + pz_hi[k]

                for i = 1:nx
                    # Load wn_window + x_window
                    for kk=1:3
                        load_z = k+kk-2 > 0 && k+kk-2<=nz;
                        for jj=1:3
                            load_y = j+jj-2 > 0 && j+jj-2<=ny;
                            for ii=1:3
                                load_x = i+ii-2 > 0 && i+ii-2<=nx;
                                if load_x & load_y & load_z
                                    @inbounds x_window[ii,jj,kk] = x[i+ii-2,j+jj-2,k+kk-2]
                                    @inbounds wn_window[ii,jj,kk] = wn[i+ii-2,j+jj-2,k+kk-2]
                                else
                                    @inbounds x_window[ii,jj,kk] = zero_x
                                    @inbounds wn_window[ii,jj,kk] = zero_w
                                end
                            end
                        end
                    end

			        coef[M,N,N] = is_deriv_mode*(cx*px_lo[i] + c1) - wm2 * wn_window[M,N,N]
			        coef[P,N,N] = is_deriv_mode*(cx*px_hi[i] + c1) - wm2 * wn_window[P,N,N]
			        coef[N,M,N] = is_deriv_mode*(xy_coef*px[i] + c2_lo) - wm2 * wn_window[N,M,N]
			        coef[N,P,N] = is_deriv_mode*(xy_coef*px[i] + c2_hi) - wm2 * wn_window[N,P,N]
			        coef[N,N,M] = is_deriv_mode*(xz_coef*px[i] + c3_lo) - wm2 * wn_window[N,N,M]
			        coef[N,N,P] = is_deriv_mode*(xz_coef*px[i] + c3_hi) - wm2 * wn_window[N,N,P]

			        coef[N,M,M] = is_deriv_mode*(w3a*px[i] + dyzLL) - wm3 * wn_window[N,M,M]
			        coef[N,M,P] = is_deriv_mode*(w3a*px[i] + dyzLH) - wm3 * wn_window[N,M,P]
			        coef[N,P,M] = is_deriv_mode*(w3a*px[i] + dyzHL) - wm3 * wn_window[N,P,M]
			        coef[N,P,P] = is_deriv_mode*(w3a*px[i] + dyzHH) - wm3 * wn_window[N,P,P]
			        coef[M,N,M] = is_deriv_mode*(-xz_coef*px_lo[i] + eyzL) - wm3 * wn_window[M,N,M]
			        coef[M,N,P] = is_deriv_mode*(-xz_coef*px_lo[i] + eyzH) - wm3 * wn_window[M,N,P]
			        coef[P,N,M] = is_deriv_mode*(-xz_coef*px_hi[i] + eyzL) - wm3 * wn_window[P,N,M]
			        coef[P,N,P] = is_deriv_mode*(-xz_coef*px_hi[i] + eyzH) - wm3 * wn_window[P,N,P]
			        coef[M,M,N] = is_deriv_mode*(-xy_coef*px_lo[i] + fyzL) - wm3 * wn_window[M,M,N]
			        coef[M,P,N] = is_deriv_mode*(-xy_coef*px_lo[i] + fyzH) - wm3 * wn_window[M,P,N]
			        coef[P,M,N] = is_deriv_mode*(-xy_coef*px_hi[i] + fyzL) - wm3 * wn_window[P,M,N]
			        coef[P,P,N] = is_deriv_mode*(-xy_coef*px_hi[i] + fyzH) - wm3 * wn_window[P,P,N]
			        coef[M,M,M] = is_deriv_mode*(-w3a*(px_lo[i] + gyzLL)) - wm4 * wn_window[M,M,M]
			        coef[M,M,P] = is_deriv_mode*(-w3a*(px_lo[i] + gyzLH)) - wm4 * wn_window[M,M,P]
			        coef[M,P,M] = is_deriv_mode*(-w3a*(px_lo[i] + gyzHL)) - wm4 * wn_window[M,P,M]
			        coef[M,P,P] = is_deriv_mode*(-w3a*(px_lo[i] + gyzHH)) - wm4 * wn_window[M,P,P]
			        coef[P,M,M] = is_deriv_mode*(-w3a*(px_hi[i] + gyzLL)) - wm4 * wn_window[P,M,M]
			        coef[P,M,P] = is_deriv_mode*(-w3a*(px_hi[i] + gyzLH)) - wm4 * wn_window[P,M,P]
			        coef[P,P,M] = is_deriv_mode*(-w3a*(px_hi[i] + gyzHL)) - wm4 * wn_window[P,P,M]
			        coef[P,P,P] = is_deriv_mode*(-w3a*(px_hi[i] + gyzHH)) - wm4 * wn_window[P,P,P]
                    coef[N,N,N] = is_deriv_mode*(-cx*px[i] - cy*py[j] - cz*pz[k]) + cNNN*wn_window[N,N,N]
                    t = complex(0.0,0.0);
                    for kk=1:3
                        for jj=1:3
                            for ii=1:3
                                @inbounds t += coef[ii,jj,kk]*x_window[ii,jj,kk];
                            end
                        end
                    end
                    y[i,j,k] = t;
                end
            end
        end
    else
    # (i+1,j+1,k+1) is the current index for the current point
    # with respect to the pml functions
        for k=1:nz
            kM,kN,kP = k,k+1,k+2
            for j=1:ny
                jM,jN,jP = j,j+1,j+2
                for i=1:nx
                    wNNN = wn[i,j,k]
                    for kk=1:3
                        load_z = k+kk-2 > 0 && k+kk-2<=nz;
                        for jj=1:3
                            load_y = j+jj-2 > 0 && j+jj-2<=ny;
                            for ii=1:3
                                load_x = i+ii-2 > 0 && i+ii-2<=nx;
                                if load_x & load_y & load_z
                                @inbounds x_window[ii,jj,kk] = x[i+ii-2,j+jj-2,k+kk-2]
                                else
                                 @inbounds x_window[ii,jj,kk] = zero_x
                                end
                            end
                        end
                    end
                    iM,iN,iP = i,i+1,i+2
                    coef[M,M,M] = -wm4*wNNN + is_deriv_mode*(-w3a*(px_hi[iM] + py_hi[jM] + pz_hi[kM]))
                    coef[N,M,M] = -wm3*wNNN + is_deriv_mode*((-yz_coef) * (pz_hi[kM] + py_hi[jM]) + w3a*px[iN])
                    coef[P,M,M] = -wm4*wNNN + is_deriv_mode*(-w3a*( px_lo[iP] + py_hi[jM] + pz_hi[kM]))
                    coef[M,N,M] = -wm3*wNNN + is_deriv_mode*(-xz_coef * (px_hi[iM] + pz_hi[kM]) + w3a*py[jN])

                    coef[N,N,M] = -wm2*wNNN + is_deriv_mode*(cz*pz_hi[kM] + yz_coef*py[jN] + xz_coef*px[iN])
                    coef[P,N,M] = -wm3*wNNN + is_deriv_mode*(-xz_coef * (pz_hi[kM] + px_lo[iP]) + w3a*py[jN])
                    coef[M,P,M] = -wm4*wNNN + is_deriv_mode*(-w3a*( px_hi[iM] + py_lo[jP] + pz_hi[kM] ));
                    coef[N,P,M] = -wm3*wNNN + is_deriv_mode*(-yz_coef * (py_lo[jP] + pz_hi[kM]) + w3a*px[iN])
                    coef[P,P,M] = -wm4*wNNN + is_deriv_mode*(-w3a*( px_lo[iP] + py_lo[jP] + pz_hi[kM]))
                    coef[M,M,N] = -wm3*wNNN + is_deriv_mode*(-xy_coef * (px_hi[iM] + py_hi[jM]) + w3a*pz[kN])
                    coef[N,M,N] = -wm2*wNNN + is_deriv_mode*(cy*py_hi[jM] + yz_coef*pz[kN] + xy_coef*px[iN])
                    coef[P,M,N] = -wm3*wNNN + is_deriv_mode*(-xy_coef * (px_lo[iP] + py_hi[jM]) + w3a*pz[kN]	)
                    coef[M,N,N] = -wm2*wNNN + is_deriv_mode*(cx*px_hi[iM] + xz_coef*pz[kN] + xy_coef*py[jN])
                    coef[N,N,N] = cNNN*wNNN - is_deriv_mode*( cx*px[iN] + cy*py[jN] + cz*pz[kN])
                    coef[P,N,N] = -wm2*wNNN + is_deriv_mode*(cx*px_lo[iP] + xz_coef*pz[kN] + xy_coef*py[jN])

                    coef[M,P,N] = -wm3*wNNN + is_deriv_mode*(-xy_coef * (px_hi[iM] + py_lo[jP]) + w3a*pz[kN])
                    coef[N,P,N] = -wm2*wNNN + is_deriv_mode*(cy*py_lo[jP] + yz_coef*pz[kN] + xy_coef*px[iN])
                    coef[P,P,N] = -wm3*wNNN + is_deriv_mode*(-xy_coef * (px_lo[iP] + py_lo[jP]) + w3a*pz[kN])
                    coef[M,M,P] = -wm4*wNNN + is_deriv_mode*(-w3a*( px_hi[iM] + py_hi[jM] + pz_lo[kP] ))
                    coef[N,M,P] = -wm3*wNNN + is_deriv_mode*(-yz_coef * (py_hi[jM] + pz_lo[kP]) + w3a*px[iN])
                    coef[P,M,P] = -wm4*wNNN + is_deriv_mode*(-w3a*( px_lo[iP] + py_hi[jM] + pz_lo[kP]))
                    coef[M,N,P] = -wm3*wNNN + is_deriv_mode*(-xz_coef * (pz_lo[kP] + px_hi[iM]) + w3a*py[jN])
                    coef[N,N,P] = -wm2*wNNN + is_deriv_mode*(cz*pz_lo[kP] + yz_coef*py[jN] + xz_coef*px[iN])
                    coef[P,N,P] = -wm3*wNNN + is_deriv_mode*(-xz_coef * (pz_lo[kP] + px_lo[iP]) + w3a*py[jN])
                    coef[M,P,P] = -wm4*wNNN + is_deriv_mode*(-w3a*( px_hi[iM] + py_lo[jP] + pz_lo[kP]))
                    coef[N,P,P] = -wm3*wNNN + is_deriv_mode*(-yz_coef * (pz_lo[kP] + py_lo[jP]) + w3a*px[iN])
                    coef[P,P,P] = -wm4*wNNN + is_deriv_mode*(-w3a*( px_lo[iP] + py_lo[jP] + pz_lo[kP]))

                    t = complex(0.0,0.0);
                    for kk=1:3
                        for jj=1:3
                            for ii=1:3
                                @inbounds t += conj(coef[ii,jj,kk])*x_window[ii,jj,kk];
                            end
                        end
                    end
                    y[i,j,k] = t;
                end
            end
        end
    end
    return vec(y)
end

function helm3d_operto_matrix(wn::Union{AbstractArray{F,3},AbstractArray{Complex{F},3}},Δ::AbstractArray{F,1},n::AbstractArray{I,1},freq::Union{F,Complex{F}},npml::AbstractArray{I,2}) where {F<:Real,I<:Integer}

    M,N,P = 1,2,3

    nx,ny,nz = n[1],n[2],n[3]
    hx,hy,hz = Δ[1]^2,Δ[2]^2,Δ[3]^2
    hyz  = hy + hz
    hxy = hy + hx
    hxz = hx + hz
    hxyz = hx + hy + hz

    npxlo,npxhi,npylo,npyhi,npzlo,npzhi = npml[1,1],npml[2,1],npml[1,2],npml[2,2],npml[1,3],npml[2,3]

    y = zeros(Complex{F},tuple(n...))

    # Create PML functions
    px_lo,px_hi = pml_func(nx,npxlo,npxhi)
    px = px_lo+px_hi
    py_lo,py_hi = pml_func(ny,npylo,npyhi)
    py = py_lo+py_hi
    pz_lo,pz_hi = pml_func(nz,npzlo,npzhi)
    pz = pz_lo+pz_hi

    # Weights
    w1  = 1.8395262e-5
    w2  = 0.29669233333333334
    w3  = 0.02747615
    wm1 = 0.49649658
    wm2 = 0.07516874999999999
    wm3 = 0.004373916666666667
    wm4 = 5.690375e-7
    w3a = (0.041214225/hxyz)::F

    cx = - (w1/hx + w2/hx + w2/hxz + w2/hxy + 8*w3a)
    cy = - (w1/hy + w2/hy + w2/hyz + w2/hxy + 8*w3a)
    cz = - (w1/hz + w2/hz + w2/hxz + w2/hyz + 8*w3a)
    cNNN = - (w1 + 3*w2 + (16*w3a*hxyz)/3 + wm1-1)

    xz_coef = w2/(2*hxz)
    xy_coef = w2/(2*hxy)
    yz_coef = w2/(2*hyz)
    wn_window = zeros(F,3,3,3)
    coef = zeros(Complex{F},3,3,3,nx,ny,nz)

    for k = 1:nz
        ks = (k > 1) ? -1 : 0
        ke = (k < nz) ? 1 : 0
        zoff = ks:ke
        for j = 1:ny
            js = (j > 1) ? -1 : 0
            je = (j < ny) ? 1 : 0
            yoff = js:je
            for i = 1:nx
                is = (i > 1) ? -1 : 0
                ie = (i < nx) ? 1 : 0
                xoff = is:ie
                # Load wn_window
			    wn_window[2+xoff,2+yoff,2+zoff] = wn[i+xoff,j+yoff,k+zoff]
                cxz = xz_coef*pz[k]
	            cxzlo = xz_coef*pz_lo[k]
	            cxzhi = xz_coef*pz_hi[k]
	            cyz = yz_coef*pz[k]
	            cyzlo = yz_coef*pz_lo[k]
	            cyzhi = yz_coef*pz_hi[k]

                c1 = xy_coef*py[j] + cxz
		        c2_lo = cy*py_lo[j] + cyz
		        c2_hi = cy*py_hi[j] + cyz
		        c3_lo = yz_coef*py[j] + cz*pz_lo[k]
		        c3_hi = yz_coef*py[j] + cz*pz_hi[k]

			    coef[M,N,N,i,j,k] = cx*px_lo[i] + c1 - wm2 * wn_window[M,N,N]
			    coef[P,N,N,i,j,k] = cx*px_hi[i] + c1 - wm2 * wn_window[P,N,N]
			    coef[N,M,N,i,j,k] = xy_coef*px[i] + c2_lo - wm2 * wn_window[N,M,N]
			    coef[N,P,N,i,j,k] = xy_coef*px[i] + c2_hi - wm2 * wn_window[N,P,N]
			    coef[N,N,M,i,j,k] = xz_coef*px[i] + c3_lo - wm2 * wn_window[N,N,M]
			    coef[N,N,P,i,j,k] = xz_coef*px[i] + c3_hi - wm2 * wn_window[N,N,P]

			    coef[N,M,M,i,j,k] = 2*w3a*px[i] - yz_coef*py_lo[j] - yz_coef*pz_lo[k] - wm3 * wn_window[N,M,M]
			    coef[N,M,P,i,j,k] = 2*w3a*px[i] - yz_coef*py_lo[j] - yz_coef*pz_hi[k] - wm3 * wn_window[N,M,P]
			    coef[N,P,M,i,j,k] = 2*w3a*px[i] - yz_coef*py_hi[j] - yz_coef*pz_lo[k] - wm3 * wn_window[N,P,M]
			    coef[N,P,P,i,j,k] = 2*w3a*px[i] - yz_coef*py_hi[j] - yz_coef*pz_hi[k] - wm3 * wn_window[N,P,P]
			    coef[M,N,M,i,j,k] = -xz_coef*px_lo[i] + 2*w3a*py[j] - xz_coef*pz_lo[k] - wm3 * wn_window[M,N,M]
			    coef[M,N,P,i,j,k] = -xz_coef*px_lo[i] + 2*w3a*py[j] - xz_coef*pz_hi[k] - wm3 * wn_window[M,N,P]
			    coef[P,N,M,i,j,k] = -xz_coef*px_hi[i] + 2*w3a*py[j] - xz_coef*pz_lo[k] - wm3 * wn_window[P,N,M]
			    coef[P,N,P,i,j,k] = -xz_coef*px_hi[i] + 2*w3a*py[j] - xz_coef*pz_hi[k] - wm3 * wn_window[P,N,P]
			    coef[M,M,N,i,j,k] = -xy_coef*px_lo[i] - xy_coef*py_lo[j] + 2*w3a*pz[k] - wm3 * wn_window[M,M,N]
			    coef[M,P,N,i,j,k] = -xy_coef*px_lo[i] - xy_coef*py_hi[j] + 2*w3a*pz[k] - wm3 * wn_window[M,P,N]
			    coef[P,M,N,i,j,k] = -xy_coef*px_hi[i] - xy_coef*py_lo[j] + 2*w3a*pz[k] - wm3 * wn_window[P,M,N]
			    coef[P,P,N,i,j,k] = -xy_coef*px_hi[i] - xy_coef*py_hi[j] + 2*w3a*pz[k] - wm3 * wn_window[P,P,N]
			    coef[M,M,M,i,j,k] = -2*w3a*(px_lo[i] + py_lo[j] + pz_lo[k]) - wm4 * wn_window[M,M,M]
			    coef[M,M,P,i,j,k] = -2*w3a*(px_lo[i] + py_lo[j] + pz_hi[k]) - wm4 * wn_window[M,M,P]
			    coef[M,P,M,i,j,k] = -2*w3a*(px_lo[i] + py_hi[j] + pz_lo[k]) - wm4 * wn_window[M,P,M]
			    coef[M,P,P,i,j,k] = -2*w3a*(px_lo[i] + py_hi[j] + pz_hi[k]) - wm4 * wn_window[M,P,P]
			    coef[P,M,M,i,j,k] = -2*w3a*(px_hi[i] + py_lo[j] + pz_lo[k]) - wm4 * wn_window[P,M,M]
			    coef[P,M,P,i,j,k] = -2*w3a*(px_hi[i] + py_lo[j] + pz_hi[k]) - wm4 * wn_window[P,M,P]
			    coef[P,P,M,i,j,k] = -2*w3a*(px_hi[i] + py_hi[j] + pz_lo[k]) - wm4 * wn_window[P,P,M]
			    coef[P,P,P,i,j,k] = -2*w3a*(px_hi[i] + py_hi[j] + pz_hi[k]) - wm4 * wn_window[P,P,P]
			    coef[N,N,N,i,j,k] = -cx*px[i] - cy*py[j] - cz*pz[k] + cNNN*wn_window[N,N,N]

            end
        end
    end
    # Zero out coefficients on the boundary layers
    for i=[M N P]
        for j=[M N P]
            for k=[M N P]
                if i==M
                    coef[i,j,k,1,:,:] = 0
                elseif i==P
                    coef[i,j,k,end,:,:] = 0
                end
                if j==M
                    coef[i,j,k,:,1,:] = 0
                elseif j==P
                    coef[i,j,k,:,end,:] = 0
                end
                if k==M
                    coef[i,j,k,:,:,1] = 0
                elseif k==P
                    coef[i,j,k,:,:,end] = 0
                end
            end
        end
    end

    coef_tup = Array{Array{Complex{F},1},1}()
    offset_tup = Array{I,1}()
    Nel = nx*ny*nz
    t = 1
    for k=[M N P]
        for j=[M N P]
            for i=[M N P]
                idx_start = 1
                idx_end = Nel
                offset = 0
                # For the coefficients before the diagonal, in lexographial organization
                # Remove the coefficient
                if t < 14
                    if i==M
                        idx_start += 1
                        offset -= 1
                    elseif i==P
                        idx_start -= 1
                        offset += 1
                    end
                    if j==M
                        idx_start += nx
                        offset -= nx
                    elseif j==P
                        idx_start -= nx
                        offset += nx
                    end
                    if k==M
                        idx_start += nx*ny
                        offset -= nx*ny
                    end
                elseif t > 14
                    if i==M
                        idx_end += 1
                        offset -= 1
                    elseif i==P
                        idx_end -= 1
                        offset += 1
                    end
                    if j==M
                        idx_end += nx
                        offset -= nx
                    elseif j==P
                        idx_end -= nx
                        offset += nx
                    end
                    if k==P
                        idx_end -= nx*ny
                        offset += nx*ny
                    end
                end
                push!(coef_tup,vec(coef[i,j,k,:,:,:])[idx_start:idx_end])
                push!(offset_tup,offset)

                t += 1
            end
        end
    end

    H = spdiagm(tuple(coef_tup...),tuple(offset_tup...),Nel,Nel)

end

function helm3d_std_7pt(v::Union{AbstractArray{F,3},AbstractArray{Complex{F},3}},Δ::AbstractArray{F,1},n::AbstractArray{I,1},freq::Union{F,Complex{F}},npml::AbstractArray{I,2},x::AbstractArray{Complex{F},3};forw_mode::Bool=true,deriv_mode::Bool=true) where {F<:Real,I<:Integer}
    nx,ny,nz = n[1],n[2],n[3]
    hx,hy,hz = Δ[1],Δ[2],Δ[3]
    hx_isq = 1/hx^2
    hy_isq = 1/hy^2
    hz_isq = 1/hz^2
    M = 1
    N = 2
    P = 3
    ω2 = (2*π*freq)^2
    npxlo,npxhi,npylo,npyhi,npzlo,npzhi = npml[1,1],npml[2,1],npml[1,2],npml[2,2],npml[1,3],npml[2,3]
    maxnpx = max(npxlo,npxhi)
    maxnpy = max(npylo,npyhi)
    maxnpz = max(npzlo,npzhi)
    ξx = x->ξ(x,nx,npxlo,npxhi,maxnpx,freq)
    ξy = y->ξ(y,ny,npylo,npyhi,maxnpy,freq)
    ξz = z->ξ(z,nz,npzlo,npzhi,maxnpz,freq)
    coef = zeros(Complex{F},3,3,3)
    x_window = zeros(Complex{F},3,3,3)
    y = zeros(Complex{F},n...)

    for k=1:nz
        ξzM = ξz(k-0.5)
        ξzN = ξz(k)
        ξzP = ξz(k+0.5)
        ξzlo = ξzM*ξzN*hz_isq
        ξzhi = ξzP*ξzN*hz_isq
        coef[N,N,M] = -ξzlo
        coef[N,N,P] = -ξzhi
        for j=1:ny
            ξyM = ξy(j-0.5)
            ξyN = ξy(j)
            ξyP = ξy(j+0.5)
            ξylo = ξyM*ξyN*hy_isq
            ξyhi = ξyP*ξyN*hy_isq
            coef[N,M,N] = -ξylo
            coef[N,P,N] = -ξyhi
            for i=1:nx
                 for kk=1:3
                    load_z = k+kk-2 > 0 && k+kk-2<=nz;
                    for jj=1:3
                        load_y = j+jj-2 > 0 && j+jj-2<=ny;
                        for ii=1:3
                            load_x = i+ii-2 > 0 && i+ii-2<=nx;
                            if load_x & load_y & load_z
                                x_window[ii,jj,kk] = x[i+ii-2,j+jj-2,k+kk-2]
                            else
                                x_window[ii,jj,kk] = complex(0,0)
                            end
                        end
                    end
                end

                ξxM = ξx(i-0.5)
                ξxN = ξx(i)
                ξxP = ξx(i+0.5)
                ξxlo = ξxM*ξxN*hx_isq
                ξxhi = ξxP*ξxN*hx_isq
                wn2 = ω2*(v[i,j,k]^(-2))
                coef[M,N,N] = -ξxlo
                coef[P,N,N] = -ξxhi
                coef[N,N,N] = ξxlo + ξxhi + ξylo + ξyhi + ξzlo + ξzhi - wn2
                y[i,j,k] = dot(vec(coef),vec(x_window))
            end
        end
    end
    return vec(y)
end

function helm3d_std_7pt_matrix(v::Union{AbstractArray{F,3},AbstractArray{Complex{F},3}},Δ::AbstractArray{F,1},n::AbstractArray{I,1},freq::Union{F,Complex{F}},npml::AbstractArray{I,2};forw_mode::Bool=true,deriv_mode::Bool=true) where {F<:Real,I<:Integer}
    nx,ny,nz = n[1],n[2],n[3]
    hx,hy,hz = Δ[1],Δ[2],Δ[3]
    hx_isq = 1/hx^2
    hy_isq = 1/hy^2
    hz_isq = 1/hz^2
    M = 1
    N = 2
    P = 3
    ω2 = (2*π*freq)^2
    npxlo,npxhi,npylo,npyhi,npzlo,npzhi = npml[1,1],npml[2,1],npml[1,2],npml[2,2],npml[1,3],npml[2,3]
    maxnpx = max(npxlo,npxhi)
    maxnpy = max(npylo,npyhi)
    maxnpz = max(npzlo,npzhi)
    ξx = x->ξ(x,nx,npxlo,npxhi,maxnpx,freq)
    ξy = y->ξ(y,ny,npylo,npyhi,maxnpy,freq)
    ξz = z->ξ(z,nz,npzlo,npzhi,maxnpz,freq)
    coef = Array{Array{Complex{F},3},1}()
    NNM = 1
    NMN = 2
    MNN = 3
    NNN = 4
    PNN = 5
    NPN = 6
    NNP = 7

    for i=1:7
        push!(coef,zeros(Complex{F},nx,ny,nz))
    end

    for k=1:nz
        ξzM = ξz(k-0.5)
        ξzN = ξz(k)
        ξzP = ξz(k+0.5)
        ξzlo = ξzM*ξzN*hz_isq
        ξzhi = ξzP*ξzN*hz_isq
        coef[NNM][:,:,k] = -ξzlo
        coef[NNP][:,:,k] = -ξzhi
        for j=1:ny
            ξyM = ξy(j-0.5)
            ξyN = ξy(j)
            ξyP = ξy(j+0.5)
            ξylo = ξyM*ξyN*hy_isq
            ξyhi = ξyP*ξyN*hy_isq
            coef[NMN][:,j,k] = -ξylo
            coef[NPN][:,j,k] = -ξyhi
            for i=1:nx
                ξxM = ξx(i-0.5)
                ξxN = ξx(i)
                ξxP = ξx(i+0.5)
                ξxlo = ξxM*ξxN*hx_isq
                ξxhi = ξxP*ξxN*hx_isq
                wn2 = ω2*(v[i,j,k]^(-2))
                coef[MNN][i,j,k] = -ξxlo
                coef[PNN][i,j,k] = -ξxhi
                coef[NNN][i,j,k] = ξxlo + ξxhi + ξylo + ξyhi + ξzlo + ξzhi - wn2
            end
        end
    end
    coef[MNN][1,:,:] = 0
    coef[PNN][end,:,:] = 0
    coef[NMN][:,1,:] = 0
    coef[NPN][:,end,:] = 0
    coef[NNM][:,:,1] = 0
    coef[NNP][:,:,end] = 0
    Nel = nx*ny*nz
    H = spdiagm((vec(coef[NNM])[nx*ny+1:end],vec(coef[NMN])[nx+1:end],vec(coef[MNN])[2:end],vec(coef[NNN]),vec(coef[PNN][1:end-1]),vec(coef[NPN][1:end-nx]),vec(coef[NNP][1:end-nx*ny])),(-nx*ny,-nx,-1,0,1,nx,nx*ny),Nel,Nel)

end

function helm3d_chen2012_27pt(v::Union{AbstractArray{F,3},AbstractArray{Complex{F},3}},Δ::AbstractArray{F,1},n::AbstractArray{I,1},freq::Union{F,Complex{F}},npml::AbstractArray{I,2},x::AbstractArray{Complex{F},3};forw_mode::Bool=true,deriv_mode::Bool=true) where {F<:Real,I<:Integer}
    nx,ny,nz = n[1],n[2],n[3]
    hx,hy,hz = Δ[1],Δ[2],Δ[3]
    hx_isq = 1/hx^2
    hy_isq = 1/hy^2
    hz_isq = 1/hz^2
    M = 1
    N = 2
    P = 3
    λmin = minimum(vec(v))/freq
    IGmin = λmin/maximum(Δ)
    ω2 = (2*π*freq)^2
    if IGmin >= 2 && IGmin <= 3.5
        γ1 = 0.5035127
        γ2 = 0.0720630
        γ3 = 0.4244243
        w1 = 0.4058413
        w2 = 0.1966284
        w3 = 0.5979158
        w4 = -0.2003855
    elseif IGmin > 3.5 && IGmin <= 5
        γ1 = 0.7617528
        γ2 = -0.0148152
        γ3 =  0.2530624
        w1 = 0.7602512
        w2 = -0.4883334
        w3 = 1.1153920
        w4 = -0.3873097
    elseif IGmin > 5 && IGmin <= 7
        γ1 = 0.8159342
        γ2 = -0.0340791
        γ3 = 0.2181449
        w1 = 1.1330134
        w2 = -1.5191327
        w3 = 2.1033335
        w4 = 0.7172142
    elseif IGmin > 7 && IGmin <= 9
        γ1 = 0.8354262
        γ2 = -0.0394517
        γ3 = 0.2040255
        w1 = 1.7177071
        w2 = -3.2400262
        w3 = 3.8084643
        w4 = -1.2861453
    elseif IGmin > 9 && IGmin <= 10
        γ1 = 0.8432810
        γ2 = -0.0414069
        γ3 = 0.1981258
        w1 = 2.4693294
        w2 = -5.4811311
        w3 = 6.0429826
        w4 = 2.0311809
    else
        γ1 = 0.8269996
        γ2 = 4.097e-7
        γ3 = 0.1729999
        w1 = 2.9473150
        w2 = -6.8805122
        w3 = 7.4116566
        w4 = -2.4784594
    end

    npxlo,npxhi,npylo,npyhi,npzlo,npzhi = npml[1,1],npml[2,1],npml[1,2],npml[2,2],npml[1,3],npml[2,3]
    maxnpx = max(npxlo,npxhi)
    maxnpy = max(npylo,npyhi)
    maxnpz = max(npzlo,npzhi)
    ξx = x->ξ(x,nx,npxlo,npxhi,maxnpx,freq)
    ξy = y->ξ(y,ny,npylo,npyhi,maxnpy,freq)
    ξz = z->ξ(z,nz,npzlo,npzhi,maxnpz,freq)

    k_sq = zeros(eltype(v),3,3,3)
    zero_t = convert(eltype(v),0)
    coef = zeros(Complex{F},3,3,3)
    x_window = zeros(Complex{F},3,3,3)
    y = zeros(Complex{F},n...)
    for k=1:nz
        ξzM = ξz(k-1)
        ξzNminus = ξz(k-0.5)
        ξzN = ξz(k)
        ξzNplus = ξz(k+0.5)
        ξzP = ξz(k+1)
        for j=1:ny
            ξyM = ξy(k-1)
            ξyNminus = ξy(k-0.5)
            ξyN = ξy(k)
            ξyNplus = ξy(k+0.5)
            ξyP = ξy(k+1)

            for i=1:nx
                ξxM = ξx(i-1)
                ξxNminus = ξx(i-0.5)
                ξxN = ξx(i)
                ξxNplus = ξx(i+0.5)
                ξxP = ξx(i+1)

                # Load 3x3x3 windows around the current point
                for kk=1:3
                    load_z = k+kk-2 > 0 && k+kk-2<=nz;
                    for jj=1:3
                        load_y = j+jj-2 > 0 && j+jj-2<=ny;
                        for ii=1:3
                            load_x = i+ii-2 > 0 && i+ii-2<=nx;
                            if load_x & load_y & load_z
                                k_sq[ii,jj,kk] = ω2*(v[i+ii-2,j+jj-2,k+kk-2])^(-2)
                                x_window[ii,jj,kk] = x[i+ii-2,j+jj-2,k+kk-2]
                            else
                                k_sq[ii,jj,kk] = zero_t
                                x_window[ii,jj,kk] = complex(0,0)
                            end

                        end
                    end
                end

                # Coefficient expressions
                coef[M,M,M] = γ3*hx_isq*ξxNminus*ξyM*ξzM/4 + γ3*hy_isq*ξxM*ξyNminus*ξzM/4 + γ3*hz_isq*ξxM*ξyM*ξzNminus/4 + k_sq[M,M,M]*w4*ξxM*ξyM*ξzM/8
                coef[M,M,N] = γ2*hx_isq*ξxNminus*ξyM*ξzN/4 + γ2*hy_isq*ξxM*ξyNminus*ξzN/4 - γ3*hz_isq*ξxM*ξyM*ξzNminus/4 - γ3*hz_isq*ξxM*ξyM*ξzNplus/4 + k_sq[M,M,N]*w3*ξxM*ξyM*ξzN/12
                coef[M,M,P] = γ3*hx_isq*ξxNminus*ξyM*ξzP/4 + γ3*hy_isq*ξxM*ξyNminus*ξzP/4 + γ3*hz_isq*ξxM*ξyM*ξzNplus/4 + k_sq[M,M,P]*w4*ξxM*ξyM*ξzP/8
                coef[M,N,M] = γ2*hx_isq*ξxNminus*ξyN*ξzM/4 + γ2*hz_isq*ξxM*ξyN*ξzNminus/4 - γ3*hy_isq*ξxM*ξyNminus*ξzM/4 - γ3*hy_isq*ξxM*ξyNplus*ξzM/4 + k_sq[M,N,M]*w3*ξxM*ξyN*ξzM/12
                coef[M,N,N] = γ1*hx_isq*ξxNminus*ξyN*ξzN - γ2*hy_isq*ξxM*ξyNminus*ξzN/4 - γ2*hy_isq*ξxM*ξyNplus*ξzN/4 - γ2*hz_isq*ξxM*ξyN*ξzNminus/4 - γ2*hz_isq*ξxM*ξyN*ξzNplus/4 + k_sq[M,N,N]*w2*ξxM*ξyN*ξzN/6
                coef[M,N,P] = γ2*hx_isq*ξxNminus*ξyN*ξzP/4 + γ2*hz_isq*ξxM*ξyN*ξzNplus/4 - γ3*hy_isq*ξxM*ξyNminus*ξzP/4 - γ3*hy_isq*ξxM*ξyNplus*ξzP/4 + k_sq[M,N,P]*w3*ξxM*ξyN*ξzP/12
                coef[M,P,M] = γ3*hx_isq*ξxNminus*ξyP*ξzM/4 + γ3*hy_isq*ξxM*ξyNplus*ξzM/4 + γ3*hz_isq*ξxM*ξyP*ξzNminus/4 + k_sq[M,P,M]*w4*ξxM*ξyP*ξzM/8
                coef[M,P,N] = γ2*hx_isq*ξxNminus*ξyP*ξzN/4 + γ2*hy_isq*ξxM*ξyNplus*ξzN/4 - γ3*hz_isq*ξxM*ξyP*ξzNminus/4 - γ3*hz_isq*ξxM*ξyP*ξzNplus/4 + k_sq[M,P,N]*w3*ξxM*ξyP*ξzN/12
                coef[M,P,P] = γ3*hx_isq*ξxNminus*ξyP*ξzP/4 + γ3*hy_isq*ξxM*ξyNplus*ξzP/4 + γ3*hz_isq*ξxM*ξyP*ξzNplus/4 + k_sq[M,P,P]*w4*ξxM*ξyP*ξzP/8
                coef[N,M,M] = γ2*hy_isq*ξxN*ξyNminus*ξzM/4 + γ2*hz_isq*ξxN*ξyM*ξzNminus/4 - γ3*hx_isq*ξxNminus*ξyM*ξzM/4 - γ3*hx_isq*ξxNplus*ξyM*ξzM/4 + k_sq[N,M,M]*w3*ξxN*ξyM*ξzM/12
                coef[N,M,N] = γ1*hy_isq*ξxN*ξyNminus*ξzN - γ2*hx_isq*ξxNminus*ξyM*ξzN/4 - γ2*hx_isq*ξxNplus*ξyM*ξzN/4 - γ2*hz_isq*ξxN*ξyM*ξzNminus/4 - γ2*hz_isq*ξxN*ξyM*ξzNplus/4 + k_sq[N,M,N]*w2*ξxN*ξyM*ξzN/6
                coef[N,M,P] = γ2*hy_isq*ξxN*ξyNminus*ξzP/4 + γ2*hz_isq*ξxN*ξyM*ξzNplus/4 - γ3*hx_isq*ξxNminus*ξyM*ξzP/4 - γ3*hx_isq*ξxNplus*ξyM*ξzP/4 + k_sq[N,M,P]*w3*ξxN*ξyM*ξzP/12
                coef[N,N,M] = γ1*hz_isq*ξxN*ξyN*ξzNminus - γ2*hx_isq*ξxNminus*ξyN*ξzM/4 - γ2*hx_isq*ξxNplus*ξyN*ξzM/4 - γ2*hy_isq*ξxN*ξyNminus*ξzM/4 - γ2*hy_isq*ξxN*ξyNplus*ξzM/4 + k_sq[N,N,M]*w2*ξxN*ξyN*ξzM/6
                coef[N,N,N] = -γ1*hx_isq*ξxNminus*ξyN*ξzN - γ1*hx_isq*ξxNplus*ξyN*ξzN - γ1*hy_isq*ξxN*ξyNminus*ξzN - γ1*hy_isq*ξxN*ξyNplus*ξzN - γ1*hz_isq*ξxN*ξyN*ξzNminus - γ1*hz_isq*ξxN*ξyN*ξzNplus + k_sq[N,N,N]*w1*ξxN*ξyN*ξzN
                coef[N,N,P] = γ1*hz_isq*ξxN*ξyN*ξzNplus - γ2*hx_isq*ξxNminus*ξyN*ξzP/4 - γ2*hx_isq*ξxNplus*ξyN*ξzP/4 - γ2*hy_isq*ξxN*ξyNminus*ξzP/4 - γ2*hy_isq*ξxN*ξyNplus*ξzP/4 + k_sq[N,N,P]*w2*ξxN*ξyN*ξzP/6
                coef[N,P,M] = γ2*hy_isq*ξxN*ξyNplus*ξzM/4 + γ2*hz_isq*ξxN*ξyP*ξzNminus/4 - γ3*hx_isq*ξxNminus*ξyP*ξzM/4 - γ3*hx_isq*ξxNplus*ξyP*ξzM/4 + k_sq[N,P,M]*w3*ξxN*ξyP*ξzM/12
                coef[N,P,N] = γ1*hy_isq*ξxN*ξyNplus*ξzN - γ2*hx_isq*ξxNminus*ξyP*ξzN/4 - γ2*hx_isq*ξxNplus*ξyP*ξzN/4 - γ2*hz_isq*ξxN*ξyP*ξzNminus/4 - γ2*hz_isq*ξxN*ξyP*ξzNplus/4 + k_sq[N,P,N]*w2*ξxN*ξyP*ξzN/6
                coef[N,P,P] = γ2*hy_isq*ξxN*ξyNplus*ξzP/4 + γ2*hz_isq*ξxN*ξyP*ξzNplus/4 - γ3*hx_isq*ξxNminus*ξyP*ξzP/4 - γ3*hx_isq*ξxNplus*ξyP*ξzP/4 + k_sq[N,P,P]*w3*ξxN*ξyP*ξzP/12
                coef[P,M,M] = γ3*hx_isq*ξxNplus*ξyM*ξzM/4 + γ3*hy_isq*ξxP*ξyNminus*ξzM/4 + γ3*hz_isq*ξxP*ξyM*ξzNminus/4 + k_sq[P,M,M]*w4*ξxP*ξyM*ξzM/8

                coef[P,M,N] = γ2*hx_isq*ξxNplus*ξyM*ξzN/4 + γ2*hy_isq*ξxP*ξyNminus*ξzN/4 - γ3*hz_isq*ξxP*ξyM*ξzNminus/4 - γ3*hz_isq*ξxP*ξyM*ξzNplus/4 + k_sq[P,M,N]*w3*ξxP*ξyM*ξzN/12
                coef[P,M,P] = γ3*hx_isq*ξxNplus*ξyM*ξzP/4 + γ3*hy_isq*ξxP*ξyNminus*ξzP/4 + γ3*hz_isq*ξxP*ξyM*ξzNplus/4 + k_sq[P,M,P]*w4*ξxP*ξyM*ξzP/8

                coef[P,N,M] = γ2*hx_isq*ξxNplus*ξyN*ξzM/4 + γ2*hz_isq*ξxP*ξyN*ξzNminus/4 - γ3*hy_isq*ξxP*ξyNminus*ξzM/4 - γ3*hy_isq*ξxP*ξyNplus*ξzM/4 + k_sq[P,N,M]*w3*ξxP*ξyN*ξzM/12

                coef[P,N,N] = γ1*hx_isq*ξxNplus*ξyN*ξzN - γ2*hy_isq*ξxP*ξyNminus*ξzN/4 - γ2*hy_isq*ξxP*ξyNplus*ξzN/4 - γ2*hz_isq*ξxP*ξyN*ξzNminus/4 - γ2*hz_isq*ξxP*ξyN*ξzNplus/4 + k_sq[P,N,N]*w2*ξxP*ξyN*ξzN/6

                coef[P,N,P] = γ2*hx_isq*ξxNplus*ξyN*ξzP/4 + γ2*hz_isq*ξxP*ξyN*ξzNplus/4 - γ3*hy_isq*ξxP*ξyNminus*ξzP/4 - γ3*hy_isq*ξxP*ξyNplus*ξzP/4 + k_sq[P,N,P]*w3*ξxP*ξyN*ξzP/12

                coef[P,P,M] = γ3*hx_isq*ξxNplus*ξyP*ξzM/4 + γ3*hy_isq*ξxP*ξyNplus*ξzM/4 + γ3*hz_isq*ξxP*ξyP*ξzNminus/4 + k_sq[P,P,M]*w4*ξxP*ξyP*ξzM/8
                coef[P,P,N] = γ2*hx_isq*ξxNplus*ξyP*ξzN/4 + γ2*hy_isq*ξxP*ξyNplus*ξzN/4 - γ3*hz_isq*ξxP*ξyP*ξzNminus/4 - γ3*hz_isq*ξxP*ξyP*ξzNplus/4 + k_sq[P,P,N]*w3*ξxP*ξyP*ξzN/12
                coef[P,P,P] = γ3*hx_isq*ξxNplus*ξyP*ξzP/4 + γ3*hy_isq*ξxP*ξyNplus*ξzP/4 + γ3*hz_isq*ξxP*ξyP*ξzNplus/4 + k_sq[P,P,P]*w4*ξxP*ξyP*ξzP/8

                y[i,j,k] = dot(vec(coef),vec(x_window))[1]
end
end
end
return vec(y)
end

function helm3d_chen2012_27pt_matrix(v::Union{AbstractArray{F,3},AbstractArray{Complex{F},3}},Δ::AbstractArray{F,1},n::AbstractArray{I,1},freq::Union{F,Complex{F}},npml::AbstractArray{I,2};deriv_mode::Bool=true) where {F<:Real,I<:Integer}
    nx,ny,nz = n[1],n[2],n[3]
    hx,hy,hz = Δ[1],Δ[2],Δ[3]
    hx_isq = 1/hx^2
    hy_isq = 1/hy^2
    hz_isq = 1/hz^2
    M = 1
    N = 2
    P = 3
    λmin = minimum(vec(v))/freq
    IGmin = λmin/maximum(Δ)
    ω2 = (2*π*freq)^2
    if IGmin >= 2 && IGmin <= 3.5
        γ1 = 0.5035127
        γ2 = 0.0720630
        γ3 = 0.4244243
        w1 = 0.4058413
        w2 = 0.1966284
        w3 = 0.5979158
        w4 = -0.2003855
    elseif IGmin > 3.5 && IGmin <= 5
        γ1 = 0.7617528
        γ2 = -0.0148152
        γ3 =  0.2530624
        w1 = 0.7602512
        w2 = -0.4883334
        w3 = 1.1153920
        w4 = -0.3873097
    elseif IGmin > 5 && IGmin <= 7
        γ1 = 0.8159342
        γ2 = -0.0340791
        γ3 = 0.2181449
        w1 = 1.1330134
        w2 = -1.5191327
        w3 = 2.1033335
        w4 = 0.7172142
    elseif IGmin > 7 && IGmin <= 9
        γ1 = 0.8354262
        γ2 = -0.0394517
        γ3 = 0.2040255
        w1 = 1.7177071
        w2 = -3.2400262
        w3 = 3.8084643
        w4 = -1.2861453
    elseif IGmin > 9 && IGmin <= 10
        γ1 = 0.8432810
        γ2 = -0.0414069
        γ3 = 0.1981258
        w1 = 2.4693294
        w2 = -5.4811311
        w3 = 6.0429826
        w4 = 2.0311809
    else
        γ1 = 0.8269996
        γ2 = 4.097e-7
        γ3 = 0.1729999
        w1 = 2.9473150
        w2 = -6.8805122
        w3 = 7.4116566
        w4 = -2.4784594
    end

    npxlo,npxhi,npylo,npyhi,npzlo,npzhi = npml[1,1],npml[2,1],npml[1,2],npml[2,2],npml[1,3],npml[2,3]
    maxnpx = max(npxlo,npxhi)
    maxnpy = max(npylo,npyhi)
    maxnpz = max(npzlo,npzhi)
    ξx = x->ξ(x,nx,npxlo,npxhi,maxnpx,freq)
    ξy = y->ξ(y,ny,npylo,npyhi,maxnpy,freq)
    ξz = z->ξ(z,nz,npzlo,npzhi,maxnpz,freq)

    k_sq = zeros(eltype(v),3,3,3)
    zero_t = convert(eltype(v),0)
    coef = zeros(Complex{F},3,3,3,n[1],n[2],n[3])
    x_window = zeros(Complex{F},3,3,3)

    for k=1:nz
        ξzM = ξz(k-1)
        ξzNminus = ξz(k-0.5)
        ξzN = ξz(k)
        ξzNplus = ξz(k+0.5)
        ξzP = ξz(k+1)
        for j=1:ny
            ξyM = ξy(k-1)
            ξyNminus = ξy(k-0.5)
            ξyN = ξy(k)
            ξyNplus = ξy(k+0.5)
            ξyP = ξy(k+1)
            for i=1:nx
                ξxM = ξx(i-1)
                ξxNminus = ξx(i-0.5)
                ξxN = ξx(i)
                ξxNplus = ξx(i+0.5)
                ξxP = ξx(i+1)

                # Load 3x3x3 windows around the current point
                for kk=1:3
                    load_z = k+kk-2 > 0 && k+kk-2<=nz;
                    for jj=1:3
                        load_y = j+jj-2 > 0 && j+jj-2<=ny;
                        for ii=1:3
                            load_x = i+ii-2 > 0 && i+ii-2<=nx;
                            if load_x & load_y & load_z
                                k_sq[ii,jj,kk] = ω2*(v[i+ii-2,j+jj-2,k+kk-2]^(-2))
                            else
                                k_sq[ii,jj,kk] = zero_t
                            end
                        end
                    end
                end

                # Coefficient expressions
                coef[M,M,M,i,j,k] = γ3*hx_isq*ξxNminus*ξyM*ξzM/4 + γ3*hy_isq*ξxM*ξyNminus*ξzM/4 + γ3*hz_isq*ξxM*ξyM*ξzNminus/4 + k_sq[M,M,M]*w4*ξxM*ξyM*ξzM/8
                coef[M,M,N,i,j,k] = γ2*hx_isq*ξxNminus*ξyM*ξzN/4 + γ2*hy_isq*ξxM*ξyNminus*ξzN/4 - γ3*hz_isq*ξxM*ξyM*ξzNminus/4 - γ3*hz_isq*ξxM*ξyM*ξzNplus/4 + k_sq[M,M,N]*w3*ξxM*ξyM*ξzN/12
                coef[M,M,P,i,j,k] = γ3*hx_isq*ξxNminus*ξyM*ξzP/4 + γ3*hy_isq*ξxM*ξyNminus*ξzP/4 + γ3*hz_isq*ξxM*ξyM*ξzNplus/4 + k_sq[M,M,P]*w4*ξxM*ξyM*ξzP/8
                coef[M,N,M,i,j,k] = γ2*hx_isq*ξxNminus*ξyN*ξzM/4 + γ2*hz_isq*ξxM*ξyN*ξzNminus/4 - γ3*hy_isq*ξxM*ξyNminus*ξzM/4 - γ3*hy_isq*ξxM*ξyNplus*ξzM/4 + k_sq[M,N,M]*w3*ξxM*ξyN*ξzM/12
                coef[M,N,N,i,j,k] = γ1*hx_isq*ξxNminus*ξyN*ξzN - γ2*hy_isq*ξxM*ξyNminus*ξzN/4 - γ2*hy_isq*ξxM*ξyNplus*ξzN/4 - γ2*hz_isq*ξxM*ξyN*ξzNminus/4 - γ2*hz_isq*ξxM*ξyN*ξzNplus/4 + k_sq[M,N,N]*w2*ξxM*ξyN*ξzN/6
                coef[M,N,P,i,j,k] = γ2*hx_isq*ξxNminus*ξyN*ξzP/4 + γ2*hz_isq*ξxM*ξyN*ξzNplus/4 - γ3*hy_isq*ξxM*ξyNminus*ξzP/4 - γ3*hy_isq*ξxM*ξyNplus*ξzP/4 + k_sq[M,N,P]*w3*ξxM*ξyN*ξzP/12
                coef[M,P,M,i,j,k] = γ3*hx_isq*ξxNminus*ξyP*ξzM/4 + γ3*hy_isq*ξxM*ξyNplus*ξzM/4 + γ3*hz_isq*ξxM*ξyP*ξzNminus/4 + k_sq[M,P,M]*w4*ξxM*ξyP*ξzM/8
                coef[M,P,N,i,j,k] = γ2*hx_isq*ξxNminus*ξyP*ξzN/4 + γ2*hy_isq*ξxM*ξyNplus*ξzN/4 - γ3*hz_isq*ξxM*ξyP*ξzNminus/4 - γ3*hz_isq*ξxM*ξyP*ξzNplus/4 + k_sq[M,P,N]*w3*ξxM*ξyP*ξzN/12
                coef[M,P,P,i,j,k] = γ3*hx_isq*ξxNminus*ξyP*ξzP/4 + γ3*hy_isq*ξxM*ξyNplus*ξzP/4 + γ3*hz_isq*ξxM*ξyP*ξzNplus/4 + k_sq[M,P,P]*w4*ξxM*ξyP*ξzP/8
                coef[N,M,M,i,j,k] = γ2*hy_isq*ξxN*ξyNminus*ξzM/4 + γ2*hz_isq*ξxN*ξyM*ξzNminus/4 - γ3*hx_isq*ξxNminus*ξyM*ξzM/4 - γ3*hx_isq*ξxNplus*ξyM*ξzM/4 + k_sq[N,M,M]*w3*ξxN*ξyM*ξzM/12
                coef[N,M,N,i,j,k] = γ1*hy_isq*ξxN*ξyNminus*ξzN - γ2*hx_isq*ξxNminus*ξyM*ξzN/4 - γ2*hx_isq*ξxNplus*ξyM*ξzN/4 - γ2*hz_isq*ξxN*ξyM*ξzNminus/4 - γ2*hz_isq*ξxN*ξyM*ξzNplus/4 + k_sq[N,M,N]*w2*ξxN*ξyM*ξzN/6
                coef[N,M,P,i,j,k] = γ2*hy_isq*ξxN*ξyNminus*ξzP/4 + γ2*hz_isq*ξxN*ξyM*ξzNplus/4 - γ3*hx_isq*ξxNminus*ξyM*ξzP/4 - γ3*hx_isq*ξxNplus*ξyM*ξzP/4 + k_sq[N,M,P]*w3*ξxN*ξyM*ξzP/12
                coef[N,N,M,i,j,k] = γ1*hz_isq*ξxN*ξyN*ξzNminus - γ2*hx_isq*ξxNminus*ξyN*ξzM/4 - γ2*hx_isq*ξxNplus*ξyN*ξzM/4 - γ2*hy_isq*ξxN*ξyNminus*ξzM/4 - γ2*hy_isq*ξxN*ξyNplus*ξzM/4 + k_sq[N,N,M]*w2*ξxN*ξyN*ξzM/6
                coef[N,N,N,i,j,k] = -γ1*hx_isq*ξxNminus*ξyN*ξzN - γ1*hx_isq*ξxNplus*ξyN*ξzN - γ1*hy_isq*ξxN*ξyNminus*ξzN - γ1*hy_isq*ξxN*ξyNplus*ξzN - γ1*hz_isq*ξxN*ξyN*ξzNminus - γ1*hz_isq*ξxN*ξyN*ξzNplus + k_sq[N,N,N]*w1*ξxN*ξyN*ξzN
                coef[N,N,P,i,j,k] = γ1*hz_isq*ξxN*ξyN*ξzNplus - γ2*hx_isq*ξxNminus*ξyN*ξzP/4 - γ2*hx_isq*ξxNplus*ξyN*ξzP/4 - γ2*hy_isq*ξxN*ξyNminus*ξzP/4 - γ2*hy_isq*ξxN*ξyNplus*ξzP/4 + k_sq[N,N,P]*w2*ξxN*ξyN*ξzP/6
                coef[N,P,M,i,j,k] = γ2*hy_isq*ξxN*ξyNplus*ξzM/4 + γ2*hz_isq*ξxN*ξyP*ξzNminus/4 - γ3*hx_isq*ξxNminus*ξyP*ξzM/4 - γ3*hx_isq*ξxNplus*ξyP*ξzM/4 + k_sq[N,P,M]*w3*ξxN*ξyP*ξzM/12
                coef[N,P,N,i,j,k] = γ1*hy_isq*ξxN*ξyNplus*ξzN - γ2*hx_isq*ξxNminus*ξyP*ξzN/4 - γ2*hx_isq*ξxNplus*ξyP*ξzN/4 - γ2*hz_isq*ξxN*ξyP*ξzNminus/4 - γ2*hz_isq*ξxN*ξyP*ξzNplus/4 + k_sq[N,P,N]*w2*ξxN*ξyP*ξzN/6
                coef[N,P,P,i,j,k] = γ2*hy_isq*ξxN*ξyNplus*ξzP/4 + γ2*hz_isq*ξxN*ξyP*ξzNplus/4 - γ3*hx_isq*ξxNminus*ξyP*ξzP/4 - γ3*hx_isq*ξxNplus*ξyP*ξzP/4 + k_sq[N,P,P]*w3*ξxN*ξyP*ξzP/12
                coef[P,M,M,i,j,k] = γ3*hx_isq*ξxNplus*ξyM*ξzM/4 + γ3*hy_isq*ξxP*ξyNminus*ξzM/4 + γ3*hz_isq*ξxP*ξyM*ξzNminus/4 + k_sq[P,M,M]*w4*ξxP*ξyM*ξzM/8

                coef[P,M,N,i,j,k] = γ2*hx_isq*ξxNplus*ξyM*ξzN/4 + γ2*hy_isq*ξxP*ξyNminus*ξzN/4 - γ3*hz_isq*ξxP*ξyM*ξzNminus/4 - γ3*hz_isq*ξxP*ξyM*ξzNplus/4 + k_sq[P,M,N]*w3*ξxP*ξyM*ξzN/12
                coef[P,M,P,i,j,k] = γ3*hx_isq*ξxNplus*ξyM*ξzP/4 + γ3*hy_isq*ξxP*ξyNminus*ξzP/4 + γ3*hz_isq*ξxP*ξyM*ξzNplus/4 + k_sq[P,M,P]*w4*ξxP*ξyM*ξzP/8

                coef[P,N,M,i,j,k] = γ2*hx_isq*ξxNplus*ξyN*ξzM/4 + γ2*hz_isq*ξxP*ξyN*ξzNminus/4 - γ3*hy_isq*ξxP*ξyNminus*ξzM/4 - γ3*hy_isq*ξxP*ξyNplus*ξzM/4 + k_sq[P,N,M]*w3*ξxP*ξyN*ξzM/12

coef[P,N,N,i,j,k] = γ1*hx_isq*ξxNplus*ξyN*ξzN - γ2*hy_isq*ξxP*ξyNminus*ξzN/4 - γ2*hy_isq*ξxP*ξyNplus*ξzN/4 - γ2*hz_isq*ξxP*ξyN*ξzNminus/4 - γ2*hz_isq*ξxP*ξyN*ξzNplus/4 + k_sq[P,N,N]*w2*ξxP*ξyN*ξzN/6

                coef[P,N,P,i,j,k] = γ2*hx_isq*ξxNplus*ξyN*ξzP/4 + γ2*hz_isq*ξxP*ξyN*ξzNplus/4 - γ3*hy_isq*ξxP*ξyNminus*ξzP/4 - γ3*hy_isq*ξxP*ξyNplus*ξzP/4 + k_sq[P,N,P]*w3*ξxP*ξyN*ξzP/12

                coef[P,P,M,i,j,k] = γ3*hx_isq*ξxNplus*ξyP*ξzM/4 + γ3*hy_isq*ξxP*ξyNplus*ξzM/4 + γ3*hz_isq*ξxP*ξyP*ξzNminus/4 + k_sq[P,P,M]*w4*ξxP*ξyP*ξzM/8
                coef[P,P,N,i,j,k] = γ2*hx_isq*ξxNplus*ξyP*ξzN/4 + γ2*hy_isq*ξxP*ξyNplus*ξzN/4 - γ3*hz_isq*ξxP*ξyP*ξzNminus/4 - γ3*hz_isq*ξxP*ξyP*ξzNplus/4 + k_sq[P,P,N]*w3*ξxP*ξyP*ξzN/12
                coef[P,P,P,i,j,k] = γ3*hx_isq*ξxNplus*ξyP*ξzP/4 + γ3*hy_isq*ξxP*ξyNplus*ξzP/4 + γ3*hz_isq*ξxP*ξyP*ξzNplus/4 + k_sq[P,P,P]*w4*ξxP*ξyP*ξzP/8

end
end
end

# Zero out coefficients on the boundary layers
for i=[M N P]
    for j=[M N P]
        for k=[M N P]
            if i==M
                coef[i,j,k,1,:,:] = 0
            elseif i==P
                coef[i,j,k,end,:,:] = 0
            end
            if j==M
                coef[i,j,k,:,1,:] = 0
            elseif j==P
                coef[i,j,k,:,end,:] = 0
            end
            if k==M
                coef[i,j,k,:,:,1] = 0
            elseif k==P
                coef[i,j,k,:,:,end] = 0
            end
        end
    end
end

coef_tup = Array{Array{Complex{F},1},1}()
offset_tup = Array{I,1}()
Nel = nx*ny*nz
t = 1
for k=[M N P]
    for j=[M N P]
        for i=[M N P]
            idx_start = 1
            idx_end = Nel
            offset = 0
            # For the coefficients before the diagonal, in lexographial organization
            # Remove the coefficient
            if t < 14
                if i==M
                    idx_start += 1
                    offset -= 1
                elseif i==P
                    idx_start -= 1
                    offset += 1
                end
                if j==M
                    idx_start += nx
                    offset -= nx
                elseif j==P
                    idx_start -= nx
                    offset += nx
                end
                if k==M
                    idx_start += nx*ny
                    offset -= nx*ny
                end
            elseif t > 14
                if i==M
                    idx_end += 1
                    offset -= 1
                elseif i==P
                    idx_end -= 1
                    offset += 1
                end
                if j==M
                    idx_end += nx
                    offset -= nx
                elseif j==P
                    idx_end -= nx
                    offset += nx
                end
                if k==P
                    idx_end -= nx*ny
                    offset += nx*ny
                end
            end
            push!(coef_tup,vec(coef[i,j,k,:,:,:])[idx_start:idx_end])
            push!(offset_tup,offset)

            t += 1
        end
    end
end

H = spdiagm(tuple(coef_tup...),tuple(offset_tup...),Nel,Nel)

end

function ξ(x,nx,nplo,nphi,maxnp,freq)
    C = 10;
    dist_from_int = (x .< (nplo+1)) .* (nplo+1-x) + (x .>= nx-nphi) .* (x-(nx-nphi));
    sigma = C * (dist_from_int/maxnp).^2;
    func = complex(1.0,sigma/freq);
end

function pml_func(nx::Int,np_lo::Int,np_hi::Int)

    s_hx = 1.0/(nx+1)
    xix = zeros(nx,2)
    Lx_hi = np_hi/nx
    Lx_lo = np_lo/nx
    gamma = zeros(nx+2,1)
    gamma[1:np_lo] = cos.((pi*(0:np_lo-1) * s_hx)/(2*Lx_lo))
    gamma[np_lo+1:nx+2-np_hi] .= 0
    gamma[nx+3-np_hi:nx+2] = cos.((pi*(1 .- (nx+2-np_hi:nx+1)*s_hx))/(2*Lx_hi))
    return 2.0 ./ ((1 .+ im*gamma[2:nx+1]).*(1 .+ im*gamma[2:nx+1] .+ 1 .+ im*gamma[1:nx] )), 2.0 ./ ((1 .+ im*gamma[2:nx+1]).*(1 .+ im*gamma[2:nx+1] .+ 1 .+ im*gamma[3:nx .+ 2] ))

end

ext: jl


text:
module TestRand

using Phylo
using DataFrames
using Random

using Test

@testset "Nonultrametric()" begin
    # Create a 10 tip tree
    nu = Nonultrametric(10)
    @test eltype(nu) == RootedTree
    @test validate!(rand(nu))
    @test Set(getleafnames(rand(nu))) == Set(getleafnames(rand(nu)))
    # Create a tree with named tips
    species = ["Dog", "Cat", "Human"]
    t = rand(Nonultrametric(species))
    @test ntrees(t) == 1
    @test validate!(t)
    @test Set(getleafnames(t)) == Set(species)

    t2 = rand(Nonultrametric{BinaryTree{OneRoot, DataFrame, Vector{Float64}}}(species))
    @test length(getnodedata(t2, species[1])) == 0
    t3 = rand(Nonultrametric{BinaryTree{OneRoot, DataFrame, Vector{String}}}(species))
    map(nodenameiter(t3)) do name
        setnodedata!(t3, name, nodehistory(t3, name))
    end
    for name in nodenameiter(t3)
        @test all(getnodedata(t3, name) .== nodehistory(t3, name))
    end
    d = rand(BrownianTrait(t, "trait", σ² = 1.0))
    @test d ≢ t
    @test species ⊆ collect(keys(d))
    @test eltype(values(d)) ≡ typeof(1.0)
    t4 = rand!(BrownianTrait(t, "trait64", σ² = 1.0), t)
    @test t4 ≡ t
    @test typeof(getnodedata(t4, species[1])["trait64"]) ≡ typeof(1.0)
    t5 = rand!(BrownianTrait(t, "trait32", 0.0f0, σ² = 1.0f0), t)
    @test typeof(getnodedata(t5, species[1])["trait32"]) ≡ typeof(1.0f0)
    @test t5 ≡ t
end

@testset "Ultrametric()" begin
    # Create a 10 tip tree
    u = Ultrametric(10)
    @test validate!(rand(u))
    @test Set(getleafnames(rand(u))) == Set(getleafnames(rand(u)))
    tree = rand(u)
    @test ntrees(tree) == 1
    heights = map(x -> getheight(tree, x), getleafnames(tree))
    @test all(h -> h ≈ heights[1], heights)
    species = ["Dog", "Cat", "Human"]
    t = rand(Ultrametric(species))
    @test validate!(t)
    @test Set(getleafnames(t)) == Set(species)

    numnodes = 50
    ul = Ultrametric{BinaryTree{OneRoot, DataFrame, Vector{Float64}}}(numnodes)
    u2 = rand(ul)
    @test eltype(ul) == BinaryTree{OneRoot, DataFrame, Vector{Float64}}
    @test length(nodefilter(isinternal, u2)) == numnodes - 2

    d = rand(BrownianTrait(t, "trait", σ² = 1.0))
    @test species ⊆ collect(keys(d))
    @test eltype(values(d)) ≡ typeof(1.0)
    t4 = rand!(BrownianTrait(t, "trait64", σ² = 1.0), t)
    @test t4 ≡ t
    @test typeof(getnodedata(t4, species[1])["trait64"]) ≡ typeof(1.0)
    t5 = rand!(BrownianTrait(t, "trait32", 0.0f0, σ² = 1.0f0), t)
    @test typeof(getnodedata(t5, species[1])["trait32"]) ≡ typeof(1.0f0)
    @test t5 ≡ t
end

@testset "TestSets" begin
    @test ntrees(rand(Ultrametric(10), 10)) == 10
    @test length(rand(Ultrametric(10), 10)) == 10
    names = ["One", "Two", "Three"]
    ts = rand(Nonultrametric(names), 20)
    @test length(ts) == 20
    @test ntrees(ts) == 20
    @test length(getleafnames(ts)) == length(names)
    @test length(getnodenames(first(gettrees(ts)))) .== length(names) * 2 - 1
end
end


text:
type GiveThisThingAType
  cxx
end

type VarVector
  cxx
end

function post_process(boxes::GiveThisThingAType, nsamples::Int)
  preimage_samples = HyperBox{Float64}[]
  @assert nsamples == icxx"$(boxes.cxx).size();"
  for j = 0:(nsamples-1)
    numdims = icxx"$(boxes.cxx)[$j].size();"
    intervals = Interval{Float64}[Interval(icxx"$(boxes.cxx)[$j][$i].lb();",icxx"$(boxes.cxx)[$j][$i].ub();") for i = 0:(numdims-1)]
    push!(preimage_samples, HyperBox(intervals))
  end
  preimage_samples
end

function pre_tlmh(lmap::LiteralMap, cnf::CMCNF, ω::IBEX.ExprSymbol, aux_vars::VarVector, box, nsamples::Int)
  println("Calling C++ pre_tlmh")
  cxx_boxes = GiveThisThingAType(@cxx sigma::pre_tlmh(lmap.cxx, cnf.cxx, ω.cxx, aux_vars.cxx, box, nsamples))
  post_process(cxx_boxes, nsamples)
end

# Construct an initial box (flattened) for all the variables
function build_init_box(Y::RandVar{Bool}, aux_vars::VarSet)
  println("Building Initial Box")
  omega_numdims = maximum(dims(Y))+1
  aux_numdims = length(aux_vars)
  numdims = omega_numdims + aux_numdims
  box = @cxx ibex::IntervalVector(numdims)
  println("Building initial box with $numdims dimensions")
  for i = 0:(numdims-1)
    # Omega is in [0,1] other variables are unbounded
    if i < omega_numdims
      icxx"$box[$i] = ibex::Interval(0,1);"
    else
      icxx"$box[$i] = ibex::Interval::ALL_REALS;"
    end
  end
  box
end

cxx"""void add_to_sp_vec(std::vector<std::shared_ptr<ibex::ExprSymbol>> &vec, ExprSymbol &x) {
  std::shared_ptr<ibex::ExprSymbol> new_ptr(&x);
  vec.push_back(new_ptr);
}
"""

function varset_to_varvec(aux_vars::VarSet)
  println("Creating Vector of Auxiary Variables")
  @show aux_vars
  varvec = VarVector(icxx"std::vector<std::shared_ptr<ibex::ExprSymbol>>();")
  for expr_symbol in aux_vars
    @cxx add_to_sp_vec(varvec.cxx,expr_symbol.cxx)
  end
  varvec
end

function pre_tlmh(Y::RandVar{Bool}, nsamples::Int, solver::Type{SigmaSolver}; args...)
  cmap, cnf, ω, aux_vars = analyze(Y)
  lmap = to_cxx_lmap(cmap)
  pre_tlmh(lmap, cnf, ω, varset_to_varvec(aux_vars), build_init_box(Y, aux_vars), nsamples)
end


text:
module XLSXasJSON

using Printf, REPL
using JSON3, StructTypes
using JSONPointer
using JSONPointer: Pointer
using XLSX
using OrderedCollections

include("index.jl")
include("jsonpointer.jl")
include("worksheet.jl")
include("workbook.jl")
include("writer.jl")

export JSONWorkbook, JSONWorksheet,
        hassheet, sheetnames,
        xlsxpath

end # module


text:
module Center_Out_Test

using Intan, FactCheck

m=Task_COT()

facts() do
    @fact m.tar --> (1.0,0.0)
    @fact m.success --> 0
    @fact m.total --> 0
    @fact m.state --> 0
    @fact m.hold_time --> 300.0
    @fact m.juice_amount --> 1.0
    @fact m.cursor_size --> 5.0
    @fact m.target_size --> 5.0
    @fact m.reach_radius --> 5.0
    @fact m.it_time --> 1.0
end

end


text:
include("splinesSerial.jl")

Su0 = larBezier(S1, [[0,0,0],[10,0,0]])
Su1 = larBezier(S1, [[0,10,0],[2.5,10,3],[5,10,-3],[7.5,10,3],[10,10,0]])
Sv0 = larBezier(S2, [[0,0,0],[0,0,3],[0,10,3],[0,10,0]])
Sv1 = larBezier(S2, [[10,0,0],[10,5,3],[10,10,0]])
coons = larCoonsPatch([Su0,Su1,Sv0,Sv1])
point = [0.2, 0.7]
coons(point)


text:
using SIIPExamples
pkgpath = dirname(dirname(pathof(SIIPExamples)))
include(joinpath(pkgpath, "test", "2_PowerSystems_examples", "parse_matpower.jl"))

FORECASTS_DIR = joinpath(base_dir, "forecasts", "5bus_ts")
fname = joinpath(FORECASTS_DIR, "timeseries_pointers_da.json")
open(fname, "r") do f
    for line in eachline(f)
        println(line)
    end
end

ts_pointers = PowerSystems.IS.read_time_series_file_metadata(fname)

add_time_series!(sys, ts_pointers)
sys

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
struct NotImplementedOperation <: AbstractGroupOperation end

struct NotImplementedManifold <: Manifold{ℝ} end

struct NotImplementedGroupDecorator{M} <: AbstractDecoratorManifold{ℝ}
    manifold::M
end


text:
if UInt64Digit[]
    const ArbDigit = UInt64
    const DoubleArbDigit = UInt128
    const SignedDoubleArbDigit = Int128
else
    const ArbDigit = UInt32
    const DoubleArbDigit = UInt64
    const SignedDoubleArbDigit = Int64
end

const BITS = UInt8(8 * sizeof(ArbDigit))
const HALF_BITS = BITS ÷ 2
const HALF = ArbDigit((1 << HALF_BITS) - 1)
const LO_MASK = DoubleArbDigit(typemax(ArbDigit))
const MAX = ArbDigit(LO_MASK)

highBits(n) = ArbDigit(n >> BITS)
lowBits(n) = ArbDigit(n & LO_MASK)

fromDoubleArbDigit(n) = (highBits(n), lowBits(n))
toDoubleArbDigit(hi, lo) = DoubleArbDigit(lo) | (DoubleArbDigit(hi) << BITS)

function u32_to_u128(a::UInt32,b::UInt32,c::UInt32,d::UInt32)
    UInt128(d) | (UInt128(c) << 32) | (UInt128(c) << 64) | UInt128(c) << 96
end

function u32_from_u128(n)
    (
        (n >> 96) % UInt32,
        (n >> 64) % UInt32,
        (n >> 32) % UInt32,
        n % UInt32
    )
end


text:
using Bezier
using Base.Test

@testset "Bezier" begin
    pts = [0.0 1.0 2.0 3.0; 0.0 -1.0 2.0 0.0]
    t = 0.0:0.25:1.0
    @test all(evalbezier(pts, t) .≈ [0.0 0.75 1.5 2.25 3.0; 0.0 -0.140625 0.375 0.703125 0.0])
    t = 0.5
    β₁, β₂ = subdivide(pts, t)
    @test all(β₁ .≈ [0.0 0.5 1.0 1.5; 0.0 -0.5 0.0 0.375])
    @test all(β₂ .≈ [1.5 2.0 2.5 3.0; 0.375 0.75 1.0 0.0])

    pts = [0.0 1.0 2.0 3.0; 0.0 -1.0 2.0 0.0; 0.0 2.0 -1.0 -1.0]
    t = 0.0:0.5:1.0
    @test all(evalbezier(pts, t) .≈ [0.0 1.5 3.0; 0.0 0.375 0.0; 0.0 0.25 -1.0])
    t = 0.25
    β₁, β₂ = subdivide(pts, t)
    @test all(β₁ .≈ [0.0 0.25 0.5 0.75; 0.0 -0.25 -0.25 -0.140625; 0.0 0.5 0.6875 0.6875])
    @test all(β₂ .≈ [0.75 1.5 2.25 3.0; -0.140625 0.1875 1.5 0.0; 0.6875 0.6875 -1.0 -1.0])
end


text:
#
# SQL fragments.
#

# Names.

sql_name(name::AbstractString) =
    "\"$(replace(name, "\"" => "\"\""))\""

_in_search_path(qname) =
    qname[1] == "pg_catalog" || qname[1] == "public"

sql_name(qname::Tuple) =
    join(sql_name.(_in_search_path(qname) ? Base.tail(qname) : qname), '.')

sql_name(names::AbstractVector) =
    join(sql_name.(names), ", ")


text:
module ANYmalRobot

using RigidBodyDynamics

packagepath() = joinpath(@__DIR__, "..")
urdfpath() = joinpath(packagepath(), "anymal_b_simple_description/urdf/anymal.urdf")

function mechanism(::Type{T} = Float64;
                   floating = true,
                   # contactmodel = default_contact_model(),
                   remove_fixed_tree_joints = true, add_flat_ground = false) where {T}
    mechanism = RigidBodyDynamics.parse_urdf(urdfpath(); scalar_type = T, floating = floating, remove_fixed_tree_joints = remove_fixed_tree_joints)

    remove_fixed_tree_joints && remove_fixed_tree_joints!(mechanism)

    return mechanism
end

end # module


text:
function ranks(v::AbstractArray)
  n     = length(v)
  place = sortperm(v)
  ord   = Array(Int64, n)

  for(i in 1:n)
    ord[place[i]] = i
  end

  ord
end

function factorial(n::Integer)
  prod(1:n)
end

function combinations(k::Integer, n::Integer)
  factorial(n) / (factorial(k) * factorial(n-k))
end


text:
TEXT_EDIT_DEFAULTS = @compat Dict{Symbol, Any}(
:Default => @compat Dict{Symbol, Any}(
	:model           => eye(Mat4),
	:screen          => ROOT_SCREEN

))
edit(text::String, style=Style(:Default); customization...) = edit(style, text, mergedefault!(style, TEXT_EDIT_DEFAULTS, customization))
edit(text::Texture{GLGlyph{Uint16}, 4, 2}, obj::RenderObject, style=Style(:Default); customization...) = edit(style, text, obj, mergedefault!(style, TEXT_EDIT_DEFAULTS, customization))

MATRIX_EDITING_DEFAULTS = @compat(Dict(
:Default => @compat(Dict(
	:offset 	     => Vec2(0),
	:backgroundcolor => rgba(0,0,0.4, 0.3),
	:color           => rgba(248.0/255.0, 248.0/255.0,242.0/255.0, 1.0),
	:maxdigits       => 5f0,
	:maxlength       => 10f0,
	:model           => eye(Mat4),
	:screen          => ROOT_SCREEN,
	:font            => getfont(),
	:newline         => -Vec3(0, getfont().props[1][2]*1.5f0, 0),
	:advance         => Vec3(getfont().props[1][1], 0, 0)
))
))
# High Level text rendering for one line or multi line text, which is decided by searching for the occurence of '\n' in text
# Low level text rendering for one line text
# Low level text rendering for multiple line text
edit{T <: Union(AbstractFixedVector, Real)}(numbers::Texture{T, 1, 2}, style=Style(:Default); customization...) = edit(style, numbers, mergedefault!(style, MATRIX_EDITING_DEFAULTS, customization))
edit{T <: Union(AbstractVector, Real)}(numbers::Input{T}, style=Style(:Default); customization...) = edit(numbers.value, style; customization...)
edit{T <: Union(AbstractVector, Real)}(numbers::T, style=Style(:Default); customization...) = edit(style, numbers, mergedefault!(style, MATRIX_EDITING_DEFAULTS, customization))

RENDEROBJECT_DEFAULTS = @compat Dict{Symbol, Any}(
:Default => @compat Dict{Symbol, Any}(
  :model           => eye(Mat4),
  :screen          => ROOT_SCREEN
))
edit(obj::RenderObject, style=Style(:Default); customizations...) = edit(style, obj, mergedefault!(style, RENDEROBJECT_DEFAULTS, customizations))

#################################################################################################################################
# Color Editing:
edit(color::AbstractAlphaColorValue, style::Style=Style(:Default); customization...) = visualize(color, style; customization...)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiNetworkingV1beta1Ingress <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiNetworkingV1beta1IngressSpec } # spec name: spec
    status::Any # spec type: Union{ Nothing, IoK8sApiNetworkingV1beta1IngressStatus } # spec name: status

    function IoK8sApiNetworkingV1beta1Ingress(;apiVersion=nothing, kind=nothing, metadata=nothing, spec=nothing, status=nothing)
        o = new()
        validate_property(IoK8sApiNetworkingV1beta1Ingress, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiNetworkingV1beta1Ingress, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiNetworkingV1beta1Ingress, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiNetworkingV1beta1Ingress, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        validate_property(IoK8sApiNetworkingV1beta1Ingress, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        o
    end
end # type IoK8sApiNetworkingV1beta1Ingress

const _property_map_IoK8sApiNetworkingV1beta1Ingress = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"), Symbol("status")=>Symbol("status"))
const _property_types_IoK8sApiNetworkingV1beta1Ingress = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiNetworkingV1beta1IngressSpec", Symbol("status")=>"IoK8sApiNetworkingV1beta1IngressStatus")
Base.propertynames(::Type{ IoK8sApiNetworkingV1beta1Ingress }) = collect(keys(_property_map_IoK8sApiNetworkingV1beta1Ingress))
Swagger.property_type(::Type{ IoK8sApiNetworkingV1beta1Ingress }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiNetworkingV1beta1Ingress[name]))}
Swagger.field_name(::Type{ IoK8sApiNetworkingV1beta1Ingress }, property_name::Symbol) =  _property_map_IoK8sApiNetworkingV1beta1Ingress[property_name]

function check_required(o::IoK8sApiNetworkingV1beta1Ingress)
    true
end

function validate_property(::Type{ IoK8sApiNetworkingV1beta1Ingress }, name::Symbol, val)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceValidation <: SwaggerModel
    openAPIV3Schema::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps } # spec name: openAPIV3Schema

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceValidation(;openAPIV3Schema=nothing)
        o = new()
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceValidation, Symbol("openAPIV3Schema"), openAPIV3Schema)
        setfield!(o, Symbol("openAPIV3Schema"), openAPIV3Schema)
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceValidation

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceValidation = Dict{Symbol,Symbol}(Symbol("openAPIV3Schema")=>Symbol("openAPIV3Schema"))
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceValidation = Dict{Symbol,String}(Symbol("openAPIV3Schema")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaProps")
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceValidation }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceValidation))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceValidation }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceValidation[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceValidation }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceValidation[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceValidation)
    true
end

function validate_property(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceValidation }, name::Symbol, val)
end


text:
struct Receiver_B0 end
struct Receiver_B1 end

# This routine ABSOLUTELY should not allocate any memory, called inside time loop.
@inbounds @fastmath function record!(it::Int64, issp::Int64, iss::Int64, pac::PFdtdc, pass::Vector{PFdtdss}, pap::PFdtdp, ::Receiver_B1)
	p=pap.p
	rinterpolatew=pass[issp].rinterpolatew
	irx1=pass[issp].irx1
	irx2=pass[issp].irx2
	irz1=pass[issp].irz1
	irz2=pass[issp].irz2

	for ipw in pac.activepw
		pw=p[ipw]
		recs=pass[issp].records[ipw]
		for (ifieldr, ifield) in enumerate(pac.irfields)
			@simd for ir = 1:pac.ageom[ipw][iss].nr
				recs[it,ir,ifieldr]=
				(
				pw[irz1[ipw][ir],irx1[ipw][ir],ifield]*
				rinterpolatew[ipw][1,ir]+
				pw[irz1[ipw][ir],irx2[ipw][ir],ifield]*
				rinterpolatew[ipw][2,ir]+
				pw[irz2[ipw][ir],irx1[ipw][ir],ifield]*
				rinterpolatew[ipw][3,ir]+
				pw[irz2[ipw][ir],irx2[ipw][ir],ifield]*
				rinterpolatew[ipw][4,ir]
				)
		end
	end
	end
end

# This routine ABSOLUTELY should not allocate any memory, called inside time loop.
@inbounds @fastmath function record!(it::Int64, issp::Int64, iss::Int64, pac::PFdtdc, pass::Vector{PFdtdss}, pap::PFdtdp, ::Receiver_B0)
	p=pap.p
	rinterpolatew=pass[issp].rinterpolatew
	irx1=pass[issp].irx1
	irz1=pass[issp].irz1

	for ipw in pac.activepw
		pw=p[ipw]
		recs=pass[issp].records[ipw]
		for (ifieldr, ifield) in enumerate(pac.irfields)
			@simd for ir = 1:pac.ageom[ipw][iss].nr
				recs[it,ir,ifieldr]=(pw[irz1[ipw][ir],irx1[ipw][ir],ifield])
		end
	end
	end
end


text:
struct Atom
    # an Atom is a repeatable unit and thus not unique in a lattice
    number::Integer
end

Base.isless(x::Atom, y::Atom) = x.number < y.number

struct LatticeDescription
    primitiveVectors::Vector{Vector{AbstractFloat}}
    basisAtoms::Dict{Atom,Vector{AbstractFloat}}
end

struct Particle
    # a Particle is a unique position in the lattice
    atom::Atom
    UC::Tuple
end

function generateLatticeShape(L::LatticeDescription, shape)
    d = length(L.primitiveVectors[1]) # dimension

    inside = zeros(d)

    res = Dict{Tuple,Vector{Atom}}()
    # you need to add the (0,0) still!
    n = 0
    # ring = collect(Iterators.flatten([x for x in Iterators.product([(-1, 0, 1) for _ in 1:d]...)]))

    inside_shape = true
    while inside_shape
        ring = [x for x in Iterators.product([(-n:n) for _ in 1:d]...) if (n in x || -n in x)]
        inside_shape = false
        for i ∈ ring
            pos = inside + sum(i .* L.primitiveVectors)
            uc_atoms = []
            for (at, at_pos) in L.basisAtoms
                if shape(at_pos .+ pos)
                    inside_shape = true
                    push!(uc_atoms, at)
                end
            end
            if !isempty(uc_atoms)
                res[i] = uc_atoms
            end
        end
        n += 1
    end
    res
end

function generateLatticeShape(L::LatticeDescription, shape, region)
    d = length(L.primitiveVectors[1]) # dimension

    inside = zeros(d)

    res = Dict{Tuple,Vector{Atom}}()
    n = 0

    # region has to be: [[x0,x1],[y0,y1]...]

    idx = collect(Iterators.product([x[1]:x[2] for x in region]...))
    for i in idx
        pos = inside + sum(i .* L.primitiveVectors)
        uc_atoms = []
        for (at, at_pos) in L.basisAtoms
            if shape(at_pos .+ pos)
                push!(uc_atoms, at)
            end
        end
        if !isempty(uc_atoms)
            res[i] = uc_atoms
        end

        n += 1
    end
    res
end

function absPositionsParticle(p::Particle, l::LatticeDescription)
    pos = sum(p.UC .* l.primitiveVectors)
    at = p.atom
    l.basisAtoms[at] .+ pos
end

function positionsLattice(unit_cells::Dict{Tuple,Vector{Atom}}, L::LatticeDescription)
    # takes a dictionary of unit cells and converts it to explicit positions
    # returns a vector of positions per atom in the system and mapping_index, which is useful if you need to know in which unit cell the atoms are, which is an unique mapping
    positions = Dict{Atom,Vector}()
    mapping = mappingAtoms(unit_cells)
    mapping_index = Dict{Atom,Vector}()
    for (partic, ind) in mapping

        pos = sum(partic.UC .* L.primitiveVectors)
        at = partic.atom
        absolute_at_pos = L.basisAtoms[at] .+ pos
        if haskey(positions, at)
            push!(positions[at], absolute_at_pos)
            push!(mapping_index[at], ind)
        else
            positions[at] =  [absolute_at_pos]
            mapping_index[at] =  [ind]
        end

    end
    positions, mapping_index
end

function mappingAtoms(unit_cells::Dict{Tuple,Vector{Atom}})
    # downfold the unit_cell into a single interaction matrix, and returns the mapping. We dont care so much about the specific ordering, as long as it is consistent
    mapping = Dict{Particle,Integer}()
    i = 1
    for (ind, uc_atoms) in unit_cells
        for at in uc_atoms
            mapping[Particle(at, ind)] = i
            i += 1
        end
    end
    mapping
end

countAtoms(unit_cells::Dict{Tuple,Vector{Atom}}) = length(mappingAtoms(unit_cells))

function coupleNN(unit_cells::Dict{Tuple,Vector{Atom}}, coupling)
    # coupling should be f(at1, at2, direction)
    d = length(first(keys(unit_cells)))
    nn_list = [x for x in Iterators.product([(-1:1) for _ in 1:d]...) if (1 in x || -1 in x || all(y -> y == 0., x))]

    mapping = mappingAtoms(unit_cells)
    N = length(mapping)
    res = zeros(N, N)

    for (ind, uc_atoms) in unit_cells
        for at1 in uc_atoms
            index1 = mapping[Particle(at1, ind)]
            for nn_pos in nn_list
                if (ind .+ nn_pos) in keys(unit_cells)
                    for at2 in unit_cells[ind .+ nn_pos]
                        index2 = mapping[Particle(at2, ind .+ nn_pos)]
                        res[index1, index2] = coupling(at1, at2, nn_pos)
                    end
                end
            end
        end
    end
    res
end

function coupleInternal(unit_cells::Dict{Tuple,Vector{Atom}}, coupling)
    # coupling should be f(at1, at2)
    mapping = mappingAtoms(unit_cells)
    N = length(mapping)
    res = zeros(N, N)

    for (ind, uc_atoms) in unit_cells
        for at1 in uc_atoms
            index1 = mapping[Particle(at1, ind)]
            for at2 in unit_cells[ind]
                if at1 != at2
                    index2 = mapping[Particle(at2, ind)]
                    res[index1, index2] = coupling(at1, at2)
                end
            end
        end
    end
    res
end

max_stars_repo_name: pgunnink/LLGTools

max_issues_repo_name: pgunnink/LLGTools

max_forks_repo_name: pgunnink/LLGTools


text:
# The action type of ToNextML policy depends on the model given
struct ToNextML{P<:Union{VDPTagProblem, DiscreteVDPTagProblem}, RNG<:AbstractRNG} <: Policy
    p::P
    rng::RNG
end

ToNextML(p::Union{VDPTagProblem, DiscreteVDPTagProblem}; rng=Random.GLOBAL_RNG) = ToNextML(p, rng)

function POMDPs.action(p::ToNextML{P}, s::TagState) where P <: VDPTagMDP
    next = next_ml_target(p.p, s.target)
    diff = next-s.agent
    return atan(diff[2], diff[1])
end

function POMDPs.action(p::ToNextML{P}, s::TagState) where P <: Union{VDPTagPOMDP, DiscreteVDPTagProblem}
    next = next_ml_target(mdp(p.p)::VDPTagMDP, s.target)
    diff = next-s.agent
    return TagAction(false, atan(diff[2], diff[1]))
end

POMDPs.action(p::ToNextML{P}, b::AbstractParticleBelief) where P <: Union{VDPTagPOMDP, DiscreteVDPTagProblem} = action(p, mode(b))
POMDPs.action(p::ToNextML{P}, b::Any) where P <: Union{VDPTagPOMDP, DiscreteVDPTagProblem} = action(p, rand(b))

struct ToNextMLSolver <: Solver
    rng::AbstractRNG
end

POMDPs.solve(s::ToNextMLSolver, p::Union{VDPTagProblem, DiscreteVDPTagProblem}) = ToNextML(p, s.rng)

# A POMDP policy which always takes a TagAction
struct ManageUncertainty <: Policy
    p::Union{VDPTagPOMDP, DiscreteVDPTagProblem}
    max_norm_std::Float64
end

function POMDPs.action(p::ManageUncertainty, b::AbstractParticleBelief)
    agent = first(particles(b)).agent
    prob_dict = ParticleFilters.probdict(b)
    if length(prob_dict) == 1
        return action(p, first(keys(prob_dict)))
    end
    target_particles = Array{Float64}(undef, 2, length(prob_dict))
    for (i, s) in enumerate(keys(prob_dict))
        target_particles[:,i] = s.target
    end
    normal_dist = fit(MvNormal, target_particles, collect(values(prob_dict))) # particles should be unique
    angle = POMDPs.action(ToNextML(mdp(p.p)), TagState(agent, mean(normal_dist)))
    return TagAction(sqrt(det(cov(normal_dist))) > p.max_norm_std, angle)
end

POMDPs.action(p::ManageUncertainty, b::Any) = TagAction(false, action(ToNextML(mdp(p.p)), rand(b)))
POMDPs.action(p::ManageUncertainty, s::TagState) = TagAction(false, action(ToNextML(mdp(p.p)), s))

mutable struct NextMLFirst{RNG<:AbstractRNG}
    rng::RNG
end

function next_action(gen::NextMLFirst, mdp::VDPTagMDP, s::TagState, snode)
    if n_children(snode) < 1
        return POMDPs.action(ToNextML(mdp, gen.rng), s)::Float64
    else
        return 2*pi*rand(gen.rng)
    end
end

next_action(gen::NextMLFirst, pomdp::Union{VDPTagPOMDP, DiscreteVDPTagProblem}, b, onode) = TagAction(false, next_action(gen, mdp(pomdp), rand(gen.rng, b), onode))


text:
# module for drawing goal flag

using Plots
pyplot()

mutable struct Goal
  pose
  radius
  value

  # init
  function Goal(x, y; radius=0.3, value=0.0)
    self = new()
    self.pose = [x, y]
    self.radius = radius
    self.value = value
    return self
  end
end

function inside(self::Goal, pose)
  if self.radius > sqrt((self.pose[1]-pose[1])^2 + (self.pose[2]-pose[2])^2)
    return true
  else
    return false
  end
end

function draw!(self::Goal)
  scatter!([self.pose[1]+0.1], [self.pose[2]+0.5],
           markershape=:rtriangle,
           markercolor=:red,
           markersize=9)
  plot!([self.pose[1], self.pose[1]], [self.pose[2], self.pose[2]+0.6],
        color=:black)
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
include("Transfer.jl")
include("Solution.jl")
include("statistics.jl")


text:
@recipe(RoundedRect) do scene
    Theme(
        rect = BBox(0, 100, 0, 100),
        cornerradius = 5,
        cornersegments = 10,
        raw = true,
        color = RGBf(0.9, 0.9, 0.9),
        strokecolor = RGBf(0, 0, 0)
    )
end

function plot!(roundrect::RoundedRect)
    @extract(roundrect, (
        rect, cornerradius, cornersegments, raw, color, strokecolor
    ))

    heightattr = roundrect.height
    widthattr = roundrect.width

    roundedrectpoints = lift(rect, cornerradius, cornersegments) do rect,
            cr, csegs

        cr = min(width(rect) / 2, height(rect) / 2, cr)

        # inner corners
        ictl = topleft(rect) .+ Point2(cr, -cr)
        ictr = topright(rect) .+ Point2(-cr, -cr)
        icbl = bottomleft(rect) .+ Point2(cr, cr)
        icbr = bottomright(rect) .+ Point2(-cr, cr)

        cstr = anglepoint.(Ref(ictr), LinRange(0, pi/2, csegs), cr)
        cstl = anglepoint.(Ref(ictl), LinRange(pi/2, pi, csegs), cr)
        csbl = anglepoint.(Ref(icbl), LinRange(pi, 3pi/2, csegs), cr)
        csbr = anglepoint.(Ref(icbr), LinRange(3pi/2, 2pi, csegs), cr)

        arr = [cstr; cstl; csbl; csbr]
    end

    poly!(roundrect, roundedrectpoints, raw = raw, color = color, strokecolor = strokecolor)
end


text:
## Grouping machinery

fast_hashed(v::AbstractVector) = isbitstype(eltype(v)) ? PooledArray(v) : v

function indices_iterator(cols)
    isempty(cols) && return (:,)
    grouping_sa = StructArray(map(refarray∘fast_hashed, cols))
    gp = GroupPerm(grouping_sa)
    return (sortperm(gp)[rg] for rg in gp)
end

function validprimaryindices(v, idxs, c::CartesianIndex)
    i = idxs === (:) ? firstindex(v, 1) : clamp(first(idxs), axes(v, 1))
    t = validtailindices(v, c)
    return (i, t...)
end

validtailindices(v, c::CartesianIndex) =
    Tuple(Broadcast.newindex(CartesianIndices(tail(axes(v))), c))

validindices(v, c::CartesianIndex) = Tuple(Broadcast.newindex(v, c))

function subgroup(e::Entry, idxs, c::CartesianIndex)
    primary = map(v -> view(v, validprimaryindices(v, idxs, c)...), e.primary)
    positional = map(v -> view(v, idxs, validtailindices(v, c)...), e.positional)
    named = map(v -> view(v, idxs, validtailindices(v, c)...), e.named)
    labels = copy(e.labels)
    map!(values(labels)) do l
        l′ = Broadcast.broadcastable(l)
        return l′[validindices(l′, c)...]
    end
    return Entry(e; primary, positional, named, labels)
end

allvariables(e::Entry) = (e.primary..., e.positional..., e.named...)
allvariables(l::Layer) = (l.positional..., l.named...)

function shape(x::Union{Entry, Layer})
    arrays = map(var -> var isa ArrayLike ? var : fill(nothing), allvariables(x))
    return Broadcast.combine_axes(arrays...)
end

splitapply(entry::Entry) = splitapply(identity, entry)

function splitapply(f, entry::Entry)
    axs = shape(entry)
    grouping = filter(v -> v isa AbstractVector, Tuple(entry.primary))
    entries = Entry[]
    foreach(indices_iterator(grouping)) do idxs
        for c in CartesianIndices(tail(axs))
            # TODO: for analyses returning several entries, rearrange in correct order.
            res = f(subgroup(entry, idxs, c))
            res isa Entry ? push!(entries, res) : append!(entries, res)
        end
    end
    return entries
end

## Transform `Layers` into list of `Entry`

assert_equal(a, b) = (@assert(a == b); a)

function unnest(arr::AbstractArray{<:AbstractArray})
    inner_size = mapreduce(size, assert_equal, arr)
    outer_size = size(arr)
    flattened = reduce(vcat, map(vec, vec(arr)))
    return reshape(flattened, inner_size..., outer_size...)
end

unnest(arr::NTuple{<:Any, <:AbstractArray}) = unnest(collect(arr))

"""
    to_entry(layer::Layer)

Convert `layer` to equivalent entry, excluding transformations.
"""
function to_entry(layer::Layer)
    labels = Dict{KeyType, Any}()
    primary_pairs, positional_list, named_pairs = [], [], []
    for c in (layer.positional, layer.named)
        for (key, val) in pairs(c)
            l, arr = getlabeledarray(layer, val)
            labels[key] = l
            if key isa Int
                push!(positional_list, arr)
            elseif !iscontinuous(arr)
                push!(primary_pairs, key => arr)
            else
                push!(named_pairs, key => arr)
            end
        end
    end
    primary = NamedTuple(primary_pairs)
    positional = Tuple(positional_list)
    named = NamedTuple(named_pairs)
    return Entry(; primary, positional, named, labels)
end

function process(layer::Layer)
    init = to_entry(layer)
    res = foldl(process, layer.transformations; init)
    return res isa Entry ? splitapply(res) : res
end

process(v::AbstractArray{Entry}, f) = map(f, v)

process(le::Entry, f) = f(le)


text:
using Pidoh
using Random

@testset "Flip bit string" begin
    @test flip(true) == false
    @test flip(false) == true
end

@testset "OneMax Fitness Function" begin
    @test fitness(BitArray([true, true, true]), OneMax(3)) == 3
    @test fitness(BitArray([false, false, false, false]), OneMax(4)) == 0
    @test fitness(BitArray([]), OneMax(0)) == 0
end

@testset "ZeroMax Fitness Function" begin
    @test fitness(BitArray([true, true, true]), ZeroMax(3)) == 0
    @test fitness(BitArray([false, false, false, false]), ZeroMax(4)) == 4
    @test fitness(BitArray([]), ZeroMax(0)) == 0
end

jump = Jump(10, 3)
@testset "ZeroMax Fitness Function" begin
    @test fitness(BitArray([0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), jump) == jump.jumpsize
    @test fitness(BitArray([1, 0, 0, 0, 0, 0, 0, 0, 0, 0]), jump) == jump.jumpsize + 1
    @test fitness(BitArray([1, 1, 1, 1, 1, 1, 1, 1, 1, 1]), jump) == jump.jumpsize + 10
    @test fitness(BitArray([1, 1, 1, 1, 1, 1, 1, 1, 1, 0]), jump) == 1
    @test fitness(BitArray([1, 1, 1, 1, 1, 1, 1, 0, 0, 0]), jump) == 10
    @test fitness(BitArray([1, 1, 1, 1, 1, 1, 1, 1, 0, 0]), jump) == 2

end

ins1 = Instance(BitArray([true, true, false]), OneMax(3))

@testset "UniformlyIndependentMutation mutation" begin
    @test mutation(ins1, UniformlyIndependentMutation(0.0)).individual == ins1.individual
    @test mutation(ins1, UniformlyIndependentMutation(1.0)).individual != ins1.individual
    @test fitness(mutation(ins1, UniformlyIndependentMutation(1.0))) == 1
end

max_stars_repo_name: DTUComputeTONIA/SDRLS