text:
using GPUifyLoops: @launch, @loop, @unroll
using Oceananigans.Utils: @loop_xy, @loop_xz, @loop_yz, launch_config

#####
##### Halo filling for value and gradient boundary conditions
#####

  @inline fill_west_halo!(c, bc::Union{VBC, GBC}, arch, grid, args...) = @launch device(arch) config=launch_config(grid, :yz)   _fill_west_halo!(c, bc, grid, args...)
  @inline fill_east_halo!(c, bc::Union{VBC, GBC}, arch, grid, args...) = @launch device(arch) config=launch_config(grid, :yz)   _fill_east_halo!(c, bc, grid, args...)
 @inline fill_south_halo!(c, bc::Union{VBC, GBC}, arch, grid, args...) = @launch device(arch) config=launch_config(grid, :xz)  _fill_south_halo!(c, bc, grid, args...)
 @inline fill_north_halo!(c, bc::Union{VBC, GBC}, arch, grid, args...) = @launch device(arch) config=launch_config(grid, :xz)  _fill_north_halo!(c, bc, grid, args...)
@inline fill_bottom_halo!(c, bc::Union{VBC, GBC}, arch, grid, args...) = @launch device(arch) config=launch_config(grid, :xy) _fill_bottom_halo!(c, bc, grid, args...)
   @inline fill_top_halo!(c, bc::Union{VBC, GBC}, arch, grid, args...) = @launch device(arch) config=launch_config(grid, :xy)    _fill_top_halo!(c, bc, grid, args...)

@inline linearly_extrapolate(c₀, ∇c, Δ) = c₀ + ∇c * Δ

@inline  left_gradient(bc::GBC, c¹, Δ, i, j, args...) = getbc(bc, i, j, args...)
@inline right_gradient(bc::GBC, cᴺ, Δ, i, j, args...) = getbc(bc, i, j, args...)

@inline  left_gradient(bc::VBC, c¹, Δ, i, j, args...) = ( c¹ - getbc(bc, i, j, args...) ) / (Δ/2)
@inline right_gradient(bc::VBC, cᴺ, Δ, i, j, args...) = ( getbc(bc, i, j, args...) - cᴺ ) / (Δ/2)

function _fill_west_halo!(c, bc::Union{VBC, GBC}, grid, args...)

           #  ↑ x ↑  interior
           #  -----  interior face
    iᴵ = 1 #    *    interior cell
    iᴮ = 1 #  =====  western boundary
    iᴴ = 0 #    *    halo cell

    @loop_yz j k grid begin
        Δ = Δx(iᴮ, j, k, grid) # Δ between first interior and first west halo point, defined at cell face.
        @inbounds ∇c = left_gradient(bc, c[iᴵ, j, k], Δ, j, k, grid, args...)
        @inbounds c[iᴴ, j, k] = linearly_extrapolate(c[iᴵ, j, k], ∇c, -Δ) # extrapolate westward in -x direction.
    end

    return nothing
end

function _fill_east_halo!(c, bc::Union{VBC, GBC}, grid, args...)

                     #  ↑ x ↑
    iᴴ = grid.Nx + 1 #    *   halo cell
    iᴮ = grid.Nx + 1 #  ===== eastern boundary
    iᴵ = grid.Nx     #    *   interior cell
                     #  ----- interior face
                     #    ↓   interior

    @loop_yz j k grid begin
        Δ = Δx(iᴮ, j, k, grid) # Δ between last interior and first east halo point, defined at cell face.
        @inbounds ∇c = right_gradient(bc, c[iᴵ, j, k], Δ, j, k, grid, args...)
        @inbounds c[iᴴ, j, k] = linearly_extrapolate(c[iᴵ, j, k], ∇c, Δ) # extrapolate eastward in +x direction.
    end

    return nothing
end

function _fill_south_halo!(c, bc::Union{VBC, GBC}, grid, args...)

           #  ↑ y ↑  interior
           #  -----  interior face
    jᴵ = 1 #    *    interior cell
    jᴮ = 1 #  =====  southern boundary
    jᴴ = 0 #    *    halo cell

    @loop_xz i k grid begin
        Δ = Δy(i, jᴮ, k, grid) # Δ between first interior and first south halo point, defined at cell face.
        @inbounds ∇c = left_gradient(bc, c[i, jᴵ, k], Δ, i, k, grid, args...)
        @inbounds c[i, jᴴ, k] = linearly_extrapolate(c[i, jᴵ, k], ∇c, -Δ) # extrapolate southward in -y direction.
    end

    return nothing
end

function _fill_north_halo!(c, bc::Union{VBC, GBC}, grid, args...)

                     #  ↑ y ↑
    jᴴ = grid.Ny + 1 #    *   halo cell
    jᴮ = grid.Ny + 1 #  ===== northern boundary
    jᴵ = grid.Ny     #    *   interior cell
                     #  ----- interior face
                     #    ↓   interior

    @loop_xz i k grid begin
        Δ = Δy(i, jᴮ, k, grid) # Δ between first interior and first north halo point, defined at cell face.
        @inbounds ∇c = right_gradient(bc, c[i, jᴵ, k], Δ, i, k, grid, args...)
        @inbounds c[i, jᴴ, k] = linearly_extrapolate(c[i, jᴵ, k], ∇c, Δ) # extrapolate northward in +y direction.
    end

    return nothing
end

function _fill_bottom_halo!(c, bc::Union{VBC, GBC}, grid, args...)

           #  ↑ z ↑  interior
           #  -----  interior face
    kᴵ = 1 #    *    interior cell
    kᴮ = 1 #  =====  bottom boundary
    kᴴ = 0 #    *    halo cell

    @loop_xy i j grid begin
        Δ = ΔzC(i, j, kᴮ, grid) # Δ between first interior and first bottom halo point, defined at cell face.
        @inbounds ∇c = left_gradient(bc, c[i, j, kᴵ], Δ, i, j, grid, args...)
        @inbounds c[i, j, kᴴ] = linearly_extrapolate(c[i, j, kᴵ], ∇c, -Δ) # extrapolate downward in -z direction.
    end

    return nothing
end

function _fill_top_halo!(c, bc::Union{VBC, GBC}, grid, args...)

                     #  ↑ z ↑
    kᴴ = grid.Nz + 1 #    *    halo cell
    kᴮ = grid.Nz + 1 #  =====  top boundary
    kᴵ = grid.Nz     #    *    interior cell
                     #  -----  interior face

    @loop_xy i j grid begin
        Δ = ΔzC(i, j, kᴮ, grid) # Δ between first interior and first top halo point, defined at cell face.
        @inbounds ∇c = right_gradient(bc, c[i, j, kᴵ], Δ, i, j, grid, args...)
        @inbounds c[i, j, kᴴ] = linearly_extrapolate(c[i, j, kᴵ], ∇c, Δ) # extrapolate upward in +z direction.
    end

    return nothing
end


text:
#!/usr/bin/env julia

usage = "
usage: julia [-switches] -- [filename] [args...]

This is a Julia script seed. Customize it to get started quickly with Julia.

Options:
  args     Specify the type of seed.
  switches -h for help, -v for version and more switches can be found on official documentation.

  For more information: http://docs.julialang.org/en/v1/manual/getting-started .
  ";

if(length(ARGS)==0)
	println("You planted a water seed.")
else
	for option in ARGS
			if option == "-h"
				println("For help use -h before the filename");
				println(usage);
			else
				println("You planted a ", option , " seed.")
			end
	end
end

print("Thank you")

ext: jl


text:
using Tmatrix
using Zygote: Zygote
using FiniteDifferences: FiniteDifferences

finite_difference_method = FiniteDifferences.central_fdm(5, 1; factor = 1e-6)

# inputs ====================================================
n_θ_points = 5
n_ϕ_points = 2
m, n, m_, n_ = 2, 3, 2, 5
n_max = 1
k1_r, k1_i, k2_r, k2_i = 1e5, 1e3, 2e5, 3e3
rotationally_symmetric = true
symmetric_about_plane_perpendicular_z = false
rx, rz = 1e-6, 1.2e-6
kind = "regular"
J_superscript = 11
Q_superscript = J_superscript

θ_1D_array = LinRange(1e-16, π, n_θ_points);
ϕ_1D_array = LinRange(1e-16, 2π, n_ϕ_points);
if rotationally_symmetric
    θ_array = collect(θ_1D_array)
    ϕ_array = zeros(size(θ_array))
else
    # eventually, this should be removed. I just keep it for sanity checks.
    θ_array, ϕ_array = Tmatrix.meshgrid(θ_1D_array, ϕ_1D_array)
end
r_array, n̂_array = Tmatrix.ellipsoid(rx, rz, θ_array);

epss = 1e-6

# test "J_mn_m_n__integrand_SeparateRealImag" ====================================================================================
J = Tmatrix.J_mn_m_n__integrand_SeparateRealImag(
    m,
    n,
    m_,
    n_,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array,
    kind,
    J_superscript,
)
J_ = Tmatrix.J_mn_m_n__integrand(
    m,
    n,
    m_,
    n_,
    complex(k1_r, k1_i) .* r_array,
    complex(k2_r, k2_i) .* r_array,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array;
    kind = kind,
    J_superscript = J_superscript,
)

println("==================================================================")
println("testing <<J_mn_m_n__integrand_SeparateRealImag>>")
println("Displaying results from _SeparateRealImag and the complex function")
println("Displaying result from _SeparateRealImag function --------------------")
display(Tmatrix.get_complex_matrix_from_concatenated_real_imag(J))
println("Displaying result from complex function ------------------------------")
display(J_)

println(
    "testing autodiff for <<J_mn_m_n__integrand_SeparateRealImag>> -----------------------------------------",
)
@time J_mn_m_n__integrand_SeparateRealImag_jacobian = Zygote.jacobian(
    Tmatrix.J_mn_m_n__integrand_SeparateRealImag,
    m,
    n,
    m_,
    n_,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array,
    kind,
    J_superscript,
)
function J_mn_m_n__integrand_SeparateRealImag_lean(
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array,
)
    Tmatrix.J_mn_m_n__integrand_SeparateRealImag(
        m,
        n,
        m_,
        n_,
        k1_r,
        k1_i,
        k2_r,
        k2_i,
        r_array,
        θ_array,
        ϕ_array,
        n̂_array,
        kind,
        J_superscript,
    )
end
@time J_mn_m_n__integrand_SeparateRealImag_lean_jacobian = Zygote.jacobian(
    J_mn_m_n__integrand_SeparateRealImag_lean,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array,
)
@time J_jacobian_finite_diff = FiniteDifferences.jacobian(
    finite_difference_method,
    J_mn_m_n__integrand_SeparateRealImag_lean,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array,
)
@time J_jacobian_manual_numerical_diff = (
    (
        J_mn_m_n__integrand_SeparateRealImag_lean(
            k1_r + epss,
            k1_i,
            k2_r,
            k2_i,
            r_array,
            θ_array,
            ϕ_array,
            n̂_array,
        ) - J_mn_m_n__integrand_SeparateRealImag_lean(
            k1_r - epss,
            k1_i,
            k2_r,
            k2_i,
            r_array,
            θ_array,
            ϕ_array,
            n̂_array,
        )
    ) / (2 * epss),
    (
        J_mn_m_n__integrand_SeparateRealImag_lean(
            k1_r,
            k1_i + epss,
            k2_r,
            k2_i,
            r_array,
            θ_array,
            ϕ_array,
            n̂_array,
        ) - J_mn_m_n__integrand_SeparateRealImag_lean(
            k1_r,
            k1_i - epss,
            k2_r,
            k2_i,
            r_array,
            θ_array,
            ϕ_array,
            n̂_array,
        )
    ) / (2 * epss),
    (
        J_mn_m_n__integrand_SeparateRealImag_lean(
            k1_r,
            k1_i,
            k2_r + epss,
            k2_i,
            r_array,
            θ_array,
            ϕ_array,
            n̂_array,
        ) - J_mn_m_n__integrand_SeparateRealImag_lean(
            k1_r,
            k1_i,
            k2_r - epss,
            k2_i,
            r_array,
            θ_array,
            ϕ_array,
            n̂_array,
        )
    ) / (2 * epss),
    (
        J_mn_m_n__integrand_SeparateRealImag_lean(
            k1_r,
            k1_i,
            k2_r,
            k2_i + epss,
            r_array,
            θ_array,
            ϕ_array,
            n̂_array,
        ) - J_mn_m_n__integrand_SeparateRealImag_lean(
            k1_r,
            k1_i,
            k2_r,
            k2_i - epss,
            r_array,
            θ_array,
            ϕ_array,
            n̂_array,
        )
    ) / (2 * epss),
    (
        J_mn_m_n__integrand_SeparateRealImag_lean(
            k1_r,
            k1_i,
            k2_r,
            k2_i,
            r_array .+ epss,
            θ_array,
            ϕ_array,
            n̂_array,
        ) - J_mn_m_n__integrand_SeparateRealImag_lean(
            k1_r,
            k1_i,
            k2_r,
            k2_i,
            r_array .- epss,
            θ_array,
            ϕ_array,
            n̂_array,
        )
    ) ./ (2 * epss),
    (
        J_mn_m_n__integrand_SeparateRealImag_lean(
            k1_r,
            k1_i,
            k2_r,
            k2_i,
            r_array,
            θ_array .+ epss,
            ϕ_array,
            n̂_array,
        ) - J_mn_m_n__integrand_SeparateRealImag_lean(
            k1_r,
            k1_i,
            k2_r,
            k2_i,
            r_array,
            θ_array .- epss,
            ϕ_array,
            n̂_array,
        )
    ) ./ (2 * epss),
)
finite_difference_method = FiniteDifferences.central_fdm(5, 1; factor = 1e-10)
@time J_jacobian_finite_diff = FiniteDifferences.jacobian(
    finite_difference_method,
    J_mn_m_n__integrand_SeparateRealImag_lean,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array,
)
hcat(J_mn_m_n__integrand_SeparateRealImag_lean_jacobian[1], J_jacobian_finite_diff[1])

# test "J_mn_m_n__SeparateRealImag" ================================================================================================
J = Tmatrix.J_mn_m_n__SeparateRealImag(
    m,
    n,
    m_,
    n_,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array,
    kind,
    J_superscript,
    rotationally_symmetric,
)
J_ = Tmatrix.J_mn_m_n_(
    m,
    n,
    m_,
    n_,
    complex(k1_r, k1_i) .* r_array,
    complex(k2_r, k2_i) .* r_array,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array;
    kind = kind,
    J_superscript = J_superscript,
    rotationally_symmetric = rotationally_symmetric,
)

println("==================================================================")
println("Displaying results from _SeparateRealImag and the complex function")
println("Displaying result from _SeparateRealImag function --------------------")
display(Tmatrix.get_complex_matrix_from_concatenated_real_imag(J))
println("Displaying result from complex function ------------------------------")
display(J_)

# the problem that Zygote is trying to differentiate "rotationally_symmetric".
# how can I prevent Zygote from differentiating with respect to a given argument?
Zygote.jacobian(
    (m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array) ->
        Tmatrix.J_mn_m_n__SeparateRealImag(
            m,
            n,
            m_,
            n_,
            k1_r,
            k1_i,
            k2_r,
            k2_i,
            r_array,
            θ_array,
            ϕ_array,
            n̂_array,
            kind,
            J_superscript,
            rotationally_symmetric,
        ),
    m,
    n,
    m_,
    n_,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array,
)

Zygote.jacobian(
    J_mn_m_n__integrand_SeparateRealImag_lean,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array,
)[5]

epss = big(1e-6);
vcat(
    (
        J_mn_m_n__integrand_SeparateRealImag_lean(
            big(k1_r),
            big(k1_i),
            big(k2_r),
            big(k2_i),
            big.(r_array .+ epss),
            big.(θ_array),
            big.(ϕ_array),
            [big.(x) for x in n̂_array],
        ) - J_mn_m_n__integrand_SeparateRealImag_lean(
            big(k1_r),
            big(k1_i),
            big(k2_r),
            big(k2_i),
            big.(r_array .- epss),
            big.(θ_array),
            big.(ϕ_array),
            [big.(x) for x in n̂_array],
        )
    ) ./ (2 * epss)...,
)

single_one_array = zeros(size(r_array));
single_one_array[1] = 1;

for epss in [big(1e-50)]
    display(
        vcat(
            (
                J_mn_m_n__integrand_SeparateRealImag_lean(
                    big(k1_r),
                    big(k1_i),
                    big(k2_r),
                    big(k2_i),
                    big.(r_array + epss * single_one_array),
                    big.(θ_array),
                    big.(ϕ_array),
                    [big.(x) for x in n̂_array],
                ) - J_mn_m_n__integrand_SeparateRealImag_lean(
                    big(k1_r),
                    big(k1_i),
                    big(k2_r),
                    big(k2_i),
                    big.(r_array - epss * single_one_array),
                    big.(θ_array),
                    big.(ϕ_array),
                    [big.(x) for x in n̂_array],
                )
            ) ./ (2 * epss)...,
        ),
    )
    println()

    display(
        vcat(
            (
                J_mn_m_n__integrand_SeparateRealImag_lean(
                    big(k1_r),
                    big(k1_i),
                    big(k2_r),
                    big(k2_i),
                    big.(r_array),
                    big.(θ_array + epss * single_one_array),
                    big.(ϕ_array),
                    [big.(x) for x in n̂_array],
                ) - J_mn_m_n__integrand_SeparateRealImag_lean(
                    big(k1_r),
                    big(k1_i),
                    big(k2_r),
                    big(k2_i),
                    big.(r_array),
                    big.(θ_array - epss * single_one_array),
                    big.(ϕ_array),
                    [big.(x) for x in n̂_array],
                )
            ) ./ (2 * epss)...,
        ),
    )
    println()
end

# test "Q_mn_m_n_SeparateRealImag"
@time Q = Tmatrix.Q_mn_m_n_SeparateRealImag(
    m,
    n,
    m_,
    n_,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array,
    kind,
    Q_superscript,
    rotationally_symmetric,
    symmetric_about_plane_perpendicular_z,
)
@time Q_ = Tmatrix.Q_mn_m_n_(
    m,
    n,
    m_,
    n_,
    complex(k1_r, k1_i),
    complex(k2_r, k2_i),
    complex(k1_r, k1_i) .* r_array,
    complex(k2_r, k2_i) .* r_array,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array;
    kind = kind,
    Q_superscript = Q_superscript,
    rotationally_symmetric = rotationally_symmetric,
)

println("==================================================================")
println("Displaying results from _SeparateRealImag and the complex function")
println("Displaying result from _SeparateRealImag function --------------------")
display(Tmatrix.get_complex_matrix_from_concatenated_real_imag(Q))
println("Displaying result from complex function ------------------------------")
display(Q_)

Zygote.jacobian(
    (m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array) ->
        Tmatrix.Q_mn_m_n_SeparateRealImag(
            m,
            n,
            m_,
            n_,
            k1_r,
            k1_i,
            k2_r,
            k2_i,
            r_array,
            θ_array,
            ϕ_array,
            n̂_array,
            kind,
            Q_superscript,
            rotationally_symmetric,
            symmetric_about_plane_perpendicular_z,
        ),
    m,
    n,
    m_,
    n_,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array,
)

# test "Q_matrix_SeparateRealImag"
@time Q = Tmatrix.Q_matrix_SeparateRealImag(
    n_max,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array,
    kind,
    rotationally_symmetric,
    symmetric_about_plane_perpendicular_z,
)
@time Q_ = Tmatrix.Q_matrix(
    n_max,
    complex(k1_r, k1_i),
    complex(k2_r, k2_i),
    complex(k1_r, k1_i) .* r_array,
    complex(k2_r, k2_i) .* r_array,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array,
    kind = kind,
    rotationally_symmetric = rotationally_symmetric,
    symmetric_about_plane_perpendicular_z = false,
    verbose = false,
)

println("==================================================================")
println("Displaying results from _SeparateRealImag and the complex function")
println("Displaying result from _SeparateRealImag function --------------------")
display(Tmatrix.get_complex_matrix_from_concatenated_real_imag(Q))
println("Displaying result from complex function ------------------------------")
display(Q_)

# The jacobian calculation is so slow. For n_max = 1, it takes 52.409485 seconds. This is slow compared to 0.015419 seconds required to evaluate the function
Zygote.jacobian(
    (n_max, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array) ->
        Tmatrix.Q_matrix_SeparateRealImag(
            n_max,
            k1_r,
            k1_i,
            k2_r,
            k2_i,
            r_array,
            θ_array,
            ϕ_array,
            n̂_array,
            kind,
            rotationally_symmetric,
            symmetric_about_plane_perpendicular_z,
        ),
    n_max,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array,
)

# test "T_matrix_SeparateRealImag"
@time T = Tmatrix.T_matrix_SeparateRealImag(
    n_max,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array,
    rotationally_symmetric,
    symmetric_about_plane_perpendicular_z,
)
@time T_ = Tmatrix.T_matrix(
    n_max,
    complex(k1_r, k1_i),
    complex(k2_r, k2_i),
    complex(k1_r, k1_i) .* r_array,
    complex(k2_r, k2_i) .* r_array,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array;
    rotationally_symmetric = rotationally_symmetric,
    symmetric_about_plane_perpendicular_z = false,
)

println("==================================================================")
println("Displaying results from _SeparateRealImag and the complex function")
println("Displaying result from _SeparateRealImag function --------------------")
display(Tmatrix.get_complex_matrix_from_concatenated_real_imag(T))
println("Displaying result from complex function ------------------------------")
display(T_)

@time T = Tmatrix.T_matrix_SeparateRealImag(
    n_max,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array,
    rotationally_symmetric,
    symmetric_about_plane_perpendicular_z,
)
@time ∂T = Zygote.jacobian(
    Tmatrix.T_matrix_SeparateRealImag,
    n_max,
    k1_r,
    k1_i,
    k2_r,
    k2_i,
    r_array,
    θ_array,
    ϕ_array,
    n̂_array,
    rotationally_symmetric,
    symmetric_about_plane_perpendicular_z,
)


text:
using NumericalMethods
using Test

@testset "NumericalMethods.jl" begin
    # Write your tests here.
end


text:
#JuliaRL/src/JuliaRL.jl

include("./BasicDQN.jl")
include("./DQN.jl")
include("./PDQN.jl")
include("./A2C.jl")

include("./display.jl")
include("./show_stats.jl")

max_stars_repo_name: FeLusiani/JuliaRL

max_issues_repo_name: FeLusiani/JuliaRL

max_forks_repo_name: FeLusiani/JuliaRL


text:
using DashBootstrapComponents, DashHtmlComponents

placement_selector = html_div(
    [
        dbc_label("Placement:"),
        dbc_radioitems(
            id = "offcanvas-placement-selector",
            options = [
                Dict("label" => "start (Default)", "value" => "start"),
                Dict("label" => "end", "value" => "end"),
                Dict("label" => "top", "value" => "top"),
                Dict("label" => "bottom", "value" => "bottom"),
            ],
            value = "start",
            inline = true,
        ),
    ],
    className = "mb-2",
)

offcanvas = html_div([
    placement_selector,
    dbc_button("Open Offcanvas", id = "open-offcanvas-placement", n_clicks = 0),
    dbc_offcanvas(
        html_p("Some offcanvas content..."),
        id = "offcanvas-placement",
        title = "Placement",
        is_open = false,
    ),
])

callback!(
    app,
    Output("offcanvas-placement", "is_open"),
    Input("open-offcanvas-placement", "n_clicks"),
    State("offcanvas-placement", "is_open"),
) do n1, is_open
    return n1 > 0 ? is_open == 0 : is_open
end;

callback!(
    app,
    Output("offcanvas-placement", "placement"),
    Input("offcanvas-placement-selector", "value"),
) do placement
    return placement
end;

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
"""
return an table of the amount of each element in every molecule
as well as what the number of the rows and columns stand for.
"""

function chemical_parse(str::String)
    elems = Dict{String, Integer}()

    eldt = [
        begin
            ## remove parentheses and brackets
            flag = true
            while flag
                flag = false
                m = match(r".*?\((.*?)\)([0-9]*)", elemdata)
                m != nothing && ((flag=true; k = m[2]=="" ? 1 : parse(Int, m[2]));
                    elemdata = replace(elemdata, r"(\(.*?\)[0-9]*)",  m[1]^k, 1))
                m = match(r".*?\[(.*?)\]([0-9]*)", elemdata)
                m != nothing && ((flag=true; k = m[2]=="" ? 1 : parse(Int, m[2]));
                    elemdata = replace(elemdata, r"(\[.*?\][0-9]*)",  m[1]^k, 1))
                m = match(r".*?\{(.*?)\}([0-9]*)", elemdata)
                m != nothing && ((flag=true; k = m[2]=="" ? 1 : parse(Int, m[2]));
                    elemdata = replace(elemdata, r"(\{.*?\}[0-9]*)",  m[1]^k, 1))
            end
            ## arrange elements' names

            len = length(elemdata)
            for i in 1:len
                '0' ≤ elemdata[i] ≤ '9' && continue
                if 'A' ≤ elemdata[i] ≤ 'Z'
                    name = i < len && 'a' ≤ elemdata[i+1] ≤ 'z' ?
                        elemdata[i:i+1] : elemdata[i:i]
                    haskey(elems, name) || (elems[name] = length(elems)+1)
                end
            end
            elemdata
        end
        for elemdata in split(str)]

    arr = zeros(Int, length(elems), length(eldt))
    ## arr[<element>, <molecule>]
    for (k, elemdata) in enumerate(eldt)
        len = length(elemdata)
        i = 1
        while i ≤ len
            @assert  'A' ≤ elemdata[i] ≤ 'Z'
            name =
                if i < len && 'a' ≤ elemdata[i+1] ≤ 'z'
                    i += 1
                    elemdata[i-1:i]
                else
                    elemdata[i:i]
                end
            i += 1
            amount_str = ""
            while i ≤ len && '0' ≤ elemdata[i] ≤ '9'
                amount_str *= elemdata[i:i]
                i += 1
            end
            arr[elems[name], k] = amount_str=="" ? 1 : parse(Int, amount_str)
        end
    end
    row = split(str)
    revelems = Dict(v=>k for (k, v) in elems)
    column = [revelems[i] for i in 1:length(revelems)]
    return arr, row, column
end

max_stars_repo_name: sclereid/collection

max_issues_repo_name: sclereid/collection

max_forks_repo_name: sclereid/collection


text:
using NODAL, Distributed, Statistics, Test

include("parameter.jl")
include("configuration.jl")
include("run.jl")
include("simulated_annealing.jl")
include("iterative_first_improvement.jl")
include("iterative_probabilistic_improvement.jl")
include("randomized_first_improvement.jl")
include("iterative_greedy_construction.jl")
include("iterated_local_search.jl")
include("unit_value.jl")


text:
using ApplicationBuilder

examples_hello = joinpath(@__DIR__, "..", "..", "examples", "hello.jl")

# Allow this file to be called either as a standalone file to build the above
# example, or from runtests.jl using a provided builddir.
@isdefined(builddir) || (builddir="builddir")

ApplicationBuilder.build_app_bundle(examples_hello;
                          verbose=true, appname="HelloWorld", builddir=builddir)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMTruss.jl/blob/master/LICENSE

""" Truss implementation for JuliaFEM. """
module FEMTruss

using FEMBase

import FEMBase: get_unknown_field_name,
                get_formulation_type,
                assemble!

"""
This truss fromulation is from
Concepts and applications of finite element analysis, forth edition
Chapter 2.4
Cook, Malkus, Plesha, Witt

# Features
- Nodal forces can be set using element `Poi1` with field
  `nodal force i`, where `i` is dof number.
- Displacements can be fixed using element `Poi1` with field
  `fixed displacement i`, where `i` is dof number.
- Temperature loads can be set by setting field `thermal expansion coefficient`
  and field `temperature difference`
"""
type Truss <: FieldProblem
end

function get_unknown_field_name(::Problem{Truss})
    return "displacement"
end

function get_formulation_type(::Problem{Truss})
    return :total
end

"""
get_truss_1d_and_Tg(element::Element{Seg2})
This function sets up the local truss 1d element needed for
the integration points to be used for the stiffness and forces
It also sets up the coordinate transformation matrix
Will need to change allocation strategy to pre-allocation later
"""
function get_truss_1d_and_Tg(element::Element{Seg2}, ndim)
    node_id1 = element.connectivity[1] #First node in element
    node_id2 = element.connectivity[2] #second node in elements
    pos = element("geometry")
    dl = pos[node_id2]-pos[node_id1];
    l = sqrt(dot(dl,dl));
    # Do the local element for calculations
    elem_1d = Element(Seg2, [1,2])
    elem_1d.id=-1 # To signal it is local
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0],
        2 => [l])
    update!(elem_1d, "geometry", X)
    if ndim == 1
        T = eye(Float64, 2)
    elseif ndim == 2
        l_theta = dl[1]/l
        m_theta = dl[2]/l
        T = [l_theta m_theta 0 0; 0 0 l_theta m_theta]
    else # All three dimensions
        l_theta = dl[1]/l
        m_theta = dl[2]/l
        n_theta = dl[3]/l
        T = [l_theta m_theta n_theta 0 0 0; 0 0 0 l_theta m_theta n_theta]
    end
    return (elem_1d,T)
end

"""
   get_truss_1d_K(element::Element{Seg2}, nnodes, time)
   This function assembles the 1d truss stiffness matrix
   Will need to change allocation strategy to pre-allocation later
   Can discuss if we really need nnodes, since that is always 2 for trusses
"""
function get_truss_1d_K(elem_1d::Element{Seg2}, nnodes, time)
    K_loc = zeros(nnodes, nnodes)
    for ip in get_integration_points(elem_1d)
        dN = elem_1d(ip, time, Val{:Grad})
        detJ = elem_1d(ip, time, Val{:detJ})
        A = elem_1d("cross section area", ip, time)
        E = elem_1d("youngs modulus", ip, time)
        K_loc += ip.weight*E*A*dN'*dN*detJ
    end
    return K_loc
end

"""
    has_temperature_load(elem_1d::Element{Seg2})
    checks if we have a temperature load on the truss element
"""
function has_temperature_load(elem_1d::Element{Seg2})
    return  haskey(elem_1d, "temperature difference")&& haskey(elem_1d, "thermal expansion coefficient")
end

"""
   get_truss_1d_f(element::Element{Seg2}, nnodes, time)
   This function assembles the 1d truss force vector
   Will need to change allocation strategy to pre-allocation later
   Can discuss if we really need nnodes, since that is always 2 for trusses
"""
function get_truss_1d_f(elem_1d::Element{Seg2}, nnodes, time)
    fe = zeros(nnodes)
    if has_temperature_load(elem_1d)
        for ip in get_integration_points(elem_1d)
            dt = elem_1d("temperature difference", ip, time)
            lambda = elem_1d("thermal expansion coefficient", ip, time)
            A = elem_1d("cross section area", ip, time)
            E = elem_1d("youngs modulus", ip, time)
            N = elem_1d(ip, time)
            detJ = elem_1d(ip, time, Val{:detJ})
            f = E*A*lambda*dt
            fe += ip.weight * N'*f # *detJ
        end
        fe[end]*=-1.0
    end
    # how do we generally set the right sides negative
    return fe
end

"""
   get_Kg_and_Tg(element::Element{Seg2}, nnodes, ndim, time)
   This function assembles the local stiffness uses global transformation matrix
   to make the global version of the local stiffnes matrix
   Will need to change allocation strategy to pre-allocation later
   Can discuss if we really need nnodes, since that is always 2 for trusses
"""

function get_Kg_and_Tg(element::Element{Seg2}, nnodes, ndim, time)
    elem_1d, T = get_truss_1d_and_Tg(element, ndim)
    # Update props for stiffness attributes
    update!(elem_1d, "youngs modulus", element("youngs modulus"))
    update!(elem_1d, "cross section area", element("cross section area"))
    K_loc = get_truss_1d_K(elem_1d, nnodes, time)
    ndofs = nnodes*ndim
    K = zeros(ndofs,ndofs)
    K = T'*K_loc*T

    return (K,T)
end

function get_F_local(element::Element{Seg2}, nnodes, ndim, time)
    elem_1d, T = get_truss_1d_and_Tg(element, ndim)
    # Update props for stiffness attributes
    update!(elem_1d, "youngs modulus", element("youngs modulus"))
    update!(elem_1d, "cross section area", element("cross section area"))
    if (has_temperature_load(element))
        update!(elem_1d, "thermal expansion coefficient", element("thermal expansion coefficient"))
        update!(elem_1d, "temperature difference", element("temperature difference"))
    end
    f = get_truss_1d_f(elem_1d, nnodes, time)
    return f
end

"""
    assemble!(assembly:Assembly, problem::Problem{Elasticity}, elements, time)
Start finite element assembly procedure for Elasticity problem.
Function groups elements to arrays by their type and assembles one element type
at time. This makes it possible to pre-allocate matrices common to same type
of elements.
"""
function assemble!(assembly::Assembly, problem::Problem{Truss},
                   element::Element{Seg2}, time)
    #Require that the number of nodes = 2 ?
    nnodes = length(element)
    ndim = get_unknown_field_dimension(problem)
    elem_1d, T = get_truss_1d_and_Tg(element, ndim)
    # setuo the things needed for local stiffness matrix
    update!(elem_1d, "youngs modulus", element("youngs modulus"))
    update!(elem_1d, "cross section area", element("cross section area"))
    K_loc = get_truss_1d_K(elem_1d, nnodes, time)
    ndofs = nnodes*ndim
    K = zeros(ndofs,ndofs)
    K = T'*K_loc*T
    gdofs = get_gdofs(problem, element)
    add!(assembly.K, gdofs, gdofs, K)
    if has_temperature_load(element)
        update!(elem_1d, "thermal expansion coefficient", element("thermal expansion coefficient"))
        update!(elem_1d, "temperature difference", element("temperature difference"))
        f = get_truss_1d_f(elem_1d, nnodes, time)*-1.0 #Needs to be negated
        fg = T'*f
        add!(assembly.f, gdofs, fg)
    end
    #add!(assembly.f, gdofs, f)
end

import FEMBase: assemble_elements!

function assemble_elements!(problem::Problem, assembly::Assembly,
                            elements::Vector{Element{Poi1}}, time::Float64)
    dim = ndofs = get_unknown_field_dimension(problem)
    Ce = zeros(ndofs,ndofs)
    ge = zeros(ndofs)
    fe = zeros(ndofs)
    for element in elements
        fill!(Ce, 0.0)
        fill!(ge, 0.0)
        fill!(fe, 0.0)
        ip = (0.0,)
        for i=1:dim
            if haskey(element, "fixed displacement $i")
                Ce[i,i] = 1.0
                ge[i] = element("fixed displacement $i", ip, time)
            end
            if haskey(element, "nodal force $i")
                fe[i] = element("nodal force $i", ip, time)
            end
        end
        gdofs = get_gdofs(problem, element)
        add!(assembly.C1, gdofs, gdofs, Ce)
        add!(assembly.C2, gdofs, gdofs, Ce)
        add!(assembly.g, gdofs, ge)
        add!(assembly.f, gdofs, fe)
    end
end

export Truss, get_Kg_and_Tg, get_F_local

end


text:
using HDF5
using Test

hf = h5open(tempname(), "w")

fv = 3.14
data = [1.,2.,3.,4.,5.,6.]
floatsize = sizeof(data[1])
h5t = datatype(data[1])
title = "lal"
name = "mym"
nfield = 2
nrec = 3
recsize = nfield * floatsize
colname = ["f1_verylongnameforfun", "f2"]
offset = [0,floatsize]
tid = [h5t.id, h5t.id]
chunk = 7
fillvalue = [3.14, 2.71]
compress = 1

HDF5.h5tb_make_table(title, hf, name, nfield, nrec, recsize, colname, offset, tid, chunk, fillvalue, compress, data)
fieldsize = [floatsize, floatsize]
HDF5.h5tb_append_records(hf, name, nrec, recsize, offset, fieldsize, data)
HDF5.h5tb_write_records(hf, name, 1, 4, recsize, offset, fieldsize, collect(1:8) .+ 20.0)
buf = fill(0.0, 100)

HDF5.h5tb_read_table(hf, name, recsize, offset, fieldsize, buf)
@test buf[1:12] == [1.0, 2.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 5.0, 6.0]
buf .= 0.0
HDF5.h5tb_read_records(hf, name, 2, 3, recsize, offset, fieldsize, buf)
@test buf[1:6] == collect(23:28)

h5_nfields, h5_nrec = HDF5.h5tb_get_table_info(hf, name)
@test h5_nfields == nfield
@test h5_nrec == 6

h5_colname, h5_fieldsize, h5_offset, h5_recsize = HDF5.h5tb_get_field_info(hf, name)

@test h5_colname == colname
@test h5_fieldsize == fieldsize
@test h5_offset == offset
@test h5_recsize == recsize

close(hf)


text:
# functions for doing complex step

import Base: flipsign

"""
  Extends Base.flipsign with method for complex numbers
"""
function flipsign(a::T, b::Number) where {T <: Complex}

  xreal = flipsign(real(a), b)
  ximag = flipsign(imag(a), b)

  return T(xreal, ximag)
end

@doc """
### EulerEquationMod.absvalue

  This function computes the absolute value in a way such that the complex
  step method works.

  There is probably a bit-twiddling way of doing this faster

"""
function absvalue(x::T) where T <: Complex

#  c_part = flipsign(imag(x), real(x))
#  return T(abs(real(x)), c_part)
  return flipsign(x, real(x))
end
#=
  x2 = complex(real(x), imag(x))

  if (real(x) < 0)
    x2 *= -1
  end

#  println("after modification, x2 = ", x2)
  return x2
end
=#
function absvalue(x)   # general fallback method
  return abs(x)
end

# array method
function absvalue(x::AbstractArray{T}) where T

  x2 = copy(x)
  for i=1:length(x)
    x2[i] = absvalue(x[i])
  end

  return x2
end

function absvalue_rev(x::Number, x_bar::Number)

#  if real(x) < 0
#    x = -x
#    x_bar = -x_bar
#  end

  # this function has the odd property that forward and reverse modes are
  # equivalent
  x2 = flipsign(x, real(x))
  x2_bar = flipsign(x_bar, real(x))

  return x2, x2_bar
end

@doc """
###Utils.absvalue_deriv

Computes the derivative of the absolute value of a variable w.r.t itself.

**Inputs**

* `val` : The variable whose derivative needs to be computed e.r.t itself

"""->

function absvalue_deriv(val::Tval) where Tval

  Tval(sign_c(val))
  #=
  if val > zero(Tval)
    return one(Tval)
  elseif val < zero(Tval)
    return -one(Tval)
  else
    return zero(Tval)
  end
  =#

end # End function absvalue_deriv

#------------------------------------------------------------------------------
# replace abs with a cubic spline fit near 0, to make it differentiable

"""
  A replacement for absvalue() that is differentiable everywhere (C^1 continuous)

  Unfortunately, it is somewhat more expensive than absvalue()
"""
function absvalue2(val::Number)

  delta = 1e-13
  if absvalue(val) > delta
    return absvalue(val)
  elseif real(val) < 0
    tmp = val + 100e-15
    return @evalpoly tmp 100e-15 -1.0 -10e12 100e24
  else  # real(val) > 0
    tmp = val
    return @evalpoly tmp 0.0 0.0 20e12 -100e24
  end
end

"""
  Derivative of abs()
"""
function absvalue2_deriv(val::T) where {T <: Number}

  delta = 1e-13
  if absvalue(val) > delta
    return sign_c(val)
  elseif real(val) < 0
    tmp = val + 100e-15
     return @evalpoly tmp -1.0 -20e12 300e24
  else  # real(val) > 0
    tmp = val
    return @evalpoly tmp 0.0 40e12 -300e24
  end
end

function max_deriv_rev(x::Tval, y::Tval, max_val_bar::Tval) where Tval

  x_bar = zero(Tval)
  y_bar = zero(Tval)

  if real(x) > real(y)
    x_bar += max_val_bar
  else
    y_bar += max_val_bar
  end

  return x_bar, y_bar
end # End function max_rev

global const ABSVALUE3_DELTA = 1e-7
"""
  Smooth abs() function based on Harten's second entropy fix
"""
function absvalue3(val::Number)

  delta = ABSVALUE3_DELTA
  val1 = absvalue(val)
  if val1 > delta
    return val1
  else
    term1 = val*val/delta
    term2 = delta
    return ((val*val)/delta + delta)/2
  end
end

function absvalue3_deriv(val::T) where {T <: Number}

  delta = ABSVALUE3_DELTA

  val1 = absvalue(val)
  if val1 > delta
    return absvalue_deriv(val)
  else
    return val/delta
  end
end

"""
  Complex-step safe version of the sign function
"""
function sign_c(a::Number)

  return sign(real(a))
end

import Base.isless

@doc """
### EulerEquationMod.isless

  This function defines the isless() function (ie. <  operator) in such a way
  that algorithmic differentiation works

"""->
function isless(x::Complex, y::Complex)
  return real(x) < real(y)
end

function isless(x::Real, y::Complex)
  return x < real(y)
end

function isless(x::Complex, y::Real)
  return real(x) < y
end
# there is no isgreater, just a reuse of isless

#------------------------------------------------------------------------------
# trig functions
# most are defind in Base, but a few are needed

import Base: atan2

"""
  atan2 for complex number.  This function protects against numerical errors
  if x -> 0, but I don't know if it is as accurate as Base.atan2
"""
function atan2(y::Complex, x::Complex)

  x2 = absvalue(x)
  y2 = absvalue(y)
  # avoid numerical problems if x -> 0
  if x2 > y2
    theta = atan(y2/x2)
  else
    theta = pi/2 - atan(x2/y2)
  end

  if x < 0
    theta = pi - theta
  end

  if y < 0
    theta = -theta
  end

  return theta
end

# if only one complex number promote both to complex
function atan2(y::Number, x::T) where {T<:Complex}
  return atan2(T(y), x)
end

function atan2(y::T, x::Number) where {T<:Complex}
  return atan2(y, T(x))
end

"""
  Reverse-mode differentiated version of atan2.

  **Inputs**

   * y
   * x
   * theta_bar

  **Outputs**

   * theta
   * y_bar
   * x_bar

"""
function atan2_rev(y::Number, x::Number, theta_bar::Number)

  x2 = absvalue(x)
  y2 = absvalue(y)
  # avoid numerical problems if x -> 0
  if x2 > y2
    theta = atan(y2/x2)
  else
    theta = pi/2 - atan(x2/y2)
  end

  if x < 0
    theta = pi - theta
  end

  if y < 0
    theta = -theta
  end

  # reverse sweep
  x_bar = zero(x)
  y_bar = zero(y)
  x2_bar = zero(x)
  y2_bar = zero(y)

  if y < 0
    theta_bar = -theta_bar
  end

  if x < 0
    theta_bar = -theta_bar
  end

  if x2 > y2
    t1 = theta_bar/(1 + y2*y2/(x2*x2))
    x2_bar += t1*-y2/(x2*x2)
    y2_bar += t1/x2
  else
    t2 = -theta_bar/(1 + x2*x2/(y2*y2))
    x2_bar += t2/y2
    y2_bar += t2*-x2/(y2*y2)
  end

  # there is no += for functions that return multiple values
  t3, t3_bar = absvalue_rev(x, x2_bar)
  x_bar += t3_bar

  t4, t4_bar = absvalue_rev(y, y2_bar)
  y_bar += t4_bar

  return theta, y_bar, x_bar
end


text:
module VisualiseTests

using Test
import MetaGraphs

using PicoQuant

@testset "Visualisation tests" begin
    @testset "Test contracting a network" begin

        qasm_str = """OPENQASM 2.0;
                    include "qelib1.inc";
                    qreg q[3];
                    h q[0];
                    cx q[0],q[1];
                    cx q[1],q[2];"""

        circ = load_qasm_as_circuit(qasm_str)
        tng = convert_qiskit_circ_to_network(circ, InteractiveBackend())
        add_input!(tng, "000")
        add_output!(tng, "000")
        g = PicoQuant.create_graph(tng)

        # Is there only one node left after the contraction?
        @test begin
            MetaGraphs.nv(g) == 9
        end

        @test begin
            MetaGraphs.ne(g) == 8
        end
    end
end
end


text:
#=
https://docs.julialang.org/en/stable/manual/interfaces/#man-interface-iteration-1

for i = I   # or  "for i in I"
    # body
end

is translated into:

state = start(I)
while !done(I, state)
    (i, state) = next(I, state)
    # body
end
=#

@inline worksheet(r::SheetRow) = r.sheet
@inline worksheet(itr::SheetRowIterator) = itr.sheet

# creates SheetRow with unpopulated rowcells
function SheetRow(ws::Worksheet, row::Int, xml_element::EzXML.Node)
    @assert EzXML.nodename(xml_element) == "row"
    return SheetRow(ws, row, xml_element, Dict{Int, Cell}(), false)
end

function populate_row_cells!(r::SheetRow)
    if !r.is_rowcells_populated
        for c in EzXML.eachelement(r.row_xml_element)
            if EzXML.nodename(c) == "c"
                cell = Cell(c)
                @assert row_number(cell) == r.row "Malformed Excel file. range_row = $(r.row), cell.ref = $(cell.ref)"
                r.rowcells[column_number(cell)] = cell
            end
        end
        r.is_rowcells_populated = true
    end
    nothing
end

Base.start(itr::SheetRowIterator) = start(itr.xml_rows_iterator)
Base.done(itr::SheetRowIterator, state) = done(itr.xml_rows_iterator, state)

#(i, state) = next(I, state)
function Base.next(itr::SheetRowIterator, state)
    xml_element, next_state = next(itr.xml_rows_iterator, state)
    row = parse(Int, xml_element["r"])
    return SheetRow(worksheet(itr), row, xml_element), next_state
end

function find_row(itr::SheetRowIterator, row::Int) :: SheetRow
    for r in itr
        if row_number(r) == row
            return r
        end
    end
    error("Row $row not found.")
end

function SheetRowIterator(ws::Worksheet)
    xroot = EzXML.root(ws.data)
    @assert EzXML.nodename(xroot) == "worksheet" "Malformed sheet $(ws.name)."

    # finds sheetData element
    for sheet_data_element in EzXML.eachelement(xroot)
        if EzXML.nodename(sheet_data_element) == "sheetData"
            return SheetRowIterator(ws, EzXML.eachelement(sheet_data_element))
        end
    end

    error("Malformed sheet $(ws.name): couldn't find sheetData element.")
end

row_number(sr::SheetRow) = sr.row

function getcell(r::SheetRow, column_index::Int) :: AbstractCell
    populate_row_cells!(r)

    if haskey(r.rowcells, column_index)
        return r.rowcells[column_index]
    else
        return EmptyCell(CellRef(row_number(r), column_index))
    end
end

function getcell(r::SheetRow, column_name::AbstractString)
    @assert is_valid_column_name(column_name) "$column_name is not a valid column name."
    return getcell(r, decode_column_number(column_name))
end

getdata(r::SheetRow, column) = celldata(worksheet(r), getcell(r, column))

Base.getindex(r::SheetRow, x) = getdata(r, x)

"""
    eachrow(sheet)

Creates a row iterator for a worksheet.

Example: Query all cells from columns 1 to 4.

```julia
left = 1  # 1st column
right = 4 # 4th column
for sheetrow in XLSX.eachrow(sheet)
    for column in left:right
        cell = XLSX.getcell(sheetrow, column)

        # do something with cell
    end
end
```
"""
eachrow(ws::Worksheet) = SheetRowIterator(ws)

#
# Table
#

function Base.isempty(sr::SheetRow)
    populate_row_cells!(sr)
    return isempty(sr.rowcells)
end

"""
   column_bounds(sr::SheetRow)

Returns a tuple with the first and last index of the columns for a `SheetRow`.
"""
function column_bounds(sr::SheetRow)

    @assert !isempty(sr) "Can't get column bounds from an empty row."

    local first_column_index::Int = first(keys(sr.rowcells))
    local last_column_index::Int = first_column_index

    for k in keys(sr.rowcells)
        if k < first_column_index
            first_column_index = k
        end

        if k > last_column_index
            last_column_index = k
        end
    end

    return (first_column_index, last_column_index)
end

# anchor_column will be the leftmost column of the column_bounds
function last_column_index(sr::SheetRow, anchor_column::Int) :: Int

    @assert !isempty(getcell(sr, anchor_column)) "Can't get column bounds based on an empty anchor cell."

    local first_column_index::Int = anchor_column
    local last_column_index::Int = first_column_index

    if length(keys(sr.rowcells)) == 1
        return anchor_column
    end

    columns = sort(collect(keys(sr.rowcells)))
    first_i = findfirst(colindex -> colindex == anchor_column, columns)
    last_column_index = anchor_column

    for i in (first_i+1):length(columns)
        if columns[i] - 1 != last_column_index
            return last_column_index
        end

        last_column_index = columns[i]
    end

    return last_column_index
end

"""
    TableRowIterator(sheet, [columns]; [first_row], [column_labels], [header])

`header` is a boolean indicating wether the first row of the table is a table header.

If `header == false` and no `names` were supplied, column names will be generated following the column names found in the Excel file.
Also, the column range will be inferred by the non-empty contiguous cells in the first row of the table.

The user can replace column names by assigning the optional `names` input variable with a `Vector{Symbol}`.

See also `gettable`.
"""
function TableRowIterator(sheet::Worksheet, cols::Union{ColumnRange, AbstractString}; first_row::Int=_find_first_row_with_data(sheet, convert(ColumnRange, cols).start), column_labels::Vector{Symbol}=Vector{Symbol}(), header::Bool=true)
    itr = SheetRowIterator(sheet)
    column_range = convert(ColumnRange, cols)

    if isempty(column_labels)
        if header
            # will use celldata to get column names
            for column_index in column_range.start:column_range.stop
                sheet_row = find_row(itr, first_row)
                cell = getcell(sheet_row, column_index)
                @assert !isempty(cell) "Header cell can't be empty."
                push!(column_labels, Symbol(celldata(sheet, cell)))
            end
        else
            # generate column_labels if there's no header information anywhere
            for c in column_range
                push!(column_labels, c)
            end
        end
    else
        # check consistency for column_range and column_labels
        @assert length(column_labels) == length(column_range) "`column_range` (length=$(length(column_range))) and `column_labels` (length=$(length(column_labels))) must have the same length."
    end

    first_data_row = header ? first_row + 1 : first_row
    return TableRowIterator(itr, Index(column_range, column_labels), first_data_row)
end

function TableRowIterator(sheet::Worksheet; first_row::Int = 1, column_labels::Vector{Symbol}=Vector{Symbol}(), header::Bool=true)
    for r in eachrow(sheet)

        # skip rows until we reach first_row
        if row_number(r) < first_row
            continue
        end

        if !isempty(r)
            columns_ordered = sort(collect(keys(r.rowcells)))
            for (ci, cn) in enumerate(columns_ordered)
                if !Missings.ismissing(getdata(r, cn))
                    # found a row with data. Will get ColumnRange from non-empty consecutive cells
                    first_row = row_number(r)
                    column_start = cn
                    column_stop = cn

                    if length(columns_ordered) == 1
                        # there's only one column
                        column_range = ColumnRange(column_start, column_stop)
                        return TableRowIterator(sheet, column_range; first_row=first_row, column_labels=column_labels, header=header)
                    else
                        # will figure out the column range
                        for ci_stop in (ci+1):length(columns_ordered)
                            cn_stop = columns_ordered[ci_stop]
                            if cn_stop - 1 != column_stop
                                column_range = ColumnRange(column_start, column_stop)
                                return TableRowIterator(sheet, column_range; first_row=first_row, column_labels=column_labels, header=header)
                            end
                            column_stop = cn_stop
                        end
                    end

                    # if got here, it's because all columns are non-empty
                    column_range = ColumnRange(column_start, column_stop)
                    return TableRowIterator(sheet, column_range; first_row=first_row, column_labels=column_labels, header=header)
                end
            end
        end
    end

    error("Couldn't find a table in sheet $(sheet.name)")
end

function _find_first_row_with_data(sheet::Worksheet, column_number::Int)
    # will find first_row
    for r in eachrow(sheet)
        if !Missings.ismissing(getdata(r, column_number))
            return row_number(r)
        end
    end
    error("Column $(encode_column_number(column_number)) has no data.")
end

@inline worksheet(tri::TableRowIterator) = tri.itr.sheet
@inline worksheet(r::TableRow) = worksheet(r.itr)

"""
Returns real sheet column numbers (based on cellref)
"""
@inline sheet_column_numbers(i::Index) = values(i.column_map)
@inline table_column_numbers(i::Index) = eachindex(i.column_labels)
@inline table_column_numbers(r::TableRow) = table_column_numbers(r.itr.index)

"""
Maps table column index (1-based) -> sheet column index (cellref based)
"""
@inline table_column_to_sheet_column_number(index::Index, table_column_number::Int) = index.column_map[table_column_number]
@inline table_columns_count(i::Index) = length(i.column_labels)
@inline table_columns_count(itr::TableRowIterator) = table_columns_count(itr.index)
@inline table_columns_count(r::TableRow) = table_columns_count(r.itr)
@inline table_row_number(r::TableRow) = r.table_row_index
@inline sheet_row_number(r::TableRow) = row_number(r.sheet_row)
@inline get_column_labels(index::Index) = index.column_labels
@inline get_column_labels(itr::TableRowIterator) = get_column_labels(itr.index)
@inline get_column_labels(r::TableRow) = get_column_labels(r.itr)
@inline get_column_label(r::TableRow, table_column_number::Int) = get_column_labels(r)[table_column_number]

# iterate into TableRow to get each column value
Base.start(r::TableRow) = start(table_column_numbers(r))
Base.done(r::TableRow, state) = done(table_column_numbers(r), state)

function Base.next(r::TableRow, state)
    (next_column_number, next_state) = next(table_column_numbers(r), state)
    return (r[next_column_number], next_state)
end

function getcell(r::TableRow, table_column_number::Int)
    table_row_iterator = r.itr
    index = table_row_iterator.index
    sheet_row_iterator = table_row_iterator.itr
    sheet = worksheet(r)
    sheet_column = table_column_to_sheet_column_number(index, table_column_number)
    sheet_row = r.sheet_row
    return getcell(sheet_row, sheet_column)
end

getdata(r::TableRow, table_column_number::Int) = celldata(worksheet(r), getcell(r, table_column_number))

function getdata(r::TableRow, column_label::Symbol)
    index = r.itr.index
    if haskey(index.lookup, column_label)
        return getindex(r, index.lookup[column_label])
    else
        error("Invalid column label: $column_label.")
    end
end

Base.getindex(r::TableRow, x) = getdata(r, x)

function Base.start(itr::TableRowIterator)
    last_state = start(itr.itr)

    # go to the first_data_row
    while !done(itr.itr, last_state)
        (sheet_row, state) = next(itr.itr, last_state)
        if row_number(sheet_row) == itr.first_data_row
            return TableRowIteratorState(state, sheet_row, 1, itr.first_data_row, false)
        end
        last_state = state
    end

    return TableRowIteratorState(state, sheet_row, 1, itr.first_data_row, true)
end

function Base.next(itr::TableRowIterator, state::TableRowIteratorState)
    sheet_row_itr_is_done = done(itr.itr, state.state)
    current_sheet_row = state.sheet_row
    current_sheet_row_number = row_number(current_sheet_row)

    if !sheet_row_itr_is_done
        next_sheet_row, next_state = next(itr.itr, state.state)
    else
        next_sheet_row, next_state = current_sheet_row, state
    end

    return TableRow(itr, current_sheet_row, state.table_row_index), TableRowIteratorState(next_state, next_sheet_row, state.table_row_index+1, current_sheet_row_number, sheet_row_itr_is_done)
end

function Base.done(itr::TableRowIterator, state::TableRowIteratorState)
    # empty table case
    if state.is_done
        return true
    end

    # check skipping rows
    if row_number(state.sheet_row) != itr.first_data_row && row_number(state.sheet_row) != (state.last_sheet_row_number + 1)
        return true
    end

    # check empty rows
    if isempty(state.sheet_row)
        return true
    end

    # check if there are any data inside column range
    for c in sheet_column_numbers(itr.index)
        if !Missings.ismissing(celldata(worksheet(itr), getcell(state.sheet_row, c)))
            return false
        end
    end

    return true
end

function infer_eltype(v::Vector{Any})
    local hasmissing::Bool = false
    local t::DataType = Any

    if isempty(v)
        return eltype(v)
    end

    for i in 1:length(v)
        if Missings.ismissing(v[i])
            hasmissing = true
        else
            if t != Any && typeof(v[i]) != t
                return Any
            else
                t = typeof(v[i])
            end
        end
    end

    if t == Any
        return Any
    else
        if hasmissing
            return Union{Missings.Missing, t}
        else
            return t
        end
    end
end

infer_eltype(v::Vector{T}) where T = T

function gettable(itr::TableRowIterator; infer_eltypes::Bool=false)
    column_labels = get_column_labels(itr)
    columns_count = table_columns_count(itr)
    data = Vector{Any}(columns_count)
    for c in 1:columns_count
        data[c] = Vector{Any}()
    end

    for r in itr # r is a TableRow
        for (ci, cv) in enumerate(r) # iterate a TableRow to get column data
            push!(data[ci], cv)
        end
    end

    if infer_eltypes
        rows = length(data[1])
        for c in 1:columns_count
            new_column_data = Vector{infer_eltype(data[c])}(rows)
            for r in 1:rows
                new_column_data[r] = data[c][r]
            end
            data[c] = new_column_data
        end
    end

    return data, column_labels
end


text:
"""
    Agent

The agent type for the party negotiation model.
"""
mutable struct Agent
    id::Integer
    party::String
    opinions::AbstractArray
end

"""
    Model

The model type for the party negotiation model.
"""
mutable struct Model
    parameter_set::ParameterSet
    agents::Array{Agent}
end

"""
    setup_model(params::ParameterSet, db::SQLite.DB)

Create a model from a given set of specifications and a database.
"""
function setup_model(params::ParameterSet, db::SQLite.DB)
    agents = create_agents(params, db)
    return Model(params, agents)
end

"""
    create_agents(params::ParameterSet, db::SQLite.DB)

Create a vector of agents with the specifications from a `ParameterSet` and a database.
"""
function create_agents(params::ParameterSet, db::SQLite.DB)
    agent_parties = reduce(vcat, fill.(params.parties, params.group_size))
    opinions = Dict(
        party => filter(p -> p.party_shorthand == party, party_opinions_view(db)).position
        for party in params.parties
    )
    agent_opinions = [Float64.(opinions[party]) for party in agent_parties]
    agents = [
        Agent(i, p, deepcopy(o))
        for (i, (p, o)) in enumerate(zip(agent_parties, agent_opinions))
    ]
    return agents
end

"""
    Meeting

When running a `Model`, each step is performed in a `Meeting`.
"""
mutable struct Meeting
    participants::Array{Agent}
    function Meeting(model::Model, comb::AbstractArray)
        participants = filter(a -> a.party in comb, model.agents)
        new(participants)
    end
end

function Base.show(io::IO, model::Model)
    print(
        """
        Model with $(length(model.agents)) agents from $(length(model.parameter_set.parties)) parties
        """
    )
end


text:
function triangle(i::Integer, parent::ArbField = RealField(64); withinterior = true)
    1 <= i <= 10 || throw(ArgumentError("give 1 ≤ i ≤ 10 not i = $i"))

    angles = [
        (3 // 4, 1 // 3, 1 // 2),
        (2 // 3, 1 // 3, 1 // 2),
        (2 // 3, 1 // 4, 1 // 2),
        (2 // 3, 1 // 3, 1 // 3),
        (3 // 4, 1 // 4, 1 // 3),
        (2 // 3, 1 // 4, 1 // 4),
        (2 // 3, 3 // 4, 3 // 4),
        (2 // 3, 2 // 3, 2 // 3),
        (1 // 2, 2 // 3, 3 // 4),
        (1 // 2, 2 // 3, 2 // 3),
    ]

    r = parent(1 // 8)
    λs =
        ball.(
            parent.([
                "12.400051652843377905",
                "13.744355213213231835",
                "20.571973537984730557",
                "21.309407630190445260",
                "24.456913796299111694",
                "49.109945263284609920",
                "4.261735",
                "5.159146",
                "6.241748",
                "6.777108",
            ]),
            r,
        )

    domain = SphericalTriangle(fmpq.(angles[i]), parent)

    if i in 1:6
        stride = i in [4, 6] ? 2 : 1
        u = SphericalVertexEigenfunction(domain, 1, stride = stride)
    elseif i == 7
        if withinterior
            us = [
                SphericalVertexEigenfunction(domain, 1, stride = 2),
                SphericalVertexEigenfunction(domain, 2, stride = 1),
                SphericalVertexEigenfunction(domain, 3, stride = 1),
                SphericalInteriorEigenfunction(domain, stride = 2),
            ]
            u = SphericalCombinedEigenfunction(domain, us, [1, 1, 1, 4])
        else
            us = [
                SphericalVertexEigenfunction(domain, 1, stride = 2),
                SphericalVertexEigenfunction(domain, 2, stride = 1),
                SphericalVertexEigenfunction(domain, 3, stride = 1),
            ]
            u = SphericalCombinedEigenfunction(domain, us, [1, 1, 1])
        end
    elseif i == 8
        if withinterior
            us = [
                [SphericalVertexEigenfunction(domain, i, stride = 2) for i = 1:3]
                SphericalInteriorEigenfunction(domain, stride = 1)
            ] #6
            u = SphericalCombinedEigenfunction(domain, us, [1, 1, 1, 4])
        else
            us = [[SphericalVertexEigenfunction(domain, i, stride = 2) for i = 1:3]]
            u = SphericalCombinedEigenfunction(domain, us, [1, 1, 1])
        end
    elseif i == 9
        if withinterior
            us = [
                SphericalVertexEigenfunction(domain, 2),
                SphericalVertexEigenfunction(domain, 3),
                SphericalInteriorEigenfunction(domain),
            ]
            u = SphericalCombinedEigenfunction(domain, us, [1, 1, 4])
        else
            us = [
                SphericalVertexEigenfunction(domain, 2),
                SphericalVertexEigenfunction(domain, 3),
            ]
            u = SphericalCombinedEigenfunction(domain, us, [1, 1])
        end
    elseif i == 10
        if withinterior
            us = [
                SphericalVertexEigenfunction(domain, 2),
                SphericalVertexEigenfunction(domain, 3),
                SphericalInteriorEigenfunction(domain, stride = 2),
            ]
            u = SphericalCombinedEigenfunction(domain, us, [1, 1, 4])
        else
            us = [
                SphericalVertexEigenfunction(domain, 2),
                SphericalVertexEigenfunction(domain, 3),
            ]
            u = SphericalCombinedEigenfunction(domain, us, [1, 1])
        end
    end

    (domain, u, λs[i])
end


text:
using ExtractPSF, OffsetArrays
using Test

bead_sigma = (1.0,2.0,3.0)
roi_radii = 2 .* ceil.(Int, bead_sigma)
bias = 100
bead_shift = (0,0,0)
noise_max = 3
bead_max = 100.0 #bead peak intensity
imgsz = (50,50,20)

#Case: first bead's roi isn't fully contained, so it should be excluded
bead_ctrs = [(9,10,6), (30, 10, 7)]

fullimg = ExtractPSF.fake_bead_sample(imgsz, bead_sigma, bead_shift, bead_ctrs, bead_max, noise_max, bias)

#rois is an array of OffsetArrays, ctrs is an array of coords (x,y,z)
rois, ctrs = find_beads(fullimg, roi_radii, bias)

@test length(rois) == 1
@test ctrs[1] == bead_ctrs[2]
@test all(rois[1] .== OffsetArray(view(fullimg, 28:32, 6:14, 1:13), -2:2, -4:4, -6:6))

#Case: bead's peak is at saturation, so should be excluded
fullimg[bead_ctrs[2]...] = typemax(eltype(fullimg))
rois, ctrs = find_beads(fullimg, roi_radii, bias)
@test length(rois) == 0

#Case: bead ROIs overlap, so both should be excluded
bead_ctrs = [(26,10,6), (30, 10, 7)]
fullimg = ExtractPSF.fake_bead_sample(imgsz, bead_sigma, bead_shift, bead_ctrs, bead_max, noise_max, bias)
rois, ctrs = find_beads(fullimg, roi_radii, bias)
@test length(rois) == 0

max_stars_repo_name: HolyLab/ExtractPSF

max_forks_repo_name: HolyLab/ExtractPSF


text:
# Version management

# NVML.driver_version() wrongly reports the forward compatible version,
# so we record the system libcuda version when we initialize the library.
const _system_version = Ref{VersionNumber}()

"""
    system_version()

Returns the latest version of CUDA supported by the system driver.
"""
function system_version()
    libcuda()   # initializes _system_version
    _system_version[]
end

"""
    version()

Returns the latest version of CUDA supported by the loaded driver.
"""
function version()
    version_ref = Ref{Cint}()
    cuDriverGetVersion(version_ref)
    major, ver = divrem(version_ref[], 1000)
    minor, patch = divrem(ver, 10)
    return VersionNumber(major, minor, patch)
end

"""
    release()

Returns the CUDA release part of the version as returned by [`version`](@ref).
"""
release() = VersionNumber(version().major, version().minor)

"""
    runtime_version()

    Returns the CUDA Runtime version.
"""
function runtime_version()
    version_ref = Ref{Cint}()
    @ccall libcudart().cudaRuntimeGetVersion(version_ref::Ptr{Cint})::CUresult
    major, ver = divrem(version_ref[], 1000)
    minor, patch = divrem(ver, 10)
    return VersionNumber(major, minor, patch)
end


text:
function areequal(theory::Vector{Rule}, exprs...;
    timeout=0, sizeout=2^14, matchlimit=5000, mod=@__MODULE__)
    G = EGraph(exprs[1])
    areequal(G, theory, exprs...;
        timeout=timeout, matchlimit=matchlimit, sizeout=sizeout, mod=mod)
end

function areequal(G::EGraph, t::Vector{Rule}, exprs...;
    timeout=0, sizeout=2^14, matchlimit=5000, mod=@__MODULE__)
    @log "Checking equality for " exprs
    if length(exprs) == 1; return true end

    ids = []
    for i ∈ exprs
        ec = addexpr!(G, i)
        push!(ids, ec.id)
    end

    # rebuild!(G)

    @log "starting saturation"

    alleq = () -> (all(x -> in_same_set(G.U, ids[1], x), ids[2:end]))

    report = saturate!(G, t; timeout=timeout, matchlimit=matchlimit,
        sizeout=sizeout, stopwhen=alleq, mod=mod)

    alleq()
end

import ..gettheory

macro areequal(theory, exprs...)
    t = gettheory(theory, __module__; compile=false)
    areequal(t, exprs...; mod=__module__)
end

macro areequalg(G, theory, exprs...)
    t = gettheory(theory, __module__; compile=false)
    areequal(getfield(__module__, G), t, exprs...; mod=__module__)
end


text:
export timeDomainFWIparam, simulateTimeDomainFWI,getTimeDomainFWIParam,getData,getRickerFunction

mutable struct timeDomainFWIparam <: ForwardProbType
    gamma					:: Vector{Float64}     # attenuation
    Sources					:: Union{Vector{Float64},SparseMatrixCSC,Array{Float64,2}}   # Sources
    Receivers				:: Union{Vector{Float64},SparseMatrixCSC,Array{Float64,2}}
	ReceiverMask			:: Array{Int8}
	Mesh      				:: RegularMesh
	ricker					:: Array
	dt						:: Float64
	T 						:: Float64
end

function getTimeDomainFWIParam(gamma,Sources,Receivers,ReceiverMask,Mesh,ricker,dt,T,workerList::Array{Int64}=[])
if workerList==[]
	ActualWorkers = workers();
else
	ActualWorkers = intersect(workerList,workers());
	if length(ActualWorkers)<length(workerList)
		warn("FWI: workerList included indices of non-existing workers.")
	end
end
numWorkers = length(ActualWorkers);
println("We have ",numWorkers," workers for time FWI");
pFor   = Array{RemoteChannel}(undef,numWorkers);
SourcesSubInd = Array{Array{Int64,1}}(undef,numWorkers);
i = 1; nextidx() = (idx=i; i+=1; idx)
nsrc  = size(Sources,2);

# send out jobs
@sync begin
	for p=ActualWorkers
		@async begin
			while true
				idx = nextidx()
				if idx > numWorkers
					break
				end
				I_k = getSourcesIndicesOfKthWorker(numWorkers,idx,nsrc)
				SourcesSubInd[idx] = I_k;
				pFor[idx] = initRemoteChannel(timeDomainFWIparam,p, gamma,Sources[:,I_k], Receivers, ReceiverMask[:,I_k], Mesh, ricker,dt,T);
				wait(pFor[idx]);
			end
		end
	end
end
return pFor, SourcesSubInd# Array of Remote Refs
end

function getRickerFunction(T,dt,fm)
	Td = 2*sqrt(6)/(pi*fm);
	nt = convert(Int64,round(T/dt));
	t  = dt:dt:nt*dt;
	ricker = (1.0.-(2*pi^2*fm^2).*((t.-3*Td).^2)).*exp.(-((t.-3*Td).^2*(pi*fm).^2));
	zeroOffset = convert(Int64,round(((3*Td)./dt) + 1));
	widthOfHalfRicker = round(Int64,(sqrt(8)/(pi*fm))/dt);
	return ricker,zeroOffset,widthOfHalfRicker;
end

function getData(m,pFor::timeDomainFWIparam,doClear::Bool=false)
Mesh       	= pFor.Mesh
gamma   	= pFor.gamma
Q       	= pFor.Sources
P       	= pFor.Receivers
dt 			= pFor.dt
T           = pFor.T;
rickFunc    = pFor.ricker
RecMask     = pFor.ReceiverMask

nrec  		= size(P,2)
nsrc  		= size(Q,2)

D = Array{Array{Float32,2}}(undef, nsrc);

nt = convert(Int64,round(T/dt));

N = prod(Mesh.n.+1);

u_ip1 = zeros(N);
u_im1 = zeros(N);
u_i   = zeros(N);
temp = zeros(N);

## m is assumed to be slowness squared here.
vsq = (dt^2)./m;
gammadt = gamma.*(0.5*dt);
tt = time_ns();
for k = 1:nsrc
	println("At source ",k," out of ",nsrc);
	q = Vector(Q[:,k]).*vsq;
	Dk = calcTrace(Mesh,N,nt,P,vsq,gammadt,q,u_ip1,u_im1,u_i,temp,rickFunc);
	D[k] = vec(RecMask[:,k]).*Dk;
	u_ip1[:] .= 0.0;
	u_im1[:] .= 0.0;
	u_i[:] .= 0.0;
end
println("time domain get data took ",(tt-time_ns())/1e+9," seconds.");

return D,pFor
end

function calcTrace(Mesh::RegularMesh,n::Int64,nt::Int64,P::SparseMatrixCSC,vsq::Array{Float64,1},gammadt::Array{Float64,1},q::Array{Float64,1},u_ip1::Array{Float64,1},u_im1::Array{Float64,1},u_i::Array{Float64,1},temp::Array{Float64,1},rickFunc::Array{Float64})
Dk    = zeros(Float32,size(P,2),nt);
# for i = 2:nt-1
	## Code with 3 u's
	# multOpNeumann!(Mesh,u_i,temp,Lap2DStencil);
	# u_ip1 = (-vsq.*temp + 2.0*u_i - u_im1 + gammadt.*u_im1+ rickFunc[i]*q)./(1.0+gammadt);
	# Dk[:,i+1] = P'*u_ip1;
	# t = u_im1;
	# u_im1 = u_i;
	# u_i = u_ip1;
	# u_ip1 = t;
# end

for i = 2:nt-1
	### Code with 3 u's
	multOpNeumann!(Mesh,u_i,temp,Lap2DStencil);
	rfi = rickFunc[i];
	@simd for j=1:n
		@inbounds temp[j] = -vsq[j]*temp[j] .+ rfi*q[j]
	end
	@simd for j=1:n
		@inbounds u_ip1[j] = (temp[j] + 2.0*u_i[j] - u_im1[j] + gammadt[j]*u_im1[j] )/(1.0+gammadt[j]);
	end
	Dk[:,i+1] = P'*u_ip1;
	t = u_im1;
	u_im1 = u_i;
	u_i = u_ip1;
	u_ip1 = t;
end

return Dk;
end


text:
function _nanremove_withnans()
    testarr = [1.0 1.0 NaN; NaN NaN 2.0; 3.0 3.0 3.0; NaN 4.0 4.0; 5.0 5.0 5.0]

    RHEOS.nanremove(testarr) == [3.0 3.0 3.0; 5.0 5.0 5.0]
end
@test _nanremove_withnans()

function _nanremove_nonans()
    testarr = [1.0 1.0 1.0; 2.0 2.0 2.0; 3.0 3.0 3.0; 4.0 4.0 4.0; 5.0 5.0 5.0]

    RHEOS.nanremove(testarr) == [1.0 1.0 1.0; 2.0 2.0 2.0; 3.0 3.0 3.0; 4.0 4.0 4.0; 5.0 5.0 5.0]
end
@test _nanremove_nonans()

function _importcsv_t_withnans()
    fildir = joinpath(@__DIR__, "testdata", "datawithnans.csv")

    rheodata1 = importcsv(fildir; t_col = 1, σ_col = 2, ϵ_col = 3)
    test1 = (rheodata1.t == [5.0, 6.0] && rheodata1.σ == [10.0, 11.0] && rheodata1.ϵ == [15.0, 16.0])

    rheodata2 = importcsv(fildir; t_col = 3, σ_col = 1, ϵ_col = 2)
    test2 = (rheodata2.t == [15.0, 16.0] && rheodata2.σ == [5.0, 6.0] && rheodata2.ϵ == [10.0, 11.0])

    rheodata3 = importcsv(fildir; t_col = 2, σ_col = 3, ϵ_col = 1)
    test3 = (rheodata3.t == [10.0, 11.0] && rheodata3.σ == [15.0, 16.0] && rheodata3.ϵ == [5.0, 6.0])

    test1 && test2 && test3
end
@test _importcsv_t_withnans()

function _importcsv_ω_withnans()
    fildir = joinpath(@__DIR__, "testdata", "datawithnans.csv")

    rheodata1 = importcsv(fildir; ω_col = 1, Gp_col = 2, Gpp_col = 3)
    test1 = (rheodata1.ω == [5.0, 6.0] && rheodata1.Gp == [10.0, 11.0] && rheodata1.Gpp == [15.0, 16.0])

    rheodata2 = importcsv(fildir; ω_col = 3, Gp_col = 1, Gpp_col = 2)
    test2 = (rheodata2.ω == [15.0, 16.0] && rheodata2.Gp == [5.0, 6.0] && rheodata2.Gpp == [10.0, 11.0])

    rheodata3 = importcsv(fildir; ω_col = 2, Gp_col = 3, Gpp_col = 1)
    test3 = (rheodata3.ω == [10.0, 11.0] && rheodata3.Gp == [15.0, 16.0] && rheodata3.Gpp == [5.0, 6.0])

    test1 && test2 && test3
end
@test _importcsv_ω_withnans()

function _importcsv_t_partial()
    fildir = joinpath(@__DIR__, "testdata", "datapartial.csv")

    rheodata1 = importcsv(fildir; t_col = 1)
    test1 = (rheodata1.t == [1.0, 2.0, 3.0, 4.0, 5.0])

    rheodata2 = importcsv(fildir; t_col = 1, σ_col = 2)
    test2 = (rheodata2.t == [1.0, 2.0, 3.0, 4.0, 5.0] && rheodata2.σ == [10.0, 20.0, 30.0, 40.0, 50.0])

    rheodata3 = importcsv(fildir; t_col = 1, ϵ_col = 1)
    test3 = (rheodata3.t == [1.0, 2.0, 3.0, 4.0, 5.0] && rheodata3.ϵ == [10.0, 20.0, 30.0, 40.0, 50.0])

    test1 && test2 && test3
end
@test _importcsv_t_withnans()

function _exportcsv_timestress()
    fildir = joinpath(@__DIR__, "testdata", "datapartial.csv")

    rheotimedataIN = importcsv(fildir; t_col = 1, σ_col = 2)

    # default column ordering for full time data is (t, σ)
    testdir = joinpath(@__DIR__, "testdata", "testloop.csv")
    exportcsv(rheotimedataIN, testdir)
    rheotimedataOUT = importcsv(testdir; t_col = 1, σ_col = 2)

    rheotimedataIN.σ == rheotimedataOUT.σ && rheotimedataIN.ϵ == rheotimedataOUT.ϵ && rheotimedataIN.t == rheotimedataOUT.t
end
@test _exportcsv_timestress()

function _exportcsv_timestrain()
    fildir = joinpath(@__DIR__, "testdata", "datapartial.csv")

    rheotimedataIN = importcsv(fildir; t_col = 1, ϵ_col = 2)

    # default column ordering for full time data is (t, ϵ)
    testdir = joinpath(@__DIR__, "testdata", "testloop.csv")
    exportcsv(rheotimedataIN, testdir)
    rheotimedataOUT = importcsv(testdir; t_col = 1, ϵ_col = 2)

    rheotimedataIN.σ == rheotimedataOUT.σ && rheotimedataIN.ϵ == rheotimedataOUT.ϵ && rheotimedataIN.t == rheotimedataOUT.t
end
@test _exportcsv_timestrain()

function _exportcsv_timefull()
    fildir = joinpath(@__DIR__, "testdata", "datanonans.csv")

    rheotimedataIN = importcsv(fildir; t_col = 1, σ_col = 2, ϵ_col = 3)

    # default column ordering for full time data is (t, σ, ϵ)
    testdir = joinpath(@__DIR__, "testdata", "testloop.csv")
    exportcsv(rheotimedataIN, testdir)
    rheotimedataOUT = importcsv(testdir; t_col = 1, ϵ_col = 3, σ_col = 2)

    rheotimedataIN.σ == rheotimedataOUT.σ && rheotimedataIN.ϵ == rheotimedataOUT.ϵ && rheotimedataIN.t == rheotimedataOUT.t
end
@test _exportcsv_timefull()

function _exportcsv_freqfull()
    fildir = joinpath(@__DIR__, "testdata", "datanonans.csv")

    rheofreqdataIN = importcsv(fildir; ω_col = 1, Gp_col = 2, Gpp_col = 3)

    # default column ordering for full time data is (ω, Gp, Gpp)
    testdir = joinpath(@__DIR__, "testdata", "testloop.csv")
    exportcsv(rheofreqdataIN, testdir)
    rheofreqdataOUT = importcsv(testdir; ω_col = 1, Gp_col = 2, Gpp_col = 3)

    rheofreqdataIN.ω == rheofreqdataOUT.ω && rheofreqdataIN.Gp == rheofreqdataOUT.Gp && rheofreqdataIN.Gpp == rheofreqdataOUT.Gpp
end
@test _exportcsv_freqfull()

# function _saveload_timefull()
#     fildir = joinpath(@__DIR__, "testdata", "datanonans.csv")

#     rheotimedataIN = importcsv(fildir; t_col = 1, σ_col = 2, ϵ_col = 3)

#     # default column ordering for full time data is (σ, ϵ, t)
#     testdir = joinpath(@__DIR__, "testdata", "testbinaryloop.jld2")
#     savedata(rheotimedataIN, testdir)
#     rheotimedataOUT = loaddata(testdir)

#     rheotimedataIN.σ == rheotimedataOUT.σ && rheotimedataIN.ϵ == rheotimedataOUT.ϵ && rheotimedataIN.t == rheotimedataOUT.t
# end
# @test _saveload_timefull()

# function _saveload_freqfull()
#     fildir = joinpath(@__DIR__, "testdata", "datanonans.csv")

#     rheofreqdataIN = importcsv(fildir; ω_col = 1, Gp_col = 2, Gpp_col = 3)

#     # default column ordering for full time data is (Gp, Gpp, ω)
#     testdir = joinpath(@__DIR__, "testdata", "testbinaryloop.jld2")
#     savedata(rheofreqdataIN, testdir)
#     rheofreqdataOUT = loaddata(testdir)

#     rheofreqdataIN.ω == rheofreqdataOUT.ω && rheofreqdataIN.Gp == rheofreqdataOUT.Gp && rheofreqdataIN.Gpp == rheofreqdataOUT.Gpp
# end
# @test _saveload_freqfull()


text:
module G

using ValidatedNumerics;
include("SubDiff.jl")
include("sti.jl")

function equation(solver)
    sti.STI(solver.system.a*solver.previous.intervals) - solver.system.b.sti
end

function initialConditions(system)
    aMatrixMid = map(x -> mid(x), system.a)
    leftMatrix = sti.constituentMatrix(aMatrixMid)
    rightVector = system.b.sti
    initialInNumbers = leftMatrix \ rightVector
    initialInNumbers
end

function subDifferential(solver)
    identityMatrix = eye(solver.system.size)
    loSubDiff = map(
        index -> calulateLoSubdifferentialRow(slicedim(solver.system.a, 1, index), solver.previous.sti, slicedim(identityMatrix, 1, index)),
        range(1, solver.system.size)
    )
    hiSubDiff = map(
        index -> calulateHiSubdifferentialRow(slicedim(solver.system.a, 1, index), solver.previous.sti, slicedim(identityMatrix, 1, index)),
        range(1, solver.system.size)
    )
    vcat(loSubDiff..., hiSubDiff...)
end

function calulateLoSubdifferentialRow{T}(aVector :: Array{Interval{T}}, xVector :: Array{T}, _ :: Array{T})
    calulateSubdifferentialRowPart(SubDiff.productLoSubDifferential, -1, aVector, xVector, _)
end

function calulateHiSubdifferentialRow{T}(aVector :: Array{Interval{T}}, xVector :: Array{T}, _ :: Array{T})
    calulateSubdifferentialRowPart(SubDiff.productHiSubDifferential, 1, aVector, xVector, _)
end

function calulateSubdifferentialRowPart{T}(partFunction, sumFactor, aVector :: Array{Interval{T}}, xVector :: Array{T}, _ :: Array{T})
    systemSize = size(xVector, 1)
    systemSizeHalf = Int(systemSize/2)
    mapped =  map(
        index -> partFunction(aVector[index], xVector[index], xVector[index + systemSizeHalf]),
        range(1, systemSizeHalf)
    )
    calculatedSum = sumFactor * mapped
    stiedCalculatedSum = sti.STI(calculatedSum)
    transpose(stiedCalculatedSum)
end
end


text:
# This file is a part of Devices.jl, licensed under the MIT License (MIT).

import Compat.Test
Test.@testset "Package Devices" begin
end


text:
# Test that https://www.google.com successfully connects
curl = curl_easy_init()
curl == C_NULL && error("curl_easy_init() failed")

# Setup the callback function to recv data
function curl_write_cb(curlbuf::Ptr{Cvoid}, s::Csize_t, n::Csize_t, p_ctxt::Ptr{Cvoid})
    sz = s * n
    data = Array{UInt8}(undef, sz)
    ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), data, curlbuf, sz)
    return sz::Csize_t
end

@testset "SSL verify" begin
    # Set up the write function to consume the curl output so we don't see it in the
    # test output
    c_curl_write_cb = @cfunction(
        curl_write_cb,
        Csize_t,
        (Ptr{Cvoid}, Csize_t, Csize_t, Ptr{Cvoid})
    )
    curl_easy_setopt(curl, CURLOPT_WRITEFUNCTION, c_curl_write_cb)

    # Set up our SSL options
    curl_easy_setopt(curl, CURLOPT_URL, "https://www.google.com")
    curl_easy_setopt(curl, CURLOPT_USE_SSL, CURLUSESSL_ALL)
    curl_easy_setopt(curl, CURLOPT_SSL_VERIFYHOST, 2)
    curl_easy_setopt(curl, CURLOPT_SSL_VERIFYPEER, 1)
    curl_easy_setopt(curl, CURLOPT_CAINFO, LibCURL.cacert)

    @testset "SSL Success" begin
        res = curl_easy_perform(curl)
        @test res == CURLE_OK
    end

end

curl_easy_cleanup(curl)


text:
module TestDocumentation

using Base.Test
using DataArrays
using DataFrames

println("\n Running documentation tests\n")

using TimeData
using Dates

fileName = joinpath(Pkg.dir("TimeData"), "data/logRet.csv")
tm = TimeData.readTimedata(fileName)[1:10, 1:4]

tm[Date(2012, 1, 4):Date(2012, 1, 10), 1:2]

tm[3:8, 2:3]

tm[3:8, [:ABT, :MMM]]

newTm = 2.*tm

mean(tm, 1)

rowmeans(tm)

vals = rand(4, 3);
dats = [Date(2013, 7, ii) for ii=1:4];
nams = [:A, :B, :C];
valsDf = composeDataFrame(vals, nams);

tm = Timematr(valsDf, dats)

td = Timedata(vals, nams, dats)
td = Timedata(vals, nams)
td = Timedata(vals, dats)
td = Timedata(vals)

typeof(valsDf[:, 1])
typeof(td[:, 1])

typeof(valsDf[1, 1])
typeof(td[1, 1])

## empty instance
typeof(td[4:3, 5:4])

## indexing by numeric indices
tmp = tm[2:3]
tmp = tm[1:3, 1:2]
tmp = tm[2, :]
tmp = tm[2]
tmp = tm[1:2, 2]
tmp = tm[3, 3]

## indexing with column names as symbols
tmp = tm[:A]
tmp = tm[2, [:A, :B]]

## logical indexing
logicCol = [true, false, true]
logicRow = repmat([true, false], 2, 1)[:]
tmp = tm[logicCol]
tmp = tm[logicRow, logicCol]
tmp = tm[logicRow, :]

## indexing with Dates
DatesToFind = [Date(2013, 7, ii) for ii=2:3]
tmp = tm[DatesToFind]
tm[Date(2013,7,1):Date(2013,7,3)]
tm[Date(2013,7,2):Date(2013,7,3), :B]
tm[Date(2013,7,3):Date(2013,7,12), [true, false, false]]

## returning the first value only
get(tm, 1, 1)

## returning all values as Array{Any,2}
kk = get(tm)
isa(kk, Array{Any})

filePath = joinpath(Pkg.dir("TimeData"), "data", "logRet.csv");
tm = readTimedata(filePath)
tm[1:5, 1:4]

tm

## str(tm) # uncomment for execution

#   writeTimedata("data/logRet2.csv", tm) # uncomment for execution

typeof(tm .+ tm)
typeof(tm .> 0.5)

tm[1:3, 1:3] .> 0.5
exp(tm[1:3, 1:3])
round(tm[1:3, 1:3], 2)

filePath = joinpath(Pkg.dir("TimeData"), "data", "logRet.csv");
tm = readTimedata(filePath)

columnsMeetingCondition = chkVars(x -> minimum(x, 1)[1, 1] .< -25, eachvar(tm)) |>
                          x -> asArr(x, Bool, false) |>
                          x -> tm[x[:]]

minimum(columnsMeetingCondition, 1)

end


text:
# build pdf version of Documenter's docs on 64-bit Linux
if Sys.ARCH === :x86_64 && Sys.KERNEL === :Linux
    cd(joinpath(@__DIR__, "..", "..")) do
        cmd = `$(Base.julia_cmd()) --project=docs/pdf/`
        @test success(`$(cmd) -e 'using Pkg; Pkg.instantiate()'`)
        @test success(`$(cmd) docs/pdf/make.jl --verbose`)
        # deploy only from Julia v1.0.X
        if VERSION.major == 1 && VERSION.minor == 1
            @test success(`$(cmd) docs/pdf/deploy.jl`)
        end
    end
end


text:
# a simplified PSNR to get rid of ImageQualityIndexes and ImageFiltering dependency
assess_psnr(ref::AbstractArray{<:Color3}, x::AbstractArray{<:Color3}) =
    assess_psnr(channelview(RGB.(ref)), channelview(RGB.(x)), 1.0)

assess_psnr(ref::AbstractArray{<:ColorTypes.Transparent3}, x::AbstractArray{<:ColorTypes.Transparent3}) =
    assess_psnr(channelview(ARGB.(ref)), channelview(ARGB.(x)), 1.0)

assess_psnr(ref::AbstractArray{<:AbstractGray}, x::AbstractArray{<:AbstractGray}) =
    assess_psnr(channelview(ref), channelview(x), 1.0)

assess_psnr(ref::AbstractArray{<:Real}, x::AbstractArray{<:Real}, peakval::Real) =
    20log10(peakval) - 10log10(_mse(float.(ref), float.(x)))

function _mse(x, y)
    @assert length(x) == length(y)
    _euclidean(x, y)/(length(x))
end
_euclidean(x, y) = sqrt(sum((x - y) .^ 2))


text:
module UnionFindAlgo

export UnionFind, union!, find, isconnected

"""
Union find with the *Weighted quick union* (WQU)
algorithm.

- `id`: the parent node's id
- `sz`: size of the tree
"""
mutable struct UnionFind
    id::Vector{Int}
    sz::Vector{Int}
    ncomp::Int
    UnionFind(n::Int) = new(collect(1:n), fill(1, n), n)
end

function UnionFind(data::NamedTuple)
    uf = UnionFind(data.nv)
    edges = data.edges
    for (p, q) in edges
        union!(uf, p, q)
    end
    uf
end

function union!(uf, p, q)
    i = find(uf, p)
    j = find(uf, q)
    i == j && return
    # weighted quick-union: always connect the root
    # of smaller tree to the root of the larger tree
    if (uf.sz[i] < uf.sz[j])
        uf.id[i] = j
        uf.sz[j] += uf.sz[i]
    else
        uf.id[j] = i
        uf.sz[i] += uf.sz[j]
    end
    uf.ncomp -= 1
    return nothing
end

"""
    find(uf, p)

Find the root of p. The component in WQUPC is defined by
the root of each node.

In theory path compression can help further flatten
the tree. There are two ways to conduct path compression:

1. find the root for `p` and then assign all nodes on the path
   from `p` to root to the newly found root. This requires running
   the `find` twice.
2. Directly change parent to grandparent, this will halve the
   height in the subtree. In code: `uf.id[i] = uf.id[uf.id[i]]`

In practice, the performance gain from path compression is not
significant. On `largeUF.txt`(V=1e6, E=2e6) approach (2) gives
about 1% improvement. The difference might be larger for larger
dataset.
"""
function find(uf, p)
    i = uf.id[p]
    while uf.id[i] != i
        # path compression: set parent to grandparent
        # this trick halves the subtree each time
        uf.id[i] = uf.id[uf.id[i]]
        # go to the parent node
        i = uf.id[i]
    end
    i
end

isconnected(uf, p, q) = find(uf, p) == find(uf, q)

end

max_stars_repo_name: Jetafull/AlgorithmsInJulia

max_issues_repo_name: Jetafull/algorithms-in-julia

max_forks_repo_name: Jetafull/algorithms-in-julia


text:
using AMQPClient
const VIRTUALHOST = "/"
const HOST = "127.0.0.1"

function send()
    # 1. Create a connection to the localhost or 127.0.0.1 of virtualhost '/'
    connection(; virtualhost=VIRTUALHOST, host=HOST) do conn
        # 2. Create a channel to send messages
        channel(conn, AMQPClient.UNUSED_CHANNEL, true) do chan
            queue = "task_queue"
            # 3. Configure the queue
            success, queue_name, message_count, consumer_count = queue_declare(chan, queue, durable=true)

            # 4. Prepare the message text
            if length(Base.ARGS) >= 1
                received = Base.ARGS[1]
            else
                received = "Hello World"
            end

            # 5. Prepare the payload
            data = convert(Vector{UInt8}, codeunits(received))
            msg = Message(data, content_type="text/plain", delivery_mode=PERSISTENT)

            # 6. Send the payload
            basic_publish(chan, msg; exchange="", routing_key=queue)
            println("Message sent: $received")
        end
    end
end

send()

max_forks_repo_name: tomashq/rabbitmq-tutorials