text:
@testset "Constructors and basic properties" begin
    # numbers
    for T in (Float64, ComplexF64), V in (BinningAnalysis.Variance, BinningAnalysis.FastVariance)
        B = LogBinner(T; accumulator=V)
        B2 = LogBinner(T; accumulator=V)

        @test length(B) == 0
        @test ndims(B) == 0
        @test isempty(B)
        @test eltype(B) == T
        @test capacity(B) == 2^32 - 1
        @test BinningAnalysis.nlevels(B) == 32
        @test B == B2
        @test !(B != B2)

        B3 = LogBinner(B, capacity=10_000)
        @test B3 == B
        @test capacity(B3) == 16383

        buffer = rand(T, 1000)
        append!(B, buffer)
        @test length(B) == 1000
        @test !isempty(B)
        @test !(B == B2)
        @test B != B2

        @test B3 != B
        @test isempty(B3)
        B3 = LogBinner(B, capacity=1000)
        @test B3 == B
        @test !isempty(B3)

        append!(B2, buffer)
        @test B == B2
        @test !(B != B2)

        append!(B, rand(T, 24))
        @test_throws OverflowError LogBinner(B, capacity=1000)

        empty!(B)
        @test length(B) == 0
        @test isempty(B)
    end

    # arrays
    for T in (Float64, ComplexF64), V in (BinningAnalysis.Variance, BinningAnalysis.FastVariance)
        B = LogBinner(zeros(T, 2, 3), accumulator=V)
        B2 = LogBinner(zeros(T, 2, 3), accumulator=V)

        @test length(B) == 0
        @test ndims(B) == 2
        @test isempty(B)
        @test eltype(B) == Array{T, 2}
        @test capacity(B) == 2^32 - 1
        @test BinningAnalysis.nlevels(B) == 32
        @test B == B2
        @test !(B != B2)

        B3 = LogBinner(B, capacity=10_000)
        @test B3 == B
        @test capacity(B3) == 16383

        buffer = [rand(T, 2, 3) for _ in 1:1000]
        append!(B, buffer)
        @test length(B) == 1000
        @test !isempty(B)
        @test !(B == B2)
        @test B != B2

        @test B3 != B
        @test isempty(B3)
        B3 = LogBinner(B, capacity=1000)
        @test B3 == B
        @test !isempty(B3)

        append!(B2, buffer)
        @test B == B2
        @test !(B != B2)

        append!(B, [rand(T, 2, 3) for _ in 1:24])
        @test_throws OverflowError LogBinner(B, capacity=1000)

        empty!(B)
        @test length(B) == 0
        @test isempty(B)
    end

    # Constructor arguments
    B = LogBinner(capacity=12345)
    @test capacity(B) == 16383
    @test_throws ArgumentError LogBinner(capacity=0)
    @test_throws ArgumentError LogBinner(capacity=-1)

    # Test error on overflow
    B = LogBinner(capacity=1)
    push!(B, 1.0)
    @test_throws OverflowError push!(B, 2.0)

    B = LogBinner(zeros(2,2), capacity=1)
    push!(B, rand(2,2))
    @test_throws OverflowError push!(B, rand(2,2))

    # time series constructor (#26)
    x = rand(10)
    B = LogBinner(x)
    @test length(B) == 10
    @test mean(B) ≈ mean(x)

    # Test equality of different sizes
    B2 = LogBinner()
    @test !(B == B2)
    @test !(B2 == B)
    @test B != B2
    @test B2 != B
    append!(B2, x)
    @test B == B2
    @test B2 == B
    @test !(B != B2)
    @test !(B2 != B)

    # Again for Array inputs
    x = [rand(2,3) for _ in 1:5]
    B = LogBinner(x)
    @test length(B) == 5
    @test mean(B) ≈ mean(x)

    B2 = LogBinner(zeros(2,3))
    @test !(B == B2)
    @test !(B2 == B)
    @test B != B2
    @test B2 != B
    append!(B2, x)
    @test B == B2
    @test B2 == B
    @test !(B != B2)
    @test !(B2 != B)
end

@testset "Checking converging data (Real)" begin
    BA = LogBinner()

    # block of maximally correlated values:
    N_corr = 16

    # number of blocks
    N_blocks = 131_072 # 2^17

    for _ in 1:N_blocks
        x = rand()
        for __ in 1:N_corr
            push!(BA, x)
        end
    end

    # BA must diverge until 16 values are binned
    for lvl in 1:4
        @test !has_converged(BA, lvl)
    end

    # Afterwards it must converge (for a while)
    for lvl in 5:8
        @test has_converged(BA, lvl)
    end
    # Later values may fluctuate due to small samples / small threshold in
    # has_converged

    # means should be consistent
    means = BinningAnalysis.all_means(BA)
    for x in means
        @test x ≈ means[1]
    end
end

@testset "Check variance for complex values" begin
    # NOTE
    # Due to the different (mathematically equivalent) versions of the variance
    # calculated here, the values are onyl approximately the same. (Float error)
    Random.seed!(1234)
    xs = rand(ComplexF64, 1_000_000)
    BA = LogBinner(ComplexF64)

    # Test small set (off by one errors are large here)
    for x in xs[1:10]; push!(BA, x) end
    @test var(BA, 1) ≈ var(xs[1:10])
    @test varN(BA, 1) ≈ var(xs[1:10])/10

    # Test full set
    for x in xs[11:end]; push!(BA, x) end
    @test var(BA, 1) ≈ var(xs)
    @test varN(BA, 1) ≈ var(xs)/1_000_000

    # all_* methods
    @test isapprox(all_vars(BA), [0.16671474067121222, 0.08324845751233179, 0.041527133392489035, 0.020847602123934883, 0.010430741538377142, 0.005111097271805531, 0.0025590988213273214, 0.001283239131297187, 0.0006322480081128456, 0.0003060164750540162, 0.00015750782337442537, 8.006142368921498e-5, 3.810111634139357e-5, 1.80535512880331e-5, 1.0002438211476061e-5, 5.505102193326117e-6, 2.8788397929968568e-6, 1.7242475507384114e-6, 7.900888818745955e-7])
    @test isapprox(all_varNs(BA), zero(all_varNs(BA)), atol=1e-6)
    @test isapprox(all_taus(BA), [0.0, -0.00065328850247931, -0.0018180968845809553, 0.00019817179932868356, 0.0005312186016332987, -0.009476150581268994, -0.008794711536776634, -0.007346737569564443, -0.014542478848703244, -0.030064159934323653, -0.01599670814224563, -0.007961042363178128, -0.03167873601168558, -0.056188229083248886, -0.008218660661725774, 0.05035147373711113, 0.0756019296606737, 0.2387501479629205, 0.289861051172009])
    @test isapprox(all_std_errors(BA), [0.0004083071646092097, 0.0004080403350462593, 0.00040756414657076514, 0.0004083880715587553, 0.0004085240073900606, 0.00040441947616030687, 0.00040470028980091994, 0.0004052963382191276, 0.00040232555162611277, 0.00039584146247874917, 0.00040172249945971054, 0.00040504357104527986, 0.00039516087376314515, 0.0003846815937765092, 0.00040493752222959487, 0.0004283729758565588, 0.00043808977717638777, 0.0004963074437049236, 0.0005131890106236348])

    @test isapprox(tau(BA), -0.008218660661725774)
    @test isapprox(std_error(BA), 0.00040493752222959487)
end

@testset "Check variance for complex vectors" begin
    Random.seed!(1234)
    xs = [rand(ComplexF64, 3) for _ in 1:1_000_000]
    BA = LogBinner(zeros(ComplexF64, 3))

    # Test small set (off by one errors are large here)
    for x in xs[1:10]; push!(BA, x) end
    @test var(BA, 1) ≈ var(xs[1:10])
    @test varN(BA, 1) ≈ var(xs[1:10])/10

    # Test full set
    for x in xs[11:end]; push!(BA, x) end
    @test var(BA, 1) ≈ var(xs)
    @test varN(BA, 1) ≈ var(xs)/1_000_000

    # all_std_errors for <:AbstractArray
    @test all(isapprox.(all_std_errors(BA), Ref(zeros(3)), atol=1e-2))

    @test all(isapprox.(tau(BA), [-0.101203, -0.0831874, -0.0112827], atol=1e-6))
    @test all(isapprox.(std_error(BA), [0.000364498, 0.00037268, 0.000403603], atol=1e-6))
end

@testset "Checking converging data (Complex)" begin
    BA = LogBinner(ComplexF64)

    # block of maximally correlated values:
    N_corr = 16

    # number of blocks
    N_blocks = 131_072 # 2^17

    for _ in 1:N_blocks
        x = rand(ComplexF64)
        for __ in 1:N_corr
            push!(BA, x)
        end
    end

    # BA must diverge until 16 values are binned
    for lvl in 1:4
        @test !has_converged(BA, lvl)
    end

    # Afterwards it must converge (for a while)
    for lvl in 5:8
        @test has_converged(BA, lvl)
    end
    # Later values may fluctuate due to small samples / small threshold in
    # has_converged

    # means should be consistent
    means = BinningAnalysis.all_means(BA)
    for x in means
        @test x ≈ means[1]
    end
end

@testset "Checking converging data (Vector)" begin
    BA = LogBinner(zeros(3))

    # block of maximally correlated values:
    N_corr = 16

    # number of blocks
    N_blocks = 131_072 # 2^17

    for _ in 1:N_blocks
        x = rand(Float64, 3)
        for __ in 1:N_corr
            push!(BA, x)
        end
    end

    # BA must diverge until 16 values are binned
    for lvl in 1:4
        @test !has_converged(BA, lvl)
    end

    # Afterwards it must converge (for a while)
    for lvl in 5:8
        @test has_converged(BA, lvl)
    end
    # Later values may fluctuate due to small samples / small threshold in
    # has_converged

    # means should be consistent
    means = BinningAnalysis.all_means(BA)
    for i in eachindex(means)
        @test means[i] ≈ means[1]
    end
end

@testset "Type promotion" begin
    Bf = LogBinner(zero(1.)) # Float64 LogBinner
    Bc = LogBinner(zero(im)) # Float64 LogBinner

    # Check that this doesn't throw (TODO: is there a better way?)
    @test (append!(Bf, rand(1:10, 10000)); true)
    @test (append!(Bc, rand(10000)); true)
end

@testset "Sum-type heuristic" begin
    # numbers
    @test typeof(LogBinner(zero(Int64))) == LogBinner{Float64, 32, BinningAnalysis.Variance{Float64}}
    @test typeof(LogBinner(zero(ComplexF16))) == LogBinner{ComplexF64, 32, BinningAnalysis.Variance{ComplexF64}}

    # arrays
    @test typeof(LogBinner(zeros(Int64, 2,2))) == LogBinner{Matrix{Float64}, 32, BinningAnalysis.Variance{Matrix{Float64}}}
    @test typeof(LogBinner(zeros(ComplexF16, 2,2))) == LogBinner{Matrix{ComplexF64}, 32, BinningAnalysis.Variance{Matrix{ComplexF64}}}
end

@testset "Indexing Bounds" begin
    BA = LogBinner(zero(Float64); capacity=1)
    for func in [:var, :varN, :mean, :tau]
        # check if func(BA, 0) throws BoundsError
        # It should as level 1 is now the initial level
        @test_throws BoundsError @eval $func($BA, 0)
        # Check that level 1 exists
        @test (@eval $func($BA, 1); true)
        # Check that level 2 throws BoundsError
        @test_throws BoundsError @eval $func($BA, 2)
    end
end

@testset "_reliable_level" begin
    BA = LogBinner()
    # Empty Binner
    @test BinningAnalysis._reliable_level(BA) == 1
    @test isnan(std_error(BA, BinningAnalysis._reliable_level(BA)))

    # One Element should still return NaN (due to 1/(n-1))
    push!(BA, rand())
    @test BinningAnalysis._reliable_level(BA) == 1
    @test isnan(std_error(BA, BinningAnalysis._reliable_level(BA)))

    # Two elements should return some value
    push!(BA, rand())
    @test BinningAnalysis._reliable_level(BA) == 1
    @test !isnan(std_error(BA, BinningAnalysis._reliable_level(BA)))

    # same behavior up to (including) 63 values (31 binned in first binned lvl)
    append!(BA, rand(61))
    @test BinningAnalysis._reliable_level(BA) == 1
    @test !isnan(std_error(BA, BinningAnalysis._reliable_level(BA)))

    # at 64 or more values, the lvl should be increasing
    push!(BA, rand())
    @test BinningAnalysis._reliable_level(BA) == 2
    @test !isnan(std_error(BA, BinningAnalysis._reliable_level(BA)))
end

@testset "Cosmetics (show, print, etc.)" begin
    B = LogBinner();
    # empty binner
    oldstdout = stdout
    (read_pipe, write_pipe) = redirect_stdout()
    println(B) # compact
    show(write_pipe, MIME"text/plain"(), B) # full
    redirect_stdout(oldstdout);
    close(write_pipe);

    # compact
    @test readline(read_pipe) == "LogBinner{Float64,32}()"
    # full
    @test readline(read_pipe) == "LogBinner{Float64,32}"
    @test readline(read_pipe) == "| Count: 0"
    @test length(readlines(read_pipe)) == 0
    close(read_pipe);

    # filled binner
    Random.seed!(1234)
    append!(B, rand(1000))
    (read_pipe, write_pipe) = redirect_stdout()
    show(write_pipe, MIME"text/plain"(), B)
    redirect_stdout(oldstdout);
    close(write_pipe);
    @test readline(read_pipe) == "LogBinner{Float64,32}"
    @test readline(read_pipe) == "| Count: 1000"
    @test readline(read_pipe) == "| Mean: 0.49685"
    @test readline(read_pipe) == "| StdError: 0.00733"
    @test length(readlines(read_pipe)) == 0
    close(read_pipe);
end


text:
using Documenter, StaticPolynomials

makedocs(
    sitename = "StaticPolynomials.jl",
    pages = [
        "Introduction" => "index.md",
        "Reference" => "reference.md",
        ]
)

deploydocs(
    repo   = "github.com/JuliaAlgebra/StaticPolynomials.jl.git",
)

ext: jl


text:
# Select a subset of columns
# TODO more general "select" with renaming, mapping, etc?

struct Project{names}
end
@inline Project(names::Tuple{Vararg{Symbol}}) = Project{names}()

@inline Project(n1::Symbol) = Project{(n1,)}()
@inline Project(n1::Symbol, n2::Symbol) = Project{(n1, n2)}()
@inline Project(n1::Symbol, n2::Symbol, n3::Symbol) = Project{(n1, n2, n3)}()
@inline Project(n1::Symbol, n2::Symbol, n3::Symbol, n4::Symbol) = Project{(n1, n2, n3, n4)}()
# TODO generalize this - constant propagation doesn't like slurping.

@inline project(x, names::Tuple{Vararg{Symbol}}) = Project{names}()(x)

@generated function (::Project{names2})(nt::NamedTuple{names}) where {names, names2}
    exprs = [:(getproperty(nt, $(Expr(:quote, n)))) for n in names2]
    return quote
        @_inline_meta
        return NamedTuple{names2}(tuple($(exprs...)))
    end
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module Processing

using ModernGL, GLFW, Colors, Tau, Images
using FreeType, FreeTypeAbstraction
using GLAbstraction, GeometryTypes, Packing

# the following isn't defined in ModernGL.jl yet.
# it is borrowed from the GLVisualize.jl package.
const GL_TEXTURE_MAX_ANISOTROPY_EXT = 0x84FE

include("constants.jl")
include("shaders.jl")

shaderBank = Dict("basicShapes" => GLuint(0),
                "texturedShapes" => GLuint(0),
                "fontDrawing" => GLuint(0),
                "drawFramebuffer" => GLuint(0))

type stateStruct
    bgCol::Array{RGB{Float64}, 1}
    fillStuff::Bool
    fillCol::Array{RGB{Float64}, 1}
    strokeStuff::Bool
    strokeCol::Array{RGB{Float64}, 1}
    tintStuff::Bool
    tintCol::Array{RGB{Float64}, 1}
    program::GLuint
    drawTexture::Bool
    aspectRatio::Float32
    preserveAspectRatio::Bool
    fbSize::Tuple{Int32,Int32}
    strokeWeight::Float32
    fontFace::AbstractString
    textSize::Float32
    height::Int
    width::Int
    left::Float32
    right::Float32
    top::Float32
    bottom::Float32
    cMode::AbstractString
    title::AbstractString
    ellipseMode::AbstractString
    rectMode::AbstractString
    shapeMode::AbstractString
    imageMode::AbstractString
    frameRate::Int
    frameCount::Int
    fontsInitialized::Bool
    alternateShader::Union{Shader,Void}
end

# need to generalize font system
state = stateStruct([RGB(0.8, 0.8, 0.8)], true, [RGB(1.0, 1.0, 1.0)], true, [RGB(0.0, 0.0, 0.0)], false, [RGB(0.0, 0.0, 0.0)], GLuint(0), false, 1., true, (Int32(0), Int32(0)), 1.0, "", 0.4, 275, 275, -1., 1., 1., -1., "RGB", "Processing.jl", "CENTER", "CORNER", "CORNER", "CORNER", 60, 0, false, nothing)

# by default, use system fonts that are known to basically always be available
@windows_only state.fontFace = "C:/Windows/Fonts/arial.ttf"
@linux_only state.fontFace = "/usr/share/fonts/DejaVu/DejaVuSansMono.ttf"
@osx_only state.fontFace = "/System/Library/Fonts/Menlo.ttc"

type vertexStruct
    shapeVertices::Array{GLfloat, 1}
    textureCoords::Array{GLfloat, 1}
    vertexStride::Int
    nVertices::Int
    shapeType::GLuint
end

shapeData = vertexStruct(GLfloat[], GLfloat[], -1, -1, GL_POINTS)

include("openglaux.jl")
include("color.jl")
include("environment.jl")
include("image.jl")
include("input.jl")
include("pixels.jl")
include("rendering.jl")
include("shapes2d.jl")
include("shapes3d.jl")
include("textures.jl")

type GLmatStruct
    currMatrix::Array{GLfloat, 2}
    matrixStack::Array{Array{GLfloat, 2}, 1}
end

# start with default identity matrix, as expected.
GLmatState = GLmatStruct(GLfloat[1 0 0 0;
                                0 1 0 0;
                                0 0 1 0;
                                0 0 0 1],
                                GLfloat[])
include("transform.jl")

include("typography.jl")

type textCharacter
    size::Array{Int, 1}
    bearing::Array{Float64, 1}
    advance::Array{Float64, 1}
    atlasOffset::Float64
end

blankChar = textCharacter(Float64[], Float64[], Float64[], 0.)

type fontStruct
    characters::Dict{Char, textCharacter}
    textAtlas::GLuint
    atlasWidth::Float64
    atlasHeight::Float64
    maxLineLength::Int
    face::Array{Ptr{FreeType.FT_FaceRec}, 1}
    fontWidth::Int
    fontHeight::Int
end

# 80 characters is choosen as the default line length, just like on Unix :)
fontState = fontStruct(Dict(' ' => blankChar), GLuint(0), 0., 0., 80, Ptr{FreeType.FT_FaceRec}[], 0, 48)

function __init__()
    FreeTypeAbstraction.init() # initialize FreeType
    fontState.face = newface(state.fontFace)
    FreeTypeAbstraction.setpixelsize(fontState.face, fontState.fontWidth, fontState.fontHeight)
end

# this will be the texture for rendering into our own framebuffer.
texColorBuffer = GLuint[0]
# this contains the coordinates for the framebuffer texture that we draw at
# the end of each render loop.
# yes, this is illegible.
quadPosCoords = GLfloat[-1.0, 1.0,
                        1.0, 1.0,
                        1.0, -1.0,
                        1.0, -1.0,
                        -1.0, -1.0,
                        -1.0, 1.0]
quadTexCoords = GLfloat[0.0, 1.0,
                        1.0, 1.0,
                        1.0, 0.0,
                        1.0, 0.0,
                        0.0, 0.0,
                        0.0, 1.0]

type GLobjs
    vaos::Array{GLuint, 1}
    posvbos::Array{GLuint, 1}
    posind::Int
    colvbos::Array{GLuint, 1}
    colind::Int
    texvbos::Array{GLuint, 1}
    texind::Int
    ebos::Array{GLuint, 1}
    fbos::Array{GLuint, 1}
    rbos::Array{GLuint, 1}
end

globjs = GLobjs(GLuint[], GLuint[], 0, GLuint[], 0, GLuint[], 0, GLuint[], GLuint[], GLuint[])

export screen, animate, endDrawing, drawingWindow

function screen(width, height; fullScreen=false, preserveAspectRatio=true, GLmajv=3, GLminv=2)
    GLFW.Init()

    GLFW.WindowHint(GLFW.VISIBLE, true)

    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, GLmajv)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, GLminv)
    GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
    GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE)

    # try to activate 10 bpc support
    GLFW.WindowHint(GLFW.RED_BITS, 10)
    GLFW.WindowHint(GLFW.GREEN_BITS, 10)
    GLFW.WindowHint(GLFW.BLUE_BITS, 10)
    GLFW.WindowHint(GLFW.ALPHA_BITS, 8)

    # anti-aliasing by default
    GLFW.WindowHint(GLFW.SAMPLES, 4)
    GLFW.WindowHint(GLFW.RESIZABLE, false)

    if state.frameRate != 60
        GLFW.WindowHint(GLFW.REFRESH_RATE, state.frameRate)
    end

    if fullScreen
        priMon = GLFW.GetPrimaryMonitor()
        vidMode = GLFW.GetVideoMode(priMon)
        state.width = vidMode.width
        state.height = vidMode.height
        window = GLFW.CreateWindow(state.width, state.height, state.title, priMon)
    else
        state.width = width
        state.height = height
        window = GLFW.CreateWindow(state.width, state.height, state.title)
    end

    state.aspectRatio = state.width/state.height
    state.preserveAspectRatio = preserveAspectRatio

    GLFW.MakeContextCurrent(window)
    GLFW.ShowWindow(window)
    GLFW.SetWindowSize(window, state.width, state.height)

    state.fbSize = GLFW.GetFramebufferSize(window)

    glViewport(0, 0, state.fbSize[1], state.fbSize[2])
    GLFW.SwapInterval(1)

    createcontextinfo()

    # framebuffer stuff

    @windows_only begin glGenFramebuffer() end
    # (at least on the one windows 7 machine with an ATI card that I tested,
    # the first call to glGenFramebuffer fails and all subsequent calls
    # succeeed)
    push!(globjs.fbos,  glGenFramebuffer())
    glBindFramebuffer(GL_FRAMEBUFFER, globjs.fbos[1])

    glActiveTexture(GL_TEXTURE0)
    glGenTextures(1, texColorBuffer)
    glBindTexture(GL_TEXTURE_2D, texColorBuffer[1])

    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA32F, state.fbSize[1], state.fbSize[2], 0, GL_RGBA, GL_FLOAT, C_NULL)

    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)

    glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D, texColorBuffer[1], 0)

    # we also need to enable a depth buffer to attach to the framebuffer
    push!(globjs.rbos,  glGenRenderbuffer())
    glBindRenderbuffer(GL_RENDERBUFFER, globjs.rbos[1])
    glRenderbufferStorage(GL_RENDERBUFFER, GL_DEPTH_COMPONENT, state.fbSize[1], state.fbSize[2])
    glFramebufferRenderbuffer(GL_FRAMEBUFFER, GL_DEPTH_ATTACHMENT, GL_RENDERBUFFER, globjs.rbos[1])

    if glCheckFramebufferStatus(GL_FRAMEBUFFER) != GL_FRAMEBUFFER_COMPLETE
        error("framebuffer not initialized.")
    end

    initShaders()

    glPointSize(2) # points are too small to see when 1-pixel big at high res

    # don't re-enable this until we have a better way of making
    # sure that vertices always wind in the correct direction.
    # currently, this is only enabled during text writing.
    # glEnable(GL_CULL_FACE)
    # glCullFace(GL_BACK)
    # glFrontFace(GL_CW)

    glEnable(GL_BLEND) # using blending by default
    glBlendFuncSeparate(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA, GL_ONE, GL_ONE)

    glEnable(GL_MULTISAMPLE) # just be double sure that multisampling is on

    # glEnable(GL_TEXTURE_2D)

    glEnable(GL_LINE_SMOOTH)
    # glEnable(GL_POINT_SMOOTH)
    # glHint(GL_POINT_SMOOTH_HINT, GL_NICEST)

    glEnable(GL_DEPTH_TEST)
    glDepthFunc(GL_LEQUAL)

    if !state.fontsInitialized
        setupFontCharacters()
    end
    state.fontsInitialized = true

    # make the color attachment of the frame buffer the default drawing
    # location
    glDrawBuffers(1, [GL_COLOR_ATTACHMENT0])
    glBindFramebuffer(GL_FRAMEBUFFER, globjs.fbos[1])
    glEnable(GL_DEPTH_TEST)
    shader("basicShapes")
    glClearColor(state.bgCol[1].r, state.bgCol[1].g, state.bgCol[1].b, 1.0)
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

    glBindFramebuffer(GL_FRAMEBUFFER, 0)
    glClearColor(1.0, 1.0, 1.0, 1.0)
    glClear(GL_COLOR_BUFFER_BIT)
    glDisable(GL_DEPTH_TEST)
    shader("drawFramebuffer")

    # draw framebuffer texture to screen
    glActiveTexture(GL_TEXTURE0)
    glBindTexture(GL_TEXTURE_2D, texColorBuffer[1])
    glDrawArrays(GL_TRIANGLES, 0, 6)

    GLFW.SwapBuffers(window)
    GLFW.PollEvents()

    glBindFramebuffer(GL_FRAMEBUFFER, globjs.fbos[1])
    glEnable(GL_DEPTH_TEST)
    shader("basicShapes")
    glClearColor(state.bgCol[1].r, state.bgCol[1].g, state.bgCol[1].b, 1.0)
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

    # for some reason, Windows 7 loses contact with the rendering context if
    # this check isn't made before every draw call
    @windows_only begin
        if GLFW.WindowShouldClose(window)
            GLFW.DestroyWindow(window)
        end
    end

    glErrorMessage()

    return window
end

screen(w; fullScreen=false, preserveAspectRatio=true, GLmajv=3, GLminv=2) = screen(w, w; fullScreen=fullScreen, preserveAspectRatio=preserveAspectRatio, GLmajv=GLmajv, GLminv=GLminv)
screen(; fullScreen=false, preserveAspectRatio=true, GLmajv=3, GLminv=2) = screen(state.width, state.height; fullScreen=fullScreen, preserveAspectRatio=preserveAspectRatio, GLmajv=GLmajv, GLminv=GLminv)

function animate(window)
    glBindFramebuffer(GL_FRAMEBUFFER, 0)
    # glClearColor(1.0, 1.0, 1.0, 1.0)
    # glClear(GL_COLOR_BUFFER_BIT)
    glDisable(GL_DEPTH_TEST)

    shader("drawFramebuffer")

    # draw framebuffer texture to screen
    glActiveTexture(GL_TEXTURE0)
    glBindTexture(GL_TEXTURE_2D, texColorBuffer[1])
    glDrawArrays(GL_TRIANGLES, 0, 6)

    GLFW.SwapBuffers(window)
    # for some reason, Windows 7 loses contact with the rendering context if
    # this command isn't run before every draw call
    @windows_only begin
        GLFW.PollEvents()
    end

    glBindFramebuffer(GL_FRAMEBUFFER, globjs.fbos[1])
    glEnable(GL_DEPTH_TEST)
    if state.alternateShader == nothing
        shader("basicShapes")
    else
        shader(state.alternateShader)
    end

    # for some reason, Windows 7 loses contact with the rendering context if
    # this check isn't made before every draw call
    @windows_only begin
        if GLFW.WindowShouldClose(window)
            GLFW.DestroyWindow(window)
        end
    end

    state.frameCount += 1
end

function animate(f::Function, window)
    while true
        f()
        animate(window)
    end
end

function endDrawing(window)
    GLFW.DestroyWindow(window)
end

function drawingWindow(window)
    GLFW.MakeContextCurrent(window)
    GLFW.ShowWindow(window)
end

end # module


text:
module CommunalHelperRailedMoveBlock

using ..Ahorn, Maple

const steeringModes = ["Horizontal", "Vertical", "Both"]

@mapdef Entity "CommunalHelper/RailedMoveBlock" RailedMoveBlock(
    x::Integer,
    y::Integer,
    width::Integer=16,
    height::Integer=16,
    nodes::Array{Tuple{Integer, Integer}, 1}=Tuple{Integer, Integer}[],
    steeringMode::String="Horizontal",
    speed::Number=120.0,
)

const placements = Ahorn.PlacementDict(
    "Railed Move Block ($steeringMode) (Communal Helper)" => Ahorn.EntityPlacement(
        RailedMoveBlock,
        "rectangle",
        Dict{String, Any}(
            "steeringMode" => steeringMode
        ),
        function(entity)
            entity.data["nodes"] = [(Int(entity.data["x"]) + Int(entity.data["width"]) + 8, Int(entity.data["y"]))]
        end,
    ) for steeringMode in steeringModes
)

Ahorn.nodeLimits(entity::RailedMoveBlock) = 1, 1

Ahorn.minimumSize(entity::RailedMoveBlock) = 16, 16
Ahorn.resizable(entity::RailedMoveBlock) = true, true

Ahorn.editingOptions(entity::RailedMoveBlock) = Dict{String, Any}(
   "steeringMode" => steeringModes
)

const midColor = (4, 3, 23) ./ 255
const highlightColor = (59, 50, 101) ./ 255
const buttonColor = (71, 64, 112, 255) ./ 255

const button = "objects/moveBlock/button"

function Ahorn.selection(entity::RailedMoveBlock)
    x, y = Ahorn.position(entity)
    nx, ny = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return [Ahorn.Rectangle(x, y, width, height), Ahorn.Rectangle(nx + floor(Int, width / 2) - 5, ny + floor(Int, height / 2) - 5, 10, 10)]
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::RailedMoveBlock, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))
    nx, ny = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    cog = "objects/zipmover/cog"
    cx, cy = x + width / 2, y + height / 2
    cnx, cny = nx + width / 2, ny + height / 2

    length = sqrt((x - nx)^2 + (y - ny)^2)
    theta = atan(cny - cy, cnx - cx)

    Ahorn.Cairo.save(ctx)

    Ahorn.translate(ctx, cx, cy)
    Ahorn.rotate(ctx, theta)

    Ahorn.setSourceColor(ctx, midColor)
    Ahorn.set_antialias(ctx, 1)
    Ahorn.set_line_width(ctx, 1);

    # Offset for rounding errors
    Ahorn.move_to(ctx, 0, 4 + (theta <= 0))
    Ahorn.line_to(ctx, length, 4 + (theta <= 0))

    Ahorn.move_to(ctx, 0, -4 - (theta > 0))
    Ahorn.line_to(ctx, length, -4 - (theta > 0))

    Ahorn.stroke(ctx)

    Ahorn.Cairo.restore(ctx)

    Ahorn.drawSprite(ctx, cog, cnx, cny)

    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    steeringMode = get(entity.data, "steeringMode", "Horizontal")
    both = true
    horizontal = false
    vertical = false

    frame = "objects/CommunalHelper/railedMoveBlock/base_both"
    icon = "objects/CommunalHelper/railedMoveBlock/both"
    if steeringMode == "Horizontal"
        frame = "objects/moveBlock/base_v"
        icon = "objects/CommunalHelper/railedMoveBlock/h"
        both = false
        horizontal = true
    elseif steeringMode == "Vertical"
        frame = "objects/moveBlock/base_h"
        icon = "objects/CommunalHelper/railedMoveBlock/v"
        both = false
        vertical = true
    end

    iconSprite = Ahorn.getSprite(icon, "Gameplay")

    Ahorn.drawRectangle(ctx, x + 2, y + 2, width - 4, height - 4, highlightColor, highlightColor)
    Ahorn.drawRectangle(ctx, x + 8, y + 8, width - 16, height - 16, midColor)

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y, 8, 0, 8, 8)
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y + height - 8, 8, 16, 8, 8)

        if both || vertical
            Ahorn.drawImage(ctx, button, x + (i - 1) * 8, y - 2, 6, 0, 8, 6, tint=buttonColor)
        end
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, x, y + (i - 1) * 8, 0, 8, 8, 8)
        Ahorn.drawImage(ctx, frame, x + width - 8, y + (i - 1) * 8, 16, 8, 8, 8)

        if both || horizontal
            Ahorn.Cairo.save(ctx)

            Ahorn.rotate(ctx, -pi / 2)
            Ahorn.drawImage(ctx, button, -y + i * 8 - height - 8, x - 2, 6, 0, 8, 6, tint=buttonColor)
            Ahorn.scale(ctx, 1, -1)
            Ahorn.drawImage(ctx, button, -y + i * 8 - height - 8, -x - 2 - width, 6, 0, 8, 6, tint=buttonColor)

            Ahorn.Cairo.restore(ctx)
        end
    end

    Ahorn.drawImage(ctx, frame, x, y, 0, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, x + width - 8, y, 16, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, x, y + height - 8, 0, 16, 8, 8)
    Ahorn.drawImage(ctx, frame, x + width - 8, y + height - 8, 16, 16, 8, 8)

    if both || vertical
        Ahorn.Cairo.save(ctx)

        Ahorn.drawImage(ctx, button, x + 2, y - 2, 0, 0, 6, 6, tint=buttonColor)
        Ahorn.scale(ctx, -1, 1)
        Ahorn.drawImage(ctx, button, -x + 2 - width, y - 2, 0, 0, 6, 6, tint=buttonColor)

        Ahorn.Cairo.restore(ctx)
    end

    if both || horizontal
        Ahorn.Cairo.save(ctx)

        Ahorn.rotate(ctx, -pi / 2)
        Ahorn.drawImage(ctx, button, -y - height + 2, x - 2, 0, 0, 8, 6, tint=buttonColor)
        Ahorn.drawImage(ctx, button, -y - 10, x - 2, 14, 0, 8, 6, tint=buttonColor)
        Ahorn.scale(ctx, 1, -1)
        Ahorn.drawImage(ctx, button, -y - height + 2, -x - 2 -width, 0, 0, 8, 6, tint=buttonColor)
        Ahorn.drawImage(ctx, button, -y - 10, -x - 2 - width, 14, 0, 8, 6, tint=buttonColor)

        Ahorn.Cairo.restore(ctx)
    end

    Ahorn.drawRectangle(ctx, x + div(width - iconSprite.width, 2) + 1, y + div(height - iconSprite.height, 2) + 1, 8, 8, highlightColor, highlightColor)
    Ahorn.drawImage(ctx, iconSprite, x + div(width - iconSprite.width, 2), y + div(height - iconSprite.height, 2))
end

end


text:
import Cambrian.get_config
using YAML
using Dates

export ESIndividual, AbstractESIndividual, get_config

abstract type AbstractESIndividual <: Cambrian.Individual end

"FloatIndividual : example Individual class using a floating point genotype in [0, 1]"
struct ESIndividual <: AbstractESIndividual
    genes::Array{Float64}
    fitness::Array{Float64}
end

function ESIndividual(cfg::NamedTuple)
    ESIndividual(rand(cfg.n_genes), -Inf*ones(cfg.d_fitness))
end

function ESIndividual(st::String)
    dict = ind_parse(st)
    ESIndividual(Float64.(dict["genes"]), Float64.(dict["fitness"]))
end

## utils

function Cambrian.get_config(cfg_file::String; kwargs...)
    cfg = YAML.load_file(cfg_file)
    for (k, v) in kwargs
        cfg[String(k)] = v
    end
    # generate id, use date if no existing id
    if ~(:id in keys(cfg))
        cfg["id"] = replace(string(Dates.now()), r"[]\.:]" => "_")
    end
    get_config(cfg)
end


text:
# using Test
# using Distributions
# using StatsFuns

# using LRMoE

# # Poisson
# param0 = rand(Distributions.Uniform(0, 1), 1)[1]
# param1 = rand(Distributions.Uniform(0.5, 8), 1)[1]

# X = rand(Distributions.Uniform(-5, 5), 20000, 7)
# α = rand(Distributions.Uniform(-1, 1), 3, 7)
# α[2,:] .= 0

# model = [PoissonExpert(param1) PoissonExpert(2*param1) PoissonExpert(3*param1);
#         PoissonExpert(3*param1) PoissonExpert(0.5*param1) PoissonExpert(1.5*param1)]

# expos = rand(Distributions.Uniform(0.1, 5), 20000)
# expos_future = ceil.(expos) .- expos
# Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

# classes = LRMoE.predict_class_prior(X, α)
# classes.prob
# classes.max_prob_idx

# classes = LRMoE.predict_class_posterior(Y_sim, X, α, model, exposure_past = expos)
# classes.prob
# classes.max_prob_idx

# m1 = LRMoE.predict_mean_prior(X, α, model, exposure_future = expos_future)
# m2 = LRMoE.predict_mean_posterior(Y_sim, X, α, model, exposure_past = expos, exposure_future = expos_future)

# minimum(m2-m1, dims = 1)
# maximum(m2-m1, dims = 1)

# v1 = LRMoE.predict_var_prior(X, α, model, exposure_future = expos_future)
# v2 = LRMoE.predict_var_posterior(Y_sim, X, α, model, exposure_past = expos, exposure_future = expos_future)

# minimum(v2-v1, dims = 1)
# maximum(v2-v1, dims = 1)

# weights = classes.prob
# model_exp = LRMoE.exposurize_model(model, exposure = expos)
# means = mean.(model_exp)

# result = fill(NaN, size(weights)[1], 2)

# c_mean = mean.(model_exp)
# g_mean = predict_mean_prior(X, α, model, exposure = expos)
# var_c_mean = fill(NaN, size(X)[1], size(model)[2])
# for i in 1:size(X)[1]
#     var_c_mean[i,:] = (c_mean[:,:,i] .- g_mean[i,:] ).^2 * weights[i,:]
# end
# c_var = var.(model_exp)
# mean_c_var = fill(NaN, size(X)[1], size(model)[2])
# for i in 1:size(X)[1]
#     mean_c_var[i,:] = c_var[:,:,i] * weights[i,:]
# end

# c_mean[:,:,1]
# g_mean[1,:]
# (c_mean[:,:,1] .- g_mean[1,:]) * weights[1,:]

# c_var[:,:,1] * weights[1,:]

# model = [LogNormalExpert(param1, 1) LogNormalExpert(2*param1, 3) LogNormalExpert(3*param1, 2);
#         LogNormalExpert(3*param1, 0.5) LogNormalExpert(0.5*param1, 1.5) LogNormalExpert(1.5*param1, 1)]

# expos = rand(Distributions.Uniform(0.1, 5), 20000)
# expos_future = ceil.(expos) .- expos
# Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

# lim = rand(Uniform(0.1, 10), size(X)[1], 2)
# tmp1 = LRMoE.predict_limit_prior(X, α, model, lim, exposure_future = expos_future)
# tmp2 = LRMoE.predict_excess_prior(X, α, model, lim, exposure_future = expos_future)
# m = LRMoE.predict_mean_prior(X, α, model, exposure_future = expos_future)
# sum(m - tmp1 - tmp2)

# tmp1 = LRMoE.predict_limit_posterior(Y_sim, X, α, model, lim, exposure_past = expos, exposure_future = expos_future)
# tmp2 = LRMoE.predict_excess_posterior(Y_sim, X, α, model, lim, exposure_past = expos, exposure_future = expos_future)
# m = LRMoE.predict_mean_posterior(Y_sim, X, α, model, exposure_past = expos, exposure_future = expos_future)
# sum(m - tmp1 - tmp2)

# qs = rand(Uniform(0, 1), size(X)[1], 2)
# tmp1 = LRMoE.predict_VaRCTE_prior(X, α, model, qs, exposure_future = expos_future)
# tmp2 = LRMoE.predict_VaRCTE_posterior(Y_sim, X, α, model, qs, exposure_past = expos, exposure_future = expos_future)

# for i in 1:size(weights)[1]
#     result[i,:] = means[:,:,i] * weights[i,:]
# end

# hcat([means[:,:,i] * weights[i,:] for i in 1:size(weights)[1]]...)

# size(means)

# means[:,:,1]

# weights[1,:]

# gate = LogitGating(α, X)

# function _exact_to_full(Y)
#     result = fill(NaN, size(Y)[1], size(Y)[2]*4)
#     for j in 1:size(Y)[2]
#         result[:,4*(j-1)+1] .= 0.0
#         result[:,4*(j-1)+2] .= Y[:,j]
#         result[:,4*(j-1)+3] .= Y[:,j]
#         result[:,4*(j-1)+4] .= Inf
#     end
#     return result
# end

# Y_full = _exact_to_full(Y_sim)
# ll_np_list = LRMoE.loglik_np(Y_full, gate, model_exp)

# z_e_obs = LRMoE.EM_E_z_obs(ll_np_list.gate_expert_ll_comp, ll_np_list.gate_expert_ll)

# tmp = mean.(model_exp)

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
module Stick

export StickEvent, readstick, sticktask

using Dates

function _stick_input_device()
    try
        for devname in readdir("/sys/class/input")
            sysfname = joinpath("/sys/class/input",devname,"device","name")
            if startswith(devname, "event") && isfile(sysfname)
                if startswith(readline(sysfname),"Raspberry Pi Sense HAT Joystick")
                    return joinpath("/dev/input",devname)
                end
            end
        end
    catch e
    end
    error("Sense Hat not found.")
end

const STICK_INPUT_DEVICE = _stick_input_device()

const EV_KEY = Cushort(0x01)

@enum Direc UP=103 LEFT=105 RIGHT=106 DOWN=108 MIDDLE=28
@enum State RELEASE=0 PRESS=1 HOLD=2

"""
    StickEventRaw

The raw data format from the input device. This is based on the struct specification from https://github.com/RPi-Distro/python-sense-hat/blob/ec8c96118d5eb58488c4ff6091cd373117e6ce08/sense_hat/stick.py#L41.
"""
struct StickEventRaw
    tv_sec::Clong
    tv_usec::Clong
    ev_type::Cushort
    ev_direc::Cushort
    ev_state::Cuint
end
Base.read(io::IO, ::Type{StickEventRaw}) = read(io, Ref{StickEventRaw}())[]

"""
    StickEvent(timestamp::DateTime, direc::Direc, state::State)

The data from a joystick event. The type contains the following fields:

* `timestamp` : `DateTime` of when the event occurred.

* `direc` : the direcion of the event, either `UP`, `DOWN`, `LEFT`, `RIGHT`, `MIDDLE` (`MIDDLE`
  occurs when the button is pressed down).

* `state` : the state of the event (`PRESS`, `RELEASE`, `HOLD`).

"""
struct StickEvent
    timestamp::DateTime
    direc::Direc
    state::State
end

Base.convert(::Type{StickEvent}, r::StickEventRaw) =
    StickEvent(Dates.unix2datetime(r.tv_sec + r.tv_usec/1_000_000),
Direc(r.ev_direc), State(r.ev_state))

"""
    readstick()

This reads a single `StickEvent` from the joystick on the Sense HAT. This will block until an event occurs.

See also `sticktask()` which may be better for reading in a loop.
"""
function readstick()
    open(STICK_INPUT_DEVICE) do dev
        fddev = RawFD(fd(dev))
        while true
            wait(fddev, readable=true)
            s = read!(dev, Ref{StickEventRaw}())[]
            if s.ev_type == EV_KEY
                return convert(StickEvent,s)
            end
        end
    end
end

"""
    sticktask()
This is a `Task` that produces `StickEvent`s when the joystick on the Sense HAT is
manipulated. It will block until new `StickEvent`s occur.
A typical usage will be to create a new task which will call this asynchronously, e.g. the
following will call the function `f(::StickEvent)` for each event:
```
@schedule for e in sticktask()
    f(e)
end
```
"""
function sticktask()
    Base.depwarn("`sticktask()` is deprecated. Use `Channel` instead to create a buffer of `StickEvent`.", :sticktask)
    @task open(STICK_INPUT_DEVICE) do dev
        fddev = RawFD(fd(dev))
        while true
            wait(fddev, readable=true)
            s = read(dev, Ref{StickEventRaw}())[]
            if s.ev_type == EV_KEY
                produce(convert(StickEvent, s))
            end
        end
    end
end

end # module


text:
using POMDPLinter
using Documenter

makedocs(;
    modules=[POMDPLinter],
    authors="Zachary Sunberg <sunbergzach@gmail.com> and contributors",
    repo="https://github.com/JuliaPOMDP/POMDPLinter.jl/blob/{commit}{path}#L{line}",
    sitename="POMDPLinter.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaPOMDP.github.io/POMDPLinter.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Requirements" => [
            "Viewing Requirements" => "requirements.md",
            "Specifying Requirements" => "specifying_requirements.md"
        ]
    ],
)

deploydocs(;
    repo="github.com/JuliaPOMDP/POMDPLinter.jl",
)


text:
using Test
using LinearAlgebra

import ClimaCore:
    ClimaCore,
    slab,
    Spaces,
    Domains,
    Meshes,
    Geometry,
    Topologies,
    Spaces,
    Fields,
    Operators

using OrdinaryDiffEq: ODEProblem, solve, SSPRK33

import Logging
import TerminalLoggers
Logging.global_logger(TerminalLoggers.TerminalLogger())

# Nonhydrostatic gravity wave
# Reference: https://climate.ucdavis.edu/pubs/UJ2012JCP.pdf Section 5.4

const R = 6.37122e6 # radius
const grav = 9.8 # gravitational constant
const Ω = 0.0 # Earth rotation (radians / sec)
const R_d = 287.058 # R dry (gas constant / mol mass dry air)
const κ = 2 / 7 # kappa
const γ = 1.4 # heat capacity ratio
const cp_d = R_d / κ # heat capacity at constant pressure
const cv_d = cp_d - R_d # heat capacity at constant volume
const T_tri = 273.16 # triple point temperature
const N = 0.01 # Brunt-Vaisala frequency
const S = grav^2 / cp_d / N^2
const T_0 = 300 # isothermal atmospheric temperature
const Δθ = 10.0 # maximum potential temperature perturbation
const R_t = R / 3 # width of the perturbation
const L_z = 20.0e3 # vertial wave length of the perturbation
const p_0 = 1.0e5 # reference pressure
const λ_c = 180.0 # center longitude of the cosine bell
const ϕ_c = 0.0 # center latitude of the cosine bell

# set up function space
function sphere_3D(
    R = 6.37122e6,
    zlim = (0, 12.0e3),
    helem = 4,
    zelem = 12,
    npoly = 4,
)
    FT = Float64
    vertdomain = Domains.IntervalDomain(
        Geometry.ZPoint{FT}(zlim[1]),
        Geometry.ZPoint{FT}(zlim[2]);
        boundary_tags = (:bottom, :top),
    )
    vertmesh = Meshes.IntervalMesh(vertdomain, nelems = zelem)
    vert_center_space = Spaces.CenterFiniteDifferenceSpace(vertmesh)

    horzdomain = Domains.SphereDomain(R)
    horzmesh = Meshes.EquiangularCubedSphere(horzdomain, helem)
    horztopology = Topologies.Topology2D(horzmesh)
    quad = Spaces.Quadratures.GLL{npoly + 1}()
    horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)

    hv_center_space =
        Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vert_center_space)
    hv_face_space = Spaces.FaceExtrudedFiniteDifferenceSpace(hv_center_space)
    return (hv_center_space, hv_face_space)
end

Φ(z) = grav * z

function pressure(ρ, e, normuvw, z)
    I = e - Φ(z) - normuvw^2 / 2
    T = I / cv_d + T_tri
    return ρ * R_d * T
end

# set up 3D domain - spherical shell
hv_center_space, hv_face_space = sphere_3D(R, (0, 10.0e3), 5, 5, 4)

# initial conditions
coords = Fields.coordinate_field(hv_center_space)
face_coords = Fields.coordinate_field(hv_face_space)

r(λ, ϕ) = R * acos(sind(ϕ_c) * sind(ϕ) + cosd(ϕ_c) * cosd(ϕ) * cosd(λ - λ_c))

function initial_condition(ϕ, λ, z)
    rd = r(λ, ϕ)
    if rd < R_t
        s = 0.5 * (1 + cos(pi * rd / R_t))
    else
        s = 0.0
    end
    p = p_0 * (1 - S / T_0 + S / T_0 * exp(-N^2 * z / grav))^(cp_d / R_d)
    θ = T_0 * exp(N^2 * z / grav) + Δθ * s * sin(2 * pi * z / L_z)
    T = θ * (p / p_0)^κ
    ρ = p / R_d / T
    e = cv_d * (T - T_tri) + grav * z
    ρe = ρ * e

    return (ρ = ρ, ρe = ρe)
end

# Coriolis
const f =
    @. Geometry.Contravariant3Vector(Geometry.WVector(2 * Ω * sind(coords.lat)))

Yc = map(coord -> initial_condition(coord.lat, coord.long, coord.z), coords)
uₕ = map(_ -> Geometry.Covariant12Vector(0.0, 0.0), coords)
w = map(_ -> Geometry.Covariant3Vector(0.0), face_coords)
Y = Fields.FieldVector(Yc = Yc, uₕ = uₕ, w = w)

function rhs!(dY, Y, _, t)
    cρ = Y.Yc.ρ # scalar on centers
    fw = Y.w # Covariant3Vector on faces
    cuₕ = Y.uₕ # Covariant12Vector on centers
    cρe = Y.Yc.ρe # scalar on centers

    dρ = dY.Yc.ρ
    dw = dY.w
    duₕ = dY.uₕ
    dρe = dY.Yc.ρe

    # 0) update w at the bottom
    # fw = -g^31 cuₕ/ g^33

    hdiv = Operators.Divergence()
    hwdiv = Operators.Divergence()
    hgrad = Operators.Gradient()
    hwgrad = Operators.Gradient()
    hcurl = Operators.Curl()
    hwcurl = Operators.Curl()

    dρ .= 0 .* cρ

    ### HYPERVISCOSITY
    # 1) compute hyperviscosity coefficients

    χe = @. dρe = hwdiv(hgrad(cρe / cρ))
    χuₕ = @. duₕ =
        hwgrad(hdiv(cuₕ)) - Geometry.Covariant12Vector(
            hwcurl(Geometry.Covariant3Vector(hcurl(cuₕ))),
        )

    Spaces.weighted_dss!(dρe)
    Spaces.weighted_dss!(duₕ)

    κ₄ = 1.0e16 # m^4/s
    @. dρe = -κ₄ * hwdiv(cρ * hgrad(χe))
    @. duₕ =
        -κ₄ * (
            hwgrad(hdiv(χuₕ)) - Geometry.Covariant12Vector(
                hwcurl(Geometry.Covariant3Vector(hcurl(χuₕ))),
            )
        )

    # 1) Mass conservation
    If2c = Operators.InterpolateF2C()
    Ic2f = Operators.InterpolateC2F(
        bottom = Operators.Extrapolate(),
        top = Operators.Extrapolate(),
    )
    cw = If2c.(fw)
    cuvw = Geometry.Covariant123Vector.(cuₕ) .+ Geometry.Covariant123Vector.(cw)

    dw .= fw .* 0

    # 1.a) horizontal divergence
    dρ .-= hdiv.(cρ .* (cuvw))

    # 1.b) vertical divergence
    vdivf2c = Operators.DivergenceF2C(
        top = Operators.SetValue(Geometry.Contravariant3Vector(0.0)),
        bottom = Operators.SetValue(Geometry.Contravariant3Vector(0.0)),
    )
    # we want the total u³ at the boundary to be zero: we can either constrain
    # both to be zero, or allow one to be non-zero and set the other to be its
    # negation

    # explicit part
    dρ .-= vdivf2c.(Ic2f.(cρ .* cuₕ))
    # implicit part
    dρ .-= vdivf2c.(Ic2f.(cρ) .* fw)

    # 2) Momentum equation

    # curl term
    hcurl = Operators.Curl()
    # effectively a homogeneous Dirichlet condition on u₁ at the boundary
    vcurlc2f = Operators.CurlC2F(
        bottom = Operators.SetCurl(Geometry.Contravariant12Vector(0.0, 0.0)),
        top = Operators.SetCurl(Geometry.Contravariant12Vector(0.0, 0.0)),
    )
    cω³ = hcurl.(cuₕ) # Contravariant3Vector
    fω¹² = hcurl.(fw) # Contravariant12Vector
    fω¹² .+= vcurlc2f.(cuₕ) # Contravariant12Vector

    # cross product
    # convert to contravariant
    # these will need to be modified with topography
    fu¹² =
        Geometry.Contravariant12Vector.(
            Geometry.Covariant123Vector.(Ic2f.(cuₕ)),
        ) # Contravariant12Vector in 3D
    fu³ = Geometry.Contravariant3Vector.(Geometry.Covariant123Vector.(fw))
    @. dw -= fω¹² × fu¹² # Covariant3Vector on faces
    @. duₕ -= If2c(fω¹² × fu³)

    # Needed for 3D:
    @. duₕ -=
        (f + cω³) ×
        Geometry.Contravariant12Vector(Geometry.Covariant123Vector(cuₕ))

    ce = @. cρe / cρ
    cp = @. pressure(cρ, ce, norm(cuvw), coords.z)

    @. duₕ -= hgrad(cp) / cρ
    vgradc2f = Operators.GradientC2F(
        bottom = Operators.SetGradient(Geometry.Covariant3Vector(0.0)),
        top = Operators.SetGradient(Geometry.Covariant3Vector(0.0)),
    )
    @. dw -= vgradc2f(cp) / Ic2f(cρ)

    cE = @. (norm(cuvw)^2) / 2 + Φ(coords.z)
    @. duₕ -= hgrad(cE)
    @. dw -= vgradc2f(cE)

    # 3) total energy

    @. dρe -= hdiv(cuvw * (cρe + cp))
    @. dρe -= vdivf2c(fw * Ic2f(cρe + cp))
    @. dρe -= vdivf2c(Ic2f(cuₕ * (cρe + cp)))

    Spaces.weighted_dss!(dY.Yc)
    Spaces.weighted_dss!(dY.uₕ)
    Spaces.weighted_dss!(dY.w)

    return dY
end

dYdt = similar(Y)
rhs!(dYdt, Y, nothing, 0.0)

# run!
using OrdinaryDiffEq
# Solve the ODE
time_end = 600
dt = 3
prob = ODEProblem(rhs!, Y, (0.0, time_end))
sol = solve(
    prob,
    SSPRK33(),
    dt = dt,
    saveat = dt,
    progress = true,
    adaptive = false,
    progress_message = (dt, u, p, t) -> t,
)

@info "Solution L₂ norm at time t = 0: ", norm(Y.Yc.ρe)
@info "Solution L₂ norm at time t = $(time_end): ", norm(sol.u[end].Yc.ρe)

# TODO: visualization artifacts

# ENV["GKSwstype"] = "nul"
# using ClimaCorePlots, Plots
# Plots.GRBackend()
# dirname = "nonhydrostatic_gravity_wave"
# path = joinpath(@__DIR__, "output", dirname)
# mkpath(path)

# function linkfig(figpath, alt = "")
#     # buildkite-agent upload figpath
#     # link figure in logs if we are running on CI
#     if get(ENV, "BUILDKITE", "") == "true"
#         artifact_url = "artifact://$figpath"
#         print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n")
#     end
# end


text:
"""
    PBox <: SimpleProbabilitySet

Probability box defined by lower and upper distributions.
"""
struct PBox <: SimpleProbabilitySet
    dlower::Union{Bernoulli, DiscreteUniform, Categorical}
    dupper::Union{Bernoulli, DiscreteUniform, Categorical}
end

"""
    PBox(d::Distribution)

Degenerate probability box of only disribution 'd'
"""
PBox(d::Union{Bernoulli, DiscreteUniform, Categorical}) = PBox(d,d) # degenerate PBox

"""
    PBox()

Default probability box of Bernoulli(0.4) and Bernoulli(0.6)
"""
PBox() = PBox(Bernoulli(0.4), Bernoulli(0.6)) # default PBox

"""
    pdfs(pb::PBox)

Tuple of lower and upper pdfs
"""
pdfs(pb::PBox) = pdf.(pb.dlower, support(pb.dlower)), pdf.(pb.dupper, support(pb.dupper))

"""
    cdfs(pb::PBox)

Tuple of lower and upper cdfs
"""
cdfs(pb::PBox) = cdf.(pb.dlower, support(pb.dlower)), cdf.(pb.dupper, support(pb.dupper))

function pints(pb::PBox)
    cdflower, cdfupper = cdfs(pb)
    n = length(cdflower)
    plower = Array{Float64}(n)
    pupper = Array{Float64}(n)
    plower[1] = cdflower[1]
    pupper[1] = cdfupper[1]
    for i = 2:n
        plower[i] = max(cdflower[i] - cdfupper[i - 1], 0)
        pupper[i] = cdfupper[i] - cdflower[i - 1]
    end
    (plower, pupper)
end


text:
module SymbolDFGs

using LightGraphs
using SparseArrays
using DocStringExtensions

import ...DistributedFactorGraphs: AbstractDFG, DFGNode, AbstractParams, NoSolverParams, DFGVariable, DFGFactor
# import DFG functions to extend
import ...DistributedFactorGraphs:  setSolverParams,
                                    getFactor,
                                    setDescription,
                                    getLabelDict,
                                    addVariable!,
                                    getVariable,
                                    getAddHistory,
                                    addFactor!,
                                    getSolverParams,
                                    exists,
                                    isVariable,
                                    isFactor,
                                    getDescription,
                                    updateVariable!,
                                    updateFactor!,
                                    deleteVariable!,
                                    deleteFactor!,
                                    getVariables,
                                    getVariableIds,
                                    ls,
                                    getFactors,
                                    getFactorIds,
                                    lsf,
                                    isFullyConnected,
                                    hasOrphans,
                                    getNeighbors,
                                    getSubgraphAroundNode,
                                    getSubgraph,
                                    getIncidenceMatrix,
                                    getIncidenceMatrixSparse,
                                    _getDuplicatedEmptyDFG

include("SymbolFactorGraphs/SymbolFactorGraphs.jl")
using .SymbolFactorGraphs

# export SymbolEdge, is_directed, has_edge
# Imports
include("entities/SymbolDFG.jl")
include("services/SymbolDFG.jl")

# Exports
export SymbolDFG

end


text:
##### Beginning of file

function __init__()::Nothing
    _import_all_on_init()
    _print_welcome_message()
    return nothing
end

##### End of file


text:
module PkgCleanup

using Pkg, TOML, REPL.TerminalMenus, REPL.Terminals

"""
    logs_directory() -> String

Guess the path of the logs directory.
"""
logs_directory() =
    first(joinpath(d, "logs") for d in Pkg.depots() if isdir(joinpath(d, "logs")))

# Adapted from BinaryBuilder.jl/src/wizard/utils.jl
"""
    yn_prompt(question::AbstractString, default = :y;
              stdin::IO=Base.stdin,
              stdout::IO=Base.stdout,
              ) -> Bool

Perform a `[Y/n]` or `[y/N]` question loop, using `default` to choose between
the prompt styles, and looping until a proper response (e.g. `"y"`, `"yes"`,
`"n"` or `"no"`) is received.  Return `true` if "yes" is selected, `false` if
"no" is selected.
"""
function yn_prompt(question::AbstractString, default = :y;
                   stdin::IO=Base.stdin,
                   stdout::IO=Base.stdout,
                   )
    default in (:y, :n) || throw(ArgumentError("Default value must be either :y or :n"))
    ynstr = default == :y ? "[Y/n]" : "[y/N]"
    while true
        print(stdout, question, " ", ynstr, ": ")
        answer = lowercase(strip(readline(stdin)))
        if isempty(answer)
            return default === :y
        elseif answer == "y" || answer == "yes"
            return true
        elseif answer == "n" || answer == "no"
            return false
        else
            println(stdout, "Unrecognized answer. Answer `y` or `n`.")
        end
    end
end

function cleanup(path::String, what::String;
                 stdin::IO=Base.stdin,
                 stdout::IO=Base.stdout,
                 stderr::IO=Base.stderr,
                 )
    toml = TOML.parse(read(path, String))
    terminal = TTYTerminal("xterm", stdin, stdout, stderr)
    entries = sort!(collect(keys(toml)))
    select_menu = MultiSelectMenu(entries)
    select_menu.selected = Set(collect(eachindex(entries)))
    choices = request(terminal, "Select the $(what) to keep in $(path)", select_menu)
    selected_entries = entries[collect(choices)]
    filter!(e -> e.first ∈ selected_entries, toml)
    deleted = setdiff(entries, selected_entries)
    if !isempty(deleted)
        println(stdout)
        println(stdout, "You are going to remove the following entries in ", path, ":")
        for el in deleted
            println(stdout, "  * ", el)
        end
        if yn_prompt("Do you confirm your choice?"; stdin, stdout)
            # open(path, "w") do io
            #     TOML.print(io, toml; sorted=true)
            # end
        end
        println(stdout)
        println(stdout, """
    To actually cleanup packages/artifacts installations, remember to run

        using Pkg
        Pkg.gc()

    or, for a more in-depth cleaning,

        using Pkg
        using Dates
        Pkg.gc(; collect_delay=Dates.Day(0))
    """)
    end
    return nothing
end

"""
    PkgCleanup.artifacts()

Interactively select the `Artifacts.toml` files to keep in the list of active artifacts.
"""
artifacts(path=joinpath(logs_directory(), "artifact_usage.toml")) = cleanup(path, "Artifacts.toml")

"""
    PkgCleanup.manifests()

Interactively select the `Manifest.toml` files to keep in the list of active manifests.
"""
manifests(path=joinpath(logs_directory(), "manifest_usage.toml")) = cleanup(path, "Manifest.toml")

end


text:
import GLPK
# Need `"presolve" => GLPK.ON` for `detect_new_linearities`, see
# https://travis-ci.org/github/JuliaPolyhedra/Polyhedra.jl/jobs/691916637#L486
lp_solver = optimizer_with_attributes(GLPK.Optimizer, "presolve" => GLPK.GLP_ON)


text:
push!(LOAD_PATH,"D:\\repo\\complex-modeler\\src_PowSysMod")
using Documenter, PowSysMod

makedocs(
    #options
    format = :html,
    sitename = "Documentation PowSysMod",
    pages = [
        "Main functions" => "main_functions.md"
        "General structures" => "general_structures.md"
        "MatpowerInput" => "indexmatpower.md"
        "GOCInput" => "indexGOC.md"
        "PolynomialOptim" => "polynomialoptim_structures.md"
    ]
)

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
using StatsFuns
using Sobol
using Random

abstract type AbstractRandomSequenceGenerator end
mutable struct PseudoRandomRSG <: AbstractRandomSequenceGenerator
  rng::MersenneTwister
  dimension::Int
  values::Vector{Float64}
  weight::Float64
end

PseudoRandomRSG(seed::Int, dimension::Int = 1, weight::Float64 = 1.0) = PseudoRandomRSG(MersenneTwister(seed), dimension, zeros(dimension), weight)

mutable struct InverseCumulativeRSG <: AbstractRandomSequenceGenerator
  rng::MersenneTwister
  dimension::Int
  values::Vector{Float64}
  weight::Float64
end

InverseCumulativeRSG(seed::Int, dimension::Int = 1, weight::Float64 = 1.0) = InverseCumulativeRSG(MersenneTwister(seed), dimension, zeros(dimension), weight)

mutable struct SobolRSG <: AbstractRandomSequenceGenerator
  rng::Sobol.SobolSeq
  dimension::Int
  values::Vector{Float64}
  weight::Float64
end

SobolRSG(dimension::Int = 1, weight::Float64 = 1.0) = SobolRSG(SobolSeq(dimension), dimension, zeros(dimension), weight)

mutable struct SobolInverseCumulativeRSG <: AbstractRandomSequenceGenerator
  rng::Sobol.SobolSeq
  dimension::Int
  values::Vector{Float64}
  weight::Float64
end

SobolInverseCumulativeRSG(dimension::Int = 1, weight::Float64 = 1.0) = SobolInverseCumulativeRSG(SobolSeq(dimension), dimension, zeros(dimension), weight)

function next_sequence!(rsg::InverseCumulativeRSG)
  # we can probably use map here with the norminvcdf
  # like this: map!(norminvcdf, rsg.values, rand(rsg, length(rsg.values)))
  @inbounds @simd for i in eachindex(rsg.values)
    x = rand(rsg.rng) # get random number
    rsg.values[i] = norminvcdf(x)
  end

  return rsg.values, rsg.weight
end

function next_sequence!(rsg::SobolRSG)
  # we can probably use map here with the norminvcdf
  # like this: map!(norminvcdf, rsg.values, rand(rsg, length(rsg.values)))
  rsg.values = next(rsg.rng)
  # for i in eachindex(rsg.values)
  #   x = next(rsg.rng) # get random number
  #   rsg.values[i] = x
  # end

  return rsg.values, rsg.weight
end

function next_sequence!(rsg::SobolInverseCumulativeRSG)
  seq = next(rsg.rng)

  @simd for i in eachindex(rsg.values)
    @inbounds rsg.values[i] = norminvcdf(seq[i])
  end

  return rsg.values, rsg.weight
end

function next_sequence!(rsg::PseudoRandomRSG)
  rsg.values = rand(rsg.rng, rsg.dimension)
  return rsg.values, rsg.weight
end

last_sequence(rsg::AbstractRandomSequenceGenerator) = rsg.values, rsg.weight

function init_sequence_generator!(rsg::SobolRSG, dimension::Int)
  rsg.dimension = dimension
  rsg.values = zeros(dimension)
  rsg.rng = SobolSeq(dimension)

  return rsg
end

function init_sequence_generator!(rsg::AbstractRandomSequenceGenerator, dimension::Int)
  # MersenneTwister already init-ed
  rsg.dimension = dimension
  rsg.values = zeros(dimension)

  return rsg
end


text:
__precompile__()

module Gasp

using Base.Threads
enter_gc_safepoint() = ccall(:jl_gc_safe_enter, Int8, ())
leave_gc_safepoint(gs) = ccall(:jl_gc_safe_leave, Void, (Int8,), gs)

import Base.length, Base.get, Base.put!, Base.flush
export Garray, GarrayMemoryHandle, Dtree, Dlog,
       ngranks, grank, affinitize, message,
       sync, distribution, access, elemsize,
       initwork, getwork, runtree

const libgasp = joinpath(dirname(@__FILE__), "..", "deps", "gasp",
        "libgasp.$(Libdl.dlext)")

const ghandle = [C_NULL]
num_garrays = 0
exiting = false

function __init__()
    global ghandle
    ccall((:gasp_init, libgasp), Int64, (Cint, Ptr{Ptr{UInt8}}, Ptr{Void}),
          length(ARGS), ARGS, pointer(ghandle, 1))
    atexit() do
        global exiting
        exiting = true
    end
end

__init__()

function __shutdown__()
    ccall((:gasp_shutdown, libgasp), Void, (Ptr{Void},), ghandle[1])
end

@inline ngranks() = ccall((:gasp_nranks, libgasp), Int64, ())
@inline grank() = ccall((:gasp_rank, libgasp), Int64, ())+1
@inline sync() = ccall((:gasp_sync, libgasp), Void, ())
@inline cpu_pause() = ccall((:cpu_pause, libgasp), Void, ())
@inline rdtsc() = ccall((:rdtsc, libgasp), Culonglong, ())
@inline start_sde_tracing() = ccall((:start_sde_tracing, libgasp), Void, ())
@inline stop_sde_tracing() = ccall((:stop_sde_tracing, libgasp), Void, ())

show_affinity() = ccall(:puts, Cint, (Cstring,), "<$(threadid())> => $(ccall(:sched_getcpu, Cint, ()))")
show_affinity_mask() = ccall((:show_affinity_mask, Gasp.libgasp), Void, (Cint,), threadid())

function affinitize(avail_cores::Int,
                    avail_threads_per_core::Int,
                    ranks_per_node::Int;
                    use_threads_per_core::Int=1,
                    show::Bool=false)
    function set_thread_affinity()
        start = (((grank() - 1) % ranks_per_node) * div(nthreads(), use_threads_per_core))
        tid = threadid()
        offset, ht = divrem(tid - 1, use_threads_per_core)
        for i = 1:ht
            offset = offset + avail_cores
        end
        target = start + offset
        ccall((:set_affinity, libgasp), Cint, (Cint,), target)
        show && show_affinity()
    end
    ccall(:jl_threading_run, Void, (Any,), Core.svec(set_thread_affinity))
end

include("Garray.jl")
include("Dtree.jl")
include("Dlog.jl")

end # module


text:
# Last run: 7/1/2020

using ScientificTypes, DataFrames, Random

function foo()
    r = rand()
    r < 0.2 && return "AAA"
    r < 0.4 && return "BBB"
    r < 0.6 && return "CCC"
    r < 0.8 && return "DDD"
    return "EEE"
end

Random.seed!(5161)
v1 = [foo() for i in 1:1_000_000];
v2 = [foo() for i in 1:1_000_000];
df = DataFrame((x1=v1,x2=v2));

using BenchmarkTools
using CategoricalArrays

@btime categorical($v1); # 28.495 ms (1000046 allocations: 34.34 MiB)
@btime coerce($v1, Multiclass);
#this branch ==> 28.755 ms (1000046 allocations: 34.34 MiB)
#master      ==> 46.380 ms (1000080 allocations: 31.48 MiB)

function foo()
    r = rand()
    r < 0.2 && return "AAA"
    r < 0.4 && return "BBB"
    r < 0.6 && return "CCC"
    r < 0.8 && return "DDD"
    return missing
end

Random.seed!(51261)
v2 = [foo() for i in 1:1_000_000];

@btime categorical($v2); # 25.051 ms (799786 allocations: 28.22 MiB)
@btime coerce($v2, Union{Missing,Multiclass});
#this branch ==> 25.712 ms (799787 allocations: 28.22 MiB)
#master      ==> 139.152 ms (1599678 allocations: 51.69 MiB)

using Random, Distributions

n = 1_115_000
Random.seed!(51611)
df  = DataFrame(
        rs  = rand(["A", "B", "C", "D", "E", "F", "G"], n),
        es  = rand(["1","3","5","7","8","9"], n),
        inc = rand([collect(5000:15000); [missing for _ ∈ 1:50]], n),
        gr  = rand([randstring(6) for _ in 1:400], n),
        year = rand(2012:2017, n),
        lt = rand([12, 24, 36, 48, 60, 72], n),
        pb = Int.(round.(rand(LogNormal(8.07, 1.116), n))),
        cb = rand(LogNormal(7.88, 1.288), n),
        ci = rand(0.03:0.01:0.17, n),
        db = rand([repeat(collect(1945:2000),100); missing], n),
        ms = rand(["1", "2", "NA"], n),
        ho = BitArray(rand([true, false], n)),
        el = Int.(round.(rand(LogNormal(1.0, 0.72), n))),
        fl = rand(n),
        ft = rand(n),
        ea = BitArray(rand([true, false], n))
        );

@btime schema($df);
#this branch ==> 296.510 μs (140 allocations: 10.78 KiB)
#master      ==> 270.863 μs (140 allocations: 10.78 KiB)

@btime coerce($df, autotype($df, :few_to_finite));
#this branch ==> 979.098 ms (11152394 allocations: 636.62 MiB)
#master      ==> 1.415 s (12273798 allocations: 609.21 MiB)


text:
module NestedMaps

export nested_map, nested_map_recursive

using ArgCheck: @argcheck

"""
    nested_map_strict(f, A, kind...)

Implements `nested_map` given the element type and size, specified by the arguments
`kind...`, which is checked. See [`nested_map_kind`](@ref).

Internal function, not part of the API.
"""
function nested_map_strict(f, A, ::Type{<: NTuple{N, Any}}) where {N}
    @argcheck all(x -> x isa NTuple{N, Any}, A) "Incompatible elements."
    ntuple(n -> f(map(x -> x[n], A)), Val{N}())
end

function nested_map_strict(f, A, ::Type{<: NamedTuple{N}}) where {N}
    @argcheck all(x -> x isa NamedTuple{N}, A) "Incompatible elements."
    NamedTuple{N}(map(n -> f(map(x -> x[n], A)), N))
end

function nested_map_strict(f, A, ::Type{Array}, axs)
    @argcheck all(x -> axes(x) == axs, A) "Incompatible elements."
    [f(map(x -> x[i], A)) for i in CartesianIndices(axs)]
end

nested_map_strict(f, A) = f(A)

"""
    nested_map_kind(elt)

Provide a `kind` argument based on an element, usually the first one, for
[`nested_map_strict`](@ref).

`kind ≡ ()` is special, suggesting no transformation.

Internal function, not part of the API.
"""
nested_map_kind(elt::NTuple{N}) where N = (NTuple{N, Any}, )

nested_map_kind(elt::NamedTuple{N}) where N = (NamedTuple{N}, )

nested_map_kind(elt::T) where {T <: AbstractArray} = (Array, axes(elt))

nested_map_kind(elt) = ()

"""
    nested_map(f, A)

Return `B` such that `B[i] == f(map(x -> x[i], A))` for all valid indexes of `i`. The latter
is determined from the first element.
"""
nested_map(f, A) = nested_map_strict(f, A, nested_map_kind(first(A))...)

"""
    nested_map_recursive(f, A)

Same as `nested_map`, but applied recursively.
"""
function nested_map_recursive(f, A)
    kind = nested_map_kind(first(A))
    if kind ≡ ()
        f(A)
    else
        nested_map_strict(x -> nested_map_recursive(f, x), A, kind...)
    end
end

end # module


text:
# Computes the logarithmic mean: (aR-aL)/(LOG(aR)-LOG(aL)) = (aR-aL)/LOG(aR/aL)
# Problem: if aL~= aR, then 0/0, but should tend to --> 0.5*(aR+aL)
#
# introduce xi=aR/aL and f=(aR-aL)/(aR+aL) = (xi-1)/(xi+1)
# => xi=(1+f)/(1-f)
# => Log(xi) = log(1+f)-log(1-f), and for small f (f^2<1.0E-02) :
#
#    Log(xi) ~=     (f - 1/2 f^2 + 1/3 f^3 - 1/4 f^4 + 1/5 f^5 - 1/6 f^6 + 1/7 f^7)
#                  +(f + 1/2 f^2 + 1/3 f^3 + 1/4 f^4 + 1/5 f^5 + 1/6 f^6 + 1/7 f^7)
#             = 2*f*(1           + 1/3 f^2           + 1/5 f^4           + 1/7 f^6)
#  (aR-aL)/Log(xi) = (aR+aL)*f/(2*f*(1 + 1/3 f^2 + 1/5 f^4 + 1/7 f^6)) = (aR+aL)/(2 + 2/3 f^2 + 2/5 f^4 + 2/7 f^6)
#  (aR-aL)/Log(xi) = 0.5*(aR+aL)*(105/ (105+35 f^2+ 21 f^4 + 15 f^6)
function ln_mean(value1::Float64, value2::Float64)
  epsilon_f2 = 1.0e-4
  ratio = value2/value1
  # f2 = f^2
  f2=(ratio*(ratio-2.)+1.)/(ratio*(ratio+2.)+1.)
  if (f2<epsilon_f2)
    return (value1+value2)*52.5/(105.0 + f2*(35.0 + f2*(21.0 +f2*15.0)))
  else
    return (value2-value1)/log(ratio)
  end
end


text:
module DistributionsAD

using PDMats,
      ForwardDiff,
      Zygote,
      Tracker,
      LinearAlgebra,
      Distributions,
      Random,
      Combinatorics

using Tracker: TrackedReal
using LinearAlgebra: copytri!
using Distributions: AbstractMvLogNormal,
                     ContinuousMultivariateDistribution

import StatsFuns: logsumexp,
                  binomlogpdf,
                  nbinomlogpdf,
                  poislogpdf,
                  nbetalogpdf
import Distributions: MvNormal,
                      MvLogNormal,
                      poissonbinomial_pdf_fft,
                      logpdf,
                      quantile,
                      PoissonBinomial

export TuringDiagNormal,
       TuringMvNormal,
       TuringMvLogNormal,
       TuringPoissonBinomial

include("common.jl")
include("univariate.jl")
include("multivariate.jl")

end


text:
module GalacticOptimJL

using Reexport, GalacticOptim, GalacticOptim.SciMLBase
@reexport using Optim
decompose_trace(trace::Optim.OptimizationTrace) = last(trace)
decompose_trace(trace::Optim.OptimizationState) = trace

function __map_optimizer_args(prob::OptimizationProblem, opt::Union{Optim.AbstractOptimizer,Optim.Fminbox,Optim.SAMIN,Optim.ConstrainedOptimizer};
    callback=nothing,
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    kwargs...)
    if !isnothing(abstol)
        @warn "common abstol is currently not used by $(opt)"
    end

    mapped_args = (; extended_trace=true, kwargs...)

    if !isnothing(callback)
        mapped_args = (; mapped_args..., callback=callback)
    end

    if !isnothing(maxiters)
        mapped_args = (; mapped_args..., iterations=maxiters)
    end

    if !isnothing(maxtime)
        mapped_args = (; mapped_args..., time_limit=maxtime)
    end

    if !isnothing(reltol)
        mapped_args = (; mapped_args..., f_tol=reltol)
    end

    return Optim.Options(; mapped_args...)
end

function SciMLBase.__solve(prob::OptimizationProblem,
    opt::Optim.AbstractOptimizer,
    data=GalacticOptim.DEFAULT_DATA;
    kwargs...)
    if !isnothing(prob.lb) | !isnothing(prob.ub)
        if !(opt isa Union{Optim.Fminbox,Optim.SAMIN,Optim.AbstractConstrainedOptimizer})
            if opt isa Optim.ParticleSwarm
                opt = Optim.ParticleSwarm(; lower=prob.lb, upper=prob.ub, n_particles=opt.n_particles)
            elseif opt isa Optim.SimulatedAnnealing
                @warn "$(opt) can currently not be wrapped in Fminbox(). The lower and upper bounds thus will be ignored. Consider using a different optimizer or open an issue with Optim.jl"
            else
                opt = Optim.Fminbox(opt)
            end
        end
    end

    return ___solve(prob, opt, data; kwargs...)
end

function ___solve(prob::OptimizationProblem, opt::Optim.AbstractOptimizer,
    data=GalacticOptim.DEFAULT_DATA;
    callback=(args...) -> (false),
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    progress=false,
    kwargs...)

    local x, cur, state

    if data != GalacticOptim.DEFAULT_DATA
        maxiters = length(data)
    end

    cur, state = iterate(data)

    function _cb(trace)
        cb_call = opt isa Optim.NelderMead ? callback(decompose_trace(trace).metadata["centroid"], x...) : callback(decompose_trace(trace).metadata["x"], x...)
        if !(typeof(cb_call) <: Bool)
            error("The callback should return a boolean `halt` for whether to stop the optimization process.")
        end
        nx_itr = iterate(data, state)
        if isnothing(nx_itr)
            true
        else
            cur, state = nx_itr
            cb_call
        end
    end

    maxiters = GalacticOptim._check_and_convert_maxiters(maxiters)
    maxtime = GalacticOptim._check_and_convert_maxtime(maxtime)

    f = GalacticOptim.instantiate_function(prob.f, prob.u0, prob.f.adtype, prob.p)

    !(opt isa Optim.ZerothOrderOptimizer) && f.grad === nothing && error("Use OptimizationFunction to pass the derivatives or automatically generate them with one of the autodiff backends")

    _loss = function (θ)
        x = f.f(θ, prob.p, cur...)
        __x = first(x)
        return prob.sense === GalacticOptim.MaxSense ? -__x : __x
    end

    fg! = function (G, θ)
        if G !== nothing
            f.grad(G, θ, cur...)
            if prob.sense === GalacticOptim.MaxSense
                G .*= false
            end
        end
        return _loss(θ)
    end

    if opt isa Optim.KrylovTrustRegion
        hv = function (H, θ, v)
            f.hv(H, θ, v, cur...)
            if prob.sense === GalacticOptim.MaxSense
                H .*= false
            end
        end
        optim_f = Optim.TwiceDifferentiableHV(_loss, fg!, hv, prob.u0)
    else
        gg = function (G, θ)
            f.grad(G, θ, cur...)
            if prob.sense === GalacticOptim.MaxSense
                G .*= false
            end
        end

        hh = function (H, θ)
            f.hess(H, θ, cur...)
            if prob.sense === GalacticOptim.MaxSense
                H .*= false
            end
        end
        optim_f = Optim.TwiceDifferentiable(_loss, gg, fg!, hh, prob.u0)
    end

    opt_args = __map_optimizer_args(prob, opt, callback=_cb, maxiters=maxiters, maxtime=maxtime, abstol=abstol, reltol=reltol; kwargs...)

    t0 = time()
    opt_res = Optim.optimize(optim_f, prob.u0, opt, opt_args)
    t1 = time()
    opt_ret = Symbol(Optim.converged(opt_res))

    SciMLBase.build_solution(prob, opt, opt_res.minimizer, prob.sense === GalacticOptim.MaxSense ? -opt_res.minimum : opt_res.minimum; original=opt_res, retcode=opt_ret)
end

function ___solve(prob::OptimizationProblem, opt::Union{Optim.Fminbox,Optim.SAMIN},
    data=GalacticOptim.DEFAULT_DATA;
    callback=(args...) -> (false),
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    progress=false,
    kwargs...)

    local x, cur, state

    if data != GalacticOptim.DEFAULT_DATA
        maxiters = length(data)
    end

    cur, state = iterate(data)

    function _cb(trace)
        cb_call = !(opt isa Optim.SAMIN) && opt.method == Optim.NelderMead() ? callback(decompose_trace(trace).metadata["centroid"], x...) : callback(decompose_trace(trace).metadata["x"], x...)
        if !(typeof(cb_call) <: Bool)
            error("The callback should return a boolean `halt` for whether to stop the optimization process.")
        end
        nx_itr = iterate(data, state)
        if isnothing(nx_itr)
            true
        else
            cur, state = nx_itr
            cb_call
        end
    end

    maxiters = GalacticOptim._check_and_convert_maxiters(maxiters)
    maxtime = GalacticOptim._check_and_convert_maxtime(maxtime)

    f = GalacticOptim.instantiate_function(prob.f, prob.u0, prob.f.adtype, prob.p)

    !(opt isa Optim.ZerothOrderOptimizer) && f.grad === nothing && error("Use OptimizationFunction to pass the derivatives or automatically generate them with one of the autodiff backends")

    _loss = function (θ)
        x = f.f(θ, prob.p, cur...)
        __x = first(x)
        return prob.sense === GalacticOptim.MaxSense ? -__x : __x
    end
    fg! = function (G, θ)
        if G !== nothing
            f.grad(G, θ, cur...)
            if prob.sense === GalacticOptim.MaxSense
                G .*= false
            end
        end
        return _loss(θ)
    end

    gg = function (G, θ)
        f.grad(G, θ, cur...)
        if prob.sense === GalacticOptim.MaxSense
            G .*= false
        end
    end
    optim_f = Optim.OnceDifferentiable(_loss, gg, fg!, prob.u0)

    opt_args = __map_optimizer_args(prob, opt, callback=_cb, maxiters=maxiters, maxtime=maxtime, abstol=abstol, reltol=reltol; kwargs...)

    t0 = time()
    opt_res = Optim.optimize(optim_f, prob.lb, prob.ub, prob.u0, opt, opt_args)
    t1 = time()
    opt_ret = Symbol(Optim.converged(opt_res))

    SciMLBase.build_solution(prob, opt, opt_res.minimizer, opt_res.minimum; original=opt_res, retcode=opt_ret)
end

function ___solve(prob::OptimizationProblem, opt::Optim.ConstrainedOptimizer,
    data=GalacticOptim.DEFAULT_DATA;
    callback=(args...) -> (false),
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    progress=false,
    kwargs...)

    local x, cur, state

    if data != GalacticOptim.DEFAULT_DATA
        maxiters = length(data)
    end

    cur, state = iterate(data)

    function _cb(trace)
        cb_call = callback(decompose_trace(trace).metadata["x"], x...)
        if !(typeof(cb_call) <: Bool)
            error("The callback should return a boolean `halt` for whether to stop the optimization process.")
        end
        nx_itr = iterate(data, state)
        if isnothing(nx_itr)
            true
        else
            cur, state = nx_itr
            cb_call
        end
    end

    maxiters = GalacticOptim._check_and_convert_maxiters(maxiters)
    maxtime = GalacticOptim._check_and_convert_maxtime(maxtime)

    f = GalacticOptim.instantiate_function(prob.f, prob.u0, prob.f.adtype, prob.p, prob.ucons === nothing ? 0 : length(prob.ucons))

    f.cons_j === nothing && error("This optimizer requires derivative definitions for nonlinear constraints. If the problem does not have nonlinear constraints, choose a different optimizer. Otherwise define the derivative for cons using OptimizationFunction either directly or automatically generate them with one of the autodiff backends")

    _loss = function (θ)
        x = f.f(θ, prob.p, cur...)
        __x = first(x)
        return prob.sense === GalacticOptim.MaxSense ? -__x : __x
    end
    fg! = function (G, θ)
        if G !== nothing
            f.grad(G, θ, cur...)
            if prob.sense === GalacticOptim.MaxSense
                G .*= false
            end
        end
        return _loss(θ)
    end
    gg = function (G, θ)
        f.grad(G, θ, cur...)
        if prob.sense === GalacticOptim.MaxSense
            G .*= false
        end
    end

    hh = function (H, θ)
        f.hess(H, θ, cur...)
        if prob.sense === GalacticOptim.MaxSense
            H .*= false
        end
    end
    optim_f = Optim.TwiceDifferentiable(_loss, gg, fg!, hh, prob.u0)

    cons! = (res, θ) -> res .= f.cons(θ)

    cons_j! = function (J, x)
        f.cons_j(J, x)
    end

    cons_hl! = function (h, θ, λ)
        res = [similar(h) for i in 1:length(λ)]
        f.cons_h(res, θ)
        for i in 1:length(λ)
            h .+= λ[i] * res[i]
        end
    end

    lb = prob.lb === nothing ? [] : prob.lb
    ub = prob.ub === nothing ? [] : prob.ub
    optim_fc = Optim.TwiceDifferentiableConstraints(cons!, cons_j!, cons_hl!, lb, ub, prob.lcons, prob.ucons)

    opt_args = __map_optimizer_args(prob, opt, callback=_cb, maxiters=maxiters, maxtime=maxtime, abstol=abstol, reltol=reltol; kwargs...)

    t0 = time()
    opt_res = Optim.optimize(optim_f, optim_fc, prob.u0, opt, opt_args)
    t1 = time()
    opt_ret = Symbol(Optim.converged(opt_res))

    SciMLBase.build_solution(prob, opt, opt_res.minimizer, opt_res.minimum; original=opt_res, retcode=opt_ret)
end

end


text:
# ---
# title: 557. Reverse Words in a String III
# id: problem557
# author: zhwang
# date: 2022-02-15
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/reverse-words-in-a-string-iii/description/>
# hidden: true
# ---
#
# Given a string, you need to reverse the order of characters in each word
# within a sentence while still preserving whitespace and initial word order.
#
# **Example 1:**
#
#
#
#     Input: "Let's take LeetCode contest"
#     Output: "s'teL ekat edoCteeL tsetnoc"
#
#
# **Note:** In the string, each word is separated by single space and there will
# not be any extra space in the string.
#
#
## @lc code=start
using LeetCode
reverse_words_557(s::String) = join([reverse(i) for i in split(s)], " ")

## @lc code=end

## @lc test=start
@testset "557.reverse-words-in-a-string-iii.jl" begin
    @test reverse_words_557("Let's take LeetCode contest") == "s'teL ekat edoCteeL tsetnoc"
    @test reverse_words_557("God Ding") =="doG gniD"
end
## @lc test=end

## @lc info=start
# link: [solution to 557](https://leetcode-cn.com/problems/reverse-words-in-a-string-iii/solution/557-fan-zhuan-zi-fu-chuan-zhong-de-dan-c-wibs/)
## @lc info=end


text:
module Enzyme

export autodiff
export Const, Active, Duplicated

using Cassette

abstract type Annotation{T} end
struct Const{T} <: Annotation{T}
    val::T
end
struct Active{T<:AbstractFloat} <: Annotation{T}
    val::T
end
Active(i::Integer) = Active(float(i))
struct Duplicated{T} <: Annotation{T}
    val::T
    dval::T
end
struct DuplicatedNoNeed{T} <: Annotation{T}
    val::T
    dval::T
end

Base.eltype(::Type{<:Annotation{T}}) where T = T

import LLVM

include("api.jl")
include("typeanalysis.jl")
include("typetree.jl")
include("utils.jl")
include("compiler.jl")

annotate() = ()
annotate(arg::Annotation, args...) = (arg, annotate(args...)...)
annotate(arg, args...) = (Const(arg), annotate(args...)...)

prepare_cc() = ()
prepare_cc(arg::Duplicated, args...) = (arg.val, arg.dval, prepare_cc(args...)...)
prepare_cc(arg::DuplicatedNoNeed, args...) = (arg.val, arg.dval, prepare_cc(args...)...)
prepare_cc(arg::Annotation, args...) = (arg.val, prepare_cc(args...)...)

function autodiff(f, args...)
    args′ =  annotate(args...)
    thunk = Compiler.thunk(f, Tuple{map(Core.Typeof, args′)...})
    thunk(prepare_cc(args′...)...)
end

import .Compiler: EnzymeCtx
# Ops that have intrinsics
for op in (sin, cos, tan, exp)
    for (T, suffix) in ((Float32, "f32"), (Float64, "f64"))
        llvmf = "llvm.$(nameof(op)).$suffix"
        @eval begin
            @inline function Cassette.overdub(::EnzymeCtx, ::typeof($op), x::$T)
                ccall($llvmf, llvmcall, $T, ($T,), x)
            end
        end
    end
end

for op in (copysign,)
    for (T, suffix) in ((Float32, "f32"), (Float64, "f64"))
        llvmf = "llvm.$(nameof(op)).$suffix"
        @eval begin
            @inline function Cassette.overdub(::EnzymeCtx, ::typeof($op), x::$T, y::$T)
                ccall($llvmf, llvmcall, $T, ($T, $T), x, y)
            end
        end
    end
end

for op in (asin,tanh)
    for (T, llvm_t) in ((Float32, "float"), (Float64, "double"))
        decl = "declare double @$(nameof(op))($llvm_t)"
        func = """
               %val = call $llvm_t @$op($llvm_t %0)
               ret $llvm_t %val
               """
       @eval begin
            @inline function Cassette.overdub(::EnzymeCtx, ::typeof($op), x::$T)
                Base.llvmcall(($decl,$func), $T, Tuple{$T}, x)
            end
        end
    end
end

@inline function pack(args...)
    ntuple(Val(length(args))) do i
        Base.@_inline_meta
        arg = args[i]
        @assert arg isa AbstractFloat
        return Duplicated(Ref(args[i]), Ref(zero(args[i])))
    end
end

@inline unpack() = ()
@inline unpack(arg) = (arg[],)
@inline unpack(arg, args...) = (arg[], unpack(args...)...)

@inline ∇unpack() = ()
@inline ∇unpack(arg::Duplicated) = (arg.dval[],)
@inline ∇unpack(arg::Duplicated, args...) = (arg.dval[], ∇unpack(args...)...)

function gradient(f, args...)
    ∇args = pack(args...)
    f′ = function (args...)
        Base.@_inline_meta
        f(unpack(args...)...)
    end
    autodiff(f′, ∇args...)
    return ∇unpack(∇args...)
end

function pullback(f, args...)
    return (c) -> begin
        ∇vals = gradient(f, args...)
        return ntuple(Val(length(∇vals))) do i
            Base.@_inline_meta
            return c*∇vals[i]
        end
    end
end

# WIP
# @inline Cassette.overdub(::EnzymeCtx, ::typeof(asin), x::Float64) = ccall(:asin, Float64, (Float64,), x)
# @inline Cassette.overdub(::EnzymeCtx, ::typeof(asin), x::Float32) = ccall(:asin, Float32, (Float32,), x)
end # module


text:
mutable struct HTTPAdapter <: HTTPAdapterBase
    server::String
    output_channel::OutputChannel
    sessionKey::String
    closed::Bool
end

HTTPAdapter() = HTTPAdapter("", OutputChannel(DEFAULT_BUFFER_SIZE), "", false)

function loop(adp::HTTPAdapter, server, qq, verifyKey; poll_interval = 1, fetch_count = 8)
    adp.server = server
    # Verify
    resp = post_restful("http://$server/verify", (), (; verifyKey))
    adp.sessionKey = resp[:session]
    # Bind
    resp = post_restful("http://$server/bind", (), (; adp.sessionKey, qq))
    resp[:code] == 0 || throw(AdapterConectionFailed(resp[:code], resp[:msg]))

    put!(adp.output_channel, resp)
    loop_body(adp, qq, poll_interval, fetch_count, Commands.response_type, SESSION_KEY_IN_HEADERS)
end


text: logis(x) = 1 ./(1 .+ exp(-x))


text:
@testset "Poisson 2D on rectangular grid" begin

    alpha = 0.05
    kx    = 0.5
    ky    = 1.0

    mesh = Mesh( 0, 2π/kx, 64, 0, 2π/ky, 128)

    fields = MeshFields( mesh )
    solutions = MeshFields( mesh)

    x = LinRange(mesh.xmin, mesh.xmax, mesh.nx+1)
    y = LinRange(mesh.ymin, mesh.ymax, mesh.ny+1)

    fields.ρ  .= - 8 .* sin.(2 .* x) .* cos.(2 .* y')

    solutions.e[1,:,:] .=   2 .* cos.(2 .* x) .* cos.(2 .* y')
    solutions.e[2,:,:] .= - 2 .* sin.(2 .* x) .* sin.(2 .* y')

    poisson! = Poisson( mesh )

    poisson!( fields )

    err = errors( fields, solutions )

    println(" error : $err ")

    @test err ≈ 0.0 atol = 1e-14

    fields.ρ .= - 4 * ( sin.(2*x) .+ cos.(2*y') )

    poisson!( fields )

    nx, ny = mesh.nx, mesh.ny

    for j in 1:ny+1, i in 1:nx+1
        solutions.e[1,i,j] =   2*cos(2*x[i])
        solutions.e[2,i,j] = - 2*sin(2*y[j])
    end

    gnuplot("test2.dat", fields)
    gnuplot("solu2.dat", solutions)

    println(" error : $err ")
    err = errors( fields, solutions )

    @test err ≈ 0.0 atol = 1e-14

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to conversion from direction cosine matrices to Euler angles.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/conversions/dcm_to_angle.jl
# =======================================

# Functions: dcm_to_angle
# -----------------------

@testset "DCM => Euler angles (Float64)" begin
    T = Float64

    # The conversion is performed by creating DCMs using the tested function
    # `create_rotation_matrix`, and then converting to Euler angles.

    # The test set is formed of three rotations angles, the rotation sequence,
    # and the comparison mode. The latter is used for the cases with
    # singularities. In those cases, we need to know if the angles a₁ and a₃
    # must be summed or subtracted due to the singularity.
    testset = [
        # ZYX
        (_rand_ang(), _rand_ang2(), _rand_ang(), :Z, :Y, :X, :none)
        (_rand_ang(), +π / 2      , _rand_ang(), :Z, :Y, :X, :sub)
        (_rand_ang(), -π / 2      , _rand_ang(), :Z, :Y, :X, :sum)
        # XYX
        (_rand_ang(), _rand_ang3(), _rand_ang(), :X, :Y, :X, :none)
        (_rand_ang(), 0           , _rand_ang(), :X, :Y, :X, :sum)
        (_rand_ang(), +π          , _rand_ang(), :X, :Y, :X, :sub)
        (_rand_ang(), -π          , _rand_ang(), :X, :Y, :X, :sub)
        # XYZ
        (_rand_ang(), _rand_ang2(), _rand_ang(), :X, :Y, :Z, :none)
        (_rand_ang(), +π / 2      , _rand_ang(), :X, :Y, :Z, :sum)
        (_rand_ang(), -π / 2      , _rand_ang(), :X, :Y, :Z, :sub)
        # XZX
        (_rand_ang(), _rand_ang3(), _rand_ang(), :X, :Z, :X, :none)
        (_rand_ang(), 0           , _rand_ang(), :X, :Z, :X, :sum)
        (_rand_ang(), +π          , _rand_ang(), :X, :Z, :X, :sub)
        (_rand_ang(), -π          , _rand_ang(), :X, :Z, :X, :sub)
        # XZY
        (_rand_ang(), _rand_ang2(), _rand_ang(), :X, :Z, :Y, :none)
        (_rand_ang(), +π / 2      , _rand_ang(), :X, :Z, :Y, :sub)
        (_rand_ang(), -π / 2      , _rand_ang(), :X, :Z, :Y, :sum)
        # YXY
        (_rand_ang(), _rand_ang3(), _rand_ang(), :Y, :X, :Y, :none)
        (_rand_ang(), 0           , _rand_ang(), :Y, :X, :Y, :sum)
        (_rand_ang(), +π          , _rand_ang(), :Y, :X, :Y, :sub)
        (_rand_ang(), -π          , _rand_ang(), :Y, :X, :Y, :sub)
        # YXZ
        (_rand_ang(), _rand_ang2(), _rand_ang(), :Y, :X, :Z, :none)
        (_rand_ang(), +π / 2      , _rand_ang(), :Y, :X, :Z, :sub)
        (_rand_ang(), -π / 2      , _rand_ang(), :Y, :X, :Z, :sum)
        # YZX
        (_rand_ang(), _rand_ang2(), _rand_ang(), :Y, :Z, :X, :none)
        (_rand_ang(), +π / 2      , _rand_ang(), :Y, :Z, :X, :sum)
        (_rand_ang(), -π / 2      , _rand_ang(), :Y, :Z, :X, :sub)
        # YZY
        (_rand_ang(), _rand_ang3(), _rand_ang(), :Y, :Z, :Y, :none)
        (_rand_ang(), 0           , _rand_ang(), :Y, :Z, :Y, :sum)
        (_rand_ang(), +π          , _rand_ang(), :Y, :Z, :Y, :sub)
        (_rand_ang(), -π          , _rand_ang(), :Y, :Z, :Y, :sub)
        # ZXY
        (_rand_ang(), _rand_ang2(), _rand_ang(), :Z, :X, :Y, :none)
        (_rand_ang(), +π / 2      , _rand_ang(), :Z, :X, :Y, :sum)
        (_rand_ang(), -π / 2      , _rand_ang(), :Z, :X, :Y, :sub)
        # ZXZ
        (_rand_ang(), _rand_ang3(), _rand_ang(), :Z, :X, :Z, :none)
        (_rand_ang(), 0           , _rand_ang(), :Z, :X, :Z, :sum)
        (_rand_ang(), +π          , _rand_ang(), :Z, :X, :Z, :sub)
        (_rand_ang(), -π          , _rand_ang(), :Z, :X, :Z, :sub)
        # ZYZ
        (_rand_ang(), _rand_ang3(), _rand_ang(), :Z, :Y, :Z, :none)
        (_rand_ang(), 0           , _rand_ang(), :Z, :Y, :Z, :sum)
        (_rand_ang(), +π          , _rand_ang(), :Z, :Y, :Z, :sub)
        (_rand_ang(), -π          , _rand_ang(), :Z, :Y, :Z, :sub)
    ]

    for test in testset
        # Unpack values in tuple.
        a₁, a₂, a₃, r₁, r₂, r₃, c = test

        D = create_rotation_matrix(a₃, r₃) *
            create_rotation_matrix(a₂, r₂) *
            create_rotation_matrix(a₁, r₁)

        rot_seq = Symbol(string(r₁) * string(r₂) * string(r₃))
        ea = dcm_to_angle(D, rot_seq)
        @test eltype(ea) === T

        # Check the rotation sequence.
        @test ea.rot_seq == rot_seq

        # Compare the representations considering the singularities.
        if c == :none
            @test ea.a1 ≈ a₁
            @test ea.a2 ≈ a₂
            @test ea.a3 ≈ a₃
        elseif c == :sum
            # Treat the singularity when a₂ is ±π.
            abs(a₂) ≈ π && (a₂ = sign(ea.a2) * π)

            @test ea.a1 ≈ _norm_ang(a₁ + a₃)
            @test ea.a2 ≈ a₂
            @test ea.a3 ≈ 0
        elseif c == :sub
            # Treat the singularity when a₂ is ±π.
            abs(a₂) ≈ π && (a₂ = sign(ea.a2) * π)

            @test ea.a1 ≈ _norm_ang(a₁ - a₃)
            @test ea.a2 ≈ a₂
            @test ea.a3 ≈ 0
        end
    end
end

@testset "DCM => Euler angles (Float32)" begin
    T = Float32

    # The conversion is performed by creating DCMs using the tested function
    # `create_rotation_matrix`, and then converting to Euler angles.

    # The test set is formed of three rotations angles, the rotation sequence,
    # and the comparison mode. The latter is used for the cases with
    # singularities. In those cases, we need to know if the angles a₁ and a₃
    # must be summed or subtracted due to the singularity.
    testset = [
        # ZYX
        (_rand_ang(T), _rand_ang2(T), _rand_ang(T), :Z, :Y, :X, :none)
        (_rand_ang(T), T(+π / 2)    , _rand_ang(T), :Z, :Y, :X, :sub)
        (_rand_ang(T), T(-π / 2)    , _rand_ang(T), :Z, :Y, :X, :sum)
        # XYX
        (_rand_ang(T), _rand_ang3(T), _rand_ang(T), :X, :Y, :X, :none)
        (_rand_ang(T), T(0)         , _rand_ang(T), :X, :Y, :X, :sum)
        (_rand_ang(T), T(+π)        , _rand_ang(T), :X, :Y, :X, :sub)
        (_rand_ang(T), T(-π)        , _rand_ang(T), :X, :Y, :X, :sub)
        # XYZ
        (_rand_ang(T), _rand_ang2(T), _rand_ang(T), :X, :Y, :Z, :none)
        (_rand_ang(T), T(+π / 2)    , _rand_ang(T), :X, :Y, :Z, :sum)
        (_rand_ang(T), T(-π / 2)    , _rand_ang(T), :X, :Y, :Z, :sub)
        # XZX
        (_rand_ang(T), _rand_ang3(T), _rand_ang(T), :X, :Z, :X, :none)
        (_rand_ang(T), T(0)         , _rand_ang(T), :X, :Z, :X, :sum)
        (_rand_ang(T), T(+π)        , _rand_ang(T), :X, :Z, :X, :sub)
        (_rand_ang(T), T(-π)        , _rand_ang(T), :X, :Z, :X, :sub)
        # YXY
        (_rand_ang(T), _rand_ang3(T), _rand_ang(T), :Y, :X, :Y, :none)
        (_rand_ang(T), T(0)         , _rand_ang(T), :Y, :X, :Y, :sum)
        (_rand_ang(T), T(+π)        , _rand_ang(T), :Y, :X, :Y, :sub)
        (_rand_ang(T), T(-π)        , _rand_ang(T), :Y, :X, :Y, :sub)
        # YXZ
        (_rand_ang(T), _rand_ang2(T), _rand_ang(T), :Y, :X, :Z, :none)
        (_rand_ang(T), T(+π / 2)    , _rand_ang(T), :Y, :X, :Z, :sub)
        (_rand_ang(T), T(-π / 2)    , _rand_ang(T), :Y, :X, :Z, :sum)
        # YZY
        (_rand_ang(T), _rand_ang3(T), _rand_ang(T), :Y, :Z, :Y, :none)
        (_rand_ang(T), T(0)         , _rand_ang(T), :Y, :Z, :Y, :sum)
        (_rand_ang(T), T(+π)        , _rand_ang(T), :Y, :Z, :Y, :sub)
        (_rand_ang(T), T(-π)        , _rand_ang(T), :Y, :Z, :Y, :sub)
        # ZXY
        (_rand_ang(T), _rand_ang2(T), _rand_ang(T), :Z, :X, :Y, :none)
        (_rand_ang(T), T(+π / 2)    , _rand_ang(T), :Z, :X, :Y, :sum)
        (_rand_ang(T), T(-π / 2)    , _rand_ang(T), :Z, :X, :Y, :sub)
        # ZXZ
        (_rand_ang(T), _rand_ang3(T), _rand_ang(T), :Z, :X, :Z, :none)
        (_rand_ang(T), T(0)         , _rand_ang(T), :Z, :X, :Z, :sum)
        (_rand_ang(T), T(+π)        , _rand_ang(T), :Z, :X, :Z, :sub)
        (_rand_ang(T), T(-π)        , _rand_ang(T), :Z, :X, :Z, :sub)
        # ZYZ
        (_rand_ang(T), _rand_ang3(T), _rand_ang(T), :Z, :Y, :Z, :none)
        (_rand_ang(T), T(0)         , _rand_ang(T), :Z, :Y, :Z, :sum)
        (_rand_ang(T), T(+π)        , _rand_ang(T), :Z, :Y, :Z, :sub)
        (_rand_ang(T), T(-π)        , _rand_ang(T), :Z, :Y, :Z, :sub)
    ]

    for test in testset
        # Unpack values in tuple.
        a₁, a₂, a₃, r₁, r₂, r₃, c = test

        D = create_rotation_matrix(a₃, r₃) *
            create_rotation_matrix(a₂, r₂) *
            create_rotation_matrix(a₁, r₁)

        rot_seq = Symbol(string(r₁) * string(r₂) * string(r₃))
        ea = dcm_to_angle(D, rot_seq)
        @test eltype(ea) === T

        # Check the rotation sequence.
        @test ea.rot_seq == rot_seq

        # Compare the representations considering the singularities.
        if c == :none
            @test ea.a1 ≈ a₁
            @test ea.a2 ≈ a₂
            @test ea.a3 ≈ a₃
        elseif c == :sum
            # Treat the singularity when a₂ is ±π.
            abs(a₂) ≈ π && (a₂ = sign(ea.a2) * π)

            @test ea.a1 ≈ _norm_ang(a₁ + a₃)
            @test ea.a2 ≈ a₂
            @test ea.a3 ≈ 0
        elseif c == :sub
            # Treat the singularity when a₂ is ±π.
            abs(a₂) ≈ π && (a₂ = sign(ea.a2) * π)

            @test ea.a1 ≈ _norm_ang(a₁ - a₃)
            @test ea.a2 ≈ a₂
            @test ea.a3 ≈ 0
        end
    end
end

@testset "DCM => Euler angles (errors)" begin
    D = DCM(I)
    @test_throws ArgumentError dcm_to_angle(D, :XXY)
end

@testset "DCM => Euler angles (internal functions)" begin
    # Float64
    # ==========================================================================
    T = Float64

    a = ReferenceFrameRotations._mod_acos(1.0 + eps(T))
    @test a isa T
    @test a == 0

    a = ReferenceFrameRotations._mod_acos(-1.0 - eps(T))
    @test a isa T
    @test a == T(π)

    a = ReferenceFrameRotations._mod_asin(1.0 + eps(T))
    @test a isa T
    @test a == T(π / 2)

    a = ReferenceFrameRotations._mod_asin(-1.0 - eps(T))
    @test a isa T
    @test a == T(-π / 2)

    # Float32
    # ==========================================================================
    T = Float32

    a = ReferenceFrameRotations._mod_acos(1.0f0 + eps(T))
    @test a isa T
    @test a == 0

    a = ReferenceFrameRotations._mod_acos(-1.0f0 - eps(T))
    @test a isa T
    @test a == T(π)

    a = ReferenceFrameRotations._mod_asin(1.0f0 + eps(T))
    @test a isa T
    @test a == T(π / 2)

    a = ReferenceFrameRotations._mod_asin(-1.0f0 - eps(T))
    @test a isa T
    @test a == T(-π / 2)
end


text:
"""
Abstract decomposition of N-mode tensor.
"""
abstract type TensorDecomposition{T<:Number, N} end

"""
Returns relative error between the re-composed and the original tensors.
"""
rel_residue(recomposed::StridedArray{T, N}, tensor::StridedArray{T, N}) where {T,N} =
    norm(recomposed .- tensor) / norm(tensor)

"""
Returns relative error between re-composed and the original tensor.
"""
rel_residue(decomp::TensorDecomposition{T, N}, tensor::StridedArray{T, N}) where {T,N} =
    rel_residue(compose(decomp), tensor)

"""
Returns relative tensor decomposition error, NaN if not available.
"""
rel_residue(decomp::TensorDecomposition) = get(decomp.props, :rel_residue, NaN)

_set_rel_residue(decomp::TensorDecomposition, error::Float64) = setindex!(decomp.props, error, :rel_residue)

_set_rel_residue(decomp::TensorDecomposition{T,N}, tensor::StridedArray{T,N}) where {T,N} =
    _set_rel_residue(decomp, rel_residue(decomp, tensor))


text:
module TestJulia
using Test
using CxxCall

@testset "tocxx" begin
    dir = mktempdir()
    lib = joinpath(dir, "libtest_tocxx")
    filepath = joinpath(dir, "test_tocxx.cxx")
    eval(cxxsetup())
    eval(cxxnewfile(filepath, ""))
    seen = Set(String[])
    types = [
        Cstring,
        Cuchar,
        Cuint,
        Cchar,
        Cdouble,
        Cfloat,
        Cvoid,
        Cwchar_t,
        Cint,
        Cptrdiff_t,
        Clong,
        Clonglong,
        Cssize_t,
        Culong,
        Csize_t,
        Cshort,
        Cwstring,
        Culonglong,
        Cushort,
       ]
    for T in types
        @test CxxCall.tocxx(T) isa String
        Tcxx = tocxx(T)
        if !(Tcxx in seen)
            @cxx lib function id(obj::Ptr{T})::Ptr{T}
                "return obj;"
            end
            push!(seen, Tcxx)
        end
    end
    cxx_write_code!()
    @test isfile(filepath)
    libpath = "$(lib).so"
    run(`g++ -shared -fPIC $(filepath) -o $(libpath)`)
    using Libdl
    Libdl.dlopen(libpath) do _
        for T in types
            obj = Ptr{T}(C_NULL)
            @test id(obj) === obj
        end
    end
end

@testset "parse_fdef" begin
    parse_fdef = CxxCall.parse_fdef

    ex1 = :(function M.f()::Res
          body
      end
    )

    def = parse_fdef(ex1)
    @test isempty(def.args)
    @test def.fun == :(M.f)
    @test def.return_type == :Res

    @test parse_fdef(:(f()::Y = y)) ==
        parse_fdef(:(function f()::Y y end))
    @test parse_fdef(:(f(x::X)::Y = y)) ==
        parse_fdef(:(function f(x::X)::Y y end))

    ex2 = :(function f(
            x::X,
            y::Y)::Z
                some
                body
            end
    )
    parse_fdef(ex2)
    def = parse_fdef(ex2)
    @test def.fun == :f
    @test def.return_type == :Z
    @test def.args == [
     (val = :x, ann = :X),
     (val = :y, ann = :Y),
    ]

    ex = :(function f(
            x::X,
            y::Y)::Z
                for i in 1:10
                    return i
                end
            end
    )
    @test_throws Exception parse_fdef(ex)

    ex = :(function f(
            y::Y)::Z where {Y}
            end
    )
    @test_throws Exception parse_fdef(ex)
end
end#module