text:
struct Elastic <: Equation end
@declare_dofs Elastic [:stxx, :styy, :stxy, :v, :w]

struct Elastic_waves <: Scenario end
struct GaussianPoint <: Scenario end
struct Elastic_source <: Scenario end
struct StiffInclusion <: Scenario end

function is_periodic_boundary(equation::Elastic, scenario::Elastic_waves)
    # Elastic_waves scenario with periodic boundary conditions
    true
end

function is_periodic_boundary(equation::Elastic, scenario::GaussianPoint)
    # Elastic_waves scenario with periodic boundary conditions
    true
end

function is_periodic_boundary(equation::Elastic, scenario::Elastic_source)
    # Elastic_waves scenario with periodic boundary conditions
    true
end

function is_periodic_boundary(equation::Elastic, scenario::StiffInclusion)
    # Elastic_waves scenario with periodic boundary conditions
    true
end
"""
# Refelcting Boundary conditions
function evaluate_boundary(eq::Elastic, scenario::Elastic_waves, face, normalidx, dofsface, dofsfaceneigh)
# dofsface and dofsfaceneigh have shape (num_2d_quadpoints, dofs)
# you need to set dofsfaceneigh
    dofsfaceneigh .= dofsface
    if normalidx == 1
        dofsfaceneigh[:, 1] = -dofsface[:, 1]
    else
        dofsfaceneigh[:, 2] = -dofsface[:, 2]
    end
end
"""

function get_initial_values(eq::Elastic, scenario::GaussianPoint, global_position; t=0.0)
    x, y = global_position
    if is_analytical_solution(eq, scenario)
        println("No analytical solution implemented")
    else
        # From Equation (22.67) from book of Leveque
        # TODO need to look fot actual values
        # some constants (Example 22.2 from book of Leveque only left side material)
        λ = 2
        μ = 1
        ρ = 1

        #if (0.2 < x < 0.35) # As we have only domain fom (0 to 1)
        if (x <= 0.5 || y <= 0.5) || (x >= 0.5 && y >= 0.5)
            #[λ+2*μ , λ , 0 , sqrt(5), 0]
            [0, 0, 0, exp(-1000 * (x - 0.5)^2 - 1000 *(y - 0.5)^2) + 1 , exp(-1000 * (x - 0.5)^2 - 1000 *(y - 0.5)^2) + 1 ]
        else
            [0 , 0 , 0 , 1 , 1]
        end

    end
end

function get_initial_values(eq::Elastic, scenario::Elastic_waves, global_position; t=0.0)
    x, y = global_position
    if is_analytical_solution(eq, scenario)
        println("No analytical solution implemented")
    else
        # From Equation (22.67) from book of Leveque
        # TODO need to look fot actual values
        # some constants (Example 22.2 from book of Leveque only left side material)
        λ = 2
        μ = 1
        ρ = 1

        if (0.2 < x < 0.35) # As we have only domain fom (0 to 1)
            [λ+2*μ , λ , 0 , sqrt(5), 0]
        else
            [0 , 0 , 0 , 0 , 0]
        end

    end
end

<<<<<<< Updated upstream
=======
function get_initial_values(eq::Elastic, scenario::StiffInclusion, global_position; t=0.0)
    x, y = global_position
    if is_analytical_solution(eq, scenario)
        println("No analytical solution implemented")
    else
        # From Equation (22.67) from book of Leveque
        # TODO need to look fot actual values
        # some constants (Example 22.2 from book of Leveque only left side material)
        λ = 4
        μ = 0.5
        ρ = 1

        if (0.6 < x < 0.65) # As we have only domain fom (0 to 1)
            [λ+2*μ , λ , 0 , sqrt(5), 0]
        else
            [0 , 0 , 0 , 0 , 0]
        end

    end
end

>>>>>>> Stashed changes
function get_initial_values(eq::Elastic, scenario::Elastic_source, global_position; t=0.0)
    x, y = global_position
    if is_analytical_solution(eq, scenario)
        #println("No analytical solution implemented")
        [( 1 * sin(6*pi*x/100) * sin(6*pi*y/100) * sin(6*pi*t/100) ) , ( 2 * sin(2*pi*x/50) * sin(2*pi*y/50) * sin(2*pi*t/50) ), ( 3 * sin(2*pi*x/100) * sin(2*pi*y/100) * sin(2*pi*t/100) ), ( 4 * sin(2*pi*x/50) * sin(2*pi*y/50) * sin(2*pi*t/10) ), ( 5 * sin(6*pi*x/100) * sin(6*pi*y/100) * sin(2*pi*t/5) ) ]
    else
        [0,0,0,0,0]
    end

end

<<<<<<< Updated upstream
function get_initial_values(eq::Elastic, scenario::StiffInclusion, global_position; t=0.0)
    x, y = global_position
    if is_analytical_solution(eq, scenario)
        println("No analytical solution implemented")
    else
        # From Equation (22.67) from book of Leveque
        # TODO need to look fot actual values
        # some constants (Example 22.2 from book of Leveque only left side material)
        λ = 2
        μ = 1
        ρ = 1

        if (0.6 < x < 0.65) # As we have only domain fom (0 to 1)
            [λ+2*μ , λ , 0 , 2, 0]
        else
            [0 , 0 , 0 , 0 , 0]
        end
    end
end

=======
>>>>>>> Stashed changes

function is_analytical_solution(equation::Elastic, scenario::Elastic_waves)
    false
end

function is_analytical_solution(equation::Elastic, scenario::GaussianPoint)
    false
end

function is_analytical_solution(equation::Elastic, scenario::Elastic_source)
    true
end

function is_analytical_solution(equation::Elastic, scenario::StiffInclusion)
    false
end

# See equation (22.39) from book of Leveque for clarification
function evaluate_flux(eq::Elastic, celldofs, cellflux)
    # size(cellflux)=(2*order*order, ndof)=(2,5)
    # size(celldofs)=(order*order, ndof)=(1,5)

    order_sq, ndof = size(celldofs)
    for j=1:order_sq
        # explicit names of dofs for clarity
        stxx = celldofs[j, 1]
        styy = celldofs[j, 2]
        stxy = celldofs[j, 3]
        v = celldofs[j, 4]
        w = celldofs[j, 5]

        #TODO need to look fot actual values
        #some constants (Example 22.2 from book of Leveque only left side material)
        λ = 2
        μ = 1
        ρ = 1

        # here flux is computed elementwise
        cellflux[j,1]= (-(λ+2*μ)) * v
        cellflux[j + order_sq, 1] = (-λ) * w

        cellflux[j,2] = (-λ) * v
        cellflux[j + order_sq, 2] = (-(λ+2*μ)) * w

        cellflux[j,3] = (-μ) * w
        cellflux[j + order_sq, 3] = (-μ) * v

        cellflux[j,4] = (-1/ρ) * stxx
        cellflux[j + order_sq, 4] = (-1/ρ) * stxy

        cellflux[j,5] = (-1/ρ) * stxy
        cellflux[j + order_sq, 5] = (-1/ρ) * styy

    end

end

function evaluate_flux_with_stiffinclusion(eq::Elastic, scenario::StiffInclusion, celldofs, cellflux, global_position)
    # size(cellflux)=(2*order*order, ndof)=(2,5)
    # size(celldofs)=(order*order, ndof)=(1,5)
    x,y = global_position
    order_sq, ndof = size(celldofs)
    for j=1:order_sq
        # explicit names of dofs for clarity
        stxx = celldofs[j, 1]
        styy = celldofs[j, 2]
        stxy = celldofs[j, 3]
        v = celldofs[j, 4]
        w = celldofs[j, 5]

        if (x <=0.5 && y <= 0.5) ||  ( x<=0.5 && y <= ((-0.5*x) + 0.75) )
            λ = 4
            μ = 0.5
            ρ = 1
        else
            λ = 2
            μ = 1
            ρ = 1
        end

        # here flux is computed elementwise
        cellflux[j,1]= (-(λ+2*μ)) * v
        cellflux[j + order_sq, 1] = (-λ) * w

        cellflux[j,2] = (-λ) * v
        cellflux[j + order_sq, 2] = (-(λ+2*μ)) * w

        cellflux[j,3] = (-μ) * w
        cellflux[j + order_sq, 3] = (-μ) * v

        cellflux[j,4] = (-1/ρ) * stxx
        cellflux[j + order_sq, 4] = (-1/ρ) * stxy

        cellflux[j,5] = (-1/ρ) * stxy
        cellflux[j + order_sq, 5] = (-1/ρ) * styy

    end

end

# ATTEMPT-1
"""
# Evaluation of flux with Source term
# See equation (57-63) from paper of Martin Kaser for clarification
function evaluate_flux_with_source(eq::Elastic, celldofs, cellflux, cell, basis)

    order_sq, ndof = size(celldofs)

    #TODO need to look fot actual values
    #some constants (Example 22.2 from book of Leveque only left side material)
    λ = 4
    μ = 0.5
    ρ = 1

    up0 = [1,2,3,4,5]
    lamdapx = [42.9,50,100,50,42.9]
    lamdapy = [42.9,50,100,50,42.9]
    timep = [42.9,50,100,10,5]

    upsol = [
            (x,y,t) -> [up0[1] * sin((2*pi*x) / lamdapx[1]) * sin((2*pi*y) / lamdapy[1]) * sin((2*pi*t) / timep[1])],
            (x,y,t) -> [up0[2] * sin((2*pi*x) / lamdapx[2]) * sin((2*pi*y) / lamdapy[2]) * sin((2*pi*t) / timep[2])],
            (x,y,t) -> [up0[3] * sin((2*pi*x) / lamdapx[3]) * sin((2*pi*y) / lamdapy[3]) * sin((2*pi*t) / timep[3])],
            (x,y,t) -> [up0[4] * sin((2*pi*x) / lamdapx[4]) * sin((2*pi*y) / lamdapy[4]) * sin((2*pi*t) / timep[4])],
            (x,y,t) -> [up0[5] * sin((2*pi*x) / lamdapx[5]) * sin((2*pi*y) / lamdapy[5]) * sin((2*pi*t) / timep[5])]
            ]

    dudtp = [
            (x,y,t) -> [((2*pi) * timep[1]) * up0[1] * sin((2*pi*x) / lamdapx[1]) * sin((2*pi*y) / lamdapy[1]) * cos((2*pi*t) / timep[1])],
            (x,y,t) -> [((2*pi) * timep[2]) * up0[2] * sin((2*pi*x) / lamdapx[2]) * sin((2*pi*y) / lamdapy[2]) * cos((2*pi*t) / timep[2])],
            (x,y,t) -> [((2*pi) * timep[3]) * up0[3] * sin((2*pi*x) / lamdapx[3]) * sin((2*pi*y) / lamdapy[3]) * cos((2*pi*t) / timep[3])],
            (x,y,t) -> [((2*pi) * timep[4]) * up0[4] * sin((2*pi*x) / lamdapx[4]) * sin((2*pi*y) / lamdapy[4]) * cos((2*pi*t) / timep[4])],
            (x,y,t) -> [((2*pi) * timep[5]) * up0[5] * sin((2*pi*x) / lamdapx[5]) * sin((2*pi*y) / lamdapy[5]) * cos((2*pi*t) / timep[5])]
            ]

    dudxp = [
            (x,y,t) -> [((2*pi) / lamdapx[1]) * up0[1] * cos((2*pi*x) / lamdapx[1]) * sin((2*pi*y) / lamdapy[1]) * sin((2*pi*t) / timep[1])],
            (x,y,t) -> [((2*pi) / lamdapx[2]) * up0[2] * cos((2*pi*x) / lamdapx[2]) * sin((2*pi*y) / lamdapy[2]) * sin((2*pi*t) / timep[2])],
            (x,y,t) -> [((2*pi) / lamdapx[3]) * up0[3] * cos((2*pi*x) / lamdapx[3]) * sin((2*pi*y) / lamdapy[3]) * sin((2*pi*t) / timep[3])],
            (x,y,t) -> [((2*pi) / lamdapx[4]) * up0[4] * cos((2*pi*x) / lamdapx[4]) * sin((2*pi*y) / lamdapy[4]) * sin((2*pi*t) / timep[4])],
            (x,y,t) -> [((2*pi) / lamdapx[5]) * up0[5] * cos((2*pi*x) / lamdapx[5]) * sin((2*pi*y) / lamdapy[5]) * sin((2*pi*t) / timep[5])]
            ]

    dudyp = [
            (x,y,t) -> [((2*pi) / lamdapy[1]) * up0[1] * sin((2*pi*x) / lamdapx[1]) * cos((2*pi*y) / lamdapy[1]) * sin((2*pi*t) / timep[1])],
            (x,y,t) -> [((2*pi) / lamdapy[2]) * up0[2] * sin((2*pi*x) / lamdapx[2]) * cos((2*pi*y) / lamdapy[2]) * sin((2*pi*t) / timep[2])],
            (x,y,t) -> [((2*pi) / lamdapy[3]) * up0[3] * sin((2*pi*x) / lamdapx[3]) * cos((2*pi*y) / lamdapy[3]) * sin((2*pi*t) / timep[3])],
            (x,y,t) -> [((2*pi) / lamdapy[4]) * up0[4] * sin((2*pi*x) / lamdapx[4]) * cos((2*pi*y) / lamdapy[4]) * sin((2*pi*t) / timep[4])],
            (x,y,t) -> [((2*pi) / lamdapy[5]) * up0[5] * sin((2*pi*x) / lamdapx[5]) * cos((2*pi*y) / lamdapy[5]) * sin((2*pi*t) / timep[5])]
            ]

    #Hence source term is according to equation (2) from same paper size of soure term is: size is (5,1)
    sp = [
        (x,y,z) -> [dudtp[1](x,y,z) + (-(λ+2*μ)) * dudxp[1](x,y,z) + (-λ) * dudyp[1](x,y,z)],
        (x,y,z) -> [dudtp[2](x,y,z) + (-λ) * dudxp[2](x,y,z) + (-(λ+2*μ)) * dudyp[2](x,y,z)],
        (x,y,z) -> [dudtp[3](x,y,z) + (-μ) * dudxp[3](x,y,z) + (-μ) * dudyp[3](x,y,z)],
        (x,y,z) -> [dudtp[4](x,y,z) + (-1/ρ) * dudxp[4](x,y,z) + (-1/ρ)  * dudyp[4](x,y,z)],
        (x,y,z) -> [dudtp[5](x,y,z) + (-1/ρ) * dudxp[5](x,y,z) + (-1/ρ)  * dudyp[5](x,y,z)]
        ]

    #Updating Flux including source terms
    linearindex = vec(collect(Iterators.product(1:basis.order, 1:basis.order))) #linear indexing

    for j=1:order_sq
        # explicit names of dofs for clarity
        stxx = celldofs[j, 1]
        styy = celldofs[j, 2]
        stxy = celldofs[j, 3]
        v = celldofs[j, 4]
        w = celldofs[j, 5]

        # Finding global Coordinates of Cell
        localx,localy = globalposition(cell, (basis.quadpoints[linearindex[j][1]], basis.quadpoints[linearindex[j][2]]))
        #@info localx
        #@info localy
        localz = 0.0 # as problem is 2 Dimensional

        # here flux is computed elementwise, one source term is divided into 2 parts equally, one along X and other along Y
        cellflux[j,1]= ( (-(λ+2*μ)) * v ) - ( ((sp[1](localx,localy,localz))[1][1])/2 )
        cellflux[j + order_sq, 1] = ( (-λ) * w )- ( ((sp[1](localx,localy,localz))[1][1])/2 )

        cellflux[j,2] = ( (-λ) * v ) - ( ((sp[2](localx,localy,localz))[1][1])/2 )
        cellflux[j + order_sq, 2] = ( (-(λ+2*μ)) * w ) - ( ((sp[2](localx,localy,localz))[1][1])/2 )

        cellflux[j,3] = ( (-μ) * w ) - ( ((sp[3](localx,localy,localz))[1][1])/2 )
        cellflux[j + order_sq, 3] = ( (-μ) * v ) - ( ((sp[3](localx,localy,localz))[1][1])/2 )

        cellflux[j,4] = ( (-1/ρ) * stxx ) - ( ((sp[4](localx,localy,localz))[1][1])/2 )
        cellflux[j + order_sq, 4] = ( (-1/ρ) * stxy ) - ( ((sp[4](localx,localy,localz))[1][1])/2 )

        cellflux[j,5] = ( (-1/ρ) * stxy ) - ( ((sp[5](localx,localy,localz))[1][1])/2 )
        cellflux[j + order_sq, 5] = ( (-1/ρ) * styy ) - ( ((sp[5](localx,localy,localz))[1][1])/2 )
    end

end
"""

# ATTEMPT 2
# Evaluation of Source term
# See equation (57-63) from paper of Martin Kaser for clarification
function evaluate_source(eq::Elastic, time, position)
<<<<<<< Updated upstream
    # Source term foe general problem (section 6.1 in same paper)
=======
    """
>>>>>>> Stashed changes
    λ = 2
    μ = 1
    ρ = 1

    up0 = [1,2,3,4,5]
    lamdapx = [42.9,50,100,50,42.9]
    lamdapy = [42.9,50,100,50,42.9]
    timep = [42.9,50,100,10,5]

    num_dofs = size(up0)
    upsol = zeros(num_dofs)
    dudtp = zeros(num_dofs)
    dudxp = zeros(num_dofs)
    dudyp = zeros(num_dofs)
    sp = zeros(num_dofs)

    upsol = [
            (x,y,t) -> [up0[1] * sin((2*pi*x) / lamdapx[1]) * sin((2*pi*y) / lamdapy[1]) * sin((2*pi*t) / timep[1])],
            (x,y,t) -> [up0[2] * sin((2*pi*x) / lamdapx[2]) * sin((2*pi*y) / lamdapy[2]) * sin((2*pi*t) / timep[2])],
            (x,y,t) -> [up0[3] * sin((2*pi*x) / lamdapx[3]) * sin((2*pi*y) / lamdapy[3]) * sin((2*pi*t) / timep[3])],
            (x,y,t) -> [up0[4] * sin((2*pi*x) / lamdapx[4]) * sin((2*pi*y) / lamdapy[4]) * sin((2*pi*t) / timep[4])],
            (x,y,t) -> [up0[5] * sin((2*pi*x) / lamdapx[5]) * sin((2*pi*y) / lamdapy[5]) * sin((2*pi*t) / timep[5])]
            ]
    for i in 1:
    dudtp = [
            (x,y,t) -> [((2*pi) * timep[1]) * up0[1] * sin((2*pi*x) / lamdapx[1]) * sin((2*pi*y) / lamdapy[1]) * cos((2*pi*t) / timep[1])],
            (x,y,t) -> [((2*pi) * timep[2]) * up0[2] * sin((2*pi*x) / lamdapx[2]) * sin((2*pi*y) / lamdapy[2]) * cos((2*pi*t) / timep[2])],
            (x,y,t) -> [((2*pi) * timep[3]) * up0[3] * sin((2*pi*x) / lamdapx[3]) * sin((2*pi*y) / lamdapy[3]) * cos((2*pi*t) / timep[3])],
            (x,y,t) -> [((2*pi) * timep[4]) * up0[4] * sin((2*pi*x) / lamdapx[4]) * sin((2*pi*y) / lamdapy[4]) * cos((2*pi*t) / timep[4])],
            (x,y,t) -> [((2*pi) * timep[5]) * up0[5] * sin((2*pi*x) / lamdapx[5]) * sin((2*pi*y) / lamdapy[5]) * cos((2*pi*t) / timep[5])]
            ]

    dudxp = [
            (x,y,t) -> [((2*pi) / lamdapx[1]) * up0[1] * cos((2*pi*x) / lamdapx[1]) * sin((2*pi*y) / lamdapy[1]) * sin((2*pi*t) / timep[1])],
            (x,y,t) -> [((2*pi) / lamdapx[2]) * up0[2] * cos((2*pi*x) / lamdapx[2]) * sin((2*pi*y) / lamdapy[2]) * sin((2*pi*t) / timep[2])],
            (x,y,t) -> [((2*pi) / lamdapx[3]) * up0[3] * cos((2*pi*x) / lamdapx[3]) * sin((2*pi*y) / lamdapy[3]) * sin((2*pi*t) / timep[3])],
            (x,y,t) -> [((2*pi) / lamdapx[4]) * up0[4] * cos((2*pi*x) / lamdapx[4]) * sin((2*pi*y) / lamdapy[4]) * sin((2*pi*t) / timep[4])],
            (x,y,t) -> [((2*pi) / lamdapx[5]) * up0[5] * cos((2*pi*x) / lamdapx[5]) * sin((2*pi*y) / lamdapy[5]) * sin((2*pi*t) / timep[5])]
            ]

    dudyp = [
            (x,y,t) -> [((2*pi) / lamdapy[1]) * up0[1] * sin((2*pi*x) / lamdapx[1]) * cos((2*pi*y) / lamdapy[1]) * sin((2*pi*t) / timep[1])],
            (x,y,t) -> [((2*pi) / lamdapy[2]) * up0[2] * sin((2*pi*x) / lamdapx[2]) * cos((2*pi*y) / lamdapy[2]) * sin((2*pi*t) / timep[2])],
            (x,y,t) -> [((2*pi) / lamdapy[3]) * up0[3] * sin((2*pi*x) / lamdapx[3]) * cos((2*pi*y) / lamdapy[3]) * sin((2*pi*t) / timep[3])],
            (x,y,t) -> [((2*pi) / lamdapy[4]) * up0[4] * sin((2*pi*x) / lamdapx[4]) * cos((2*pi*y) / lamdapy[4]) * sin((2*pi*t) / timep[4])],
            (x,y,t) -> [((2*pi) / lamdapy[5]) * up0[5] * sin((2*pi*x) / lamdapx[5]) * cos((2*pi*y) / lamdapy[5]) * sin((2*pi*t) / timep[5])]
            ]

    #Hence source term is according to equation (2) from same paper size of soure term is: size is (5,1)
    sp = [
        (x,y,t) -> [dudtp[1](x,y,t) + (-(λ+2*μ)) * dudxp[1](x,y,t) + (-λ) * dudyp[1](x,y,t)],
        (x,y,t) -> [dudtp[2](x,y,t) + (-λ) * dudxp[2](x,y,t) + (-(λ+2*μ)) * dudyp[2](x,y,t)],
        (x,y,t) -> [dudtp[3](x,y,t) + (-μ) * dudxp[3](x,y,t) + (-μ) * dudyp[3](x,y,t)],
        (x,y,t) -> [dudtp[4](x,y,t) + (-1/ρ) * dudxp[4](x,y,t) + (-1/ρ)  * dudyp[4](x,y,t)],
        (x,y,t) -> [dudtp[5](x,y,t) + (-1/ρ) * dudxp[5](x,y,t) + (-1/ρ)  * dudyp[5](x,y,t)]
        ]
    """

    # Source term from Lamb's problem:Ricker welvet (section 6.1 in same paper)
    rhosolid = 2200
    cp=3200
    cs=1847
    xs=[0.5,0.5] #Source application
    fc=14.5 #Ricker welvet central frequesncy
    a1=-2000
    a2=-(pi*pi*fc*fc)
    td=0.08 #source delay
    px,py = position
    y = [px-0.5, py-1]
    y_norm = y[1]^2 + y[2]^2
    dirac_delta = 5*(1/sqrt(pi))*exp(-25*(y_norm))

    g= a1 * (0.5 + (a2 * (time-td)^2 )) * (exp( a2 * (time -td)^2 ))

    sp = [(x,y,z) -> [0],
          (x,y,z) -> [0],
          (x,y,z) -> [0],
          (x,y,z) -> [0],
          (x,y,z) -> [(g/rhosolid)*dirac_delta]
        ]
<<<<<<< Updated upstream
    """
=======

>>>>>>> Stashed changes

    #@info "from function"
    #@info sp[1](1,2,3)
    #Returning sp
    sp

end

# See equation (22.17 and 22.18) from book of Leveque
function max_eigenval(eq::Elastic, celldata, normalidx)
    maxeigenval = 0.

    #TODO need to look fot actual values
    #some constants (Example 22.2 from book of Leveque only left side material)
    λ = 4
    μ = 0.5
    ρ = 1

    cp = sqrt((λ+2*μ)/ρ)
    cs = sqrt(μ/ρ)

    if (cp > maxeigenval && cp > cs)
        maxeigenval = cp
    elseif (cs > maxeigenval && cs > cp)
        maxeigenval = cs
    end
    maxeigenval
end


text:
using Cubature
using Distributions
using Gadfly
using PMTK4BasicModels

y1 = 90
n1 = 100
y2 = 2
n2 = 2
β1 = Beta(y1+1, n1-y1+1)
β2 = Beta(y2+1, n2-y2+1)

"""
Analysis of binomial differences.
"""
function unpairedTestHandednessDemo()
    θ1, θ2, δ, dens = contingencyTableUnpairedDiffPostMc(n1, y1, n2, y2)
    diff = θ1-θ2

    qq = quantile(EmpiricalUnivariateDistribution(diff), [.025, .975])

    pp = plot(x=δ, y=dens, Geom.line, Theme(default_color=colorant"black", line_width=3px),
              xintercept=qq, Geom.vline(color=colorant"blue", size=3px))
end

"""
Compute a probability whether buying from the first seller is better than from the second.
First value is given by a numerical integration method, the second - by MCMC.
"""
function amazonsSellerDelta()
    integrand(θ::Vector{Float64}) = pdf(β1, θ[1]) .* pdf(β2, θ[2]) .* (θ[1] > θ[2])
    pGreaterExact, err = hcubature(integrand, [0, 0], [1, 1], reltol=1e-4)

    θ1, θ2, δ, dens = contingencyTableUnpairedDiffPostMc(n1, y1, n2, y2)
    pGreaterMC = mean(θ1.>θ2)

    [pGreaterExact, pGreaterMC]
end

"""
Exact posteriors p(θ|D)
"""
function amazonSellersTheta()
    xs = linspace(0.001, 0.999, 40)
    ps1 = pdf(β1, xs)
    ps2 = pdf(β2, xs)

    pdf1_l = layer(x=xs, y=ps1, Geom.line, Theme(default_color=colorant"red", line_width=3px))
    pdf2_l = layer(x=xs, y=ps2, Geom.point, Theme(default_color=colorant"green", line_width=3px))

    pp = plot(pdf1_l, pdf2_l)
end

ext: jl


text:
@testset "Bitmask flags" begin
    test(BitmaskDefinition, bitmask_by_name, :VkQueryPipelineStatisticFlagBits, :(
        @bitmask_flag QueryPipelineStatisticFlag::UInt32 begin
            QUERY_PIPELINE_STATISTIC_INPUT_ASSEMBLY_VERTICES_BIT = 1
            QUERY_PIPELINE_STATISTIC_INPUT_ASSEMBLY_PRIMITIVES_BIT = 2
            QUERY_PIPELINE_STATISTIC_VERTEX_SHADER_INVOCATIONS_BIT = 4
            QUERY_PIPELINE_STATISTIC_GEOMETRY_SHADER_INVOCATIONS_BIT = 8
            QUERY_PIPELINE_STATISTIC_GEOMETRY_SHADER_PRIMITIVES_BIT = 16
            QUERY_PIPELINE_STATISTIC_CLIPPING_INVOCATIONS_BIT = 32
            QUERY_PIPELINE_STATISTIC_CLIPPING_PRIMITIVES_BIT = 64
            QUERY_PIPELINE_STATISTIC_FRAGMENT_SHADER_INVOCATIONS_BIT = 128
            QUERY_PIPELINE_STATISTIC_TESSELLATION_CONTROL_SHADER_PATCHES_BIT = 256
            QUERY_PIPELINE_STATISTIC_TESSELLATION_EVALUATION_SHADER_INVOCATIONS_BIT = 512
            QUERY_PIPELINE_STATISTIC_COMPUTE_SHADER_INVOCATIONS_BIT = 1024
        end
    ))

    test(BitmaskDefinition, bitmask_by_name, :VkSparseMemoryBindFlagBits, :(
        @bitmask_flag SparseMemoryBindFlag::UInt32 begin
            SPARSE_MEMORY_BIND_METADATA_BIT = 1
        end
    ))

    test(BitmaskDefinition, bitmask_by_name, :VkShaderStageFlagBits, :(
        @bitmask_flag ShaderStageFlag::UInt32 begin
            SHADER_STAGE_VERTEX_BIT = 1
            SHADER_STAGE_TESSELLATION_CONTROL_BIT = 2
            SHADER_STAGE_TESSELLATION_EVALUATION_BIT = 4
            SHADER_STAGE_GEOMETRY_BIT = 8
            SHADER_STAGE_FRAGMENT_BIT = 16
            SHADER_STAGE_COMPUTE_BIT = 32
            SHADER_STAGE_RAYGEN_BIT_KHR = 256
            SHADER_STAGE_ANY_HIT_BIT_KHR = 512
            SHADER_STAGE_CLOSEST_HIT_BIT_KHR = 1024
            SHADER_STAGE_MISS_BIT_KHR = 2048
            SHADER_STAGE_INTERSECTION_BIT_KHR = 4096
            SHADER_STAGE_CALLABLE_BIT_KHR = 8192
            SHADER_STAGE_TASK_BIT_NV = 64
            SHADER_STAGE_MESH_BIT_NV = 128
            SHADER_STAGE_SUBPASS_SHADING_BIT_HUAWEI = 16384
            SHADER_STAGE_ALL_GRAPHICS = $(Int(0x0000001f))
            SHADER_STAGE_ALL = $(Int(0x7fffffff))
        end
    ))

    test(BitmaskDefinition, bitmask_by_name, :VkCullModeFlagBits, :(
        @bitmask_flag CullModeFlag::UInt32 begin
            CULL_MODE_FRONT_BIT = 1
            CULL_MODE_BACK_BIT = 2
            CULL_MODE_NONE = 0
            CULL_MODE_FRONT_AND_BACK = 3
        end
    ))
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
module adapters_mysql_execute_result_test

using Test # @test
using Octo.Adapters.MySQL # Repo Schema Raw USE INSERT

Repo.debug_sql()

include("options.jl")

Repo.connect(;
    adapter = Octo.Adapters.MySQL,
    Options.for_mysql...
)

struct Employee
end
Schema.model(Employee, table_name="Employee", primary_key="ID")

result = Repo.execute([DROP TABLE IF EXISTS Employee])
@test result === nothing

result = Repo.execute(Raw("""CREATE TABLE IF NOT EXISTS Employee
                 (
                     ID INT NOT NULL AUTO_INCREMENT,
                     Name VARCHAR(255),
                     Salary FLOAT(7,2),
                     JoinDate DATE,
                     LastLogin DATETIME,
                     LunchTime TIME,
                     OfficeNo TINYINT,
                     JobType ENUM('HR', 'Management', 'Accounts'),
                     Senior BIT(1),
                     empno SMALLINT,
                     PRIMARY KEY (ID)
                 );"""))
@test result === nothing

result = Repo.execute(Raw("""INSERT INTO Employee (Name, Salary, JoinDate, LastLogin, LunchTime, OfficeNo, JobType, Senior, empno)
                 VALUES
                 ('John', 10000.50, '2015-8-3', '2015-9-5 12:31:30', '12:00:00', 1, 'HR', b'1', 1301),
                 ('Tom', 20000.25, '2015-8-4', '2015-10-12 13:12:14', '13:00:00', 12, 'HR', b'1', 1422),
                 ('Jim', 30000.00, '2015-6-2', '2015-9-5 10:05:10', '12:30:00', 45, 'Management', b'0', 1567);
              """))
@test result.num_affected_rows == 3
inserted = Repo.execute_result(INSERT)
@test inserted.id == 1

changes = (Name="Tim", Salary=15000.50, JoinDate="2015-7-25", LastLogin="2015-10-10 12:12:25", LunchTime="12:30:00", OfficeNo=56, JobType="Accounts", empno=3200)
inserted = Repo.insert!(Employee, changes)
@test inserted.id == 4
@test inserted.num_affected_rows == 1

result = Repo.delete!(Employee, 1:5)
@test result.num_affected_rows == 4

result = Repo.execute("update Employee set Name = 'New Name' where ID > 100")
@test result.num_affected_rows == 0

result = Repo.execute("delete from Employee")
@test result.num_affected_rows == 0

Repo.disconnect()

end # module adapters_mysql_execute_result_test


text:
using GLAbstraction, GLWindow, GeometryTypes, ColorTypes, GLFW, Reactive, ModernGL
const vert = vert"""
{{GLSL_VERSION}}
in vec2 pos;
out vec2 g_position;

void main() {
    g_position = pos;
}
"""

const frag = frag"""
{{GLSL_VERSION}}

out vec4 outColor;

void main() {
    outColor = vec4(1,0,0,1);
}
"""

const geom = geom"""

{{GLSL_VERSION}}

layout(points) in;
layout(triangle_strip, max_vertices = 4) out;
in vec2 g_position[];

void main(void)
{
  // get the four vertices passed to the shader:
  vec2 p0 = g_position[0];   // start of previous segment

  gl_Position = vec4(p0, 0, 1);
  EmitVertex();

  gl_Position = vec4(p0+vec2(0,0.1), 0, 1);
  EmitVertex();

  gl_Position = vec4(p0+vec2(0.1,0), 0, 1);
  EmitVertex();

  gl_Position = vec4(p0+vec2(0.1), 0, 1);
  EmitVertex();

  EndPrimitive();
}
"""
const window = GLWindow.create_glcontext("Geometry Shader")

const b = Point2f0[(-0.5,0),(0.0, 0.0),(0.4, 0.3)]

data = Dict{Symbol, Any}(
    :pos => GLBuffer(b),
)

program = GLAbstraction.LazyShader(vert, geom, frag)
robj = std_renderobject(data, program, Signal(AABB(Vec3f0(0), Vec3f0(1))), GL_POINTS)
println(robj)
glClearColor(0,0,0,1)
glClearColor(0, 0, 0, 1)

while isopen(window)
  	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    render(robj)
  	swapbuffers(window)
  	pollevents()
end


text:
module layers
using LinearAlgebra
export halfspace,layer_plain,layer_patterned, modes_freespace,Layer,Halfspace

"""
    Layer(thi,V,W,X,q)
    Computed Eigenmodes and propagation properties of a layer
    thi: layer thickness
    V: coordinate transform between eigenmode amplitude and magnetic field
    W: coordinate transform between eigenmode amplitude and electric field
    X: transmission of amplitude vector through the layer
    q: effective propagation constant of the eigenmodes
"""
struct Layer
    thi::Float64
    V::AbstractArray{Complex{Float64},2}
    W::AbstractArray{Complex{Float64},2}
    X::AbstractArray{Complex{Float64},2}
    q::AbstractArray{Complex{Float64},2}
end
"""
Halfspace(V,W,Kz)
Computed Eigenmodes of an unpatterned, infinitely thin halfspace
V: coordinate transform between eigenmode amplitude and magnetic field
W: coordinate transform between eigenmode amplitude and electric field
Kz: z-axis component of the propagation vector
"""
struct Halfspace
    V::AbstractArray{Complex{Float64},2}
    W::AbstractArray{Complex{Float64},2}
    Kz::AbstractArray{Complex{Float64},2}
end
"""
modes_freespace(Kx,Ky)
Computes the eigenmodes of propagation through free space, for normalization
Kx: x-axis component of the propagation vector
Ky: y-axis component of the propagation vector
returns
V0: coordinate transform between free space eigenmode amplitude and magnetic field
Kz0: z-axis component of the propagation vector in free space
"""
function modes_freespace(Kx,Ky)
    #just because |k|=1
    Kz0=sqrt.(Complex.(I-Kx*Kx-Ky*Ky))
    #P0 is identity
    Q0=[Kx*Ky I-Kx*Kx;Ky*Ky-I -Ky*Kx]
    #propagation
    q0=1im*Kz0
    q0=[q0 q0*0;0*q0 q0]
    #Free space, so W is identity
    #W0=I+0*Q0
    V0=Q0/Diagonal(q0)
    return V0,Kz0
end
"""
modes_patterned(Kx,Ky,k0,thi,epsilon)
Computes the eigenmodes of propagation through a patterned layer
Kx: x-axis component of the propagation vector
Ky: y-axis component of the propagation vector
k0: is 2pi/lambda, to normalize the thickness
thi: layer thickness
epsilon: Permittivity of layer in reciprocal space
returns
V: coordinate transform between eigenmode amplitude and magnetic field in the layer
W: coordinate transform between eigenmode amplitude and electric field in the layer
X: transmission of amplitude vector through the layer
q: effective propagation constant of the eigenmodes
"""
function modes_patterned(Kx,Ky,k0,thi,epsilon)
    eta=I/epsilon
    P=[Kx*eta*Ky I-Kx*eta*Kx;Ky*eta*Ky-I -Ky*eta*Kx]
    Q=[Kx*Ky epsilon-Kx*Kx;Ky*Ky-epsilon -Ky*Kx]
    #eigenmodes
    ev=eigen(Matrix(P*Q))
    q=Diagonal(sqrt.(Complex.(ev.values)))
    q[real.(q).>0].*=-1
    #W is transform between amplitude vector and E-Field
    W=ev.vectors
    #V is transform between amplitude vector and H-Field
    V=Q*W/Diagonal(q)

    #X the factor applied to the amplitudes when propagatin through the layer
    X=exp(q*k0*thi)
    return V,W,X,q
end
"""
layer_patterned(Kx,Ky,k0,thi,epsilon)
Computes the eigenmodes of propagation through a patterned layer in a Layer object
Kx: x-axis component of the propagation vector
Ky: y-axis component of the propagation vector
k0: is 2pi/lambda, to normalize the thickness
thi: layer thickness
epsilon: Permittivity of layer in reciprocal space
returns
patterned layer object
"""
function layer_patterned(Kx,Ky,k0,thi,epsilon)
    V,W,X,q=modes_patterned(Kx,Ky,k0,thi,epsilon)
    return Layer(thi,V,W,X,q)
end
"""
modes_plain(Kx,Kyk0,thi,epsilon)
Computes the eigenmodes of propagation through a plain layer
Kx: x-axis component of the propagation vector
Ky: y-axis component of the propagation vector
k0: is 2pi/lambda, to normalize the thickness
thi: layer thickness
epsilon: Bulk permittivity of the layer
returns
V: coordinate transform between eigenmode amplitude and magnetic field in the layer
W: coordinate transform between eigenmode amplitude and electric field in the layer
X: transmission of amplitude vector through the layer
q: effective propagation constant of the eigenmodes
"""
function modes_plain(Kx,Ky,k0,thi,epsilon)
    #Kz=P*Q
    Kz=sqrt.(Complex.(epsilon*I-Kx*Kx-Ky*Ky))
    Q=[Kx*Ky epsilon*I-Kx*Kx;Ky*Ky-epsilon*I -Ky*Kx]
    q=[1im*Kz zeros(size(Kz));zeros(size(Kz)) 1im*Kz]
    q[real.(q).>0].*=-1
    #W is identity
    W=I
    V=Q/Diagonal(q)
    X=exp(Matrix(q*k0*thi))
    return V,W,X,q
end
"""
layer_plain(Kx,Kyk0,thi,epsilon)
Computes the eigenmodes of propagation through a plain layer, in a Layer object
Kx: x-axis component of the propagation vector
Ky: y-axis component of the propagation vector
k0: is 2pi/lambda, to normalize the thickness
thi: layer thickness
epsilon: Bulk permittivity of the layer
returns
layer object
"""
function layer_plain(Kx,Ky,k0,thi,epsilon)
    V,W,X,q=modes_plain(Kx,Ky,k0,thi,epsilon)
    return Layer(thi,V,V*0+W,X,q)
end
"""
modes_halfspace(Kx,Ky,epsilon)
Computes the eigenmodes of an infinitely thin halfspace
Kx: x-axis component of the propagation vector
Ky: y-axis component of the propagation vector
epsilon: bulk permittivity of the halfspace
returns
V: coordinate transform between eigenmode amplitude and magnetic field in the halfspace
W: coordinate transform between eigenmode amplitude and electric field in the halfspace
Kz: z-axis component of the propagation vector in the halfspace
"""
function modes_halfspace(Kx,Ky,epsilon)
    Kz=sqrt.(Complex.(epsilon*I-Kx*Kx-Ky*Ky))
    Q=[Kx*Ky epsilon*I-Kx*Kx;Ky*Ky-epsilon*I -Ky*Kx]
    q0=[1im*Kz zeros(size(Kz));zeros(size(Kz)) 1im*Kz]
    W=I
    V=Q/Diagonal(q0)
    return V,W+0*V,Kz
end
"""
halfspace(Kx,Ky,epsilon)
Computes the eigenmodes of an infinitely thin halfspace in a Halfspace object
Kx: x-axis component of the propagation vector
Ky: y-axis component of the propagation vector
epsilon: bulk permittivity of the halfspace
returns
halfspace object
"""
function halfspace(Kx,Ky,epsilon)
    V,W,Kz=modes_halfspace(Kx,Ky,epsilon)
    return Halfspace(V,W,Kz)
end

end


text:
@testset "Reduction #3 per Harris" begin
    num_samples = 2500
    num_ants = 1
    num_correlators = 3
    correlator = EarlyPromptLateCorrelator(NumAnts(num_ants), NumAccumulators(num_correlators))
    correlator_sample_shifts = get_correlator_sample_shifts(GPSL1(), correlator, 2.5e6Hz, 0.5)
    input = StructArray{ComplexF32}(
        (
            CUDA.ones(Float32, (num_samples, num_ants, num_correlators)),
            CUDA.zeros(Float32, (num_samples, num_ants, num_correlators))
        )
    )
    threads_per_block = 256
    blocks_per_grid = cld(num_samples, threads_per_block)
    accum = StructArray{ComplexF32}(
        (
            CUDA.zeros(Float32, (blocks_per_grid, num_ants, num_correlators)),
            CUDA.zeros(Float32, (blocks_per_grid, num_ants, num_correlators))
        )
    )
    shmem_size = sizeof(Float32) * threads_per_block
    for corr_idx = 1:num_correlators
        # re samples
        @cuda threads=threads_per_block blocks=blocks_per_grid shmem=shmem_size reduce_3(
            view(accum.re, :, :, corr_idx),
            view(input.re, :, :, corr_idx),
            num_samples
        )
        # im samples
        @cuda threads=threads_per_block blocks=blocks_per_grid shmem=shmem_size reduce_3(
            view(accum.im, :, :, corr_idx),
            view(input.im, :, :, corr_idx),
            num_samples
        )
    end
    for corr_idx = 1:num_correlators
        # re samples
        @cuda threads=threads_per_block blocks=1 shmem=shmem_size reduce_3(
            view(accum.re, :, :, corr_idx),
            view(accum.re, :, :, corr_idx),
            size(accum, 1)
        )
        # im samples
        @cuda threads=threads_per_block blocks=1 shmem=shmem_size reduce_3(
            view(accum.im, :, :, corr_idx),
            view(accum.im, :, :, corr_idx),
            size(accum, 1)
        )
    end
    CUDA.@allowscalar begin
        accum_true = ComplexF32[num_samples num_samples num_samples]
        @test Array(accum)[1, :, :,] ≈ accum_true
    end
    # @benchmark CUDA.@sync begin
    #     for corr_idx = 1:$num_correlators
    #         # re samples
    #         @cuda threads=$threads_per_block blocks=$blocks_per_grid shmem=$shmem_size $reduce_3(
    #             $view(accum.re, :, :, corr_idx),
    #             $view(input.re, :, :, corr_idx),
    #             $num_samples
    #         )
    #         # im samples
    #         @cuda threads=$threads_per_block blocks=$blocks_per_grid shmem=$shmem_size $reduce_3(
    #             $view(accum.im, :, :, corr_idx),
    #             $view(input.im, :, :, corr_idx),
    #             $num_samples
    #         )
    #     end
    #     for corr_idx = 1:$num_correlators
    #         # re samples
    #         @cuda threads=$threads_per_block blocks=1 shmem=$shmem_size $reduce_3(
    #             $view(accum.re, :, :, corr_idx),
    #             $view(accum.re, :, :, corr_idx),
    #             $size(accum, 1)
    #         )
    #         # im samples
    #         @cuda threads=$threads_per_block blocks=1 shmem=$shmem_size $reduce_3(
    #             $view(accum.im, :, :, corr_idx),
    #             $view(accum.im, :, :, corr_idx),
    #             $size(accum, 1)
    #         )
    #     end
    # end
end

@testset "Complex Reduction #3 per Harris" begin
    num_samples = 2500
    num_ants = 1
    num_correlators = 3
    correlator = EarlyPromptLateCorrelator(NumAnts(num_ants), NumAccumulators(num_correlators))
    correlator_sample_shifts = get_correlator_sample_shifts(GPSL1(), correlator, 2.5e6Hz, 0.5)
    input = StructArray{ComplexF32}(
        (
            CUDA.ones(Float32, (num_samples, num_ants, num_correlators)),
            CUDA.zeros(Float32, (num_samples, num_ants, num_correlators))
        )
    )
    threads_per_block = 256
    blocks_per_grid = cld(num_samples, threads_per_block)
    accum = StructArray{ComplexF32}(
        (
            CUDA.zeros(Float32, (blocks_per_grid, num_ants, num_correlators)),
            CUDA.zeros(Float32, (blocks_per_grid, num_ants, num_correlators))
        )
    )
    shmem_size = sizeof(ComplexF32) * threads_per_block
    for corr_idx = 1:num_correlators
        @cuda threads=threads_per_block blocks=blocks_per_grid shmem=shmem_size reduce_cplx_3(
            view(accum.re, :, :, corr_idx),
            view(accum.im, :, :, corr_idx),
            view(input.re, :, :, corr_idx),
            view(input.im, :, :, corr_idx),
            num_samples
        )
    end
    for corr_idx = 1:num_correlators
        @cuda threads=threads_per_block blocks=1 shmem=shmem_size reduce_cplx_3(
            view(accum.re, :, :, corr_idx),
            view(accum.im, :, :, corr_idx),
            view(accum.re, :, :, corr_idx),
            view(accum.im, :, :, corr_idx),
            size(accum, 1)
        )
    end
    CUDA.@allowscalar begin
        accum_true = ComplexF32[num_samples num_samples num_samples]
        @test Array(accum)[1, :, :,] ≈ accum_true
    end
    # @benchmark CUDA.@sync begin
    #     for corr_idx = 1:$num_correlators
    #         @cuda threads=$threads_per_block blocks=$blocks_per_grid shmem=$shmem_size $reduce_cplx_3(
    #             $view(accum.re, :, :, corr_idx),
    #             $view(accum.im, :, :, corr_idx),
    #             $view(input.re, :, :, corr_idx),
    #             $view(input.im, :, :, corr_idx),
    #             $num_samples
    #         )
    #     end
    #     for corr_idx = 1:$num_correlators
    #         @cuda threads=$threads_per_block blocks=1 shmem=$shmem_size $reduce_cplx_3(
    #             $view(accum.re, :, :, corr_idx),
    #             $view(accum.im, :, :, corr_idx),
    #             $view(accum.re, :, :, corr_idx),
    #             $view(accum.im, :, :, corr_idx),
    #             $size(accum, 1)
    #         )
    #     end
    # end
end

@testset "Complex Multi Reduction #3 per Harris" begin
    num_samples = 2500
    num_ants = 1
    num_correlators = 3
    correlator = EarlyPromptLateCorrelator(NumAnts(num_ants), NumAccumulators(num_correlators))
    correlator_sample_shifts = get_correlator_sample_shifts(GPSL1(), correlator, 2.5e6Hz, 0.5)
    input = StructArray{ComplexF32}(
        (
            CUDA.ones(Float32, (num_samples, num_ants, num_correlators)),
            CUDA.zeros(Float32, (num_samples, num_ants, num_correlators))
        )
    )
    threads_per_block = 256
    blocks_per_grid = cld(num_samples, threads_per_block)
    accum = StructArray{ComplexF32}(
        (
            CUDA.zeros(Float32, (blocks_per_grid, num_ants, num_correlators)),
            CUDA.zeros(Float32, (blocks_per_grid, num_ants, num_correlators))
        )
    )
    shmem_size = sizeof(ComplexF32) * threads_per_block * num_ants * num_correlators
    @cuda threads=threads_per_block blocks=blocks_per_grid shmem=shmem_size reduce_cplx_multi_31(
        accum.re,
        accum.im,
        input.re,
        input.im,
        num_samples,
        NumAnts(num_ants),
        correlator_sample_shifts
    )
    @cuda threads=threads_per_block blocks=1 shmem=shmem_size reduce_cplx_multi_31(
        accum.re,
        accum.im,
        accum.re,
        accum.im,
        size(accum, 1),
        NumAnts(num_ants),
        correlator_sample_shifts
    )
    CUDA.@allowscalar begin
        accum_true = ComplexF32[num_samples num_samples num_samples]
        @test Array(accum)[1, :, :,] ≈ accum_true
    end
    # @benchmark CUDA.@sync begin
    #     @cuda threads=$threads_per_block blocks=$blocks_per_grid shmem=$shmem_size $reduce_cplx_multi_3(
    #         $accum.re,
    #         $accum.im,
    #         $input.re,
    #         $input.im,
    #         $num_samples,
    #         $NumAnts(num_ants),
    #         $correlator_sample_shifts
    #     )
    #     @cuda threads=$threads_per_block blocks=1 shmem=$shmem_size $reduce_cplx_multi_3(
    #         $accum.re,
    #         $accum.im,
    #         $accum.re,
    #         $accum.im,
    #         $size(accum, 1),
    #         $NumAnts(num_ants),
    #         $correlator_sample_shifts
    #     )
    # end
end

@testset "Reduction #4 per Harris" begin
    num_samples = 2500
    num_ants = 1
    num_correlators = 3
    correlator = EarlyPromptLateCorrelator(NumAnts(num_ants), NumAccumulators(num_correlators))
    correlator_sample_shifts = get_correlator_sample_shifts(GPSL1(), correlator, 2.5e6Hz, 0.5)
    input = StructArray{ComplexF32}(
        (
            CUDA.ones(Float32, (num_samples, num_ants, num_correlators)),
            CUDA.zeros(Float32, (num_samples, num_ants, num_correlators))
        )
    )
    threads_per_block = 256
    # only half the grid size for reduce_4
    blocks_per_grid = cld(num_samples, threads_per_block) ÷ 2
    accum = StructArray{ComplexF32}(
        (
            CUDA.zeros(Float32, (blocks_per_grid, num_ants, num_correlators)),
            CUDA.zeros(Float32, (blocks_per_grid, num_ants, num_correlators))
        )
    )
    shmem_size = sizeof(Float32) * threads_per_block
    for corr_idx = 1:num_correlators
        # re samples
        @cuda threads=threads_per_block blocks=blocks_per_grid shmem=shmem_size reduce_4(
            view(accum.re, :, :, corr_idx),
            view(input.re, :, :, corr_idx),
            num_samples
        )
        # im samples
        @cuda threads=threads_per_block blocks=blocks_per_grid shmem=shmem_size reduce_4(
            view(accum.im, :, :, corr_idx),
            view(input.im, :, :, corr_idx),
            num_samples
        )
    end
    for corr_idx = 1:num_correlators
        # re samples
        @cuda threads=threads_per_block÷2 blocks=1 shmem=shmem_size reduce_4(
            view(accum.re, :, :, corr_idx),
            view(accum.re, :, :, corr_idx),
            size(accum, 1)
        )
        # im samples
        @cuda threads=threads_per_block÷2 blocks=1 shmem=shmem_size reduce_4(
            view(accum.im, :, :, corr_idx),
            view(accum.im, :, :, corr_idx),
            size(accum, 1)
        )
    end
    CUDA.@allowscalar begin
        accum_true = ComplexF32[num_samples num_samples num_samples]
        @test Array(accum)[1, :, :,] ≈ accum_true
    end
    # @benchmark CUDA.@sync begin
    #     for corr_idx = 1:$num_correlators
    #         # re samples
    #         @cuda threads=$threads_per_block blocks=$blocks_per_grid shmem=$shmem_size $reduce_4(
    #             $view(accum.re, :, :, corr_idx),
    #             $view(input.re, :, :, corr_idx),
    #             $num_samples
    #         )
    #         # im samples
    #         @cuda threads=$threads_per_block blocks=$blocks_per_grid shmem=$shmem_size $reduce_4(
    #             $view(accum.im, :, :, corr_idx),
    #             $view(input.im, :, :, corr_idx),
    #             $num_samples
    #         )
    #     end
    #     for corr_idx = 1:$num_correlators
    #         # re samples
    #         @cuda threads=$threads_per_block÷2 blocks=1 shmem=$shmem_size $reduce_4(
    #             $view(accum.re, :, :, corr_idx),
    #             $view(accum.re, :, :, corr_idx),
    #             $size(accum, 1)
    #         )
    #         # im samples
    #         @cuda threads=$threads_per_block÷2 blocks=1 shmem=$shmem_size $reduce_4(
    #             $view(accum.im, :, :, corr_idx),
    #             $view(accum.im, :, :, corr_idx),
    #             $size(accum, 1)
    #         )
    #     end
    # end
end

@testset "Complex Reduction #4 per Harris" begin
    num_samples = 2500
    num_ants = 1
    num_correlators = 3
    correlator = EarlyPromptLateCorrelator(NumAnts(num_ants), NumAccumulators(num_correlators))
    correlator_sample_shifts = get_correlator_sample_shifts(GPSL1(), correlator, 2.5e6Hz, 0.5)
    input = StructArray{ComplexF32}(
        (
            CUDA.ones(Float32, (num_samples, num_ants, num_correlators)),
            CUDA.zeros(Float32, (num_samples, num_ants, num_correlators))
        )
    )
    threads_per_block = 256
    # only half the grid size for reduce_4
    blocks_per_grid = cld(num_samples, threads_per_block) ÷ 2
    accum = StructArray{ComplexF32}(
        (
            CUDA.zeros(Float32, (blocks_per_grid, num_ants, num_correlators)),
            CUDA.zeros(Float32, (blocks_per_grid, num_ants, num_correlators))
        )
    )
    shmem_size = sizeof(ComplexF32) * threads_per_block
    for corr_idx = 1:num_correlators
        @cuda threads=threads_per_block blocks=blocks_per_grid shmem=shmem_size reduce_cplx_4(
            view(accum.re, :, :, corr_idx),
            view(accum.im, :, :, corr_idx),
            view(input.re, :, :, corr_idx),
            view(input.im, :, :, corr_idx),
        )
    end
    for corr_idx = 1:num_correlators
        @cuda threads=threads_per_block blocks=1 shmem=shmem_size reduce_cplx_4(
            view(accum.re, :, :, corr_idx),
            view(accum.im, :, :, corr_idx),
            view(accum.re, :, :, corr_idx),
            view(accum.im, :, :, corr_idx),
        )
    end
    CUDA.@allowscalar begin
        accum_true = ComplexF32[num_samples num_samples num_samples]
        @test Array(accum)[1, :, :,] ≈ accum_true
    end
    # @benchmark CUDA.@sync begin
    #     for corr_idx = 1:$num_correlators
    #         @cuda threads=$threads_per_block blocks=$blocks_per_grid shmem=$shmem_size $reduce_cplx_4(
    #             $view(accum.re, :, :, corr_idx),
    #             $view(accum.im, :, :, corr_idx),
    #             $view(input.re, :, :, corr_idx),
    #             $view(input.im, :, :, corr_idx),
    #         )
    #     end
    #     for corr_idx = 1:$num_correlators
    #         @cuda threads=$threads_per_block blocks=$1 shmem=$shmem_size $reduce_cplx_4(
    #             $view(accum.re, :, :, corr_idx),
    #             $view(accum.im, :, :, corr_idx),
    #             $view(accum.re, :, :, corr_idx),
    #             $view(accum.im, :, :, corr_idx),
    #         )
    #     end
    # end
end

@testset "Complex Multi Reduction #4 per Harris" begin
    num_samples = 2500
    num_ants = 1
    num_correlators = 3
    correlator = EarlyPromptLateCorrelator(NumAnts(num_ants), NumAccumulators(num_correlators))
    correlator_sample_shifts = get_correlator_sample_shifts(GPSL1(), correlator, 2.5e6Hz, 0.5)
    input = StructArray{ComplexF32}(
        (
            CUDA.ones(Float32, (num_samples, num_ants, num_correlators)),
            CUDA.zeros(Float32, (num_samples, num_ants, num_correlators))
        )
    )
    threads_per_block = 256
    # only half the grid size for reduce_4
    blocks_per_grid = cld(num_samples, threads_per_block) ÷ 2
    accum = StructArray{ComplexF32}(
        (
            CUDA.zeros(Float32, (blocks_per_grid, num_ants, num_correlators)),
            CUDA.zeros(Float32, (blocks_per_grid, num_ants, num_correlators))
        )
    )
    shmem_size = sizeof(ComplexF32) * threads_per_block * num_ants * num_correlators
    @cuda threads=threads_per_block blocks=blocks_per_grid shmem=shmem_size reduce_cplx_multi_4(
        accum.re,
        accum.im,
        input.re,
        input.im,
        num_samples,
        NumAnts(num_ants),
        correlator_sample_shifts
    )
    @cuda threads=threads_per_block blocks=1 shmem=shmem_size reduce_cplx_multi_4(
        accum.re,
        accum.im,
        accum.re,
        accum.im,
        blocks_per_grid,
        NumAnts(num_ants),
        correlator_sample_shifts
    )
    CUDA.@allowscalar begin
        accum_true = ComplexF32[num_samples num_samples num_samples]
        @test Array(accum)[1, :, :,] ≈ accum_true
    end
    # @benchmark CUDA.@sync begin
    #     @cuda threads=$threads_per_block blocks=$blocks_per_grid shmem=$shmem_size reduce_cplx_multi_4(
    #         $accum.re,
    #         $accum.im,
    #         $input.re,
    #         $input.im,
    #         $num_samples,
    #         $NumAnts(num_ants),
    #         $correlator_sample_shifts
    #     )
    #     @cuda threads=$threads_per_block blocks=1 shmem=$shmem_size reduce_cplx_multi_4(
    #         $accum.re,
    #         $accum.im,
    #         $accum.re,
    #         $accum.im,
    #         $blocks_per_grid,
    #         $NumAnts(num_ants),
    #         $correlator_sample_shifts
    #     )
    # end
end

@testset "Reduction #5 per Harris" begin
    num_samples = 2048
    num_ants = 1
    num_correlators = 3
    signal_duration = 0.001
    sampling_frequency = (num_samples/signal_duration)Hz
    correlator = EarlyPromptLateCorrelator(NumAnts(num_ants), NumAccumulators(num_correlators))
    correlator_sample_shifts = get_correlator_sample_shifts(GPSL1(), correlator, sampling_frequency, 0.5)
    input = StructArray{ComplexF32}(
        (
            CUDA.ones(Float32, (num_samples, num_ants, num_correlators)),
            CUDA.zeros(Float32, (num_samples, num_ants, num_correlators))
        )
    )
    threads_per_block = 256
    # only half the grid size for reduce_4
    blocks_per_grid = cld(cld(num_samples, threads_per_block), 2)
    accum = StructArray{ComplexF32}(
        (
            CUDA.zeros(Float32, (blocks_per_grid, num_ants, num_correlators)),
            CUDA.zeros(Float32, (blocks_per_grid, num_ants, num_correlators))
        )
    )
    shmem_size = sizeof(Float32) * threads_per_block
    for corr_idx = 1:num_correlators
        # re samples
        @cuda threads=threads_per_block blocks=blocks_per_grid shmem=shmem_size reduce_5(
            view(accum.re, :, :, corr_idx),
            view(input.re, :, :, corr_idx),
            num_samples
        )
        # im samples
        @cuda threads=threads_per_block blocks=blocks_per_grid shmem=shmem_size reduce_5(
            view(accum.im, :, :, corr_idx),
            view(input.im, :, :, corr_idx),
            num_samples
        )
    end
    for corr_idx = 1:num_correlators
        # re samples
        @cuda threads=threads_per_block blocks=1 shmem=shmem_size reduce_5(
            view(accum.re, :, :, corr_idx),
            view(accum.re, :, :, corr_idx),
            size(accum, 1)
        )
        # im samples
        @cuda threads=threads_per_block blocks=1 shmem=shmem_size reduce_5(
            view(accum.im, :, :, corr_idx),
            view(accum.im, :, :, corr_idx),
            size(accum, 1)
        )
    end
    CUDA.@allowscalar begin
        accum_true = ComplexF32[num_samples num_samples num_samples]
        @test Array(accum)[1, :, :,] ≈ accum_true
    end
    # @benchmark CUDA.@sync begin
    #     for corr_idx = 1:$num_correlators
    #         # re samples
    #         @cuda threads=$threads_per_block blocks=$blocks_per_grid shmem=$shmem_size $reduce_5(
    #             $view(accum.re, :, :, corr_idx),
    #             $view(input.re, :, :, corr_idx),
    #             $num_samples
    #         )
    #         # im samples
    #         @cuda threads=$threads_per_block blocks=$blocks_per_grid shmem=$shmem_size $reduce_5(
    #             $view(accum.im, :, :, corr_idx),
    #             $view(input.im, :, :, corr_idx),
    #             $num_samples
    #         )
    #     end
    #     for corr_idx = 1:$num_correlators
    #         # re samples
    #         @cuda threads=$threads_per_block blocks=1 shmem=$shmem_size $reduce_5(
    #             $view(accum.re, :, :, corr_idx),
    #             $view(accum.re, :, :, corr_idx),
    #             $size(accum, 1)
    #         )
    #         # im samples
    #         @cuda threads=$threads_per_block blocks=1 shmem=$shmem_size $reduce_5(
    #             $view(accum.im, :, :, corr_idx),
    #             $view(accum.im, :, :, corr_idx),
    #             $size(accum, 1)
    #         )
    #     end
    # end
end

@testset "Complex Reduction #5 per Harris" begin
    num_samples = 2048
    num_ants = 1
    num_correlators = 3
    signal_duration = 0.001
    sampling_frequency = (num_samples/signal_duration)Hz
    correlator = EarlyPromptLateCorrelator(NumAnts(num_ants), NumAccumulators(num_correlators))
    correlator_sample_shifts = get_correlator_sample_shifts(GPSL1(), correlator, sampling_frequency, 0.5)
    input = StructArray{ComplexF32}(
        (
            CUDA.ones(Float32, (num_samples, num_ants, num_correlators)),
            CUDA.zeros(Float32, (num_samples, num_ants, num_correlators))
        )
    )
    threads_per_block = 256
    # only half the grid size for reduce_5
    blocks_per_grid = cld(cld(num_samples, threads_per_block), 2)
    accum = StructArray{ComplexF32}(
        (
            CUDA.zeros(Float32, (blocks_per_grid, num_ants, num_correlators)),
            CUDA.zeros(Float32, (blocks_per_grid, num_ants, num_correlators))
        )
    )
    shmem_size = sizeof(ComplexF32) * threads_per_block
    for corr_idx = 1:num_correlators
        @cuda threads=threads_per_block blocks=blocks_per_grid shmem=shmem_size reduce_cplx_5(
            view(accum.re, :, :, corr_idx),
            view(accum.im, :, :, corr_idx),
            view(input.re, :, :, corr_idx),
            view(input.im, :, :, corr_idx),
            num_samples
        )
    end
    Array(accum)
    for corr_idx = 1:num_correlators
        @cuda threads=threads_per_block blocks=1 shmem=shmem_size reduce_cplx_5(
            view(accum.re, :, :, corr_idx),
            view(accum.im, :, :, corr_idx),
            view(accum.re, :, :, corr_idx),
            view(accum.im, :, :, corr_idx),
            blocks_per_grid
        )
    end
    Array(accum)
    CUDA.@allowscalar begin
        accum_true = ComplexF32[num_samples num_samples num_samples]
        @test Array(accum)[1, :, :,] ≈ accum_true
    end
    @benchmark CUDA.@sync begin
        for corr_idx = 1:$num_correlators
            @cuda threads=$threads_per_block blocks=$blocks_per_grid shmem=$shmem_size $reduce_cplx_5(
                $view($accum.re, :, :, corr_idx),
                $view($accum.im, :, :, corr_idx),
                $view($input.re, :, :, corr_idx),
                $view($input.im, :, :, corr_idx),
                $num_samples
            )
        end
        for corr_idx = 1:$num_correlators
            @cuda threads=$threads_per_block blocks=$1 shmem=$shmem_size $reduce_cplx_5(
                $view($accum.re, :, :, corr_idx),
                $view($accum.im, :, :, corr_idx),
                $view($accum.re, :, :, corr_idx),
                $view($accum.im, :, :, corr_idx),
                $blocks_per_grid
            )
        end
    end
end

@testset "Complex Multi Reduction #5 per Harris" begin
    num_samples = 2048
    num_ants = 1
    num_correlators = 3
    signal_duration = 0.001
    sampling_frequency = (num_samples/signal_duration)Hz
    correlator = EarlyPromptLateCorrelator(NumAnts(num_ants), NumAccumulators(num_correlators))
    correlator_sample_shifts = get_correlator_sample_shifts(GPSL1(), correlator, sampling_frequency, 0.5)
    input = StructArray{ComplexF32}(
        (
            CUDA.ones(Float32, (num_samples, num_ants, num_correlators)),
            CUDA.zeros(Float32, (num_samples, num_ants, num_correlators))
        )
    )
    threads_per_block = 256
    # only half the grid size for reduce_5
    blocks_per_grid = cld(cld(num_samples, threads_per_block), 2)
    accum = StructArray{ComplexF32}(
        (
            CUDA.zeros(Float32, (blocks_per_grid, num_ants, num_correlators)),
            CUDA.zeros(Float32, (blocks_per_grid, num_ants, num_correlators))
        )
    )
    shmem_size = sizeof(ComplexF32) * threads_per_block * num_ants * num_correlators
    @cuda threads=threads_per_block blocks=blocks_per_grid shmem=shmem_size reduce_cplx_multi_5(
        accum.re,
        accum.im,
        input.re,
        input.im,
        num_samples,
        NumAnts(num_ants),
        correlator_sample_shifts
    )
    Array(accum)
    @cuda threads=threads_per_block blocks=1 shmem=shmem_size reduce_cplx_multi_5(
        accum.re,
        accum.im,
        accum.re,
        accum.im,
        blocks_per_grid,
        NumAnts(num_ants),
        correlator_sample_shifts
    )
    Array(accum)
    CUDA.@allowscalar begin
        accum_true = ComplexF32[num_samples num_samples num_samples]
        @test Array(accum)[1, :, :,] ≈ accum_true
    end
    # @benchmark CUDA.@sync begin
    #     @cuda threads=$threads_per_block blocks=$blocks_per_grid shmem=$shmem_size reduce_cplx_multi_4(
    #         $accum.re,
    #         $accum.im,
    #         $input.re,
    #         $input.im,
    #         $num_samples,
    #         $NumAnts(num_ants),
    #         $correlator_sample_shifts
    #     )
    #     @cuda threads=$threads_per_block blocks=1 shmem=$shmem_size reduce_cplx_multi_4(
    #         $accum.re,
    #         $accum.im,
    #         $accum.re,
    #         $accum.im,
    #         $blocks_per_grid,
    #         $NumAnts(num_ants),
    #         $correlator_sample_shifts
    #     )
    # end
end

# @testset "Complex Multi NANT Reduction #5 per Harris" begin
    num_samples = 2 ^ 15
    num_ants = 16
    num_correlators = 1
    signal_duration = 0.001
    sampling_frequency = (num_samples/signal_duration)Hz
    # correlator = EarlyPromptLateCorrelator(NumAnts(num_ants), NumAccumulators(num_correlators))
    # correlator_sample_shifts = get_correlator_sample_shifts(GPSL1(), correlator, sampling_frequency, 0.5)
    input = StructArray{ComplexF32}(
        (
            CUDA.ones(Float32, (num_samples, num_ants)),
            CUDA.zeros(Float32, (num_samples, num_ants))
        )
    )
    threads_per_block = 1024
    # only half the grid size for reduce_5
    blocks_per_grid_x = cld(cld(num_samples, threads_per_block), 2)
    blocks_per_grid_y = num_ants
    blocks_per_grid = (blocks_per_grid_x, blocks_per_grid_y)
    accum = StructArray{ComplexF32}(
        (
            CUDA.zeros(Float32, (blocks_per_grid_x, blocks_per_grid_y)),
            CUDA.zeros(Float32, (blocks_per_grid_x, blocks_per_grid_y))
        )
    )
    shmem_size = sizeof(ComplexF32) * threads_per_block

    @cuda threads=threads_per_block blocks=blocks_per_grid shmem=shmem_size reduce_cplx_multi_nant_5(
        accum.re,
        accum.im,
        input.re,
        input.im,
        num_samples,
        NumAnts(num_ants),
    )
    Array(accum)
    @cuda threads=threads_per_block blocks=(1, blocks_per_grid_y) shmem=shmem_size reduce_cplx_multi_nant_5(
        accum.re,
        accum.im,
        accum.re,
        accum.im,
        blocks_per_grid_x,
        NumAnts(num_ants)
    )
    Array(accum)
    CUDA.@allowscalar begin
        accum_true = ComplexF32[num_samples num_samples num_samples]
        @test Array(accum)[1, :, :,] ≈ accum_true
    end
    @benchmark CUDA.@sync begin
        @cuda threads=$threads_per_block blocks=$blocks_per_grid shmem=$shmem_size $reduce_cplx_multi_nant_5(
            $accum.re,
            $accum.im,
            $input.re,
            $input.im,
            $num_samples,
            $NumAnts(num_ants),
        )
        @cuda threads=$threads_per_block blocks=(1, $blocks_per_grid_y) shmem=$shmem_size $reduce_cplx_multi_nant_5(
            $accum.re,
            $accum.im,
            $accum.re,
            $accum.im,
            $blocks_per_grid_x,
            $NumAnts(num_ants)
        )
    end
end

max_stars_repo_name: mfkiwl/GPUAcceleratedTracking

max_issues_repo_name: mfkiwl/GPUAcceleratedTracking

max_forks_repo_name: mfkiwl/GPUAcceleratedTracking


text:
## weights
function getweight(P, i, j, P2, TEs, M, maxmag, flags) # Phase, index, neighbor, ...
    weight = 1.0
    if flags[1] weight *= (0.1 + 0.9phasecoherence(P, i, j)) end
    if flags[2] weight *= (0.1 + 0.9phasegradientcoherence(P, P2, TEs, i, j)) end
    if flags[3] weight *= (0.1 + 0.9phaselinearity(P, i, j)) end
    if !isnothing(M)
        small, big = minmax(M[i], M[j])
        if flags[4] weight *= (0.1 + 0.9magcoherence(small, big)) end
        if flags[5] weight *= (0.1 + 0.9magweight(small, maxmag)) end
        if flags[6] weight *= (0.1 + 0.9magweight2(big, maxmag)) end
    end
    return weight
end

phasecoherence(P, i, j) = 1 - abs(γ(P[i] - P[j]) / π)
phasegradientcoherence(P, P2, TEs, i, j) = max(0, 1 - abs(γ(P[i] - P[j]) - γ(P2[i] - P2[j]) * TEs[1] / TEs[2]))
magcoherence(small, big) = (small / big) ^ 2
magweight(small, maxmag) = 0.5 + 0.5min(1, small / (0.5 * maxmag))
magweight2(big, maxmag) = 0.5 + 0.5min(1, (0.5 * maxmag) / big) # too high magnitude is not good either (flow artifact)

phaselinearity(P, i, j, k) = max(0, 1 - abs(rem2pi(P[i] - 2P[j] + P[k], RoundNearest)/2))
function phaselinearity(P, i, j)
    neighbor = j - i
    h = i - neighbor
    k = j + neighbor
    if 0 < h && k <= length(P)
        return phaselinearity(P, h, i, j) * phaselinearity(P, i, j, k)
    else
        return 0.9
    end
end

"""
    calculateweights(wrapped; weights=:romeo, kwargs...)

Calculates weights for all edges.
size(weights) == [3, size(wrapped)...]

###  Optional keyword arguments:

- `weights`: Options are [`:romeo`] | `:romeo2` | `:romeo3` | `:bestpath`.
- `mag`: Additional mag weights are used.
- `mask`: Unwrapping is only performed inside the mask.
- `phase2`: A second reference phase image (possibly with different echo time).
   It is used for calculating the phasecoherence weight.
- `TEs`: The echo times of the phase and the phase2 images as a tuple (eg. (5, 10) or [5, 10]).

"""
function calculateweights(wrapped; weights=:romeo, kwargs...)
    weights = if weights isa AbstractArray{<:Real,4}
        weights
    elseif weights == :bestpath
        calculateweights_bestpath(wrapped; kwargs...)
    else
        calculateweights_romeo(wrapped, weights; kwargs...)
    end
    # these edges do not exist
    weights[1,end,:,:] .= 0
    weights[2,:,end,:] .= 0
    weights[3,:,:,end] .= 0
    return weights
end

calculateweights_romeo(wrapped, weights::AbstractArray{T,4}; kw...) where T = UInt8.(weights)
# defines default weights
function calculateweights_romeo(wrapped, weights::Symbol; kwargs...)
    flags = falses(6)
    if weights == :romeo
        flags[1:4] .= true
    elseif weights == :romeo2
        flags[[1,4]] .= true # phasecoherence, magcoherence
    elseif weights == :romeo3
        flags[[1,2,4]] .= true # phasecoherence, phasegradientcoherence, magcoherence
    elseif weights == :romeo4
        flags[1:4] .= true
    elseif weights == :romeo6
        flags[1:6] .= true # additional magweight, magweight2
    else
        throw(ArgumentError("Weight '$weights' not defined!"))
    end
    return calculateweights_romeo(wrapped, flags; kwargs...)
end

function calculateweights_romeo(wrapped, flags::AbstractArray{Bool,1}; type::Type{T}=UInt8, rescale=rescale, kwargs...) where T
    mask, P2, TEs, M, maxmag = parsekwargs(kwargs, wrapped)
    updateflags!(flags, wrapped, P2, TEs, M)
    stridelist = getdimoffsets(wrapped)
    weights = zeros(T, 3, size(wrapped)...)
    for dim in 1:3
        neighbor = stridelist[dim]
        for I in LinearIndices(wrapped)
            J = I + neighbor
            if mask[I] && checkbounds(Bool, wrapped, J)
                w = getweight(wrapped, I, J, P2, TEs, M, maxmag, flags)
                weights[dim + (I-1)*3] = rescale(w)
            end
        end
    end
    return weights
end

## utility functions
function parsekwargs(kwargs, wrapped)
    getval(key) = if haskey(kwargs, key) kwargs[key] else nothing end
    mag = getval(:mag)
    mask = getval(:mask)
    if !isnothing(mask)
        if !isnothing(mag)
            mag = mag .* mask
        end
    else
        mask = trues(size(wrapped))
    end
    maxmag = if !isnothing(mag)
        maximum(mag[isfinite.(mag)]) else nothing
    end
    return mask, getval(:phase2), getval(:TEs), mag, maxmag
end

function updateflags!(flags, P, P2, TEs, M)
    if isnothing(M)
        flags[4:6] .= false
    end
    if isnothing(P2) || isnothing(TEs)
        flags[2] = false
    end
end

# from: 1 is best and 0 worst
# to: 1 is best, NBINS is worst, 0 is not valid (not added to queue)
function rescale(w)
    if 0 ≤ w ≤ 1
        max(round(Int, (1 - w) * (NBINS - 1)), 1)
    else
        0
    end
end

## best path weights
# Abdul-Rahamn https://doi.org/10.1364/AO.46.006623

function calculateweights_bestpath(wrapped; kwargs...)
    scale(w) = UInt8.(min(max(round((1 - (w / 10)) * (NBINS - 1)), 1), 255))
    weights = scale.(getbestpathweight(wrapped))
    if haskey(kwargs, :mask) # apply mask to weights
        mask = kwargs[:mask]
        weights .*= reshape(mask, 1, size(mask)...)
    end
    weights
end

function getbestpathweight(φ)
    R = 1 ./ getD(φ)
    weight = zeros(3, size(R)...)
    for idim in 1:3
        n = getdimoffsets(R)[idim]
        @inbounds for i in 1:length(R)-n
            weight[idim + 3i] = R[i] + R[i+n]
        end
    end
    return weight
end

function getD(φ)
    directions = Iterators.product(-1:1,-1:1,-1:1)
    neighbors = unique(abs(sum(getdimoffsets(φ) .* d)) for d in directions if d != (0,0,0))
    D2 = zeros(size(φ))
    @inbounds for n in neighbors, i in 1+n:length(φ)-n
        D2[i] += (γ(φ[i-n] - φ[i]) - γ(φ[i] - φ[i+n])) ^ 2
    end
    return sqrt.(D2)
end


text:
using Gadfly

p1 = plot(x=[1,1,1], y=[1,2,3], size=[1mm,2mm,4mm], Geom.point,
    label=["qwerQWER","asdfASDF","zxcvZXCV"], Geom.label(position=:dynamic, hide_overlaps=false));
p2 = plot(x=[1,1,1], y=[1,2,3], Geom.point,
    label=["qwerQWER","asdfASDF","zxcvZXCV"], Geom.label(position=:dynamic, hide_overlaps=false));
p3 = plot(x=[1,1,1], y=[1,2,3], size=[1mm,2mm,4mm], Geom.point,
    label=["qwerQWER","asdfASDF","zxcvZXCV"], Geom.label(position=:right, hide_overlaps=false));
p4 = plot(x=[1,1,1], y=[1,2,3], Geom.point,
    label=["qwerQWER","asdfASDF","zxcvZXCV"], Geom.label(position=:right, hide_overlaps=false));
gridstack([p1 p2; p3 p4])


text:
@testset "dat" begin
    d = dat(2003, 6, 1, 0.0)
    @test d ≈ 32.0 atol=1e-9
    d = dat(2008, 1, 17, 0.0)
    @test d ≈ 33.0 atol=1e-9
    @test_logs (:warn,) dat(2100, 1, 17, 0.0)
    @test_throws ERFAException dat(1e9, 1, 1, 0.0)
    @test_throws ERFAException dat(-1e9, 1, 1, 0.0)
    @test_throws ERFAException dat(2000, 13, 1, 0.0)
    @test_throws ERFAException dat(2000, 1, 32, 0.0)
    @test_throws ERFAException dat(2000, 1, 1, 3.0)
end

@testset "d2dtf" begin
    y, m, d, H, M, S, F = d2dtf("UTC", 5, 2400000.5, 49533.99999)
    @test (y, m, d, H, M, S, F) == (1994, 6, 30, 23, 59, 60, 13599)
    @test_logs (:warn,) d2dtf("UTC", 5, 0, 0)
    @test_throws ERFAException d2dtf("TAI", 5, 2e9, 0.0)
end

@testset "d2tf" begin
    @test d2tf(4, -0.987654321) == ('-', 23, 42, 13, 3333)
end

@testset "dtdb" begin
    d = dtdb(2448939.5, 0.123, 0.76543, 5.0123, 5525.242, 3190.0)
    @test d ≈ -0.1280368005936998991e-2 atol=1e-17
end

@testset "dtf2d" begin
    jd1, jd2 = dtf2d("UTC", 1994, 6, 30, 23, 59, 60.13599)
    @test jd1 + jd2 ≈ 2449534.49999 atol=1e-6
    @test_logs (:warn,) dtf2d("UTC", 1950, 6, 30, 23, 59, 59.13599)
    @test_logs (:warn,) dtf2d("UTC", 2021, 6, 30, 23, 59, 60.13599)
    @test_logs (:warn,) dtf2d("UTC", 1950, 6, 30, 23, 59, 60.13599)
    @test_throws ERFAException dtf2d("TAI", -2e9, 6, 30, 23, 59, 59.0)
    @test_throws ERFAException dtf2d("TAI", 2021, 13, 30, 23, 59, 59.0)
    @test_throws ERFAException dtf2d("TAI", 2021, 12, 32, 23, 59, 59.0)
    @test_throws ERFAException dtf2d("TAI", 2021, 12, 31, 24, 59, 59.0)
    @test_throws ERFAException dtf2d("TAI", 2021, 12, 31, 23, 60, 59.0)
    @test_throws ERFAException dtf2d("TAI", 2021, 12, 31, 23, 59, -59.0)
end


text:
using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates

threads=(32,5)
blocks =8
mainArrDims= (516,523,826)
dataBdim = (43,21,17)
metaData =view(MetaDataUtils.allocateMetadata(mainArrDims,dataBdim),1:4,1:4,1:4,:);
#metaData = view(MetaDataUtils.allocateMetadata(mainArrDims,dataBdim),1:9,2:3,4:6,: );
metaDataDims=size(metaData)
iterThrougWarNumb = cld(14,threads[2])
resShmem = CuArray(falses(dataBdim[1]+2, dataBdim[2]+2, dataBdim[3]+2 ))
totalFp,totalFn = 100000,100000
resList,resListIndicies,maxResListIndex= allocateResultLists(totalFp,totalFn)

loopXMeta= fld(metaDataDims[1],threads[1])
loopYZMeta= fld(metaDataDims[2]*metaDataDims[3],blocks )

#we set some offsets to make it simple for evaluation we will keeep it  each separated by 50
offset = -49
for metaX in 1:3, metaY in 1:3, metaZ in 1:3
   for quueueNumb in 1:14
      offset+=350
      #set offset
      metaData[metaX,metaY,metaZ,getResOffsetsBeg()+quueueNumb]=offset
      #set counter
   end#for
end #for

#setting number of fp ... so it should be validated as fit to be validated
offset = -49
   for quueueNumb in 1:14
         offset+=350
         #... so it should be validated as fit to be validated
           metaData[2,2,2,getBeginingOfFpFNcounts()+quueueNumb]+=8 # so all will remain active
           metaData[2,2,2,getNewCountersBeg()+quueueNumb]+=6
           metaData[2,2,2,getOldCountersBeg()+quueueNumb]+=5
         # so it should not be validated
         metaData[3,3,3,getBeginingOfFpFNcounts()+quueueNumb]+=1 # so all will be inactive
         metaData[3,3,3,getNewCountersBeg()+quueueNumb]+=1

end#for

#we are simulating some results in some of the result queues
offset = -49
   for quueueNumb in 1:14
         offset+=350
      for j in 1:(quueueNumb*10)
            # if(quueueNumb ==2 )
            #    println("list $([innerJ+1,innerJ+2,innerJ+3,1,1,1]) ind $(offset+j+innerJ*quueueNumb)   ")
            # end
        #the bigger the number the more repetitions and more non repeated elements
           resList[offset+(j+1),:]= [mod(j,quueueNumb*5)+1,mod(j,quueueNumb*5)+2,mod(j,quueueNumb*5)+3,1,1,1,1]

           metaData[1,1,1,getNewCountersBeg()+quueueNumb]+=1
           metaData[1,1,1,getBeginingOfFpFNcounts()+quueueNumb]+=2 # so all will remain active
   end#for
end#for

#should be first queue in block 3,3,3
resList[127701+1,:] = [1,1,1,1,1,1,1]
resList[127701+2,:] = [1,1,1,0,1,1,1]
resList[127701+3,:] = [1,1,1,1,1,1,1]#repeat
resList[127701+4,:] = [1,2,1,1,1,1,1]
resList[127701+5,:] = [1,1,2,1,1,1,1]
resList[127701+6,:] = [2,1,1,1,1,1,1]

metaData[3,3,3,getNewCountersBeg()+1]+=6
####### check is test written well (testing the test)
@test Int64(metaData[3,3,3,getNewCountersBeg()+1]) ==7

   #setting res list indicies
   for i in 1:(totalFp+totalFn)
      if(resList[i,1]!=0)
         resListIndicies[i] = UInt32(getResLinIndex(resList[i,1],resList[i,2],resList[i,3],resList[i,4],mainArrDims) )
      end
   end

 offsetForTwelve= Int64(metaData[1,1,1,getResOffsetsBeg()+12])
 notShuffl= resListIndicies[(offsetForTwelve+1):(offsetForTwelve+40)]= shuffle( resListIndicies[(offsetForTwelve+1):(offsetForTwelve+40)])

 offset = -49
for quueueNumb in 1:14
      offset+=350
     #the bigger the number the more repetitions and more non repeated elements
      @test Int64.(metaData[1,1,1,getNewCountersBeg()+quueueNumb])==quueueNumb*10
      #we test  entries
      tempList= resListIndicies[offset:offset+350]
      # println("**************************  $(quueueNumb)")
      # println(Int64.(filter(el->el>0,tempList)))
      # println(Int64.(unique(tempList)))

      @test length(filter(el->el>0,tempList))==quueueNumb*10
      @test length(unique(tempList))==quueueNumb*5+1

end#for quueueNumb
maxResListIndex= length(resListIndicies)

globalCurrentFnCount,globalCurrentFpCount= CUDA.zeros(UInt32,1),CUDA.zeros(UInt32,1)
inBlockLoopXZIterWithPadding,shmemblockDataLoop,shmemblockDataLenght,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength = calculateLoopsIter(dataBdim,threads[1],threads[2],metaDataDims,blocks)

function loadAndSanForDuplKernel(shmemblockDataLoop,shmemblockDataLenght,loopWarpMeta,metaDataLength,globalCurrentFnCount,globalCurrentFpCount,maxResListIndex,resListIndicies,metaData,iterThrougWarNumb,mainArrDims ,metaDataDims,loopXMeta,loopYZMeta,resList,dataBdim)
   locArr= UInt32(0)
   localOffset= UInt32(0)
   offsetIter= UInt16(0)
   shmemSum =  @cuStaticSharedMem(UInt32,(36,14)) # we need this additional spots
   alreadyCoveredInQueues =@cuStaticSharedMem(UInt32,(15))
   shmemblockData = @cuStaticSharedMem(UInt32,(37, 32))

   sync_threads()

   @iterateLinearly shmemblockDataLoop shmemblockDataLenght shmemblockData[i]=0

   for i in 1:14
      shmemSum[(threadIdxX()),i]= 0
   end

   MetadataAnalyzePass.@metaDataWarpIter( metaDataDims,loopWarpMeta,metaDataLength,
       begin

        # @exOnWarp 16 CUDA.@cuprint "idX $(threadIdxX())  xMeta $(xMeta) yMeta $(yMeta+1) zMeta $(zMeta+1)  \n "

        @loadAndScanForDuplicates(iterThrougWarNumb,locArr,offsetIter,localOffset)
sync_threads()

 for i in 1:14
    shmemSum[(threadIdxX()),i]= 0
    shmemSum[33,i]= 0
    shmemSum[34,i]= 0
    shmemSum[35,i]= 0
    shmemSum[36,i]= 0
 end
# for i in 1:14
#   @ifX i alreadyCoveredInQueues[i]=0
# end

 sync_threads()
       end)

       sync_threads()
       @ifXY 1 1 begin
         if(blockIdxX()==1)
         summA=0
         suumB=0
         for i in 1:12
            summA+=alreadyCoveredInQueues[i]
            suumB+=i*5
         end
         for i in 13:14
            summA+=alreadyCoveredInQueues[i]
            suumB+=i*10
         end
         # CUDA.@cuprint " alreadyCoveredInQueues[1] $(alreadyCoveredInQueues[1])   should be $(1*5)\n "
         # CUDA.@cuprint " alreadyCoveredInQueues[2] $(alreadyCoveredInQueues[2]) should be $(2*5) \n "
         # CUDA.@cuprint " alreadyCoveredInQueues[3] $(alreadyCoveredInQueues[3]) should be $(3*5) \n "
         # CUDA.@cuprint " alreadyCoveredInQueues[4] $(alreadyCoveredInQueues[4]) should be $(4*5) \n "
         # CUDA.@cuprint " alreadyCoveredInQueues[5] $(alreadyCoveredInQueues[5]) should be $(5*5) \n "
         # CUDA.@cuprint " alreadyCoveredInQueues[6] $(alreadyCoveredInQueues[6]) should be $(6*5) \n "
         # CUDA.@cuprint " alreadyCoveredInQueues[7] $(alreadyCoveredInQueues[7]) should be $(7*5) \n "
         # CUDA.@cuprint " alreadyCoveredInQueues[8] $(alreadyCoveredInQueues[8]) should be $(8*5) \n "
         # CUDA.@cuprint " alreadyCoveredInQueues[9] $(alreadyCoveredInQueues[9]) should be $(9*5) \n "
         # CUDA.@cuprint " alreadyCoveredInQueues[10] $(alreadyCoveredInQueues[10]) should be $(10*5) \n "
         # CUDA.@cuprint " alreadyCoveredInQueues[11] $(alreadyCoveredInQueues[11]) should be $(11*5) \n "
         # CUDA.@cuprint " alreadyCoveredInQueues[12] $(alreadyCoveredInQueues[12]) should be $(12*5) \n "
         # CUDA.@cuprint " alreadyCoveredInQueues[13] $(alreadyCoveredInQueues[13]) should be $(13*10) \n "
         # CUDA.@cuprint " alreadyCoveredInQueues[14] $(alreadyCoveredInQueues[14]) should be $(14*10) \n "
         # CUDA.@cuprint "summA  $(summA) should be  $(suumB) \n "

      end
         #   CUDA.@cuprint """  valuee fp $(alreadyCoveredInQueues[1]+ alreadyCoveredInQueues[3]+ alreadyCoveredInQueues[5]+ alreadyCoveredInQueues[7]+ alreadyCoveredInQueues[9]+ alreadyCoveredInQueues[11]+ alreadyCoveredInQueues[13])
         #      alreadyCoveredInQueues[1] $(alreadyCoveredInQueues[1]) alreadyCoveredInQueues[3] $(alreadyCoveredInQueues[3]) alreadyCoveredInQueues[5] $(alreadyCoveredInQueues[5]) alreadyCoveredInQueues[7] $(alreadyCoveredInQueues[7]) alreadyCoveredInQueues[9] $(alreadyCoveredInQueues[9]) alreadyCoveredInQueues[11] $(alreadyCoveredInQueues[11]) alreadyCoveredInQueues[13] $(alreadyCoveredInQueues[13])

         #      """
          atomicAdd(globalCurrentFpCount, alreadyCoveredInQueues[1]+ alreadyCoveredInQueues[3]+ alreadyCoveredInQueues[5]+ alreadyCoveredInQueues[7]+ alreadyCoveredInQueues[9]+ alreadyCoveredInQueues[11]+ alreadyCoveredInQueues[13])
       end
          @ifXY 2 1 atomicAdd(globalCurrentFnCount, alreadyCoveredInQueues[2]+ alreadyCoveredInQueues[4]+ alreadyCoveredInQueues[6]+ alreadyCoveredInQueues[8]+ alreadyCoveredInQueues[10]+ alreadyCoveredInQueues[12]+ alreadyCoveredInQueues[14])

       sync_threads()
   return
end

numbOfFp = 0
numbOfFn = 0

for i in [1,3,5,7,9,11]
   numbOfFp+=i*5
end
numbOfFp+=13*10
for i in [2,4,6,8,10,12]
   numbOfFn+=i*5
end
numbOfFn+=14*10
numbOfFp+=7*2
numbOfFn+=7*2

numbOfFp+=5

@cuda threads=threads blocks=blocks loadAndSanForDuplKernel(shmemblockDataLoop,shmemblockDataLenght,loopWarpMeta,metaDataLength,globalCurrentFnCount,globalCurrentFpCount,maxResListIndex,resListIndicies,metaData,iterThrougWarNumb,mainArrDims ,metaDataDims,loopXMeta,loopYZMeta,resList,dataBdim)

Int64(globalCurrentFnCount[1]+globalCurrentFpCount[1])== 660+5
offset = -49
   # for quueueNumb in 1:14
   for quueueNumb in 1:12
         offset+=350
        #the bigger the number the more repetitions and more non repeated elements
         @test Int64(metaData[1,1,1,getNewCountersBeg()+quueueNumb])==quueueNumb*5#j instead of quueueNumb*quueueNumb
         tempList= Array(resListIndicies[offset:offset+350])
         @test length(filter(el->el>0,tempList))==quueueNumb*5#+1 becouse of 0's entry
  end#for

@testset "scanning and checking is to be validated" begin
    quueueNumb =12
   offset = -49
   offset+=350*(quueueNumb)
   @test Int64(metaData[1,1,1,getNewCountersBeg()+quueueNumb])==quueueNumb*5#j instead of quueueNumb*quueueNumb
   tempList= Array(resListIndicies)[offset:offset+350];
   @test length(filter(el->el>0,tempList))==quueueNumb*5#+1 becouse of 0's entry
   @test Int64(metaData[3,3,3,getNewCountersBeg()+1]) ==6

@test Int64(globalCurrentFnCount[1]) ==numbOfFn
@test Int64(globalCurrentFpCount[1]) ==numbOfFp

   quueueNumb =2
   offset = -49
   offset+=350*(quueueNumb)
   @test Int64(metaData[1,1,1,getNewCountersBeg()+quueueNumb])==quueueNumb*5#j instead of quueueNumb*quueueNumb
   tempList= Array(resListIndicies)[offset:offset+350];
   @test length(filter(el->el>0,tempList))==quueueNumb*5#+1 becouse of 0's entry
   Int64.(filter(el->el>0,tempList))
   Int64.(unique(filter(el->el>0,tempList)))

  ############### testing is to be analyzed
  #main block
  @test metaData[1,1,1,getIsToBeAnalyzedNumb()+15 ]==1
  @test metaData[1,1,1,getIsToBeAnalyzedNumb()+16 ]==1
  #right
  @test metaData[2,1,1,getIsToBeAnalyzedNumb()+3 ]==1
  @test metaData[2,1,1,getIsToBeAnalyzedNumb()+4 ]==1
  #bottom
  @test metaData[1,1,2,getIsToBeAnalyzedNumb()+11 ]==1
  @test metaData[1,1,2,getIsToBeAnalyzedNumb()+12 ]==1
  #anterior
  @test metaData[1,2,1,getIsToBeAnalyzedNumb()+7 ]==1
  @test metaData[1,2,1,getIsToBeAnalyzedNumb()+8 ]==1

  ### for block  2,2,2 all queaueas are active  and should activate surrounding blocks

  #main block
  metaX = 2
  metaY = 2
  metaZ = 2

  @test metaData[metaX,metaY,metaZ,getIsToBeAnalyzedNumb()+15 ]==1
  @test metaData[metaX,metaY,metaZ,getIsToBeAnalyzedNumb()+16 ]==1

  @test metaData[metaX-1,metaY,metaZ,getIsToBeAnalyzedNumb()+1 ]==1
  @test metaData[metaX-1,metaY,metaZ,getIsToBeAnalyzedNumb()+2 ]==1
  @test metaData[metaX+1,metaY,metaZ,getIsToBeAnalyzedNumb()+3 ]==1
  @test metaData[metaX+1,metaY,metaZ,getIsToBeAnalyzedNumb()+4 ]==1
  @test metaData[metaX,metaY-1,metaZ,getIsToBeAnalyzedNumb()+5 ]==1
  @test metaData[metaX,metaY-1,metaZ,getIsToBeAnalyzedNumb()+6 ]==1
  @test metaData[metaX,metaY+1,metaZ,getIsToBeAnalyzedNumb()+7 ]==1
  @test metaData[metaX,metaY+1,metaZ,getIsToBeAnalyzedNumb()+8 ]==1
  @test metaData[metaX,metaY,metaZ-1,getIsToBeAnalyzedNumb()+9 ]==1
  @test metaData[metaX,metaY,metaZ-1,getIsToBeAnalyzedNumb()+10 ]==1
  @test metaData[metaX,metaY,metaZ+1,getIsToBeAnalyzedNumb()+11 ]==1
  @test metaData[metaX,metaY,metaZ+1,getIsToBeAnalyzedNumb()+12 ]==1

  metaX = 3
  metaY = 3
  metaZ = 3

  @test metaData[metaX,metaY,metaZ,getIsToBeAnalyzedNumb()+15 ]==1
  @test metaData[metaX,metaY,metaZ,getIsToBeAnalyzedNumb()+16 ]==0

  @test metaData[metaX-1,metaY,metaZ,getIsToBeAnalyzedNumb()+1 ]==1
  @test metaData[metaX-1,metaY,metaZ,getIsToBeAnalyzedNumb()+2 ]==0
  @test metaData[metaX+1,metaY,metaZ,getIsToBeAnalyzedNumb()+3 ]==0
  @test metaData[metaX+1,metaY,metaZ,getIsToBeAnalyzedNumb()+4 ]==0
  @test metaData[metaX,metaY-1,metaZ,getIsToBeAnalyzedNumb()+5 ]==0
  @test metaData[metaX,metaY-1,metaZ,getIsToBeAnalyzedNumb()+6 ]==0
  @test metaData[metaX,metaY+1,metaZ,getIsToBeAnalyzedNumb()+7 ]==0
  @test metaData[metaX,metaY+1,metaZ,getIsToBeAnalyzedNumb()+8 ]==0
  @test metaData[metaX,metaY,metaZ-1,getIsToBeAnalyzedNumb()+9 ]==0
  @test metaData[metaX,metaY,metaZ-1,getIsToBeAnalyzedNumb()+10 ]==0
  @test metaData[metaX,metaY,metaZ+1,getIsToBeAnalyzedNumb()+11 ]==0
  @test metaData[metaX,metaY,metaZ+1,getIsToBeAnalyzedNumb()+12 ]==0

end

  ### for block 3,3,3 only queaue 1 should be active  and rest should not activate surrounding blocks

#mixed entries ...

#    zzzzzzzzzzz    shmemSum[36,innerWarpNumb]  4  innerWarpNumb 1    idX 1  idY 1  xMeta 1 yMeta 0 zMeta 0
#    zzzzzzzzzzz    shmemSum[36,innerWarpNumb]  9  innerWarpNumb 2    idX 1  idY 2  xMeta 1 yMeta 0 zMeta 0
#    zzzzzzzzzzz    shmemSum[36,innerWarpNumb]  14  innerWarpNumb 3    idX 1  idY 3  xMeta 1 yMeta 0 zMeta 0
#    zzzzzzzzzzz    shmemSum[36,innerWarpNumb]  47  innerWarpNumb 4    idX 1  idY 4  xMeta 1 yMeta 0 zMeta 0
#    zzzzzzzzzzz    shmemSum[36,innerWarpNumb]  177  innerWarpNumb 5    idX 1  idY 5  xMeta 1 yMeta 0 zMeta 0
#    zzzzzzzzzzz    shmemSum[36,innerWarpNumb]  1  innerWarpNumb 1    idX 3  idY 1  xMeta 3 yMeta 2 zMeta 2

#    offset = -49+(350*8)
#    Int64.(resList[offset:offset+10,:])

#    indd 1102501  xMeta 1 yMeta 1 zMeta 1 innerWarpNumb 2 tempCount 656
#    indd 1102501  xMeta 1 yMeta 1 zMeta 1 innerWarpNumb 2 tempCount 653

#    Int64.(unique(resListIndicies[offset:offset+350]))
#    Int64(resListIndicies[656 ]) #0
#    Int64(resListIndicies[653 ]) #223716...
#    Int64(resListIndicies[652 ]) #223716...
#    Int64(metaData[1102501]) #4

#    list [1, 2, 3, 1, 1, 1] ind 652
#    list [2, 3, 4, 1, 1, 1] ind 654
#    list [1, 2, 3, 1, 1, 1] ind 653
#    list [2, 3, 4, 1, 1, 1] ind 655

# Int64(resListIndicies[127701] )
# Int64(resListIndicies[127701+1] )
# Int64(resListIndicies[127701+2] )
# Int64(resListIndicies[127701+3] )
# Int64(resListIndicies[127701+4] )
# Int64(resListIndicies[127701+5] )

# Int64(resListIndicies[127701] ) ==223181869
# Int64(resListIndicies[127701+1] ) ==223181869
# Int64(resListIndicies[127701+2] ) ==223181869
# Int64(resListIndicies[127701+3] ) ==223181869
# Int64(resListIndicies[127701+4] ) ==223181869
# Int64(resListIndicies[127701+5] ) ==223181869

# 223451221

# 223181869 -223181869

# for i in  0:0
#    println("111")
# end

# using BenchmarkTools
# BenchmarkTools.DEFAULT_PARAMETERS.samples = 100
# BenchmarkTools.DEFAULT_PARAMETERS.seconds =60
# BenchmarkTools.DEFAULT_PARAMETERS.gcsample = true

# function toBench()
#    CUDA.@sync @cuda threads=threads blocks=blocks loadAndSanForDuplKernel(resListIndicies,metaData,iterThrougWarNumb,mainArrDims ,metaDataDims,loopXMeta,loopYZMeta,resList,dataBdim)
# end
# bb2 = @benchmark toBench()


text:
import FireflyModule.Firefly

function light_intensity(source_node::Firefly, destination_node::Firefly, λ::Float)::Float
    """Calculates the light intensity for the source firefly.
    """
    I_0 = source_node.cost
    r, _ = hamming_distance(source_node, destination_node)
    return I_0 * ℯ ^ (-λ * (r^2))
end

max_stars_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_issues_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_forks_repo_name: zahidkizmaz/firefly-algorithm-tsp


text:
module JLD2
using DataStructures
import Base.sizeof
export jldopen

const OBJECT_HEADER_SIGNATURE = reinterpret(UInt32, UInt8['O', 'H', 'D', 'R'])[1]

# Currently we specify that all offsets and lengths are 8 bytes
const Length = UInt64

# Currently we specify a 512 byte header
const FILE_HEADER_LENGTH = 512
const FILE_HEADER = "Julia data file (HDF5), version "
const CURRENT_VERSION = v"0.2"

struct UnsupportedVersionException <: Exception end
struct UnsupportedFeatureException <: Exception end
struct InvalidDataException <: Exception end

include("Lookup3.jl")
include("mmapio.jl")
include("misc.jl")

# RelOffset represents an HDF5 relative offset. It differs from a file offset (used
# elsewhere) in that it is relative to the superblock base address. In practice,
# this means that FILE_HEADER_LENGTH has been subtracted. `fileoffset` and
# `h5offset` convert between RelOffsets and file offsets
struct RelOffset
    offset::UInt64
end
define_packed(RelOffset)
Base.:(==)(x::RelOffset, y::RelOffset) = x === y
Base.hash(x::RelOffset) = hash(x.offset)

const UNDEFINED_ADDRESS = RelOffset(0xffffffffffffffff)
const NULL_REFERENCE = RelOffset(0)

struct JLDWriteSession{T<:Union{Dict{UInt,RelOffset},Union{}}}
    h5offset::T
    objects::Vector{Any}

    JLDWriteSession{T}() where T = new()
    JLDWriteSession{T}(h5offset, objects) where T = new(h5offset, objects)
end
JLDWriteSession() = JLDWriteSession{Dict{UInt,RelOffset}}(Dict{UInt,RelOffset}(), Any[])

mutable struct GlobalHeap
    offset::Int64
    length::Length
    free::Length
    objects::Vector{Int64}
end

abstract type H5Datatype end

struct CommittedDatatype <: H5Datatype
    header_offset::RelOffset
    index::Int
end

struct ReadRepresentation{T,ODR} end
struct CustomSerialization{T,S} end

symbol_length(x::Symbol) = ccall(:strlen, Int, (Cstring,), x)

mutable struct JLDFile{T<:IO}
    io::T
    writable::Bool
    written::Bool
    created::Bool
    datatype_locations::OrderedDict{RelOffset,CommittedDatatype}
    datatypes::Vector{H5Datatype}
    datatype_wsession::JLDWriteSession{Dict{UInt,RelOffset}}
    datasets::OrderedDict{String,RelOffset}
    jlh5type::ObjectIdDict
    h5jltype::ObjectIdDict
    jloffset::Dict{RelOffset,WeakRef}
    end_of_data::Int64
    global_heaps::Dict{RelOffset,GlobalHeap}
    global_heap::GlobalHeap
end
JLDFile(io::IO, writable::Bool, written::Bool, created::Bool) =
    JLDFile(io, writable, written, created, OrderedDict{RelOffset,CommittedDatatype}(), H5Datatype[],
            JLDWriteSession(), OrderedDict{String,RelOffset}(), ObjectIdDict(),
            ObjectIdDict(), Dict{RelOffset,WeakRef}(),
            Int64(FILE_HEADER_LENGTH + sizeof(Superblock)), Dict{RelOffset,GlobalHeap}(),
            GlobalHeap(0, 0, 0, Int64[]))

fileoffset(f::JLDFile, x::RelOffset) = Int64(x.offset + FILE_HEADER_LENGTH)
h5offset(f::JLDFile, x::Int64) = RelOffset(x - FILE_HEADER_LENGTH)

#
# File
#

function jldopen(fname::AbstractString, wr::Bool, create::Bool, truncate::Bool)
    exists = isfile(fname)
    io = MmapIO(fname, wr, create, truncate)
    f = JLDFile(io, wr, truncate, !exists || truncate)

    if !truncate
        if String(read(io, UInt8, length(FILE_HEADER))) != FILE_HEADER
            throw(ArgumentError(string('"', fname, "\" is not a JLD file")))
        end

        ver = convert(VersionNumber, read_bytestring(io))
        if ver < v"0.2"
            throw(ArgumentError("only JLD2 files are presently supported"))
        elseif ver > CURRENT_VERSION
            warn('"', fname, "\" was written in JLD file format version ", ver,
                 ", but this version of JLD supports only JLD file format ", CURRENT_VERSION,
                 ". Some or all data in the file may not be readable")
        end
    end

    if wr
        seek(io, 0)
        # Write JLD header
        write(io, FILE_HEADER)
        print(io, CURRENT_VERSION)
    end

    if !truncate
        seek(io, FILE_HEADER_LENGTH)
        superblock = read(io, Superblock)
        f.end_of_data = superblock.end_of_file_address

        root_group_offset = fileoffset(f, superblock.root_group_object_header_address)
        seek(io, root_group_offset)
        root_group = read(io, Group)
        if position(io) == f.end_of_data
            f.end_of_data = root_group_offset
        end

        for i = 1:length(root_group.names)
            name = root_group.names[i]
            offset = root_group.offsets[i]
            if name == "_types"
                types_group_offset = fileoffset(f, offset)
                seek(io, types_group_offset)
                types_group = read(io, Group)
                if position(io) == f.end_of_data
                    f.end_of_data = types_group_offset
                end

                for i = 1:length(types_group.offsets)
                    f.datatype_locations[types_group.offsets[i]] = CommittedDatatype(types_group.offsets[i], i)
                end
                resize!(f.datatypes, length(types_group.offsets))
            else
                f.datasets[name] = offset
            end
        end
    end

    f
end

function jldopen(fname::AbstractString, mode::AbstractString="r")
    mode == "r"  ? jldopen(fname, false, false, false) :
    mode == "r+" ? jldopen(fname, true, false, false) :
    mode == "a" || mode == "a+" ? jldopen(fname, true, true, false) :
    mode == "w" || mode == "w+" ? jldopen(fname, true, true, true) :
    throw(ArgumentError("invalid open mode: $mode"))
end

function Base.read(f::JLDFile, name::AbstractString)
    f.end_of_data == 0 && throw(ArgumentError("file is closed"))
    haskey(f.datasets, name) || throw(ArgumentError("file has no dataset $name"))
    read_dataset(f, f.datasets[name])
end

# Populate f.datatypes and f.jlh5types with all of the committed datatypes from
# a file. We need to do this before writing to make sure we reuse written
# datatypes.
function load_datatypes(f::JLDFile)
    dts = f.datatypes
    cdts = f.datatype_locations
    @assert length(dts) == length(cdts)
    i = 1
    for cdt in values(cdts)
        !isassigned(dts, i) && jltype(f, cdt)
        i += 1
    end
end

function Base.write(f::JLDFile, name::AbstractString, obj, wsession::JLDWriteSession=JLDWriteSession())
    f.end_of_data == 0 && throw(ArgumentError("file is closed"))
    !f.writable && throw(ArgumentError("file was opened read-only"))
    !f.written && !f.created && load_datatypes(f)
    f.written = true

    io = f.io
    seek(io, f.end_of_data)
    header_offset = write_dataset(f, obj, wsession)
    f.datasets[name] = header_offset
    nothing
end

function Base.close(f::JLDFile)
    io = f.io
    if f.written
        seek(io, f.end_of_data)

        names = String[]
        sizehint!(names, length(f.datasets)+1)
        offsets = RelOffset[]
        sizehint!(offsets, length(f.datasets)+1)

        # Write types group
        if !isempty(f.datatypes)
            push!(names, "_types")
            push!(offsets, h5offset(f, position(io)))
            write(io, Group(String[@sprintf("%08d", i) for i = 1:length(f.datatypes)],
                            collect(keys(f.datatype_locations))))
        end

        # Write root group
        root_group_object_header_address = h5offset(f, position(io))
        for (k, v) in f.datasets
            push!(names, k)
            push!(offsets, v)
        end
        write(io, Group(names, offsets))

        eof_position = position(io)
        truncate(io, eof_position)
        seek(io, FILE_HEADER_LENGTH)
        write(io, Superblock(0, FILE_HEADER_LENGTH, UNDEFINED_ADDRESS,
              eof_position, root_group_object_header_address))
    end
    f.end_of_data = 0
    close(io)
    nothing
end

include("superblock.jl")
include("object_headers.jl")
include("groups.jl")
include("dataspaces.jl")
include("attributes.jl")
include("datatypes.jl")
include("datasets.jl")
include("global_heaps.jl")
include("data.jl")

end # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule wannier90_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("wannier90")
JLLWrappers.@generate_main_file("wannier90", UUID("c5400fa0-8d08-52c2-913f-1e3f656c1ce9"))
end  # module wannier90_jll


text:
module BenchmarkToolsMPI
"""
Make BenchmarkTools MPI-aware

Overrides a few functions from BenchmarkTools to avoid hanging when running with MPI.
"""

using Reexport
@reexport using BenchmarkTools
using MPI

import BenchmarkTools
using BenchmarkTools: Benchmark, gcscrub, Parameters, RESOLUTION, Trial

function __init__()
    MPI.Init()
end

function BenchmarkTools._lineartrial(b::Benchmark, p::Parameters = b.params;
                                     maxevals = RESOLUTION, kwargs...)
    params = Parameters(p; kwargs...)
    estimates = zeros(maxevals)
    completed = 0
    params.gctrial && gcscrub()
    start_time = time()
    for evals in eachindex(estimates)
        params.gcsample && gcscrub()
        params.evals = evals
        estimates[evals] = first(b.samplefunc(params))
        completed += 1
        run_time = MPI.bcast(time() - start_time, 0, MPI.COMM_WORLD)
        (run_time > params.seconds) && break
    end
    return estimates[1:completed]
end

function BenchmarkTools._run(b::Benchmark, p::Parameters;
                             verbose = false, pad = "", kwargs...)
    params = Parameters(p; kwargs...)
    @assert params.seconds > 0.0 "time limit must be greater than 0.0"
    params.gctrial && gcscrub()
    start_time = Base.time()
    trial = Trial(params)
    params.gcsample && gcscrub()
    s = b.samplefunc(params)
    push!(trial, s[1:end-1]...)
    return_val = s[end]
    iters = 2
    run_time = start_time
    while (run_time - start_time) < params.seconds && iters ≤ params.samples
        params.gcsample && gcscrub()
        push!(trial, b.samplefunc(params)[1:end-1]...)
        iters += 1
        run_time = MPI.bcast(time() - start_time, 0, MPI.COMM_WORLD)
    end
    return trial, return_val
end

end # BenchmarkToolsMPI


text:
"""
    sparsemax(x; dims = 1)

[Sparsemax](https://arxiv.org/abs/1602.02068) turns input array `x`
into sparse probability distributions that sum to 1 along the dimensions specified by `dims`.

Similar to softmax, each dimension is considered independent. For a matrix input `x` it will
by default (`dims = 1`) treat it as a batch of vectors, with each column independent. Keyword
`dims = 2` will instead treat rows independently.

# Examples

```jldoctest
julia> sparsemax(reshape(1:15, 5,3) ./ 2) # dims=1
5×3 Matrix{Float64}:
 0.0   0.0   0.0
 0.0   0.0   0.0
 0.0   0.0   0.0
 0.25  0.25  0.25
 0.75  0.75  0.75

julia> sparsemax(reshape(1:15, 5, 3) ./ 20) # Scale matters
5×3 Matrix{Float64}:
 0.1   0.1   0.1
 0.15  0.15  0.15
 0.2   0.2   0.2
 0.25  0.25  0.25
 0.3   0.3   0.3

 julia> sparsemax(reshape(1:15, 5, 3) ./ 20 .+ [-20 0 20])
 5×3 Matrix{Float64}:
  0.1   0.1   0.1
  0.15  0.15  0.15
  0.2   0.2   0.2
  0.25  0.25  0.25
  0.3   0.3   0.3
```
"""
function sparsemax(x::AbstractArray; dims::Integer=1)
    n_dims = length(size(x))
    if n_dims < dims
        error_msg = "Dims ($dims) is larger than the array dimensions ($n_dims)."
        throw(DimensionMismatch(error_msg))
    end

    if isa(x, AbstractVector)
        z = _sort_vector(x; dims=dims)
    else
        z = sort(float(x); dims=dims, rev=true)
    end
    mask = _sparsemax_mask(z, dims)
    tausum = sum(z .* mask; dims)
    kay = sum(mask; dims)
    z = _relu.(x  .- (tausum .- 1) ./ kay)
end

function _sort_vector(x::AbstractVector; dims::Integer=1, rev=true)
    if dims == 1
        z = sort(float(x); rev=rev)
    elseif dims == 2
        z = x
    end
    return z
end

function _sparsemax_mask(z::AbstractArray, dim::Integer)
    acc = cumsum(z; dims=dim)
    if dim == 1
        acc .= 1 .+ axes(z,1) .* z .> acc
    elseif dim == 2
        acc .= 1 .+ axes(z,2)' .* z .> acc
    else
        # This isn't type-stable. Writing into `acc` ensures the whole function still is:
        cnt = reshape(axes(x, dim), ntuple(_->1, dim-1)..., :)
        acc .= 1 .+ cnt .* z .> acc
    end
    acc
end

_ifelse(p, x, y) = ifelse(p, promote(x, y)...)
_relu(x) = _ifelse(x>0, x, false)  # different gradient at zero

function ∇sparsemax(Δ::AbstractArray, y::AbstractArray; dims = 1)
    nonzeros = Δ.!=0.0
    out = Δ .* y
    total = sum(out .* nonzeros; dims=dims) / sum(nonzeros; dims=dims)
    out = nonzeros .* (out .- total)
 end

function rrule(::typeof(sparsemax), xs; dims=1)
    y = sparsemax(xs; dims=dims)
    sparsemax_pullback(Δ) = (NoTangent(), ∇sparsemax(unthunk(Δ), y; dims = dims))
    return y, sparsemax_pullback
 end


text:
"""
# Usage

    solve(FDProb::FractionalDifferenceProblem, T, h, PECEDifference())

Use the PECE algorithm to solve fractional difference equations

### References

```tex
@article{陈福来2011分数阶微分差分方程的,
  title={分数阶微分差分方程的 Matlab 求解},
  author={陈福来 and 王华 and 李势丰},
  journal={湘南学院学报},
  volume={32},
  number={5},
  pages={1--4},
  year={2011}
}
```
"""
struct PECEDifference <: FractionalDiffEqAlgorithm end

function solve(FDProb::FractionalDifferenceProblem, T, h, ::PECEDifference)
    @unpack fun, α, u0 = FDProb
    N = round(Int, T/h)

    # Initialize
    b=zeros(N); a=zeros(N); x=zeros(N); t=zeros(N)
    t1=zeros(N+1); y=zeros(N+1)

    @fastmath @inbounds @simd for i=1:N
        k = 1
        @fastmath @inbounds @simd for j=1:N-i
            k = k*j
        end
        b[i] = gamma(N-i+α)/k
        a[N+1-i] = b[i]
    end
    x[1] = u0+a[1]*fun(u0)/gamma(α)
    t[1] = 1

    @fastmath @inbounds @simd for i=2:N
        temp = 0
        @fastmath @inbounds @simd for j=1:i-1
            temp = a[i+1-j]*fun(x[j])
        end
        temp = temp/gamma(α)
        x[i] = u0 + temp+a[1]*(u0+temp)/gamma(α)
        t[i] = i
    end

    y[1] = u0
    t1[1] = 0
    @fastmath @inbounds @simd for i = 1:N
        y[i+1] = x[i]
        t1[i+1] = t[i]
    end
    return FDifferenceSolution(t1, y)
end


text:
# This file contains functions for creating spatial weights from points data using a k nearest neighbors

"""
    knearneigh(X, Y; k)
Build a spatial weights object from a set of coordinates using `k` nearest neighbors.
"""
function knearneigh(X::Vector{Float64}, Y::Vector{Float64}; k::Int)::SpatialWeights

    n, ny = length(X), length(Y)

    n == ny || throw(DimensionMismatch("dimensions must match: X has length ($(n)), Y has length ($ny)"))

    # Get nearest neighbors using a KDTree
    cpoints = vcat(X', Y')

    kdtree = KDTree(cpoints)
    idxs,  = knn(kdtree, cpoints, k + 1, true)

    # Check neighbours
    neighs = copy.(fill(Int[], n))
    weights = copy.(fill(Float64[], n))
    nneighs = zeros(Int,n)

    for i in 1:n
        neighsi = idxs[i]
        neighs[i] = sort!(neighsi[neighsi .!= i])
        nneighs[i] = length(neighs[i])
        weights[i] = ones(nneighs[i]) ./ nneighs[i]
    end

    SpatialWeights(n, neighs, weights, nneighs, :row)

end

"""
    knearneigh(P; k)
Build a spatial weights object from a vector of points using ``k`` nearest neighbors.
"""
function knearneigh(P::Vector{T} where T<:Union{Missing,AbstractPoint}; k::Int)::SpatialWeights

    cpoints = reduce(hcat, map(a -> coordinates(a), P))

    knearneigh(cpoints[1,:], cpoints[2,:]; k)
end

"""
    knearneigh(A; k)
Build a spatial weights object from a table A that constains a points geometry column using ``k`` nearest neighbors..
"""
function knearneigh(A::Any; k::Int)::SpatialWeights
    istable(A) || throw(ArgumentError("Unknown geometry or not points geometry"))

    (:geometry in propertynames(A)) || throw(ArgumentError("table does not have :geometry information"))

    knearneigh(A.geometry; k)
end