text:
module SkillRating

    const q = log(10)/400

    g(rd) = 1/sqrt(1 + 3q^2 * (rd^2)/π^2)
    Es(r, r2, rd2) = 1/(1+10^(-g(rd2) * (r - r2) / 400))
    d2(r, rjs, RDjs) = 1/((q^2) * sum(g.(RDjs).^2 .* Es.(r, rjs, RDjs).* (1 .- Es.(r, rjs, RDjs))))

    """
        glicko_rating(r, RD, rjs::Vector, RDjs::Vector, sjs::Vector)
        glicko_rating(r, RD, rj::T, RDj::T, sj::T) where T <: Real

    # Arguments:
    * `r`: Rating of player
    * `RD`: Rating deviation of player
    * `rjs`: Ratings of oppenents within rating period.
    * `RDjs`: Rating devations of oppenents within rating period.
    * `sjs`: Game outcomes from the perspective of the "player": 1 - won, 0.5 - draw, 0 - lost.

    Calculate a new glicko rating for games within a rating period or for a single game.

    See: http://www.glicko.net/glicko/glicko.pdf and Glickman, Mark E., "Parameter estimation
    in large dynamic paired comparison experiments" (1999) Applied Statistics, 48, 377-394
    (http://www.glicko.net/research/glicko.pdf)
    """
    glicko_rating(r, RD, rjs::Vector, RDjs::Vector, sjs::Vector) = r + q/(1/RD^2 + 1/d2(r, rjs, RDjs)) * sum(g.(RDjs) .* (sjs .- Es.(r, rjs, RDjs)))

    """
        glicko_RD(r, RD, rjs::Vector, RDjs::Vector, sjs::Vector)
        glicko_RD(r, RD, rj::T, RDj::T, sj::T) where T <: Real

    # Arguments:
    * `r`: Rating of player
    * `RD`: Rating deviation of player
    * `rjs`: Ratings of oppenents within rating period.
    * `RDjs`: Rating devations of oppenents within rating period.
    * `sjs`: Game outcomes from the perspective of the "player": 1 - won, 0.5 - draw, 0 - lost.

    Calculate a new glicko rating deviation for games within a rating period  or for a single game.

    See: http://www.glicko.net/glicko/glicko.pdf and
    Glickman, Mark E., "Parameter estimation in large dynamic paired comparison experiments" (1999) Applied Statistics, 48, 377-394
    (http://www.glicko.net/research/glicko.pdf)
    """
    glicko_RD(r, RD, rjs::Vector, RDjs::Vector, sjs::Vector) = sqrt(1/(1/RD^2 + 1/d2(r, rjs, RDjs)))

    glicko_rating(r, RD, rj::T, RDj::T, sj::T) where T <: Real = glicko_rating(r, RD, [rj], [RDj], [sj])
    glicko_RD(r, RD, rj::T, RDj::T, sj::T) where T <: Real = glicko_RD(r, RD, [rj], [RDj], [sj])

    """
        glicko_RD_increase(rd; c = 63.2, maxrd=350)

    # Arguments:
    * `RD`: Rating deviation at end of rating period.
    * `c`: Increase in uncertainty between rating periods.
    * `maxRD`: Rating deviation of a player who does not play.

    Calculate rating deviation after the end of a rating period.

    See: http://www.glicko.net/glicko/glicko.pdf and
    Glickman, Mark E., "Parameter estimation in large dynamic paired comparison experiments" (1999) Applied Statistics, 48, 377-394
    (http://www.glicko.net/research/glicko.pdf)
    """
    glicko_RD_increase(RD; c = 63.2, maxRD=350) = min(sqrt(RD^2+c^2), maxRD)

    """
        glicko_c(rating_perods_to_max, maxrd, medianrd)

    Calculate c, the increase in uncertainty between rating periods.

    See: http://www.glicko.net/glicko/glicko.pdf and
    Glickman, Mark E., "Parameter estimation in large dynamic paired comparison experiments" (1999) Applied Statistics, 48, 377-394
    (http://www.glicko.net/research/glicko.pdf)
    """
    glicko_c(rating_perods_to_max, maxrd, medianrd) = sqrt((maxrd^2 - medianrd^2)/rating_perods_to_max)

    """
    elo_expected(RA, RB)

    # Arguments:
    * `RA`: Rating of player A
    * `RB`: Rating of player B

    Calculate the expected outcome using glicko.

    See: http://www.glicko.net/glicko/glicko.pdf and
    Glickman, Mark E., "Parameter estimation in large dynamic paired comparison experiments" (1999) Applied Statistics, 48, 377-394
    (http://www.glicko.net/research/glicko.pdf)
    """
    glicko_expected(r1, RD1, r2, RD2) = 1 / (1 + 10^(-g(sqrt(RD1^2 + RD2^2)) * (r1 - r2)/400))

    """
        elo_rating(r, rjs::Vector, sjs::Vector; k = 32)

    # Arguments:
    * `r`: Rating of player
    * `rjs`: Ratings of oppenents within rating period.
    * `sjs`: Game outcomes from the perspective of the "player": 1 - won, 0.5 - draw, 0 - lost.

    Calculate the elo rating change.
    See: https://en.wikipedia.org/wiki/Elo_rating_system#Mathematical_details
    """
    elo_rating(r, rjs::Vector, sjs::Vector; k = 32) = r + k * (sum(sjs) - sum(elo_expected.(r, rjs)))

    """
    elo_expected(RA, RB)

    # Arguments:
    * `RA`: Rating of player A
    * `RB`: Rating of player B

    Calculate the expected outcome using elo.
    See: https://en.wikipedia.org/wiki/Elo_rating_system#Mathematical_details
    """
    elo_expected(RA, RB) = 1 / (1 + 10 ^ ((RB - RA) / 400))

    Eg2(μ, μj, ϕj) = 1/(1+exp(-gg2(ϕj)*(μ - μj)))
    gg2(ϕ) = 1/sqrt(1 + 3*(ϕ^2)/π^2)
    compute_v(μ, μjs, ϕjs) = 1/sum(gg2.(ϕjs).^2 .* Eg2.(μ, μjs, ϕjs) .* (1 .- Eg2.(μ, μjs, ϕjs)))
    compute_Δ(μ, μjs, ϕjs, sjs, v) = v * sum(gg2.(ϕjs) .* (sjs .- Eg2.(μ, μjs, ϕjs)))
    ϕincrease(ϕ, σ) = sqrt(ϕ^2 + σ^2)

    """
        glicko2(μ::Real, ϕ::Real, σ::Real, μjs::Vector{Real}, ϕjs::Vector{Real}, sjs::Vector{Real}; τ=0.05, ϵ=10^-6)

    # Arguments:
    * `μ`: Rating of player
    * `ϕ`: Rating deviation of player
    * `σ`: Rating volatility
    * `μjs`: Ratings of oppenents within rating period
    * `μjs`: Rating deviations within rating period
    * `sjs`: Game outcomes from the perspective of the "player": 1 - won, 0.5 - draw, 0 - lost.
    * `τ=0.5`: smaller values constrain the change in volatility over time, reasonable range between 0.3 and 1.2
    * `ϵ=10^-6`: Convergence bound

    Returns tuple `(μ, ϕ, σ)`: glicko2 rating, rating deviation and rating volatility for games in one rating period.

    See: http://www.glicko.net/glicko/glicko2.pdf and
    Glickman, Mark E., "Dynamic paired comparison models with stochastic variances" (2001),
    Journal of Applied Statistics, 28, 673-689.
    (http://www.glicko.net/research/dpcmsv.pdf)
    """
    function glicko2(μ::Real, ϕ::Real, σ::Real, μjs::Vector{R}, ϕjs::Vector{R}, sjs::Vector{S}; τ=0.5, ϵ=10^-6) where {R <: Real, S <: Real}
        v = compute_v(μ, μjs, ϕjs)
        Δ = compute_Δ(μ, μjs, ϕjs, sjs, v)
        a = log(σ^2)
        f(x) = exp(x)*(Δ^2 - ϕ^2  - v - exp(x)) / (2 * (ϕ^2 + v + exp(x))^2) - (x - a)/τ^2
        A = a
        if Δ^2 > ϕ^2 + v
            B = log(Δ^2 - ϕ^2  - v)
        else
            k = 1
            while (f(a - k*τ) < 0)
                k += 1
            end
            B = a - k*τ
        end
        fA = f(A)
        fB = f(B)
        while abs(B - A) > ϵ
            C = A + (A-B) * fA/ (fB-fA)
            fC = f(C)
            if fC * fB < 0
                A = B
                fA = fB
            else
                fA /= 2
            end
            B = C
            fB = fC
        end
        σnew = exp(A/2)
        ϕincreased = sqrt(ϕ^2 + σnew^2)
        ϕnew = 1/sqrt(1/ϕincreased^2 + 1/v)
        μnew = μ + ϕnew^2 * sum(gg2.(ϕjs) .* (sjs .- Eg2.(μ, μjs, ϕjs)))
        μnew, ϕnew, σnew
    end

    """
        glicko1_to_glicko2(r::real, RD::Real)

    Convert ratings `r` and rating deviations `RD` from glicko1 (Elo-like) scale to glicko 2 scale (μ and ϕ).
    """
    glicko1_to_glicko2(r::Real, RD::Real) = ((r-1500)/173.7178, RD/173.7178)

    export glicko_rating, glicko_RD, glicko_RD_increase, glicko_c, glicko_expected, elo_rating, elo_expected, glicko2
end

max_stars_repo_name: altre/SkillRating

max_issues_repo_name: altre/SkillRating

max_forks_repo_name: altre/SkillRating


text:
"""
Selection of random objects on the screen is realized by rendering an
object id + plus an arbitrary index into the framebuffer.
The index can be used for e.g. instanced geometries.
"""
struct SelectionID{T <: Integer} <: FieldVector{2, T}
    id::T
    index::T
end

function draw_fullscreen(vao_id)
    glBindVertexArray(vao_id)
    glDrawArrays(GL_TRIANGLES, 0, 3)
    glBindVertexArray(0)
    return
end

struct PostprocessPrerender end

function (sp::PostprocessPrerender)()
    glDepthMask(GL_TRUE)
    glDisable(GL_DEPTH_TEST)
    glDisable(GL_BLEND)
    glDisable(GL_CULL_FACE)
    return
end

const PostProcessROBJ = RenderObject{PostprocessPrerender}

mutable struct GLFramebuffer
    resolution::Node{NTuple{2, Int}}
    id::NTuple{2, GLuint}

    color::Texture{RGBA{N0f8}, 2}
    objectid::Texture{Vec{2, GLuint}, 2}
    depth::Texture{GLAbstraction.DepthStencil_24_8, 2}
    position::Texture{Vec4f0, 2}
    normal_occlusion::Texture{Vec4f0, 2}
    ssao_noise::Texture{Vec2f0, 2}

    color_luma::Texture{RGBA{N0f8}, 2}

    postprocess::NTuple{5, PostProcessROBJ}
end

Base.size(fb::GLFramebuffer) = size(fb.color) # it's guaranteed, that they all have the same size

loadshader(name) = joinpath(@__DIR__, "GLVisualize", "assets", "shader", name)

rcpframe(x) = 1f0 ./ Vec2f0(x[1], x[2])

"""
Creates a postprocessing render object.
This will transfer the pixels from the color texture of the Framebuffer
to the screen and while at it, it can do some postprocessing (not doing it right now):
E.g fxaa anti aliasing, color correction etc.
"""
function postprocess(
        color, position, normal_occlusion, ssao_noise,
        objectid, color_luma,
        framebuffer_size
    )

    # SSAO setup
    N_samples = 64
    lerp_min = 0.1f0
    lerp_max = 1.0f0
    kernel = map(1:N_samples) do i
        n = normalize([2.0rand() .- 1.0, 2.0rand() .- 1.0, rand()])
        scale = lerp_min + (lerp_max - lerp_min) * (i / N_samples)^2
        v = Vec3f0(scale * rand() * n)
    end

    # compute occlusion
    shader1 = LazyShader(
        loadshader("postprocessing/fullscreen.vert"),
        loadshader("postprocessing/SSAO.frag"),
        view = Dict(
            "N_samples" => "$N_samples"
        )
    )
    data1 = Dict{Symbol, Any}(
        :position_buffer => position,
        :normal_occlusion_buffer => normal_occlusion,
        :kernel => kernel,
        :noise => ssao_noise,
        :noise_scale => map(s -> Vec2f0(s ./ 4.0), framebuffer_size),
        :projection => Node(Mat4f0(I)),
        :bias => Node(0.025f0),
        :radius => Node(0.5f0)
    )
    pass1 = RenderObject(data1, shader1, PostprocessPrerender(), nothing)
    pass1.postrenderfunction = () -> draw_fullscreen(pass1.vertexarray.id)

    # blur occlusion and combine with color
    shader2 = LazyShader(
        loadshader("postprocessing/fullscreen.vert"),
        loadshader("postprocessing/SSAO_blur.frag")
    )
    data2 = Dict{Symbol, Any}(
        :normal_occlusion => normal_occlusion,
        :color_texture => color,
        :ids => objectid,
        :inv_texel_size => map(t -> Vec2f0(1f0/t[1], 1f0/t[2]), framebuffer_size),
        :blur_range => Node(Int32(2))
    )
    pass2 = RenderObject(data2, shader2, PostprocessPrerender(), nothing)
    pass2.postrenderfunction = () -> draw_fullscreen(pass2.vertexarray.id)

    # calculate luma for FXAA
    shader3 = LazyShader(
        loadshader("postprocessing/fullscreen.vert"),
        loadshader("postprocessing/postprocess.frag")
    )
    data3 = Dict{Symbol, Any}(
        :color_texture => color
    )
    pass3 = RenderObject(data3, shader3, PostprocessPrerender(), nothing)
    pass3.postrenderfunction = () -> draw_fullscreen(pass3.vertexarray.id)

    # perform FXAA
    shader4 = LazyShader(
        loadshader("postprocessing/fullscreen.vert"),
        loadshader("postprocessing/fxaa.frag")
    )
    data4 = Dict{Symbol, Any}(
        :color_texture => color_luma,
        :RCPFrame => lift(rcpframe, framebuffer_size),
    )
    pass4 = RenderObject(data4, shader4, PostprocessPrerender(), nothing)
    pass4.postrenderfunction = () -> draw_fullscreen(pass4.vertexarray.id)

    # draw color buffer
    shader5 = LazyShader(
        loadshader("postprocessing/fullscreen.vert"),
        loadshader("postprocessing/copy.frag")
    )
    data5 = Dict{Symbol, Any}(
        :color_texture => color
    )
    pass5 = RenderObject(data5, shader5, PostprocessPrerender(), nothing)
    pass5.postrenderfunction = () -> draw_fullscreen(pass5.vertexarray.id)

    return (pass1, pass2, pass3, pass4, pass5)
end

function attach_framebuffer(t::Texture{T, 2}, attachment) where T
    glFramebufferTexture2D(GL_FRAMEBUFFER, attachment, GL_TEXTURE_2D, t.id, 0)
end

function GLFramebuffer(fb_size::NTuple{2, Int})
    # First Framebuffer
    render_framebuffer = glGenFramebuffers()
    glBindFramebuffer(GL_FRAMEBUFFER, render_framebuffer)

    color_buffer = Texture(RGBA{N0f8}, fb_size, minfilter = :nearest, x_repeat = :clamp_to_edge)
    objectid_buffer = Texture(Vec{2, GLuint}, fb_size, minfilter = :nearest, x_repeat = :clamp_to_edge)
    position_buffer = Texture(Vec4f0, fb_size, minfilter = :nearest, x_repeat = :clamp_to_edge)
    normal_occlusion_buffer = Texture(Vec4f0, fb_size, minfilter = :nearest, x_repeat = :clamp_to_edge)

    depth_buffer = Texture(
        Ptr{GLAbstraction.DepthStencil_24_8}(C_NULL), fb_size,
        minfilter = :nearest, x_repeat = :clamp_to_edge,
        internalformat = GL_DEPTH24_STENCIL8,
        format = GL_DEPTH_STENCIL
    )

    ssao_noise = Texture(
        [normalize(Vec2f0(2.0rand(2) .- 1.0)) for _ in 1:4, __ in 1:4],
        minfilter = :nearest, x_repeat = :repeat
    )

    attach_framebuffer(color_buffer, GL_COLOR_ATTACHMENT0)
    attach_framebuffer(objectid_buffer, GL_COLOR_ATTACHMENT1)
    attach_framebuffer(position_buffer, GL_COLOR_ATTACHMENT2)
    attach_framebuffer(normal_occlusion_buffer, GL_COLOR_ATTACHMENT3)
    attach_framebuffer(depth_buffer, GL_DEPTH_ATTACHMENT)
    attach_framebuffer(depth_buffer, GL_STENCIL_ATTACHMENT)

    status = glCheckFramebufferStatus(GL_FRAMEBUFFER)
    @assert status == GL_FRAMEBUFFER_COMPLETE

    # Second Framebuffer
    color_luma_framebuffer = glGenFramebuffers()
    glBindFramebuffer(GL_FRAMEBUFFER, color_luma_framebuffer)

    color_luma = Texture(RGBA{N0f8}, fb_size, minfilter=:linear, x_repeat=:clamp_to_edge)
    attach_framebuffer(color_luma, GL_COLOR_ATTACHMENT0)

    @assert status == GL_FRAMEBUFFER_COMPLETE

    glBindFramebuffer(GL_FRAMEBUFFER, 0)
    fb_size_node = Node(fb_size)

    p = postprocess(
        color_buffer, position_buffer, normal_occlusion_buffer, ssao_noise, #occlusion,
        objectid_buffer, color_luma,
        fb_size_node
    )

    return GLFramebuffer(
        fb_size_node,
        (render_framebuffer, color_luma_framebuffer),
        color_buffer, objectid_buffer, depth_buffer,
        position_buffer, normal_occlusion_buffer, ssao_noise,# occlusion,
        color_luma,
        p
    )
end

function Base.resize!(fb::GLFramebuffer, window_size)
    ws = Int.((window_size[1], window_size[2]))
    if ws != size(fb) && all(x-> x > 0, window_size)
        resize_nocopy!(fb.color, ws)
        resize_nocopy!(fb.objectid, ws)
        resize_nocopy!(fb.depth, ws)
        resize_nocopy!(fb.position, ws)
        resize_nocopy!(fb.normal_occlusion, ws)
        resize_nocopy!(fb.color_luma, ws)
        fb.resolution[] = ws
    end
    nothing
end

struct MonitorProperties
    name::String
    isprimary::Bool
    position::Vec{2, Int}
    physicalsize::Vec{2, Int}
    videomode::GLFW.VidMode
    videomode_supported::Vector{GLFW.VidMode}
    dpi::Vec{2, Float64}
    monitor::GLFW.Monitor
end

function MonitorProperties(monitor::GLFW.Monitor)
    name = GLFW.GetMonitorName(monitor)
    isprimary = GLFW.GetPrimaryMonitor() == monitor
    position = Vec{2, Int}(GLFW.GetMonitorPos(monitor)...)
    physicalsize = Vec{2, Int}(GLFW.GetMonitorPhysicalSize(monitor)...)
    videomode = GLFW.GetVideoMode(monitor)
    sfactor = Sys.isapple() ? 2.0 : 1.0
    dpi = Vec(videomode.width * 25.4, videomode.height * 25.4) * sfactor ./ Vec{2, Float64}(physicalsize)
    videomode_supported = GLFW.GetVideoModes(monitor)

    MonitorProperties(name, isprimary, position, physicalsize, videomode, videomode_supported, dpi, monitor)
end

was_destroyed(nw::GLFW.Window) = nw.handle == C_NULL

function GLContext()
    context = GLFW.GetCurrentContext()
    version = opengl_version_number()
    glsl_version = glsl_version_number()
    return GLContext(context, version, glsl_version, unique_context_counter())
end

function ShaderAbstractions.native_switch_context!(x::GLFW.Window)
    GLFW.MakeContextCurrent(x)
end

function ShaderAbstractions.native_context_alive(x::GLFW.Window)
    GLFW.is_initialized() && !was_destroyed(x)
end

function destroy!(nw::GLFW.Window)
    was_current = ShaderAbstractions.is_current_context(nw)
    if !was_destroyed(nw)
        GLFW.SetWindowShouldClose(nw, true)
        GLFW.PollEvents()
        GLFW.DestroyWindow(nw)
        nw.handle = C_NULL
    end
    was_current && ShaderAbstractions.switch_context!()
end

function windowsize(nw::GLFW.Window)
    was_destroyed(nw) && return (0, 0)
    size = GLFW.GetFramebufferSize(nw)
    return (size.width, size.height)
end

function Base.isopen(window::GLFW.Window)
    was_destroyed(window) && return false
    try
        return !GLFW.WindowShouldClose(window)
    catch e
        # can't be open if GLFW is already terminated
        e.code == GLFW.NOT_INITIALIZED && return false
        rethrow(e)
    end
end

ext: jl


text:
module SurfaceWaves

export
    ∂t_uˢ,
    ∂t_vˢ,
    ∂t_wˢ,
    x_curl_Uˢ_cross_U,
    y_curl_Uˢ_cross_U,
    z_curl_Uˢ_cross_U

using Oceananigans.Grids: AbstractGrid

using Oceananigans.Fields
using Oceananigans.Operators

"""
    abstract type AbstractStokesDrift end

Parent type for parameter structs for Stokes drift fields
associated with surface waves
"""
abstract type AbstractStokesDrift end

#####
##### Functions for "no surface waves"
#####

@inline ∂t_uˢ(i, j, k, grid::AbstractGrid{FT}, ::Nothing, time) where FT = zero(FT)
@inline ∂t_vˢ(i, j, k, grid::AbstractGrid{FT}, ::Nothing, time) where FT = zero(FT)
@inline ∂t_wˢ(i, j, k, grid::AbstractGrid{FT}, ::Nothing, time) where FT = zero(FT)

@inline x_curl_Uˢ_cross_U(i, j, k, grid::AbstractGrid{FT}, ::Nothing, U, time) where FT = zero(FT)
@inline y_curl_Uˢ_cross_U(i, j, k, grid::AbstractGrid{FT}, ::Nothing, U, time) where FT = zero(FT)
@inline z_curl_Uˢ_cross_U(i, j, k, grid::AbstractGrid{FT}, ::Nothing, U, time) where FT = zero(FT)

#####
##### Uniform surface waves
#####

"""
    UniformStokesDrift{UZ, VZ, UT, VT} <: AbstractStokesDrift

Parameter struct for Stokes drift fields associated with surface waves.
"""
struct UniformStokesDrift{UZ, VZ, UT, VT} <: AbstractStokesDrift
    ∂z_uˢ :: UZ
    ∂z_vˢ :: VZ
    ∂t_uˢ :: UT
    ∂t_vˢ :: VT
end

addzero(args...) = 0

"""
    UniformStokesDrift(; ∂z_uˢ=addzero, ∂z_vˢ=addzero,
                                    ∂t_uˢ=addzero, ∂t_vˢ=addzero)

Construct a set of functions that describes the Stokes drift field beneath
a uniform surface gravity wave field.
"""
UniformStokesDrift(; ∂z_uˢ=addzero, ∂z_vˢ=addzero, ∂t_uˢ=addzero, ∂t_vˢ=addzero) =
    UniformStokesDrift(∂z_uˢ, ∂z_vˢ, ∂t_uˢ, ∂t_vˢ)

const USD = UniformStokesDrift

@inline ∂t_uˢ(i, j, k, grid, sw::USD, time) = sw.∂t_uˢ(znode(Cell, k, grid), time)
@inline ∂t_vˢ(i, j, k, grid, sw::USD, time) = sw.∂t_vˢ(znode(Cell, k, grid), time)

@inline ∂t_wˢ(i, j, k, grid::AbstractGrid{FT}, sw::USD, time) where FT = zero(FT)

@inline x_curl_Uˢ_cross_U(i, j, k, grid, sw::USD, U, time) =
    @inbounds ℑxzᶠᵃᶜ(i, j, k, grid, U.w) * sw.∂z_uˢ(znode(Cell, k, grid), time)

@inline y_curl_Uˢ_cross_U(i, j, k, grid, sw::USD, U, time) =
    @inbounds ℑyzᵃᶠᶜ(i, j, k, grid, U.w) * sw.∂z_vˢ(znode(Cell, k, grid), time)

@inline z_curl_Uˢ_cross_U(i, j, k, grid, sw::USD, U, time) = @inbounds begin (
    - ℑxzᶜᵃᶠ(i, j, k, grid, U.u) * sw.∂z_uˢ(znode(Face, k, grid), time)
    - ℑyzᵃᶜᶠ(i, j, k, grid, U.v) * sw.∂z_vˢ(znode(Face, k, grid), time) )
end

end # module


text:
using Test, PyCall
using PyCall: f_contiguous, PyBUF_ND_CONTIGUOUS, array_format, npy_initialized,
NoCopyArray, isbuftype, setdata!

pyutf8(s::PyObject) = pycall(s."encode", PyObject, "utf-8")
pyutf8(s::String) = pyutf8(PyObject(s))

@testset "PyBuffer" begin
    @testset "String Buffers" begin
        b = PyCall.PyBuffer(pyutf8("test string"))
        @test ndims(b) == 1
        @test (length(b),) == (length("test string"),) == (size(b, 1),) == size(b)
        @test stride(b, 1) == 1
        @test PyCall.iscontiguous(b) == true
    end

    if !npy_initialized
        println(stderr, "Warning: skipping array related buffer tests since NumPy not available")
    else
        np = pyimport("numpy")
        listpy = pybuiltin("list")
        arrpyo(args...; kwargs...) =
            pycall(np."array", PyObject, args...; kwargs...)
        listpyo(args...) = pycall(listpy, PyObject, args...)
        pytestarray(sz::Int...; order="C") =
            pycall(arrpyo(1.0:prod(sz), "d")."reshape", PyObject, sz, order=order)

        @testset "Non-native-endian" begin
            wrong_endian_str = ENDIAN_BOM == 0x01020304 ? "<" : ">"
            wrong_endian_arr =
                pycall(np."ndarray", PyObject, 2; buffer=UInt8[0,1,3,2],
                                                   dtype=wrong_endian_str*"i2")
            # Not supported, so throws
            @test_throws ArgumentError NoCopyArray(wrong_endian_arr)
            @test_throws ArgumentError PyArray(wrong_endian_arr)
        end

        @testset "dtype should match eltype" begin
            npy2jl = Dict("int32"=>Int32, "float32"=>Float32,
                          "int64"=>Int64, "float64"=>Float64)
            for (pytype, jltype) in npy2jl
                @testset "dtype $pytype should match eltype $jltype" begin
                    jlarr = jltype[1:10;]
                    nparr = arrpyo(jlarr, dtype=pytype)
                    @test pystr(nparr."dtype") == pytype
                    jlarr2 = convert(PyAny, nparr)
                    @test eltype(jlarr2) == jltype
                    @test jlarr2 == jlarr
                end
            end
        end

        # f_contiguous(T, sz, st)
        @testset "f_contiguous 1D" begin
            # contiguous case: stride == sizeof(T)
            @test f_contiguous(Float64, (4,), (8,)) == true
            # non-contiguous case: stride != sizeof(T)
            @test f_contiguous(Float64, (4,), (16,)) == false
        end

        @testset "f_contiguous 2D" begin
            # contiguous: st[1] == sizeof(T), st[2] == st[1]*sz[1]
            @test f_contiguous(Float64, (4, 2), (8, 32)) == true
            # non-contiguous: stride != sizeof(T), but st[2] == st[1]*sz[1]
            @test f_contiguous(Float64, (4, 2), (16, 64)) == false
            # non-contiguous: stride == sizeof(T), but st[2] != st[1]*sz[1]
            @test f_contiguous(Float64, (4, 2), (8, 64)) == false
        end

        @testset "copy f_contig 1d" begin
            apyo = arrpyo(1.0:10.0, "d")
            pyarr = PyArray(apyo)
            jlcopy = copy(pyarr)
            @test pyarr.f_contig == true
            @test pyarr.c_contig == true
            @test all(jlcopy .== pyarr)
            # check it's not aliasing the same data
            jlcopy[1] = -1.0
            @test pyarr[1] == 1.0
        end

        @testset "copy c_contig 2d" begin
            apyo = pytestarray(2,3) # arrpyo([[1,2,3],[4,5,6]], "d")
            pyarr = PyArray(apyo)
            jlcopy = copy(pyarr)
            @test pyarr.c_contig == true
            @test pyarr.f_contig == false
            # check all is in order
            for i in 1:size(pyarr, 1)
                for j in 1:size(pyarr, 1)
                    @test jlcopy[i,j] == pyarr[i,j]
                end
            end
            # check it's not aliasing the same data
            jlcopy[1,1] = -1.0
            @test pyarr[1,1] == 1.0
        end

        @testset "Non contiguous PyArrays" begin
            @testset "1d non-contiguous" begin
                # create an array of four Int32s, with stride 8
                nparr = pycall(np."ndarray", PyObject, 4,
                                buffer=UInt32[1,0,1,0,1,0,1,0],
                                dtype="i4", strides=(8,))
                pyarr = PyArray(nparr)

                # The convert goes via a PyArray then a `copy`
                @test convert(PyAny, nparr) == [1, 1, 1, 1]

                @test eltype(pyarr) == Int32
                @test sizeof(eltype(pyarr)) == 4
                @test pyarr.info.st == (8,)
                # not f_contig because not contiguous
                @test pyarr.f_contig == false
                @test copy(pyarr) == Int32[1, 1, 1, 1]
            end

            @testset "2d non-contiguous" begin
                nparr = pycall(np."ndarray", PyObject,
                                buffer=UInt32[1,0,2,0,1,0,2,0,
                                              1,0,2,0,1,0,2,0], order="f",
                                dtype="i4", shape=(2, 4), strides=(8,16))
                pyarr = PyArray(nparr)

                # The convert goes via a PyArray then a `copy`
                @test convert(PyAny, nparr) == [1 1 1 1; 2 2 2 2]
                pyarr = convert(PyArray, nparr)
                @test eltype(pyarr) == Int32
                @test pyarr.info.st == (8, 16)
                # not f_contig because not contiguous
                @test pyarr.f_contig == false
                @test copy(pyarr) == Int32[1 1 1 1; 2 2 2 2]
            end
        end

        @testset "NoCopyArray 1d" begin
            ao = arrpyo(1.0:10.0, "d")
            pybuf = PyBuffer(ao, PyBUF_ND_CONTIGUOUS)
            T, native_byteorder = array_format(pybuf)
            @test T == Float64
            @test native_byteorder == true
            @test size(pybuf) == (10,)
            @test strides(pybuf) == (1,) .* sizeof(T)
            nca = NoCopyArray(ao)
            @test !(nca isa PermutedDimsArray)
            @test nca isa Array
            @test nca[3] == get(ao,2)
            @test nca[4] == get(ao,3)
        end

        @testset "NoCopyArray 2d f-contig" begin
            ao = arrpyo(reshape(1.0:12.0, (3,4)) |> collect, "d", order="F")
            pybuf = PyBuffer(ao, PyBUF_ND_CONTIGUOUS)
            T, native_byteorder = array_format(pybuf)
            @test T == Float64
            @test native_byteorder == true
            @test size(pybuf) == (3,4)
            @test strides(pybuf) == (1, 3) .* sizeof(T)
            nca = NoCopyArray(ao)
            @test !(nca isa PermutedDimsArray)
            # @show typeof(nca) (nca isa Array)
            @test nca isa Array
            @test size(nca) == (3,4)
            @test strides(nca) == (1,3)
            @test nca[3,2] == get(ao, (2,1))
            @test nca[2,3] == get(ao, (1,2))
        end

        @testset "NoCopyArray 3d c-contig" begin
            ao = pytestarray(3,4,5)
            pybuf = PyBuffer(ao, PyBUF_ND_CONTIGUOUS)
            T, native_byteorder = array_format(pybuf)
            @test T == Float64
            @test native_byteorder == true
            @test size(pybuf) == (3,4,5)
            @test strides(pybuf) == (20,5,1) .* sizeof(T)
            nca = NoCopyArray(ao)
            @test nca isa PermutedDimsArray
            @test !(nca isa Array)
            @test size(nca) == (3,4,5)
            @test strides(nca) == (20,5,1)
            @test nca[2,3,4] == get(ao, (1,2,3))
            @test nca[3,2,4] == get(ao, (2,1,3))
        end

        @testset "isbuftype" begin
            @test isbuftype(PyObject(0)) == false
            @test isbuftype(listpyo((1.0:10.0...,))) == false
            @test isbuftype(arrpyo(1.0:10.0, "d")) == true
            @test isbuftype(PyObject([1:10...])) == true
        end

        # TODO maybe move these to a test_pyarray.jl
        @testset "setdata!" begin
            ao1 = arrpyo(1.0:10.0, "d")
            pyarr = convert(PyArray, ao1)
            ao2 = arrpyo(11.0:20.0, "d")
            setdata!(pyarr, ao2)
            @test all(pyarr[1:10] .== 11.0:20.0)
        end

        @testset "similar on PyArray PyVec getindex" begin
            jlarr1 = [1:10;]
            jlarr2 = hcat([1:10;], [1:10;])
            pyarr1 = pycall(np."array", PyArray, jlarr1)
            pyarr2 = pycall(np."array", PyArray, jlarr2)
            @test all(pyarr1[1:10]    .== jlarr1[1:10])
            @test all(pyarr2[1:10, 2] .== jlarr2[1:10, 2])
            @test all(pyarr2[1:10, 1:2] .== jlarr2)
        end
    end
end


text:
using ConnectSDK
using JSON
using HTTPClient.HTTPC
using Base.Test

# write your own tests here
connect = ConnectClient("55b078611f2ddd01a84c28db", "15765BE8FF529140640431906EE1D749-79F100C25D90EB54EACE1376058FD10D32ABEEBE0C8353157C7708FD623EC5093BA915C70F4FE2B878C2DEAD3743E60A")
temp = JSON.parse(bytestring(HTTPC.get("api.openweathermap.org/data/2.5/weather?id=7839562&units=metric&APPID=159d95e992dd299342cd8190e6e1f964").body))["main"]["temp"]
testdata = {"client" => "Julia", "temp" => temp}

r = push(connect, "tests", testdata)
@test r.http_code == 200


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ e6e24892-e2b7-472f-9a5c-7d1cbe34eeb6
html"<button onclick=present()> Presente </button>"

# ╔═╡ f92f9d30-b368-11eb-2e66-a7bbf5a4a943
md"""
# Plano da disciplina de Mecânica dos Sólidos II

## Ementa da Disciplina:

1. Diagramas de Força Axial, Cortante e de Momentos.
2. Tensão e Deformação.
3. Torção.
4. Tensão de Flexão e de Cisalhamento em Vigas.
5. Tensões Compostas.
6. Transformação de Tensão.
7. Critério de Escoamento e de Fratura.
8. Vasos de Pressão.
9. Deflexão de Vigas.
10. Métodos de Energia.
11. Flambagem de Colunas.
12. Método dos Elementos Finitos.

## Resultados de aprendizagem esperados

#### 🥇 Identificar as teorias e equações que fundamentam a Mecânica dos Sólidos

#### 🥈 Identificar e calcular os esforços aplicados que se relacionam com a tensão e deformação nos materias.

## Metodologia

### Estratégia

Com base no livro de ELMÔR FILHO, Gabriel et al. $(html"Uma nova sala de aula é possível: aprendizagem ativa na educação em engenharia.") Rio de Janeiro: LTC, 2019.

Será aplicada a **Sala de aula invertida** com utilização da Estratégia *Co-op co-op* que é uma estratégia cooperativa de aprendizagem ativa.

🧑🎓 Os estudantes trabalham juntos, em pequenas equipes (4 ou 5 membros), para investigar um tópico ou resolver um problema, e produzir um produto de grupo que eles compartilham com toda a classe.

🧑🎓 Os estudantes cooperam dentro de suas equipes para produzir algo que irá beneficiar o grande grupo (Estudos de caso que ficarão disponíveis na internet até o fim da disciplina).

👨🔬 Mediação do professor em todas as etapas de desenvolvimento da estratégia.

## Roteiro

#### Momento Pré

🕐 1. O professor anuncia por meio de uma orientação por escrito (Sala):

   - A atividade a ser desenvolvida.
   - O respectivo tema e subtópicos.
   - Uma síntese da estratégia e
   - O cronograma incluindo o período para as avaliações da atividade.

🕑 2. Os estudantes têm a tarefa de ler sobre o assunto e postar no Classroom resumo de uma página do que estudou (Resumo pré) (Extra sala).

## Roteiro

#### Momento Aula

🕒 3. Rápida verificação do entendimento do tema pelo professor (Sala).

🕓 4. Promoção de discussão sobre os interesses no tema que será estudado. O que querem aprender e experienciar durante a unidade? (Sala).

🕔 5. Seleção do tópico que cada grupo irá estudar. Os estudantes discutem seus interesses e selecionam um tópico para seu grupo (Aula).

🕕 6. Cada grupo divide seu tópico em subtópicos para dividir o trabalho dentro do grupo (Aula).

🕖 7. Preparação dos subtópicos pelos membros do grupo (estudo individual).

🕗 8. Os membros do grupo discutem o tema, como um painel de especialistas, apresentando os conceitos e questões chaves (reunião do grupo).

🕘 9. Apresentação formais dos subtópicos em cada grupo (reunião do grupo).

🕙 10. Apresentações de subtópicos pelos grupos (seleção por sorteio a critério do professsor) (sala).

## Roteiro

#### Momento Pós

🕚 11. Preparação das apresentações de cada um dos grupos (reunião do grupo).

🕛 12. Envio via Classroom do relatório da atividade para postagem na internet.

🕐 13. Apresentações de cada um dos grupos (sala).

🕑 14. Avaliações (Extra sala).

🕒 15. Feedback (sala).

## Recursos e Materiais Didáticos

1. Google Classroom;
2. Google Meet para encontros;
3. Jamboard como quadro branco;
4. Google Forms para formulários e avaliações;
5. Google Drive para arquivos (Docs, Slides, Planilhas, Desenhos);
6. Os recurso $(html"Vscode") e $(html"Linguagem Julia") e $(html"Pluto.jl") para Cálculos.
7. Página do professor para registro dos casos.

## Indicações de Referências

### Bibliografia principal:

1. BEER, Ferdinand et al; Mecânica dos Materiais, 5ª ed. 2011.

### Bibliografia Complementar:

2. BEER, Ferdinand Pierre; et al; Resistência dos materiais.

3. HIBBELER, R. C. Resistência dos materiais.

4. POPOV, Egor P. Introdução à mecânica dos sólidos.

5. NORTON, Robert L. Projeto de máquinas: uma abordagem integrada.

6. BUDYNAS, Richard G ; NISBETT, J. Keith. Elementos de máquinas de Shigley.

## Procedimento de Avaliação da Apredizagem

Avaliação será:

- Contínua
- Individual
- Por atividade a ser desenvolvida

## Procedimento de Avaliação da Apredizagem

### Regras:

- As avaliações serão feitas via $(html"Formulário Eletrônico no Google Forms (Acesse aqui)"), com excessão do Resumo pré que deve ser entregue via Classroom.
- Só serão aceitos emails da Unifesspa.
- No caso de duplicidade de avaliação será considerada sempre a última avaliação enviada.
- Serão consideradas apenas as avaliações realizadas no período estipulado, ou até quando estiverem disponíveis no formulário.
- Cada avaliação receberá uma pontuação de 0 a 100.
- A nota para o conceito final será a média das notas das atividades propostas.

## Procedimento de Avaliação da Apredizage

### Composição das Notas

A seguir são enumeradas os itens que serão avaliados.

Para cada atividade serão avaliados:

1. Apresentações dos grupos (peso 1).
2. Contribuições individuais (peso 1).
3. Relatório em arquivo html para postagem (peso 1).

## Procedimento de Avaliação da Apredizagem

### Composição das Notas

#### 1. Apresentações dos grupos (peso 1)

Avalições que serão consideradas neste item:

- Avaliação por pares (peso 1)
- Avaliação do professor (peso 1)

## Procedimento de Avaliação da Apredizagem

### Composição das Notas

#### 2. Contribuições individuais (peso 1)

Avaliações que serão consideradas:

- Avaliação do resumo pré (peso 1)
- Avaliação por pares (peso 1)
- Avaliação do professor (peso 1)

## Procedimento de Avaliação da Apredizagem

### Composição das Notas

#### 3. Relatório em arquivo html para postagem (peso 1).

- Autoavaliação (peso 1)
- Avaliação por pares (peso 1)
- Avaliação do professor (peso 1)

"""

# ╔═╡ Cell order:
# ╟─e6e24892-e2b7-472f-9a5c-7d1cbe34eeb6
# ╟─f92f9d30-b368-11eb-2e66-a7bbf5a4a943

max_stars_repo_path: Plano da disciplina.jl

max_stars_repo_name: jefferds/Solid-mechanics-II

max_issues_repo_path: Plano da disciplina.jl

max_issues_repo_name: jefferds/Solid-mechanics-II

max_forks_repo_path: Plano da disciplina.jl

max_forks_repo_name: jefferds/Solid-mechanics-II


text:
using DelimitedFiles
using LinearAlgebra

function convolve_linear(signal::Array{T, 2}, filter::Array{T, 2},
                         output_size) where {T <: Number}

    # convolutional output
    out = Array{Float64,2}(undef, output_size)
    sum = 0

    for i = 1:output_size[1]
        for j = 1:output_size[2]
            for k = max(1, i-size(filter)[1]):i
                for l = max(1, j-size(filter)[2]):j
                    if k <= size(signal)[1] && i-k+1 <= size(filter)[1] &&
                       l <= size(signal)[2] && j-l+1 <= size(filter)[2]
                        sum += signal[k,l] * filter[i-k+1, j-l+1]
                    end
                end
            end

            out[i,j] = sum
            sum = 0
        end
    end

    return out
end

function create_gaussian_kernel(kernel_size)

    kernel = zeros(kernel_size, kernel_size)

    # The center must be offset by 0.5 to find the correct index
    center = kernel_size * 0.5 + 0.5

    sigma = sqrt(0.1*kernel_size)

    for i = 1:kernel_size
        for j = 1:kernel_size
            kernel[i,j] = exp(-((i-center)^2 + (j-center)^2) / (2*sigma^2))
        end
    end

    return normalize(kernel)

end

function create_sobel_operators()
    Sx = [1.0, 2.0, 1.0]*[-1.0 0.0 1.0] / 9
    Sy = [-1.0, 0.0, 1.0]*[1.0 2.0 1.0] / 9

    return Sx, Sy
end

function compute_sobel(signal)
    Sx, Sy = create_sobel_operators()

    Gx = convolve_linear(signal, Sx, size(signal) .+ size(Sx))
    Gy = convolve_linear(signal, Sy, size(signal) .+ size(Sy))

    return sqrt.(Gx.^2 .+ Gy.^2)
end

# Simple function to create a square grid with a circle embedded inside of it
function create_circle(image_resolution, grid_extents, radius)
    out = zeros(image_resolution, image_resolution)

    for i = 1:image_resolution
        x_position = ((i-1)*grid_extents/image_resolution)-0.5*grid_extents
        for j = 1:image_resolution
            y_position = ((j-1)*grid_extents/image_resolution)-0.5*grid_extents
            if x_position^2 + y_position^2 <= radius^2
                out[i,j] = 1.0
            end
        end
    end

    return out
end

function main()

    # Random distribution in x
    x = rand(100, 100)

    # Gaussian signals
    y = [exp(-(((i-50)/100)^2 + ((j-50)/100)^2)/.01) for i = 1:100, j=1:100]

    # Normalization is not strictly necessary, but good practice
    normalize!(x)
    normalize!(y)

    # full convolution, output will be the size of x + y
    full_linear_output = convolve_linear(x, y, size(x) .+ size(y))

    # simple boundaries
    simple_linear_output = convolve_linear(x, y, size(x))

    # outputting convolutions to different files for plotting in external code
    writedlm("full_linear.dat", full_linear_output)
    writedlm("simple_linear.dat", simple_linear_output)

    # creating simple circle and 2 different Gaussian kernels
    circle = create_circle(50,2,0.5)

    normalize!(circle)

    small_kernel = create_gaussian_kernel(3)
    large_kernel = create_gaussian_kernel(25)

    small_kernel_output = convolve_linear(circle, small_kernel,
                                          size(circle).+size(small_kernel))
    large_kernel_output = convolve_linear(circle, large_kernel,
                                          size(circle).+size(large_kernel))

    writedlm("small_kernel.dat", small_kernel_output)
    writedlm("large_kernel.dat", large_kernel_output)

    # Using the circle for Sobel operations as well
    sobel_output = compute_sobel(circle)

    writedlm("sobel_output.dat", sobel_output)

end


text:
# Dynamic Time Warping (DTW)

using VoiceConversion.DTWs

println("testing: Dynamic time warping")

let
    # template sequence [1,2,3,4]
    v1 = [1. 2 3;1 2 4;1 8 5;10 3 6]'
    # sequence to be aligned [1 2 2 3 4]
    v2 = [1. 2 3;1 2 4;1 2 5; 1 8 5;10 3 6]'

    # Perform alignment
    d = DTW(bstep=1, fstep=0)
    indices = fit!(d, v1, v2)

    expected = Int[1, 2, 2, 3, 4]
    @test indices == expected
end

let
    v1 = [0.; 1; 2; 3; 4; 5]'
    v2 = [0.; 0; 1; 2; 3; 4; 4; 5;]'

    # Perform alignment
    d = DTW(bstep=1, fstep=0)
    indices = fit!(d, v1, v2)

    expected = Int[1, 1, 2, 3, 4, 5, 5, 6]
    @test indices == expected
end


text:
#---------------------------------------
# types for a single EB sample
#---------------------------------------
abstract type EBayesSample{T<:Number} end
abstract type AbstractNormalSample{T<:Number} <: EBayesSample{T} end

"""
    NormalSample(Z,σ)

A observed sample ``Z`` drawn from a Normal distribution with known variance ``\\sigma^2 > 0``.

```math
Z \\sim \\mathcal{N}(\\mu, \\sigma^2)
```

``\\mu`` is assumed unknown. The type above is used when the sample ``Z`` is to be used for estimation or inference of ``\\mu``.

```julia
NormalSample(0.5, 1.0)          #Z=0.5, σ=1
```
"""
struct NormalSample{T <: Number} <: AbstractNormalSample{T}
    Z::T
    σ::T
end

function NormalSample(Z::T) where {T<:Number}
    NormalSample(Z, one(T))
end

NormalSample(Z::Number, σ::Number) = NormalSample(promote(Z, σ)...)

response(s::NormalSample) = s.Z
Statistics.var(s::NormalSample) = s.σ^2

eltype(s::NormalSample{T}) where T = T
zero(s::NormalSample{T}) where T = zero(T)
support(ss::NormalSample) = RealInterval(-Inf, +Inf)

#---------------------------------------
# types for a single EB sample
#---------------------------------------
const EBayesSamples = AbstractArray{EBS} where EBS <: EBayesSample

const NormalSamples{T} = StructArray{NormalSample{T}} where T

function NormalSamples(Zs::AbstractVector{T}, σs::AbstractVector{T}) where T
    NormalSamples{T}((Zs, σs))
end

response(ss::NormalSamples) = ss.Z
Statistics.var(ss::NormalSamples) = ss.σ .^ 2

"""
    StandardNormalSample(Z)

A observed sample ``Z`` drawn from a Normal distribution with known variance ``\\sigma^2 =1``.

```math
Z \\sim \\mathcal{N}(\\mu, 1)
```

``\\mu`` is assumed unknown. The type above is used when the sample ``Z`` is to be used for estimation or inference of ``\\mu``.

```julia
StandardNormalSample(0.5)          #Z=0.5
```
"""
struct StandardNormalSample{T <: Number} <: AbstractNormalSample{T}
    Z::T
end

response(s::StandardNormalSample) = s.Z
Statistics.var(s::StandardNormalSample) = 1

eltype(s::StandardNormalSample{T}) where T = T
zero(s::StandardNormalSample{T}) where T = zero(T)
support(ss::StandardNormalSample) = RealInterval(-Inf, +Inf)

const AbstractNormalSamples = AbstractArray{NS} where NS <: Union{NormalSample,
                                                                  StandardNormalSample}

response(ss::AbstractNormalSamples) = response.(ss)
Statistics.var(ss::AbstractNormalSamples) = var.(ss)
zeros(ss::AbstractNormalSamples) = zeros(eltype(response(ss)), length(ss))

# Poisson

# Binomial

# Replicated , ReplicatedArray


text:
export register_hashpipe_thread, find_hashpipe_thread

# Hashpipe thread abstract type
"""
    HashpipeThread

Abstract type describing a Hashpipe thread.

    name - A string containing the thread's name. Used to match command line
            thread spcifiers to thread metadata so that the pipeline can be
            constructed as specified on the command line.
    skey - A string containing the thread's status buffer "status" key.
            Typically 8 characters or less, uppercase, and ends with "STAT".
            If it is non-NULL and non-empty, HASHPIPE will automatically
            store/update this key in the status buffer with the thread's
            status at initialization ("init") and exit ("exit").
    init - A pointer to the thread's initialization function.
    run  - A pointer to the thread's run function.
    ibuf - A structure describing the thread's input data buffer (if any).
    obuf - A structure describing the thread's output data buffer (if any).

## Description of Hashpipe threads
A hashpipe_thread structure encapsulates metadata and functionality for one
or more threads that can be used in a processing pipeline.  The hashpipe
executable dynamically assembles a pipeline at runtime consisting of
multiple hashpipe threads.

Hashpipe threads must register themselves with the hashpipe executable via a
call to register_hashpipe_thread().  This is typically performed from a
static C function with the constructor attribute in the hashpipe thread's
source file.

Hashpipe threads are identified by their name.  The hashpipe executable
finds (registered) hashpipe threads by their name.  A hashpipe thread can be
input-only, output-only, or both input and output.  An input thread has an
associated output data buffer into which it writes data.  An output thread
has an associated input data buffer from which it reads data.  An
input/output thread has both.

Input-only threads source data into the pipeline.  They do not get their
input data from a shared memory ring buffer.  They get their data from
external sources (e.g.  files or the network) or generate it internally
(e.g.  for test vectors).  Input-only threads have an output data buffer,
but no input data buffer (their input does not come from a shared memory
ring buffer).

Output-only threads sink data from the pipeline.  Thy do not put their
output data into a shared memory ring buffer.  They send their data to
external sinks (e.g. files or the network) of consume it internally (e.g.
comparing against expected output).  Output-only threads have an input data
buffer, but no output data buffer (their output data does not go the a
shared memory ring buffer).

Input/output threads get their input data from one shared memory region
(their input data buffer), process it, and store the output data in another
shared memory region (their output data buffer).
"""
abstract type HashpipeThread end

function init(thread::HashpipeThread)
    @error "Using default Hashpipe thread init function. Define a specific init function."
end

function run(thread::HashpipeThread)
    @error "Using default Hashpipe thread run function. Define a specific run function."
end

"""
    THREAD_OK

OK status from run.
"""
const global THREAD_OK  = 0
"""
    THREAD_ERROR

Error status from run.
"""
const global THREAD_ERROR = -1
"""
    MAX_HASHPIPE_THREADS

Maximum number of threads that be defined by plugins.
"""
const global MAX_HASHPIPE_THREADS = 1024

"""
    thread_desc_t

Metadata describing a Hashpipe thread (usually defined by a plugin per thread).

The thread_desc structure is used to store metadata describing a hashpipe thread.
Typically a hashpipe plugin will define one of these hashpipe thread
descriptors per hashpipe thread.
"""
struct thread_desc_t
    name::Cstring # TODO: Double check on NULL terminated assumption with Dave
    skey::Cstring # ^^
    init::Ptr{Cvoid}
    run::Ptr{Cvoid}
    ibuf_desc::databuf_desc_t
    obuf_desc::databuf_desc_t
end

"""
    thread_args_t

Implementation data for a Hashpipe thread.
"""
mutable struct thread_args_t
    thread_desc::Ptr{thread_desc_t}
    instance_id::Cint
    input_buffer::Cint
    output_buffer::Cint
    cpu_mask::UInt32
    finished::Cint
    finished_c::Ptr{Cvoid} # TODO: Change to mimic pthread_cond_t
    finished_m::Ptr{Cvoid} # TODO: Change to mimic pthread_mutex_t
    st::status_t
    ibuf::Ptr{databuf_t}
    obuf::Ptr{databuf_t}
    user_data::Ptr{Cvoid}
end

"""
    register_hashpipe_thread(ptm)

Register a Hashpipe thread with the pipeline executable.
"""
function register_hashpipe_thread(ptm)
    ccall((:register_hashpipe_thread, libhashpipe), Cint, (Ptr{thread_desc_t},), ptm)
end

"""
    find_hashpipe_thread(name)

Find a Hashpipe thread with the given name.

Generally used only by the hashpipe executable.  Returns a pointer to its
hashpipe_thread_desc_t structure or NULL if a test with the given name is not found.

Names are case-sensitive.
"""
function find_hashpipe_thread(name)
    ccall((:find_hashpipe_thread, libhashpipe), Ptr{thread_desc_t}, (Cstring,), name)
end

"""
    run_threads()

Function that threads use to determine whether to keep running.
"""
function run_threads()
    ccall((:run_threads, libhashpipe), Cint, ())
end

"""
    get_cpu_affinity()
"""
function get_cpu_affinity()
    ccall((:get_cpu_affinity, libhashpipe), UInt32, ())
end


text:
# B-Spline Basis function for Bezier Bernsetirn splines
@doc """
### basisFunctions

Computes the basis functions needed to compute B-Splines. The calculation starts
from scratch and the entire triangular array of B-spline values of orders
1,2,...,jhigh is generated order by order (column by column).

**Inputs**

*  `U`     : Knot vector
*  `order` : Order of B-spline
*  `u`     : Location where basis functions are to be evaluated
*  `span`  : Knot span index
*  `N`     : Basis Functions

**Outputs**

*  None

SOURCE: "A practical guide to splines" bt C. de Boor Page 112 BSPLVB

"""->
function basisFunctions(U::AbstractVector{T}, order, u::T2, span, N::AbstractVector{T}) where {T, T2}

  # (U, order, u, span, N)
  #order = map.order[di]
  # U = map.edge_knot[di]

  Tffd = promote_type(T, T2)
  jmax = 20 # arbitrarily imposed
  delta_l = zeros(Tffd, jmax)
  delta_r = zeros(Tffd, jmax)

  N[1] = 1.0

  for j = 1:(order-1)
    delta_r[j] = U[span+j] - u
    delta_l[j] = u - U[span+1-j]
    saved = 0.0
    for i = 1:j
      term = N[i]/(delta_r[i] + delta_l[j+1-i])
      N[i] = saved + delta_r[i]*term
      saved = delta_l[j+1-i]*term
    end  # end for i = 1:r
    N[j+1] = saved
  end  # end for j = 1:(order-1)

  return nothing
end

@doc """
### derivBasisFunctions

Computes the basis functions and its first & second (method 2) derivative at a
particular point u.

**Inputs**

*  `U`     : Knot vector
*  `order` : Order of B-spline
*  `u`     : Location where basis functions are to be evaluated
*  `span`  : Knot span index
*  `N`     : Basis Functions
*  `Nderiv`: Derivative of the basis function

**Outputs**

*  None

"""->

function derivBasisFunctions(u::T, U::AbstractVector{T2},
                             order, span, N, Nderiv) where {T, T2}

  Tffd = promote_type(T, T2)

  # order = map.order[di]
  N[1] = 1.0
  Nderiv[1] = 0.0

  # dr =sview(map.dr, :, :)
  # dl =sview(map.dl, :, :)

  dr = zeros(Tffd, order-1)
  dl = zeros(Tffd, order-1)

  if order > 1
    for k = 1:order-1
      kp1 = k + 1
      dr[k] = U[span+k] - u
      dl[k] = u - U[span+1-k]
      saved = 0.0
      dsaved = 0.0
      for i = 1:k
        temp = 1/( dr[i] + dl[k+1-i] )
        dtemp = Nderiv[i]*temp
        temp = N[i]*temp

        N[i] = saved + dr[i]*temp
        Nderiv[i] = dsaved - temp + dr[i]*dtemp

        saved = dl[k+1-i]*temp
        dsaved = temp + dl[k+1-i]*dtemp
      end  # end for i = 1:k
      N[k+1] = saved
      Nderiv[k+1] = dsaved
    end  # end for k = 1:order-1
  end  # end if order > 1

  return nothing
end

function derivBasisFunctions(map, N, Nderiv, N2deriv, di, u, span)

  order = map.order[di]
  N[1] = 1.0
  Nderiv[1] = 0.0
  N2deriv[1] = 0.0

  dr =sview(map.dr, :, :)
  dl =sview(map.dl, :, :)

  if order > 1
    for k = 1:order-1
      dr[k,di] = map.knot[span+k,di] - u
      dl[k,di] = u - map.knot[span+1-k, di]
      saved = 0.0
      dsaved = 0.0
      d2saved = 0.0
      for i = 1:k
        temp = 1/( dr[i,di] + dl[k+1-i, di] )
        dtemp = Nderiv[i]*temp
        d2temp = N2deriv[i]*temp
        temp = bval[i]*temp

        N[i] = saved + dr[i,di]*temp
        Nderiv[i] = dsaved - temp + dr[i,di]*dtemp
        N2deriv[i] = d2saved - 2*dtemp + dr[i,di]*d2temp

        saved = dl[k+1-i,di]*temp
        dsaved = temp + dl[k+1-i,di]*dtemp
        d2saved = 2*dtemp + dl[k+1-i,di]*d2temp
      end  # End for i = 1:k
      N[k+1] = saved
      Nderiv[k+1] = dsaved
      N2deriv[k+1] = d2saved
    end # End for k = 1:order-1
  end  # End if order > 1

  return nothing
end


text:
include("smearing_functions.jl")
abstract type Smear end

function get_smear_type(potential::Symbol)
    if potential == :ttm3
        return Smear_TTM3()
    elseif potential == :ttm21
        return Smear_TTM21()
    elseif potential == :ttm22
        return Smear_TTM22()
    elseif potential == :ttm4
        return Smear_TTM4()
    else
        @assert false "Didn't receive a valid potential. Pass as a symbol :ttm3, :ttm21, :ttm22, or :ttm4."
    end
end

struct Smear_TTM3 <: Smear
    aCC::Float64
    aCD::Float64
    aDD_intermolecular::Float64
    aDD_intramolecular::Float64
    smear1::Function
    smear2::Function
    smear3::Function
    Smear_TTM3() = new(0.175, 0.175, 0.175, 0.175, smear1_ttm3, smear2_ttm3, smear3_ttm3)
end

struct Smear_TTM21 <: Smear
    aCC::Float64
    aCD::Float64
    aDD_intermolecular::Float64
    aDD_intramolecular::Float64
    smear1::Function
    smear2::Function
    smear3::Function
    Smear_TTM21() = new(0.2, 0.2, 0.3, 0.3, smear1_ttm3, smear2_ttm3, smear3_ttm3)
end

struct Smear_TTM22 <: Smear
    aCC::Float64
    aCD::Float64
    aDD_intermolecular::Float64
    aDD_intramolecular::Float64
    smear1::Function
    smear2::Function
    smear3::Function
    Smear_TTM22() = new(0.2, 0.2, 0.3, 0.3, smear1_ttm3, smear2_ttm3, smear3_ttm3)
end

struct Smear_TTM4 <: Smear
    aCC::Float64
    aCD::Float64
    aDD_intermolecular::Float64
    aDD_intramolecular::Float64
    smear1::Function
    smear2::Function
    smear3::Function
    Smear_TTM4() = new(0.4, 0.4, 0.055, 0.626, smear1_ttm4, smear2_ttm4, smear3_ttm4)
end

max_stars_repo_name: heindelj/MolecularPotentials

max_issues_repo_name: heindelj/MolecularPotentials

max_forks_repo_name: heindelj/MolecularPotentials


text:
####################################
# The AbstractBlockArray interface #
####################################

"""
    abstract AbstractBlockArray{T, N} <: AbstractArray{T, N}

The abstract type that represents a blocked array. Types that implement
the `AbstractBlockArray` interface should subtype from this type.

** Typealiases **

* `AbstractBlockMatrix{T}` -> `AbstractBlockArray{T, 2}`

* `AbstractBlockVector{T}` -> `AbstractBlockArray{T, 1}`

* `AbstractBlockVecOrMat{T}` -> `Union{AbstractBlockMatrix{T}, AbstractBlockVector{T}}`
"""
abstract type AbstractBlockArray{T, N} <: LayoutArray{T, N} end
const AbstractBlockMatrix{T} = AbstractBlockArray{T, 2}
const AbstractBlockVector{T} = AbstractBlockArray{T, 1}
const AbstractBlockVecOrMat{T} = Union{AbstractBlockMatrix{T}, AbstractBlockVector{T}}

block2string(b, s) = string(join(map(string,b), '×'), "-blocked ", Base.dims2string(s))
_block_summary(a) = string(block2string(blocksize(a), size(a)), " ", typeof(a))
Base.summary(a::AbstractBlockArray) = _block_summary(a)
_show_typeof(io, a) = show(io, typeof(a))
function _block_summary(io, a)
    print(io, block2string(blocksize(a), size(a)))
    print(io, ' ')
    _show_typeof(io, a)
end
Base.summary(io::IO, a::AbstractBlockArray) = _block_summary(io, a)

# avoid to_shape which complicates axes
Base.similar(a::AbstractBlockArray{T}) where {T}                             = similar(a, T)
Base.similar(a::AbstractBlockArray, ::Type{T}) where {T}                     = similar(a, T, axes(a))
Base.similar(a::AbstractBlockArray{T}, dims::Tuple) where {T}                = similar(a, T, dims)

# If all we know is size, just return an Array which conforms to BlockArray interface
Base.similar(::Type{<:AbstractBlockArray{T,N}}, dims::Dims) where {T,N} = similar(Array{T,N}, dims)

Base.IndexStyle(::Type{<:AbstractBlockArray}) = IndexCartesian()

# need to overload axes to return BlockAxis
@inline size(block_array::AbstractBlockArray) = map(length, axes(block_array))
@inline axes(block_array::AbstractBlockArray) = throw(error("axes for ", typeof(block_array), " is not implemented"))

"""
    BlockBoundsError([A], [inds...])

Thrown when a block indexing operation into a block array, `A`, tried to access an out-of-bounds block, `inds`.
"""
struct BlockBoundsError <: Exception
    a::Any
    i::Any
    BlockBoundsError() = new()
    BlockBoundsError(a::AbstractArray) = new(a)
    BlockBoundsError(a::AbstractArray, @nospecialize(i)) = new(a,i)
end

BlockBoundsError(a::AbstractArray, I::Block) = BlockBoundsError(a, I.n)

function Base.showerror(io::IO, ex::BlockBoundsError)
    print(io, "BlockBoundsError")
    if isdefined(ex, :a)
        print(io, ": attempt to access ")
            print(io, summary(ex.a))
        if isdefined(ex, :i)
            print(io, " at block index [")
            join(io, ex.i, ',')
            print(io, ']')
        end
    end
end

"""
    blockcheckbounds(A, inds...)

Throw a `BlockBoundsError` if the specified block indexes are not in bounds for the given block array.
Subtypes of `AbstractBlockArray` should
specialize this method if they need to provide custom block bounds checking behaviors.

```jldoctest; setup = quote using BlockArrays end
julia> A = BlockArray(rand(2,3), [1,1], [2,1]);

julia> blockcheckbounds(A, 3, 2)
ERROR: BlockBoundsError: attempt to access 2×2-blocked 2×3 BlockMatrix{Float64, Matrix{Matrix{Float64}}, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}} at block index [3,2]
[...]
```
"""
@inline function blockcheckbounds(A::AbstractArray{T, N}, i::Vararg{Integer, N}) where {T,N}
    if blockcheckbounds(Bool, A, i...)
        return
    else
        throw(BlockBoundsError(A, i))
    end
end

@inline function blockcheckbounds(::Type{Bool}, A::AbstractArray{T, N}, i::Vararg{Integer, N}) where {T,N}
    n = blockaxes(A)
    k = 0
    for idx in 1:N # using enumerate here will allocate
        k += 1
        @inbounds _i = i[idx]
        Block(_i) in n[k] || return false
    end
    for idx = N+1:length(i)
        isone(i[idx]) || return false
    end
    return true
end

blockcheckbounds(A::AbstractArray{T, N}, i::Block{N}) where {T,N} = blockcheckbounds(A, i.n...)
blockcheckbounds(A::AbstractArray{T, N}, i::Vararg{Block{1},N}) where {T,N} = blockcheckbounds(A, Int.(i)...)
blockcheckbounds(A::AbstractVector{T}, i::Block{1}) where {T,N} = blockcheckbounds(A, Int(i))

@propagate_inbounds Base.setindex!(block_arr::AbstractBlockArray{T,N}, v, block::Block{N}) where {T,N} =
    setindex!(block_arr, v, Block.(block.n)...)
@inline @propagate_inbounds function Base.setindex!(block_arr::AbstractBlockArray{T,N}, v, block::Vararg{Block{1}, N}) where {T,N}
    blockcheckbounds(block_arr, block...)
    dest = view(block_arr, block...)
    size(dest) == size(v) || throw(DimensionMismatch(string("tried to assign $(size(v)) array to $(size(dest)) block")))
    copyto!(dest, v)
    block_arr
end

@inline @propagate_inbounds Base.setindex!(block_arr::AbstractBlockArray{T,N}, v, blockindex::BlockIndex{N}) where {T,N} =
    view(block_arr, block(blockindex))[blockindex.α...] = v
@inline @propagate_inbounds Base.setindex!(block_arr::AbstractBlockVector{T}, v, blockindex::BlockIndex{1}) where {T} =
    view(block_arr, block(blockindex))[blockindex.α...] = v
@inline @propagate_inbounds Base.setindex!(block_arr::AbstractBlockArray{T,N}, v, blockindex::Vararg{BlockIndex{1},N}) where {T,N} =
    block_arr[BlockIndex(blockindex)] = v

viewblock(block_arr, block) = Base.invoke(view, Tuple{AbstractArray, Any}, block_arr, block)
@inline Base.view(block_arr::AbstractBlockArray{<:Any,N}, block::Block{N}) where N = viewblock(block_arr, block)
@inline Base.view(block_arr::AbstractBlockArray{<:Any,N}, block::Block{1}) where N = view(block_arr, Block(Base._ind2sub(map(Base.OneTo, blocksize(block_arr)), Int(block))...))
@inline Base.view(block_arr::AbstractBlockVector, block::Block{1}) = viewblock(block_arr, block)
@inline @propagate_inbounds Base.view(block_arr::AbstractBlockArray, block::Block{1}...) = view(block_arr, Block(block))

"""
    eachblock(A::AbstractBlockArray)

Create a generator that iterates over each block of an `AbstractBlockArray`
returning views.

```jldoctest; setup = quote using BlockArrays end
julia> v = Array(reshape(1:6, (2, 3)))
2×3 Matrix{Int64}:
 1  3  5
 2  4  6

julia> A = BlockArray(v, [1,1], [2,1])
2×2-blocked 2×3 BlockMatrix{Int64}:
 1  3  │  5
 ──────┼───
 2  4  │  6

julia> Matrix.(collect(eachblock(A)))
2×2 Matrix{Matrix{Int64}}:
 [1 3]  [5;;]
 [2 4]  [6;;]

julia> sum.(eachblock(A))
2×2 Matrix{Int64}:
 4  5
 6  6
```
"""
function eachblock(A::AbstractBlockArray)
    # blockinds = CartesianIndices(blocksize(A))
    blockinds = CartesianIndices(axes.(blocklasts.(axes(A)),1))
    (view(A, Block(Tuple(I))) for I in blockinds)
end

# Use memory laout for sub-blocks
@inline Base.getindex(A::AbstractMatrix, kr::Colon, jr::Block{1}) = ArrayLayouts.layout_getindex(A, kr, jr)
@inline Base.getindex(A::AbstractMatrix, kr::Block{1}, jr::Colon) = ArrayLayouts.layout_getindex(A, kr, jr)
@inline Base.getindex(A::AbstractMatrix, kr::Block{1}, jr::AbstractVector) = ArrayLayouts.layout_getindex(A, kr, jr)
@inline Base.getindex(A::AbstractArray{T,N}, kr::Block{1}, jrs...) where {T,N} = ArrayLayouts.layout_getindex(A, kr, jrs...)
@inline Base.getindex(A::AbstractArray{T,N}, block::Block{N}) where {T,N} = ArrayLayouts.layout_getindex(A, block)
@inline Base.getindex(A::AbstractMatrix, kr::AbstractVector, jr::Block) = ArrayLayouts.layout_getindex(A, kr, jr)
@inline Base.getindex(A::AbstractMatrix, kr::BlockRange{1}, jr::BlockRange{1}) = ArrayLayouts.layout_getindex(A, kr, jr)
@inline Base.getindex(A::LayoutMatrix, kr::BlockRange{1}, jr::BlockRange{1}) = ArrayLayouts.layout_getindex(A, kr, jr)
for Typ in (:AbstractTriangular, :Adjoint, :Transpose, :Symmetric, :Hermitian)
    @eval @inline Base.getindex(A::$Typ{<:Any,<:LayoutMatrix}, kr::BlockRange{1}, jr::BlockRange{1}) = ArrayLayouts.layout_getindex(A, kr, jr)
end

###
# permutedims
#
# just use transpose for now
###

Base.permutedims(A::AbstractBlockVector{<:Number}) = transpose(A)
Base.permutedims(A::AbstractBlockMatrix{<:Number}) = transpose(A)


text:
using Test
using OpenQuantumSystems
using Random, SparseArrays, LinearAlgebra, StableRNGs
import QuantumOpticsBase

@testset "memory kernel" begin

    Random.seed!(StableRNG(0), 1)

    D(op1::Array, op2::Array) = abs(norm(op1 - op2))
    D(x1::StateVector, x2::StateVector) = norm(x2 - x1)
    D(op1::AbstractOperator, op2::AbstractOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))
    D(op1::AbstractSuperOperator, op2::AbstractSuperOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))

    # TODO: change to macro
    mode1 = Mode(0.2, 1.0)
    mode2 = Mode(0.3, 2.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1], 2, [2.0, 200.0])
    mol2 = Molecule([mode2], 2, [3.0, 300.0])
    aggCore = AggregateCore([mol1, mol2])
    aggCore.coupling[2, 3] = 50
    aggCore.coupling[3, 2] = 50
    agg = setupAggregate(aggCore)
    aggTools = agg.tools
    aggOperators = agg.operators

    Ham_I = agg.operators.Ham_I
    Ham_0 = agg.operators.Ham_0
    Ham = agg.operators.Ham

    basis = agg.tools.basis
    indicesLen = agg.tools.bSize
    indices = agg.tools.indices
    indicesMap = agg.tools.indicesMap
    FCFact = agg.tools.FCfactors
    FCProd = agg.tools.FCproduct

    Ham_0_lambda, Ham_0_S = eigen(Ham_0.data)
    Ham_0_Sinv = inv(Ham_0_S)
    Ham_0_lambda = diagm(Ham_0_lambda)

    # T = 300
    # W0 = thermal_state_composite(T, [0.0, 0.8, 0.2], Ham, indices; diagonalize=false, diagonal=true)
    data = Matrix(Hermitian(rand(ComplexF64, indicesLen, indicesLen)))
    W0 = DenseOperator(basis, basis, data)
    normalize!(W0)
    W0_bath = get_rho_bath(
        W0,
        aggCore,
        aggTools;
        justCopy = true,
    )
    rho0 = trace_bath(W0.data, aggCore, aggTools)

    W_ab = take_el_part(W0.data, 1, 1, indicesMap)
    W_ab_len = size(W_ab, 1)
    @test 1e-15 > D(W_ab, W0.data[1:W_ab_len, 1:W_ab_len])

    tspan = [0.0:0.1:0.3;]
    for t_i = 1:length(tspan)
        for s_i = 1:t_i
            t = tspan[t_i]
            s = tspan[s_i]
            # println(t, " ", s)
            H_II_t = getInteractionHamIPictureA(Ham_I.data, Ham_0_lambda, Ham_0_S, Ham_0_Sinv, t)
            H_II_s = getInteractionHamIPictureA(Ham_I.data, Ham_0_lambda, Ham_0_S, Ham_0_Sinv, s)

            MemoryKernel_1 = MemoryKernel_1_traced(
                H_II_t,
                H_II_s,
                W0_bath.data,
                aggCore,
                aggTools
            )
            rho = zero(complex(rho0))
            for a = 1:size(rho, 1),
                b = 1:size(rho, 2),
                c = 1:size(rho, 1),
                d = 1:size(rho, 2)

                rho[a, b] += MemoryKernel_1[a, b, c, d] * rho0[c, d]
            end
            MemoryKernel_1_ref = H_II_t * H_II_s * W0.data
            rho_ref = trace_bath(
                MemoryKernel_1_ref,
                aggCore,
                aggTools
            )
            @test 1e-6 > D(rho, rho_ref)
            # println(D(rho, rho_ref))

            MemoryKernel_2 = MemoryKernel_2_traced(
                H_II_t,
                H_II_s,
                W0_bath.data,
                aggCore,
                aggTools
            )
            rho = zero(complex(rho0))
            for a = 1:size(rho, 1),
                b = 1:size(rho, 2),
                c = 1:size(rho, 1),
                d = 1:size(rho, 2)

                rho[a, b] += MemoryKernel_2[a, b, c, d] * rho0[c, d]
            end
            MemoryKernel_2_ref = H_II_t * W0.data * H_II_s
            rho_ref = trace_bath(
                MemoryKernel_2_ref,
                aggCore,
                aggTools
            )
            @test 1e-6 > D(rho, rho_ref)
            # println(D(rho, rho_ref))

            MemoryKernel_3 = MemoryKernel_3_traced(
                H_II_t,
                H_II_s,
                W0_bath.data,
                aggCore,
                aggTools
            )
            rho = zero(complex(rho0))
            for a = 1:size(rho, 1),
                b = 1:size(rho, 2),
                c = 1:size(rho, 1),
                d = 1:size(rho, 2)

                rho[a, b] += MemoryKernel_3[a, b, c, d] * rho0[c, d]
            end
            MemoryKernel_3_ref = H_II_s * W0.data * H_II_t
            rho_ref = trace_bath(
                MemoryKernel_3_ref,
                aggCore,
                aggTools
            )
            @test 1e-6 > D(rho, rho_ref)
            # println(D(rho, rho_ref))

            MemoryKernel_4 = MemoryKernel_4_traced(
                H_II_t,
                H_II_s,
                W0_bath.data,
                aggCore,
                aggTools
            )
            rho = zero(complex(rho0))
            for a = 1:size(rho, 1),
                b = 1:size(rho, 2),
                c = 1:size(rho, 1),
                d = 1:size(rho, 2)

                rho[a, b] += MemoryKernel_4[a, b, c, d] * rho0[c, d]
            end
            MemoryKernel_4_ref = W0.data * H_II_s * H_II_t
            rho_ref = trace_bath(
                MemoryKernel_4_ref,
                aggCore,
                aggTools
            )
            @test 1e-6 > D(rho, rho_ref)
            # println(D(rho, rho_ref))
        end
    end

    tspan = [0.0:0.1:0.3;]
    for t_i = 1:length(tspan)
        for s_i = 1:t_i
            t = tspan[t_i]
            s = tspan[s_i]
            H_II_t = getInteractionHamIPictureA(Ham_I.data, Ham_0_lambda, Ham_0_S, Ham_0_Sinv, t)
            H_II_s = getInteractionHamIPictureA(Ham_I.data, Ham_0_lambda, Ham_0_S, Ham_0_Sinv, s)

            MemoryKernel = MemoryKernel_traced(
                H_II_t,
                H_II_s,
                W0_bath.data,
                aggCore,
                aggTools
            )
            rho = zero(complex(rho0))
            for a = 1:size(rho, 1),
                b = 1:size(rho, 2),
                c = 1:size(rho, 1),
                d = 1:size(rho, 2)

                rho[a, b] += MemoryKernel[a, b, c, d] * rho0[c, d]
            end
            MemoryKernel_ref = H_II_s * W0.data - W0.data * H_II_s
            MemoryKernel_ref = H_II_t * MemoryKernel_ref - MemoryKernel_ref * H_II_t
            rho_ref = trace_bath(
                MemoryKernel_ref,
                aggCore,
                aggTools
            )
            @test 1e-6 > D(rho, rho_ref)
        end
    end

end # testset

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
module DiskTracer

# package code goes here
export constants,
    model

# These statements just straight up dump the source code directly here, making DiskJockey.jl
# act as one giant file.
include("constants.jl")
include("model.jl")
include("geometry.jl")
include("trace.jl")
include("trace_simple.jl")

end # module


text:
#=
ArithematicAnalysis:
- Julia version: 1.4.2
- Author: akhil
- Date: 2020-05-28
=#
module ArithematicAnalysis
    include("Bisection.jl")
    include("GaussianElimination.jl")
    include("StaticEquilibrium.jl")
end


text:
using Test

using AutomotiveSimulator
using EzXML
using AutoUrban

include("test_convert2xodr.jl")
# broken
# include("test_drivermodel.jl")

max_stars_repo_name: sisl/UrbanDrivingSimulation

max_issues_repo_name: sisl/UrbanDrivingSimulation

max_forks_repo_name: sisl/UrbanDrivingSimulation


text:
mutable struct SchellingAgent <: AbstractAgent
    id::Int # The identifier number of the agent
    pos::Dims{2} # The x, y location of the agent on a 2D grid
    mood::Bool # whether the agent is happy in its position. (true = happy)
    group::Int # The group of the agent,  determines mood as it interacts with neighbors
end

"""
``` julia
schelling(;
    numagents = 320,
    griddims = (20, 20),
    min_to_be_happy = 3
)
```
Same as in [Schelling's segregation model](@ref).
"""
function schelling(; numagents = 320, griddims = (20, 20), min_to_be_happy = 3)
    @assert numagents < prod(griddims)
    space = GridSpace(griddims, periodic = false)
    properties = Dict(:min_to_be_happy => min_to_be_happy)
    model =
        ABM(SchellingAgent, space; properties = properties, scheduler = random_activation)
    for n in 1:numagents
        agent = SchellingAgent(n, (1, 1), false, n < numagents / 2 ? 1 : 2)
        add_agent_single!(agent, model)
    end
    return model, schelling_agent_step!, dummystep
end

function schelling_agent_step!(agent, model)
    agent.mood == true && return # do nothing if already happy
    count_neighbors_same_group = 0
    for neighbor in nearby_agents(agent, model)
        if agent.group == neighbor.group
            count_neighbors_same_group += 1
        end
    end
    if count_neighbors_same_group ≥ model.min_to_be_happy
        agent.mood = true
    else
        move_agent_single!(agent, model)
    end
end


text:
# Adaptation of BayesLogit C code by Nicholas Polson, James Scott, Jesse Windle, 2012-2019

# double rtinvchi2(double h, double trunc, RNG& r)
#
#   double h2 = h * h
#   double R = trunc / h2
#   double X = 0.0
#   # I need to consider using a different truncated normal sampler.
#   double E1 = r.expon_rate(1.0) double E2 = r.expon_rate(1.0)
#   while ( (E1*E1) > (2 * E2 / R))
#     # printf("E %g %g %g %g\n", E1, E2, E1*E1, 2*E2/R)
#     E1 = r.expon_rate(1.0) E2 = r.expon_rate(1.0)
#   end
#   # printf("E %g %g \n", E1, E2)
#   X = 1 + E1 * R
#   X = R / (X * X)
#   X = h2 * X
#   return X
# end

const PISQ = __PI^2

const trunc_schedule = [
    0.64,
    0.68,
    0.72,
    0.75,
    0.78,
    0.8,
    0.83,
    0.85,
    0.87,
    0.89,
    0.91,
    0.93,
    0.95,
    0.96,
    0.98,
    1,
    1.01,
    1.03,
    1.04,
    1.06,
    1.07,
    1.09,
    1.1,
    1.12,
    1.13,
    1.15,
    1.16,
    1.17,
    1.19,
    1.2,
    1.21,
    1.23,
    1.24,
    1.25,
    1.26,
    1.28,
    1.29,
    1.3,
    1.32,
    1.33,
    1.34,
    1.35,
    1.36,
    1.38,
    1.39,
    1.4,
    1.41,
    1.42,
    1.44,
    1.45,
    1.46,
    1.47,
    1.48,
    1.5,
    1.51,
    1.52,
    1.53,
    1.54,
    1.55,
    1.57,
    1.58,
    1.59,
    1.6,
    1.61,
    1.62,
    1.63,
    1.65,
    1.66,
    1.67,
    1.68,
    1.69,
    1.7,
    1.71,
    1.72,
    1.74,
    1.75,
    1.76,
    1.77,
    1.78,
    1.79,
    1.8,
    1.81,
    1.82,
    1.84,
    1.85,
    1.86,
    1.87,
    1.88,
    1.89,
    1.9,
    1.91,
    1.92,
    1.93,
    1.95,
    1.96,
    1.97,
    1.98,
    1.99,
    2,
    2.01,
    2.02,
    2.03,
    2.04,
    2.05,
    2.07,
    2.08,
    2.09,
    2.1,
    2.11,
    2.12,
    2.13,
    2.14,
    2.15,
    2.16,
    2.17,
    2.18,
    2.19,
    2.21,
    2.22,
    2.23,
    2.24,
    2.25,
    2.26,
    2.27,
    2.28,
    2.29,
    2.3,
    2.31,
    2.32,
    2.33,
    2.35,
    2.36,
    2.37,
    2.38,
    2.39,
    2.4,
    2.41,
    2.42,
    2.43,
    2.44,
    2.45,
    2.46,
    2.47,
    2.48,
    2.49,
    2.51,
    2.52,
    2.53,
    2.54,
    2.55,
    2.56,
    2.57,
    2.58,
    2.59,
    2.6,
    2.61,
    2.62,
    2.63,
    2.64,
    2.65,
    2.66,
    2.68,
    2.69,
    2.7,
    2.71,
    2.72,
    2.73,
    2.74,
    2.75,
    2.76,
    2.77,
    2.78,
    2.79,
    2.8,
    2.81,
    2.82,
    2.83,
    2.84,
    2.85,
    2.87,
    2.88,
    2.89,
    2.9,
    2.91,
    2.92,
    2.93,
    2.94,
    2.95,
    2.96,
    2.97,
    2.98,
    2.99,
    3,
    3.01,
    3.02,
    3.03,
    3.04,
    3.06,
    3.07,
    3.08,
    3.09,
    3.1,
    3.11,
    3.12,
    3.13,
    3.14,
    3.15,
    3.16,
    3.17,
    3.18,
    3.19,
    3.2,
    3.21,
    3.22,
    3.23,
    3.24,
    3.25,
    3.27,
    3.28,
    3.29,
    3.3,
    3.31,
    3.32,
    3.33,
    3.34,
    3.35,
    3.36,
    3.37,
    3.38,
    3.39,
    3.4,
    3.41,
    3.42,
    3.43,
    3.44,
    3.45,
    3.46,
    3.47,
    3.49,
    3.5,
    3.51,
    3.52,
    3.53,
    3.54,
    3.55,
    3.56,
    3.57,
    3.58,
    3.59,
    3.6,
    3.61,
    3.62,
    3.63,
    3.64,
    3.65,
    3.66,
    3.67,
    3.68,
    3.69,
    3.71,
    3.72,
    3.73,
    3.74,
    3.75,
    3.76,
    3.77,
    3.78,
    3.79,
    3.8,
    3.81,
    3.82,
    3.83,
    3.84,
    3.85,
    3.86,
    3.87,
    3.88,
    3.89,
    3.9,
    3.91,
    3.92,
    3.93,
    3.95,
    3.96,
    3.97,
    3.98,
    3.99,
    4,
    4.01,
    4.02,
    4.03,
    4.04,
    4.05,
    4.06,
    4.07,
    4.08,
    4.09,
    4.1,
    4.11,
    4.12,
    4.13,
]

struct PolyaGammaAltSampler <:
       Distributions.Sampleable{Distributions.Univariate,Distributions.Continuous}
    h::Float64
    z::Float64
    maxinner::Int64
end

function a_coef(n::Int64, x::Float64, h::Float64)
    d_n = 2.0 * n + h
    log_out =
        h * _log_c(2.0) - loggamma(h) + loggamma(n + h) - loggamma(n + 1) + _log_c(d_n) -
        0.5 * _log_c(2.0 * __PI * (x^3)) - 0.5 * d_n * d_n / x
    out = _exp_c(log_out)
    return out::Float64
end

function a_coef_recursive(
    n::Float64,
    x::Float64,
    h::Float64,
    coef_h::Float64,
    gnh_over_gn1_gh::Float64,
)
    d_n = 2.0 * n + h
    # gamma_nh_over_n *= (n + h - 1) / n  # Can speed up further by separate function for a0 and an, n > 0.
    if (n != 0)
        gnh_over_gn1_gh *= (n + h - 1) / n
    else
        gnh_over_gn1_gh = 1.0
    end
    coef = coef_h * gnh_over_gn1_gh
    log_kernel = -0.5 * (_log_c(x * x * x) + d_n * d_n / x) + (d_n)
    return coef * _exp_c(log_kernel)::Float64
end

function w_left(x::Float64, h::Float64, z::Float64)
    if (z != 0)
        # p.13 l(x|h)
        return (1 + _exp_c(-2 * abs(z)))^h *
               Distributions.pdf(Distributions.InverseGaussian(h / z, h * h), x)
    else
        # p.11 l(x|h)
        return 2^h * Distributions.pdf(Distributions.InverseGamma(0.5, 0.5 * h * h), x)
    end
end

function w_right(x::Float64, h::Float64, z::Float64)
    if (z != 0)
        # p.13 r(x|h)
        lambda_z = (PISQ + 4 * (z^2)) / 8
        return (__PI / (2 * lambda_z))^h *
               Distributions.pdf(Distributions.Gamma(h, 1 / lambda_z), x)
    else
        # p.11 r(x|h)
        return (4 / __PI)^h * Distributions.pdf(Distributions.Gamma(h, 1 / (PISQ / 8)), x)
    end
    #lambda_z = __PI^2 * 0.125 + 0.5 * z * z
    #p = _exp_c(h * _log_c((__PI/2) / lambda_z)) * (1.0-p_gamma_rate(trunc, h, lambda_z, false))
end

function g_tilde(x::Float64, h::Float64, trunc::Float64)
    if (x > trunc)
        return _exp_c(
            h * _log_c(0.5 * __PI) + (h - 1) * _log_c(x) - PISQ * 0.125 * x - loggamma(h),
        )
    else
        return h * _exp_c(
            h * _log_c(2.0) - 0.5 * _log_c(2.0 * __PI * x * x * x) - 0.5 * h * h / x,
        )
    end
    # out = h * pow(2, h) * pow(2 * __PI * pow(x,3), -0.5) * _exp_c(-0.5 * pow(h,2) / x)
end

########################################
# Sample #
########################################

function draw_abridged(
    rng::Distributions.AbstractRNG,
    h::Float64,
    z::Float64,
    max_inner::Int64,
)
    # p.132 Windle(2013)
    if (h < 1 || h > 4)
        println("draw h = $h must be in [1,4]")
        return 0
    end

    # Change the parameter.
    z = abs(z) * 0.5

    trunc = trunc_schedule[Int64(floor((h - 1.0) * 100))]

    # Now sample 0.25 * J^*(1, z := z/2).

    lambda_z = (PISQ + 4 * (z^2)) / 8
    weight_left = w_left(trunc, h, z)
    weight_right = w_right(trunc, h, z)
    prob_right = weight_right / (weight_right + weight_left)

    # printf("prob_right: %g\n", prob_right)

    coef1_h = _exp_c(h * _log_c(2.0) - 0.5 * _log_c(2.0 * __PI))
    # double gamma_nh_over_n = RNG::Gamma(h)
    gnh_over_gn1_gh = 1.0 # Will fill in value on first call to a_coef_recursive.

    num_trials = 0
    total_iter = 0

    while (num_trials < 10000)
        num_trials += 1

        X = 0.0
        Y = 0.0

        # if (Distributions.rand(rng) < p/(p+q))
        uu = Distributions.rand(rng) #runif(1)
        if (uu < prob_right)
            X = Distributions.rand(rng, TruncatedGamma(h, 1 / lambda_z, trunc, Inf))
        else
            X = Distributions.rand(rng, TruncatedInverseGaussian(h / z, h^2, 0.0, trunc))
        end
        # double S  = a_coef(0, X, h)
        S = a_coef_recursive(0.0, X, h, coef1_h, gnh_over_gn1_gh)
        a_n = copy(S)
        gt = g_tilde(X, h, trunc)
        Y = Distributions.rand(rng) * gt
        decreasing = false

        n = 0
        go = true
        # Cap the number of iterations?
        while (go && n <= max_inner)
            total_iter += 1
            n += 1
            # Break infinite loop.  Put first so it always checks n==0.
            prev = copy(a_n)
            # a_n  = a_coef(n, X, h)
            a_n = a_coef_recursive(Float64(n), X, h, coef1_h, gnh_over_gn1_gh)
            # printf("a_n=%g, a_n2=%g\n", a_n, a_n2)
            decreasing = a_n <= prev

            if ((n % 2) == 1)
                S = S - a_n
                if (Y <= S && decreasing)
                    return X
                end
            else
                S = S + a_n
            end
            if (Y > S && decreasing)
                go = false
            end
        end
        # Need Y <= S in event that Y = S, e.g. when X = 0.
    end
    # We should never get here.
    return -1.0
end # draw

function Distributions.rand(rng::Distributions.AbstractRNG, s::PolyaGammaAltSampler)
    if (s.h < 1)
        println("draw h = ", s.h, " must be >= 1")
        return 0
    end
    n = floor((s.h - 1.0) / 4.0)
    remain = s.h - 4.0 * n
    x = 0.0
    for i = 1:n
        x += draw_abridged(rng, 4.0, s.z, s.maxinner)
        if (remain > 4.0)
            x += (
                draw_abridged(rng, 0.5 * remain, s.z, s.maxinner) +
                draw_abridged(rng, 0.5 * remain, s.z, s.maxinner)
            )
        else
            x += draw_abridged(rng, remain, s.z, s.maxinner)
        end
    end
    return x
end

########################################
# APPENDIX #
########################################

# We should only have to calculate Gamma(h) once.  We can then get Gamma(n+h)
# from the recursion Gamma(z+1) = z Gamma(z).  Not sure how that is in terms of
# stability, but that should save us quite a few computations.  This affects
# a_coef and g_tilde.


text:
using Random, Test, StructLH

mutable struct MoaTest
    x
    y
    z
    a
end

function merge_object_arrays_test()
    @testset "Merge object arrays" begin
        rng = MersenneTwister(123);
        # First dimension must be the same for all arrays
        n = 17;
        x = rand(rng, n, 2);
        y = rand(rng, n, 3);
        z = rand(rng, n, 3, 2);
        a = rand(rng, n);
        oTg = MoaTest(x, y, z, a);

        idxV = [3,4,8];
        nSrc = length(idxV);
        x2 = rand(rng, nSrc, 2);
        y2 = 9.3;
        z2 = rand(rng, nSrc, 3, 2);
        a2 = rand(rng, nSrc);
        oSrc = MoaTest(x2, y2, z2, a2);

        merge_object_arrays!(oSrc, oTg, idxV, false, dbg = true);

        @test oTg.x[idxV,:] ≈ x2
        @test oTg.y ≈ y  # Skipped field; not of array type in source
        @test oTg.z[idxV,:,:] ≈ z2
        @test oTg.a[idxV] ≈ a2
    end
end

@testset "Merge" begin
    merge_object_arrays_test()
end

# -------------


text:
#-------------binary_trees.jl--------------------------------------------------#
#
# Purpose: This implements a basica binary tree, DFS, and a Priority Queue
#
#------------------------------------------------------------------------------#

# Now we need a typealias for the tree
abstract BT

# Type of Binary Tree
type BTree <: BT
    weight::Float64
    bitpath::String
    left::BT
    right::BT
end

# First, we need to create a type for empty nodes
type Empty <: BT end
const et = Empty()

# Depth-First Search
function DFS(node::BTree)
    println(node.bitpath)
    if !isa(node.right, Empty)
        #println(0)
        DFS(node.right)
    end

    if !isa(node.left, Empty)
        #println(1)
        DFS(node.left)
    end

    if isa(node.left, Empty) && isa(node.right, Empty)
        println(node.bitpath)
    end
end

# Function to create binary tree with a certain number of generations
function create_BTree(node::BTree, gen::Int64, max_gen::Int64, bitpath::String)

    node.weight = rand()
    node.bitpath = bitpath

    if gen < max_gen
        node.right = BTree(0, "", et, et)
        node.left = BTree(0, "", et, et)
        create_BTree(node.right, gen+1, max_gen, string(bitpath,1))
        create_BTree(node.left, gen+1, max_gen, string(bitpath,0))
    end

end

function main()

    tree = BTree(0, "", et, et,)
    create_BTree(tree, 0, 3, "")
    DFS(tree)

end

main()

max_stars_repo_name: vchuravy/skillpill-julia

max_issues_repo_name: vchuravy/skillpill-julia

max_forks_repo_name: vchuravy/skillpill-julia


text:
export GradientDescent, SecondMomentMatching, GramSchmidt, QRFact, SVDFact, orthonormalize!
# QR, SVD not exported due to conflict with LinearAlgebra identifiers.

struct GradientDescent
    maxiter::Int
    μ::Float64
    atol::Float64
    rtol::Float64
    function GradientDescent(;maxiter = 100, μ = 1.0, atol = 1e-8, rtol = 1e-8)
        return new(maxiter, μ, atol, rtol)
    end
end

struct QRFact end

struct SVDFact end

struct SecondMomentMatching end

struct GramSchmidt end

function orthonormalize!(U::AbstractMatrix, alg::GradientDescent)
    @unpack μ, atol, rtol, maxiter = alg
    r = size(U,2)
    K = U'*U
    A = Matrix{eltype(K)}(I, r, r)
    dA = similar(A)
    iter = 0
    ϵ = norm(A'*K*A - I)^2
    while ϵ > atol && ϵ > r*rtol && iter < maxiter
        dA .= - K*A*(A'*K*A - I)
        A .+= μ*dA
        iter += 1
        ϵ = norm(A'*K*A - I)^2
    end
    if iter == maxiter
        @warn "Gradient flow orthonormalization did not converge.
               Iterations exceeded maxiters = $maxiter.
               Primal residual: $(norm(A'*K*A - I)^2)"
    end
    U .= U*A
end

function orthonormalize!(U::AbstractMatrix, Z::AbstractMatrix, alg::GradientDescent)
    @unpack μ, atol, rtol, maxiter = alg
    K = U'*U
    r = size(U,2)
    A = Matrix{eltype(K)}(I, r, r)
    Ainv = Matrix{eltype(K)}(I, r, r)
    dA = similar(A)
    iter = 0
    ϵ = norm(A'*K*A - I)
    while ϵ > atol && ϵ > r*rtol && iter < maxiter
        dA .= - K*A*(A'*K*A - I)
        A .+= μ*dA
        Ainv .-= μ*Ainv*dA*Ainv
        iter += 1
        ϵ = norm(A'*K*A - I)
    end
    if iter == maxiter
        @warn "Gradient flow orthonormalization did not converge.
               Iterations exceeded maxiters = $maxiter.
               Primal residual: $(norm(A'*K*A - I)^2)"
    end
    U .= U*A
    Z .= Z*Ainv'
end

function orthonormalize!(U::AbstractMatrix, ::SecondMomentMatching)
    P = Symmetric(U'*U)
    Λ, V = eigen(P)
    U .= U*V*Diagonal(1 ./ sqrt.(Λ)) * V'
end

function orthonormalize!(U::AbstractMatrix, Z::AbstractMatrix, ::SecondMomentMatching)
    P = Symmetric(U'*U)
    Λ, V = eigen(P)
    U .= U*V*Diagonal(1 ./ sqrt.(Λ)) * V'
    Z .= Z*V*Diagonal(sqrt.(Λ)) * V'
end

function orthonormalize!(U::AbstractMatrix, ::QRFact)
    Q, _ = qr!(U)
    U .= Matrix(Q)
end

function orthonormalize!(U::AbstractMatrix, Z::AbstractMatrix, ::QRFact)
    Q, R = qr!(U)
    U .= Matrix(Q)
    Z .= Z*R'
end

function orthonormalize!(U::AbstractMatrix, ::SVDFact)
    Q, _, _ = svd!(U)
    U .= Q
end

function orthonormalize!(U::AbstractMatrix, Z::AbstractMatrix, ::SVDFact)
    Q, S, V = svd!(U)
    U .= Q
    Z .= Z*V*Diagonal(S)'
end

function orthonormalize!(U::AbstractMatrix, ::GramSchmidt)
    u = @view U[:,1]
    u ./= norm(u)
    for i in 2:size(U,2)
        u = @view U[:,i]
        subU = @view U[:,1:i-1]
        ips = subU'*u
        mul!(u,subU,ips,-1,1)
        u ./= norm(u)
    end
end

function orthonormalize!(U::AbstractMatrix, Z::AbstractMatrix, ::GramSchmidt)
    u = @view U[:,1]
    z = @view Z[:,1]
    r = norm(u)
    u ./= r
    z .*= r
    for i in 2:size(U,2)
        u = @view U[:,i]
        z = @view Z[:,i]
        subU = @view U[:,1:i-1]
        subZ = @view Z[:,1:i-1]
        ips = subU'*u
        mul!(u, subU, ips, -1, 1)
        mul!(subZ, z, ips', 1, 1)
        r = norm(u)
        u ./= r
        z .*= r
    end
end

# this implementation could return in the inverse of the GramSchmidt procedure in R
# function orthonormalize!(U, Z, ::GramSchmidt)
#     v = @view U[:,1]
#     R = UpperDiagonal(zeros(size(U,2), size(U,2)))
#     R[1,1] = norm(v)
#     v ./= R[1,1]
#     for i in 2:size(U,2)
#         v = @view U[:,i]
#         subU = @view U[:,1:i-1]
#         u = subU'*v
#         R[1:i-1,i] = u
#         R[i,i] = norm(v)
#         mul!(v,subU,u,-1,1)
#         v .*= R[i,i]
#     end
#     Z .= Z*R'
# end

# ToDo:
#  1. add general minimum second order moment matching orthonormalization:
#         https://link.springer.com/content/pdf/10.1007/s00211-021-01178-8.pdf


text:
module Bias

using LinearAlgebra
using Statistics
using Random
using SparseArrays

include("word_sets.jl")

export WEAT_WORD_SETS

function get_weat_idx_set(word_set::NamedTuple, vocab::Dict)
    words2indices(words) = [vocab[w].index for w in words if w in keys(vocab)]
    return (S=words2indices(word_set.S), T=words2indices(word_set.T),
            A=words2indices(word_set.A), B=words2indices(word_set.B))
end

function normalize_rows(X::AbstractArray)
    return mapslices(normalize, X, dims=2)
end

function normalize_rows(X::SparseMatrixCSC)
    return mapslices(normalize, X, dims=1)'
end

function effect_size(S::AbstractArray, T::AbstractArray, A::AbstractArray,
        B::AbstractArray)
    Ŝ = normalize_rows(S)
    T̂ = normalize_rows(T)
    Â = normalize_rows(A)
    B̂ = normalize_rows(B)

    μSA = mean(Ŝ * Â', dims=2)
    μSB = mean(Ŝ * B̂', dims=2)
    μTA = mean(T̂ * Â', dims=2)
    μTB = mean(T̂ * B̂', dims=2)

    dS = μSA - μSB
    dT = μTA - μTB
    return (mean(dS) - mean(dT)) / std(vcat(dS, dT))
end

# Helper to grab word vecs for you
function effect_size(W, weat_idx_set::NamedTuple)
    S = W[weat_idx_set.S, :]
    T = W[weat_idx_set.T, :]
    A = W[weat_idx_set.A, :]
    B = W[weat_idx_set.B, :]
    return effect_size(S, T, A, B)
end

function effect_size(X::SparseMatrixCSC, weat_idx_set::NamedTuple)
    S = X[:, weat_idx_set.S]
    T = X[:, weat_idx_set.T]
    A = X[:, weat_idx_set.A]
    B = X[:, weat_idx_set.B]
    return effect_size(S, T, A, B)
end

# Helper to compute effect size after changes to the embedding
function effect_size(W, weat_idx_set::NamedTuple, deltas::Dict)
    weat_vec_set = []
    delta_indices = keys(deltas) # the indices that have changes
    for indices in weat_idx_set
        vecs = W[indices, :]
        for (idx, pos) = zip(delta_indices, indexin(delta_indices, indices))
            # idx: word index of changed vectors
            # pos: relative position of that index in the "vecs" matrix
            if pos != nothing
                vecs[pos, :] += deltas[idx]
            end
        end
        push!(weat_vec_set, vecs)
    end
    return effect_size(weat_vec_set...)
end

# Compute p-value of the bias
function p_value(W, word_indices, N=10_000)
    es = effect_size(W, word_indices)
    ST = vcat(word_indices.S, word_indices.T)
    boundary = length(word_indices.S)
    trials = zeros(N)
    for i in 1:N
        perm = randperm(length(ST))
        S = ST[perm[1:boundary]]
        T = ST[perm[boundary+1:end]]
        trials[i] = effect_size(W, (S=S, T=T, A=word_indices.A, B=word_indices.B))
    end
    return mean(trials .> es)
end

end


text:
import Base: @deprecate

# to be deprecated
const FINALIZE_ATEXIT = Ref(true)
@deprecate finalize_atexit() true false

@deprecate Info(::Cint) Info() false
@deprecate Info_set(info::Info,key::AbstractString,value::AbstractString) info[Symbol(key)] = value false
@deprecate Info_get(info::Info,key::AbstractString) info[Symbol(key)] false
@deprecate Info_delete(info::Info,key::AbstractString) delete!(info, Symbol(key)) false
@deprecate Info_free(info::Info) free(info) false

@deprecate Comm_free(comm::Comm) free(comm) false
@deprecate Type_Commit(dtyp) Type_Commit!(dtyp) false

@deprecate Cart_coords!(comm::Comm, maxdims::Integer) Cart_coords(comm, maxdims) false

# window functions
@deprecate(Win_create(base::Array{T}, info::Info, comm::Comm, win::Win) where {T},
           (win = Win_create(base, comm; info...)), false)

@deprecate(Win_create_dynamic(info::Info, comm::Comm, win::Win),
           (win = Win_create_dynamic(comm; info...)), false)

@deprecate(Win_allocate_shared(::Type{T}, len::Int, info::Info, comm::Comm, win::Win) where T,
           ((win, ptr) = Win_allocate_shared(T, len, comm; info...); ptr), false)

@deprecate Win_free(win::Win) free(win) false