text:
"""
    TimeStep!(model::PlanktonModel, ΔT::Int64, diags::PlanktonDiagnostics, resultpath::String)
Update physiology processes and nutrient field of `PlanktonModel` one time step forward.

Keyword Arguments
=================
- `model`: `PlanktonModel` to be updated one time step forward.
- `ΔT`: The length of a time step.
- `diags`: `PlanktonDiagnostics` to be updated.
- `resultpath` (optional): The file path to store model output.
"""
function TimeStep!(model::PlanktonModel, ΔT, diags::PlanktonDiagnostics)
    model.t = model.t+ΔT
    model.iteration = model.iteration+1

    @inbounds model.timestepper.vel½.u.data .= (model.timestepper.vel₀.u.data .+ model.timestepper.vel₁.u.data) .* 0.5
    @inbounds model.timestepper.vel½.v.data .= (model.timestepper.vel₀.v.data .+ model.timestepper.vel₁.v.data) .* 0.5
    @inbounds model.timestepper.vel½.w.data .= (model.timestepper.vel₀.w.data .+ model.timestepper.vel₁.w.data) .* 0.5

    zero_fields!(model.timestepper.plk)
    @inbounds model.timestepper.Chl .= 0.0
    @inbounds model.timestepper.pop .= 0.0  # may add an option of self grazing, besides shared grazing
    ##### plankton advection and diffusion
    for sp in keys(model.individuals.phytos)
        ##### RK4
        plankton_advection!(model.individuals.phytos[sp].data, model.timestepper.velos, model.grid,
                               model.timestepper.vel₀, model.timestepper.vel½, model.timestepper.vel₁, ΔT, model.arch)
        ##### Diffusion
        plankton_diffusion!(model.individuals.phytos[sp].data, model.timestepper.rnd,
                            model.bgc_params["κhP"], ΔT, model.grid, model.arch)

        #### calculate accumulated Chla quantity (not concentration)
        find_inds!(model.individuals.phytos[sp].data, model.grid, model.arch)
        acc_counts!(model.timestepper.Chl, model.timestepper.pop,
                    model.individuals.phytos[sp].data.Chl, model.individuals.phytos[sp].data.ac,
                    model.individuals.phytos[sp].data.xi, model.individuals.phytos[sp].data.yi,
                    model.individuals.phytos[sp].data.zi, model.arch)
    end

    ##### calculate PAR
    calc_par!(model.timestepper.par, model.arch, model.timestepper.Chl, model.timestepper.PARF,
              model.grid, model.bgc_params["kc"], model.bgc_params["kw"])

    ##### plankton physiology
    for sp in keys(model.individuals.phytos)
        find_NPT!(model.timestepper.nuts, model.individuals.phytos[sp].data.xi,
                  model.individuals.phytos[sp].data.yi, model.individuals.phytos[sp].data.zi,
                  model.individuals.phytos[sp].data.ac, model.nutrients.NH4.data,
                  model.nutrients.NO3.data, model.nutrients.PO4.data, model.nutrients.DOC.data,
                  model.timestepper.par, model.timestepper.temp, model.timestepper.pop, model.arch)

        plankton_update!(model.individuals.phytos[sp].data, model.timestepper.nuts,
                             model.timestepper.rnd, model.individuals.phytos[sp].p,
                             model.timestepper.plk, diags.plankton[sp], ΔT, model.t, model.arch, model.mode)
    end

    ##### diagnostics for nutrients
    @inbounds diags.tracer.PAR .+= model.timestepper.par
    for key in keys(diags.tracer)
        if key in keys(model.nutrients)
            @inbounds diags.tracer[key] .+= model.nutrients[key].data
        end
    end

    nut_update!(model.nutrients, model.timestepper.Gcs, model.timestepper.nut_temp, model.arch,
                model.grid, model.bgc_params, model.timestepper.vel₁, model.timestepper.plk, ΔT, model.iteration)

    @inbounds model.timestepper.vel₀.u.data .= model.timestepper.vel₁.u.data
    @inbounds model.timestepper.vel₀.v.data .= model.timestepper.vel₁.v.data
    @inbounds model.timestepper.vel₀.w.data .= model.timestepper.vel₁.w.data

    return nothing
end


text:
using SimpleDirectMediaLayer
const SDL2 = SimpleDirectMediaLayer

using SimpleDirectMediaLayer:
    Renderer

export
    RenderContext,
    save!, restore!, post, pre, apply!,
    render_line, render_lines, render_aa_rectangle, render_outlined_polygon

using .Rendering
using ..Math
using ..Geometry

import ..Math.transform!

@enum RenderStyle FILLED OUTLINE BOTH
@enum PolygonState OPEN CLOSED

mutable struct RenderContext
    renderer::Ptr{Renderer}

    raster_font::RasterFont
    vector_font::VectorFont

    # State management
    state::Array{State}
    stack_top::Integer

    clear_color::Palette
    draw_color::Palette

    # Device/Window dimensions
    width::Int32
    height::Int32

    # The current transform based on top
    current::Math.AffineTransform

    post::Math.AffineTransform   # A preallocated cache

    # Scratch working variables
    v1::Geometry.Point{Float64}
    v2::Geometry.Point{Float64}

    function RenderContext(world::Ranger.World, raster_font::RasterFont, vector_font::Rendering.VectorFont)
        o = new()

        o.renderer = world.renderer
        o.raster_font = raster_font
        o.vector_font = vector_font
        o.state = []
        o.stack_top = 1
        o.clear_color = Orange()
        o.draw_color = White()
        o.width = world.window_width
        o.height = world.window_height

        o.current = AffineTransform{Float64}()
        o.post = AffineTransform{Float64}()

        o.v1 = Geometry.Point{Float64}()
        o.v2 = Geometry.Point{Float64}()

        o
    end
end

# const SAMPLESIZE = 16

function initialize(context::RenderContext, world::Ranger.World)
    # SDL2.GL_SetAttribute(SDL2.GL_MULTISAMPLEBUFFERS, SAMPLESIZE)
    # SDL2.GL_SetAttribute(SDL2.GL_MULTISAMPLESAMPLES, SAMPLESIZE)
    # SDL2.SetHint(SDL2.HINT_RENDER_SCALE_QUALITY, "1")

    SDL2.SetRenderDrawBlendMode(context.renderer, SDL2.BLENDMODE_BLEND);

    state = State()

    for _ in 0:STATE_STACK_DEPTH
        push!(context.state, copy(state))
    end

    # Apply centered view-space matrix
    apply!(context, world.view_space);
end

function apply!(context::RenderContext, aft::AffineTransform)
    # Concat this transform onto the current transform but don't push it.
    # Use post multiply
    Math.multiply!(aft, context.current, context.post)
    Math.set!(context.current, context.post)
end

function pre(context::RenderContext)
    c = context.clear_color
    SDL2.SetRenderDrawColor(context.renderer, c.r, c.g, c.b, c.a)
    SDL2.RenderClear(context.renderer);

    # Draw checkered board as an clear indicator for debugging
    # NOTE: disable this code for release builds
    # draw_checkerboard(context);
end

# Push the current state onto the stack
function save!(context::RenderContext)
    top = context.state[context.stack_top]

    top.clear_color = context.clear_color
    top.draw_color = context.draw_color
    Math.set!(top.current, context.current)

    context.stack_top += 1;

    # debugging
    # rect = SDL2.Rect(0, 0, 1, 1)
    # rect.x = 100
    # rect.y = 100
    # rect.w = 400
    # rect.h = 400

    # SDL2.SetRenderDrawColor(context.renderer, 255, 127, 0, 255)
    # draw_text(context, 100, 100, "Ranger", 3, 2, false)
    # draw_filled_rectangle(context, rect)
    # SDL2.SetRenderDrawColor(context.renderer, 255, 255, 255, 255)
    # draw_outlined_rectangle(context, rect)
end

# Pop the current state
function restore!(context::RenderContext)
    context.stack_top -= 1

    top = context.state[context.stack_top]

    context.clear_color = top.clear_color
    context.draw_color = top.draw_color
    Math.set!(context.current, top.current)

    c = context.clear_color
    SDL2.SetRenderDrawColor(context.renderer, c.r, c.g, c.b, c.a)
end

function post(context::RenderContext)
    SDL2.RenderPresent(context.renderer);
end

# Transform a Point
function transform!(context::RenderContext, p::Geometry.Point{Float64}, out::Geometry.Point{Float64})
    Math.transform!(context.current, p, out)
end

# Transform a line or rectangle corners
function transform!(context::RenderContext, p1::Geometry.Point{Float64}, p2::Geometry.Point{Float64}, out1::Geometry.Point{Float64}, out2::Geometry.Point{Float64})
    Math.transform!(context.current, p1, out1)
    Math.transform!(context.current, p2, out2)
end

function transform!(context::RenderContext, vertices::Array{Geometry.Point{Float64},1}, bucket::Array{Point{Float64},1})
    for (idx, vertex) in enumerate(vertices)
        Math.transform!(context.current, vertex, bucket[idx])
    end
end

function transform!(context::RenderContext, mesh::Geometry.Mesh)
    for (idx, vertex) in enumerate(mesh.vertices)
        Math.transform!(context.current, vertex, mesh.bucket[idx])
    end
end

function transform!(context::RenderContext, polygon::Math.Polygon)
    for (idx, vertex) in enumerate(polygon.mesh.vertices)
        Math.transform!(context.current, vertex, polygon.mesh.bucket[idx])
    end
end

function transform!(context::RenderContext, line::Math.Line)
    Math.transform!(context.current, line.p1, line.bucket[1])
    Math.transform!(context.current, line.p2, line.bucket[2])
end

# ,--.,--.,--.,--.,--.,--.,--.,--.,--.,--.,--.,--.,--.,--.,--.,--.,--.
# Draw functions render directly to the device.
# ,__.,__.,__.,__.,__.,__.,__.,__.,__.,__.,__.,__.,__.,__.,__.,__.,__.
function set_draw_color(context::RenderContext, color::Palette)
    context.draw_color = color
    SDL2.SetRenderDrawColor(context.renderer, color.r, color.g, color.b, color.a)
end

function draw_point(context::RenderContext, x::Int32, y::Int32)
    SDL2.RenderDrawPoint(context.renderer, x, y);
end

function draw_line(context::RenderContext, x1::Float64, y1::Float64, x2::Float64, y2::Float64)
    SDL2.RenderDrawLine(context.renderer, Int32(round(x1)), Int32(round(y1)), Int32(round(x2)), Int32(round(y2)));
end

function draw_outlined_rectangle(context::RenderContext, rect::SDL2.Rect)
    SDL2.RenderDrawRect(context.renderer, pointer_from_objref(rect));
end

function draw_filled_rectangle(context::RenderContext, rect::SDL2.Rect)
    SDL2.RenderFillRect(context.renderer, pointer_from_objref(rect));
end

rect_draw = SDL2.Rect(0, 0, 1, 1)

function draw_outlined_rectangle(context::RenderContext, minx::Int32, miny::Int32, maxx::Int32, maxy::Int32)
    rect_draw.x = minx
    rect_draw.y = miny
    rect_draw.w = maxx - minx
    rect_draw.h = maxy - miny
    SDL2.RenderDrawRect(context.renderer, pointer_from_objref(rect_draw));
end

function draw_filled_rectangle(context::RenderContext, minx::Int32, miny::Int32, maxx::Int32, maxy::Int32)
    rect_draw.x = minx
    rect_draw.y = miny
    rect_draw.w = maxx - minx
    rect_draw.h = maxy - miny
    SDL2.RenderFillRect(context.renderer, pointer_from_objref(rect_draw));
end

function draw_horz_line(context::RenderContext, x1::Int32, x2::Int32, y::Int32)
    SDL2.RenderDrawLine(context.renderer, x1, y, x2, y);
end

function draw_horz_line(context::RenderContext, x1::AbstractFloat, x2::AbstractFloat, y::AbstractFloat)
    SDL2.RenderDrawLine(context.renderer, Int32(round(x1)), Int32(round(y)), Int32(round(x2)), Int32(round(y)));
end

function draw_vert_line(context::RenderContext, x::Int32, y1::Int32, y2::Int32)
    SDL2.RenderDrawLine(context.renderer, x, y1, x, y2);
end

function draw_vert_line(context::RenderContext, x::AbstractFloat, y1::AbstractFloat, y2::AbstractFloat)
    SDL2.RenderDrawLine(context.renderer, Int32(round(x)), Int32(round(y1)), Int32(round(x)), Int32(round(y2)));
end

const shifts = [0,1,2,3,4,5,6,7]

# This type of font doesn't support transforms. It will always be axis aligned.
function draw_text(context::RenderContext, x::Integer, y::Integer, text::String, scale::Integer, fill::Integer, invert::Bool)
    cx = Int32(x)
    s = Int32(scale)
    row_width = Int32(get_glyph_width())

    # Is the text colored or the space around it (aka inverted)
    bit_invert = invert ? 0 : 1

    for char in text
        if char == ' '
            cx += row_width * s # move to next column/char/glyph
            continue
        end

        gy = Int32(y) # move y back to the "top" for each char
        glyph = get_glyph(context.raster_font, char)
        for g in glyph
            gx = cx # set to current column
            for shift in shifts # scan each pixel in the glyph
                bit = (g >> shift) & 1

                if bit == bit_invert
                    if scale == 1
                        SDL2.RenderDrawPoint(context.renderer, gx, gy)
                    else
                        fillet = fill
                        if fill > scale
                            fillet = 0
                        end
                        for xl in 0:(scale - fillet)
                            for yl in 0:(scale - fillet)
                                SDL2.RenderDrawPoint(context.renderer, Int32(gx + xl), Int32(gy + yl));
                            end
                        end
                    end
                end
                gx += s
            end
            gy += s # move to next pixel-row in char
        end
        cx += row_width * s # move to next column/char/glyph
    end
end

function draw_checkerboard(context::RenderContext)
    flip = false
    size = 200
    col = 0
    row = 0

    while row < context.height
        while col < context.width
            if flip
                SDL2.SetRenderDrawColor(context.renderer, 100, 100, 100, 255)
            else
                SDL2.SetRenderDrawColor(context.renderer, 80, 80, 80, 255)
            end

            rect_draw.x = col
            rect_draw.y = row
            rect_draw.w = col + size
            rect_draw.h = row + size

            SDL2.RenderFillRect(context.renderer, pointer_from_objref(rect_draw))

            flip = !flip;
            col += size;
        end

        flip = !flip;
        col = 0;
        row += size;
    end
end

# ,--.,--.,--.,--.,--.,--.,--.,--.,--.,--.,--.,--.,--.,--.,--.,--.,--.
# Render functions render based on transformed vertices
# The Render functions use the Draw functions.
# ,__.,__.,__.,__.,__.,__.,__.,__.,__.,__.,__.,__.,__.,__.,__.,__.,__.
function render_line(context::RenderContext, x1::Float64, y1::Float64, x2::Float64, y2::Float64)
    draw_line(context, Int32(round(x1)), Int32(round(y1)), Int32(round(x2)), Int32(round(y2)));
end

function render_lines(context::RenderContext, mesh::Geometry.Mesh)
    first = true

    for v in mesh.bucket
        if first
            context.v1 = v
            first = false
            continue
        else
            context.v2 = v
            first = true
        end

        draw_line(context,
            context.v1.x, context.v1.y,
            context.v2.x, context.v2.y);
    end
end

function render_checkerboard(context::RenderContext, mesh::Geometry.Mesh, oddColor::Palette, evenColor::Palette)
    # render a grid of rectangles defined by min/max points
    flip = false
    vertices = mesh.bucket
    build = true

    # TODO detect negative change in X so that the flip value
    # alternates correctly for even grid sizes. To lazy to fix at the moment.
    for (idx, vertex) in enumerate(mesh.bucket)
        if build
            Geometry.set!(context.v1, vertex)
            build = false
            continue
        else
            Geometry.set!(context.v2, vertex)
            build = true
        end

        if flip
            SDL2.SetRenderDrawColor(context.renderer, oddColor.r, oddColor.g, oddColor.b, oddColor.a)
        else
            SDL2.SetRenderDrawColor(context.renderer, evenColor.r, evenColor.g, evenColor.b, evenColor.a)
        end

        # upper-left
        minx = Int32(round(context.v1.x))
        miny = Int32(round(context.v1.y))

        # bottom-right
        maxx = Int32(round(context.v2.x))
        maxy = Int32(round(context.v2.y))

        rect_draw.x = minx
        rect_draw.y = miny
        rect_draw.w = maxx - minx
        rect_draw.h = maxy - miny

        SDL2.RenderFillRect(context.renderer, pointer_from_objref(rect_draw))

        flip = !flip;
    end
end

# Render an axis aligned rectangle. Rotating any of the vertices
# will cause strange rendering behaviours
function render_aa_rectangle(context::RenderContext, min::Geometry.Point{Float64}, max::Geometry.Point{Float64}, fillStyle::RenderStyle)
    # upper-left
    minx = Int32(round(min.x))
    miny = Int32(round(min.y))

    # bottom-right
    maxx = Int32(round(max.x))
    maxy = Int32(round(max.y))

    if fillStyle == FILLED
        draw_filled_rectangle(context, minx, miny, maxx, maxy)
    elseif fillStyle == OUTLINE
        draw_outlined_rectangle(context, minx, miny, maxx, maxy)
    else
        draw_filled_rectangle(context, minx, miny, maxx, maxy)
        draw_outlined_rectangle(context, minx, miny, maxx, maxy)
    end
end

function render_aabb_rectangle(context::RenderContext, aabb::Geometry.AABB)
    # upper-left
    minx = Int32(round(aabb.min.x))
    miny = Int32(round(aabb.min.y))

    # bottom-right
    maxx = Int32(round(aabb.max.x))
    maxy = Int32(round(aabb.max.y))

    draw_outlined_rectangle(context, minx, miny, maxx, maxy)
end

function render_outlined_polygon(context::RenderContext, polygon::Geometry.Polygon, state::PolygonState)
    bucs = polygon.mesh.bucket

    for idx in 1:(length(bucs) - 1)
        SDL2.RenderDrawLine(context.renderer,
            Int32(round(bucs[idx].x)), Int32(round(bucs[idx].y)),
            Int32(round(bucs[idx + 1].x)), Int32(round(bucs[idx + 1].y)))
    end

    if state == CLOSED
        SDL2.RenderDrawLine(context.renderer,
            Int32(round(bucs[length(bucs)].x)), Int32(round(bucs[length(bucs)].y)),
            Int32(round(bucs[1].x)), Int32(round(bucs[1].y)))
    end
end

ext: jl

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
using Test
using StaticArrays, IntervalSets
import ClimaCore.DataLayouts: IJFH
import ClimaCore:
    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry
using StaticArrays, IntervalSets, LinearAlgebra

FT = Float64

radius = FT(1)
domain = Domains.SphereDomain(radius)

scurl = Operators.Curl()
wcurl = Operators.WeakCurl()

@testset "curl" begin
    Ne = 7
    Nq = 6

    mesh = Meshes.EquiangularCubedSphere(domain, Ne)
    grid_topology = Topologies.Topology2D(mesh)

    quad = Spaces.Quadratures.GLL{Nq}()
    space = Spaces.SpectralElementSpace2D(grid_topology, quad)
    coords = Fields.coordinate_field(space)

    u_local = @. Geometry.UVVector(
        sind(coords.long) * sind(coords.lat),
        cosd(coords.long),
    )

    u = Geometry.transform.(Ref(Geometry.Covariant12Axis()), u_local)

    curl_us = scurl.(u)
    Spaces.weighted_dss!(curl_us)

    curl_uw = wcurl.(u)
    Spaces.weighted_dss!(curl_uw)

    curl_exact = @. Geometry.Contravariant3Vector(
        -2 * sind(coords.long) * cosd(coords.lat) / radius,
    )

    @test curl_us ≈ curl_exact rtol = 1e-2
    @test curl_uw ≈ curl_exact rtol = 1e-2

end

"""
    convergence_rate(err, Δh)
Estimate convergence rate given vectors `err` and `Δh`
    err = C Δh^p+ H.O.T
    err_k ≈ C Δh_k^p
    err_k/err_m ≈ Δh_k^p/Δh_m^p
    log(err_k/err_m) ≈ log((Δh_k/Δh_m)^p)
    log(err_k/err_m) ≈ p*log(Δh_k/Δh_m)
    log(err_k/err_m)/log(Δh_k/Δh_m) ≈ p
"""
convergence_rate(err, Δh) =
    [log(err[i] / err[i - 1]) / log(Δh[i] / Δh[i - 1]) for i in 2:length(Δh)]

@testset "convergence tests for the curl operator on the sphere" begin
    Nes = [3, 9, 27]
    Nqs = [4, 6]

    for (Iq, Nq) in enumerate(Nqs)
        err_scurl = Array{FT}(undef, length(Nes))
        err_wcurl = Array{FT}(undef, length(Nes))
        Δh = Array{FT}(undef, length(Nes))

        for (Ie, Ne) in enumerate(Nes)
            mesh = Meshes.EquiangularCubedSphere(domain, Ne)
            grid_topology = Topologies.Topology2D(mesh)

            quad = Spaces.Quadratures.GLL{Nq}()
            space = Spaces.SpectralElementSpace2D(grid_topology, quad)
            coords = Fields.coordinate_field(space)

            u_local = @. Geometry.UVVector(
                sind(coords.long) * sind(coords.lat),
                cosd(coords.long),
            )

            u = Geometry.transform.(Ref(Geometry.Covariant12Axis()), u_local)

            curl_us = scurl.(u)
            Spaces.weighted_dss!(curl_us)

            curl_uw = wcurl.(u)
            Spaces.weighted_dss!(curl_uw)

            curl_exact = @. Geometry.Contravariant3Vector(
                -2 * sind(coords.long) * cosd(coords.lat) / radius,
            )

            err_scurl[Ie] = norm(curl_us .- curl_exact)
            err_wcurl[Ie] = norm(curl_uw .- curl_exact)
            Δh[Ie] = 1 / Ne

        end

        convergence_rate_scurl = convergence_rate(err_scurl, Δh)
        convergence_rate_wcurl = convergence_rate(err_scurl, Δh)
        for Ie in range(1, length = length(Nes) - 1)
            @test convergence_rate_scurl[Ie] ≈ (Nq - 1) atol = 0.1
            @test convergence_rate_wcurl[Ie] ≈ (Nq - 1) atol = 0.1
        end
    end
end


text:
mutable struct TCPStreamEndpoint
    addr::Sockets.InetAddr{IPv4}
    seq::Union{Nothing, UInt32}
    ack::Union{Nothing, UInt32}
    fin_sent::Bool
end

function tcpstreamid(e1::Sockets.InetAddr{IPv4}, e2::Sockets.InetAddr{IPv4})
    h = zero(UInt)
    if e1.host < e2.host
        h = hash(Int(e1.host), h)
        h = hash(Int(e2.host), h)
    else
        h = hash(Int(e2.host), h)
        h = hash(Int(e1.host), h)
    end
    if e1.port < e2.port
        h = hash(e1.port, h)
        h = hash(e2.port, h)
    else
        h = hash(e2.port, h)
        h = hash(e1.port, h)
    end
    h
end

@enum TCPEvent begin
    TCPE_OPEN           # -> On SYN
    TCPE_CONNECTED      # -> After SYN & SYN-ACK & ACK
    TCPE_CLOSE          # -> After FIN
    TCPE_DISCONNECTED   # -> After RST | FIN & ACK & FIN & ACK
    TCPE_ERROR
    TCPE_RETRANSMIT
end

# TODO deal with ip fragmentation
const TCPHandler = FunctionWrapper{Nothing, Tuple{Sockets.InetAddr{IPv4}, Sockets.InetAddr{IPv4}, Union{TCPEvent, TCPPacket}}}
mutable struct TCPStream
    endpoint1::TCPStreamEndpoint
    endpoint2::TCPStreamEndpoint
    handler::TCPHandler

    function TCPStream(handler, src::Sockets.InetAddr{IPv4}, dst::Sockets.InetAddr{IPv4})
        new(
            TCPStreamEndpoint(src, nothing, nothing, false),
            TCPStreamEndpoint(dst, nothing, nothing, false),
            TCPHandler(handler),
        )
    end
end

function TCPStream(handler, iheader::IPv4Header, theader::TCPHeader)
    TCPStream(
        handler,
        Sockets.InetAddr(iheader.src_ip, theader.src_port),
        Sockets.InetAddr(iheader.dst_ip, theader.dst_port),
    )
end

function get_endpoint_(stream::TCPStream, ip::IPv4, port::UInt16)
    addr = Sockets.InetAddr(ip, port)
    stream.endpoint1.addr == addr && return stream.endpoint1
    stream.endpoint2.addr == addr && return stream.endpoint2
    error("Unknown endpoint $ip $port")
end

function Base.push!(stream::TCPStream, iheader::IPv4Header, packet::TCPPacket)
    theader = packet.header

    src_endpoint = get_endpoint_(stream, iheader.src_ip, theader.src_port)
    dst_endpoint = get_endpoint_(stream, iheader.dst_ip, theader.dst_port)

    if theader.RST
        @assert !theader.SYN
        stream.handler(src_endpoint.addr, dst_endpoint.addr, TCPE_DISCONNECTED)
        return
    end

    if theader.SYN
        if src_endpoint.seq !== nothing
            @error "Seq set before SYN" src_endpoint dst_endpoint iheader packet
            stream.handler(src_endpoint.addr, dst_endpoint.addr, TCPE_ERROR)
            return
        end
        src_endpoint.seq = theader.seq_num + 1
        if theader.ACK
            src_endpoint.ack !== nothing && @error "Ack set before SYN" src_endpoint dst_endpoint iheader packet
            if dst_endpoint.seq === nothing
                @error "SYN-ACK was not preceded by SYN" src_endpoint dst_endpoint iheader packet
                stream.handler(src_endpoint.addr, dst_endpoint.addr, TCPE_ERROR)
                return
            elseif dst_endpoint.seq != theader.ack_num
                @error "SYN-ACK with unexpected ack_num" src_endpoint dst_endpoint iheader packet
                stream.handler(src_endpoint.addr, dst_endpoint.addr, TCPE_ERROR)
                return
            else
                @debug "SYN-ACK received" src_endpoint dst_endpoint iheader packet
            end
            src_endpoint.ack = theader.ack_num
            stream.handler(src_endpoint.addr, dst_endpoint.addr, TCPE_CONNECTED)
        else
            @debug "SYN received" src_endpoint dst_endpoint iheader packet
            stream.handler(src_endpoint.addr, dst_endpoint.addr, TCPE_OPEN)
        end
    elseif theader.FIN
        if !isempty(packet.payload)
            @error "TCP FIN packet with non-empty payload" src_endpoint dst_endpoint iheader packet
            stream.handler(src_endpoint.addr, dst_endpoint.addr, TCPE_ERROR)
            return
        end
        if src_endpoint.fin_sent
            @error "TCP double FIN" src_endpoint dst_endpoint iheader packet
            stream.handler(src_endpoint.addr, dst_endpoint.addr, TCPE_ERROR)
            return
        end

        src_endpoint.fin_sent = true
        src_endpoint.ack = theader.ack_num
        src_endpoint.seq = theader.seq_num + 1
        if dst_endpoint.fin_sent
            stream.handler(src_endpoint.addr, dst_endpoint.addr, TCPE_DISCONNECTED)
        else
            stream.handler(src_endpoint.addr, dst_endpoint.addr, TCPE_CLOSE)
        end
    else
        @assert theader.ACK
        if src_endpoint.seq === nothing || dst_endpoint.seq === nothing || dst_endpoint.ack === nothing
            @error "TCP stream init missed" src_endpoint dst_endpoint iheader packet
            # NOTE: We've seen this connection mid-way, we don't support this at this point
            stream.handler(src_endpoint.addr, dst_endpoint.addr, TCPE_ERROR)
            return
        end

        if theader.seq_num > src_endpoint.seq
            @error "TCP stream seq gap" src_endpoint dst_endpoint iheader packet
            stream.handler(src_endpoint.addr, dst_endpoint.addr, TCPE_ERROR)
            return
        elseif theader.seq_num < src_endpoint.seq
            @warn "TCP retransmission" src_endpoint dst_endpoint iheader packet
            stream.handler(src_endpoint.addr, dst_endpoint.addr, TCPE_RETRANSMIT)
            if theader.seq_num + length(packet.payload) != src_endpoint.seq
                @error "TCP retransmission with different data length" src_endpoint dst_endpoint iheader packet
                stream.handler(src_endpoint.addr, dst_endpoint.addr, TCPE_ERROR)
            end
            return
        end

        if src_endpoint.ack === nothing
            @debug "TCP connection initialised" src_endpoint dst_endpoint
        end
        src_endpoint.ack = theader.ack_num
        src_endpoint.seq = theader.seq_num + length(packet.payload)

        if !isempty(packet.payload)
            stream.handler(src_endpoint.addr, dst_endpoint.addr, packet)
        end
    end
    nothing
end

max_stars_repo_name: ancapdev/NetworkProtocols

max_issues_repo_name: ancapdev/NetworkProtocols

max_forks_repo_name: ancapdev/NetworkProtocols


text:
module JLLUtils

include("exceptions.jl");
include("units.jl");
include("stats.jl");
include("arrays.jl");
include("datemath.jl");
include("plotting.jl");

import .JLLUnits: parse_unit_string, strip_units;
import .JLLStats: normalize_to, normalize_to!;
import .JLLPlots: plotas, plotas!;

export parse_unit_string, strip_units, normalize_to, normalize_to!, plotas, plotas!;

end # module

max_stars_repo_name: joshua-laughner/Julia-JLLUtils

max_issues_repo_name: joshua-laughner/Julia-JLLUtils

max_forks_repo_name: joshua-laughner/Julia-JLLUtils


text:
function asap_colmem(num_links, A, dim::Tuple{Int64,Int64}, mem_cols, vert_spacing)
    arch = TopLevel{A,2}("asap_colmem")

    # Extra parameters
    num_fifos = 2
    row = dim[1]
    col = dim[2]

    # generate memory locations according to num of mem cols and vertical spacing
    hori_spacing = floor((col-2)/(mem_cols))
    mem_rows = floor(row/(vert_spacing+2))
    memories = Array{MemoryLocation{2},1}()
    offsets = [(0,-1),(0,1)]
    prev_col = 0
    for r = 1:mem_rows
        for c = 1:mem_cols
            if (prev_col == (Int64((hori_spacing*c+1))-2)) && (prev_col != 0)
                error("$mem_cols memory columns cannot fit inside row:
                       $row and column: $col architecture.")
            end
            address = (Int64((vert_spacing+2)*r), Int64((hori_spacing*c+2)))
            other_addresses = (Int64(((vert_spacing+2)*r)+1),
                               Int64((hori_spacing*c+2)))
            m = MemoryLocation(address, other_addresses, offsets)
            push!(memories,m)
            prev_col = Int64((hori_spacing*c+1))
        end
    end

    mem_addrs = mem_addresses(memories)
    mems = mem_locations(memories)
    mem_neighbor_addrs = mem_neighbors(memories)

    ####################
    # Normal Processor #
    ####################
    processor = build_processor_tile(num_links)
    for r in 2:row+1, c in 2:col+1
        addr = Address(r,c)
        # Skip if this is an address occupied by a memory or memory neighbor.
        if in(addr, mem_addrs) || in(addr, mem_neighbor_addrs)
            continue
        end
        add_child(arch, processor, addr)
    end

    ####################
    # Memory Processor #
    ####################
    # Instantiate a memory processor at every address neighboring a memory.
    memory_processor = build_processor_tile(num_links, include_memory = true)
    for mem_neighbor_addr in mem_neighbor_addrs
        add_child(arch, memory_processor, mem_neighbor_addr)
    end

    #################
    # 2 Port Memory #
    #################
    memory = build_memory(2)
    for mem_addr in mems
        add_child(arch, memory, mem_addr)
    end

    handler_spacing = floor(row/4)
    handler_row = (Int64(handler_spacing), Int64(handler_spacing*2),
                   Int64(handler_spacing*3), Int64(handler_spacing*4))
    #################
    # Input Handler #
    #################
    input_handler = build_input_handler(1)
    for (r,c) ∈ zip(handler_row, (1, col+2, 1, col+2))
        add_child(arch, input_handler, Address(r,c))
    end

    ##################
    # Output Handler #
    ##################
    output_handler = build_output_handler(1)
    for (r,c) ∈ zip(handler_row, (col+2, 1, col+2, 1))
        add_child(arch, output_handler, Address(r,c))
    end

    connect_processors(arch, num_links)
    connect_io(arch, num_links)
    connect_memories_colmem(arch, memories)
    connect_processors_colmem(arch, num_links)

    return arch
end

function connect_memories_colmem(arch, memories)
    # Create metadata dictionary for the memory links.
    request_metadata = Dict(
        "cost"          => 1.0,
        "capacity"      => 1,
        "link_class"    => "memory_request_link",
   )

    response_metadata = Dict(
        "cost"          => 1.0,
        "capacity"      => 1,
        "link_class"    => "memory_response_link",
   )

    ###########################
    # Connect 1 port memories #
    ###########################
    proc_rule = x -> search_metadata!(x, "attributes", "memory_processor", in)
    mem_rule = x -> search_metadata!(x, "attributes", "memory_2port", in)

    for mem in memories
        address = mem.address
        for (i, offset) in enumerate(mem.offsets)
            # Build Request Link.
            offset_rule = [(-offset, "memory_out", "in[$(i-1)]")]
            connection_rule(arch, offset_rule, proc_rule, mem_rule,
                            metadata = request_metadata,
                            valid_addresses = (address + offset,))

            # Build response link.
            offset_rule = [(offset, "out[$(i-1)]", "memory_in")]
            connection_rule(arch, offset_rule, mem_rule, proc_rule,
                            metadata = response_metadata,
                            valid_addresses = (address,))
        end
    end

    return nothing
end

function connect_processors_colmem(tl, num_links)
    fn = x -> search_metadata!(x, "attributes", "processor", in)
    src_rule = fn
    dst_rule = fn

    # Build metadata dictionary for capacity and cost
    metadata = Dict(
        "cost"          => 1.0,
        "capacity"      => 1,
        "link_class"    => "circuit_link"
    )

    offsets = [CartesianIndex(0,2), CartesianIndex(0,-2)]
    src_dirs = ("east", "west")
    dst_dirs = ("west", "east")
    src_ports = [["$(src)_out[$i]" for i in 0:num_links-1] for src in src_dirs]
    dst_ports = [["$(dst)_in[$i]" for i in 0:num_links-1] for dst in dst_dirs]

    offset_rules = []
    for (o,s,d) in zip(offsets, src_ports, dst_ports)
        rules = [(o,i,j) for (i,j) in zip(s,d)]
        append!(offset_rules, rules)
    end

    connection_rule(tl, offset_rules, src_rule, dst_rule, metadata = metadata)
end


text:
module settings

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
#
# spg2.jl --
#
# Implements Spectral Projected Gradient Method (Version 2: "continuous
# projected gradient direction") to find the local minimizers of a given
# function with convex constraints, described in:
#
# [1] E. G. Birgin, J. M. Martinez, and M. Raydan, "Nonmonotone spectral
#     projected gradient methods on convex sets", SIAM Journal on Optimization
#     10, pp. 1196-1211 (2000).
#
# [2] E. G. Birgin, J. M. Martinez, and M. Raydan, "SPG: software for
#     convex-constrained optimization", ACM Transactions on Mathematical
#     Software (TOMS) 27, pp. 340-349 (2001).
#
# ----------------------------------------------------------------------------
#
# This file is part of OptimPack.jl which is licensed under the MIT
# "Expat" License:
#
# Copyright (C) 2014: Éric Thiébaut.
#
# ----------------------------------------------------------------------------

type SpgResult{T,N}
    x::DenseArray{T,N}
    xbest::DenseArray{T,N}
    f::Float64
    fbest::Float64
    pgtwon::Float64
    pginfn::Float64
    fcnt::Int
    pcnt::Int
    iter::Int
    status::Symbol
    SpgResult(x::DenseArray{T,N}, xbest::DenseArray{T,N}) = new(x, xbest, 0.0, 0.0, 0.0, 0.0, 0, 0, 0, :WORK_IN_PROGRESS)
end

function spg2{T,N}(fg!::Function, prj!::Function, x0::DenseArray{T,N}, m::Integer;
                   maxit::Integer=-1, maxfc::Integer=-1,
                   eps1::Real=1e-6, eps2::Real=1e-6,
                   eta::Real=1.0, ftol::Real=1e-4,
                   lmin::Real=1e-30, lmax::Real=1e+30,
                   amin::Real=0.1, amax::Real=0.9,
                   printer=nothing, verb::Bool=false)

    # Allocate workspace.
    x0 = copy(x0)
    dims = size(x0)
    space = DenseVectorSpace(T, dims)
    d = Array(T, dims)
    g = Array(T, dims)
    x = Array(T, dims)
    g = Array(T, dims)
    g0 = Array(T, dims)
    xbest = Array(T, dims)
    _x = wrap(space, x)
    _g = wrap(space, g)
    _d = wrap(space, d)
    _x0 = wrap(space, x0)
    _g0 = wrap(space, g0)
    _xbest = wrap(space, xbest)

    # Initialization.
    ws = SpgResult(x, xbest)
    local sty, sts, f0, f
    if m > 1
        lastfv = Array(T, m)
        fill!(lastfv, inf(T))
    else
        lastfv = nothing
    end

    # Project initial guess.
    prj!(x, x0)
    ws.pcnt += 1

    # Evaluate function and gradient.
    f = fg!(x, g)
    ws.fcnt += 1

    # Initialize best solution and best function value.
    ws.fbest = f
    copy!(_xbest, _x)

    # Main loop.
    while true

        # Compute continuous projected gradient (and its norms).
        axpby!(_d, 1.0, _x, eta, _g)
        prj!(_d, _d)
        axpby!(_d, 1.0/eta, _x, -1.0/eta, _d)
        ws.pcnt += 1
        ws.pgtwon = norm2(_d)
        ws.pginfn = norminf(_d)

        # Print iteration information
        if printer <: Function
            printer(ws)
        end
        if verb
            @printf("ITER = %-5d  EVAL = %-5d  PROJ = %-5d  F(%s) =%24.17e  ||PG||oo = %17.10e\n",
                    ws.iter, ws.fcnt, ws.pcnt, (f <= ws.fbest ? "+" : "-"), f, ws.pginfn)
        end

        # Test stopping criteria.
        if ws.pginfn <= eps1
            # Gradient infinite-norm stopping criterion satisfied, stop.
            ws.status = :CONVERGENCE_WITH_INFNORM
            return ws
        end
        if ws.pgtwon <= eps2
            # Gradient 2-norm stopping criterion satisfied, stop.
            ws.status = :CONVERGENCE_WITH_TWONORM
            return ws
        end
        if maxit >= 0 && ws.iter >= maxit
            # Maximum number of iterations exceeded, stop.
            ws.status = :TOO_MANY_ITERATIONS
            return ws
        end
        if maxfc >= 0 && ws.fcnt >= maxfc
            # Maximum number of function evaluations exceeded, stop.
            ws.status = :TOO_MANY_EVALUATIONS
            return ws
        end

        # Store function value for the nonmonotone line search and
        # find maximum function value since m last calls.
        if m > 1
            lastfv[(ws.iter%m) + 1] = f
            fmax = maximum(lastfv)
        else
            fmax = f
        end

        # Compute spectral steplength.
        if iter == 0
            # Initial steplength.
            lambda = min(lmax, max(lmin, 1.0/ws.pginfn))
        else
            _s = _x0 # alias
            _y = _g0 # alias
            axpby!(_s, 1.0, _x, -1.0, _x0)
            axpby!(_y, 1.0, _g, -1.0, _g0)
            sty = dot(_s, _y)
            if sty > 0.0
                # Safeguarded Barzilai & Borwein spectral steplength.
                sts = dot(_s, _s)
                lambda = min(lmax, max(lmin, sts/sty))
            else
                lambda = lmax
            end
        end

        # Save current point.
        copy!(_x0, _x)
        copy!(_g0, _g)
        f0 = f

        # Compute the spectral projected gradient direction and <G,D>
        axpby!(_x, 1.0, _x, -lambda, _g)
        prj(_x, _x)
        ws.pcnt += 1
        axpby!(_d, 1.0, _x, -1.0, _x0)
        delta = dot(_g0, _d)

        # Nonmonotone line search.
        stp = 1.0 # Step length for first trial.
        while true
            # Evaluate function and gradient at trial point.
            f = fg!(x, g)
            ws.fcnt += 1

            # Compare the new function value against the best function
            # value and, if smaller, update the best function value and the
            # corresponding best point.
            if f < ws.fbest
                ws.fbest = f
                copy!(_xbest, _x)
            end

            # Test stopping criteria.
            if f <= fmax + stp*ftol*delta
                # Nonmonotone Armijo-like stopping criterion satisfied, stop.
                break
            end
            if maxfc >= 0 && ws.fcnt >= maxfc
                # Maximum number of function evaluations exceeded, stop.
                ws.status = :TOO_MANY_EVALUATIONS
                return ws
            end

            # Safeguarded quadratic interpolation.
            q = -delta*(stp*stp)
            r = (f - f0 - stp*delta)*2.0
            if r > 0.0 && amin*r <= q <= amax*stp*r
                stp = q/r
            else
                stp /= 2.0
            end

            # Compute trial point.
            axpby!(_x, 1.0, _x0, stp, _d)
        end

        if ws.status != :WORK_IN_PROGRESS
            # The number of function evaluations was exceeded inside the line
            # search.
            return ws
        end

        # Proceed with next iteration.
        ws.iter += 1

    end
end


text:
export status_t, status_exists, status_attach, status_detach, status_lock, status_unlock,
        status_clear, status_chkinit, update_status, status_buf_lock_unlock

# Status constants
"""
    const global STATUS_TOTAL_SIZE = 184320

The maximum size of a HASHPIPE status buffer - 184,320 bytes.
"""
const global STATUS_TOTAL_SIZE = 184320 # 2880 * 64
"""
    const global STATUS_RECORD_SIZE = 80

The size of a status record - 80 characters.
"""
const global STATUS_RECORD_SIZE = 80

"""
    Hashpipe Status struct

Data representing a Hashpipe status buffer.

Need to create empty status struct before trying to attaching
to existing status buffer.
Example:
'''
    instance_id = 0
    status = status_t(0,0,0,0)
    status_attach(instance_id, Ref(status))
'''
"""
mutable struct status_t
    "Instance ID of this status buffer. DO NOT SET/CHANGE!"
    instance_id::Cint
    "Shared memory segment id."
    shmid::Cint
    "POSIX semaphore descriptor for locking "
    p_lock::Ptr{UInt8}
    "Pointer to data area."
    p_buf::Ptr{UInt8}
end

"""
    status_exists(instance_id)

Check whether or not the Hashpipe status buffer exists for the given Hashpipe instance.

Returns non-zero if the status buffer for instance already exists.
"""
function status_exists(instance_id::Int)
    exists::Int8 = ccall((:hashpipe_status_exists,
                libhashpipestatus),
                Int8, (Int8,), instance_id)
    return exists
end

"""
    status_attach(instance_id::Int, p_status::Ref{status_t})

Populate the Hashpipe status pointed to by p_status with the status values of the Hashpipe
instance given (created if doesn't already exist).

Attach/create lock semaphore as well.  Return nonzero on error.
"""
function status_attach(instance_id::Int, p_status::Ref{status_t})
    error::Int8 = ccall((:hashpipe_status_attach, libhashpipestatus),
                    Int, (Int, Ref{status_t}), instance_id, p_status)
    return error
end

function status_detach(p_status::Ref{status_t})
    error::Int8 = ccall((:hashpipe_status_detach, libhashpipestatus),
                    Int, (Ref{status_t}, ), p_status)
    return error
end

"""
    status_lock(p_status::Ref{status_t})

Lock the status pointed to by p_status (probably for updating status values).

If locked, will sleep while waiting for the buffer to become unlocked.
"""
function status_lock(p_status::Ref{status_t})
    error::Int8 = ccall((:hashpipe_status_lock, libhashpipestatus),
                    Int, (Ref{status_t},), p_status)
    return error
end

"""
    status_unlock(p_status::Ref{status_t})

Unlock the status pointed to by p_status (probably after updating status values).

If unlocked, will sleep while waiting for the buffer to become locked.
"""
function status_unlock(p_status::Ref{status_t})
    error::Int8 = ccall((:hashpipe_status_unlock, libhashpipestatus),
                    Int, (Ref{status_t},), p_status)
    return error
end

"""
    status_buf_lock_unlock(f::Function, r_status::Ref{status_t})

Safely lock and unlock a shared status buffer for updating its values. This must be done
so that the status buffer values aren't changed by multiple processes at the same time.

Example:

"""
function status_buf_lock_unlock(f::Function, r_status::Ref{status_t})
        try
            status_lock(r_status)
            f() # or f(st) TODO: test which of these is better
        catch e
            @error "Error locking hashpipe status buffer - Error: $e"
        finally
            status_unlock(r_status)
        end
end

"""
    status_clear(p_status::Ref{status_t})

Clear the status values of the status buffer shared memory.
"""
function status_clear(p_status::Ref{status_t})
    ccall((:hashpipe_status_clear, libhashpipestatus),
            Int, (Ref{status_t},), p_status)
    return nothing
end

"""
    status_chkinit(p_status::Ref{status_t})

Check the status buffer for appropriate formatting (existence of "END").
If not found, zero it out and add END.
"""
function status_chkinit(p_status::Ref{status_t})
    ccall((:hashpipe_status_chkinit, libhashpipestatus),
            Int, (Ref{status_t},), p_status)
    return nothing
end

"""
    Base.display(s::status_t)

Display a Hashpipe status nicely in REPL.
"""
function Base.display(s::status_t)
    # Calculate max number of header records
    BUFFER_MAX_RECORDS = Int(STATUS_TOTAL_SIZE / STATUS_RECORD_SIZE)
    println("Instance ID: $(s.instance_id)")
    println("shmid: $(s.shmid)")

    # Check to see if valid pointer
    if s.p_lock != C_NULL
        lock = unsafe_wrap(Array, s.p_lock, (1))[1]
        println("Lock: $lock")
    else
        println("Lock: NULL")
    end

    # Check to see if valid pointer
    if s.p_buf != C_NULL
        println("Buffer:")
        string_array = unsafe_wrap(Array, s.p_buf, (STATUS_RECORD_SIZE, BUFFER_MAX_RECORDS))
        # Iterate through records and display them
        for record in 1:size(string_array, 2)
            record_string = String(string_array[:, record])
            println("\t", record_string)
            if record_string[1:3] == "END"
                return nothing
            end
        end
    else
        println("Buffer: NULL")
    end
    return nothing
end

"""
    Base.display(r::Ref{status_t})

Display Hashpipe status from reference nicely in REPL.
"""
function Base.display(r::Ref{status_t})
    display(r[])
    return nothing
end

# Todo Auto-generate all update_status possibilities

function update_status(status::status_t, key::String, value::String)::Int8
    key = Cstring(pointer(key)) # Need to convert for hput functions
    error::Int8 = ccall((:hputs, libhashpipestatus),
                    Int, (Ptr{UInt8}, Cstring, Cstring),
                    status.p_buf, Cstring(pointer(key)), Cstring(pointer(value)))
    return error
end

function update_status(status::status_t, key::String, value::Int)::Int8
    key = Cstring(pointer(key)) # Need to convert for hput functions
    error::Int8 = ccall((:hputi4, libhashpipestatus),
                    Int, (Ptr{UInt8}, Cstring, Cint),
                    status.p_buf, Cstring(pointer(key)), Cint(value))
    return error
end


text:
using Yao, YaoExtensions
using Test, BitBasis

function toffli(b::BitStr)
    t = @inbounds b[1] ⊻ (b[3] & b[2])
    return @inbounds bit_literal(t, b[2], b[3])
end

"""
    pshift(n::Int) -> Function

return a peridoc shift function.
"""
pshift(n::Int) = function (b::BitStr{N}) where N mod(b+n, 1<<N) end

@testset "test toffli" begin
    g = YaoExtensions.mathgate(3, toffli)
    check_truth(b1, b2) = apply!(ArrayReg(b1), g) == ArrayReg(b2)
    @test check_truth(bit"000", bit"000")
    @test check_truth(bit"001", bit"001")
    @test check_truth(bit"010", bit"010")
    @test check_truth(bit"011", bit"011")
    @test check_truth(bit"100", bit"100")
    @test check_truth(bit"101", bit"101")
    @test check_truth(bit"110", bit"111")
    @test check_truth(bit"111", bit"110")
end

@testset "test pshift" begin
    # bint
    nbits = 5
    ab = YaoExtensions.mathgate(nbits, pshift(3))
    mb = YaoExtensions.mathgate(nbits, pshift(-3))
    @test apply!(zero_state(nbits), ab) == product_state(nbits, 3)
    @test apply!(zero_state(nbits), mb) == product_state(nbits, 1<<nbits-3)
    @test applymatrix(ab) ≈ mat(ab)
    @test isunitary(ab)

    ab = YaoExtensions.mathgate(nbits, pshift(3))
    @test apply!(zero_state(nbits), ab) == product_state(nbits, 0b00011)
    @test isunitary(ab)
end

max_stars_repo_name: QuantumBFS/YaoExtensions

max_issues_repo_name: QuantumBFS/YaoExtensions

max_forks_repo_name: QuantumBFS/YaoExtensions


text:
mutable struct RobustMaximinModel <: AbstractModel
    settings::Settings
    constraints::Vector{Constraint}
    obj::RobustMaximinObjective
    output::Output
    RobustMaximinModel(settings, constraints, obj, output) =
        new(settings, constraints, obj, output)
    RobustMaximinModel() = new(Settings(), Constraint[], RobustMaximinObjective(), Output())
end


text:
# Compute the starting x & v positions:

using DelimitedFiles
using JLD2

path = "/Users/ericagol/Software/TRAPPIST1_Spitzer"
# Relevant code:
include(string(path,"/src/regress.jl"))
include(string(path,"/src/NbodyGradient/src/ttv.jl"))

# Load in parameters and results from the likelihood profile:
#grid_file=string(path,"/data/T1_likelihood_profile_student_planetc_3.0sig.jld2")
@load "/Users/ericagol/Software/TRAPPIST1_Spitzer//data/T1_likelihood_profile_student_planetc_3.0sig.jld2"

# Choose the maximum likelihood elements:
elements = readdlm(string(path,"/data/elements_noprior_students.txt"),',')
x,v = init_nbody(elements,t0,n)

println("x: ",x)
println("v: ",v)

writedlm("T1_maxlike_xv.txt",[x;v])

open("T1_maxlike_tt.txt","w") do io
for i=2:8
  for j=1:count1[i]
    println(io,string(i," ",j," ",tt_grid[1,i,j,16]))
  end
end
end


text:
module TestSetMacro

module Clone
using Setfield: setmacro, lensmacro

macro lens(ex)
    lensmacro(identity, ex)
end

macro set(ex)
    setmacro(identity, ex)
end

end#module Clone

using Setfield: Setfield
using Test
using .Clone: Clone

using StaticArrays: @SMatrix

@testset "setmacro, lensmacro isolation" begin

    # test that no symbols like `IndexLens` are needed:
    @test Clone.@lens(_                                   ) isa Setfield.Lens
    @test Clone.@lens(_.a                                 ) isa Setfield.Lens
    @test Clone.@lens(_[1]                                ) isa Setfield.Lens
    @test Clone.@lens(first(_)                            ) isa Setfield.Lens
    @test Clone.@lens(_[end]                              ) isa Setfield.Lens
    @test Clone.@lens(_[$1]                           ) isa Setfield.Lens
    @test Clone.@lens(_.a[1][end, end-2].b[$1, $1]) isa Setfield.Lens

    @test Setfield.@lens(_.a) === Clone.@lens(_.a)
    @test Setfield.@lens(_.a.b) === Clone.@lens(_.a.b)
    @test Setfield.@lens(_.a.b[1,2]) === Clone.@lens(_.a.b[1,2])

    o = (a=1, b=2)
    @test Clone.@set(o.a = 2) === Setfield.@set(o.a = 2)
    @test Clone.@set(o.a += 2) === Setfield.@set(o.a += 2)

    m = @SMatrix [0 0; 0 0]
    m2 = Clone.@set m[end-1, end] = 1
    @test m2 === @SMatrix [0 1; 0 0]
    m3 = Clone.@set(first(m) = 1)
    @test m3 === @SMatrix[1 0; 0 0]
end

end#module


text:
using PyPlot, Revise, HiQGA.transD_GP, Random
## set up
nfreqsperdecade     = 20
ntimesperdecade     = 20
modelprimary        = true
provideddt          = true
doconvramp          = false
nkᵣeval             = 200
times               = 10 .^LinRange(-5,-1.1, 60)
lowpassfcs          = [1e6]
freqhigh = 1e6
freqlow = 1e-8
## model
zfixed   = [-1e5,   0]
rho      = [1e12,   100]
nmax = 200
##  geometry
rRx = 100.
zRx = 0.
zTx = 0.
modelprimary = true
##
F = transD_GP.AEM_VMD_HMD.HFieldDHT(freqlow = freqlow,
                      freqhigh = freqhigh,
                      zTx    = zTx,
                      rRx    = rRx,
                      times = times,
                      nfreqsperdecade = nfreqsperdecade,
                      ntimesperdecade = ntimesperdecade,
                      zRx    = zRx,
                      nkᵣeval = nkᵣeval,
                      modelprimary = modelprimary,
                      provideddt = provideddt,
                      lowpassfcs = lowpassfcs)
transD_GP.AEM_VMD_HMD.getfieldTD!(F, zfixed, rho)
## plot FD
figure()
loglog(F.freqs, abs.(real.(F.HFD_z)))
loglog(F.freqs, abs.(imag.(F.HFD_z)))
xlim(1e-1,1e5)
ylim(1e-12, 1e-6)
xlabel("Hz")
title("Compare with W&H 4.2")
## plot time domain
figure(figsize=(4,7))
# first dhzdt
F.provideddt = true
transD_GP.AEM_VMD_HMD.getfieldTD!(F, zfixed, rho)
loglog(F.interptimes*1e3, abs.(F.HTD_z_interp), label="dhz/dt")
xlim(1e-5,1e3)
ylim(1e-11,1e-1)
grid(true, which="both")
# then dhdt
F.provideddt = false
transD_GP.AEM_VMD_HMD.getfieldTD!(F, zfixed, rho)
loglog(F.interptimes*1e3, abs.(F.HTD_z_interp), label="hz")
legend()
xlabel("time in ms")
title("Compare with W&H Fig 4.4")
plt.tight_layout()
## Radial fields
figure(figsize=(4,7))
F.useprimary = 0.
F.getradialH = true
F.provideddt = true
transD_GP.AEM_VMD_HMD.getfieldTD!(F, zfixed, rho)
loglog(F.interptimes*1e3, abs.(F.HTD_r_interp), label="dhr/dt")
ylim(1e-11,1e-1)
legend()
ylabel("dh/dt")
# grid(true, which="both")
xlabel("time in ms")
F.provideddt = false
transD_GP.AEM_VMD_HMD.getfieldTD!(F, zfixed, rho)
ax = twinx(gca())
ax.loglog(F.interptimes*1e3, abs.(F.HTD_r_interp), label="hr", color="k")
ax.set_ylabel("h")
ax.set_xlim(1e-3,10^3.5)
ax.set_ylim(10^-7.05,10^-16.9)
ax.invert_yaxis()
# ax.grid(true, which="both")
title("Compare with W&H Fig 4.5")
legend(loc="lower left")
plt.tight_layout()
## investigative function
function plotkernelfunc(G::typeof(F); decfactor=10, J0=true)
    figure()
    for ifreq = 1:decfactor:length(G.freqs)
        freq = G.freqs[ifreq]
        if J0
            loglog(G.interpkᵣ, abs.(real(G.J0_kernel_v[:,ifreq])), label="$freq Hz", "--")
            loglog(G.interpkᵣ, abs.(imag(G.J0_kernel_v[:,ifreq])), "-")
        else
            loglog(G.interpkᵣ, abs.(real(G.J1_kernel_v[:,ifreq])), label="$freq Hz", "--")
            loglog(G.interpkᵣ, abs.(imag(G.J1_kernel_v[:,ifreq])), "-")
        end
    end
    J0 ? title("J0 v") : title("J1 v")
    legend()
    figure()
    s1 = subplot(211)
    if J0
        pcolormesh(G.interpkᵣ, G.freqs, log10.(abs.(real(G.J0_kernel_v)))')
    else
        pcolormesh(G.interpkᵣ, G.freqs, log10.(abs.(real(G.J1_kernel_v)))')
    end
    colorbar()
    J0 ? title("real J0 kernel v") : title("real J1 kernel v")
    ylabel("Hz")
    subplot(212, sharex=s1, sharey=s1)
    if J0
        pcolormesh(G.interpkᵣ, G.freqs, log10.(abs.(imag(G.J0_kernel_v)))')
    else
        pcolormesh(G.interpkᵣ, G.freqs, log10.(abs.(imag(G.J1_kernel_v)))')
    end
    colorbar()
    J0 ? title("imag J0 kernel v") : title("imag J1 kernel v")
    gca().set_yscale("log")
    gca().set_xscale("log")
    xlabel("kᵣ 1/m"); ylabel("Hz")
    plt.tight_layout()
end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
T=9
tsname="joOuterProd"
@testset "$tsname" begin
for t=T-2:T # start test loop
    m=t
    n=m+rand(-2:2:2)
    ccA=joOuterProd(rand(joComplex,m,2),rand(joComplex,n,2);DDT=joComplex)
    crA=joOuterProd(rand(joComplex,m,2),rand(joFloat,n,2);DDT=joComplex)
    rcA=joOuterProd(rand(joFloat,m),rand(joComplex,n);DDT=joComplex)
    rrA=joOuterProd(rand(joFloat,m),rand(joFloat,n);DDT=joFloat)

    verbose && println("$tsname ($m x $n)")
    @testset "$m x $n" begin
        @test isadjoint(ccA)[1]
        @test islinear(ccA)[1]
        @test isadjoint(crA)[1]
        @test islinear(crA)[1]
        @test isadjoint(rcA)[1]
        @test islinear(rcA)[1]
        @test isadjoint(rrA)[1]
        @test islinear(rrA)[1]
    end

end # end test loop
end


text:
export ReservoirArrayBuffer

using Random
using ElasticArrays
using MacroTools: @forward

mutable struct ReservoirArrayBuffer{T,N,B<:ElasticArray{T,N},R<:AbstractRNG} <:
               AbstractArray{T,N}
    buffer::B
    n::Int
    capacity::Int
    rng::R
end

ReservoirArrayBuffer{T}(dims::Int...; rng = Random.GLOBAL_RNG) where {T} =
    ReservoirArrayBuffer(ElasticArray{T}(undef, dims[1:end-1]..., 0), 0, dims[end], rng)

@forward ReservoirArrayBuffer.buffer Base.size,
Base.getindex,
Base.length,
Base.sizeof,
Base.IndexStyle

# TODO: rename all push! to append!

function Base.push!(b::ReservoirArrayBuffer{T,N}, x) where {T,N}
    b.n += 1
    if b.n <= b.capacity
        append!(b.buffer, x)
    else
        i = rand(b.rng, 1:b.n)
        if i <= b.capacity
            stride = b.buffer.kernel_length.divisor
            b.buffer.data[(stride*(i-1)+1):stride*i] .= x
        end
    end
end


text:
using Dalinar
using Test

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


text:
using PtFEM, Test

data = Dict(
  # GenericSolid(ndim, nst, nels, nn, nip, finite_element(nod, nodof), axisymmentrix)
  :struc_el => GenericSolid(3, 6, 6, 8, 1, Tetrahedron(4, 3), false),
  :properties => [100.0 0.3 0.0;],
  :g_coord => [
    0.0 0.0 0.0; 1.0 0.0 0.0; 0.0 0.0 -1.0; 1.0 0.0 -1.0;
    0.0 1.0 0.0; 1.0 1.0 0.0; 0.0 1.0 -1.0; 1.0 1.0 -1.0],
  :g_num => [
    1, 3, 4, 7,
    1, 4, 2, 7,
    1, 2, 5, 7,
    6, 4, 8, 7,
    6, 2, 4, 7,
    6, 5, 2, 7
  ],
  :support => [
    (1, [0 0 1]), (2, [1 0 1]), (3, [0 0 0]), (4, [1 0 0]),
    (5, [0 1 1]), (7, [0 1 0]), (8, [1 1 0])
  ],
  :loaded_nodes => [
    ( 1, [0.0 0.0 -1/6]),
    ( 2, [0.0 0.0 -1/3]),
    ( 5, [0.0 0.0 -1/3]),
    ( 6, [0.0 0.0 -1/6])
  ]
)

@time m, dis_df, fm_df = p54(data)

@test m.loads[0:12] ≈ [0.0,-0.010000000000000004,0.003000000000000002,
  -0.010000000000000004,0.003000000000000001,0.003000000000000001,
  -0.010000000000000004,0.0030000000000000014,0.0030000000000000014,
  -0.010000000000000004,0.003000000000000001,0.0030000000000000022,
  0.003000000000000002]


text:
module ScalarTests

using ReverseDiff, ForwardDiff, Test, DiffRules, SpecialFunctions, NaNMath

include(joinpath(dirname(@__FILE__), "../utils.jl"))

x, a, b = rand(3)
tp = InstructionTape()
int_range = 1:10

function test_forward(f, x, tp::InstructionTape, is_domain_err_func::Bool)
    xt = ReverseDiff.TrackedReal(x, zero(x), tp)
    y = f(x)

    # record
    yt = f(xt)
    @test yt == y
    @test length(tp) == 1

    # reverse
    ReverseDiff.seed!(yt)
    ReverseDiff.reverse_pass!(tp)
    @test deriv(xt) == ForwardDiff.derivative(f, x)

    # forward
    x2 = is_domain_err_func ? rand() + 1 : rand()
    ReverseDiff.value!(xt, x2)
    ReverseDiff.forward_pass!(tp)
    @test value(yt) == f(x2)
    ReverseDiff.value!(xt, x)

    empty!(tp)
end

function test_forward(f, a, b, tp)
    at = ReverseDiff.TrackedReal(a, 0.0, tp)
    bt = ReverseDiff.TrackedReal(b, 0.0, tp)
    c = f(a, b)

    ########################################

    # record
    ct = f(at, b)
    @test ct == c
    @test length(tp) == 1

    # reverse
    ReverseDiff.seed!(ct)
    ReverseDiff.reverse_pass!(tp)
    test_approx(deriv(at), ForwardDiff.derivative(x -> f(x, b), a))
    ReverseDiff.unseed!(at)

    # forward
    a2 = isa(a, Int) ? rand(int_range) : rand()
    ReverseDiff.value!(at, a2)
    ReverseDiff.forward_pass!(tp)
    @test value(ct) == f(a2, b)
    ReverseDiff.value!(at, a)

    empty!(tp)

    ########################################

    # record
    ct = f(a, bt)
    @test ct == c
    @test length(tp) == 1

    # reverse
    ReverseDiff.seed!(ct)
    ReverseDiff.reverse_pass!(tp)
    test_approx(deriv(bt), ForwardDiff.derivative(x -> f(a, x), b))
    ReverseDiff.unseed!(bt)

    # forward
    b2 = isa(b, Int) ? rand(int_range) : rand()
    ReverseDiff.value!(bt, b2)
    ReverseDiff.forward_pass!(tp)
    @test value(ct) == f(a, b2)
    ReverseDiff.value!(bt, b)

    empty!(tp)

    ########################################

    # record
    ct = f(at, bt)
    @test ct == c
    @test length(tp) == 1

    # reverse
    ReverseDiff.seed!(ct)
    ReverseDiff.reverse_pass!(tp)
    test_approx(deriv(at), ForwardDiff.derivative(x -> f(x, b), a))
    test_approx(deriv(bt), ForwardDiff.derivative(x -> f(a, x), b))

    # forward
    a2 = isa(a, Int) ? rand(int_range) : rand()
    b2 = isa(b, Int) ? rand(int_range) : rand()
    ReverseDiff.value!(at, a2)
    ReverseDiff.value!(bt, b2)
    ReverseDiff.forward_pass!(tp)
    @test value(ct) == f(a2, b2)
    ReverseDiff.value!(bt, a)
    ReverseDiff.value!(bt, b)

    empty!(tp)
end

function test_skip(f, x, tp)
    xt = ReverseDiff.TrackedReal(x, zero(x), tp)
    y = f(x)
    yt = f(xt)
    @test yt == y
    @test isempty(tp)
end

function test_skip(f, a, b, tp)
    at = ReverseDiff.TrackedReal(a, zero(a), tp)
    bt = ReverseDiff.TrackedReal(b, zero(b), tp)
    c = f(a, b)

    ct = f(at, b)
    @test ct == c
    @test isempty(tp)

    ct = f(a, bt)
    @test ct == c
    @test isempty(tp)

    ct = f(at, bt)
    @test ct == c
    @test isempty(tp)
end

DOMAIN_ERR_FUNCS = (:asec, :acsc, :asecd, :acscd, :acoth, :acosh)

for (M, f, arity) in DiffRules.diffrules()
    f === :rem2pi && continue
    if arity == 1
        test_println("forward-mode unary scalar functions", string(M, ".", f))
        is_domain_err_func = in(f, DOMAIN_ERR_FUNCS)
        n = is_domain_err_func ? x + 1 : x
        test_forward(eval(:($M.$f)), n, tp, is_domain_err_func)
    elseif arity == 2
        in(f, SKIPPED_BINARY_SCALAR_TESTS) && continue
        test_println("forward-mode binary scalar functions", f)
        test_forward(eval(:($M.$f)), a, b, tp)
    end
end

INT_ONLY_FUNCS = (:iseven, :isodd)

for f in ReverseDiff.SKIPPED_UNARY_SCALAR_FUNCS
    test_println("SKIPPED_UNARY_SCALAR_FUNCS", f)
    n = in(f, DOMAIN_ERR_FUNCS) ? x + 1 : x
    n = in(f, INT_ONLY_FUNCS) ? ceil(Int, n) : n
    test_skip(eval(f), n, tp)
end

for f in ReverseDiff.SKIPPED_BINARY_SCALAR_FUNCS
    test_println("SKIPPED_BINARY_SCALAR_FUNCS", f)
    test_skip(eval(f), a, b, tp)
end

end # module


text:
@testset "Interactions" begin
    c1 = SVector(1.0, 1.0, 1.0)u"nm"
    c2 = SVector(1.3, 1.0, 1.0)u"nm"
    c3 = SVector(1.4, 1.0, 1.0)u"nm"
    a1 = Atom(charge=1.0, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1")
    box_size = SVector(2.0, 2.0, 2.0)u"nm"
    dr12 = vector(c1, c2, box_size)
    dr13 = vector(c1, c3, box_size)

    for inter in (LennardJones(), Mie(m=6, n=12))
        @test isapprox(
            force(inter, dr12, c1, c2, a1, a1, box_size),
            SVector(16.0, 0.0, 0.0)u"kJ * mol^-1 * nm^-1",
            atol=1e-9u"kJ * mol^-1 * nm^-1",
        )
        @test isapprox(
            force(inter, dr13, c1, c3, a1, a1, box_size),
            SVector(-1.375509739, 0.0, 0.0)u"kJ * mol^-1 * nm^-1",
            atol=1e-9u"kJ * mol^-1 * nm^-1",
        )
        @test isapprox(
            potential_energy(inter, dr12, c1, c2, a1, a1, box_size),
            0.0u"kJ * mol^-1",
            atol=1e-9u"kJ * mol^-1",
        )
        @test isapprox(
            potential_energy(inter, dr13, c1, c3, a1, a1, box_size),
            -0.1170417309u"kJ * mol^-1",
            atol=1e-9u"kJ * mol^-1",
        )
    end

    inter = SoftSphere()
    @test isapprox(
        force(inter, dr12, c1, c2, a1, a1, box_size),
        SVector(32.0, 0.0, 0.0)u"kJ * mol^-1 * nm^-1",
        atol=1e-9u"kJ * mol^-1 * nm^-1",
    )
    @test isapprox(
        force(inter, dr13, c1, c3, a1, a1, box_size),
        SVector(0.7602324486, 0.0, 0.0)u"kJ * mol^-1 * nm^-1",
        atol=1e-9u"kJ * mol^-1 * nm^-1",
    )
    @test isapprox(
        potential_energy(inter, dr12, c1, c2, a1, a1, box_size),
        0.8u"kJ * mol^-1",
        atol=1e-9u"kJ * mol^-1",
    )
    @test isapprox(
        potential_energy(inter, dr13, c1, c3, a1, a1, box_size),
        0.0253410816u"kJ * mol^-1",
        atol=1e-9u"kJ * mol^-1",
    )

    inter = Coulomb()
    @test isapprox(
        force(inter, dr12, c1, c2, a1, a1, box_size),
        SVector(1543.727311, 0.0, 0.0)u"kJ * mol^-1 * nm^-1",
        atol=1e-5u"kJ * mol^-1 * nm^-1",
    )
    @test isapprox(
        force(inter, dr13, c1, c3, a1, a1, box_size),
        SVector(868.3466125, 0.0, 0.0)u"kJ * mol^-1 * nm^-1",
        atol=1e-5u"kJ * mol^-1 * nm^-1",
    )
    @test isapprox(
        potential_energy(inter, dr12, c1, c2, a1, a1, box_size),
        463.1181933u"kJ * mol^-1",
        atol=1e-5u"kJ * mol^-1",
    )
    @test isapprox(
        potential_energy(inter, dr13, c1, c3, a1, a1, box_size),
        347.338645u"kJ * mol^-1",
        atol=1e-5u"kJ * mol^-1",
    )

    c1_grav = SVector(1.0, 1.0, 1.0)u"m"
    c2_grav = SVector(6.0, 1.0, 1.0)u"m"
    a1_grav, a2_grav = Atom(mass=1e6u"kg"), Atom(mass=1e5u"kg")
    box_size_grav = SVector(20.0, 20.0, 20.0)u"m"
    dr12_grav = vector(c1_grav, c2_grav, box_size_grav)
    inter = Gravity()
    @test isapprox(
        force(inter, dr12_grav, c1_grav, c2_grav, a1_grav, a2_grav, box_size_grav),
        SVector(-0.266972, 0.0, 0.0)u"kg * m * s^-2",
        atol=1e-9u"kg * m * s^-2",
    )
    @test isapprox(
        potential_energy(inter, dr12_grav, c1_grav, c2_grav,
                         a1_grav, a2_grav, box_size_grav),
        -1.33486u"kg * m^2 * s^-2",
        atol=1e-9u"kg * m^2 * s^-2",
    )

    b1 = HarmonicBond(b0=0.2u"nm", kb=300_000.0u"kJ * mol^-1 * nm^-2")
    b2 = HarmonicBond(b0=0.6u"nm", kb=100_000.0u"kJ * mol^-1 * nm^-2")
    fs = force(b1, c1, c2, box_size)
    @test isapprox(
        fs.f1,
        SVector(30000.0, 0.0, 0.0)u"kJ * mol^-1 * nm^-1",
        atol=1e-9u"kJ * mol^-1 * nm^-1",
    )
    @test isapprox(
        fs.f2,
        SVector(-30000.0, 0.0, 0.0)u"kJ * mol^-1 * nm^-1",
        atol=1e-9u"kJ * mol^-1 * nm^-1",
    )
    fs = force(b2, c1, c3, box_size)
    @test isapprox(
        fs.f1,
        SVector(-20000.0, 0.0, 0.0)u"kJ * mol^-1 * nm^-1",
        atol=1e-9u"kJ * mol^-1 * nm^-1",
    )
    @test isapprox(
        fs.f2,
        SVector(20000.0, 0.0, 0.0)u"kJ * mol^-1 * nm^-1",
        atol=1e-9u"kJ * mol^-1 * nm^-1",
    )
    @test isapprox(
        potential_energy(b1, c1, c2, box_size),
        1500.0u"kJ * mol^-1",
        atol=1e-9u"kJ * mol^-1",
    )
    @test isapprox(
        potential_energy(b2, c1, c3, box_size),
        2000.0u"kJ * mol^-1",
        atol=1e-9u"kJ * mol^-1",
    )

    b1 = MorseBond(D=100.0u"kJ * mol^-1", α=10.0u"nm^-1", r0=0.2u"nm")
    b2 = MorseBond(D=200.0u"kJ * mol^-1", α=5.0u"nm^-1" , r0=0.6u"nm")
    fs = force(b1, c1, c2, box_size)
    @test isapprox(
        fs.f1,
        SVector(465.0883158697, 0.0, 0.0)u"kJ * mol^-1 * nm^-1",
        atol=1e-9u"kJ * mol^-1 * nm^-1",
    )
    @test isapprox(
        fs.f2,
        SVector(-465.0883158697, 0.0, 0.0)u"kJ * mol^-1 * nm^-1",
        atol=1e-9u"kJ * mol^-1 * nm^-1",
    )
    fs = force(b2, c1, c3, box_size)
    @test isapprox(
        fs.f1,
        SVector(-9341.5485409432, 0.0, 0.0)u"kJ * mol^-1 * nm^-1",
        atol=1e-9u"kJ * mol^-1 * nm^-1",
    )
    @test isapprox(
        fs.f2,
        SVector(9341.5485409432, 0.0, 0.0)u"kJ * mol^-1 * nm^-1",
        atol=1e-9u"kJ * mol^-1 * nm^-1",
    )
    @test isapprox(
        potential_energy(b1, c1, c2, box_size),
        39.9576400894u"kJ * mol^-1",
        atol=1e-9u"kJ * mol^-1",
    )
    @test isapprox(
        potential_energy(b2, c1, c3, box_size),
        590.4984884025u"kJ * mol^-1",
        atol=1e-9u"kJ * mol^-1",
    )
end

@testset "Spatial" begin
    @test vector1D(4.0, 6.0, 10.0) ==  2.0
    @test vector1D(1.0, 9.0, 10.0) == -2.0
    @test vector1D(6.0, 4.0, 10.0) == -2.0
    @test vector1D(9.0, 1.0, 10.0) ==  2.0

    @test vector1D(4.0u"nm", 6.0u"nm", 10.0u"nm") ==  2.0u"nm"
    @test vector1D(1.0u"m" , 9.0u"m" , 10.0u"m" ) == -2.0u"m"
    @test_throws Unitful.DimensionError vector1D(6.0u"nm", 4.0u"nm", 10.0)

    @test vector(SVector(4.0, 1.0, 6.0), SVector(6.0, 9.0, 4.0),
                    SVector(10.0, 10.0, 10.0)) == SVector(2.0, -2.0, -2.0)
    @test vector(SVector(4.0, 1.0, 1.0), SVector(6.0, 4.0, 3.0),
                    SVector(10.0, 5.0, 3.5)) == SVector(2.0, -2.0, -1.5)
    @test vector(SVector(4.0, 1.0), SVector(6.0, 9.0),
                    SVector(10.0, 10.0)) == SVector(2.0, -2.0)
    @test vector(SVector(4.0, 1.0, 6.0)u"nm", SVector(6.0, 9.0, 4.0)u"nm",
                    SVector(10.0, 10.0, 10.0)u"nm") == SVector(2.0, -2.0, -2.0)u"nm"

    @test wrap_coords(8.0 , 10.0) == 8.0
    @test wrap_coords(12.0, 10.0) == 2.0
    @test wrap_coords(-2.0, 10.0) == 8.0

    @test wrap_coords(8.0u"nm" , 10.0u"nm") == 8.0u"nm"
    @test wrap_coords(12.0u"m" , 10.0u"m" ) == 2.0u"m"
    @test_throws ErrorException wrap_coords(-2.0u"nm", 10.0)

    vels_units   = [maxwell_boltzmann(12.0u"u", 300.0u"K") for _ in 1:1_000]
    vels_nounits = [maxwell_boltzmann(12.0    , 300.0    ) for _ in 1:1_000]
    @test 0.35u"nm * ps^-1" < std(vels_units) < 0.55u"nm * ps^-1"
    @test 0.35 < std(vels_nounits) < 0.55
end

@testset "Neighbor lists" begin
    for neighbor_finder in (DistanceNeighborFinder, TreeNeighborFinder, CellListMapNeighborFinder)
        nf = neighbor_finder(nb_matrix=trues(3, 3), n_steps=10, dist_cutoff=2.0u"nm")
        s = System(
            atoms=[Atom(), Atom(), Atom()],
            coords=[SVector(1.0, 1.0, 1.0)u"nm", SVector(2.0, 2.0, 2.0)u"nm",
                    SVector(5.0, 5.0, 5.0)u"nm"],
            box_size=SVector(10.0, 10.0, 10.0)u"nm",
            neighbor_finder=nf,
        )
        neighbors = find_neighbors(s, s.neighbor_finder; parallel=false)
        @test neighbors.list == [(2, 1, false)] || neighbors.list == [(1, 2, false)]
        if run_parallel_tests
            neighbors = find_neighbors(s, s.neighbor_finder; parallel=true)
            @test neighbors.list == [(2, 1, false)] || neighbors.list == [(1, 2, false)]
        end
        show(devnull, nf)
    end

    # Test passing the box_size and coordinates as keyword arguments to CellListMapNeighborFinder
    coords = [SVector(1.0, 1.0, 1.0)u"nm", SVector(2.0, 2.0, 2.0)u"nm", SVector(5.0, 5.0, 5.0)u"nm"]
    box_size = SVector(10.0, 10.0, 10.0)u"nm"
    neighbor_finder=CellListMapNeighborFinder(
        nb_matrix=trues(3, 3), n_steps=10, x0=coords,
        unit_cell=box_size, dist_cutoff=2.0u"nm",
    )
    s = System(
        atoms=[Atom(), Atom(), Atom()],
        coords=coords, box_size=box_size,
        neighbor_finder=neighbor_finder,
    )
    neighbors = find_neighbors(s, s.neighbor_finder; parallel=false)
    @test neighbors.list == [(2, 1, false)] || neighbors.list == [(1, 2, false)]
    if run_parallel_tests
        neighbors = find_neighbors(s, s.neighbor_finder; parallel=true)
        @test neighbors.list == [(2, 1, false)] || neighbors.list == [(1, 2, false)]
    end
end

@testset "Analysis" begin
    pdb_path = joinpath(data_dir, "1ssu.pdb")
    struc = read(pdb_path, BioStructures.PDB)
    cm_1 = BioStructures.coordarray(struc[1], BioStructures.calphaselector)
    cm_2 = BioStructures.coordarray(struc[2], BioStructures.calphaselector)
    coords_1 = SVector{3, Float64}.(eachcol(cm_1)) / 10 * u"nm"
    coords_2 = SVector{3, Float64}.(eachcol(cm_2)) / 10 * u"nm"
    @test rmsd(coords_1, coords_2) ≈ 2.54859467758795u"Å"
    if run_gpu_tests
        @test rmsd(cu(coords_1), cu(coords_2)) ≈ 2.54859467758795u"Å"
    end

    bb_atoms = BioStructures.collectatoms(struc[1], BioStructures.backboneselector)
    coords = SVector{3, Float64}.(eachcol(BioStructures.coordarray(bb_atoms))) / 10 * u"nm"
    bb_to_mass = Dict("C" => 12.011u"u", "N" => 14.007u"u", "O" => 15.999u"u")
    atoms = [Atom(mass=bb_to_mass[BioStructures.element(bb_atoms[i])]) for i in 1:length(bb_atoms)]
    @test isapprox(radius_gyration(coords, atoms), 11.51225678195222u"Å", atol=1e-6u"nm")
end


text:
using Flux, NNlib, GraphNeuralNetworks, Graphs, SparseArrays
using Documenter

DocMeta.setdocmeta!(GraphNeuralNetworks, :DocTestSetup,
    :(using GraphNeuralNetworks, Graphs, SparseArrays, NNlib, Flux);
    recursive=true)

makedocs(;
    modules=[GraphNeuralNetworks, NNlib, Flux, Graphs, SparseArrays],
    doctest=false, clean=true,
    sitename = "GraphNeuralNetworks.jl",
    pages = ["Home" => "index.md",
             "Graphs" => "gnngraph.md",
             "Message Passing" => "messagepassing.md",
             "Model Building" => "models.md",
             "Datasets" => "datasets.md",
             "API Reference" =>
               [
                "GNNGraph" => "api/gnngraph.md",
                "Basic Layers" => "api/basic.md",
                "Convolutional Layers" => "api/conv.md",
                "Pooling Layers" => "api/pool.md",
                "Message Passing" => "api/messagepassing.md",
                "Utils" => "api/utils.md",
               ],
              "Developer Notes" => "dev.md",
              "Summer Of Code" => "gsoc.md",
            ],
)

deploydocs(repo="github.com/CarloLucibello/GraphNeuralNetworks.jl.git")


text:
module Router

import HTTP
import ..State

include("HealthHandler.jl")

function HttpRouter(state::State.SharedState)
    router = HTTP.Router()
    HTTP.@register(router, "GET", "/health", HealthHandler.Handler(state))

    return router
end

end # module

max_stars_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia

max_issues_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia

max_forks_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia


text:
using MultivariateStats
using LinearAlgebra
using Test
import Statistics: mean, cov, var
import Random

Random.seed!(34568)

## FA with zero mean

X = randn(5, 10)
Y = randn(3, 10)

W = qr(randn(5, 5)).Q[:, 1:3]
Ψ = fill(0.1, 5)
M = FactorAnalysis(Float64[], W, Ψ)

@test indim(M) == 5
@test outdim(M) == 3
@test mean(M) == zeros(5)
@test loadings(M) == W
@test var(M) == Ψ

T = inv(I+W'*diagm(0 => 1 ./ var(M))*W)*W'*diagm(0 => 1 ./ var(M))
@test transform(M, X[:,1]) ≈ T * X[:,1]
@test transform(M, X) ≈ T * X

R = cov(M)*W*inv(W'W)
@test reconstruct(M, Y[:,1]) ≈ R * Y[:,1]
@test reconstruct(M, Y) ≈ R * Y

## PCA with non-zero mean

mval = rand(5)
M = FactorAnalysis(mval, W, Ψ)

@test indim(M) == 5
@test outdim(M) == 3
@test mean(M) == mval
@test loadings(M) == W
@test var(M) == Ψ

@test transform(M, X[:,1]) ≈ T * (X[:,1] .- mval)
@test transform(M, X) ≈ T * (X .- mval)

@test reconstruct(M, Y[:,1]) ≈ R * Y[:,1] .+ mval
@test reconstruct(M, Y) ≈ R * Y .+ mval

## prepare training data

d = 5
n = 1000

R = collect(qr(randn(d, d)).Q)
@test R'R ≈ Matrix(I, 5, 5)
rmul!(R, Diagonal(sqrt.([0.5, 0.3, 0.1, 0.05, 0.05])))

X = R'randn(5, n) .+ randn(5)
mval = vec(mean(X, dims=2))
Z = X .- mval

# facm (default) & faem
fa_methods = [:cm, :em]

fas = FactorAnalysis[]
for method in fa_methods
    let M = fit(FactorAnalysis, X, method=method, maxiter=5000),
        P = projection(M),
        W = loadings(M)
    push!(fas,M)

    @test indim(M) == 5
    @test outdim(M) == 4
    @test mean(M) == mval
    @test P'P ≈ Matrix(I, 4, 4)
    @test all(isapprox.(cov(M), cov(X, dims=2), atol=1e-3))

    M = fit(FactorAnalysis, X; mean=mval, method=method)
    @test loadings(M) ≈ W

    M = fit(FactorAnalysis, Z; mean=0, method=method)
    @test loadings(M) ≈ W

    M = fit(FactorAnalysis, X; maxoutdim=3, method=method)
    P = projection(M)

    @test indim(M) == 5
    @test outdim(M) == 3
    @test P'P ≈ Matrix(I, 3, 3)
    end
end

# compare two algorithms
M1, M2 = fas
@test all(isapprox.(cov(M1), cov(M2), atol=1e-3)) # noise
LL(m, x) = (-size(x,2)/2)*(size(x,1)*log(2π) + log(det(cov(m))) + tr(inv(cov(m))*cov(x, dims=2)))
@test LL(M1, X) ≈ LL(M2, X) # log likelihood

# test that fit works with Float32 values
X2 = convert(Array{Float32,2}, X)
# Float32 input
M = fit(FactorAnalysis, X2; method=:cm, maxoutdim=3)
M = fit(FactorAnalysis, X2; method=:em, maxoutdim=3)

# views
M = fit(FactorAnalysis, view(X2, :, 1:100), method=:cm, maxoutdim=3)
M = fit(FactorAnalysis, view(X2, :, 1:100), method=:em, maxoutdim=3)


text:
@testset "Optimal tranport" begin
    using STMO.OptimalTransport

    C = [1 1 0;
        0 1 1;
        1 0 1]

    @testset "Monge" begin

        @test monge_brute_force(C) == ([3, 1, 2], 0)
        @test monge_brute_force(1.0C) == ([3, 1, 2], 0.0)
    end

    @testset "Sinkhorn" begin
        a, b = [1, 1, 1] / 3, [0.25, 0.6, 0.15]
        Ph = sinkhorn(C, a, b, λ=10, ϵ=1e-10)
        Pl = sinkhorn(C, a, b, λ=10, ϵ=1e-10)

        @test sum(Ph, dims=1)[:] ≈ b
        @test sum(Ph, dims=2)[:] ≈ a
        @test sum(Pl, dims=1)[:] ≈ b
        @test sum(Pl, dims=2)[:] ≈ a

        @test sum(Ph .* C) ≤ sum(Pl .* C)  # lower cost

        h(p) = -sum(p .* log.(p))
        @test h(Ph) ≤ h(Pl)
    end
end


text:
import Pkg

const GROUP = get(ENV, "GROUP", "All")

function dev_subpkg(subpkg)
    subpkg_path = joinpath("lib", subpkg)
    Pkg.develop(Pkg.PackageSpec(path=subpkg_path))
end

function activate_subpkg_env(subpkg)
    subpkg_path = joinpath("lib", subpkg)
    Pkg.activate(subpkg_path)
    Pkg.develop(Pkg.PackageSpec(path=subpkg_path))
    Pkg.instantiate()
end

Pkg.update()

# All packages need the core
dev_subpkg("ArrayInterfaceCore")
Pkg.develop(Pkg.PackageSpec(path="."))

Pkg.test("ArrayInterface"; coverage=true)


text:
println("Master node pid: $(getpid())")

# Load libraries on master node.
# NOTE: In theory, this step is not needed. In practice, this step
# guarantees that libraries are compiled and cached before they
# are used by the workers. If this is not done, workers may
# try to compile the same library simultaneously and write the compiled cache
# to the same locatoin.
println("Compile libraries on master node ..."); flush(stdout)
include("imports.jl")
println("Finished compiling libraries on master node ..."); flush(stdout)

using Distributed

# Number of processors to use.
if length(ARGS) == 0
  println("Running in debug mode! This does not reproduce results in paper!")
  flush(stdout)
  num_proc = 3
  NMCMC = 20
  NBURN = 10
  Kmcmcs = [3, 4, 5]
  mm = 0
else
  num_proc = parse(Int, ARGS[1])
  Kmcmcs = parse.(Int, split(ARGS[2], ','))  # collect(3:3:33)
  mm = parse(Int, ARGS[3])  # 0, 1, or 2. missing mechanism.
  NMCMC = 6000
  NBURN = 10000
end

# Remove old workers
rmprocs(filter(w -> w > 1, workers()))

# Add processeors
println("Adding $(num_proc) workers ..."); flush(stdout)
addprocs(num_proc)

# Print worker process ids
worker_pids = pmap(id -> getpid(), workers())
println("Worker pids: $(join(worker_pids, " "))"); flush(stdout)

# Load the same libraries on worker nodes.
println("Loading libraries on worker nodes ..."); flush(stdout)
@everywhere include("imports.jl")
println("Finished loading libraries on workers node ..."); flush(stdout)

# Create model data object
dat = let
  function separatesamples(Y)
    sample_ids = sort(unique(Y.sample_id))
    return Matrix.(DF.select(Y[Y.sample_id .== i, :], DF.Not(:sample_id))
                   for i in sample_ids)
  end

  path_to_data = "../data/cb_transformed_reduced.csv"
  Y = coalesce.(CSV.read(path_to_data), NaN)
  y = separatesamples(Y)

  CytofResearch.Model.Data(y)
end

@everywhere function fit(Kmcmc::Int, nmcmc::Int, nburn::Int, dat, mm::Int)
  # Directory to store output.
  output_dir = "results/cb-runs/mm_$(mm)/Kmcmc_$(Kmcmc)"

  # y quantiles for missing mechanism
  if mm == 0
    yQuantiles =  [0., 0.25, 0.5]
  elseif mm == 1
    yQuantiles =  [0., 0.2, 0.4]
  else
    yQuantiles =  [0., 0.15, 0.3]
  end

  # Create output directory if needed.
  mkpath(output_dir)

  # Redirect output to a file.
  CytofResearch.Model.redirect_stdout_to_file("$(output_dir)/log.txt") do
    println("Output directory: $(output_dir)"); flush(stdout)
    println("Worker node pid: $(getpid())"); flush(stdout)

    # Converts number of samples to thinning factor.
    nsamps_to_thin(nsamps::Int, nmcmc::Int) = max(1, div(nmcmc, nsamps))

    # Specification of L (as in model)
    Lmcmc = Dict(0 => 5, 1 => 3)

    # Set random seed for reproducibility.
    Random.seed!(0)

    # Make model constants
    @time c = CytofResearch.Model.defaultConstants(
      dat, Kmcmc, Lmcmc,
      tau0=10.0, tau1=10.0,
      sig2_prior=InverseGamma(3.0, 2.0),
      alpha_prior=Gamma(0.1, 10.0),
      yQuantiles=yQuantiles,
      pBounds=[.05, .80, .05],
      similarity_Z=CytofResearch.Model.sim_fn_abs(10000),
      probFlip_Z=2.0 / (dat.J * Kmcmc),
      noisyDist=Normal(0.0, 10.0))

    # Print model constant for sanity.
    CytofResearch.Model.printConstants(c); flush(stdout)

    # Initialize MCMC smartly
    println("Smart initializing ..."); flush(stdout)
    @time init = CytofResearch.Model.smartInit(c, dat)

    # TODO: print probMissing curves

    # NOTE: This is the core of this script. It fits the FAM.
    # Everything else is setup / for parallelization.
    results = CytofResearch.Model.cytof_fit(
      init, c, dat;
      nmcmc=nmcmc, nburn=nburn,
      monitors=[[:Z, :lam, :W, :sig2, :delta, :alpha, :v, :eta, :eps],
                [:y_imputed, :gam]],
      thins=[2, nsamps_to_thin(10, nmcmc)],
      printFreq=10,
      computeDIC=true, computeLPML=true,
      computedden=true,
      sb_ibp=false, use_repulsive=false,
      joint_update_Z=true)

    # Write to disk
    BSON.bson("$(output_dir)/output.bson", results)

    println("DONE!")
  end # end of redirect
end

# Run the thing in parallel
numjobs = length(Kmcmcs)
status = pmap(fit,
              Kmcmcs, fill(NMCMC, numjobs),
              fill(NBURN, numjobs), fill(dat, numjobs),
              fill(mm, numjobs),
              on_error=identity)
status_sanitized = map(s -> s == nothing ? "success" : s, status)

# Printing success / failure status of runs.
for (K, s) in zip(Kmcmcs, status_sanitized)
  println("Kmcmc: $(K) => status: $(s)")
end

println("ALL DONE!")

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
using Flux: onehotbatch
using Distributions: Uniform
import Base.insert!

mutable struct UDRLMemory{S<:Real, A<:Real, R<:Real}
    episodes::Vector{Tuple{R, Vector{Tuple{Vector{S}, A, R}}}}
    obss::Matrix{S}
    dhs::Matrix{S}
    drs::Matrix{S}
    targets::Vector{A}
    maxreward::S
    maxlen::Int64
    minibatchsize::Int64
    function UDRLMemory{S, A, R}(nO, maxreward, maxlen, minibatchsize) where {S, A, R}
        new{S, A, R}(
            Vector{Tuple{R, Vector{Tuple{Vector{S}, Vector{A}, R}}}}(),
            Matrix{S}(undef, nO, minibatchsize),
            Matrix{S}(undef, 1, minibatchsize),
            Matrix{S}(undef, 1, minibatchsize),
            Vector{A}(undef, minibatchsize),
            maxreward,
            maxlen,
            minibatchsize
        )
    end
end

function insert!(memory::UDRLMemory{S, A, R}, r::R, e::Vector{Tuple{Vector{S}, A, R}}) where {S, A, R}
    len = length(memory.episodes)
    i = findfirst(x -> x[1] < r, memory.episodes) |> x -> x === nothing ? len + 1 : x
    if i <= memory.maxlen
        Base.insert!(memory.episodes, i, (r, e))
        len + 1 > memory.maxlen && deleteat!(memory.episodes, memory.maxlen)
    end
end

function command(memory::UDRLMemory{S, A, R}, best) where {S, A, R}
    length(memory.episodes) < best && return (one(S), memory.maxreward)
    h = [length(memory.episodes[i][2]) for i in 1:best]
    r = [memory.episodes[i][1] for i in 1:best]
    rmean = mean(r)
    rstd = std(r, mean=rmean) |> x -> isnan(x) ? zero(R) : x
    return S(mean(h)), S(rand(Uniform(rmean, rmean+rstd+1f-5)))
end

function sample(memory::UDRLMemory, nA)
    @inbounds for i in 1:memory.minibatchsize
        traj = memory.episodes[rand(1:length(memory.episodes))][2]
        t0 = rand(1:length(traj))
        t1 = rand(t0:length(traj))
        memory.obss[:, i] = traj[t0][1]
        memory.dhs[1, i] = t1 + 1 - t0
        memory.drs[1, i] = reduce((x, y) -> x + y[3], traj[t0:t1], init=zero(eltype(memory.drs)))
        memory.targets[i] = traj[t0][2]
    end
    return memory.obss, memory.dhs, memory.drs, onehotbatch(memory.targets, 1:nA)
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
module LifeGame

using Plots
using ProgressMeter

struct Board{T}
    state::T
    tmp::T
end
Board(state) = Board(state, similar(state))
function fivexfive(n=200)
    state = zeros(Int8, n, n)
    m = n ÷ 2
    state[m-2:m+2, m-2:m+2] .= [1 1 1 0 1; 1 0 0 0 0; 0 0 0 1 1; 0 1 1 0 1; 1 0 1 0 1]
    Board(state)
end
randboard(n=200) = Board(rand(Int8[0, 1], n, n))

P(i, m) = ifelse(i == m, 1, i+1)
Q(i, m) = ifelse(i == 1, m, i-1)
function _update!(v, u)
    m, n = size(u)
    @inbounds for j in 1:n, i in 1:m
        i₊, i₋, j₊, j₋ = P(i, m), Q(i, m), P(j, n), Q(j, n)
        N = u[i₋,j₋] + u[i,j₋] + u[i₊,j₋] + u[i₋,j] + u[i₊,j] + u[i₋,j₊] + u[i,j₊] + u[i₊,j₊]
        v[i,j] = N == 3 ? 1 : (N == 2 && !iszero(u[i,j])) ? 1 : 0
    end
end

function update!(board::Board, niters=1)
    state, tmp = board.state, board.tmp
    for _ in 1:niters÷2
        _update!(tmp, state)
        _update!(state, tmp)
    end
    if isodd(niters)
        _update!(tmp, state)
        state .= tmp
    end
end

function gif!(board::Board, niters, nskips=1; gifname="life.gif", fps=20, size=(240, 240))
    prog = Progress(niters, 0)
    anim = @animate for t in 1:niters
        heatmap(board.state; size, colorbar=false, ticks=false, axis=false, frame=false)
        update!(board, nskips)
        next!(prog)
    end
    gif(anim, gifname; fps)
end

end

# %%
boardrandom = LifeGame.randboard()
LifeGame.gif!(boardrandom, 500; gifname="liferandom.gif")

# %%
board5x5 = LifeGame.fivexfive()
LifeGame.gif!(board5x5, 2000; gifname="life5x5.gif", fps=60)

# %%
state = Int8[
    0 0 0 0 0 0 0
    0 0 0 0 0 0 0
    0 0 0 0 0 0 0
    0 0 0 0 0 0 0
    0 1 1 0 0 0 0
    1 0 1 0 0 0 0
    0 0 1 0 0 0 0
]
board = LifeGame.Board(state)
board.state

# %%
LifeGame.update!(board, 12)
board.state

# %%
LifeGame.gif!(board, 27; gifname="lifeglider.gif", fps=5)

# %%

max_stars_repo_path: 0014/Game of Life.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0014/Game of Life.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0014/Game of Life.jl

max_forks_repo_name: genkuroki/public


text:
function MCMCDiagnosticTools.gewekediag(
    chains::Chains;
    sections = _default_sections(chains),
    kwargs...
)
    # Subset the chain.
    _chains = Chains(chains, _clean_sections(chains, sections))

    # Compute statistics for each chain.
    chains_array = _chains.value.data
    results = map(axes(chains_array, 3)) do chain
        vec_of_namedtuples = map(axes(chains_array, 2)) do param
            return MCMCDiagnosticTools.gewekediag(
                cskip(@view(_chains.value.data[:, param, chain])); kwargs...
            )
        end
        namedtuple_of_vecs = Tables.columntable(vec_of_namedtuples)
        return namedtuple_of_vecs
    end

    # Create data frames.
    parameters = (parameters = names(_chains),)
    dfs = [
        ChainDataFrame("Geweke diagnostic - Chain $i", merge(parameters, result))
        for (i, result) in enumerate(results)
    ]

    return dfs
end


text:
import LightXML

#
# Table
#
Base.showable(::MIME"text/html", table::Table) = true

function Base.show(io::IO, a::MIME"text/html", table::Table)
	show(io, a, table.potential)
end

#
# Vector{Table}
#
Base.showable(::MIME"text/html", tables::Vector{Table}) = true

function Base.show(io::IO, a::MIME"text/html", tables::Vector{Table})
	for table in tables
		show(io, a, table)
	end
end

"""
    readxdsl( filename::AbstractString )
Return a DiscreteBayesNet read from the xdsl file
"""
function readxdsl( filename::AbstractString )

    # This currently assumes that the states are Integers
    # Thus, all <state id="???"/> must be integers

    splitext(filename)[2] == ".xdsl" || error("readxdsl only supports .xdsl format")

    xdoc  = LightXML.parse_file(filename)
    xroot = LightXML.root(xdoc)
    ces   = LightXML.get_elements_by_tagname(xroot, "nodes")[1]
    cpts  = collect(LightXML.child_elements(ces))

    varnames = Array{Symbol}(undef, length(cpts))
    for (i,e) in enumerate(cpts)
        id = LightXML.attribute(e, "id")
        varnames[i] = Symbol(id)
    end

    bn = DiscreteBayesNet()

    for (i,e) in enumerate(cpts)

        node_sym = varnames[i]

        for s in LightXML.get_elements_by_tagname(e, "state")
            attr = convert(String, LightXML.attribute(s, "id"))
            @assert(!isa(match(r"\d", attr), Nothing), "All state ids must be integers")
        end

        # set the node's domain
        states   = [parse(Int, match(r"\d", convert(String, LightXML.attribute(s, "id"))).match) for s in LightXML.get_elements_by_tagname(e, "state")]

        probs = Float64[parse(Float64, s) for s in split(LightXML.content(LightXML.find_element(e, "probabilities")))]

        # set any parents & populate probability table
        parents_elem = LightXML.get_elements_by_tagname(e, "parents")
        if !isempty(parents_elem)
            parents = NodeName[Symbol(s) for s in split(LightXML.content(parents_elem[1]))]

            # populate probability table
            reverse!(parents) # because SMILE varies first parent least quickly
            parental_ncategories = _get_parental_ncategories(bn, parents)
            k = length(states)
            Q = prod(parental_ncategories)
            distributions = Array{Categorical}(undef, Q)
            for q in 1:Q
                hi = k*q
                lo = hi - k + 1
                distributions[q] = Categorical{Float64,Vector{Float64}}(probs[lo:hi])
            end

            push!(bn, DiscreteCPD(node_sym, parents, parental_ncategories, distributions))
        else
            # no parents
            push!(bn, DiscreteCPD(node_sym, probs))
        end
    end

    bn
end

"""
    write(io, text/plain, bn)

Writes a text file containing the sufficient statistics for a discrete Bayesian network.
This was inspired by the format listed in Appendix A of
"Correlated Encounter Model for Cooperative Aircraft in the National Airspace System Version 1.0"
by Mykel Kochenderfer.

The text file contains the following parameters:
- variable labels: A space-delimited list specifies the variable labels, which are symbols.
                   The ordering of the variables in this list determines the ordering of the variables
                   in the other tables. Note that the ordering of the variable labels is not
                   necessarily topological.
- graphical structure: A binary matrix is used to represent the graphical structure of the Bayesian
                   network. A 1 in the ith row and jth column means that there is a directed edge
                   from the ith varible to the jth variable in the Bayesian network. The ordering
                   of the variables are as defined in the variable labels section of the file.
                   The entries are 0 or 1 and are not delimited.
- variable instantiations: A list of integers specifying the number of instantiations for each variable.
                   The list is space-delimited.
- sufficient statistics: A list of space-delimited integers Pₐⱼₖ  which specifies the sufficient statistics.
                   The array is ordered first by increasing k, then increasing j, then increasing i.
                   The variable ordering is defined in the variable labels section of the file.
                   The list is a flattened matrices, where each matrix is rₐ × qₐ where rₐ is the number of
                   instantiations of variable a and qₐ is the number of instantiations of the parents of
                   variable a. The ordering is the same as the ordering of the distributions vector in
                   the CategoricalCPD type.
                   The entires in Pₐⱼₖ are floating point probability values.

For example, the network Success -> Forecast
with Success ∈ [1, 2] and P(1) = 0.2, P(2) = 0.8
and Forecast ∈ [1, 2, 3] with
    P(1 | 1) = 0.4, P(2 | 1) = 0.4, P(3 | 1) = 0.2
    P(1 | 2) = 0.1, P(2 | 2) = 0.3, P(3 | 2) = 0.6

Is output as:

Success Forecast
01
00
2 3
2 4 4 1 3
"""
function Base.write(io::IO, mime::MIME"text/plain", bn::DiscreteBayesNet)

    n = length(bn)
    arr_names = names(bn)

    # variable labels
    for (i,name) in enumerate(arr_names)
        print(io, name, i != n ? " " : "\n")
    end

    # graphical structure
    for i in arr_names
        for j in arr_names
            print(io, has_edge(bn, i, j) ? "1" : "0")
        end
        print(io, "\n")
    end

    # variable instantiations
    for (i,name) in enumerate(arr_names)
        cpd = get(bn, name)
        print(io, ncategories(cpd), i != n ? " " : "\n")
    end

    # sufficient statistics
    space = false
    for name in arr_names
        cpd = get(bn, name)
        for D in cpd.distributions
            for p in probs(D)[1:end-1]
                str = @sprintf("%.16g", p)
                print(io, space ? " " : "" , str)
                space = true
            end
        end
    end
    print(io, "\n")
end

function Base.read(io::IO, mime::MIME"text/plain", ::Type{DiscreteBayesNet})

    # variable labels
    arr_names = [Symbol(s) for s in split(readline(io))]
    n = length(arr_names)

    # graphical structure
    adj = falses(n,n)
    for i in 1 : n
        for (j,b) in enumerate(readline(io)[1:n])
            adj[i,j] = b == '1'
        end
    end

    # variable instantiations
    rs = [parse(Int, s) for s in split(readline(io))]

    # sufficient statistics
    stats = split(readline(io)) # strings for now

    # build DBN
    idx = 0

    function pull_next_prob_vector(r)
        probs = Array{Float64}(undef, r)
        for j in 1 : r-1
            probs[j] = parse(Float64, stats[idx += 1])
        end
        probs[end] = 1.0 - sum(probs[1:end-1])
        if probs[end] < 0.0
            probs[end] = 0.0
            probs ./= sum(probs)
        end
        return probs
    end

    cpds = Array{DiscreteCPD}(undef, n)
    for i in 1 : n
        name = arr_names[i]
        r = rs[i]
        parents = findall(adj[:,i])
        if isempty(parents)
            probs = pull_next_prob_vector(r)
            cpds[i] = DiscreteCPD(name, probs)
        else
            parent_names = arr_names[parents]
            parental_ncategories = rs[parents]
            Q = prod(parental_ncategories)
            distributions = Array{Categorical}(undef, Q)
            for q in 1 : Q
                probs = pull_next_prob_vector(r)
                distributions[q] = Categorical(probs)
            end
            cpds[i] = DiscreteCPD(name, parent_names, parental_ncategories, distributions)
        end
    end

    return BayesNet(cpds)
end


text:
#####
##### Game trace
#####

"""
    MuTrace{State}

An object that collects all states visited during a game, along with the
rewards obtained at each step and the successive player policies to be used
as targets for the neural network.

# Constructor

    MuTrace(initial_state)

"""
struct MuTrace{State}
  states :: Vector{State}
  policies :: Vector{Vector{Float64}}
  rewards :: Vector{Float64}
  actions :: Vector{Int}
  rootvalues :: Vector{Float64}
  function MuTrace(init_state)
    return new{typeof(init_state)}([init_state], [], [], [], [])
  end
end

function valid_trace(t::MuTrace)
  return length(t.policies) == length(t.rewards) == length(t.actions) == length(t.states)-1
end

"""
    Base.push!(t::Trace, π, r, s, a, v)

Add a (target policy, reward, new state, selected action, rootvalue) to a trace.
"""
function Base.push!(t::MuTrace, π, r, s, a, v)
  push!(t.states, s)
  push!(t.policies, π)
  push!(t.rewards, r)
  push!(t.actions, a)
  push!(t.rootvalues, v)
end

function Base.length(t::MuTrace)
  return length(t.rewards)
end

function total_reward(t::MuTrace, gamma=1.)
  return sum(gamma^(i-1) * r for (i, r) in enumerate(t.rewards))
end

# function total_reward(t::Trace, gamma=1.)
#   return sum([gamma^(i-1) * r for (i, r) in enumerate(t.rewards)])
# end

# function debug_trace(gspec::AbstractGameSpec, t::Trace)
#   n = length(t)
#   for i in 1:n
#     println("Transition $i:")
#     game = GI.init(gspec, t.states[i])
#     GI.render(game)
#     for (a, p) in zip(GI.available_actions(game),  t.policies[i])
#       print("$(GI.action_string(gspec, a)): $(fmt(".3f", p))  ")
#     end
#     println("")
#     println("Obtained reward of: $(t.rewards[i]).")
#     println("")
#   end
#   println("Showing final state:")
#   GI.render(GI.init(gspec, t.states[n + 1]))
# end


text:
# Julia version of Boyd et al's Matlab script available on:
#   https://web.stanford.edu/~boyd/papers/admm/nonconvex/regressor_sel.html

# regressor_sel  Solve lasso problem via ADMM
#
#  x, history = regressor_sel(A, b, K, rho, alpha)
#
# Attempts to solve the following problem via ADMM:
#
#   minimize || Ax - b ||_2^2
#   subject to card(x) <= K
#
# where card() is the number of nonzero entries.
#
# The solution is returned in the vector x.
#
# history is a structure that contains the objective value, the primal and
# dual residual norms, and the tolerances for the primal and dual residual
# norms at each iteration.
#
# rho is the augmented Lagrangian parameter.
#
# alpha is the over-relaxation parameter (typical values for alpha are
# between 1.0 and 1.8).
#
# More information can be found in the paper linked at:
# http://www.stanford.edu/~boyd/papers/distr_opt_stat_learning_admm.html
#
function regressor_sel(A, b, K, rho;
    MaxIterations = 1000, Quiet = true, ABSTOL = 1e-4, RELTOL = 1e-2)

    # Data preprocessing and setup
    K = iceil(K)
    m, n = size(A)
    Atb = A' * b            # save a matrix-vector multiply

    local r_norm::Float64, s_norm::Float64, eps_pri::Float64, eps_dual::Float64
    x = zeros(n)
    z = zeros(n)
    u = zeros(n)
    L, U = factor(A, rho)   # cache the factorization

    for k in 1:MaxIterations
        # x-update
        q = Atb + rho*(z - u)    # temporary value
        if ( m >= n )    # if skinny
            x = U \ (L \ q)
        else            # if fat
            x = q/rho - (A'*(U \ ( L \ (A*q) )))/rho^2
        end
        x = vec(x)

        # z-update with relaxation
        zold = z
        z = keep_largest!(x + u, K)

        # u-update
        u = u + (x - z)

        # Check convergence
        r_norm  = norm(x - z)
        s_norm  = norm(-rho*(z - zold))
        eps_pri = sqrt(n) * ABSTOL + RELTOL * max(norm(x), norm(-z))
        eps_dual= sqrt(n) * ABSTOL + RELTOL * norm(rho*u)

        if !Quiet
            println(@sprintf("%3d\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.2f\n", k,
                 r_norm, eps_pri, s_norm, eps_dual, objective(A,b,x)))
        end

        if (r_norm < eps_pri) && (s_norm < eps_dual)
            break
        end
    end

    return z
end

function objective(A, b, x)
    sum( (A*x - b).^2 )
end

function keep_largest!(z, K)
    pos = sortperm(abs(z), rev=true)
    z[pos[K+1:end]] = 0
    z
end

import Base.factor
function factor(A, rho)
    m, n = size(A)
    if ( m >= n )  # if skinny
        L = chol( A'*A + rho*speye(n), :L )
    else            # if fat
        L = chol( speye(m) + 1/rho*(A*A'), :L )
    end

    return L, L'
end